第7章功率譜估計7.1功率譜估計及其分析方法

7.2經典譜估計法

7.3改進的非參數化方法

7.4模型參數化方法

7.5自相關矩陣的本征分析法

7.6小結

7.1功率譜估計及其分析方法

如果用φxx(m)表示隨機信號x(n)的自相關函數，Pxx(ω)表示它的功率譜密度，則有

(7.1)對于平穩隨機過程，根據各態遍歷假設，集合的平均可以用時間的平均代替，于是上式可寫成(7.2)在N→∞的極限情況下，Pxx(ω)是不可能收斂的，這是因為對于無限時域的隨機信號，它的傅里葉變換是不存在的。由于實際的隨機信號只能是它的一個實現或一個樣本序列的片段，因此，根據有限個樣本序列估計的信號功率譜密度為(7.3)式中，為(7.4)故(7.5)式中，XN(ω)是有限長序列XN(n)(n=0,1,…,N－1)的傅里葉變換。功率譜估計已被廣泛地應用于各種信號處理中。在信號處理的許多場所，要求預先知道信號的功率譜密度。功率譜估計也常被用來得到線性系統的參數估計。由于白噪聲的自功率譜PSD為一常數，即σω2，于是有

Pyy(ω)=σω2|H(ω)|2

(7.6)

因此，通過估計輸出信號的自功率譜(PSD)，可以估計出系統的頻率特性|H(ω)|。

7.2經典譜估計法

經典功率譜估計法可分成兩種:一種是自相關函數估計法，或稱為BT法;另一種是基于DFT的周期圖法，或稱為直接法。

7.2.1自相關函數估計法

設觀察到N個樣本序列{xn}的值為x(0)，x(1)，…，x(N－1)，現需由此N個數據來估計自相關函數φxx(m)。由于xn只能觀察到0≤n≤N－1的N個值，而n<0與n>N－1時的xn值是未知的，一般只能假定為零。根據自相關函數的定義得到

(7.7)由于x(n)只有N個觀測值，因此對于每一個固定延遲m，可以利用的數據只有(N－|m|－1)個，且在[0，N－1]范圍內，所以實際計算為(7.8)

考慮乘積項的長度，自相關序列的估計為(7.9)式中，m取絕對值是因為φxx(m)=φxx(－m)，m為負值時上式仍適用。規定的求和上下限的原則是保持充分利用全部數據。

【例7.1】

采用自相關函數估計法，求帶有白噪聲干擾的頻率為10Hz的正弦信號的功率譜。

MATLAB程序如下:

%MATLABPROGRAM7-1

clf;f=10;

N=1000;Fs=500;%數據長度和采樣頻率

n=[0:N－1];t=n/Fs;%時間序列

Lag=100;%延遲樣本點數

randn(′state′，0);%設置產生隨機數的初始狀態

x=sin(2*pi*f*t)+0.6*randn(1，length(t));%原始信號

[c，lags]=xcorr(x，Lag，′unbiased′);%對原始信號進行無偏自相關估計

subplot(311)，plot(t，x);%繪原始信號x

xlabel(′t/s′);ylabel(′x(t)′);grid;

legend(′含噪聲的信號x(t)′);

subplot(312);plot(lags/Fs,c);%繪x信號自相關，lags/Fs為時間序列

xlabel(′t/s′);ylabel(′Rxx(t)′);

legend(′信號的自相關Rxx′);grid;

Pxx=fft(c，length(lags));%利用FFT變換計算信號的功率譜

fp=(0:length(Pxx)－1)′*Fs/length(Pxx);%求功率譜的橫坐標f

Pxmag=abs(Pxx);%求幅值

subplot(313);

plot(fp(1:length(Pxx)/2),Pxmag(1:length(Pxx)/2));%繪制功率譜曲線

xlabel(′f/Hz′);ylabel(′|Pxx|′);grid;

legend(′信號的功率譜Pxx′);

程序運行結果如圖7.1所示。圖7.1自相關函數估計間接法求功率譜密度7.2.2周期圖法

周期圖法是直接將離散信號x(n)進行傅里葉變換來求取功率譜估計。假設x(n)為有限長隨機信號序列，其功率譜估計可表示為

(7.10)式中，xN(n)與xN(n+m)的下標N表示它們是有限長序列。令l=n+m，則(7.11)式中，(7.12)即XN(ω)是有限長序列x(n)的傅里葉變換。顯然，XN(ω)是周期性的。直接將XN(ω)的模的平方除以N求得功率譜的估計稱為周期圖，并用IN(ω)表示為(7.13)

【例7.2】

利用FFT算法求信號x(t)=sin(2πf1t)+cos(2πf2t)+u(t)的功率譜。其中，f1=40Hz，f2=80Hz，u(t)為白噪聲，采樣頻率fs=1kHz。信號長度取256和1024。

MATLAB程序如下:

%MATLABPROGRAM7-2

clf;

Fs=1000;

f1=40;f2=80;

%Case1:N=256

N=256;Nfft=256;

n=[0:N-1];

t=n/Fs;

xn=sin(2*pi*f1*t)+cos(2*pi*f2*t)+randn(1，N);

Pxx=10*log10(abs(fft(xn，Nfft).^2)/(N+1));

f=(0:length(Pxx)－1)*Fs/length(Pxx);

subplot(211);

plot(f，Pxx);

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′N=256′);grid;

%Case2:N=1024

N=1024;Nfft=1024;

n=[0:N－1];

t=n/Fs;

xn=sin(2*pi*f1*t)+cos(2*pi*f2*t)+randn(1，N);

Pxx=10*log10(abs(fft(xn，Nfft).^2)/N);

f=(0:length(Pxx)－1)*Fs/length(Pxx);

subplot(212);

plot(f，Pxx);

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′N=1024′);grid;

程序運行結果如圖7.2所示。由圖可以看出，在頻率40Hz和80Hz處功率譜有兩個峰值，說明信號中含有40Hz和80Hz的周期成分。但是功率譜密度都在很大范圍內波動，而且并沒有因信號取樣點數N的增加而有明顯改進。圖7.2DFT周期圖法功率譜估計周期圖的數學期望值可表示為(7.14)式中，RN代表矩形序列。令(7.15)ω(m)是兩個矩形函數的卷積，為一個三角函數，稱為巴特利特(Bartlett)窗函數。

7.3改進的非參數化方法

7.3.1分段平均周期圖法

如果x1，x2，…，xL是不相關的隨機變量，每一個具有期望值μ，方差σ2，則它們的數學平均的期望值等于μ，數學平均的方差等于σ2/L，即(7.16)(7.17)當L→∞時，Var[

]→0，可達到一致譜估計的目的。因而降低估計量方差的一種有效方法是將若干個獨立估計值進行平均。把這種方法應用于譜估計就是Bartlett平均周期圖法。

Bartlett平均周期圖法是將序列x(n)(0≤n≤N－1)分段求周期圖再平均。設將x(n)分成L段，每段有M個樣本，因而N=LM，第i段樣本序列可寫成

xi(n)=x(n+iM－M)0≤n≤M－1,1≤i≤L

(7.18)

第i段的周期圖為(7.19)一般可假定各段的周期圖IiM(ω)是相互獨立的。譜估計可定義為L段周期圖的平均，即(7.20)它的數學期望值為

(7.21)由此得到(7.22)式中，M=N/L，因此Bartlett估計的期望值是真實譜Pxx(ω)與三角窗函數的卷積。由于三角窗函數不等于δ函數，所以Bartlett估計也是有偏估計，但當N→∞時是無偏估計。假定各段周期圖是相互獨立的，可得到周期圖的方差為(7.23)

【例7.3】利用分段平均周期圖法求信號x(t)=sin(2πf1t)+sin(2πf2t)+u(t)的功率譜。其中，f1=60Hz，f2=120Hz,u(t)為白噪聲,采樣頻率fs=1kHz，信號長度為1024。

MATLAB程序如下:

%MATLABPROGRAM7-3

clf;

Fs=1000;

f1=60;f2=120;

N=1024;Nsec=256;%分四段

n=[0:N-1];

t=n/Fs;

xn=sin(2*pi*f1*t)+sin(2*pi*f2*t)+randn(1，N);

Pxx1=abs(fft(xn(1:256)，Nsec).^2)/Nsec;

Pxx2=abs(fft(xn(257:512)，Nsec).^2)/Nsec;

Pxx3=abs(fft(xn(513:768)，Nsec).^2)/Nsec;

Pxx4=abs(fft(xn(769:1024)，Nsec).^2)/Nsec;

Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);

f=(0:length(Pxx)-1)*Fs/length(Pxx);

subplot(211);

plot(f，Pxx);

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′N=256*4′);grid;

程序運行結果如圖7.3所示。圖7.3分段平均周期圖法功率譜估計7.3.2加窗平均周期圖法

采用一合適的窗函數W(ejω)與周期圖IN(ω)卷積，可得到譜估計的另一種常用方法——平滑周期圖。(7.24)(7.25)式中，與ω(m)分別是IN(ω)與W(ejω)的傅里葉反變換。是一個實偶非負函數，必有窗函數ω(m)是一個非矩形窗的偶序列，其長度為(2M－1)，且其傅里葉變換非負。

【例7.4】

利用加窗平均周期圖法求例7.3的隨機信號序列的功率譜。

MATLAB程序如下:

%MATLABPROGRAM7-4

clf;

Fs=1000;

f1=60;f2=120;

N=1024;Nsec=256;%分四段

n=[0:N－1];

t=n/Fs;

xn=sin(2*pi*f1*t)+sin(2*pi*f2*t)+randn(1，N);

w=triang(256)′;%采用三角窗加窗處理

Pxx1=abs(fft(w.*xn(1:256)，Nsec).^2)/norm(w)^2;

Pxx2=abs(fft(w.*xn(257:512)，Nsec).^2)/norm(w)^2;

Pxx3=abs(fft(w.*xn(513:768)，Nsec).^2)/norm(w)^2;

Pxx4=abs(fft(w.*xn(769:1024)，Nsec).^2)/norm(w)^2;

Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);

f=(0:length(Pxx)-1)*Fs/length(Pxx);

subplot(211);

plot(f，Pxx);

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′N=256*4′);grid;

程序運行結果如圖7.4所示。相比分段平均周期圖法，采用加窗處理后，譜峰加寬，而噪聲譜均在0dB附近，更為平坦。圖7.4加窗平均周期圖法功率譜估計7.3.3Welch法

Welch法采用信號重疊分段、加窗函數以及FFT等算法來計算一個信號序列的自功率譜(PSD)和兩個信號序列的互功率譜(CSD)。在重疊分段的基礎上，選擇適當的窗函數ω(n)，加進周期圖計算中，得到每一段的周期圖(7.26)式中，為歸一化因子，這使得無論什么樣的窗函數均可使譜估計非負。分段重疊處理會使方差減小，一般取重疊部分長度為分段度度的二分之一。

【例7.5】

對一個256點的白噪聲序列，取每一分段長度為64，重疊點數為32，利用Welch法求其功率譜估計。繪制該功率譜曲線，并與256點一次求周期圖的功率譜曲線進行比較。

MATLAB程序如下:

%MATLABPROGRAM7-5

clc;

N=256;Fs=1;

u=rand(1，N);

u=u-mean(u);

%用Welch法平均估計實驗數據的功率譜

[xpsd，w]=pwelch(u，hamming(64)，32，4096，1);

mmax=max(xpsd);

xpsd=xpsd/mmax;

xpsd=10*log10(xpsd+0.000001);

subplot(211);

plot(w，xpsd);

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′重疊分段，分段值為64′);grid;

%256點一次求周期圖的功率譜

[xpsd，w]=pwelch(u，hamming(256)，0，4096，1);

mmax=max(xpsd);

xpsd=xpsd/mmax;

xpsd=10*log10(xpsd+0.000001);

subplot(212);

plot(w，xpsd);

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′全周期N=256′);grid;

程序運行結果如圖7.5所示。從圖中可以看出，Welch法比周期圖法的曲線更平滑，功率譜估計效果最好。圖7.5Welch法的功率譜估計

MATLAB還可以利用函數psd和csd來實現Welch法的功率譜估計。

(1)函數psd利用Welch法來

估計一個信號的自功率譜，調用格式為

Pxx=psd(x)

Pxx=psd(x，Nfft)

Pxx=psd(x，Nfft，Fs，windows，Noverlap)

[Pxx，F]=psd(x，Nfft，Fs，windows，Noverlap，dflag)

psd(x，…，dflag)

其中，x為信號序列;windows定義窗函數和x分段序列的長度，窗函數長度必須小于或等于Nfft，否則會出錯;dflag為取出信號趨勢分量的選擇項:linear—去除直線趨勢分量，mean—去除均值分量，none—不作去除趨勢處理。其余參數與函數pwelch中的相同。

【例7.6】

求取隨機信號序列x(n)=sin(2πnf)+e－0.2n+u(n)，n=0，1，…，N－1的功率譜估計，其中，u(n)為白噪聲序列，歸一化頻率i=0.3，N=1024。

MATLAB程序如下:

%MATLABPROGRAM7-6

clf;

Fs=1;f=0.3;

N=1024;Nfft=256;n=[0:N－1];

xn=sin(2*pi*f*n)+exp(－0.2*n)+randn(1，N);

windows=hanning(256);

noverlap=128;

dflag=′none′;

[Pxx，F]=psd(xn，Nfft，Fs，windows，noverlap，dflag);

plot(F，10*log10(Pxx));

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′WelchMethod(N=1024)′);grid;

程序運行結果如圖7.6所示。由圖可知，采用Welch法求得的功率譜與前面所述方法相比，效果最好，使信號得以突出，從而抑制了噪聲。圖7.6函數psd實現Welch法功率譜估計

(2)函數csd利用Welch法估計兩個隨機信號的互功率譜密度，調用格式為

Pxy=csd(x，y)

Pxy=csd(x，y，Nfft)

Pxy=csd(x，y，Nfft，Fs，windows，Noverlap)

[Pxy，F]=csd(x，Nfft，Fs，windows，Noverlap，dflag)

csd(x，y，…，dflag)

其中，x和y分別為隨機信號;Pxy為x和y的互功率譜估計。其他參數同函數psd。

【例7.7】

產生兩個長度為8000的白噪聲信號和一個帶有白噪聲的1kHz的周期信號，求兩個白噪聲信號及白噪聲與帶有噪聲的周期信號的互功率譜，并繪制互功率譜曲線。FFT所采用的長度為1024，采用500個點的三角窗，并且沒有重疊，采樣頻率fs=6kHz。

MATLAB程序如下:

%MATLABPROGRAM7-7

clf;

Fs=6000;

%信號采樣頻率

randn(′state′，0);%設置隨機狀態初始值

x=randn(8000，1);%產生第一個白噪聲信號

y=randn(8000，1);%產生第二個白噪聲信號

z=2*sin(2*pi*1000*[1:8000]′/Fs)+randn(8000，1);

%產生含周期成分的噪聲信號

[Pxy，f]=csd(x，y，1024，Fs，triang(500)，0);

%白噪聲信號互功率譜，無重疊三角窗

[Pxz，f]=csd(x，z，1024，Fs，triang(500)，0);

%噪聲與周期信號互功率譜，無重疊三角窗

subplot(211);plot(f，10*log10(Pxy));grid;%繪制x，y信號互功率譜

xlabel(′f/Hz′);ylabel(′互相關譜/dB′);title('噪聲信號的互功率譜');

subplot(212);plot(f，10*log10(Pxz));grid;%繪制x，z信號的互相關譜

xlabel(′f/Hz′);ylabel(′互相關譜/dB′);title('帶噪聲周期信號的互相關譜');

程序運行結果如圖7.7所示。圖7.7兩個噪聲信號及噪聲信號和含噪聲周期信號的互功率譜7.3.4多窗口法

多窗口法(MultitaperMethod，MTM法)利用多個正交窗口(tapers)獲得各自獨立的近似功率譜估計，然后綜合這些估計最終得到一個序列的功率譜估計。

MATLAB信號處理工具箱提供函數pmtm實現MTM法估計功率譜密度，調用格式為

Pxx=pmtm(x)

Pxx=pmtm(x，nw，Nfft)

[Pxx，F]=pmtm(x，nw，Nfft，Fs)

其中，x為信號序列;nw為時間帶寬積，缺省值為4，通常可取2、5/2、3、7/2等;Nfft為FFT的長度;Fs為采樣頻率。

【例7.8】

采用多窗口(MTM)法求取例7.3中隨機信號序列的功率譜估計。

MATLAB程序如下:

%MATLABPROGRAM7-8

clf;

Fs=1000;

f1=60;f2=120;

N=1024;Nfft=256;n=[0:N－1];t=n/Fs;%數據長度，分段數據長度，時間序列

randn(′state′，0);

%設置產生隨機數的初始狀態

xn=sin(2*pi*f1*t)+sin(2*pi*f2*t)+randn(1，N);

[Pxx1，F]=pmtm(xn,4,Nfft,Fs);%用多窗口法(NW=4)估計功率譜

plot(F，10*log10(Pxx1));

%繪制功率譜

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′MTM(NW=4)′);grid;

[Pxx2，F]=pmtm(xn，2，Nfft，Fs);

%用多窗口法(NW=2)估計功率譜

plot(F，10*log10(Pxx2));%繪制功率譜

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′MTM(NW=2)′);grid;

程序運行結果如圖7.8所示。由圖可知，多窗口法得到了較好的功率譜估計。當帶寬NW=4時，譜估計采用較多的窗口，使得每個窗口的數據點數均減少了，因此頻率域分辨率降低了，且功率譜的波動較小。當NW=2時，譜估計采用較少的窗口，每個窗口的數據點數均增多了，使得頻率域的分辨率提高了，但噪聲譜的分辨率也隨之增高了，因而功率譜具有相對較大的波動。圖7.8多窗口(MTM)法的功率譜估計

7.4模型參數化方法

7.4.1AR模型法

在實際中所遇到的隨機過程，常常可以用一個具有有理分式的傳遞函數的模型來表示，因此，可以用一個線性差分方程作為產生隨機序列x(n)的系統模型

(7.27)式中，ω(n)表示白噪聲，將上式變換到z域則有(7.28)該模型的傳遞函數為(7.29)式中，。

當輸入的白噪聲的功率譜密度為Pωω(z)=σω2時，輸出的功率譜密度為(7.30)將z=ejω代入式(7.30)，得(7.31)

設a0=b0=1，所有其余的bl均為零，則(7.32)式(7.32)的形式使它被稱為p階自回歸模型，簡稱AR模型。將式(7.32)進行Z變換，可得AR模型的傳遞函數為

(7.33)自回歸模型的H(z)只有極點，沒有零點，因此又稱為全極點模型。當采用自回歸模型時，輸出的功率譜密度為(7.34)此時，只要能求得σ2ω及所有的ak參量，就可求得Pxx(ω)。由此可見，用模型法作功率譜估計，實際上要解決的是模型的參數估計問題。為了得到AR模型的功率譜估計Pxx(ω)，必須求取參數a1，a2，a3，…，ap及σ2ω。AR模型參數的估計可以通過自相關法(Yule-Walker法)來求解。自相關法的計算簡單，但譜估計的分辨率相對較差。AR參數與自相關函數φxx(m)之間可用Yule-Walker方程來表示:(7.35)由式(7.32)可知，x(n)只與ω(n)相關，而與ω(n+m)(m≥1)無關，故式(7.35)中的第二項為(7.36)代入式(7.35)得(7.37)將m=1，2，…，p分別代入式(7.37)并寫成矩陣形式，得Yule-Walker方程為(7.38)將式(7.30)與式(7.37)合在一起寫成歸一化的正規矩陣形式，得(7.39)式(7.39)的系數矩陣用[φ]p+1表示，稱為自相關矩陣，是對稱矩陣，也稱(p+1)×(p+1)的托普尼茲(Toeplitz)矩陣，其與主對角線平行的斜對角線上的元素都是相同的，因此存在高效算法，其中應用廣泛的有Levinson-Durbin算法。

【例7.9】

已知隨機序列x(t)=3sin(2πf1t)+

2sin(2πf2t)+u(t)，f1=20Hz，f2=100Hz，u(t)為白噪聲，采樣間隔為0.002s，長度N=2048。當AR模型的階次分別為8和14時，用自相關法求解AR模型的系數a以及功率譜估計。

MATLAB程序如下:

%MATLABPROGRAM7-9

clf;

dalt=0.002;Fs=1/dalt;%采樣頻率

N=2048;

t=[0:dalt:dalt*(N－1)];

f1=20;f2=100;

x=3*sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+randn(1，N);%生成輸入信號

subplot(311);plot(t，x，′k′);

xlabel(′t/s′);ylabel(′x(t)′);legend(′x(t)′);grid;

%用自相關法求AR模型的系數a和輸入白噪聲的功率譜E

[a，E]=aryule(x，8);%AR模型的階次為8

%根據定義計算信號的功率譜密度

fpx=abs(fft(a，N)).^2;Pxx=E./fpx;

f=(0:length(Pxx)－1)*Fs/length(Pxx);%計算各點對應的頻率值

subplot(312);plot(f，10*log10(Pxx));

xlabel(′f/Hz′);ylabel(′Pxx/dB′);

legend(′ARorder=8′);grid;

[a，E]=aryule(x，14);%AR模型的階次為14

fpx=abs(fft(a，N)).^2;Pxx=E./fpx;

f=(0:length(Pxx)－1)*Fs/length(Pxx);

subplot(313);plot(f，10*log10(Pxx));

xlabel(′f/Hz′);ylabel(′Pxx/dB′);

legend(′ARorder=14′);grid;

程序運行結果如圖7.9所示。由圖可知，AR模型法的功率譜估計很好，曲線更平滑。當AR模型的階數取值更大時，效果更佳。圖7.9AR模型的自相關法求功率譜估計由線性預測器的定義可知，參數ak與σ2ω分別等于x(n)的p階線性預測器的系數apk與最小均方誤差E[e2(n)]min，故AR模型法又稱為線性預測AR模型法。因此，這種估計功率譜密度的方法也可歸結為求在最小均方誤差準則下的線性預測器的系數apk的問題。

線性預測誤差為

(7.40)

將上式進行Z變換，得(7.41)于是，有

(7.42)由式(7.42)可知，He(z)是以x(n)作為輸入信號、誤差e(n)作為輸出信號的濾波器的傳遞函數，該濾波器稱為預測誤差濾波器。當apk按最小均方誤差準則時，有apk=ak

，故(7.43)即預測誤差濾波器是x(n)的形成系統的逆濾波器,因此可得到

e(n)=u(n)

(7.44)即將x(n)送入預測誤差濾波器，其輸出e(n)等于系統的激勵信號白噪聲u(n)。7.4.2最大熵功率譜估計法

熵代表一種不定度，最大熵為最大不定度，即它的時間序列最隨機，而它的PSD應最平伏。按熵的定義，當隨機信號x的取值為離散時，熵H定義為(7.45)式中,pi為出現狀態i的概率。當x的取值為連續時,熵定義為(7.46)式中，p(x)為概率密度函數。對于時間序列x0，x1，…，xN，概率密度應由聯合概率密度函數p(x0，x1，…，xN)代替。最大熵功率譜估計法假定隨機過程是平穩高斯過程。假設{xn}為零均值且呈高斯分布的隨機過程，則它的一維高斯分布為(7.47)式中，(7.48)

N維高斯分布為(7.49)式中，(7.50)又因為(7.51)則一維高斯分布的熵為(7.52)

N維高斯分布信號的熵為(7.53)式中，detΦ代表矩陣Φ的行列式，若要使熵H最大，就要求detΦ最大。若已知φxx(0)，φxx(1)，…，φxx(N)，由于自相關函數的矩陣必是正定的，故矩陣Φ(N+1)的行列必大于零，即

(7.54)

為了得到最大熵，用φxx(N+1)對式(7.54)求導，使得，由此得到det[Φ(N+1)]最大的φxx(N+1)，滿足下列方程:(7.55)式(7.55)是φxx(N+1)的一次函數，由此式可解出φxx(N+1)。同理可以類推φxx(N+2)。因此，若每步都按最大熵的原則外推一個自相關序列的值，則可以外推到任意多個數值而不必認為它們是零，這就是最大熵功率譜估計法的基本思想。

【例7.10】

用最大熵功率譜估計法估計例7.6中的隨機信號序列功率譜。

MATLAB程序如下:

%MATLABPROGRAM7-10

clf;

f=0.3;

N=1024;Nfft=256;

n=[0:N－1];

randn(′state′，0);

xn=sin(2*pi*f*n)+exp(－0.2*n)+randn(1，N);

[Pxx，F]=pmem(xn，14，Nfft，Fs);

plot(F，10*log10(Pxx));

xlabel(′f/Hz′);ylabel(′功率譜/dB′);

title(′MaxmumEntropyMethod(order=14)′);

grid;

程序執行結果如圖7.10所示。圖7.10最大熵功率譜估計最大熵功率譜估計法和Welch功率譜估計法相比，最大熵功率譜曲線較光滑。本例中選定的AR模型的階數order必須大于10。7.4.3Levinson-Durbin遞推算法

最大熵功率譜估計法與線性預測AR模型法是等價的，它們都可歸結為求解Yule-Walker方程中各AR系數ak(k=1，2，…，p)的問題。若直接從Yule-Walker方程求解參數，需要作求逆矩陣的運算，當p很大時，運算量也很大，而且當模型階數增加一階時，矩陣會增大一維，需要全部重新計算。Levinson-Durbin算法對Yule-Walker方程提供了一個高效率的解法。

Levinson-Durbin算法是一種遞推求解算法。遞推算法先從一階AR模型求得一階參數a11及σ21，再從二階AR模型的矩陣方程解得a22、a21、σ2，最后求得所要求的p階模型參數{ap1，ap2，…，app，σ2p}。其中，參數a的第一個下標是指AR模型的階數。

遞推公式為

(7.56)(7.57)(7.58)一般階數是未知的，當遞推到第k階時，σk2滿足所允許的值，就可選階數p=k。若信號的模型是p階的AR模型，則應有

akk=0，σk2=σ2p

(7.59)

式(7.59)說明σ2p已達最小均方誤差值。由于σ2k>0，故對于任何k必有

|akk|<1k=1，2，…，p

(7.60)

式(7.60)正是He(z)的所有極點均在單位圓內(即穩定性)的充分必要條件，同時也是Φ為正定矩陣的充分必要條件。

【例7.11】設信號x(n)是一個帶白噪聲的4階IIR濾波器的脈沖響應，IIR濾波器的傳遞函數無零點，其分母多項式系數為a=[10.10.10.10.1]，用線性預測法建立其AR模型。

MATLAB程序如下:

%MATLABPROGRAM7-11

clc;

randn(′state′，0);

%設置隨機函數的狀態

a=[1，0.1，0.1，0.1，0.1];

%濾波器分母多項式系數

x=impz(1,a,10)+randn(10,1)/20;%求得帶噪聲的濾波器脈沖響應

r=xcorr(x);%求信號x的自相關序列rxx

r(1:length(x)－1)=[];%該序列的前面部分受到邊界的影響，剔除

aa=levinson(r，4)%采用Levison-Durbin遞推算法求出AR模型系數aa程序運行結果為

aa=

1.00000.18490.12790.11140.1839

若信號不包含隨機噪聲，則用上面程序求得的信號AR模型和所采用的全極點濾波器模型完全相同。即使加入隨機噪聲，得到的結果也與原模型參數a接近。7.4.4Burg遞推算法

Levinson-Durbin遞推算法求解Yule-Walker方程中的AR系數，雖可以簡化計算，但需要已知自相關序列φxx(m)。實際上，自相關序列φxx(m)只能從時間序列x(n)的有限個數據中得到它的估計值。當時間序列較短時，的估計誤差很大，這將給AR參數ak的計算引入很大誤差，導致功率譜估計出現譜線分裂與譜峰頻率偏移等現象。

最大熵功率譜估計法只說明了如何從已知的N個φxx(m)外推N個以外的φxx(m)值，從而得到高分辨率的功率譜，但并未涉及如何從有限時間序列來得到這N個φxx(m)的問題。Burg提出了一種直接由時間序列計算AR模型參數的方法，稱為Burg算法。

Burg遞推算法的優點是不需要估計自相關函數，可以直接從已知的x(n)序列求得參數Kp，并可保證滿足穩定性的充要條件:|Kp|<1。Burg算法是以前向均方誤差與后向均方誤差之和最小為準則求得Kp的。

Burg遞推算法對于短的時間序列x(n)能得到較正確的估計，因此得到普遍應用。但Burg算法受Levinson關系式的約束，因此不能完全克服Levinson-Durbin算法中的缺點，仍存在某些譜線分裂與頻率偏移的現象。

【例7.12】

用Burg遞推算法估計例7.9中AR模型的功率譜。

MATLAB程序如下:

%MATLABPROGRAM7-12

clf;

dalt=0.002;Fs=1/dalt;%采樣頻率

N=2048;

t=[0:dalt:dalt*(N－1)];

f1=20;f2=100;

x=3*sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+randn(1，N);%生成輸入信號

subplot(311);plot(t，x，′k′);

xlabel(′t/s′);ylabel(′x(t)′);legend(′x(t)′);grid;

%使用Burg算法得到功率譜估計

[xpsd，F]=pburg(x，8，N，Fs);%AR模型階數為8

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

subplot(312);

plot(F，xpsd);xlabel(′f/Hz′);ylabel(′Pxx/dB′);

legend(′ARorder=8′);grid;

[xpsd，F]=pburg(x，14，N，Fs);%AR模型階數為14

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

subplot(313);

plot(F，xpsd);xlabel(′f/Hz′);ylabel(′Pxx/dB′);

legend(′ARorder=14′);grid;程序運行結果如圖7.11所示。由圖可知，Burg遞推算法的功率譜估計曲線更平滑。當AR模型的階數取值較大時，功率譜估計更接近于實際值。圖7.11Burg算法實現功率譜估計

7.5自相關矩陣的本征分析法

在功率譜估計的諸多實現方法中，基于自相關矩陣特征值Vk的多信號分類(MultipeSignalClassification，MUSIC)法和特征向量(EigenVector，EV)法是兩種重要的方法。對于含有復正弦信號和白噪聲的系統，系統自相關矩陣由信號自相關矩陣和噪聲自相關矩陣兩部分組成，即系統自相關矩陣R包含有兩個子空間信息:信號子空間和噪聲子空間。

MUSIC法和EV法原理完全相同。MUSIC法計算功率譜密度為

(7.61)

EV法采用特征值加權法計算功率譜密度，為(7.62)式(7.61)和式(7.62)中，e(f)為復正弦向量;V為輸入信號自相關矩陣的特征向量;H為共軛轉置運算;λk為特征值函數。

【例7.13】

試用MUSIC法求信號x(t)=2sin(2πf1t)+cos(2πf2t)+u(t)的功率譜估計。其中，f1=50Hz，f2=100Hz,u(t)為白噪聲，采樣頻率fs=0.5kHz，信號長度為256。

MATLAB程序如下:

%MATLABPROGRAM7-13

clf;

N=256;Nfft=256;

Fs=500;f1=50;f2=100;

n=[0:N-1];t=n/Fs;

randn(′state′，0);

xn=2*sin(2*pi*f1*t)+cos(2*pi*f2*t)+randn(1，N);

[Pxx，f]=pmusic(xn，[7，1.1]，Nfft，Fs，32，16);

Pxx=10*log10(Pxx);

plot(f，Pxx);

xlabel(′f/Hz′);ylabel(′功率譜/dB′);title(′MusicMethod′);

grid;

程序運行結果如圖7.12所示。圖7.12MUSIC法實現功率譜估計

【例7.14】

已知信號x(n)=0.5ejπn/2+2ejπn/5+3ejπn/10+

u(n)，n=1，2，…，N，其中，u(n)為白噪聲，信號長度N=1024。試用EV法求信號的功率譜估計。

M