实验 2 随机过程的计算机模拟 一、实验目的 1、给定功率谱（相关函数）和概率分布，通过计算机模拟分析产生相应的随机过程； 2、通过该随即过程的实际功率谱（相关函数）和概率分布验证该实验的有效性； 3、学会运用 Matlab 函数对随机过程进行模拟。 二、实验原理 1、标准正态分布随机序列的产生方法：利用随机变量函数变换的方法。 设 r1，r2 为两个相互独立的（0，1）均匀分布的随机数，如果要产生服从均值为 m,方差为 正态分布的随机数 x，则可以按如下变换关系产生： 2、正态随机矢量的模拟：设有一 N 维正态随机矢量 ，其概率密度为 为协方差矩阵，且是正定的。 3、具有有理谱的正态随机序列的模拟 根据随机过程通过线性系统的理论，白噪声通过线性系统后，输出是正态的，且输出功 率谱只与系统的传递函数有关。利用这一性质，我们可以产生正态随机过程。 如上图所示，输入 W(n)为白噪声，假定功率谱密度为 G (z) = 1 W ，通过离散线性系统 后，输出 X (n)是正态随机序列， 由于要求模拟的随机序列具有有理谱，则 G (z) X 可表示为： 其中，G (z) X+ 表示有理谱部分，即所有的零极点在单位圆之内，G (z) X? 表示非有理谱部分，即所有零极点在单位圆之外。 4、满足一定相关函数的平稳正态随机过程的模拟，当已知平稳随机过程的相关函数而要确定该随机过程的模拟算法。很显然，只要我们设计一个合适的滤波器，使得该白噪声通过滤波器后，输出的功率谱满足上述相关函数的傅里叶变换，就可以模拟得到该随机过程。 三、实验内容 1、产生两组相互独立的（0，1）均匀分布的随机数（随机数个数：500） 程序及图形如下: clear; x=randn(1,500); y= randn(1,500); subplot(2,1,1); plot(x); title(第二组); subplot(2,1,2); plot(y); title(第一组) 2、按照实验原理中的方法产生一组均值为 1，方差为 1 的正态分布的随机序列（序列长度：500） 程序及图形如下： clear; y=1+sqrt(1)*randn(1,500); plot(y); title(正态分布，均值方差都为 1) 3、画出功率谱密度为 G(w)=1/(1.25+cosw) 的功率谱图（一个周期内），采用均匀采样方法，采样点数为 500 程序及图形如下： clear; w=rand(1,500); M=1.25+cos(w); N=1; G=N./M; plot(G); title(均匀采样功率频谱); 5、模拟产生一个功率谱为 G(w)=1/(1.25+cosw) 的正态随机序列 程序及图形如下: clear; w=randn(1,500); M=1.25+cos(w); N=1; G=N./M; plot(G); title(均匀采样功率频谱); 4、实验中所遇到问题及解决方法 问题 1、对 Matlab 软件很生疏、编程也不熟悉。 解决方法：我和同学利用学校资源和网络查找和参看了许多有关 Matlab 的资料，对其中一些基本知识和有关随机信号处理的章节作了详细了解和练习。同时也向老师同学请教，经过实验课上和平时的练习，渐渐地有了一些好的进展，这个过程很漫长，但是很值得我们花时间和精力去了解。 问题 2、对 Matlab 中与统计信号处理随机过程中的某些函数的运用很有困难。 解决方法：也是查找 Matlab 的书籍中有关内容，然后在 Matlab 上学会用help 中的相关辅助、查找功能。 问题 3、运用 Matlab 编程时编写的程序经过运行之后有错，而且很难发现其中错误。 解决方法：充分利用 Matlab 运行出错之后的英文提示进行分析和改正，然后也要比较 Matlab 语言与 C 语言的差异和共同点，这样比较学习有利于我们更好地了解这门语言。 5、实验总结及心得体会 实验总结： 本实验是运用 Matlab 作为工具来对随机过程中的功率谱（相关函数）和概率分布函数进行模拟和验证，由于随机过程中涉及的数据和运运算往往比较繁多和复杂，运用 Matlab 这个软件的强大的数据、运算和图像处理功能可以很好的解决随机过程中的一些问题。 此外，这个实验也从实际动手的角度加深了我们对随机过程特征估计的理解，运用 Matlab 处理的图形也可以很深刻的帮助我们理解相关知识。 心得体会： 这次实验对于同学和我来说刚开始时是很不容易的，但是经过和同学的协作、查找参看一些相关资料，我们对 Matlab 的实际操作之后，我们还是有一些收获的，我们对 Matlab 有了进一步认识，对于随机过程这一重要内容也有所了解，对于统计信号分析的一些知识也不仅仅只是再停留在理论方面了，这次实验让我们以实际动手的方式去认知感受统计信号的知识。更重要的是我觉得，这次实验也在一定程度上，锻炼、提高我们通信工程专业学生的根据理论分析和实验工具来设计分析实验的思维和能力。 此外，我们发现做实验时理论知识也是很重要的，只有对理论知识有了很深的理解，这样才有可能运用课堂所学内容去指导实践，也只有这样的实践才会加强我们对知识的掌握程度。 实验 3 随机过程的特征估计 一、实验目的 1、了解随机过程特征估计的基本概念和方法； 2、学会运用 Matlab 函数对随机过程进行特征估计； 3、通过实验了解不同估计方法所估计出来结果之间的差异。 二、实验原理 设随机序列 X(n)、Y(n)为各态历经过程，样本分别为 x(n)、y(n)（n=0,1,.N-1）。 1、均值的估计 2、方差的估计 方差估计有两种情况，如果均值 mx 已知，则 如果均值未知，那么 3、相关函数估计 4、功率谱估计 功率谱的估计有几种方法， （1）自相关法:先求相关函数的估计， 然后对估计的相关函数做傅立叶变换， （2）周期图法:先对序列 x(n)做傅立叶变换， 则功率谱估计为 周期图法是一种非参数谱估计方法，另外还有一种修正的周期图方法，也叫Welch 法， MATLAB 有周期图和 Welch 法的谱估计函数。 （3）现代谱估计技术 现代谱估计主要有参数谱估计和子空间谱估计。参数谱估计法是假定待估计功率谱的信号是白噪声驱动线性系统的输出，常用的方法有基于最大墒估计的伯格算法和 Yuler-Walk 自回归(AR)方法，这些方法是估计线性系统的参数，通常会得到比经典谱估计方法更好的估计。子空间法也称为高分辨率谱估计或超分辨率谱估计，常用的方法有 MUSIC 法和特征矢量法，这些方法是根据相关矩阵的特征分析或特征分解得到对信号的频率分量的估计，特别适合于线谱（即正弦信号）的估计，是低信噪比环境下检测正弦信号的有效方法。 MATLAB 有许多估计数字特征的统计函数： 1.均值与方差 mean(A)，返回序列的均值，序列用矢量 A 表示。VAR(X)，返回序列 X 的方差。 2.互相关函数估计 xcorr 用法： c = xcorr(x,y) c = xcorr(x) c = xcorr(x,y,option) c = xcorr(x,option) xcorr(x,y)计算 X 与 Y 的互相关，矢量 X 表示序列 x(n)，矢量 Y 表示序列y(n)。xcorr(x)计算 X 的自相关。option 选项是： biased：有偏估计 unbiased:：无偏估计 coeff：m=0 的相关函数值归一化为 1。 none：不作归一化处理。 3.功率谱估计：MATLAB 提供了许多功率铺估计的函数： 三、实验内容 1、产生一组均值为 1，方差为 4 的正态分布的随机序列（1000 个样本），估计该序列的均值与方差。 程序及图形为： clc,clear N=1000;alp=1;sig=1;delt=1;a=1; mm=zeros(1,N);x1=rand(1,N);x2=rand(1,N); x3=a.*sqrt(-2*log(x1).*cos(2*pi*x2)+mm; %产生随高斯分布的随机数 y(1)=sig*x3(1); for n=2:N y(n)=exp(-alp)*y(n-1)+sig*sqrt(1-exp(-2*alp*delt)*x3(n); end i=1:N; plot(i,y); hold on; plot(i,mm,-);title(正态分布随机序列) M=0;for i=1:N M=M+y(i);endM=M/ND=0;for i=1:N D=D+(y(i)-M)2; end D=D/N for m=1:N %计算自相关函数正半轴% for n=1:N-m+1 rr(n)=y(n)*y(n+m-1); end r2(m)=sum(rr)/N;end M =-0.0061 D = 0.8837 2、按如下模型产生一组随机序列： x(n)=0.8x(n-1)+w(n) 其中 w(n)为均值为 1，方差为 4 的正态分布白噪声序列。估计过程的自相关函数与功率谱。程序及图形为： N=500; u=randn(N,1); w=1+2.*u; x(1)=w(1); for i=2:N; x(i)=0.8*x(i-1)+w(i); end subplot(2,2,1); plot(x); subplot(2,2,3); R=xcorr(x,coeff); plot(R); subplot(2,2,4); P=periodogram(x); plot(P); 3、设信号为 x(n)=sin(2f1n)+2cos(2f2n)+w(n)，n=1,2,.,N，其中f1=0.05，f2=0.12，w(n)为正态白噪声，试在 N=356 和 1024 点时，分别产生随机序列 x(n)、画出 x(n)的波形并估计 x(n)的相关函数和功率谱。 程序及图形为： Fs=1000; n=0:1/Fs:1; xn=sin(2*pi*0.05*n)+2*cos(2*pi*0.12*n)+randn(size(n);