《MATLAB輔助現(xiàn)代工程數(shù)字信號處理》課件第7章

第7章功率譜估計(jì)7.1功率譜估計(jì)及其分析方法

7.2經(jīng)典譜估計(jì)法

7.3改進(jìn)的非參數(shù)化方法

7.4模型參數(shù)化方法

7.5自相關(guān)矩陣的本征分析法

7.6小結(jié)

7.1功率譜估計(jì)及其分析方法

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

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

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

(7.6)

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

7.2經(jīng)典譜估計(jì)法

經(jīng)典功率譜估計(jì)法可分成兩種:一種是自相關(guān)函數(shù)估計(jì)法，或稱為BT法;另一種是基于DFT的周期圖法，或稱為直接法。

7.2.1自相關(guān)函數(shù)估計(jì)法

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

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

考慮乘積項(xiàng)的長度，自相關(guān)序列的估計(jì)為(7.9)式中，m取絕對值是因?yàn)棣誼x(m)=φxx(－m)，m為負(fù)值時(shí)上式仍適用。規(guī)定的求和上下限的原則是保持充分利用全部數(shù)據(jù)。

【例7.1】

采用自相關(guān)函數(shù)估計(jì)法，求帶有白噪聲干擾的頻率為10Hz的正弦信號的功率譜。

MATLAB程序如下:

%MATLABPROGRAM7-1

clf;f=10;

N=1000;Fs=500;%數(shù)據(jù)長度和采樣頻率

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

Lag=100;%延遲樣本點(diǎn)數(shù)

randn(′state′，0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)

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

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

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

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

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

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

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

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

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

fp=(0:length(Pxx)－1)′*Fs/length(Pxx);%求功率譜的橫坐標(biāo)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′);

程序運(yùn)行結(jié)果如圖7.1所示。圖7.1自相關(guān)函數(shù)估計(jì)間接法求功率譜密度7.2.2周期圖法

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

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

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

7.3改進(jìn)的非參數(shù)化方法

7.3.1分段平均周期圖法

如果x1，x2，…，xL是不相關(guān)的隨機(jī)變量，每一個(gè)具有期望值μ，方差σ2，則它們的數(shù)學(xué)平均的期望值等于μ，數(shù)學(xué)平均的方差等于σ2/L，即(7.16)(7.17)當(dāng)L→∞時(shí)，Var[

]→0，可達(dá)到一致譜估計(jì)的目的。因而降低估計(jì)量方差的一種有效方法是將若干個(gè)獨(dú)立估計(jì)值進(jìn)行平均。把這種方法應(yīng)用于譜估計(jì)就是Bartlett平均周期圖法。

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

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

(7.18)

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

(7.21)由此得到(7.22)式中，M=N/L，因此Bartlett估計(jì)的期望值是真實(shí)譜Pxx(ω)與三角窗函數(shù)的卷積。由于三角窗函數(shù)不等于δ函數(shù)，所以Bartlett估計(jì)也是有偏估計(jì)，但當(dāng)N→∞時(shí)是無偏估計(jì)。假定各段周期圖是相互獨(dú)立的，可得到周期圖的方差為(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;

程序運(yùn)行結(jié)果如圖7.3所示。圖7.3分段平均周期圖法功率譜估計(jì)7.3.2加窗平均周期圖法

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

【例7.4】

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

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;

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

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

【例7.5】

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

MATLAB程序如下:

%MATLABPROGRAM7-5

clc;

N=256;Fs=1;

u=rand(1，N);

u=u-mean(u);

%用Welch法平均估計(jì)實(shí)驗(yàn)數(shù)據(jù)的功率譜

[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點(diǎn)一次求周期圖的功率譜

[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;

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

MATLAB還可以利用函數(shù)psd和csd來實(shí)現(xiàn)Welch法的功率譜估計(jì)。

(1)函數(shù)psd利用Welch法來

估計(jì)一個(gè)信號的自功率譜，調(diào)用格式為

Pxx=psd(x)

Pxx=psd(x，Nfft)

Pxx=psd(x，Nfft，F(xiàn)s，windows，Noverlap)

[Pxx，F(xiàn)]=psd(x，Nfft，F(xiàn)s，windows，Noverlap，dflag)

psd(x，…，dflag)

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

【例7.6】

求取隨機(jī)信號序列x(n)=sin(2πnf)+e－0.2n+u(n)，n=0，1，…，N－1的功率譜估計(jì)，其中，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(xiàn)]=psd(xn，Nfft，F(xiàn)s，windows，noverlap，dflag);

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

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

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

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

(2)函數(shù)csd利用Welch法估計(jì)兩個(gè)隨機(jī)信號的互功率譜密度，調(diào)用格式為

Pxy=csd(x，y)

Pxy=csd(x，y，Nfft)

Pxy=csd(x，y，Nfft，F(xiàn)s，windows，Noverlap)

[Pxy，F(xiàn)]=csd(x，Nfft，F(xiàn)s，windows，Noverlap，dflag)

csd(x，y，…，dflag)

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

【例7.7】

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

MATLAB程序如下:

%MATLABPROGRAM7-7

clf;

Fs=6000;

%信號采樣頻率

randn(′state′，0);%設(shè)置隨機(jī)狀態(tài)初始值

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

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

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

%產(chǎn)生含周期成分的噪聲信號

[Pxy，f]=csd(x，y，1024，F(xiàn)s，triang(500)，0);

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

[Pxz，f]=csd(x，z，1024，F(xiàn)s，triang(500)，0);

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

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

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

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

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

程序運(yùn)行結(jié)果如圖7.7所示。圖7.7兩個(gè)噪聲信號及噪聲信號和含噪聲周期信號的互功率譜7.3.4多窗口法

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

MATLAB信號處理工具箱提供函數(shù)pmtm實(shí)現(xiàn)MTM法估計(jì)功率譜密度，調(diào)用格式為

Pxx=pmtm(x)

Pxx=pmtm(x，nw，Nfft)

[Pxx，F(xiàn)]=pmtm(x，nw，Nfft，F(xiàn)s)

其中，x為信號序列;nw為時(shí)間帶寬積，缺省值為4，通?？扇?、5/2、3、7/2等;Nfft為FFT的長度;Fs為采樣頻率。

【例7.8】

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

MATLAB程序如下:

%MATLABPROGRAM7-8

clf;

Fs=1000;

f1=60;f2=120;

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

randn(′state′，0);

%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)

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

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

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

%繪制功率譜

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

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

[Pxx2，F(xiàn)]=pmtm(xn，2，Nfft，F(xiàn)s);

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

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

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

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

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

7.4模型參數(shù)化方法

7.4.1AR模型法

在實(shí)際中所遇到的隨機(jī)過程，常?？梢杂靡粋€(gè)具有有理分式的傳遞函數(shù)的模型來表示，因此，可以用一個(gè)線性差分方程作為產(chǎn)生隨機(jī)序列x(n)的系統(tǒng)模型

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

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

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

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

【例7.9】

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

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

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;

%用自相關(guān)法求AR模型的系數(shù)a和輸入白噪聲的功率譜E

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

%根據(jù)定義計(jì)算信號的功率譜密度

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

f=(0:length(Pxx)－1)*Fs/length(Pxx);%計(jì)算各點(diǎn)對應(yīng)的頻率值

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;

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

線性預(yù)測誤差為

(7.40)

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

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

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

e(n)=u(n)

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

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

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

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

(7.54)

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

【例7.10】

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

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(xiàn)]=pmem(xn，14，Nfft，F(xiàn)s);

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

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

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

grid;

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

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

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

遞推公式為

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

akk=0，σk2=σ2p

(7.59)

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

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

(7.60)

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

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

MATLAB程序如下:

%MATLABPROGRAM7-11

clc;

randn(′state′，0);

%設(shè)置隨機(jī)函數(shù)的狀態(tài)

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

%濾波器分母多項(xiàng)式系數(shù)

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

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

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

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

aa=

1.00000.18490.12790.11140.1839

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

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

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

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

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

【例7.12】

用Burg遞推算法估計(jì)例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算法得到功率譜估計(jì)

[xpsd，F(xiàn)]=pburg(x，8，N，F(xiàn)s);%AR模型階數(shù)為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(xiàn)]=pburg(x，14，N，F(xiàn)s);%AR模型階數(shù)為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;程序運(yùn)行結(jié)果如圖7.11所示。由圖可知，Burg遞推算法的功率譜估計(jì)曲線更平滑。當(dāng)AR模型的階數(shù)取值較大時(shí)，功率譜估計(jì)更接近于實(shí)際值。圖7.11Burg算法實(shí)現(xiàn)功率譜估計(jì)

7.5自相關(guān)矩陣的本征分析法

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

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

(7.61)

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

【例7.13】

試用MUSIC法求信號x(t)=2sin(2πf1t)+cos(2πf2t)+u(t)的功率譜估計(jì)。其中，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，F(xiàn)s，32，16);

Pxx=10*log10(Pxx);

plot(f，Pxx);

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

grid;

程序運(yùn)行結(jié)果如圖7.12所示。圖7.12MUSIC法實(shí)現(xiàn)功率譜估計(jì)

【例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法求信號的功率譜估計(jì)。

M