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

王锋

实验五：多种功率谱估计的比较

一、实验目的

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:512

z1(i)=-a1.*z1(i-1)+e1(i);

end

f1=0.08;f2=0.38;f3=0.40;

for i=1:512

xn11(i)=2.*cos(2.*pi.*f1.*i)+2.*cos(2.*pi.*f2.*i)+2.*cos(2.*pi.*f3.*i)+z1(i);

end

xn1(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:512

z2(i)=(-a21).*z2(i-1)+(-a22).*z2(i-2)+(-a23).*z2(i-3)+(-a24).*z2(i-4)+e2(i);

end

xn2(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:512

z3(i)=(-a31).*z3(i-1)+(-a32).*z3(i-2)+(-a33).*z3(i-3)+(-a34).*z3(i-4)+e3(i);

end

xn3(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和50

for i=1:m

X11(:,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法

end

px111=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:m

px15(:,i)=pyulear(xn1,15,512);

px16(:,i)=pburg(xn1,15,512);

end

px115=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和50

for i=1:m

X21(:,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法

end

px211=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:m

px25(:,i)=pyulear(xn2,15,512);

px26(:,i)=pburg(xn2,15,512);

end

px215=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和50

for i=1:m

X31(:,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法

end

pX311=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:m

px35(:,i)=pyulear(xn3,15,512);

px36(:,i)=pburg(xn3,15,512);

end

pX315=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 相关的操作。

建议：课上可以多教授一些关于MATLAB的相关函数的用法，以及常见的误区。