现代信号分析
3.2功率谱计算
由上式的输入输出关系可以看出这是一个AR模型，由已知数据得到：
H(z)?11?2.2377z?1?3.7476z?2?2.6293z?3?0.9224z?4
本题中的随机过程是由单位方差白噪声激励一个差分系统方程而产生。首先利用Matlab产生该信号源：
n=1:1:1024; wn=randn(1,length(n)); %产生随机单位白噪声序列 sn=filter([1],[1,-2.2377,3.7476,-2.6293,0.9224],wn); subplot; plot(sn,'r','LineWidth',1.4); title('s(n)','fontsize',12); 产生的信号源如下图3-1所示。对信号源信号求FFT变换，得到信号源的频谱图，作为估计算法的参考标准，如图3-2所示。
图3-1 Matlab产生的信号源 图3-2 信号源的频谱图
1）经典谱估计
下面首先采用直接法（周期图法）进行功率谱估计，程序如下： % nfft=128 window=boxcar(length(sn)); nfft=128; [Pxx,f]=periodogram(sn,window,nfft); subplot(211);plot(f,10*log10(Pxx)); xlabel('频率');ylabel('相对功率谱密度dB'); title('周期图法 点数128'); % nfft=1024 window=boxcar(length(sn)); % 矩形窗 nfft=1024; [Pxx,f]=periodogram(sn,window,nfft); subplot(212);;plot(f,10*log10(Pxx)); xlabel('频率'); ylabel('相对功率谱密度dB'); title('周期图法 点数1024'); 改变采样点数，获得的周期图法谱估计如下图所示： - 11 -
图3-3 不同序列长度周期法谱估计
运用周期图法通过对采样数据长度的调整，即在采样数据为128、1024时比较周期图法的效果，可知采集数据序列长度越大，谱分辨率越高，但谱曲线起伏也越大。因此需要改进，这里的改进指的是数据方差的改进。当前主要有两种改进方法：一种是周期图的平滑，即采用间接法估计功率谱，即BT法；另一种是平均法，它将长度为N的数据X(n)分成L段，分别求出每段的功率谱，然后加以平均，如Barlett法和Welch法。这里我们就采用以上3种方法对直接法进行改进。Matlab程序如下：
%BT法 nfft=1024; %以点数1024为例 cxn=xcorr(sn,'unbiased'); %求序列自相关函数 P1=fft(cxn,nfft); subplot(311); plot(f,10*log10(P1(1:257)),'r','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB ','fontsize',12); title('BT法','fontsize',12) % Bartlett法 window=boxcar(length(sn)); %矩形窗 noverlap=0; %数据无重叠 p=0.9; %置信概率 p2=psd(sn,nfft,1,window,noverlap,p); %计算PSD subplot(312);plot(f,10*log10(p2(1:257)),'r','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('Bartlett法','fontsize',12); %Welch法 window=boxcar(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(313);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('Welch法','fontsize',12); 各个经典功率谱估计方法的结果如图3-4所示。 - 12 -
现代信号分析
图3-4 经典谱估计方法的估计结果
为了了解使用不同种窗函数对功率谱的影响，以Welth法为例，进行了应用各种常见窗函数的仿真实验。程序如下:
%矩形窗 window=boxcar(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(221);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('矩形窗','fontsize',12); %汉明窗 window=hamming(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(222);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('汉明窗','fontsize',12); %布莱克曼窗 window=blackman(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(223);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('布莱克曼窗','fontsize',12); %三角窗 window=triang(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(224);plot(f,10*log10(p3),'b','LineWidth',1.4); - 13 -
xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('三角窗','fontsize',12); 估计结果如图3-5所示：
图3-5 不同窗函数对功率谱估计结果的影响
综合以上的波形，可以对经典功率谱估计作如下总结：经典功率谱估计，不论是周期图法还是BT法都可以用FFT快速计算，且计算量较小，物理概念明确，因而仍是目前较常用的谱估计方法。功率谱的分辨率正比于数据长度，当数据长度较大时，谱的曲线起伏加剧。BT法是对周期图法的改进，但其谱的平滑是以牺牲分辨率为代价的。Welth方法是对BT法德改进，它能使用更多的窗函数，并能更好的控制功率谱密度的估计的方差特性。在实际应用中，在方差、偏差和分辨率之间存在着矛盾，只能根据需要作折中处理。
2）现代谱估计
下面，根据AR模型法和burg算法对信号源进行现代功率谱估计。Matlab程序如下：
%------------现代功率谱分析---------% n=1:1:1024; wn=randn(1,length(n)); %产生随机单位白噪声序列 sn=filter([1],[1,-2.7377,3.7476,-2.6293,0.9224],wn); [Pxx,f]=periodogram(sn,window,nfft); nfft=1024; %AR模型法 [p4,w]=pyulear(sn,100,nfft); figure(3);subplot(2,1,1); plot(f,10*log10(abs(p4)),'b','LineWidth',1.4); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('AR模型法','fontsize',12); %burg算法 [p5,w]=pburg(sn,100,nfft); subplot(2,1,2); plot(f,10*log10(abs(p5)),'b','LineWidth',1.4); xlabel('频率','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('burg算法','fontsize',12); - 14 -
现代信号分析
估计结果如图3-6所示：
图3-6 现代功率谱估计方法的估计结果
为了分析现代功率谱中模型阶数对谱估计的影响，现利用AR模型法，改变模型阶数。分别用10、20、100作为模型阶数在采样点数分别为128和1024下进行分析，结果如下： %---------=-------现代功率谱分析---=---------% %AR模型法 n=1:1:1024; wn=randn(1,length(n)); %产生随机单位白噪声序列 sn=filter([1],[ 1,-2.7377,3.7476,-2.6293,0.9224],wn); [P1,w4]=pmcov(sn,10,128); [P2,w5]=pmcov(sn,50,128); [P3,w6]=pmcov(sn,100,128); [P4,w1]=pmcov(sn,10,1024); [P5,w2]=pmcov(sn,50,1024); [P6,w3]=pmcov(sn,100,1024); subplot(3,2,1); plot(w4/pi,10*log10(P1)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('128 10阶AR模型功率谱图');grid on; subplot(3,2,3); plot(w5/pi,10*log10(P2)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('128 50阶AR模型功率谱图');grid on; subplot(3,2,5); plot(w6/pi,10*log10(P3)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('128 100阶AR模型功率谱图');grid on; subplot(3,2,2); plot(w1/pi,10*log10(P4)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('1024 10阶AR模型功率谱图');grid on; subplot(3,2,4); plot(w2/pi,10*log10(P5)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('1024 50阶AR模型功率谱图');grid on; subplot(3,2,6); plot(w3/pi,10*log10(P6)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('1024 100阶AR模型功率谱图');grid on; - 15 -