企业官网建设流程全解析

简介：这份PDF实验报告面向高校随机信号分析与处理课程的学习者，聚焦离散时间随机过程仿真这一核心难点。报告以正态随机过程为例，完整呈现从均匀分布随机数生成、白噪声正态序列构造，到具有特定相关函数与功率谱的随机过程x(n)仿真的全过程，并给出MATLAB代码与运行结果。读者可据此掌握rand函数、Box-Muller变换、集合统计法计算均值方差与相关函数、区间积分概率验证等关键方法，同时对照理论值分析误差来源。资源包共1个PDF文件，约525KB，内容涵盖实验目的、要求、程序代码、结果与体会五大模块，结构完整，适合作为课程实验参考或自学仿真练习的对照材料。目前已有183人学习下载，便于快速理解随机过程数值特征与编程实现之间的对应关系。

1. 从一份实验报告说起：正态随机过程仿真到底在算什么

很多人第一次拿到「相关正态随机过程的仿真」这个题目，第一反应是打开 MATLAB 敲randn，画个直方图就交差。但如果你真把这份实验报告里的五个任务拆开看，会发现它其实在逼你回答一个更底层的问题：给定一个相关函数和功率谱，怎么用均匀分布随机数一步步造出一个满足该统计特性的离散随机过程，并且用集合统计去验证它。这不是调库能糊弄过去的，因为实验明确要求从rand出发，经过 Box-Muller 变换生成白噪声，再通过一阶 AR 递推构造有色噪声，最后用时间平均代替集合平均去估计均值、方差和相关函数。适合正在学随机信号分析与处理、需要交实验报告或者想搞懂仿真链路每一环的从业者。下面我按任务顺序把代码、参数和踩坑点全部拆开。

2. 均匀分布到正态分布：Box-Muller 变换的代码实现与参数校验

2.1 为什么不能直接调 randn

实验要求第一条就写死了：利用计算机语言的 [0,1] 区间均匀分布随机数产生函数生成两个相互独立的序列。这意味着你不能用randn一步到位，必须自己实现从均匀分布到正态分布的映射。常见做法是 Box-Muller 变换，它把两个独立均匀分布变量 u1、u2 映射成两个独立标准正态变量。MATLAB 里rand产生的是 [0,1) 区间的均匀分布，注意是左闭右开，理论上 u1 可以取到 0，而log(0)是负无穷，这就是第一个坑。

% 任务1：生成两个相互独立的均匀分布序列并画直方图 clc; clear; N = 100000; u1 = rand(N,1); % 列向量，[0,1) 均匀分布 u2 = rand(N,1); % 与 u1 独立 subplot(121); histogram(u1); title('u1 分布'); subplot(122); histogram(u2); title('u2 分布');

这段代码本身没有难度，但要注意rand(N,1)和rand(1,N)在后续矩阵运算里的维度差异。任务 2 用的是行向量，任务 1 用的是列向量，混用会导致.*报错。参数 N 取 100000 是实验给定的，样本量足够大时直方图才能看出均匀性。如果你把 N 降到 1000，直方图会像锯齿一样，别急着怀疑代码，先加样本量。

2.2 Box-Muller 公式的代码落地

任务 2 要求生成均值 m=0、根方差 σ=1 的白色正态分布序列。实验报告里给出的公式是en = sqrt(-2*log(u1)).*cos(2*pi*u2)，这是 Box-Muller 的一种形式，只取了 cos 分支。严格来说 Box-Muller 一次产生两个独立正态变量，另一个是sqrt(-2*log(u1)).*sin(2*pi*u2)，但实验只要求一个序列，所以只用 cos 分支没问题。

% 任务2：Box-Muller 变换生成标准正态白噪声 clc; clear; u1 = rand(1,100000); u2 = rand(1,100000); en = sqrt(-2*log(u1)).*cos(2*pi*u2); % 标准正态分布 e(n) histogram(en,100); % 100 个 bin 看分布形状 title('e(n) 直方图');

逻辑说明：-2*log(u1)把均匀分布映射成指数分布，开根号后得到瑞利分布，再乘以cos(2*pi*u2)相当于给瑞利分布加上均匀相位，结果就是标准正态。参数说明：log是自然对数，u1不能为 0，实际运行时rand返回 0 的概率极低但并非不可能，稳妥做法是加一个极小量u1 = max(u1, eps)。histogram(en,100)的 100 是 bin 数量，bin 太少看不出正态钟形，bin 太多会毛刺，100 是实验报告里的经验值。

提示：如果你用 Python 复现，numpy.random.rand同样返回 [0,1)，Box-Muller 写法一致，但要注意np.log对 0 会返回-inf并触发警告，建议先做截断。

3. 一阶 AR 递推生成有色噪声：相关函数与功率谱的对应关系

3.1 递推公式的推导与初始条件

任务 3 是整个实验的核心。实验给出离散随机过程 x(n) 服从均值 mx=0、根方差 σx=2，相关函数和功率谱有特定形式。报告里给出的生成方法是x(n+1) = a*x(n) + 2*sqrt(1-a*a)*en(n+1)，初始条件x(1) = 2*sqrt(1-a*a)*en(1)。这里 a=0.6 是 AR(1) 模型的系数，决定了相邻样本的相关性。

% 任务3：一阶 AR 递推生成相关正态随机过程 clc; clear; u1 = rand(1,100000); u2 = rand(1,100000); en = sqrt(-2*log(u1)).*cos(2*pi*u2); % 白噪声 a = 0.6; x = zeros(1,100000); x(1) = 2*sqrt(1-a*a)*en(1); % 初始条件，保证方差从第一点就稳定 for n = 1:100000-1 x(n+1) = a*x(n) + 2*sqrt(1-a*a)*en(n+1); end histogram(x,100); title('x(n) 直方图');

逻辑说明：这个递推式本质是 AR(1) 过程，a是自回归系数，2*sqrt(1-a*a)是驱动噪声的增益。为什么系数是 2 而不是 1？因为实验要求根方差 σx=2，而en的根方差是 1，所以驱动项要放大 2 倍。sqrt(1-a*a)是为了让稳态方差归一化，推导一下：如果 x(n) = ax(n-1) + be(n)，稳态方差 Var(x) = b²/(1-a²) * Var(e)，代入 Var(e)=1、Var(x)=4，得 b² = 4(1-a²)，所以 b = 2*sqrt(1-a²)。初始条件用同样的增益是为了避免前几个样本方差偏小，如果你直接令 x(1)=0，前几十个点的方差会明显低于理论值，集合统计时均值可能没问题但方差会偏小。

参数说明：a=0.6对应相关函数 r(k) = σx² * a^|k| = 4 * 0.6^k，功率谱是典型的低通形状。a 越接近 1，相关性越强，功率谱越窄；a 越接近 0，越接近白噪声。实验报告里任务 4 计算的相关函数理论值q(p) = ax*ax*0.6.^p正是这个公式，其中 ax 是估计的标准差。

3.2 集合统计验证：均值、方差与相关函数

任务 4 要求用集合统计的方法计算均值、标准差和相关函数，并与理论值比较。这里有个概念要澄清：严格意义上的集合统计需要多条独立实现的样本轨道，但实验里只有一条长序列，所以实际上是用时间平均代替集合平均。对于各态历经的平稳随机过程，时间平均收敛到集合平均，AR(1) 过程满足这个条件。

% 任务4：时间平均估计均值、标准差和相关函数 sum_x = 0; for i = 1:100000 sum_x = sum_x + x(i); end mx = sum_x / 100000; % 均值估计 sum_x2 = 0; for i = 1:100000 sum_x2 = sum_x2 + x(i)*x(i); end ax = sqrt(sum_x2 / 100000); % 标准差估计 r = zeros(1,4); for k = 1:4 sum_r = 0; for j = 1:100000-k sum_r = sum_r + x(j)*x(j+k); end r(k) = sum_r / (100000-k); % 相关函数估计 end q = zeros(1,4); for p = 1:4 q(p) = ax*ax * 0.6^p; % 理论相关函数 end k = 1:4; hold on; plot(k, r, 'o-'); % 估计值 plot(k, q, 'x--'); % 理论值 legend('估计相关函数','理论相关函数');

逻辑说明：均值估计用sum/100000，标准差估计用sqrt(sum(x²)/100000)，注意这里除以 N 而不是 N-1，因为实验报告里就是这么写的，而且 N=100000 时两者差异可以忽略。相关函数估计r(k) = sum(x(j)*x(j+k))/(N-k)，分母是 N-k 而不是 N，这是无偏估计的常见做法。理论值q(p) = ax² * 0.6^p里的 ax 用的是估计标准差而不是理论值 2，这样比较更公平，能看出估计误差。

参数说明：k 取 1 到 4 是实验报告给定的范围，实际做的时候可以扩展到 10 或 20，观察相关函数是否按指数衰减。如果你发现 r(k) 在 k 较大时波动很大，那是因为 N-k 变小导致估计方差增大，这是正常现象，不是代码 bug。

注意：实验报告里任务 4 的代码有个笔误，sum=0之后在第二个循环里没有重置 sum，导致标准差计算错误。上面代码里我用了sum_x2单独变量，避免这个坑。如果你直接抄报告里的代码，标准差会算成累加值，结果完全不对。

4. 区间积分与概率验证：从直方图到数值积分的排查

4.1 四个区间的比例计算

任务 5 要求根据生成的 x(n)，在 100000 个数据中计算 (-∞,-2)、[-2,0]、(0,2]、(2,∞) 四个区间上数据出现的比例 P1、P2、P3、P4，并与理想值比较。因为 x(n) 服从均值 0、根方差 2 的正态分布，标准化后 Z = x/2 服从标准正态。区间 (-∞,-2) 对应 Z < -1，概率是 0.5 - Φ(1) ≈ 0.1587；[-2,0] 对应 -1 ≤ Z ≤ 0，概率是 Φ(1) - 0.5 ≈ 0.3413；后面两个区间对称。

% 任务5：统计四个区间的数据比例 num1 = 0; num2 = 0; num3 = 0; num4 = 0; for i = 1:100000 if x(i) < -2 num1 = num1 + 1; elseif x(i) >= -2 && x(i) <= 0 num2 = num2 + 1; elseif x(i) > 0 && x(i) <= 2 num3 = num3 + 1; else num4 = num4 + 1; end end p1 = num1/100000; p2 = num2/100000; p3 = num3/100000; p4 = num4/100000; disp('实验值：'); disp([p1,p2,p3,p4]); % 数值积分计算理论值 p2_theory = 0; for i = 1:200000 z = i * 0.00001; p2_theory = p2_theory + 1/(sqrt(2*pi)*2)*exp(-z*z/(2*2*2))*0.00001; end p3_theory = p2_theory; p1_theory = (1 - 2*p2_theory)/2; p4_theory = p1_theory; disp('理论值：'); disp([p1_theory,p2_theory,p3_theory,p4_theory]);

逻辑说明：统计部分用 if-elseif 链，注意边界条件，x(i) >= -2 && x(i) <= 0包含了 -2 和 0，而x(i) < -2不包含 -2，这样四个区间不重叠也不遗漏。数值积分部分，实验报告里用步长 0.00001 从 0 积到 2，被积函数是正态分布密度函数1/(sqrt(2*pi)*2)*exp(-z²/(2*2²))，注意这里的 2 是根方差 σx=2，不是方差。积分上限 200000*0.00001=2，正好覆盖 [0,2] 区间。

参数说明：步长 0.00001 对应 200000 个积分点，精度足够。如果你把步长改成 0.001，积分点降到 2000，p2 的误差会到小数点后三位，可能影响比较结论。实验报告里 p2 和 p3 理论值相同是因为正态分布关于均值对称，p1 和 p4 相同也是对称性。

4.2 结果差异的排查思路

实验报告给出的结果是均值 -0.0099、方差 1.9963，与理论值 0 和 4 差异很小。但如果你自己跑出来均值偏差超过 0.05，先别怀疑算法，按下面顺序排查：

第一，检查en的生成是否用了同一组 u1、u2。如果你在任务 3 里重新生成了 u1、u2，那 x(n) 和任务 2 的 en 就不是同一组白噪声，但这不影响统计特性，只影响具体数值。第二，检查 AR 递推的循环范围，for n = 1:100000-1生成 x(2) 到 x(100000)，共 99999 个点，加上 x(1) 正好 100000 个。如果你写成for n = 1:100000，会访问 x(100001) 导致数组越界。第三，检查相关函数估计的分母，用 N-k 还是 N，两者在 k 较小时差异不大，但 k 接近 N 时差异显著。

提示：如果你用 Python 的numpy.random.rand和numpy.log，注意np.log(0)会返回-inf并触发 RuntimeWarning，虽然结果里-inf会被后续运算吞掉，但建议加u1 = np.clip(u1, 1e-12, None)保持干净。

5. 避坑与常见问题：从代码报错到统计偏差的排查清单

5.1 现象：直方图不是正态钟形，而是尖峰或双峰

原因：Box-Muller 变换里u1或u2的维度不对，导致.*广播出错。比如u1是 1×100000 行向量，u2是 100000×1 列向量，MATLAB 会做隐式扩展生成 100000×100000 矩阵，内存爆炸且结果完全错误。解决：统一用行向量或列向量，推荐全部用rand(1,N)。

5.2 现象：均值接近 0 但方差明显偏小

原因：AR 递推的初始条件写成了x(1) = 0或者x(1) = en(1)，没有乘2*sqrt(1-a*a)。前几十个样本的方差会从 0 逐渐爬升到稳态值，如果序列不够长，整体方差就被拉低了。解决：初始条件必须和递推式的驱动项增益一致，即x(1) = 2*sqrt(1-a*a)*en(1)。

5.3 现象：相关函数估计值在 k=1 时就偏离理论值很远

原因：r(k)的分母用了 N 而不是 N-k，或者循环里sum_r没有在每次 k 循环开始时清零。实验报告里的代码在任务 4 就有这个笔误，sum变量被复用导致累加。解决：每个 k 循环内部重新初始化累加变量，分母用N-k。

5.4 现象：区间比例 P1 和 P4 不相等

原因：统计区间边界写错了，比如x(i) <= -2和x(i) < -2混用，导致 -2 这个点被重复计数或遗漏。虽然 -2 恰好等于边界值的概率极低，但逻辑上要严谨。解决：统一用左闭右开或左开右闭，四个区间互斥且完备。

5.5 现象：数值积分理论值与查表值对不上

原因：积分步长太大，或者被积函数里的 σ 写成了 σ²。正态分布密度函数是1/(sqrt(2*pi)*σ) * exp(-(x-m)²/(2σ²))，分母是 σ 不是 σ²，指数里是 2σ²。实验报告里 σ=2，所以分母是sqrt(2*pi)*2，指数里是2*2*2=8。解决：对照公式逐项检查，步长建议不超过 0.0001。

6. 进阶技巧：用向量化替代循环，把 100000 点仿真压到毫秒级

上面所有代码都用了 for 循环，因为实验报告是这么写的，方便对照。但实际做项目时，100000 点的 AR 递推用循环在 MATLAB 里要跑零点几秒，如果做蒙特卡洛需要几千次独立实现，循环就成了瓶颈。我一般会改成向量化或者用filter函数。

% 向量化生成 AR(1) 过程：用 filter 替代 for 循环 clc; clear; N = 100000; a = 0.6; u1 = rand(1,N); u2 = rand(1,N); en = sqrt(-2*log(max(u1,eps))).*cos(2*pi*u2); b = 2*sqrt(1-a*a); % filter 的分子是 [b]，分母是 [1, -a]，对应差分方程 x(n) = a*x(n-1) + b*e(n) x = filter(b, [1, -a], en); % 注意 filter 默认初始条件为 0，前几个点有暂态，可以丢弃前 100 个点 x = x(101:end); histogram(x,100); title('向量化生成的 x(n)');

逻辑说明：filter(b, [1, -a], en)实现的是差分方程x(n) - a*x(n-1) = b*e(n)，即x(n) = a*x(n-1) + b*e(n)，和 for 循环完全等价。但filter默认零初始条件，所以前几个样本有暂态过程，方差从 0 爬升。丢弃前 100 个点后，剩余序列的统计特性与稳态一致。参数说明：b是驱动增益，[1, -a]是分母系数，注意符号，filter的分母系数对应x(n) - a*x(n-1)，所以是-a不是a。

验证方法：跑一遍向量化代码，再跑一遍 for 循环代码，比较两者的均值和方差，差异应该在 0.01 以内。如果差异大，检查filter的系数符号是否写反。另外，filter的暂态长度取决于 a，a=0.6 时大约 20 个点就衰减完了，丢弃 100 个点足够安全。

从那以后我每次做随机过程仿真，都会先用 for 循环写一版对照实验报告，再用filter或向量化写一版用于批量跑数据，两版结果对不上就先查初始条件和系数符号。这个习惯帮我省了很多次返工。希望帮到你。

本文还有配套的精品资源，点击获取