古典法功率譜估計(jì)

1、古典法功率譜估計(jì)古典法功率譜估計(jì)一、 信號(hào)的產(chǎn)生（一）信號(hào)組成在本實(shí)驗(yàn)中，需要事先產(chǎn)生待估計(jì)的信號(hào)，為了使實(shí)驗(yàn)結(jié)果較為明顯，我產(chǎn)生了由兩個(gè)不同頻率的正弦信號(hào) （頻率差相對(duì)較大） 和加性高斯白噪聲組成的信號(hào)。（二）程序xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為0，方差為1 的AWGN信號(hào)figure(1)plot(n,xn);title(a)兩個(gè)正弦信號(hào)與白噪聲疊加的時(shí)域波形)（三）信號(hào)波形兩個(gè)正弦信號(hào)與白噪聲疊加的時(shí)域波形1086420-2-4-6-80100200300400500600（二）運(yùn)算簡(jiǎn)要框圖X(n

2、)快速加窗（2M輸出矩形窗截?cái)嘞嚓P(guān)法譜估計(jì)運(yùn)算簡(jiǎn)要框圖圖中快速相關(guān)的輸出時(shí)從- （N-1）到（ N-1）的 2N-1 點(diǎn)，加窗后截取的是 - （ M-1）到（ M-1）的，最后做（2M-1）點(diǎn)FFT，即可得到結(jié)果。（三）程序示例程序的主要思路就是按照運(yùn)算框圖一步一步進(jìn)行計(jì)算，下面附程序并進(jìn)行簡(jiǎn)要解釋：N=512,n=0:N-1; %N 是 FFT的變換區(qū)間 xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為0，方差為 1 的 AWGN的信號(hào)Xk=fft(xn,1024);%進(jìn)行 2N-1 點(diǎn) FFT,系統(tǒng)會(huì)自動(dòng)補(bǔ) 0Sk

3、=abs(Xk).*(abs(Xk)./N; %取頻譜幅度的平方，并除以N，以此作為對(duì)xn 真實(shí)功率譜的估計(jì)Rn=ifft(Sk);Sk1=fft(Rn,512);figure(2)subplot(2,1,1);plot(n/N,Sk1);ylabel(Sk)title(b)相關(guān)法估計(jì)功率譜密度 )Sk2=10*log(Sk1);%對(duì)估計(jì)出的 Sk 取對(duì)數(shù)，使畫出的圖更加突出特點(diǎn)subplot(2,1,2);plot(n/N,Sk2);ylabel(10log(PSD)（四）結(jié)果分析下面是程序運(yùn)行后的結(jié)果相關(guān)法估計(jì)功率譜密度200150kS10050000.10.20.30.40.50.60.

4、70.80.9160)40P(gol0201000.10.20.30.40.50.60.70.80.91從上圖中我們可以較為明顯的看到信號(hào)中有兩個(gè)頻率分量，一個(gè)在 0.2 處，一個(gè)在 0.4 處，與產(chǎn)生的信號(hào)相一致。但是我們不難看出，估計(jì)出的功率譜譜線非常不平坦，有很多起伏。三、 周期圖法譜估計(jì)（一）算法原理簡(jiǎn)介周期圖法又稱直接法。它是從隨機(jī)信號(hào)x(n) 中截取 N 長(zhǎng)的一段，把它視為能量有限 x(n) 真實(shí)功率譜的估計(jì)的抽樣。其具體步驟如下：第一步：由獲得的N點(diǎn)數(shù)據(jù)構(gòu)成有限長(zhǎng)序列直接求傅里葉變換，得頻譜。第二步：取頻譜幅度的平方，并除以N，以此作為對(duì)x(n) 真實(shí)功率譜的估計(jì)。事實(shí)上，周期圖

5、法譜估計(jì)與自相關(guān)法譜估計(jì)的差異只是估計(jì)自相關(guān)函數(shù)的方法不同。（二）運(yùn)算簡(jiǎn)要框圖矩形窗（長(zhǎng)度N）截?cái)郮(n)N 點(diǎn)FFT圖中用 FFT來代替傅里葉變換（三）程序示例N=512,n=0:N-1;xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的 AWGN的信號(hào)Xk1=fft(xn,512);%進(jìn)行 N點(diǎn) FFTSk3=abs(Xk1).*(abs(Xk1)./N;% 取頻譜幅度的平方，并除以 N，以此作為對(duì) xn 真實(shí)功率譜的估計(jì)Sk4=10*log(Sk3);%對(duì)估計(jì)出的Sk 取對(duì)數(shù)，使畫出的圖更加突出

6、特點(diǎn)figure(3)subplot(2,1,1);plot(n/N,Sk3);title(c)周期圖法估計(jì)功率譜密度)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk4);ylabel(10log(PSD)（四）結(jié)果分析下面是程序運(yùn)行后的結(jié)果周期圖法估計(jì)功率譜密度300200kS100000.10.20.30.40.50.60.70.80.91100)50DSP(0go01-50-10000.10.20.30.40.50.60.70.80.91從上圖中我們同樣可以看到信號(hào)中有兩個(gè)頻率分量， 一個(gè)在 0.2 處，一個(gè)在 0.4 處，與產(chǎn)生的信號(hào)相一致。但是我們不難看出，

7、估計(jì)出的功率譜譜線與相關(guān)法功率譜估計(jì)一樣非常不平坦，有很多起伏。四、 Bartlett法功率譜估計(jì)（一）算法原理簡(jiǎn)介當(dāng)我們用相關(guān)法或者周期圖法對(duì)信號(hào)的功率譜進(jìn)行估計(jì)時(shí)，都不是對(duì)的一致估計(jì)，主要原因是方差大。于是就產(chǎn)生了周期圖法的改進(jìn)。改進(jìn)的主要途徑是平滑和平均。平滑是用一個(gè)適當(dāng)?shù)拇昂瘮?shù)與計(jì)算的功率譜進(jìn)行卷積，是譜線平滑。這種方法的出的譜估計(jì)是無偏的，方差也小，但分辨率下降。平均就是將截取的數(shù)據(jù)段再分成 L 個(gè)小段，分別計(jì)算功率譜后取功率譜的平均。因?yàn)?L 個(gè)平均的方差比隨機(jī)變量的單獨(dú)方差小 L 倍，所以當(dāng) L 趨于無窮時(shí)，L 個(gè)平均的方差趨于零，可以達(dá)到一致估計(jì)的目的。（二）運(yùn)算簡(jiǎn)要框圖矩形窗

8、截?cái)郮(n)分成 L 小段周 期圖輸出對(duì) 求（三）程序示例%L=2時(shí) bartlett法N=256,n=0:255;x1n=2*cos(2*pi*0.2*n)+cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的 AWGN的信號(hào)的前半段Xk1=fft(x1n,N);%進(jìn)行 N點(diǎn) FFTSk5=abs(Xk1).2./N;%取頻譜幅度的平方，并除以N，以此作為對(duì) xn 真實(shí)功率譜的估計(jì)n=256:511;x2n=2*cos(2*pi*0.2*n)+cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的

9、AWGN的信號(hào)后半段Xk2=fft(x2n,N);%進(jìn)行 N點(diǎn) FFTSk6=abs(Xk2).2./N;%取頻譜幅度的平方，并除以N，以此作為對(duì) xn 真實(shí)功率譜的估計(jì)Sk7=(Sk5+Sk6)/2;%相加求平均Sk8=10*log(Sk7);n=0:255;figure(4)subplot(2,1,1);plot(n/N,Sk7);title(d) Bartlett法估計(jì)功率譜密度 L=2)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk8);ylabel(10log(PSD)%L=4時(shí) bartlett法N=128,n=0:127;x1n=2*cos(2*pi*0

10、.2*n)+cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的 AWGN的信號(hào)的 1/4 段 Xk1=fft(x1n,N);% 進(jìn)行 N點(diǎn) FFTSk1=abs(Xk1).2./N;% 取頻譜幅度的平方，并除以 N，以此作為對(duì) xn 真實(shí)功率譜的估計(jì)n=128:255;x2n=2*cos(2*pi*0.2*n)+cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的 AWGN的信號(hào)的 1/4 段 Xk2=fft(x2n,N);% 進(jìn)行 N點(diǎn) FFTSk2=abs(Xk2).*(abs(Xk2)./N

11、;% 取頻譜幅度的平方，并除以 N，以此作為對(duì) xn真實(shí)功率譜的估計(jì)N=128,n=256:383;x3n=2*cos(2*pi*0.2*n)+cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的 AWGN的信號(hào)的 1/4 段 Xk3=fft(x3n,N);% 進(jìn)行 N點(diǎn) FFTSk3=abs(Xk3).2./N;% 取頻譜幅度的平方，并除以 N，以此作為對(duì) xn 真實(shí)功率譜的估計(jì)n=384:511;x4n=2*cos(2*pi*0.2*n)+cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的 A

12、WGN的信號(hào)的 1/4 段 Xk4=fft(x4n,N);% 進(jìn)行 N點(diǎn) FFTSk4=abs(Xk4).*(abs(Xk4)./N;% 取頻譜幅度的平方，并除以 N，以此作為對(duì) xn真實(shí)功率譜的估計(jì)Sk5=(Sk1+Sk2+Sk3+Sk4)/4; % 相加求平均Sk6=10*log(Sk5);n=0:127;figure(5)subplot(2,1,1);plot(n/N,Sk5);title(e) Bartlett法估計(jì)功率譜密度L=4)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk6);ylabel(10log(PSD)（四）結(jié)果分析下面是程序運(yùn)行后的結(jié)果(d)

13、 Bartlett 法 估 計(jì) 功 率 譜 密 度 L=2300200kS100000.10.20.30.40.50.60.70.80.9160)40DSP(g20ol010-2000.10.20.30.40.50.60.70.80.91kS)DSP(gol01(e) Bartlett 法 估 計(jì) 功 率 譜 密 度 L=480604020000.10.20.30.40.50.60.70.80.916040200-2000.10.20.30.40.50.60.70.80.91上圖分別是 L=2 和 L=4 時(shí)用 bartlett法進(jìn)行信號(hào)功率譜估計(jì)的波形。從上圖中我們同樣可以看到信號(hào)中有兩個(gè)頻

14、率分量，一個(gè)在0.2 處，一個(gè)在 0.4 處，與產(chǎn)生的信號(hào)相一致。但是我們不難看出，估計(jì)出的功率譜譜線與之前相關(guān)法功率譜估計(jì)和周期圖法功率譜估計(jì)相比，波形相對(duì)平坦了一些。 L=2 和 L=4 時(shí)用 bartlett 法進(jìn)行信號(hào)功率譜估計(jì)的波形相比我們可以很明顯的看出L 大的那個(gè)波形更加平坦，這與之前在算法原理中介紹的一樣， L 越大，平均后的方差就越小，越能達(dá)到一致估計(jì)的目的。五、Welch 法功率譜估計(jì)（一）算法原理簡(jiǎn)介現(xiàn)在比較常用的功率譜估計(jì)改進(jìn)方法是 Welch 法，又叫加權(quán)交疊平均法。這種方法以加窗（加權(quán)）求取平滑，以分段重疊求得平均，因此集平均與平滑的優(yōu)點(diǎn)于一體，同時(shí)也不可避免的帶有

15、兩者的缺點(diǎn)。其主要步驟如下：第一步：將 N長(zhǎng)的數(shù)據(jù)段分成L 個(gè)小段，每小段M點(diǎn)，相鄰小段見交疊M/2點(diǎn)。第二步：對(duì)個(gè)小段加同樣的品掛窗后求傅里葉變換第三步：求個(gè)小段功率譜的平均，得這里（二）運(yùn)算簡(jiǎn)要框圖矩形窗截?cái)郮 （ n）分成L小段輸出交疊相（三）程序示例N=128;n=0:127;x1n=2*cos(2*pi*0.2*n)+cos(2*pi*0.4*n)+2*randn(size(n);%產(chǎn)生加有均值為 0，方差為 1 的 AWGN的信號(hào)一部分 wn=hanning(128);L=7;U=61.7942;%U 是 128 長(zhǎng)窗函數(shù)的能量 , 那個(gè)函數(shù)不會(huì)寫， 這個(gè)數(shù)是自己加出來的。x11n

16、=x1n.*wn;Xk1=fft(x11n,128);Sk1=abs(Xk1).2.*1/U;N=128;n=64:191;x2n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x22n=x2n.*wn;Xk2=fft(x22n,128);Sk2=abs(Xk2).2.*(1/U);N=128,n=128:255;x3n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n); x33n=x3n.*wn;Xk3=fft(x33n,128);Sk3=abs(Xk3).2.*(1/U);N=128

17、,n=192:319;x4n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x44n=x4n.*wn;Xk4=fft(x44n,128);Sk4=abs(Xk4).2.*(1/U);N=128,n=256:383;x5n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n); x55n=x5n.*wn;Xk5=fft(x55n,128);Sk5=abs(Xk5).2.*(1/U);N=128,n=320:447;x6n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2

18、*randn(size(n);x66n=x6n.*wn;Xk6=fft(x66n,128);Sk6=abs(Xk6).2.*(1/U);N=128,n=384:511;x7n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x77n=x7n.*wn;Xk7=fft(x77n,128);Sk7=abs(Xk7).2.*(1/U);n=0:127;Sk=(Sk1+Sk2+Sk3+Sk4+Sk5+Sk6+Sk7)./L;figure(6)title(f) Welch法（漢寧窗， L=7,64 點(diǎn)交疊） )subplot(2,1,1);plot(n/N,Sk);ylabel(Sk)Skk=10*log(Sk);subplot(2,1,2);plot(n/N,Skk);ylabel(10log(PSD)（四）結(jié)果分析下面是程序運(yùn)行后的結(jié)果8060kS4020000.10.20.30.40.50.60.70.80.9160)40DSP(g20ol01