经典功率谱估计及Matlab仿真
?原文出自：http://www.cnblogs.com/jacklu/p/5140913.html
功率譜估計(jì)在分析平穩(wěn)各態(tài)遍歷隨機(jī)信號(hào)頻率成分領(lǐng)域被廣泛使用，并且已被成功應(yīng)用到雷達(dá)信號(hào)處理、故障診斷等實(shí)際工程中。本文給出了經(jīng)典功率譜估計(jì)的幾類方法，并通過Matlab的實(shí)驗(yàn)仿真對(duì)經(jīng)典功率譜估計(jì)方法性能進(jìn)行了分析，最后說明了經(jīng)典功率譜估計(jì)法的局限性和造成這種局限性的原因。
1.引言
給定一個(gè)標(biāo)準(zhǔn)的正弦信號(hào)，我們可以通過傅里葉變換來(lái)分析它的頻率成分。然而，實(shí)際工程應(yīng)用中，由于存在著各種干擾、噪聲，我們得到的信號(hào)往往不是理想的，如圖1-1這種信號(hào)，具有不確定性，幅度不能預(yù)知，非周期，但往往服從一定的統(tǒng)計(jì)特性，這種信號(hào)叫作隨機(jī)信號(hào)。需要注意的是，本文所說的隨機(jī)信號(hào)是指平穩(wěn)各態(tài)遍歷的隨機(jī)信號(hào)，關(guān)于非平穩(wěn)隨機(jī)信號(hào)的分析方法[1]本文不予討論。
圖1-1 一種隨機(jī)信號(hào)時(shí)域形式
?
對(duì)于圖1-1的隨機(jī)信號(hào)，我們可以通過功率譜來(lái)分析它的頻率成分，如圖1-2所示為圖1-1隨機(jī)信號(hào)的功率譜。實(shí)際過程中，我們只能獲得隨機(jī)信號(hào)的一些離散數(shù)據(jù)點(diǎn)（假設(shè)為N個(gè)），本文將討論如何利用這N個(gè)數(shù)據(jù)點(diǎn)，來(lái)得到一個(gè)"非精確"的功率譜來(lái)對(duì)真實(shí)隨機(jī)信號(hào)的功率譜進(jìn)行估計(jì)，并討論如何更好的估計(jì)，即在下一章要講述的幾個(gè)經(jīng)典的功率譜估計(jì)法。
圖1-2 上圖所示的隨機(jī)信號(hào)功率譜
2.經(jīng)典功率譜估計(jì)法
上一章我們已經(jīng)知道功率譜估計(jì)法是通過利用已經(jīng)獲得的N個(gè)數(shù)據(jù)點(diǎn)，來(lái)得到一個(gè)"非精確"的功率譜對(duì)真實(shí)隨機(jī)信號(hào)的功率譜進(jìn)行估計(jì)，所以在給出具體的方法之前，如何來(lái)評(píng)價(jià)我們得出的這個(gè)"非精確"的功率譜的好壞呢？
評(píng)價(jià)功率譜性能好壞的標(biāo)準(zhǔn)有很多，本文只給出兩個(gè)影響最大的標(biāo)準(zhǔn)：分辨率和方差。分辨率即功率譜上能夠區(qū)分的最小相鄰頻率成分，分辨率越高，我們觀察信號(hào)的頻率成分越清晰；方差大小則反映到功率譜波動(dòng)性的大小，如果方差太大，功率譜波動(dòng)性大，則很容易造成有用的頻率成分被噪聲淹沒。所以，我們希望得到的這個(gè)"非精確"的功率譜，分辨率越高越好，方差越小越好。
同時(shí)，我們給出概率論與數(shù)理統(tǒng)計(jì)中所學(xué)的一致估計(jì)和非一致估計(jì)的概念，假定真實(shí)信號(hào)的功率譜為，估計(jì)得到的"非精確"功率譜為，如果滿足公式（2-1），則稱為一致估計(jì)。
???????????????????????????????????????? （2-1）
2.1周期圖法
2.1.1周期圖法原理
我們已知N個(gè)離散的數(shù)據(jù)點(diǎn)，對(duì)這些數(shù)據(jù)點(diǎn)進(jìn)行傅里葉變換，得到式（2-2）： (2-2)
再對(duì)（2-2）式取模的平方，除以N，即得到一個(gè)"非精確"的譜，如式（2-3）所示，這就是周期圖法的原理。
（2-3）
2.1.2周期圖法性能（Matlab仿真）
上一小節(jié)我們已經(jīng)給出了周期圖法的原理。本節(jié)將通過Matlab仿真給出數(shù)據(jù)點(diǎn)數(shù)N對(duì)功率譜性能好壞的影響，正如上文所述，將通過對(duì)所得功率譜的分辨率和方差兩方面進(jìn)行分析。
我們?cè)贛atlab中通過三個(gè)正弦函數(shù)和白噪聲疊加，構(gòu)造了一個(gè)隨機(jī)信號(hào)。其數(shù)學(xué)形式如式（2-4），其中頻率、、，幅值、、，相位為相互獨(dú)立在上服從均勻分布隨機(jī)相位，為均值為0，方差為1的實(shí)值高斯白噪聲，采樣頻率為1000。信號(hào)的時(shí)域形式如圖2-1所示。
（2-4）
圖2-1 實(shí)驗(yàn)所用的隨機(jī)信號(hào)
?
當(dāng)數(shù)據(jù)點(diǎn)數(shù)N分別為128、256、512和1024時(shí)，得到的功率譜分別如圖2-2、圖2-3、圖2-4和圖2-5所示。分辨率能夠直觀的通過功率譜圖形看出，方差的數(shù)值由表2-1給出。
?
圖2-2 N=128時(shí)周期圖法得到的功率譜
?
圖2-3 N=256時(shí)周期圖法得到的功率譜
?
圖2-4 N=512時(shí)周期圖法得到的功率譜
?
圖2-5 N=1024時(shí)周期圖法得到的功率譜
?
表2-1 不同N值得到功率譜的方差值
| N | 128 | 256 | 512 | 1024 |
| 方差 | 92.7108 | 130.9109 | 160.9187 | 483.5894 |
?
通過上面實(shí)驗(yàn)結(jié)果的比較，我們很容易發(fā)現(xiàn)，周期圖法得到的功率譜隨著數(shù)據(jù)點(diǎn)數(shù)N的增大，分辨率變大、方差變也大。
2.1.3平均周期圖法
周期圖法得到的功率譜與我們所期望的"分辨率大、方差小"是矛盾的。為了進(jìn)一步降低方差，將N個(gè)觀測(cè)樣本數(shù)據(jù)點(diǎn)分為L(zhǎng)段，每段數(shù)據(jù)長(zhǎng)度為M, 分別對(duì)每段數(shù)據(jù)求周期圖功率譜估計(jì)，然后求平均值，這種方法稱平均周期圖法。
那么這種方法會(huì)如何改善方差呢？我們給出證明： （2-5）
其中：
（2-6）
由式（2-5）我們可以看出，平均周期圖法將原來(lái)的方差變?yōu)樵瓉?lái)的，L為分段數(shù)。
2.1.4平均周期圖法性能（Matlab仿真）
當(dāng)數(shù)據(jù)點(diǎn)數(shù)N為1024，分段數(shù)分別為8、4、2時(shí)，平均周期圖法得到的功率譜分別如圖2-6、圖2-7、圖2-8所示。分辨率能夠直觀的通過功率譜圖形看出，方差的數(shù)值由表2-2給出。
?
圖2-6 L=8時(shí)平均周期圖法得到的功率譜
?
圖2-7 L=4時(shí)平均周期圖法得到的功率譜
?
圖2-8 L=2時(shí)平均周期圖法得到的功率譜
?
表2-2 不同L值得到功率譜的方差值
| L | 8 | 4 | 2 | 1 |
| 方差 | 96.3756 | 190.9647 | 400.6464 | 483.5894 |
?
L=1時(shí)，平均周期圖法退化為周期圖法。通過上面實(shí)驗(yàn)結(jié)果的比較，我們很容易發(fā)現(xiàn)，平均周期圖法得到的功率譜隨著分段數(shù)L變大，方差變小，但分辨率變小。
當(dāng)觀測(cè)樣本序列數(shù)據(jù)個(gè)數(shù)N固定時(shí)，要降低方差需要增加分段數(shù)L。當(dāng)N不大時(shí)分段長(zhǎng)度M取值較小，則功率譜分辨率降低到較低的水平。若分段數(shù)L固定時(shí)，增加分辨率需要分段長(zhǎng)度M，則需要采集到更長(zhǎng)的檢測(cè)數(shù)據(jù)序列。實(shí)際中恰恰是檢測(cè)樣本序列長(zhǎng)度不足。
2.1.5修正的平均周期圖法
上一節(jié)已經(jīng)提到實(shí)際中檢測(cè)樣本序列長(zhǎng)度是有限的。對(duì)現(xiàn)有數(shù)據(jù)長(zhǎng)度N，如果能獲得更多的段數(shù)分割，將會(huì)得到更小的方差。允許數(shù)據(jù)段間有重疊部分，來(lái)得到更多的段數(shù)。對(duì)段間重疊長(zhǎng)度的選取，最簡(jiǎn)單是取為段長(zhǎng)度M的一半。由式（2-5）可知更多的段數(shù)可以進(jìn)一步降低方差。
數(shù)據(jù)截?cái)嗟倪^程中相當(dāng)于數(shù)據(jù)加矩形窗，矩形窗幅度較大的旁瓣會(huì)造成"頻譜泄漏"。我們分段時(shí)采取的窗函數(shù)更為多樣(三角窗,海明窗等), 以減小截?cái)鄶?shù)據(jù)（加矩形窗）窗函數(shù)帶來(lái)的影響[2]
2.1.6修正的平均周期圖法性能（Matlab仿真）
利用修正平均周期圖法，分別使用矩形窗、Blackman窗和Hamming窗得到的功率譜如圖2-9所示。
?
圖2-9 不同窗函數(shù)的修正平均周期圖法得到的功率譜
?
可以發(fā)現(xiàn)，矩形窗的分辨率最高，但是方差也最大，這是由于矩形窗頻譜主瓣最窄，分辨率因此最高，旁瓣也高，導(dǎo)致頻譜泄漏最嚴(yán)重，方差最大。
2.1.7總結(jié)
周期圖法獲得的功率譜隨著樣本點(diǎn)數(shù)越多，分辨率越大、方差越大；平均周期圖法以犧牲分辨率來(lái)進(jìn)一步改善方差；修正的平均周期圖法允許段的重疊來(lái)進(jìn)一步增大分段數(shù)、或者分段數(shù)相同，每段樣本點(diǎn)數(shù)變多。無(wú)論是哪種方法都沒有徹底結(jié)局方差與分辨率之間的矛盾。
2.2相關(guān)功率譜估計(jì)法-BT法
正如我們之前介紹的，要提高功率譜估計(jì)的分辨率，必須增加數(shù)據(jù)序列的長(zhǎng)度N，但是較長(zhǎng)的數(shù)據(jù)序列，由噪聲引起的隨機(jī)性得到更為充分的體現(xiàn)－較大的方差。事實(shí)上，當(dāng)N無(wú)窮大時(shí)，方差為一非零常數(shù)。即周期圖法無(wú)法實(shí)現(xiàn)功率譜的一致估計(jì)。而這節(jié)講述的相關(guān)功率譜估計(jì)法（下文稱作BT法），是一致估計(jì)。
2.2.1 BT法的原理
維納辛欽定理指出，隨機(jī)信號(hào)的相關(guān)函數(shù)與它的功率譜是一對(duì)傅里葉變換對(duì)。BT法就是基于這個(gè)原理。先由觀測(cè)數(shù)據(jù)估計(jì)出自相關(guān)函數(shù)，然后求自相關(guān)函數(shù)的傅立葉變換，以此變換作為對(duì)功率譜的估計(jì)，也稱為間接法。BT法要求信號(hào)長(zhǎng)度N以外的信號(hào)為零，這也造成BT法的局限性。
的自相關(guān)函數(shù)定義如式（2-7）所示，得到的功率譜記為，則BT法可以表述為式（2-8）。
（2-7）
（2-8）
2.2.2 BT法的性能（Matlab仿真）
數(shù)據(jù)點(diǎn)數(shù)N分別為128、256、512和1024的BT法，得到的功率譜如圖2-10、圖2-11、圖2-12和圖2-13所示。
?
圖2-10 N=128時(shí)，BT法得到的功率譜
?
圖2-11 N=256時(shí)，BT法得到的功率譜
?
圖2-12 N=512時(shí)，BT法得到的功率譜
?
圖2-13 N=1024時(shí)，BT法得到的功率譜
?
由上面實(shí)驗(yàn)可以發(fā)現(xiàn)，M隨著N的增大而增大時(shí)，分辨率提高，方差變大。BT法仍然沒有解決分辨率與方差之間的矛盾，但是BT法得到的功率譜當(dāng)N為無(wú)窮大時(shí)，方差會(huì)趨向于零，即為一致估計(jì)[2]。
2.2.3 周期圖法與BT法的關(guān)系
相關(guān)函數(shù)可以寫成如式（2-9）的卷積形式： (2-9）
設(shè)序列的傅立葉變換為，則當(dāng)M=N-1時(shí)，功率譜的估計(jì)可表示為式（2-10）的形式。可以看出周期圖法可以看作BT法在取M=N-1時(shí)的特例。
?（2-10）
結(jié) 論
本文通過Matlab仿真，以一個(gè)具體的隨機(jī)信號(hào)為例，簡(jiǎn)單介紹了周期圖法、平均周期圖法、修正的平均周期圖法以及BT法的基本原理，并對(duì)這些方法的性能進(jìn)行分析。可以看出，無(wú)論是周期圖法及其改進(jìn)算法還是BT法都沒有從根本上解決分辨率與方差的矛盾。經(jīng)典功率譜估計(jì)是利用傅里葉變換估計(jì)功率譜，而我們之前分析隨機(jī)信號(hào)不滿足傅里葉變換的條件，所以經(jīng)典功率譜估計(jì)方法不得不從無(wú)限長(zhǎng)數(shù)據(jù)點(diǎn)截取有限長(zhǎng)數(shù)據(jù)點(diǎn)，加入限制條件（周期圖法實(shí)際上假定N點(diǎn)外數(shù)據(jù)周期重復(fù)、BT法假定N點(diǎn)外數(shù)據(jù)為零）來(lái)"強(qiáng)制"作傅里葉變換，這也是造成它局限性的原因。
參考資料
[1] 朱哲,鐘宏偉. 非平穩(wěn)隨機(jī)信號(hào)分析處理方法研究[J] 安徽電子信息技術(shù)學(xué)院學(xué)報(bào) 2008.6:28-28
[2] 皇甫堪. 現(xiàn)代數(shù)字信號(hào)處理[M]. 電子工業(yè)出版社
?部分matlab程序代碼：
周期圖法：（by宋同學(xué)）
1 Fs=1000; 2 f1=50; 3 f2=125; 4 f3=135; 5 N=128; 6 Nfft=N; 7 n=0:N-1; 8 t=n/Fs;%采用的時(shí)間序列 9 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 10 figure; 11 plot(n,xn);grid on;title('時(shí)域信號(hào)'); 12 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 13 f=(0:length(P)-1)*Fs/length(P);%給出頻率序列 14 figure; 15 plot(f(1:N/2),P(1:N/2));grid on;title('功率譜（dB圖）');ylabel('功率譜/dB'); 16 xlabel('頻率/Hz'); 17 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 18 f=(0:length(Pxx)-1)*Fs/length(Pxx);%給出頻率序列 19 figure; 20 plot(f(1:N/2),Pxx(1:N/2));%繪制功率譜曲線 21 xlabel('頻率/Hz');ylabel('功率譜'); 22 title('周期圖 N=128');grid on; 23 std(Pxx)^2 24 25 N=256; 26 Nfft=N; 27 n=0:N-1; 28 t=n/Fs;%采用的時(shí)間序列 29 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 30 figure; 31 plot(n,xn);grid on;title('時(shí)域信號(hào)'); 32 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 33 f=(0:length(P)-1)*Fs/length(P);%給出頻率序列 34 figure; 35 plot(f(1:N/2),P(1:N/2));grid on;title('功率譜（dB圖）');ylabel('功率譜/dB'); 36 xlabel('頻率/Hz'); 37 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 38 f=(0:length(Pxx)-1)*Fs/length(Pxx);%給出頻率序列 39 figure; 40 plot(f(1:N/2),Pxx(1:N/2));%繪制功率譜曲線 41 xlabel('頻率/Hz');ylabel('功率譜'); 42 title('周期圖 N=256');grid on; 43 std(Pxx)^2 44 45 N=512; 46 Nfft=N; 47 n=0:N-1; 48 t=n/Fs;%采用的時(shí)間序列 49 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 50 figure; 51 plot(n,xn);grid on;title('時(shí)域信號(hào)'); 52 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 53 f=(0:length(P)-1)*Fs/length(P);%給出頻率序列 54 figure; 55 plot(f(1:N/2),P(1:N/2));grid on;title('功率譜（dB圖）');ylabel('功率譜/dB'); 56 xlabel('頻率/Hz'); 57 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 58 f=(0:length(Pxx)-1)*Fs/length(Pxx);%給出頻率序列 59 figure; 60 plot(f(1:N/2),Pxx(1:N/2));%繪制功率譜曲線 61 xlabel('頻率/Hz');ylabel('功率譜'); 62 title('周期圖 N=512');grid on; 63 std(Pxx)^2 64 65 N=1024; 66 Nfft=N; 67 n=0:N-1; 68 t=n/Fs;%采用的時(shí)間序列 69 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 70 figure; 71 plot(n(1:1000),xn(1:1000));grid on;title('時(shí)域信號(hào)'); 72 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 73 f=(0:length(P)-1)*Fs/length(P);%給出頻率序列 74 figure; 75 plot(f(1:N/2),P(1:N/2));grid on;title('功率譜（dB圖）');ylabel('功率譜/dB'); 76 xlabel('頻率/Hz'); 77 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 78 f=(0:length(Pxx)-1)*Fs/length(Pxx);%給出頻率序列 79 figure; 80 plot(f(1:N/2),Pxx(1:N/2));%繪制功率譜曲線 81 xlabel('頻率/Hz');ylabel('功率譜'); 82 title('周期圖 N=1024');grid on; 83 std(Pxx)^2
平均周期圖法
修正的平均周期圖法（by周同學(xué)）
1 clear;%求1024點(diǎn)功率譜以及方差 2 Fs=1000; 3 f1=50; 4 f2=125; 5 f3=135; 6 n=0:1/Fs:1; 7 xn=cos(2*pi*f1*n)+1.5*cos(2*pi*f2*n)+cos(2*pi*f3*n)+1.5*randn(size(n)); 8 M=512;N=1024;%數(shù)據(jù)總點(diǎn)數(shù)1024，每段長(zhǎng)度M 9 window=hamming(M); 10 Pxxtotal=0; 11 L=N/(M/2)-1; 12 for i=1:1:(L-1) 13 m=abs(fft(window'.*(xn((M/2+M/2*i-M+1):(M/2+M/2*i))),M).^2)/M; 14 Pxxtotal=Pxxtotal+m; 15 end 16 window=hamming(Fs-(N-M+1)+1); 17 mend=abs(fft(window'.*(xn((N-M+1):Fs)),M).^2)/M; 18 Pxxtotal=(Pxxtotal+mend)/L; 19 Pxx=10*log10((Pxxtotal));%Fourier振幅譜平方的平均值，并轉(zhuǎn)化為dB 20 f=(0:length(Pxx)-1)*Fs/length(Pxx);%給出頻率序列 21 w1=fft(window,N); 22 w10=abs(fftshift(w1)); 23 24 plot(f(1:M/2),Pxx(1:M/2));%繪制功率譜曲線 25 xlabel('頻率/Hz');ylabel('功率譜/dB'); 26 grid on; 27 B=var(Pxxtotal)BT法（by陽(yáng)同學(xué)）
1 Fs=1000; 2 n=0:1/Fs:1; 3 x=cos(2*pi*50*n)+1.5*cos(2*pi*125*n)+cos(2*pi*135*n)+1.5*randn(size(n)); 4 nfft=1024; 5 ncxk=3*nfft/4; 6 cxn=xcorr(x,'unbiased'); 7 CXk=fft(cxn,ncxk); 8 Pxx=abs(CXk); 9 index=0:round(ncxk/2-1); 10 k=index*Fs/ncxk; 11 C=Pxx(index+1); 12 P=(Pxx(index+1)); 13 plot(k,P);grid on 14 var(C) 15 title('BT法功率譜估計(jì) N=1024'); 16 xlabel('頻率 Hz'); 17 ylabel('功率');?
轉(zhuǎn)載于:https://www.cnblogs.com/jacklu/p/5140913.html
總結(jié)
以上是生活随笔為你收集整理的经典功率谱估计及Matlab仿真的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 配置java编译环境
- 下一篇： ZZULI 1876: 蛤玮的项链 Ha