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

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

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

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

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

5、數(shù)據(jù)的長(zhǎng)度和FFT所用的數(shù)據(jù)長(zhǎng)度nfft=1024;N=1024;%每段長(zhǎng)度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;%直接法求功率譜%將隨機(jī)信號(hào)x(n)的N點(diǎn)樣本值看作xN(n)能量有限信號(hào)，取其傅立葉變換，得到X ；然后再取其幅值的平方，并除以N作為x(n)%的真實(shí)功率譜P的估計(jì)。%計(jì)算序列的DFTXN=fft(xn,nfft);%對(duì)序列取絕對(duì)值后平方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;%間接法求功率譜%又稱(chēng)為自相關(guān)法或BT法，是由隨機(jī)信號(hào)N個(gè)觀察值估計(jì)出自相關(guān)函數(shù)R(m),然后再求R(m)傅立葉變換作為功率譜的估計(jì)PBT%計(jì)算序列的自相關(guān)函數(shù)RmRm=xcorr(xn,'unbiased'); %計(jì)算自相關(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點(diǎn)的有限長(zhǎng)序列x(n)分段，對(duì)各段用周期圖法求解功率后再平均。%加矩形窗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);%繪制功率譜曲線(xiàn)figure(3);plot(Fpbar,PBAR);xlabel('頻率/Hz');ylabel('功率譜/dB');title('bartlett法 4*256');grid;%welch法求功率譜%Welch法對(duì)Bartlett法進(jìn)行了兩方面的修正，一是選擇適當(dāng)?shù)拇昂瘮?shù)w(n)，并再周期圖計(jì)算前直接加進(jìn)去，加窗的優(yōu)點(diǎn)是無(wú)論什么樣的窗函數(shù)均可使譜估計(jì)非負(fù)。%二是在分段時(shí)，可使各段之間有重疊，這樣會(huì)使方差減小。%加漢寧窗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