实验一随机过程地模拟与特征估计.docx

1、实验一随机过程地模拟与特征估计实验一：随机过程的模拟与特征估计一、实验目的了解随机过程特征估计的基本概念和方法，学会运用MATLAB软件产生各种随机过程，对随机过程的特征进行估计，并通过实验了解不同估计方法所估计出来的结果之间的差异。二、实验原理（1）高斯白噪声的产生利用MATLAB函数randn产生（2）自相关函数的估计MATLAB自带的函数为xcorr()，阐述xcorr的用法 R=xcorr(x,y)或R=xcorr(x,y,option) 用来求序列x(n)与y(n)的互相关函数R=xcorr(x)或R=xcorr(x,option) 用来求序列x(n)的自相关函数option选项是:

2、biased有偏估计， unbiased无偏估计， coeff m=0 的相关函数值归一化为1none不作归一化处理（3）功率谱的估计利用周期图方法估计功率谱，提示：MATLAB自带的函数为periodogram()，阐述periodogram()的用法; 阐述其它谱估计方法的用法。Pxx,w=periodgram(x) Pxx为对应频率w的功率谱密度值。Pxx,w=periodgram(x,window) window=boxcar(n)矩形窗（RectangleWindow）window=triang(n)三角窗（TriangularWindow）window=hanning(n)汉宁窗（

3、HanningWindow）window=hamming(n)海明窗（HammingWindow）window=blackman(n)布拉克曼窗（BlackmanWindow）window=kaiser(n,beta)恺撒窗（KaiserWindow）Window代表与x等长度的窗序列，对数据进行加窗。其它谱估计方法：相关函数法(BT法)该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。r=xccor(x);R=fft(r);Pxx=abs(R);平均周期图法和平滑平均周期图法对于周期图的功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若

4、N太小,谱的分辨率又不好,因此需要改进。两种改进的估计法是平均周期图法和平滑平均周期图法。Bartlett法:Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。Welch法:Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w(n),并在周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠,这样会使方差减小。现代谱估计-非参数法 MTM法pmtm：使用正交窗口来截取获得相互独立的改进周期图法功率谱估计，然后再把这些估计结果结合得到最终的估计。随着NW的增大，窗的个数增多，会有更多的谱估计，从

5、而谱估计的方差得到减小，但同时带来谱泄露的增大，而且正的谱估计的结果将会有更大的偏差。MUSIC法pmusic：基于矩阵特征分解的谱估计非参数方法，它把相关数据矩阵中的信息分类，把信息分配到信号的子空间或噪声的子空间。它适合于普遍情况下的正弦信号参数估计的方法，是多信号分类法的简称。特征向量法peig：也是一种基于矩阵特征分解的谱估计非参数方法，它主要适用于混有噪声的正弦信号的功率谱估计，此方法利用相关矩阵的特征值来对MUSIC法公式中的求和进行加权得到的。Matlab函数（4）均值的估计 MATLAB自带的函数为mean()（5）方差的估计MATLAB自带的函数为var()(6)AR(1)模

6、型的理论自相关函数和理论功率谱对于AR(1)模型，自相关函数为 ，其功率谱为。三、实验容（带*为选作）1. 相关高斯随机序列的产生按如下模型产生一组随机序列，其中为均值为1，方差为4的正态分布白噪声序列。（1）产生并画出a=0.8和a=0.2的x(n)的波形；（2）估计x(n)的均值和方差；（3）估计x(n)的自相关函数，并画出相关函数的图形当a=0.8时程序：clc,clear all;a=0.8;w =1+sqrt(4)* randn(1,1000);%w（n） 均值为1 方差为4的正态分布白噪声x(1)=w(1)/(sqrt(1-a2);%x(1)初值条件for n=2:1000 x(n

7、)=a*x(n-1) +w(n);endsubplot(2,1,1);plot (x);%画出x（n）的波形title(x（n）波形图);axis(0 1000 -10 20);%设定坐标区间mean=mean(x)%估计x（n）的均值var=var(x)%估计x（n）的方差R=xcorr(x,coeff);%估计x（n）自相关函数,归一化subplot(2,1,2);plot(R);%画出自相关函数title(x(n)自相关函数);axis(0 2000 -0.3 1);结果：估计x（n）的均值及方差mean = 4.6889var = 10.8210X（n）的波形及其估计自相关函数图形：

8、当a=0.2时程序：clc,clear all;a=0.2;w =1+sqrt(4)* randn(1,1000);%w（n） 均值为1 方差为4的正态分布白噪声x(1)=w(1)/(sqrt(1-a2);%x(1)初值条件for n=2:1000 x(n)=a*x(n-1) +w(n);endsubplot(2,1,1);plot (x);%画出x（n）的波形title(x（n）波形图);axis(0 1000 -6 10);%设定坐标区间mean=mean(x)%估计x（n）的均值var=var(x)%估计x（n）的方差R=xcorr(x,coeff);%估计x（n）自相关函数,归一化su

9、bplot(2,1,2);plot(R);%画出自相关函数title(x(n)自相关函数);axis(0 2000 -0.3 1)结果：估计x（n）的均值及方差mean = 1.2322var = 3.9102X（n）的波形及其估计自相关函数图形： 对比a=0.8和a=0.2时的情况，当a0时，当a越大，x（n）变化幅度越大，变化越快，自相关函数下降越缓慢。 2. 两个具有不同频率的正弦信号的识别设信号为，其中为零均值正态白噪声，方差为。（1）假定，针对,和, 两种情况，使用周期图periodogram()的方法估计功率谱。取N=1000，由于本题中f0.5,故周期图函数使用默认采样频率Fs=

10、1Hz程序：clc,clear all;f1=0.05;f2=0.08;w =sqrt(1)* randn(1,1000);%w（n） 均值为0 方差为1的正态分布白噪声n=1:1000;x(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);subplot(2,1,1);periodogram(x);title(f1=0.05,f2=0.08时功率谱)f1=0.05;f2=0.20;y(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);subplot(2,1,2);periodogram(y);title(f1=0.05,f2=0.2

11、0时功率谱)f与图像上角频率对应关系从两幅图中可以看出，在f1和f2对应的频率上功率能量明显较大，第一幅图w1=0.1pi，w2=0.16pi，第二幅图w1=0.1pi，w2=0.2pi，且w2对应能量幅度大于w1.（2）假定,，针对和两种情况，用周期图periodogram()的方法估计功率谱程序clc,clear all;f1=0.05;f2=0.08;w1 =sqrt(1)* randn(1,1000);%w（n） 均值为0方差为1的正态分布白噪声w2 =sqrt(4)* randn(1,1000);%w（n） 均值为0方差为4的正态分布白噪声n=1:1000;x(n)=sin(2*pi

12、*f1*n)+2*cos(2*pi*f2*n)+w1(n);subplot(2,1,1);periodogram(x);title(噪声方差为1时功率谱)y(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w2(n);subplot(2,1,2);periodogram(y);title(噪声方差为4时功率谱) 对比两幅图像，可以看出，噪声方差越大能量越大。*(3) 假定, 选用不同的谱估计方法进行估计，并进行比较。clc,clear all;f1=0.05;f2=0.08;w =sqrt(4)* randn(1,1000);%w（n） 均值为0方差为4的正态分布白噪声n

13、=1:1000;x(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%自相关法R=xcorr(x);subplot(5,1,1);Pw=fft(R/2000);f=2*(0:length(Pw)-1)/length(Pw);plot(f,10*log10(abs(Pw);axis(0 1 -50 50 );title(理论x（n）功率谱);grid on;title(自相关法)subplot(5,1,2);periodogram(x);xlabel( );ylabel( );title(周期图法)subplot(5,1,3);window=hann(1000);

14、periodogram(x,window);xlabel( );title(汉宁周期图法)subplot(5,1,4)pwelch(x,512);xlabel( );ylabel( );title(平滑平均周期图法)subplot(5,1,5)pyulear(x,999);xlabel( );ylabel( );title(Yule-Walker方法)周期图法比自相关方产生图像更为平滑，汉宁窗周期图法产生的波形比无窗函数周期法减少了大量高次谐波，而平滑平均周期图法比汉宁窗周期图法产生波形更为平滑。从图中看出，f1和f2 对应频点能量明显高。3. 理论值与估计值的对比分析设有AR(1)模型，W(

15、n)是零均值正态白噪声，方差为4。 用MATLAB模拟产生X(n)的500个样本，并估计它的均值和方差； 画出X(n)的理论的自相关函数和功率谱； 估计X(n)的自相关函数和功率谱。clc,clear all;hold off;a=-0.8;w =1+sqrt(4)* randn(1,500);%w（n） 均值为1 方差为4的正态分布白噪声x(1)=w(1)/(sqrt(1-a2);%x(1)初值条件for n=2:500 x(n)=a*x(n-1) +w(n);endsubplot(111)plot (x);%画出x（n）的波形 title(x（n）波形图);axis(0 500 -20 2

16、0);%设定坐标区间mean=mean(x)%估计x（n）的均值var=var(x)%估计x（n）的方差for i=1:500; R(i)=4*(ai)/(1-a2);end;for j=1:999 if j500 R1(j)=R(j-499); end endendR2=xcorr(x);%估计x（n）自相关函数,归一化figure;subplot(2,1,1);i=1:999;plot(i-500,R1(i);%画出自相关函数title(理论x(n)自相关函数);subplot(2,1,2);plot(R2);%画出自相关函数title(估计x(n)自相关函数);figure;subplo

17、t(2,1,1);Pw=fft(R/50000);f=2*(0:length(Pw)-1)/length(Pw);plot(f,10*log10(abs(Pw);title(理论x（n）功率谱);subplot(2,1,2);periodogram(x,twosided);title(估计x（n）功率谱);axis(0 2 -40 40 )mean = 0.5657var = 11.5339可以看出估计的自相关函数与理论值大致吻合，但高频上有些噪声误差。但随着样本数的增多，估计曲线接近理论值。若采样点数为50004. 随机信号通过线性系统分析考虑图示系统假定w为正态分布的随机序列和均匀分布的随

18、机序列，分别估计输出的概率密度。Xn=wn+0.9wn-1-0.1wn-2w为均匀分布的随机序列clc,clear all;w=rand(1,1000000);x(1)=w(1);x(2)=w(2)+0.9*w(1);for n=3:1000000 x(n)=w(n)+0.9*w(n-1)-0.1*w(n-2);endsubplot(2,1,1);ksdensity(x);title(输出概率密度曲线);subplot(212);m=-0.5:0.1:2hist(x,m);title(输出序列直方图);w为正态分布的随机序列可以看出，当输入随机信号服从正态分布，通过线性系统后，输出的概率密度曲

19、线仍为正态分布；而输入随机信号服从均匀分布，输出概率密度曲线为三角形。四、实验思考题（1）自相关函数R(m)最大值应该在n=0，用MATLAB估计得到的结果与理论的结果相同吗？为什么？MATLAB 中，自相关函数最大值在序列中点。这是因为matlab中数组是从n=1开始定义的。实际上与理论结果一样，matlab计算出的R（m）相当于实际的函数向右平移一半的序列点数。（2）对两个正弦信号的识别，识别的效果如何？信噪比对识别效果有何影响？在第2问中，我们可以在估计功率谱看到两个明显凸起的频点，与题目中的正弦信号对应。两个频点相距越近时，需要一定量的采样点才能确保识别。信噪比越高，识别信号难度减低，所需采用点减少。5、心得体会 通过此次实验，掌握了使用matlab进行随机序列数字特征估计的方法。对之前的知识理解加深。