测试信号分析与处理作业实验五

王锋实验五：多种功率谱估计的比较一、实验目的a.了解功率谱估计在信号分析中的作用；b.掌握随机信号分析的基础理论，掌握参数模型描述形式下的随机信号的功率谱的计算方法；c.掌握在计算机上产生随机信号的方法；d.了解不同的功率谱估计方法的优缺点。二、实验准备有三个信号源，分别代表三种随机信号（序列）。信号源1：123()2cos(2)2cos(2)2cos(2)()xnfnfnfnzn其中，1230.08,=0.38,0.40fffz(n)是一个一阶AR过程，满足方程：()(1)(1)()znaznen(1)0.823321ae(n)是一高斯分布的实白噪声序列，方差20.1信号源2和信号源3：都是4阶的AR过程，它们分别是一个宽带和一个窄带过程，满足方程：()(1)(1)(2)(2)(3)(3)(4)(4)()xnaxnaxnaxnaxnene(n)是一高斯分布的实白噪声序列，方差2，参数如下：参数信号源a(1)a(2)a(3)a(4)2信号源2－1.3001.200－0.6000.2501.000信号源3－2.7503.799－2.6500.9281.000三、实验内容a.描绘出这三个实验信号的真实功率谱波形。b.在计算机上分别产生这个三个信号，令所得到的数据长度N=256。注意：产生信号的时候注意避开起始瞬态点。例如，可以产生长度为512的信号序列，然后取后面256个点作为实验数据。c.分别用如下的谱估计方法，对三个信号序列进行谱估计。1、经典谱估计周期图法自相关法平均周期图法（Bartlett法）Welch法（可选每段64点，重叠32点，用Hamming窗）2、现代谱估计Yule-Walker方程（自相关法）最小二乘法注：阶次p可在3－20之间，由自己给定。四、实验结果分析生成的信号源进行一次估计时的结果：信号1的经典谱估计：信号1的现代谱估计：信号2的经典谱估计：信号2的现代谱估计：信号3的经典谱估计：信号3的现代谱估计：运行50次求平均的结果：信号1的经典谱估计：信号1的现代谱估计：信号2的经典谱估计：信号2的现代谱估计：信号3的经典谱估计：信号3的现代谱估计：运行1次估计结果分析：可见现代谱估计的分辨率和光滑性都比较好。而经典谱估计中，周期图法和自相关法的分辨率较好，但是光滑性不好，起伏剧烈。Bartlett法和Welch法光滑性比较好，但是是分辨率降低了很多。运行50次结果分析：可以看到运行50次求平均后，周期图法以及自相关法的光滑性都得到了改善，并且没有降低分辨率。这是因为50次求平均对于平稳信号来说相当于进行了数据长度更长的welch法和平均周期图法，因此在没有牺牲分辨率的情况下，改善了光滑性。实验还发现，burg算法的结果很不稳定，分辨率有时极差，这可能与随机信号的采样点数太少有关。需要进步一研究。五、原程序清单%=======生成信号源1=======e1=wgn(1,512,0.1);%生成高斯白噪声a1=-0.823321;z1(1)=0;fori=2:512z1(i)=-a1.*z1(i-1)+e1(i);endf1=0.08;f2=0.38;f3=0.40;fori=1:512xn11(i)=2.*cos(2.*pi.*f1.*i)+2.*cos(2.*pi.*f2.*i)+2.*cos(2.*pi.*f3.*i)+z1(i);endxn1(1:256)=xn11(257:512);%取后半段256点subplot(3,1,1);plot(1:256,xn1)title('信号源1');axis([0256-1010])%=======生成信号源2=======e2=wgn(1,512,1);%生成高斯白噪声a21=-1.300;a22=1.200;a23=-0.600;a24=0.250;z2(1)=0;z2(2)=0;z2(3)=0;z2(4)=0;fori=5:512z2(i)=(-a21).*z2(i-1)+(-a22).*z2(i-2)+(-a23).*z2(i-3)+(-a24).*z2(i-4)+e2(i);endxn2(1:256)=z2(257:512);%取后半段256点subplot(3,1,2);plot(1:256,xn2)title('信号源2');axis([0256-1010])%=======生成信号源3=======e3=wgn(1,512,1);%生成高斯白噪声a31=-2.750;a32=3.799;a33=-2.650;a34=0.928;z3(1)=0;z3(2)=0;z3(3)=0;z3(4)=0;fori=5:512z3(i)=(-a31).*z3(i-1)+(-a32).*z3(i-2)+(-a33).*z3(i-3)+(-a34).*z3(i-4)+e3(i);endxn3(1:256)=z3(257:512);%取后半段256点subplot(3,1,3);plot(1:256,xn3)title('信号源3');axis([0256-100100])%====信号1的经典谱估计===m=50;%实验次数m取1和50fori=1:mX11(:,i)=fft(xn1,512);px11(:,i)=1/256*abs(X11(:,i)).^2;%周期图法Rx12(:,i)=xcorr(xn1);px12(:,i)=1/256*abs(fft(Rx12(:,i)));%自相关法xnk1(1:64)=xn1(1:64);xnk2(1:64)=xn1(65:128);xnk3(1:64)=xn1(129:192);xnk4(1:64)=xn1(193:256);Xnk1(:,i)=fft(xnk1,128);pxnk1(:,i)=1/64*abs(Xnk1(:,i)).^2;Xnk2(:,i)=fft(xnk2,128);pxnk2(:,i)=1/64*abs(Xnk2(:,i)).^2;Xnk3(:,i)=fft(xnk3,128);pxnk3(:,i)=1/64*abs(Xnk3(:,i)).^2;Xnk4(:,i)=fft(xnk4,128);pxnk4(:,i)=1/64*abs(Xnk4(:,i)).^2;px13(:,i)=1/4.*(pxnk1(:,i)+pxnk2(:,i)+pxnk3(:,i)+pxnk4(:,i));%4段平均周期图法px14(:,i)=pwelch(xn1,[],64,512);%welch法endpx111=mean(px11,2);px112=mean(px12,2);px113=mean(px13,2);px114=mean(px14,2);figure(4)subplot(4,1,1);plot(1:512,px111);title('信号1-周期图法');axis([02560300]);subplot(4,1,2);plot(1:511,px112);title('信号1-自相关法');axis([02560300]);subplot(4,1,3);plot(1:length(px113),px113);title('信号1-4段平均周期图法');axis([0640100]);subplot(4,1,4);plot(1:length(px114),px114);title('信号1-Welch法');axis([0256040]);%====信号1的现代谱估计===fori=1:mpx15(:,i)=pyulear(xn1,15,512);px16(:,i)=pburg(xn1,15,512);endpx115=mean(px15,2);px116=mean(px16,2);figure(5)subplot(2,1,1);plot(px115);title('信号1-Yule-Walker自相关法');subplot(2,1,2);plot(px116);title('信号1-Burg算法');%====信号2的经典谱估计===m=50;%实验次数m取1和50fori=1:mX21(:,i)=fft(xn2,512);px21(:,i)=1/256*abs(X21(:,i)).^2;%周期图法Rx22(:,i)=xcorr(xn2);px22(:,i)=1/256*abs(fft(Rx22(:,i)));%自相关法xnk21(1:64)=xn2(1:64);xnk22(1:64)=xn2(65:128);xnk23(1:64)=xn2(129:192);xnk24(1:64)=xn2(193:256);xnk221(:,i)=fft(xnk21,128);pxnk21(:,i)=1/64*abs(xnk221(:,i)).^2;xnk222(:,i)=fft(xnk22,128);pxnk22(:,i)=1/64*abs(xnk22(:,i)).^2;xnk223(:,i)=fft(xnk23,128);pxnk23(:,i)=1/64*abs(xnk223(:,i)).^2;xnk224(:,i)=fft(xnk24,128);pxnk24(:,i)=1/64*abs(xnk24(:,i)).^2;px23(:,i)=1/4.*(pxnk21(:,i)+pxnk22(:,i)+pxnk23(:,i)+pxnk24(:,i));%4段平均周期图法px24(:,i)=pwelch(xn2,[],64,512);%welch法endpx211=mean(px21,2);px212=mean(px22,2);px213=mean(px23,2);px214=mean(px24,2);figure(6)subplot(4,1,1);plot(1:512,px211);title('信号2-周期图法');axis([02560100]);subplot(4,1,2);plot(1:511,px212);title('信号2-自相关法');axis([02560100]);subplot(4,1,3);plot(1:length(px213),px213);title('信号2-4段平均周期图法');axis([064050]);subplot(4,1,4);plot(1:length(px214),px214);title('信号2-Welch法');axis([0256010]);%====信号2的现代谱估计===fori=1:mpx25(:,i)=pyulear(xn2,15,512);px26(:,i)=pburg(xn2,15,512);endpx215=mean(px25,2);px216=mean(px26,2);figure(7)subplot(2,1,1);plot(px215);title('信号2-Yule-Walker自相关法');subplot(2,1,2);plot(px216);title('信号2-Burg算法');%====信号3的经典谱估计===m=50;%实验次数m取1和50fori=1:mX31(:,i)=fft(xn3,512);pX31(:,i)=1/256*abs(X31(:,i)).^2;%周期图法Rx32(:,i)=xcorr(xn3);px32(:,i)=1/256*abs(fft(Rx32(:,i)));%自相关法xnk31(1:64)=xn3(1:64);xnk32(1:64)=xn3(65:128);xnk33(1:64)=xn3(129:192);xnk34(1:64)=xn3(193:256);xnk321(:,i)=fft(xnk31,128);pxnk31(:,i)=1/64*abs(xnk321(:,i)).^2;xnk322(:,i)=fft(xnk32,128);pxnk32(:,i)=1/64*abs(xnk32(:,i)).^2;xnk323(:,i)=fft(xnk3