古典法功率譜估計

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

2、)快速加窗（2M輸出矩形窗截斷相關(guān)法譜估計運算簡要框圖圖中快速相關(guān)的輸出時從- （N-1）到（ N-1）的 2N-1 點，加窗后截取的是 - （ M-1）到（ M-1）的，最后做（2M-1）點FFT，即可得到結(jié)果。（三）程序示例程序的主要思路就是按照運算框圖一步一步進行計算，下面附程序并進行簡要解釋：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的信號Xk=fft(xn,1024);%進行 2N-1 點 FFT,系統(tǒng)會自動補 0Sk

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

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

5、法譜估計與自相關(guān)法譜估計的差異只是估計自相關(guān)函數(shù)的方法不同。（二）運算簡要框圖矩形窗（長度N）截斷X(n)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的信號Xk1=fft(xn,512);%進行 N點 FFTSk3=abs(Xk1).*(abs(Xk1)./N;% 取頻譜幅度的平方，并除以 N，以此作為對 xn 真實功率譜的估計Sk4=10*log(Sk3);%對估計出的Sk 取對數(shù)，使畫出的圖更加突出

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

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

8、截斷X(n)分成 L 小段周 期圖輸出對 求（三）程序示例%L=2時 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的信號的前半段Xk1=fft(x1n,N);%進行 N點 FFTSk5=abs(Xk1).2./N;%取頻譜幅度的平方，并除以N，以此作為對 xn 真實功率譜的估計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的信號后半段Xk2=fft(x2n,N);%進行 N點 FFTSk6=abs(Xk2).2./N;%取頻譜幅度的平方，并除以N，以此作為對 xn 真實功率譜的估計Sk7=(Sk5+Sk6)/2;%相加求平均Sk8=10*log(Sk7);n=0:255;figure(4)subplot(2,1,1);plot(n/N,Sk7);title(d) Bartlett法估計功率譜密度 L=2)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk8);ylabel(10log(PSD)%L=4時 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的信號的 1/4 段 Xk1=fft(x1n,N);% 進行 N點 FFTSk1=abs(Xk1).2./N;% 取頻譜幅度的平方，并除以 N，以此作為對 xn 真實功率譜的估計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的信號的 1/4 段 Xk2=fft(x2n,N);% 進行 N點 FFTSk2=abs(Xk2).*(abs(Xk2)./N

11、;% 取頻譜幅度的平方，并除以 N，以此作為對 xn真實功率譜的估計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的信號的 1/4 段 Xk3=fft(x3n,N);% 進行 N點 FFTSk3=abs(Xk3).2./N;% 取頻譜幅度的平方，并除以 N，以此作為對 xn 真實功率譜的估計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的信號的 1/4 段 Xk4=fft(x4n,N);% 進行 N點 FFTSk4=abs(Xk4).*(abs(Xk4)./N;% 取頻譜幅度的平方，并除以 N，以此作為對 xn真實功率譜的估計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法估計功率譜密度L=4)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk6);ylabel(10log(PSD)（四）結(jié)果分析下面是程序運行后的結(jié)果(d)

13、 Bartlett 法 估 計 功 率 譜 密 度 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 法 估 計 功 率 譜 密 度 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 時用 bartlett法進行信號功率譜估計的波形。從上圖中我們同樣可以看到信號中有兩個頻

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

15、兩者的缺點。其主要步驟如下：第一步：將 N長的數(shù)據(jù)段分成L 個小段，每小段M點，相鄰小段見交疊M/2點。第二步：對個小段加同樣的品掛窗后求傅里葉變換第三步：求個小段功率譜的平均，得這里（二）運算簡要框圖矩形窗截斷X （ 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的信號一部分 wn=hanning(128);L=7;U=61.7942;%U 是 128 長窗函數(shù)的能量 , 那個函數(shù)不會寫， 這個數(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 點交疊） )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é)果分析下面是程序運行后的結(jié)果8060kS4020000.10.20.30.40.50.60.70.80.9160)40DSP(g20ol01