王锋实验五：多种功率谱估计的比较一、实验目的a.了解功率谱估计在信号分析中的作用；b.掌握随机信号分析的基础理论，掌握参数模型描述形式下的随机信 号的功率谱的计算方法；c.掌握在计算机上产生随机信号的方法；d.了解不同的功率谱估计方法的优缺点。

二、实验准备有三个信号源，分别代表三种随机信号（序列）。

信号源1：123()2cos(2)2cos(2)2cos(2)()x n f n f n f n z n πππ=+++其中，1230.08,=0.38,0.40f f f == z(n)是一个一阶 AR 过程，满足方程： ()(1)(1)()z n a z n e n =--+ (1)0.823321a =-e(n)是一高斯分布的实白噪声序列，方差20.1σ=信号源2和信号源3：都是4阶的AR 过程，它们分别是一个宽带和一个窄带过程，满足方程： ()(1)(1)(2)(2)(3)(3)(4)(4)()x n a x n a x n a x n a x n e n =--------+ e(n)是一高斯分布的实白噪声序列，方差2σ，参数如下：三、实验内容a. 描绘出这三个实验信号的真实功率谱波形。

b. 在计算机上分别产生这个三个信号，令所得到的数据长度 N= 256 。

注意：产生信号的时候注意避开起始瞬态点。

例如，可以产生长度为512 的信号序列，然后取后面256 个点作为实验数据。

c. 分别用如下的谱估计方法，对三个信号序列进行谱估计。

1、经典谱估计 周期图法 自相关法平均周期图法（Bartlett 法）Welch法（可选每段64 点，重叠32 点，用Hamming 窗）2、现代谱估计Yule - Walker方程（自相关法）最小二乘法注：阶次p可在3－20之间，由自己给定。

四、实验结果分析生成的信号源进行一次估计时的结果：信号1的经典谱估计：信号1的现代谱估计：信号2的现代谱估计：信号3的现代谱估计：运行50次求平均的结果：信号1的经典谱估计：信号1的现代谱估计：信号2的现代谱估计：信号3的现代谱估计：运行1次估计结果分析：可见现代谱估计的分辨率和光滑性都比较好。

而经典谱估计中，周期图法和自相关法的分辨率较好，但是光滑性不好，起伏剧烈。

Bartlett法和Welch法光滑性比较好，但是是分辨率降低了很多。

运行50次结果分析：可以看到运行50次求平均后，周期图法以及自相关法的光滑性都得到了改善，并且没有降低分辨率。

这是因为50次求平均对于平稳信号来说相当于进行了数据长度更长的welch 法和平均周期图法，因此在没有牺牲分辨率的情况下，改善了光滑性。

实验还发现，burg算法的结果很不稳定，分辨率有时极差，这可能与随机信号的采样点数太少有关。

需要进步一研究。

五、原程序清单%=======生成信号源1=======e1=wgn(1,512,0.1);%生成高斯白噪声a1=-0.823321;z1(1)=0;for i=2:512z1(i)=-a1.*z1(i-1)+e1(i);endf1=0.08;f2=0.38;f3=0.40;for i=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([0 256 -10 10])%=======生成信号源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;for i=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([0 256 -10 10])%=======生成信号源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;for i=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([0 256 -100 100])%====信号1的经典谱估计===m=50;%实验次数m取1和50for i=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*a bs(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*a bs(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([0 256 0 300]);subplot(4,1,2);plot(1:511,px112);title('信号1-自相关法');axis([0 256 0 300]);subplot(4,1,3);plot(1:length(px113),px113);title('信号1-4段平均周期图法');axis([0 64 0 100]); subplot(4,1,4);plot(1:length(px114),px114);title('信号1-Welch法');axis([0 256 0 40]);%====信号1的现代谱估计===for i=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和50for i=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([0 256 0 100]);subplot(4,1,2);plot(1:511,px212);title('信号2-自相关法');axis([0 256 0 100]);subplot(4,1,3);plot(1:length(px213),px213);title('信号2-4段平均周期图法');axis([0 64 0 50]); subplot(4,1,4);plot(1:length(px214),px214);title('信号2-Welch法');axis([0 256 0 10]);%====信号2的现代谱估计===for i=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和50for i=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(xnk33,128);pxnk33(:,i)=1/64*abs(xnk323(:,i)).^2;xnk324(:,i)=fft(xnk34,128);pxnk34 (:,i)=1/64*abs(xnk34(:,i)).^2;px33(:,i)=1/4.*(pxnk31(:,i)+pxnk32(:,i)+pxnk33(:,i)+pxnk34(:,i));%4段平均周期图法px34(:,i)=pwelch(xn3,[],64,512);%welch法endpX311=mean(pX31,2);pX312=mean(px32,2);pX313=mean(px33,2);pX314=mean(px34,2);figure(8)subplot(4,1,1);plot(1:512,pX311);title('信号3-周期图法');axis([0 256 0 40000]);subplot(4,1,2);plot(1:511,pX312);title('信号3-自相关法');axis([0 256 0 40000]);subplot(4,1,3);plot(1:length(pX313),pX313);title('信号3-4段平均周期图法');axis([0 64 0 20000]); subplot(4,1,4);plot(1:length(pX314),pX314);title('信号3-Welch法');axis([0 256 0 6000]);%====信号3的现代谱估计===for i=1:mpx35(:,i)=pyulear(xn3,15,512);px36(:,i)=pburg(xn3,15,512);endpX315=mean(px35,2);pX316=mean(px36,2);figure(9)subplot(2,1,1);plot(pX315);title('信号3-Yule-Walker自相关法');subplot(2,1,2);plot(pX316);title('信号3-Burg算法');六、实验后的体会和建议通过此次实验，了解不同的功率谱估计方法的优缺点，更加熟悉了MATLAB 相关的操作。