經(jīng)典功率譜設(shè)計

1、DSP實驗報告實驗題目: 實驗 功率譜估計一、實驗要求：（1）理解功率譜估計的基本概念；（2）掌握經(jīng)典功率譜估計方法直接法和間接法；（3）掌握改進的經(jīng)典功率譜估計方法，例如Welch法。二、實驗內(nèi)容與原理：功率譜估計就是基于有限的數(shù)據(jù)尋找信號、隨機過程或系統(tǒng)的頻率成分。它表示隨機信號頻率域的統(tǒng)計特性。隨機信號是無始無終具有無限能量的，所以其傅立葉變換并不存在，因為它不滿足絕對可積的條件。因此需要研究其在頻率域上的功率分布情況，即功率譜密度或功率譜。根據(jù)實驗要求，完成該實驗首先要正確的生成被估計信號。數(shù)據(jù)的長度和FFT所用的數(shù)據(jù)長度都設(shè)為1024。1.周期圖法：直接法，即周期圖法，是由傅立葉變換

2、得到的：將隨機信號的N點樣本值看作能量有限信號，取其傅立葉變換，得到；然后再取其幅值的平方，并除以N作為的真實功率譜的估計，即實驗中，將隨機信號x(n)的N點樣本值看作xN(n)能量有限信號，取其傅立葉變換，得到X，然后再取其幅值的平方，并除以N作為x(n)的真實功率譜P的估計。2.間接法：間接法，又稱為自相關(guān)法或BT法，是由隨機信號N個觀察值，估計出自相關(guān)函數(shù)，然后再求的傅立葉變換作為功率譜的估計： 即如下計算： 實驗中由隨機信號N個觀察值估計出自相關(guān)函數(shù)R(m),然后再求R(m)傅立葉變換作為功率譜的估計PBT。直接法和間接法的方差性能很差，而且當(dāng)數(shù)據(jù)長度太大時，譜曲線起伏加劇；若數(shù)據(jù)長度

3、太小，則譜的分辨率又不好，所以需要改進。改進的直接譜估計方法由Bartlett法和Welch法。3.BARTLETT算法Bartlett法將采樣數(shù)據(jù)分成L段，每段的長度都是M，即N=LM，對每段數(shù)據(jù)加矩形窗，再計算其各自的功率譜，把對應(yīng)相加，再取平均，得到平均周期圖。即如下過程：首先將觀測數(shù)據(jù)分為L段，每段長M，分段后的數(shù)據(jù)可以用式(1)表示： ； (1)其中 平均的周期圖為： BARTLETT法是周期圖算法的一種改進。由概率論的知識知道，如果 是N個不相關(guān)的隨機變量，每個隨機變量的期望值為，方差為，那么將這N個隨機變量求平均，它的期望仍為，方差變?yōu)椤ARTLETT法即是受此啟發(fā)，將觀測數(shù)據(jù)

4、分段，先求每段數(shù)據(jù)的周期圖，再求平均的周期圖，當(dāng)分段較多時，估計出的功率譜較平滑，頻率分辨率較差；當(dāng)分段較少時，估計出的功率譜起伏較大，頻率分辨率較好。4.WELCH算法Welch法是對Bartlett法的改進：一，在對分段時，可允許每段數(shù)據(jù)有部分重疊；二，每段數(shù)據(jù)窗口可以不是矩形窗口，例如使用漢寧窗或哈明窗，記為。然后按Bartlett法求每一段的功率譜，記為=，其中。平均后的功率譜為：（段數(shù)） 加窗的優(yōu)點是使得無論對于什么樣的窗函數(shù)均可以譜估計為非負值；二是在分段時，各段之間有重疊，這樣會使方差減小。三、實驗?zāi)康模悍謩e用四種不同的發(fā)放進行功率譜估計，并對比結(jié)果。四、實驗程序：clear;%

5、數(shù)據(jù)的長度和FFT所用的數(shù)據(jù)長度nfft=1024;N=1024;%每段長度Ns=256;%產(chǎn)生含有噪聲的序列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;%并轉(zhuǎn)化為dBPERdb=10*log10(PER);%給出頻

6、率序列f=(0:length(PERdb)-1)/length(PERdb);%繪制功率譜圖形figure(1);plot(f,PERdb);xlabel('頻率/Hz');ylabel('功率譜/dB');title('直接法 N=1024');grid;%間接法求功率譜%又稱為自相關(guān)法或BT法，是由隨機信號N個觀察值估計出自相關(guān)函數(shù)R(m),然后再求R(m)傅立葉變換作為功率譜的估計PBT%計算序列的自相關(guān)函數(shù)RmRm=xcorr(xn,'unbiased'); %計算自相關(guān)函數(shù)Rm的DTFTPBT=fft(Rm,nfft);

7、%把PBT轉(zhuǎn)化為dBPBTdb=10*log10(abs(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);P

8、BAR1=abs(fft(window.*xn(1:256),Ns).2)/normvalue2;%第一段功率譜PBAR2=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振幅譜的平均值，并轉(zhuǎn)化為dBPBAR=10*log10(PBAR1+PBAR2+PBAR3+PBAR4)/4);%

9、給出頻率序列Fpbar=(0:length(PBAR)-1)/length(PBAR);%繪制功率譜曲線figure(3);plot(Fpbar,PBAR);xlabel('頻率/Hz');ylabel('功率譜/dB');title('bartlett法 4*256');grid;%welch法求功率譜%Welch法對Bartlett法進行了兩方面的修正，一是選擇適當(dāng)?shù)拇昂瘮?shù)w(n)，并再周期圖計算前直接加進去，加窗的優(yōu)點是無論什么樣的窗函數(shù)均可使譜估計非負。%二是在分段時，可使各段之間有重疊，這樣會使方差減小。%加漢寧窗hanningwind

10、ow=hanning(Ns)'normvalue=norm(window);PWel1=abs(fft(hanningwindow.*xn(1:256),Ns).2)/normvalue2;%第一段功率譜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振幅譜的平均值，并轉(zhuǎn)化為dBPWel=10*log10(PWel1+PWel2+PWel3+PWel4+PWel5+PWel6+PWel7)/7