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

1、word古典法功率譜估計(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)波形二、相關(guān)法功率譜估計(jì)一算法原理簡(jiǎn)介此方法以相關(guān)函數(shù)為媒介來(lái)計(jì)算功率譜，所以又叫間接法。它是1958年由Blackman 和Tu

2、key提出。這種方法的具體步驟是： 第一步：從無(wú)限長(zhǎng)隨機(jī)序列x(n)中截取長(zhǎng)度N的有限長(zhǎng)序列 第二步：由N長(zhǎng)序列求2M-1點(diǎn)的自相關(guān)函數(shù)序列。 第三步：由相關(guān)函數(shù)的傅式變換求功率譜。 以上過(guò)程中經(jīng)歷了兩次截?cái)?，一次是將x(n)截成N長(zhǎng)，稱為加數(shù)據(jù)窗，一次是將想x(n)截成2M-1長(zhǎng)，稱為加延遲窗。因此所得的功率譜僅是近似值，也叫譜估計(jì)。一般取M<<N，因?yàn)橹挥挟?dāng)M較小時(shí)，序列傅式變換的點(diǎn)數(shù)才較小，功率譜的計(jì)算量才不至于大到難以實(shí)現(xiàn)，而且譜估計(jì)質(zhì)量也較好。因此，在FFT問(wèn)世之前，相關(guān)法是最常用的譜估計(jì)方法。當(dāng)FFT問(wèn)世后，情況有所變化。因?yàn)榻財(cái)嗪蟮?(nxN可視作能量信號(hào)，由相關(guān)卷積

3、定理可得這就將相關(guān)化為線性卷積，而線性卷積又可以用快速卷積來(lái)實(shí)現(xiàn)。我們可對(duì)上式兩邊取2N-1點(diǎn)DFT，那么有于是將時(shí)域卷積變?yōu)轭l域乘積，用快速相關(guān)求自相關(guān)函數(shù)估值的完整方案如下： 1. 對(duì)N長(zhǎng))(nxN的充N-1個(gè)零，成為2N-1長(zhǎng)的。 2. 求2N-1點(diǎn)的FFT，得3、求 4、 求2N-1點(diǎn)的IFFT二運(yùn)算簡(jiǎn)要框圖X(n)快速相關(guān)加窗截?cái)?M-1點(diǎn)FFT輸出矩形窗截?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)要

4、解釋：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=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)法

5、估計(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é)果從上圖中我們可以較為明顯的看到信號(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ù)

6、構(gòu)成有限長(zhǎng)序列直接求傅里葉變換，得頻譜。第二步：取頻譜幅度的平方，并除以N，以此作為對(duì)x(n)真實(shí)功率譜的估計(jì)。事實(shí)上，周期圖法譜估計(jì)與自相關(guān)法譜估計(jì)的差異只是估計(jì)自相關(guān)函數(shù)的方法不同。二運(yùn)算簡(jiǎn)要框圖 矩形窗長(zhǎng)度N截?cái)郚點(diǎn)FFTX(n) 圖中用FFT來(lái)代替傅里葉變換三程序例如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ì)

7、xn真實(shí)功率譜的估計(jì)Sk4=10*log(Sk3); %對(duì)估計(jì)出的Sk取對(duì)數(shù)，使畫出的圖更加突出特點(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é)果從上圖中我們同樣可以看到信號(hào)中有兩個(gè)頻率分量，一個(gè)在0.2處，一個(gè)在0.4處，與產(chǎn)生的信號(hào)相一致。但是我們不難看出，估計(jì)出的功率譜譜線與相關(guān)法功率譜估計(jì)一樣非常不平坦，有很多

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

9、時(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的AWGN的信號(hào)后半段Xk2=fft(x2n,N); %進(jìn)行N點(diǎn)FFTSk6=abs(X

10、k2).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.2*n)+ cos(2*p

11、i*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;%取頻譜幅度的平方，并除以N，以此作為對(duì)xn真實(shí)功率譜的估計(jì)N=128

12、,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的AWGN的信號(hào)的1/4段Xk4=fft(x4n,N);%進(jìn)行N點(diǎn)FFTSk4=abs(Xk4).*(abs(Xk

13、4)./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é)果上圖分別是L=2和L=4時(shí)用bartlett法進(jìn)行信號(hào)功率譜估計(jì)的波形。從

14、上圖中我們同樣可以看到信號(hào)中有兩個(gè)頻率分量，一個(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ì)改良方法是Welch法，又叫加權(quán)交疊平均法。這種方法以加窗加權(quán)求取平滑，以分段重疊求得平均，因此集平均與平滑的優(yōu)點(diǎn)于一體，同時(shí)也不可防止的帶有

15、兩者的缺點(diǎn)。其主要步驟如下：第一步：將N長(zhǎng)的數(shù)據(jù)段分成L個(gè)小段，每小段M點(diǎn)，相鄰小段見(jiàn)交疊M/2點(diǎn)。第二步：對(duì)個(gè)小段加同樣的品掛窗后求傅里葉變換第三步：求個(gè)小段功率譜的平均，得這里 二運(yùn)算簡(jiǎn)要框圖矩形窗截?cái)嘞嗉忧笃骄化B進(jìn)行功率譜估計(jì) Xn分成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ù)是自己加出來(lái)的。x11n=x1n.*

16、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);

17、N=128,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

18、*pi*0.4*n)+2*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)