随机过程模拟及特征估计

《随机过程模拟及特征估计》由会员分享，可在线阅读，更多相关《随机过程模拟及特征估计（20页珍藏版）》请在金锄头文库上搜索。

1、 本本 科科 实实 验验 报报 告告 实验名称： 随机过程的模拟与特征估计 一、实验目的 了解随机过程特征估计的基本概念和方法， 学会运用 MATLAB 软件产生各种随机过程， 对随机过程的特征进行估计，并通过实验了解不同估计方法所估计出来的结果之间的差异。 二、实验内容和原理 （一（一）实验实验原理原理 （1）高斯白噪声的产生 利用 MATLAB 函数 randn 产生： （2）自相关函数的估计 1 0 1 0 1 () ( ) ( ) 1 ( ) Nm n x Nm xn mn n x nm x n N R m R mxx Nm 对有偏估计 对无偏估计 MATLAB 中可利用 xcorr

2、函数估计随机过程的互相关序列。 c = xcorr(x,y)返回一个长度 2N-1 的向量 C （x， y 均为长度为 N 的向量， 若二者不等长， 则短的序列补零取齐） 。令 y=x，则可以求得随机序列的自相关序列。 通常相关函数需要标准化来产生精确估计。 c = xcorr(x,y,option)指定了标准化选项: a) biased有偏估计 b) unbiased无偏估计，无偏估计是系统误差为 0 的估计 c) coeff归一化，返回的向量最大值为 1 d) none未标准化 c = xcorr(x,y,maxlags,option) 中的 maxlags 可以控制延时长度最大为 max

3、lags （3）功率谱的估计 利用周期图方法估计功率谱， 21 ( )( ) x GX N Matlab 中的 periodogram()函数可以用来估计功率谱。 perdiodogram （x） 用来计算信号 x 的功率谱估计 （PSD estimate） 。 默认使用的是矩形窗。 如果 x 是实信号，则为双边功率谱估计。 periodogram(x,window)可以使用不同的窗（如哈明窗，汉宁窗等）进行平滑。 其它谱估计方法： 1.相关法谱估计（BT 功率谱） 以相关函数为媒介进行功率谱计算， 1 0 1 ( ) Nm xn mn n R mxx Nm ， 当仅使用长度为 2M-1（-M

4、+1 到 M-1）的自相关序列时，对其进行傅里叶变换： 1 +1 ( )( ) M jwm xx mM GR m e 这里用 fft 进行傅里叶变换 2.相关法谱估计 （4）均值的估计 1 1 1 ( ) N x n mx n N MATLAB 自带的函数为 mean()。 （5）方差的估计 1 22 1 1 ( ) N x n x nx N MATLAB 自带的函数为 var()。 (6) AR(1)模型的理论自相关函数和理论功率谱 对于 AR(1)模型( )(1)( )X naX nW n， 自相关函数为 2| | 2 ( ) 1 m X a Rm a ， 其功率谱 为 2 2 ( ) (

5、1) X j G ae 。 （二二）实验实验内容内容 1. 1. 相关高斯随机序列的产生相关高斯随机序列的产生 按如下模型产生一组随机序列( )(1)( )x nax nw n，其中( )w n为均值为 1，方差为 4 的正态分布白噪声序列。 （1）产生并画出 a=0.8 和 a=0.2 的 x(n)的波形； （2）估计 x(n)的均值和方差； （3）估计 x(n)的自相关函数，并画出相关函数的图形。 2. 两个具有不同频率的正弦信号的识别两个具有不同频率的正弦信号的识别 设信号为 12 ( )sin(2)2cos(2)( )x nf nf nw n，1,2,nN， 其中( )w n为零均 值

6、正态白噪声，方差为 2 。 （1）假定 2 1 ，针对 1 0.05f , 2 0.08f 和 1 0.05f , 2 0.20f 两种情况， 使用周期图 periodogram()的方法估计功率谱。 （2） 假定 1 0.05f , 2 0.08f ， 针对 2 1 和 2 4 两种情况， 用周期图 periodogram() 的方法估计功率谱 *(3) 假定 1 0.05f , 2 0.08f , 2 4, 选用不同的谱估计方法进行估计，并进行比较。 3. 理论值与估计值的对比分析理论值与估计值的对比分析 设有 AR(1)模型， ( )0.8(1)( )X nX nW n ， W(n)是零

7、均值正态白噪声，方差为 4。 用 MATLAB 模拟产生 X(n)的 500 个样本，并估计它的均值和方差； 画出 X(n)的理论的自相关函数和功率谱； 估计 X(n)的自相关函数和功率谱。 4. 4. 随机信号通过线性系统分析随机信号通过线性系统分析 考虑图示系统 其中 w 为均匀分布的随机序列，画出输出端的概率密度和直方图。 三、实验结果与分析 （一）（一）. . 相关高斯随机序列的产生相关高斯随机序列的产生 （1）画出（1）产生并画出 a=0.8 和 a=0.2 的 x(n)的波形，代码如下： %code1_1.m a=0.8;b=0.2;N=500; u=randn(N,1);v=ra

8、ndn(N,1); z-1z-1 + -0.1 0.9 wnxn z-1z-1 + -0.1 0.9 wnxn x(1)=2.*u(1)./sqrt(1-a2);%根据 AR 模型确定第一项 y(1)=2.*v(1)./sqrt(1-a2); w(1)=1+2.*u(1); for i=2:N w(i)=1+2.*u(i); x(i)=a*x(i-1)+w(i); y(i)=a*y(i-1)+w(i);%计算 x(n) end subplot(2,1,1)%画出 a=0.8 时 x(n)波形 plot(1:N,x) title(a=0.8 时 x(n) ylabel(x(n) xlabel(n

9、) 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.9837 a=0.2 时均值为 5.1801，方差为 12.2964 （3）估计 x(n)的自相关函数，并画出相关函数的图形。 以 a=0.8 为例，用 xcorr 函数计算自相关函数，产生序列并作图。代码如下： 010020030040050

10、0 -10 0 10 20 a=0.8时 x(n) Rx(n) n 0100200300400500 -10 0 10 20 a=0.2时 x(n) Rx(n) n %code1_3.m r1=xcorr(x); plot(t,r1) title(a=0.8 时 x(n)的自相关函数) ylabel(Rx(n) xlabel(n) 计算得 4 (0)max( 1)1.46 10Rxr 这与 222 (0)+=9.984.4029.26 XX Rxm矛盾 这是由于 xcorr 缺省计算的是无归一化的自相关序列。 对无偏估计求自相关函数： %code1_3.m t=-250:250; r2=xco

11、rr(x,250,unbiased); plot(t,r2) title(a=0.8 时 x(n)的自相关函数) ylabel(Rx(n) xlabel(n) Rxmax=29.2428 -500-400-300-200-1000100200300400500 0 5000 10000 15000 a=0.8时 x(n)的 自 相 关 函 数 Rx(n) n -250-200-150-100-50050100150200250 10 15 20 25 30 a=0.8时 x(n)的 自 相 关 函 数 Rx(n) n 可以看到此时跟自相关函数 R(0)相等。 （二二）. .两个具有不同频率的正

12、弦信号的识别两个具有不同频率的正弦信号的识别 (1) 不同频点的不同频点的功率谱功率谱估计估计 利用 periodogram()函数进行功率谱估计，代码如下： %code2_1.m N=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)

13、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,

14、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) 可以看出明显区分

15、出来两个频点，当10.05Hz,20.08Hzff时信号在 0.05Hz,0.08Hz 处信号功率谱密度达到最大。当10.05Hz,20.20Hzff时 信号在0.05Hz,0.20Hz 处信号功率谱密度达到最大。 (2) 不同方差的不同方差的功率谱功率谱估计估计 %code2_2.m N=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); w