经典功率谱设计(共9页)

1、精选优质文档-倾情为你奉上DSP实验报告实验题目: 实验 功率谱估计一、实验要求：（1）理解功率谱估计的基本概念；（2）掌握经典功率谱估计方法直接法和间接法；（3）掌握改进的经典功率谱估计方法，例如Welch法。二、实验内容与原理：功率谱估计就是基于有限的数据寻找信号、随机过程或系统的频率成分。它表示随机信号频率域的统计特性。随机信号是无始无终具有无限能量的，所以其傅立叶变换并不存在，因为它不满足绝对可积的条件。因此需要研究其在频率域上的功率分布情况，即功率谱密度或功率谱。根据实验要求，完成该实验首先要正确的生成被估计信号。数据的长度和FFT所用的数据长度都设为1024。1.周期图法：直接法，

2、即周期图法，是由傅立叶变换得到的：将随机信号的N点样本值看作能量有限信号，取其傅立叶变换，得到；然后再取其幅值的平方，并除以N作为的真实功率谱的估计，即实验中，将随机信号x(n)的N点样本值看作xN(n)能量有限信号，取其傅立叶变换，得到X，然后再取其幅值的平方，并除以N作为x(n)的真实功率谱P的估计。2.间接法：间接法，又称为自相关法或BT法，是由随机信号N个观察值，估计出自相关函数，然后再求的傅立叶变换作为功率谱的估计： 即如下计算： 实验中由随机信号N个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT。直接法和间接法的方差性能很差，而且当数据长度太大时，

3、谱曲线起伏加剧；若数据长度太小，则谱的分辨率又不好，所以需要改进。改进的直接谱估计方法由Bartlett法和Welch法。3.BARTLETT算法Bartlett法将采样数据分成L段，每段的长度都是M，即N=LM，对每段数据加矩形窗，再计算其各自的功率谱，把对应相加，再取平均，得到平均周期图。即如下过程：首先将观测数据分为L段，每段长M，分段后的数据可以用式(1)表示： ； (1)其中 平均的周期图为： BARTLETT法是周期图算法的一种改进。由概率论的知识知道，如果 是N个不相关的随机变量，每个随机变量的期望值为，方差为，那么将这N个随机变量求平均，它的期望仍为，方差变为。BARTLETT

4、法即是受此启发，将观测数据分段，先求每段数据的周期图，再求平均的周期图，当分段较多时，估计出的功率谱较平滑，频率分辨率较差；当分段较少时，估计出的功率谱起伏较大，频率分辨率较好。4.WELCH算法Welch法是对Bartlett法的改进：一，在对分段时，可允许每段数据有部分重叠；二，每段数据窗口可以不是矩形窗口，例如使用汉宁窗或哈明窗，记为。然后按Bartlett法求每一段的功率谱，记为=，其中。平均后的功率谱为：（段数） 加窗的优点是使得无论对于什么样的窗函数均可以谱估计为非负值；二是在分段时，各段之间有重叠，这样会使方差减小。三、实验目的：分别用四种不同的发放进行功率谱估计，并对比结果。四

5、、实验程序：clear;%数据的长度和FFT所用的数据长度nfft=1024;N=1024;%每段长度Ns=256;%产生含有噪声的序列xnn=0:N-1;w1=2*pi*0.02;w2=2*pi*0.28;wn=randn(1,N);xn=sin(w1*n)+2*cos(w2*n)+wn;%直接法求功率谱%将随机信号x(n)的N点样本值看作xN(n)能量有限信号，取其傅立叶变换，得到X ；然后再取其幅值的平方，并除以N作为x(n)%的真实功率谱P的估计。%计算序列的DFTXN=fft(xn,nfft);%对序列取绝对值后平方PER=abs(XN).2/N;%并转化为dBPERdb=10*lo

6、g10(PER);%给出频率序列f=(0:length(PERdb)-1)/length(PERdb);%绘制功率谱图形figure(1);plot(f,PERdb);xlabel(频率/Hz);ylabel(功率谱/dB);title(直接法 N=1024);grid;%间接法求功率谱%又称为自相关法或BT法，是由随机信号N个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT%计算序列的自相关函数RmRm=xcorr(xn,unbiased); %计算自相关函数Rm的DTFTPBT=fft(Rm,nfft);%把PBT转化为dBPBTdb=10*log10(ab

7、s(PBT);Fbt=(0:length(PBTdb)-1)/length(PBTdb);%绘制功率谱图形figure(2);plot(Fbt,PBTdb);xlabel(频率/Hz);ylabel(功率谱/dB);title(间接法 N=1024);grid;%bartlett法求功率谱%Bartlett平均周期图的方法是将N点的有限长序列x(n)分段，对各段用周期图法求解功率后再平均。%加矩形窗window=ones(1,Ns);normvalue=norm(window);PBAR1=abs(fft(window.*xn(1:256),Ns).2)/normvalue2;%第一段功率谱P

8、BAR2=abs(fft(window.*xn(257:512),Ns).2)/normvalue2;%第二段功率谱PBAR3=abs(fft(window.*xn(513:768),Ns).2)/normvalue2;%第三段功率谱PBAR4=abs(fft(window.*xn(769:1024),Ns).2)/normvalue2;%第四段功率谱%求Fourier振幅谱的平均值，并转化为dBPBAR=10*log10(PBAR1+PBAR2+PBAR3+PBAR4)/4);%给出频率序列Fpbar=(0:length(PBAR)-1)/length(PBAR);%绘制功率谱曲线figur

9、e(3);plot(Fpbar,PBAR);xlabel(频率/Hz);ylabel(功率谱/dB);title(bartlett法 4*256);grid;%welch法求功率谱%Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。%二是在分段时，可使各段之间有重叠，这样会使方差减小。%加汉宁窗hanningwindow=hanning(Ns);normvalue=norm(window);PWel1=abs(fft(hanningwindow.*xn(1:256),Ns).2)/nor

10、mvalue2;%第一段功率谱PWel2=abs(fft(hanningwindow.*xn(129:384),Ns).2)/normvalue2;%第二段功率谱PWel3=abs(fft(hanningwindow.*xn(257:512),Ns).2)/normvalue2;%第三段功率谱PWel4=abs(fft(hanningwindow.*xn(385:640),Ns).2)/normvalue2;%第四段功率谱PWel5=abs(fft(hanningwindow.*xn(513:768),Ns).2)/normvalue2;%第五段功率谱PWel6=abs(fft(hanningwindow.*xn(641:896),Ns).2)/normvalue2;%第六段功率谱PWel7=abs(fft(hanningwindow.*xn(769:1024),Ns).2)/normvalue2;%第七段功率谱%求Fourier振幅谱的平均值，并转化为dBPWel=10*log10(PWel1+PWel2+PWel3+PWel4+PWe