简介本资源是湖南大学《随机信号分析与处理》课程配套实验报告面向电子信息、通信工程等专业本科生及信号处理初学者聚焦正态随机过程的建模仿真与统计特性验证。报告完整覆盖实验目的、要求、MATLAB实现代码含均匀/正态序列生成、自相关过程x(n)递推构造、均值方差计算、区间概率统计等5大任务、结果直方图与理论对比分析以及深入的实验体会助力读者理解随机过程核心概念并提升编程实践能力。压缩包为单个PDF文件525KB内容结构清晰含目录、原理说明、分步代码、运行结果图示及误差分析便于对照学习与复现。已有183人下载学习适合作为课程作业参考、实验预习材料或随机信号仿真入门范例。1. 正态随机过程仿真为什么不是“生成几个高斯数”就完事——从HNU实验一暴露的建模断层说起在某高校信号与系统类课程的实验一中“正态随机过程的仿真”常被学生简化为调用np.random.normal()生成一列数、画个直方图、算个均值方差就交作业。但真正卡住多数人的是实验指导书里那句轻描淡写的“观察样本函数的遍历性”——你生成了100条长度为1024的样本曲线却无法回答哪一条能代表整个过程的统计特性时间平均和集合平均为何不等自相关函数R(τ)的衰减形态怎么和功率谱密度S(f)对应上这些不是数学题而是信号建模的底层逻辑断层。本实验本质是搭建一个可验证、可复现、可拆解的正态随机过程黑匣子它必须同时满足严平稳性统计特性不随时间平移改变、各态历经性单一样本足够推断整体、且其功率谱需可控如白噪声、带限噪声、指数衰减相关。本文不讲定义只带你用Python从零搭起这个黑匣子跑通从时域样本生成→统计验证→频域检验→相关性调控的全链路并把HNU实验报告中最常被扣分的3个隐性陷阱直接摊开在代码注释里。2. 用确定性算法生成“真随机”正态随机过程的三重构造法正态随机过程不是“随机数序列”而是一族满足特定联合分布的随机变量集合。实验一要求的“仿真”核心在于构造出符合严平稳各态历经正态分布的样本函数族。常见做法是三层嵌套先设计功率谱密度决定相关性再通过滤波器实现决定物理可实现性最后驱动白噪声提供随机源。这比直接采样高斯分布更贴近工程实际——比如通信中接收机前端噪声就是带限高斯过程而非无限带宽白噪声。2.1 白噪声种子为什么必须用np.random.Generator而非np.random.randn旧式写法np.random.randn(N)已被弃用且无法复现相同随机序列。HNU实验要求多次运行结果一致便于批改对比必须使用现代NumPy的Generator接口import numpy as np # ✅ 正确可复现、线程安全、支持多种分布 rng np.random.default_rng(seed42) # seed固定所有后续调用可复现 white_noise rng.normal(loc0.0, scale1.0, size8192) # 标准正态白噪声 # ❌ 错误全局状态难追踪seed设置后仍可能被其他库污染 # np.random.seed(42) # white_noise np.random.randn(8192)参数说明loc0.0强制均值为0严平稳必要条件scale1.0对应单位方差后续滤波会缩放size8192取足够长序列以降低估计偏差HNU实验建议≥4096。seed42是教学常用固定值确保全班结果一致——这不是玄学是实验可重复性的底线。2.2 功率谱设计从“想要什么相关性”反推滤波器正态随机过程的统计特性由其自相关函数R(τ)完全决定而R(τ)与功率谱密度S(f)构成傅里叶变换对。实验一常要求实现指数衰减型自相关对应一阶低通滤波器输出或sinc型自相关对应矩形带通滤波器输出。我们以最典型的指数型为例$$ R(\tau) \sigma^2 e^{-\alpha|\tau|} \quad \Leftrightarrow \quad S(f) \frac{2\sigma^2\alpha}{\alpha^2 (2\pi f)^2} $$该S(f)物理可实现对应连续时间一阶RC低通系统。离散化后可用IIR滤波器逼近。关键参数alpha控制相关时间尺度alpha越大R(τ)衰减越快过程越“不相关”alpha越小记忆性越强相邻点越相似。from scipy import signal import numpy as np def design_exponential_filter(fs, alpha, order1): 设计离散化一阶低通滤波器使输出过程具有指数自相关 fs: 采样率HzHNU实验默认fs1000 alpha: 连续时间衰减系数rad/s需满足 alpha pi*fs 避免混叠 返回: b, a 系数用于scipy.signal.lfilter # 双线性变换s 2*fs*(z-1)/(z1) # 连续传递函数 H(s) alpha / (s alpha) # 离散化后 H(z) (alpha * (1z^-1)) / (2*fs alpha (2*fs - alpha)*z^-1) b0 alpha b1 alpha a0 2 * fs alpha a1 2 * fs - alpha return [b0, b1], [a0, a1] # HNU实验典型参数fs1000Hz, alpha200 rad/s → 相关时间约5ms fs 1000 alpha 200 b, a design_exponential_filter(fs, alpha) # 生成白噪声并滤波 rng np.random.default_rng(seed42) white_noise rng.normal(0, 1, 8192) gaussian_process signal.lfilter(b, a, white_noise) # 归一化方差至1便于后续统计验证 gaussian_process gaussian_process / np.std(gaussian_process)逻辑说明signal.lfilter是因果、稳定、实时可实现的IIR滤波器完美匹配“物理系统输出随机过程”的建模需求。b,a系数由双线性变换严格推导非经验拟合——这是HNU实验报告中区分“照抄代码”和“理解建模”的分水岭。gaussian_process此时已是严平稳正态过程每个时刻服从N(0,1)任意有限维联合分布为多元正态且R(τ)严格指数衰减。2.3 样本函数族生成一次滤波 ≠ 一个过程一个随机过程由无穷多个样本函数组成。实验一要求生成多条样本曲线如M50条每条长度N1024。关键陷阱不能对同一白噪声滤波50次必须为每条样本独立生成白噪声M 50 # 样本函数数量HNU实验要求≥20 N 1024 # 每条样本长度HNU实验要求≥512 # ✅ 正确每条样本使用独立随机种子 sample_functions np.zeros((M, N)) for i in range(M): rng_i np.random.default_rng(seed42 i) # 每条样本不同seed但可复现 white_i rng_i.normal(0, 1, N) sample_functions[i, :] signal.lfilter(b, a, white_i) sample_functions[i, :] / np.std(sample_functions[i, :]) # 逐条归一化 # ❌ 错误所有样本共享同一白噪声失去独立性 # white_common rng.normal(0, 1, N) # for i in range(M): # sample_functions[i, :] signal.lfilter(b, a, white_common)参数说明seed42i确保每条样本独立且可复现逐条归一化保证每条样本方差为1避免滤波器增益差异导致统计失真。此步生成的sample_functions是真正的“过程”——后续所有统计量均值、方差、自相关都需在此矩阵上计算。3. 统计验证遍历性不是口号是三个可计算的数值指标HNU实验报告扣分重灾区学生画出10条曲线就宣称“观察到遍历性”却未给出任何量化证据。遍历性Ergodicity指时间平均等于集合平均必须通过三组数值对比验证一阶矩遍历性单条样本的时间均值 ≈ 所有样本在t时刻的集合均值二阶矩遍历性单条样本的时间自相关 ≈ 所有样本在τ延迟下的集合自相关分布遍历性单条样本的直方图 ≈ 所有样本在任一时刻的集合直方图以下代码一次性完成全部验证并输出HNU实验报告要求的对比表格。import numpy as np from scipy import signal import matplotlib.pyplot as plt # 假设已生成 sample_functions (M50, N1024) M, N sample_functions.shape # 1. 一阶矩遍历性验证 time_mean np.mean(sample_functions[0, :]) # 第一条样本的时间均值 ensemble_mean np.mean(sample_functions[:, 0]) # t0时刻的集合均值任意t均可 # 更严谨计算所有t时刻集合均值的均值 ensemble_mean_all_t np.mean(sample_functions, axis0).mean() # 2. 二阶矩遍历性验证 # 时间自相关第一条样本的自相关无偏估计 time_acf np.correlate(sample_functions[0, :], sample_functions[0, :], modefull) time_acf time_acf[len(time_acf)//2:] / len(sample_functions[0, :]) # 只取τ≥0归一化 time_acf time_acf[:N//2] # 截断到合理延迟 # 集合自相关所有样本在延迟τ下的平均 ensemble_acf np.zeros(N//2) for tau in range(N//2): # 计算所有样本的x(t)x(ttau)平均 products [] for i in range(M): if tau N: prod sample_functions[i, :-tau] * sample_functions[i, tau:] products.extend(prod) ensemble_acf[tau] np.mean(products) # 3. 分布遍历性验证 # 单条样本直方图第一条 hist_time, bins np.histogram(sample_functions[0, :], bins50, densityTrue) # 集合直方图所有样本所有点合并 all_points sample_functions.flatten() hist_ensemble, _ np.histogram(all_points, binsbins, densityTrue) # 输出验证结果HNU实验报告核心表格 print( HNU实验一遍历性验证结果 ) print(f{指标:12} {时间平均:12} {集合平均:12} {相对误差(%):12}) print(- * 50) print(f{均值:12} {time_mean:12.4f} {ensemble_mean_all_t:12.4f} {abs(time_mean-ensemble_mean_all_t)/abs(ensemble_mean_all_t)*100:12.2f}) print(f{方差:12} {np.var(sample_functions[0,:]):12.4f} {np.var(sample_functions):12.4f} {abs(np.var(sample_functions[0,:])-np.var(sample_functions))/np.var(sample_functions)*100:12.2f}) # 计算ACF在τ0,1,10处的误差 tau_points [0, 1, 10] for tau in tau_points: if tau len(time_acf): err abs(time_acf[tau] - ensemble_acf[tau]) / (abs(ensemble_acf[tau]) 1e-8) print(f{ACF(τstr(tau)):12} {time_acf[tau]:12.4f} {ensemble_acf[tau]:12.4f} {err*100:12.2f})关键提示HNU实验报告明确要求“列出τ0,1,10处的ACF对比”。此处time_acf[0]即方差应≈1time_acf[1]反映相邻点相关性指数型应≈e⁻αᵀˢTₛ1/fstime_acf[10]体现长程记忆应显著衰减。若相对误差5%说明样本长度N不足或M太小——这是实验一最常见的翻车点根源常在2.3节的种子复用错误。4. 频域检验功率谱密度如何成为“过程指纹”正态随机过程的功率谱密度S(f)是其最本质的特征如同DNA。HNU实验一虽未明说但所有滤波器设计都隐含S(f)目标。验证S(f)是否符合预期是判断仿真是否成功的终极标准。这里不用plt.psd()这种黑盒函数而是用Welch法手动计算并与理论S(f)曲线叠加对比。4.1 Welch法实操窗口、重叠、FFT点数的三重权衡Welch法通过分段加窗、平均来降低功率谱估计方差。HNU实验要求采样率fs1000Hz故频率轴范围为[0,500]Hz。关键参数选择窗口长度取N_w256对应0.256秒平衡频率分辨率Δffs/N_w≈3.9Hz与方差抑制重叠率50%即每段移动128点提升统计稳定性FFT点数N_fft1024零填充提高谱线密度不影响分辨率from scipy import signal import numpy as np def welch_psd(x, fs, nperseg256, noverlap128, nfft1024): 手动实现Welch法返回f, Sxx f, Pxx signal.welch( x, fsfs, windowhann, # 汉宁窗抑制频谱泄漏 npersegnperseg, # 每段长度 noverlapnoverlap, # 重叠点数 nfftnfft, # FFT点数 detrendconstant, # 去直流分量严平稳要求均值为0 return_onesidedTrue ) return f, Pxx # 对第一条样本计算PSD f_welch, Sxx_welch welch_psd(sample_functions[0, :], fs1000) # 理论功率谱密度指数相关型 def theoretical_sxx_exponential(f, alpha, sigma21.0): 理论S(f) 2*sigma2*alpha / (alpha^2 (2*pi*f)^2) return 2 * sigma2 * alpha / (alpha**2 (2*np.pi*f)**2) Sxx_theory theoretical_sxx_exponential(f_welch, alpha200)4.2 可视化与误差分析为什么你的PSD图总“毛刺”plt.figure(figsize(10, 6)) plt.semilogy(f_welch, Sxx_welch, b-, labelWelch估计 (样本1), linewidth1.5) plt.semilogy(f_welch, Sxx_theory, r--, label理论S(f), linewidth2) plt.xlabel(频率 f (Hz)) plt.ylabel(功率谱密度 Sxx(f)) plt.title(正态随机过程功率谱密度验证) plt.legend() plt.grid(True, whichboth, ls-) plt.xlim(0, 200) # 聚焦低频区理论曲线在此最显著 plt.show() # 计算PSD相对误差HNU实验报告要求 # 取f∈[0,100]Hz区间计算均方相对误差 mask (f_welch 0) (f_welch 100) err_psd np.sqrt(np.mean(((Sxx_welch[mask] - Sxx_theory[mask]) / (Sxx_theory[mask] 1e-8))**2)) print(fPSD估计相对误差 (0-100Hz): {err_psd*100:.2f}%)避坑指南若PSD图出现剧烈毛刺尤其高频区90%概率是nperseg太小128或noverlap太低50%。HNU实验环境常因内存限制被迫减小nperseg此时必须同步增加M样本数来补偿——因为Welch法本质是用时间分集替代集合分集。这是学生最容易忽略的底层逻辑当无法获取足够长的单条样本时多条短样本的PSD平均等效于单条长样本的Welch估计。5. 避坑HNU实验一报告里最常被扣分的5个隐形陷阱学生常因“步骤正确但结果异常”被扣分根源在于未理解实验背后的信号建模约束。以下是某高校助教整理的5个血泪经验每条按“现象→原因→解决”结构呈现直击HNU实验报告评分细则。5.1 现象自相关函数R(τ)在τ0处值≠1且随τ增大不衰减原因未对滤波后样本进行逐条方差归一化。IIR滤波器增益不为1不同alpha导致输出方差差异巨大若仅对整个sample_functions矩阵归一化会破坏各条样本的独立性。解决在2.3节循环内对每条sample_functions[i, :]单独执行/ np.std(...)。验证np.std(sample_functions[0, :])必须≈1.0。5.2 现象Welch功率谱在f0处出现尖峰DC分量原因样本存在非零均值违反严平稳前提。np.random.normal(0,1)理论上均值为0但有限长度下存在估计偏差滤波过程可能引入微小直流偏移。解决在welch_psd函数中强制detrendconstant已包含或预处理x x - np.mean(x)。HNU实验明确要求“去除直流分量”。5.3 现象多条样本曲线在时域图中完全重叠无法分辨原因绘图时未设置alpha透明度或zorder且未添加图例。Matplotlib默认将多条线绘制在同一图层粗线覆盖细线视觉上变成单条。解决绘图时显式设置alpha0.3和zorderi例如for i in range(min(10, M)): # 只画前10条 plt.plot(sample_functions[i, :100], alpha0.3, zorderi, labelf样本{i1}) plt.legend(ncol2)5.4 现象集合自相关ensemble_acf[tau]在τ较大时震荡发散原因计算ensemble_acf时对每个τ使用了不同数量的样本点。例如τ100时每条样本只剩N-100个有效点若M50则总点数仅50×(N-100)统计不足。解决统一截断所有样本到相同长度。在生成sample_functions后执行# 截断为N_eff512确保大τ时仍有足够点数 N_eff 512 sample_functions sample_functions[:, :N_eff]5.5 现象实验报告结论写“过程具有遍历性”但未给出任何数值误差原因HNU实验评分标准第3条明确“遍历性验证需提供定量误差误差5%视为未通过”。学生常只画图不计算或只计算均值误差却忽略ACF误差。解决必须输出4.3节中的完整表格包含均值、方差、ACF(τ0,1,10)、PSD四类误差且所有误差需5%。若5%立即检查seed是否固定、M是否≥50、N是否≥1024。6. 进阶技巧用“相关时间”反推滤波器参数让实验报告脱颖而出HNU实验一的高分秘诀不在于复现标准流程而在于展示对物理意义的掌控力。例如指导书要求“生成相关时间为5ms的过程”但未给alpha值。此时你需要从定义出发相关时间τ_c定义为R(τ)衰减至R(0)/e的τ值。对指数型R(τ)σ²e⁻ᵃᵗ显然τ_c1/α。6.1 从物理指标反推数字滤波器参数假设HNU实验要求相关时间τ_c 5 ms 0.005 s采样率fs 1000 Hz → 采样间隔Tₛ 0.001 s则连续时间衰减系数$$ \alpha \frac{1}{\tau_c} 200 \text{ rad/s} $$但直接代入2.2节公式会因双线性变换引入频率畸变。更精确的做法是将期望的离散时间自相关R[k] e⁻ᵏᵀˢ/ᵀᶜ 与滤波器脉冲响应h[k]关联。一阶IIR滤波器h[k] β(1-β)ᵏk≥0其自相关为R_h[k] β²(1-β)ᵏ/(1(1-β)²)令其匹配R[k]即可解出β。def alpha_to_beta(alpha, fs, methodmatch_tau): 将连续时间alpha转换为离散IIR滤波器系数beta methodmatch_tau: 匹配相关时间推荐HNU实验适用 methodbilinear: 双线性变换理论严格但高频畸变 if method match_tau: # 离散相关时间 tau_c_discrete -T_s / ln(1-beta) # 令其等于连续tau_c则 beta 1 - exp(-T_s / tau_c) Ts 1 / fs beta 1 - np.exp(-Ts / (1/alpha)) # 1/alpha tau_c return beta else: # bilinear return alpha / (alpha 2*fs) # HNU实验tau_c5ms → beta 1 - exp(-0.001/0.005) 0.181 beta alpha_to_beta(alpha200, fs1000, methodmatch_tau) print(f匹配相关时间的IIR系数 beta {beta:.3f}) # 构造IIR滤波器H(z) beta / (1 - (1-beta)z^-1) b_iir [beta] a_iir [1, -(1-beta)]6.2 实验报告加分项动态调节相关时间并可视化在报告末尾添加一段代码展示tau_c从1ms到20ms变化时样本曲线形态如何从“剧烈抖动”变为“缓慢漂移”并同步更新ACF和PSD图。这证明你不仅会做更能解释“为什么”。taus_ms [1, 3, 5, 10, 20] fig, axes plt.subplots(2, 3, figsize(15, 10)) axes axes.flatten() for i, tau_ms in enumerate(taus_ms): alpha_i 1000 / tau_ms # tau_ms单位ms转为1/s beta_i alpha_to_beta(alpha_i, fs1000, methodmatch_tau) b_iir, a_iir [beta_i], [1, -(1-beta_i)] # 生成单条样本 rng_i np.random.default_rng(seed42i) x signal.lfilter(b_iir, a_iir, rng_i.normal(0,1,512)) x / np.std(x) # 画时域 axes[i].plot(x, b-, linewidth0.8) axes[i].set_title(fτ_c {tau_ms} ms) axes[i].grid(True) # 第6个子图理论ACF对比 tau_vec np.arange(0, 50) for tau_ms in taus_ms[:3]: # 只画3条避免混乱 alpha_i 1000 / tau_ms acf_theory np.exp(-tau_vec * 0.001 / (tau_ms/1000)) # T_s0.001s axes[5].plot(tau_vec, acf_theory, labelfτ_c{tau_ms}ms) axes[5].set_xlabel(延迟 k (采样点)) axes[5].set_ylabel(R[k]) axes[5].legend() axes[5].grid(True) plt.tight_layout() plt.show()这段代码的价值在于它把抽象的“相关时间”变成了可感知的波形特征。HNU实验报告中助教看到这张图会立刻认定你掌握了随机过程的核心物理意义——这比堆砌10行公式更有说服力。我带过的某跨平台系统项目中正是靠这类“物理可解释性”设计让算法在客户现场一次通过验收。希望帮到你。本文还有配套的精品资源点击获取