随机过程的模拟与特征估计

本科实验报告实验名称：随机过程的模拟与特征估计一、实验目的了解随机过程特征估计的基本概念和方法，学会运用MATLAB软件产生各种随机过程，对随机过程的特征进行估计，并通过实验了解不同估计方法所估计出来的结果之间的差异。二、实验内容和原理（一）实验原理（1）高斯白噪声的产生利用MATLAB函数randn产生：（2）自相关函数的估计10101()()ˆ()1ˆ()NmnxNmxnmnnxnmxnNRmRmxxNm对有偏估计对无偏估计MATLAB中可利用xcorr函数估计随机过程的互相关序列。c=xcorr(x,y)返回一个长度2N-1的向量C（x，y均为长度为N的向量，若二者不等长，则短的序列补零取齐）。令y=x，则可以求得随机序列的自相关序列。通常相关函数需要标准化来产生精确估计。c=xcorr(x,y,'option')指定了标准化选项:a)‘biased’有偏估计b)‘unbiased’无偏估计，无偏估计是系统误差为0的估计c)‘coeff’归一化，返回的向量最大值为1d)‘none’未标准化c=xcorr(x,y,maxlags,'option')中的maxlags可以控制延时长度最大为maxlags（3）功率谱的估计利用周期图方法估计功率谱，21ˆ()()xGXNMatlab中的periodogram()函数可以用来估计功率谱。perdiodogram（x）用来计算信号x的功率谱估计（PSDestimate）。默认使用的是矩形窗。如果x是实信号，则为双边功率谱估计。periodogram(x,window)可以使用不同的窗（如哈明窗，汉宁窗等）进行平滑。其它谱估计方法：1.相关法谱估计（BT功率谱）以相关函数为媒介进行功率谱计算，101ˆ()NmxnmnnRmxxNm，当仅使用长度为2M-1（-M+1到M-1）的自相关序列时，对其进行傅里叶变换：1+1ˆˆ()()MjwmxxmMGRme这里用fft进行傅里叶变换2.相关法谱估计（4）均值的估计111ˆ()NxnmxnNMATLAB自带的函数为mean()。（5）方差的估计12211ˆ[()]NxnxnxNMATLAB自带的函数为var()。(6)AR(1)模型的理论自相关函数和理论功率谱对于AR(1)模型()(1)()XnaXnWn，自相关函数为2||2()1mXaRma，其功率谱为22()(1)XjGae。（二）实验内容1.相关高斯随机序列的产生按如下模型产生一组随机序列()(1)()xnaxnwn，其中()wn为均值为1，方差为4的正态分布白噪声序列。（1）产生并画出a=0.8和a=0.2的x(n)的波形；（2）估计x(n)的均值和方差；（3）估计x(n)的自相关函数，并画出相关函数的图形。2.两个具有不同频率的正弦信号的识别设信号为12()sin(2)2cos(2)()xnfnfnwn，1,2,,nN，其中()wn为零均值正态白噪声，方差为2。（1）假定21，针对10.05f,20.08f和10.05f,20.20f两种情况，使用周期图periodogram()的方法估计功率谱。（2）假定10.05f,20.08f，针对21和24两种情况，用周期图periodogram()的方法估计功率谱*(3)假定10.05f,20.08f,24,选用不同的谱估计方法进行估计，并进行比较。3.理论值与估计值的对比分析设有AR(1)模型，()0.8(1)()XnXnWn，W(n)是零均值正态白噪声，方差为4。用MATLAB模拟产生X(n)的500个样本，并估计它的均值和方差；画出X(n)的理论的自相关函数和功率谱；估计X(n)的自相关函数和功率谱。4.随机信号通过线性系统分析考虑图示系统其中w为均匀分布的随机序列，画出输出端的概率密度和直方图。三、实验结果与分析（一）.相关高斯随机序列的产生（1）画出（1）产生并画出a=0.8和a=0.2的x(n)的波形，代码如下：%code1_1.ma=0.8;b=0.2;N=500;u=randn(N,1);v=randn(N,1);z-1z-1+-0.10.9w[n]x[n]z-1z-1+-0.10.9w[n]x[n]x(1)=2.*u(1)./sqrt(1-a^2);%根据AR模型确定第一项y(1)=2.*v(1)./sqrt(1-a^2);w(1)=1+2.*u(1);fori=2:Nw(i)=1+2.*u(i);x(i)=a*x(i-1)+w(i);y(i)=a*y(i-1)+w(i);%计算x(n)endsubplot(2,1,1)%画出a=0.8时x(n)波形plot(1:N,x)title('a=0.8时x(n)')ylabel('x(n)')xlabel('n')mean(x);var(x);subplot(2,1,2);%画出a=0.2时x(n)波形plot(1:N,y)title('a=0.2时x(n)')ylabel('x(n)')xlabel('n')（2）估计x(n)的均值和方差；用mean(),var()可以直接求得序列的均值和方差a=0.8时均值为4.3908，方差为9.9837a=0.2时均值为5.1801，方差为12.2964（3）估计x(n)的自相关函数，并画出相关函数的图形。以a=0.8为例，用xcorr函数计算自相关函数，产生序列并作图。代码如下：0100200300400500-1001020a=0.8时x(n)Rx(n)n0100200300400500-1001020a=0.2时x(n)Rx(n)n%code1_3.mr1=xcorr(x);plot(t,r1)title('a=0.8时x(n)的自相关函数')ylabel('Rx(n)')xlabel('n')计算得4(0)max(1)1.4610Rxr这与222(0)+=9.984.4029.26XXRxm矛盾这是由于xcorr缺省计算的是无归一化的自相关序列。对无偏估计求自相关函数：%code1_3.mt=-250:250;r2=xcorr(x,250,'unbiased');plot(t,r2)title('a=0.8时x(n)的自相关函数')ylabel('Rx(n)')xlabel('n')Rxmax=29.2428-500-400-300-200-1000100200300400500050001000015000a=0.8时x(n)的自相关函数Rx(n)n-250-200-150-100-500501001502002501015202530a=0.8时x(n)的自相关函数Rx(n)n可以看到此时跟自相关函数R(0)相等。（二）.两个具有不同频率的正弦信号的识别(1)不同频点的功率谱估计利用periodogram()函数进行功率谱估计，代码如下：%code2_1.mN=500;n=[1:500];f1=0.05;f2=0.08;f3=0.2;u=randn(N,1);w=u';x=sin(2*pi*f1*n)+2*cos(2*f2*pi*n)+w;y=sin(2*pi*f1*n)+2*cos(2*f3*pi*n)+w;subplot(2,3,1)plot(n,x)title('x(n)(f1=0.05,f2=0.08,δ=1)')ylabel('x(n)');xlabel('n')subplot(2,3,2)[Pww1,w1]=periodogram(x);w11=w1/2/pi;plot(w11,Pww1);title('x(n)功率谱(f1=0.05,f2=0.08,δ=1)')ylabel('Gx(w)(w/Hz)');xlabel('w(Hz)')subplot(2,3,3)Px1=10*log10(Pww1./w11);plot(w11,Px1);title('x(n)功率谱(f1=0.05,f2=0.08,δ=1)')ylabel('Gx(w)(dB/Hz)');xlabel('w(Hz)')subplot(2,3,4)plot(n,y)title('x(n)(f1=0.05,f2=0.20,δ=1)')ylabel('x(n)');xlabel('n')subplot(2,3,5);[Pww2,w2]=periodogram(y);w21=w2/2/pi;plot(w21,Pww2);title('x(n)功率谱(f1=0.05,f2=0.20,δ=1)')ylabel('Gx(w)(w/Hz)');xlabel('w(Hz)')subplot(2,3,6);Px2=10*log10(Pww2./w21);plot(w21,Px2);title('x(n)功率谱(f1=0.05,f2=0.20,δ=1)')ylabel('Gx(w)(dB/Hz)');xlabel('w(Hz)')可以看出明显区分出来两个频点，当10.05Hz,20.08Hzff时信号在0.05Hz,0.08Hz处信号功率谱密度达到最大。当10.05Hz,20.20Hzff时信号在0.05Hz,0.20Hz处信号功率谱密度达到最大。(2)不同方差的功率谱估计%code2_2.mN=500;n=[1:500];f1=0.05;f2=0.08;sigma1=1;sigma2=2;u=randn(N,1);v=randn(N,1);w1=sigma1*u';w2=sigma2*v';x=sin(2*pi*f1*n)+2*cos(2*f2*pi*n)+w1;y=sin(2*pi*f1*n)+2*cos(2*f2*pi*n)+w2;subplot(2,3,1)%x(n)波形plot(n,x)title('x(n)(f1=0.05,f2=0.08,δ=1)')ylabel('x(n)');xlabel('n')subplot(2,3,2)%x(n)功率谱[Pww1,w1]=periodogram(x);w11=w1/2/pi;plot(w11,Pww1);title('x(n)功率谱(f1=0.05,f2=0.08,δ=1)')ylabel('Gx(w)(w/Hz)');xlabel('w(Hz)')0200400600-505x(n)(f1=0.05,f2=0.08,δ=1)x(n)n00.20.40.60.8050100150x(n)功率谱(f1=0.05,f2=0.08,δ=1)Gx(w)(w/Hz)w(Hz)00.20.40.60.8-40-2002040x(n)功率谱(f1=0.05,f2=0.08,δ=1)Gx(w)(dB/Hz)w(Hz)0200400600-10-505x(n)(f1=0.05,f2=0.20,δ=1)x(n)n00.20.40.60.8050100x(n)功率谱(f1=0.05,f2=0.20,δ=1)Gx(w)(w/Hz)w(Hz)00.20.40.60.8-40-2002040x(n)功率谱(f1=0.05,f2=0.20,δ=1)Gx(w)(dB/Hz)w(Hz)subplot(2,3,3)%x(n)功率谱Px1=10*log10(Pww1);plot(w11,Px1);title('x(n)功率谱(f1=0.05,f2=0.08,δ=1)')ylabel('Gx(w)(dB/Hz)');xlabel('w(Hz)')subplot(2,3,4)%x(n)波形plot(n,y)title('x(n)(f1=0.05,f2=0.08,δ=2)')ylabel('x(n)');xlabel('n')subplot(2,3,5);%x(n)功率谱[Pww2,w2]=periodogram(y);w21=w2/2/pi;plot(w21,Pww2);title('x(n)功率谱(f1=0.05,f