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

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

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

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

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

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

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

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

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