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进行随机序列数字特征估计的方法。对之前的知识理解加深。
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1