简介：本资源是湖南大学《随机信号分析与处理》课程配套实验报告，面向电子信息、通信工程等专业本科生及信号处理初学者，聚焦正态随机过程的建模仿真与统计特性验证。报告完整覆盖实验目的、要求、MATLAB实现代码（含均匀/正态序列生成、AR模型构建x(n)、均值方差与自相关函数计算、区间概率统计等5大任务）、结果直方图与数值对比分析，以及深入的实验体会，助力读者理解随机过程核心概念并提升编程实践能力。资源为单文件PDF，大小525KB，内容结构清晰，含目录、代码分段注释、图表与理论值比对表格，便于精读与复现。已有183人学习下载，适合作为课程作业参考、实验预习材料或随机信号仿真入门实操范本。
1. 正态随机过程仿真：为什么“随机”不是乱来，而是可复现、可验证、可工程化的信号行为？
在某高校《随机与信号》实验课中，学生第一次接触“正态随机过程”时，常误以为只要调用np.random.normal()就算完成——结果跑出的样本路径波动剧烈、自相关函数不衰减、功率谱密度不成形，实验报告里画出的曲线被批注：“像噪声，不像过程”。这恰恰暴露了核心误区：正态随机过程 ≠ 一堆独立正态数的拼接，而是具有严格统计结构（均值函数、自相关函数、遍历性）的连续/离散时间函数族。本实验聚焦“相关正态随机过程”的仿真，即构造满足指定自相关函数 $R_X(\tau)$ 的宽平稳高斯过程样本，其本质是将理论中的“协方差矩阵开方”操作落地为可执行的数值生成流程。它直接支撑通信系统信道建模、传感器噪声注入、雷达回波仿真等真实场景。适合刚学完随机过程基础、手握 Python 环境、需要把课本公式（如维纳-辛钦定理、Cholesky 分解、谱表示法）真正跑通的本科生与入门工程师。别怕“随机”二字——这里每一步都有确定性算法，每一次“翻车”都对应一个可定位的数学假设偏差。
2. 从理论到代码：用三种主流方法生成相关正态随机过程
生成相关正态随机过程，本质是解决同一个问题：给定目标自相关函数 $R_X(\tau)$（或等价的功率谱密度 $S_X(f)$），构造长度为 $N$ 的离散时间序列 ${x[n]}_{n=0}^{N-1}$，使其经验统计量逼近理论要求。常见做法是三条技术路径：时域滤波法（最直观）、频域谱表示法（最高效）、协方差矩阵分解法（最通用）。三者数学等价，但实现细节、计算开销、适用边界差异显著。我一般会先用滤波法快速验证思路，再用谱法处理长序列，最后用矩阵法调试非平稳或特殊约束场景。下面逐个拆解，每种方法都给出最小可运行代码、关键参数说明及选择依据。
2.1 时域滤波法：用白噪声“染色”，让独立样本产生相关性
这是最符合直觉的方法：先生成零均值单位方差的独立同分布（i.i.d.）高斯白噪声 $w[n]$，再通过一个线性时不变（LTI）滤波器 $h[n]$，输出 $x[n] = w[n] * h[n]$。根据线性系统理论，输出过程的功率谱密度为 $S_X(f) = |H(f)|^2 S_W(f) = |H(f)|^2$（因 $S_W(f)=1$），故只需设计 $h[n]$ 使其频率响应幅值平方等于目标 $S_X(f)$。对指数衰减自相关 $R_X(\tau) = \sigma^2 e^{-\alpha|\tau|}$，其对应功率谱为洛伦兹型 $S_X(f) = \frac{2\sigma^2\alpha}{\alpha^2 + (2\pi f)^2}$，此时可用一阶 IIR 滤波器实现。
import numpy as np import matplotlib.pyplot as plt def generate_by_filtering(N, sigma2=1.0, alpha=0.2, fs=1.0): """ 用一阶IIR滤波器生成指数自相关的正态过程 :param N: 样本点数 :param sigma2: 目标过程方差 :param alpha: 自相关衰减系数（单位：1/s） :param fs: 采样率（Hz），决定时间分辨率 :return: 长度为N的数组x[n] """ # 生成白噪声 w = np.random.normal(0, 1, N) # 设计一阶IIR滤波器系数：y[n] = a1*y[n-1] + b0*w[n] # 理论推导得：a1 = exp(-alpha/fs), b0 = sqrt(sigma2 * (1 - a1**2)) a1 = np.exp(-alpha / fs) b0 = np.sqrt(sigma2 * (1 - a1**2)) # 初始化输出数组 x = np.zeros(N) # 初始状态设为0（忽略暂态影响，或预热足够长） for n in range(1, N): x[n] = a1 * x[n-1] + b0 * w[n] return x # 生成示例 N = 4096 x_filt = generate_by_filtering(N, sigma2=4.0, alpha=0.5, fs=10.0)逻辑说明与参数说明：
a1是滤波器极点位置，直接控制自相关衰减速率；alpha越大，a1越小，过程越“粗糙”（记忆性越弱）。b0的推导来自稳态方差匹配：$\mathrm{Var}(x) = \frac{b_0^2}{1-a_1^2} = \sigma^2$，故 $b_0 = \sigma \sqrt{1-a_1^2}$。若忽略此步，输出方差会严重偏离设定值。- 为何用 IIR 而非 FIR？因指数自相关对应无限长冲激响应，FIR 近似需极高阶数，而 IIR 用 2 个参数即可精确实现。但注意：IIR 滤波器相位非线性，若需严格线性相位（如某些通信场景），需改用 FIR 或全通补偿。
2.2 频域谱表示法：在频域“雕刻”功率，再逆变换回时域
当目标功率谱 $S_X(f)$ 形式复杂（如多峰、带限、有尖锐谐振），时域滤波器设计困难时，谱表示法成为首选。其核心是：对白噪声做傅里叶变换 → 在频域按 $\sqrt{S_X(f)}$ 缩放幅度 → 保持原相位 → 逆傅里叶变换。该方法天然保证输出为高斯过程（因线性变换不改变高斯性），且计算效率高（FFT 复杂度 $O(N\log N)$），特别适合生成万点以上长序列。
def generate_by_spectral(N, S_func, fs=1.0, sigma2_target=None): """ 用谱表示法生成相关正态过程 :param N: 样本点数（建议为2的幂次以加速FFT） :param S_func: 功率谱密度函数，输入为频率f（Hz），输出为S_X(f) :param fs: 采样率（Hz） :param sigma2_target: 若指定，则对输出进行方差归一化（因数值积分误差） :return: 长度为N的实数数组x[n] """ # 生成白噪声并FFT w = np.random.normal(0, 1, N) W = np.fft.rfft(w) # 实数FFT，返回N//2+1个复数 # 计算对应频率轴（rfft的频率点） freqs = np.fft.rfftfreq(N, d=1/fs) # 单位：Hz # 计算目标幅度谱（sqrt(S_X(f))），注意：rfft输出含直流和Nyquist分量 amplitude = np.sqrt(S_func(freqs)) # 缩放W的幅度（保持相位不变） W_scaled = W * amplitude # 逆变换并取实部（因输入为实数，IFFT后虚部应为0，取实部防浮点误差） x = np.fft.irfft(W_scaled, n=N).real # 方差归一化（可选） if sigma2_target is not None: current_var = np.var(x) if current_var > 1e-10: x = x * np.sqrt(sigma2_target / current_var) return x # 定义目标功率谱：带限白噪声（矩形谱） def S_rect(f, B=2.0, sigma2=4.0): """带宽B Hz的矩形功率谱，总功率=σ²""" return np.where(np.abs(f) <= B, sigma2 / (2*B), 0.0) # 生成示例 x_spec = generate_by_spectral(8192, lambda f: S_rect(f, B=1.5, sigma2=4.0), fs=10.0)逻辑说明与参数说明：
S_func必须定义在rfftfreq返回的频率点上，且需处理f=0（直流）和f=f_s/2（奈奎斯特）的特殊权重（rfft中直流和奈奎斯特分量只出现一次，其余频率成对共轭）。本例中S_rect已正确处理。sigma2_target是重要安全阀：因数值 FFT 和插值存在积分误差，直接缩放后的输出方差可能偏离理论值 5–10%。开启归一化可确保统计量达标，尤其在教学实验中避免因数值误差被质疑“结果不对”。- 为何用
rfft而非fft？因输入信号为实数，rfft自动利用共轭对称性，内存和计算量减半，且输出频率范围明确（0 到 $f_s/2$），避免手动处理负频率的繁琐。
2.3 协方差矩阵分解法：当自相关函数任意给定，且需精确控制有限维统计量
前两种方法依赖于 $R_X(\tau)$ 具有解析表达式（能写出 $S_X(f)$ 或设计滤波器）。但实验中常遇到：老师给一个离散表格形式的 $R_X[k]$（$k=0,\pm1,\pm2,...,\pm M$），或要求生成满足特定协方差结构（如空间相关、非均匀采样）的过程。此时，直接构造 $N\times N$ 协方差矩阵 $\mathbf{C}$，再对其 Cholesky 分解 $\mathbf{C} = \mathbf{L}\mathbf{L}^T$，最后用 $\mathbf{x} = \mathbf{L}\mathbf{w}$ 生成，是最普适的方案。它不依赖频域性质，但计算复杂度为 $O(N^3)$，仅适用于 $N<5000$ 的中小规模。
def generate_by_cholesky(R_func, N, sigma2=1.0): """ 用Cholesky分解生成相关正态过程 :param R_func: 自相关函数，输入为时延k（整数），输出为R_X[k] :param N: 样本点数 :param sigma2: 目标方差（用于校准R_func(0)） :return: 长度为N的数组x[n] """ # 构造N×N协方差矩阵C C = np.zeros((N, N)) for i in range(N): for j in range(N): k = abs(i - j) # 时延 C[i, j] = R_func(k) # 强制C(0,0) = sigma2，确保方差正确 C[0, 0] = sigma2 # Cholesky分解（要求C正定） try: L = np.linalg.cholesky(C) # L为下三角阵 except np.linalg.LinAlgError: # 若C非正定（如R_func定义不当），添加微小扰动 eps = 1e-8 * np.eye(N) L = np.linalg.cholesky(C + eps) # 生成白噪声向量 w = np.random.normal(0, 1, N) # 生成相关过程 x = L @ w # 矩阵乘法 return x # 定义自相关函数：三角窗（简单示例） def R_tri(k, M=5, sigma2=4.0): """三角自相关：R[k] = σ² * max(0, 1 - |k|/M)""" if k <= M: return sigma2 * max(0, 1 - k/M) else: return 0.0 # 生成示例（小规模） x_chol = generate_by_cholesky(R_tri, N=1024, sigma2=4.0)逻辑说明与参数说明：
R_func(k)输入必须是非负整数时延k，因为协方差矩阵是对称 Toeplitz 矩阵，C[i,j]只依赖|i-j|。若给定的是负时延值，需先转换为k=|i-j|。np.linalg.cholesky要求矩阵严格正定。实际中，当R_func截断过早或N过大时，C可能接近奇异。代码中加入eps扰动是工程常用技巧，其量级（1e-8）远小于主对角线元素，不影响统计特性。- 为何不用
eigendecomposition？虽然特征分解（$\mathbf{C}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^T$）更稳定（可处理半正定），但计算更慢且需处理负特征值截断。Cholesky 在正定时更快，且教学实验中N较小，稳定性足够。
3. 验证三板斧：如何证明你生成的真是“相关正态过程”？
生成只是第一步，验证才是实验报告得分的关键。很多同学只画一条样本路径就交差，却忽略了“过程”的统计本质——它是一族随机变量，需通过多次独立实现来估计统计量。验证必须包含三个维度：分布形态（是否正态）、时域相关性（是否匹配 $R_X(\tau)$）、频域特性（是否匹配 $S_X(f)$）。下面给出可直接复用的验证脚本框架与判据阈值。
3.1 分布验证：Q-Q 图比直方图更可靠，K-S 检验给定量结论
正态性检验不能只看直方图是否“钟形”——小样本下直方图波动大，且无法区分偏态与峰态。Q-Q 图（分位数-分位数图）是黄金标准：若样本分位数与理论正态分位数严格共线，则分布为正态。辅以 Kolmogorov-Smirnov（K-S）检验给出 p 值，p > 0.05 接受原假设（服从正态分布）。
from scipy import stats def validate_distribution(x, title=""): """ 验证单次实现x的正态性 :param x: 一维数组 :param title: 图标题前缀 """ fig, axes = plt.subplots(1, 2, figsize=(12, 4)) # Q-Q 图 stats.probplot(x, dist="norm", plot=axes[0]) axes[0].set_title(f"{title} Q-Q Plot") # K-S 检验 ks_stat, ks_p = stats.kstest(x, 'norm', args=(np.mean(x), np.std(x, ddof=1))) axes[0].text(0.05, 0.95, f'KS p-value = {ks_p:.4f}', transform=axes[0].transAxes, verticalalignment='top') # 直方图叠加正态拟合曲线 axes[1].hist(x, bins=50, density=True, alpha=0.7, label='Empirical') mu, std = np.mean(x), np.std(x, ddof=1) x_norm = np.linspace(mu-4*std, mu+4*std, 100) axes[1].plot(x_norm, stats.norm.pdf(x_norm, mu, std), 'r-', lw=2, label='Fitted Normal') axes[1].set_title(f"{title} Histogram & Fit") axes[1].legend() plt.tight_layout() plt.show() print(f"[{title}] Mean={mu:.4f}, Std={std:.4f}, KS p-value={ks_p:.4f}") # 对三种方法生成的样本分别验证 validate_distribution(x_filt, "Filtering Method") validate_distribution(x_spec, "Spectral Method") validate_distribution(x_chol, "Cholesky Method")判据解读：
- Q-Q 图中点越贴近直线，正态性越好；边缘点（尾部）偏离预示厚尾或薄尾。
- K-S p 值 > 0.05 是基本门槛；若 p < 0.01，大概率非正态（如滤波器初始状态未清零导致暂态污染）。
- 注意：单次实现的
np.mean(x)和np.std(x)是估计值，必然与理论值有偏差。验证重点是分布形状，而非均值/方差的绝对精度。
3.2 自相关函数验证：用np.correlate计算，但必须归一化与去偏
np.correlate(x,x,mode='full')输出的是未归一化、含偏置的互相关，直接绘图会因R_X(0)值巨大而掩盖其他时延的细节。正确做法是：除以R_X(0)归一化，并用scipy.signal.correlate的method='direct'确保精度（避免 FFT 方法在短序列下的泄漏）。
from scipy import signal def estimate_autocorr(x, max_lag=100): """ 估计自相关函数（归一化，无偏估计） :param x: 输入序列 :param max_lag: 最大时延（避免计算过长） :return: lag_arr, R_est（时延数组和估计值数组） """ # 使用无偏估计：除以 (N - |k|) 而非 N N = len(x) R_full = signal.correlate(x, x, mode='full', method='direct') # 提取中心部分（对应时延 -N+1 到 N-1） R_centered = R_full[N-1-max_lag : N+max_lag] # 长度 2*max_lag+1 lags = np.arange(-max_lag, max_lag+1) # 归一化：R_est[k] = R_full[k+N-1] / ((N - |k|) * var(x)) var_x = np.var(x, ddof=1) R_est = np.zeros_like(R_centered, dtype=float) for i, lag in enumerate(lags): denom = (N - abs(lag)) * var_x if denom > 1e-10: R_est[i] = R_centered[i] / denom else: R_est[i] = 0.0 return lags, R_est # 估计并绘图 lags_filt, R_filt = estimate_autocorr(x_filt, max_lag=50) lags_spec, R_spec = estimate_autocorr(x_spec, max_lag=50) lags_chol, R_chol = estimate_autocorr(x_chol, max_lag=50) plt.figure(figsize=(10, 6)) plt.plot(lags_filt, R_filt, 'b-', label='Filtering (Est.)', alpha=0.8) plt.plot(lags_spec, R_spec, 'g--', label='Spectral (Est.)', alpha=0.8) plt.plot(lags_chol, R_chol, 'r-.', label='Cholesky (Est.)', alpha=0.8) # 绘制理论曲线（以指数为例） tau_theory = np.arange(-50, 51) R_theory = np.exp(-0.5 * np.abs(tau_theory)) # alpha=0.5, sigma2=1 plt.plot(tau_theory, R_theory, 'k:', label='Theory: exp(-0.5|τ|)', lw=2) plt.xlabel('Lag τ') plt.ylabel('Autocorrelation R_X(τ)') plt.title('Estimated vs Theoretical Autocorrelation') plt.legend() plt.grid(True) plt.show()关键细节：
- 无偏估计：分母用
(N - |k|)而非N，因时延k下只有N-|k|对数据点可用于计算。这是教材强调但常被忽略的要点。method='direct'避免 FFT 法在max_lag << N时的频谱泄漏，保证时域估计精度。- 理论曲线必须与估计使用相同参数（如
alpha=0.5），否则对比无意义。实验报告中务必注明所有参数。
3.3 功率谱密度验证：Welch 方法是工业标准，参数设置决定成败
周期图法（直接对|X(f)|^2平均）方差大，Welch 方法通过分段、加窗、平均，显著降低方差，是 MATLABpwelch和 Pythonscipy.signal.welch的默认策略。但窗口长度、重叠率、FFT 点数设置不当，会导致分辨率与方差矛盾。
from scipy.signal import welch def validate_psd(x, fs=1.0, title=""): """ 用Welch法估计PSD并绘图 :param x: 输入序列 :param fs: 采样率 :param title: 标题 """ # Welch参数：nperseg=256（平衡分辨率与方差），noverlap=128（50%重叠），nfft=1024（零填充提升频率分辨率） f, Pxx = welch(x, fs=fs, nperseg=256, noverlap=128, nfft=1024, window='hann', scaling='density') plt.figure(figsize=(10, 5)) plt.semilogy(f, Pxx, label='Estimated PSD') # 绘制理论PSD（以指数自相关对应的洛伦兹谱为例） S_theory = (2 * 0.5 * 1.0) / (0.5**2 + (2*np.pi*f)**2) # alpha=0.5, sigma2=1.0 plt.semilogy(f, S_theory, 'r--', label='Theory: Lorentzian', lw=2) plt.xlabel('Frequency (Hz)') plt.ylabel('PSD (V²/Hz)') plt.title(f'{title} Power Spectral Density') plt.legend() plt.grid(True) plt.show() # 验证 validate_psd(x_filt, fs=10.0, title="Filtering Method")Welch 参数指南：
nperseg: 每段长度。越大则频率分辨率越高（$\Delta f = f_s / nperseg$），但段数减少，方差增大。推荐nperseg ≈ N/10到N/4。noverlap: 重叠点数。50% 重叠（noverlap=nperseg//2）是经典选择，兼顾平滑性与计算量。nfft: FFT 点数。大于nperseg时自动零填充，提升频率轴密度（插值），不提高真实分辨率，但使曲线更平滑易读。window: 汉宁窗（'hann'）是默认推荐，抑制频谱泄漏；矩形窗（'boxcar'）分辨率最高但泄漏严重，仅用于理论验证。
4. 避坑指南：这五个血泪经验，让我少改三次实验报告
在指导某高校《随机与信号》实验课的三年中，我整理出学生在本实验中最常踩的五个坑。每个坑都对应一个具体现象、根本原因和可立即执行的解决方案。它们不是玄学，而是数学假设与工程实现之间的缝隙。
4.1 现象：生成的样本路径看起来“太光滑”或“太毛刺”，自相关函数在零时延处有尖锐峰值，但衰减过快或过慢
原因：滤波器设计时，a1和b0的推导未考虑采样率fs的尺度效应。例如，将连续时间指数衰减率α（单位：1/s）直接代入离散时间公式a1 = e^{-α}，而忽略了a1 = e^{-α T_s} = e^{-α / f_s}。当fs=1时侥幸正确，fs=10时a1被错误放大 10 倍，导致过程记忆性过强。
解决：所有连续时间参数（α,B等）必须显式除以fs再代入离散公式。在代码中用注释标明单位，如# alpha: 0.5 s⁻¹, so discrete pole at exp(-0.5/10) for fs=10Hz。
4.2 现象：Q-Q 图显示明显 S 形弯曲（尾部上翘），K-S 检验 p 值 < 0.001，但直方图看似正常
原因：生成过程中使用了未初始化的滤波器状态。例如，在generate_by_filtering中，x[0]被设为 0，但x[1] = a1*x[0] + b0*w[1]仍受初始零状态影响，前 10–20 点存在暂态，污染整体分布。尤其当a1接近 1（慢衰减）时，暂态持续更久。
解决：预热（warm-up）滤波器。生成远超N的序列（如2*N），丢弃前N//2点，再取后N点。或改用scipy.signal.lfilter并设置zi初始条件为稳态方差对应值。
4.3 现象：谱表示法生成的序列，其估计 PSD 在高频端出现异常凸起，或总功率（积分）远大于sigma2
原因：rfft的频率轴freqs与功率谱函数S_func的定义域不匹配。例如，S_func被写成lambda f: 1.0 if abs(f)<1 else 0.0，但rfft的freqs包含f=0和f=f_s/2，而f_s/2处的功率应为S_X(f_s/2)/2（因奈奎斯特分量只计一次），但代码中未做此修正。
解决：对rfft输出的功率谱进行加权。在generate_by_spectral中，对freqs[0]（直流）和freqs[-1]（奈奎斯特，当N为偶数时存在）的amplitude乘以1/sqrt(2)，其余点保持不变。这是scipy.signal.periodogram的内部处理逻辑。
4.4 现象：Cholesky 分解报错LinAlgError: Matrix is not positive definite，尤其当N>2000或R_func截断较早时
原因：由离散R_func构造的 Toeplitz 矩阵C在数值上接近奇异。理论中，若R_func是某个平稳过程的自相关，则C必正定；但截断后的R_func（如只定义k=0..10，其余为 0）破坏了正定性，因它隐含了“过程在k>10后完全不相关”的强假设，与连续谱矛盾。
解决：用正则化（regularization）替代简单报错。在generate_by_cholesky中，当cholesky失败时，不抛出异常，而是添加eps * eye(N)并重试。eps取1e-8 * C[0,0]量级，既保证分解成功，又不显著改变统计特性。这是数值计算中的标准实践。
4.5 现象：多次运行同一段代码，生成的样本路径完全不同，但实验报告要求“可复现”
原因：np.random的全局随机种子未固定。每次运行np.random.normal都使用不同种子，导致结果不可复现，无法对比不同方法或调试问题。
解决：在脚本开头强制设置全局种子。添加np.random.seed(42)（或任何固定整数）。更佳实践是创建独立Generator：rng = np.random.default_rng(42)，然后用rng.normal()。这样避免干扰其他库的随机行为，且default_rng是 NumPy 1.17+ 推荐方式。
5. 进阶技巧：用“遍历性验证”把实验报告升格为小型研究
教科书说“宽平稳过程若具各态历经性，则时间平均等于集合平均”，但学生常止步于背诵。本实验的终极价值，是让你亲手验证这一抽象概念。下面这个技巧，能让报告脱颖而出：用单次长序列的时间平均，去逼近理论上的集合平均统计量，并量化误差。它不增加代码量，却体现对随机过程本质的深刻理解。
5.1 时间平均 vs 集合平均：构造一个可量化的验证闭环
集合平均（ensemble average）需无穷多次独立实现，现实中不可行。但遍历性允许我们用单次实现的长时间观测来估计。例如：
- 均值函数：理论为常数 $\mu_X$，时间平均为 $\frac{1}{N}\sum_{n=0}^{N-1} x[n]$
- 自相关函数：理论为 $R_X(\tau)$，时间平均为 $\frac{1}{N-|\tau|}\sum_{n=0}^{N-1-|\tau|} x[n]x[n+|\tau|]$
- 功率谱密度：理论为 $S_X(f)$，时间平均为 Welch 法估计的 PSD
关键在于：比较时间平均与理论值的相对误差，并观察其随 $N$ 增大的收敛趋势。这比单次验证更有说服力。
def ergodicity_test(x, R_theory_func, fs=1.0, N_list=[1024, 4096, 16384]): """ 遍历性验证：测试不同长度N下，时间平均统计量对理论值的逼近程度 :param x: 超长序列（长度 >= max(N_list)） :param R_theory_func: 理论自相关函数 R_X(k) :param N_list: 测试的序列长度列表 :return: 字典，含各N下的误差 """ results = {} for N in N_list: x_short = x[:N] # 时间平均均值 mu_time = np.mean(x_short) mu_theory = 0.0 # 假设零均值 mu_error = abs(mu_time - mu_theory) / (abs(mu_theory) + 1e-6) # 时间平均自相关（在tau=10处） tau_test = 10 if N > tau_test: R_time = np.mean(x_short[:-tau_test] * x_short[tau_test:]) R_theory = R_theory_func(tau_test) R_error = abs(R_time - R_theory) / (abs(R_theory) + 1e-6) else: R_error = np.nan # Welch PSD 在中心频率f=0.5Hz处的误差 f, Pxx = welch(x_short, fs=fs, nperseg=min(256, N//4), noverlap=128, nfft=1024) idx_f05 = np.argmin(np.abs(f - 0.5)) S_time = Pxx[idx_f05] if idx_f05 < len(Pxx) else np.nan S_theory = (2 * 0.5 * 1.0) / (0.5**2 + (2*np.pi*0.5)**2) # 同前 S_error = abs(S_time - S_theory) / (abs(S_theory) + 1e-6) if not np.isnan(S_time) else np.nan results[N] = { 'mu_error': mu_error, 'R_error_at_10': R_error, 'S_error_at_0p5Hz': S_error } return results # 生成超长序列（10万点） x_long = generate_by_spectral(100000, lambda f: S_rect(f, B=1.5, sigma2=4.0), fs=10.0) # 执行遍历性测试 ergo_results = ergodicity_test(x_long, R_tri, fs=10.0) # 打印结果表格 print("遍历性验证结果（相对误差 %）:") print(f"{'N':<8} {'Mean Error':<12} {'R(10) Error':<12} {'S(0.5Hz) Error':<15}") print("-" * 55) for N, errs in ergo_results.items(): print(f"{N:<8} {errs['mu_error']*100:<12. <p> <a href="https://download.csdn.net/download/m0_66182318/88532430" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>
