功率譜估計(jì)性能分析及其MATLAB實(shí)現(xiàn)

1、功率譜估計(jì)性能分析及其MATLAB實(shí)現(xiàn)一、 經(jīng)典功率譜估計(jì)分類簡(jiǎn)介1 間接法根據(jù)維納-辛欽定理，1958年Blackman和Turkey給出了這一方法的具體實(shí)現(xiàn)，即先由N個(gè)觀察值xN(n)，估計(jì)出自相關(guān)函數(shù)rx(m)，求自相關(guān)函數(shù)傅里葉變換，以此變換結(jié)果作為對(duì)功率譜Px(w)的估計(jì)。2 直接法直接法功率譜估計(jì)是間接法功率譜估計(jì)的一個(gè)特例，又稱為周期圖法，它是把隨機(jī)信號(hào)的N個(gè)觀察值xN(n)直接進(jìn)行傅里葉變換，得到XN(ejw)，然后取其幅值的平方，再除以N，作為對(duì)功率譜Px(w)的估計(jì)。3 改進(jìn)的周期圖法將N點(diǎn)的觀察值分成L個(gè)數(shù)據(jù)段，每段的數(shù)據(jù)為M，然后計(jì)算L個(gè)數(shù)據(jù)段的周期圖的平均Pper(w

2、)，作為功率譜的估計(jì)，以此來(lái)改善用N點(diǎn)觀察數(shù)據(jù)直接計(jì)算的周期圖Pper(w)的方差特性。根據(jù)分段方法的不同，又可以分為Welch法和Bartlett法。Welch法所分的數(shù)據(jù)段可以互相重疊，選用的數(shù)據(jù)窗可以是任意窗。Bartlett法所分的數(shù)據(jù)段互不重疊，選用的數(shù)據(jù)窗是矩形窗。二、 經(jīng)典功率譜估計(jì)的性能比較1 仿真結(jié)果為了比較經(jīng)典功率譜估計(jì)的性能，本文采用的信號(hào)是高斯白噪聲加兩個(gè)正弦信號(hào)，采樣率Fs=1000Hz，兩個(gè)正弦信號(hào)的頻率分別為f1=200Hz，f2=210Hz。所用數(shù)據(jù)長(zhǎng)度N=400.仿真結(jié)果如下：(a)(b)(c)(d)(e)(f)Figure1 經(jīng)典功率譜估計(jì)的仿真結(jié)果Figu

3、re1(a)示出了待估計(jì)信號(hào)的時(shí)域波形；Figure2(b)示出了用該數(shù)據(jù)段直接求出的周期圖，所用的數(shù)據(jù)窗為矩形窗；Figure2(c)是用BT法（間接法）求出的功率譜曲線，對(duì)自相關(guān)函數(shù)用的平滑窗為矩形窗，長(zhǎng)度M=128，數(shù)據(jù)沒(méi)有加窗；Figure2(d)是用BT法（間接法）求出的功率譜曲線，對(duì)自相關(guān)函數(shù)用的平滑窗為Hamming窗，長(zhǎng)度M=64，數(shù)據(jù)沒(méi)有加窗；Figure2(e)是用Welch平均法求出的功率譜曲線，每段數(shù)據(jù)的長(zhǎng)度為64點(diǎn)，重疊32點(diǎn)，使用的Hamming窗；Figure2(f)是用Welch平均法求出的功率譜曲線，每段數(shù)據(jù)的長(zhǎng)度為100點(diǎn)，重疊48點(diǎn)，使用的Hamming窗

4、；2 性能比較1) 直接法得到的功率譜分辨率最高，但是方差性能最差，功率譜起伏劇烈，容易出現(xiàn)虛假譜峰；2) 間接法由于使用了平滑窗對(duì)直接法估計(jì)的功率譜進(jìn)行了平滑，因此方差性能比直接法好，功率譜比直接法估計(jì)的要平滑，但其分辨率比直接法低。3) Welch平均周期圖法是三種經(jīng)典功率譜估計(jì)方法中方差性能最好的，估計(jì)的功率譜也最為平滑，但這是以分辨率的下降及偏差的增大為代價(jià)的。3 關(guān)于經(jīng)典功率譜估計(jì)的總結(jié)1) 功率譜估計(jì)，不論是直接法還是間接法都可以用FFT快速計(jì)算，且物理概念明確，因而仍是目前較常用的譜估計(jì)方法。2) 譜的分辨率較低，它正比于2/N，N是所使用的數(shù)據(jù)長(zhǎng)度。3) 方差性能不好，不是真實(shí)

5、功率譜的一致估計(jì)，且N增大時(shí)，功率譜起伏加劇。4) 周期圖的平滑和平均是和窗函數(shù)的使用緊密關(guān)聯(lián)的，平滑和平均主要是用來(lái)改善周期圖的方差性能，但往往又減小了分辨率和增加了偏差，沒(méi)有一個(gè)窗函數(shù)能使估計(jì)的功率譜在方差、偏差和分辨率各個(gè)方面都得到改善，因此使用窗函數(shù)只是改進(jìn)估計(jì)質(zhì)量的一個(gè)技巧問(wèn)題，并不能從根本上解決問(wèn)題。三、 AR模型功率譜估計(jì)1 AR模型功率譜估計(jì)簡(jiǎn)介AR模型功率譜估計(jì)是現(xiàn)代譜估計(jì)中最常用的一種方法，這是因?yàn)锳R模型參數(shù)的精確估計(jì)可以用解一組線性方程（Yule-Walker方程）的方法求得。其核心思想是：將信號(hào)看成是一個(gè)p階AR過(guò)程，通過(guò)建立Yule-Walker方程求解AR模型的參

6、數(shù)，從而得到功率譜的估計(jì)。由于已知的僅僅是長(zhǎng)度有限的觀測(cè)數(shù)據(jù)，因此AR模型參數(shù)的求得，通常是首先通過(guò)某種算法求得自相關(guān)函數(shù)的估計(jì)值，進(jìn)而求得AR模型參數(shù)的估計(jì)值。常用的幾種AR模型參數(shù)提取方法有：1) 自相關(guān)法假定觀測(cè)數(shù)據(jù)區(qū)間之外的數(shù)據(jù)為0，在均方誤差意義下使得數(shù)據(jù)的前向預(yù)測(cè)誤差最小。由此估計(jì)的自相關(guān)矩陣式正定的，且具有Toeplitz性，可以用Levison-Durbin算法求解。2) 協(xié)方差法不作觀測(cè)數(shù)據(jù)區(qū)間之外的數(shù)據(jù)為0的假設(shè)，在均方誤差意義下使得數(shù)據(jù)的前向預(yù)測(cè)誤差最小。由此估計(jì)的自相關(guān)矩陣式半正定的，且不具有Toeplitz性，得到的AR模型可能不穩(wěn)定。3) 修正的協(xié)方差法不作觀測(cè)數(shù)據(jù)

7、區(qū)間之外的數(shù)據(jù)為0的假設(shè)，在均方誤差意義下使得數(shù)據(jù)的前向預(yù)測(cè)誤差與后向預(yù)測(cè)誤差之和最小。由此估計(jì)的自相關(guān)矩陣式半正定的，且不具有Toeplitz性，得到的AR模型可能不穩(wěn)定。但得到的一階AR模型是穩(wěn)定的。4) Burg法在約束AR模型的參數(shù)滿足Levison-Durbin遞歸條件的前提下，在均方誤差意義下使得數(shù)據(jù)的前向預(yù)測(cè)誤差與后向預(yù)測(cè)誤差之和最小。得到的AR模型是穩(wěn)定的，但有時(shí)可能出現(xiàn)譜線分裂現(xiàn)象。仍然用前面的仿真信號(hào)，取AR模型的階數(shù)p=48，用上述各種AR模型參數(shù)提取方法估計(jì)的功率譜如figure2所示。Figure2 AR模型功率譜估計(jì)的仿真結(jié)果2 AR模型功率譜估計(jì)與經(jīng)典功率譜估計(jì)性

8、能比較分別采用直接法和AR模型功率譜估計(jì)中的自相關(guān)法得到的上述信號(hào)的功率譜估計(jì)如figure3所示：Figure3 AR模型功率譜估計(jì)與經(jīng)典功率譜估計(jì)比較小結(jié)：1) 由于AR模型是一個(gè)有理分式，因而估計(jì)出的譜要比經(jīng)典法的譜平滑。2) AR譜估計(jì)的一些方法隱含著數(shù)據(jù)和自相關(guān)函數(shù)的外推，可獲得高的分辨率。3 關(guān)于AR模型功率譜估計(jì)總結(jié)Figure4給出了階數(shù)對(duì)AR模型功率譜估計(jì)的影響，采用的是AR模型功率譜估計(jì)中的Burg法。Figure4 不同階數(shù)的AR模型功率譜估計(jì)小結(jié)：階數(shù)越高，得到的譜的分辨率也越高，但方差也越大，將會(huì)產(chǎn)生更多的虛假譜峰；四、 附錄本文所用來(lái)仿真的MATLAB程序如下：%

9、Research On Classic PSD Estimate Methods% Author: Chen Feiqiang% Date: 2010-12-14% Generate Signalclear all; close all; clc;randn(state,0);Fs = 1000; % sample frequencyt = 0:1/Fs:.3;sigma = 1; % noise variancex = cos(2*pi*t*200) + cos(2*pi*t*210) + sigma*randn(size(t); % generate signal: % cosine +

10、noisefigure;plot(t,x), xlabel(t), ylabel(x);title(Signal In Time Domain);grid on;% Estimate PSD using periodogram methodwindow = rectwin(length(x); % window function used xn = x.*window; index = nextpow2(length(x);N = pow2(index);Xw = fft(xn, N); % do N-points FFTPw = Xw.*conj(Xw)/N; % Calculate PSD

11、k = 0:floor(N/2 - 1);figure;plot(k*Fs/N, 10*log10(Pw(k+1)/max(Pw);title(Periodogram PSD Estimate);xlabel(Frequency(Hz);ylabel(Normalized PSD(dB),grid on;hold on% Estimate PSD using BT methodwindow_a = rectwin(length(x); % window function for data x(n) xn = x.*window_a; Rx = xcorr(xn); % auto-correla

12、tion function estimateN = length(Rx);M = floor(N/4); % the length of smooth window% M = 100;window_v = rectwin(M); % smooth window for BT methodRxWin = Rx(1:M).*window_v; % smooth window multiply auto-correlation function Pw = abs(fft(RxWin, N);k = 0:floor(N/2 - 1);figure;plot(k*Fs/N, 10*log10(Pw(k+1)/max(Pw),r);title(BT Method PSD Estimate);xlabel(Frequency(Hz);ylabel(Normalized PSD(dB),grid on;% Estimate PSD using Welch methodwindow = 32; % the length of each segmentnoverlap = 8; % overlap number for each segmentnfft = pow2(nextpow2(length(x); % nfft-points FFT for ea