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

第10章隨機(jī)信號(hào)的高階譜分析10.1高階累積量與高階譜

10.2累積量與雙譜的性質(zhì)

10.3雙譜估計(jì)

10.4高階譜分析的應(yīng)用

10.5小結(jié)

10.1高階累積量與高階譜

10.1.1累積量

設(shè)X表示有限階矩的隨機(jī)變量，定義X的矩生成函數(shù)或特征參數(shù)為

(10.1)Φx(ω)=E[exp(jωX)]考慮由k個(gè)實(shí)隨機(jī)變量{xk}組成一個(gè)隨機(jī)向量X=[x1，x2，…，xk]T，設(shè)ω=[ω1，ω2，…，ωk]T，定義{xk}的聯(lián)合特征函數(shù)為(10.2)或(10.3)定義序列{xk}的k階累積量生成函數(shù)為(10.4)

考慮各態(tài)遍歷序列的情況，定義k個(gè)實(shí)隨機(jī)變量(x1，x2，…，xk)的r階聯(lián)合累積量為(10.5)式中，r=r1+r2+…+rk。定義隨機(jī)變量(x1，x2，…，xk)的r階聯(lián)合矩為(10.6)因此，隨機(jī)變量的聯(lián)合累積量可以通過(guò)它的聯(lián)合矩來(lái)表示?？紤]r1=r2=…=rk=1的情況，對(duì)于具有零均值的實(shí)隨機(jī)變量，其二階、三階和四階累積量分別為(10.7)(10.8)(10.9)由此可知，三階以下的累積量和矩是等價(jià)的，但四階以上是不同的。四階以上累積量包含四階矩及互相關(guān)函數(shù)。由累積量確定矩的關(guān)系為

(10.10)(10.11)假定序列{x(t)}為一零均值的k階實(shí)平穩(wěn)隨機(jī)過(guò)程，滿足各態(tài)遍歷性，其k階累積量定義為隨機(jī)變量x(t)，x(t+τ1)，…，x(t+τk－1)的聯(lián)合k階累積量(10.12)

根據(jù)平穩(wěn)性的假設(shè)，可知(10.13)(10.14)

(10.15)為了更明確地說(shuō)明累積量的物理意義，這里介紹累積量的另一種定義。設(shè){x(t)}為一零均值的k階實(shí)平穩(wěn)隨機(jī)過(guò)程，設(shè){g(t)}為一高斯隨機(jī)過(guò)程，若{x(t)}和{g(t)}具有相同的自相關(guān)函數(shù)，則定義{x(t)}的k階累積量為(10.16)

【例10.1】

利用MATLAB產(chǎn)生兩個(gè)遠(yuǎn)程的發(fā)送信號(hào)，并采用歸一化處理使其中一個(gè)信號(hào)有時(shí)延;然后利用MATLAB產(chǎn)生兩個(gè)相互獨(dú)立、互不相關(guān)的噪聲，且噪聲滿足零均值;再將噪聲加入到發(fā)送信號(hào)上，形成兩個(gè)基站的接收信號(hào);最后求兩個(gè)接收信號(hào)的高階累積量。

MATLAB程序如下:

%MATLABPROGRAM10-1

%常量初始化********************************

N=2048;%取2048個(gè)時(shí)間點(diǎn)

w0=pi/4;

%主頻

d1=25;

%時(shí)間延遲，用來(lái)對(duì)信號(hào)作歸一化處理

d2=10;m=50;p=25;

%信號(hào)源**************************

signal1=zeros(1，N);%1行N列的值都為零的矩陣

signal2=zeros(1，N);

st=signalSOI(N);%2048個(gè)值為1或者-1的一維矩陣

fori=1:N

signal1=st*cos(w0*i);

end

fori=1:N－d1－d2

%歸一化處理使信號(hào)2(signal2)在信號(hào)1的時(shí)延基礎(chǔ)上產(chǎn)生

signal2(1，i)=signal1(1，i+d1)+signal1(1，i+d2);

end

%平穩(wěn)隨機(jī)過(guò)程的噪聲**********************

noise1=normrnd(0，0.7，1，N);%產(chǎn)生均值為0，方差為0.7，1行N列的隨機(jī)噪聲

noise2=normrnd(0，0.7，1，N);

%產(chǎn)生均值為0，方差為0.7，1行N列的隨機(jī)噪聲

%形成接收信號(hào)**********************

x=signal1+noise1;

y=signal2+noise2;

%求高階累積量**********************

Cyx=zeros(2*m+1，1);

b=0.25;%cyclefrequency

fort=m:－1:－m

ift>=0

forn=1:N－t

Cyx(m+1－t，1)=Cyx(m+1－t，1)+y(1，n+t)*y(1，n+t)*x(1，n)*x(1，n)*exp(－j*2*pi*b*(n+t/2));

end

Cyx(m+1－t，1)=Cyx(m+1－t，1)/N;

else

forn=－t+1:N

Cyx(m+1－t，1)=Cyx(m+1－t，1)+y(1，n+t)*y(1，n+t)*

x(1，n)*x(1，n)*exp(－j*2*pi*b*(n+t/2));

end

Cyx(m+1－t，1)=Cyx(m+1－t，1)/N;

end

end

Cx=zeros(2*m+1，2*p+1);

forl=m-p:－1:－m－p

fork=m－p:m+p

t=k－(m－p－l);

ift>=0

forn=1:N－t

Cx(m－p－l+1,k+1－(m－p))=Cx(m－p－l+1,k+1－(m－p))

+x(1,n+t)*x(1,n+t)*x(1,n)*x(1,n)*exp(－j*2*pi*b*(n+t/2));

end

Cx(m－p－l+1，k+1－(m－p))=Cx(m－p－l+1，k+1－(m－p))/N;

else

forn=－t+1:N

Cx(m－p－l+1，k+1－(m－p))=Cx(m－p－l+1，k+1－(m－

p))+x(1,n+t)*x(1,n+t)*x(1,n)*x(1,n)*exp(－j*2*pi*b*(n+t/2));

end

Cx(m－p－l+1，k+1－(m－p))=Cx(m－p－l+1，k+1－(m－p))/N;

end

end

end

A=zeros(1，2*p+1);

forl=p:－1:－p

A(1，p+1－l)=exp(－j*pi*b*l);

end

D=diag(A);

B=Cx*D;

A=inv((B′*B))*B′*Cyx;

x=[-p:p];y=A;

plot(x，y);

其中，自定義功能函數(shù)signalSOI完成信號(hào)的原始數(shù)據(jù)的接收，該子程序?yàn)?/p>

functionst=signalSOI(t)

BaudRate=4;%1/Tbisbaudrate

r=rand(1，ceil(t/BaudRate));

fori=1:ceil(t/BaudRate)

ifr(i)>=0.5judge(i)=1;

else judge(i)=－1;

end

end

fori=1:t/BaudRate

forn=1:BaudRate

st(BaudRate*(i－1)+n)=judge(i);

end

end

程序運(yùn)行結(jié)果如圖10.1所示。圖10.1高階累積量估計(jì)曲線10.1.2高階譜

設(shè)實(shí)隨機(jī)序列{xk}的均值為0，k階平穩(wěn)，其k階累積量Ck,x(τ1，τ2，…，τk－1)存在，則定義{xk}的k階譜為(10.17)式中，ω=[ω1,ω2,…,ωk－1]T，τ=[τ1,τ2,…,τk－1]T。

k階譜又稱為多譜或累量譜，一般情況下為復(fù)數(shù)，其存在的充分條件之一是k階累積量滿足絕對(duì)可積。不難看出，維納-辛欽定理是k=2時(shí)的特例。在高階譜分析中，k=3時(shí)的高階譜稱為雙譜(Bispectrum)，其應(yīng)用非常廣泛。

10.2累積量與雙譜的性質(zhì)

10.2.1累積量的性質(zhì)

(1)若λi(1，2，…，k)為常數(shù)，{xk}為隨機(jī)變量，則有

(10.18)(2)累積量對(duì)所有變量對(duì)稱，即(10.19)式中，(i，j，…，p)是(1，2，…，k)的一個(gè)排列。

(3)累積量對(duì)變量具有可加性。設(shè)x0，y0是兩個(gè)不同的隨機(jī)變量，則(10.20)

(4)如果λ為常數(shù)，則(10.21)

(5)隨機(jī)變量{xk}與{yk}彼此統(tǒng)計(jì)獨(dú)立，則有

cum(x1+y1,x2+y2,…,xk+yk)=cum(x1,x2,…,xk)+cum(y1,y2,…,yk)(10.22)利用該性質(zhì)可以進(jìn)行高斯背景中的非高斯信號(hào)檢測(cè)。

(6)如果隨機(jī)變量{xk}中的一部分與其他部分統(tǒng)計(jì)獨(dú)立，則有

cum(x1，x2，…，xk)=0

(10.23)這說(shuō)明獨(dú)立同分布的非高斯隨機(jī)過(guò)程的累積量是一種沖激函數(shù)形式。10.2.2雙譜及其性質(zhì)

1.雙譜的定義

雙譜(Bispectrum)是高階譜在R=3時(shí)的特例，采用符號(hào)Bx(ω1，ω2)表示，為一個(gè)二維復(fù)數(shù)。(10.24)雙譜是三階自相關(guān)函數(shù)Rx(m，n)的二維傅里葉變換。設(shè)隨機(jī)序列{x(R)}的三階自相關(guān)函數(shù)Rx(m，n)為(10.26)Rx(m，n)=E[x(R)x(R+m)x(R+n)](10.25)則隨機(jī)序列的雙譜Bx(ω1，ω2)為

2.雙譜的性質(zhì)

設(shè)實(shí)平穩(wěn)隨機(jī)系列{x(k)}均值為0，根據(jù)隨機(jī)序列三階自相關(guān)定義可得以下對(duì)稱關(guān)系:

Rx(n,m)=Rx(m,n)=Rx(－m,n－m)=Rx(n－m,－m)

=Rx(－n,m－n)=Rx(m－n，－n)

(10.27)

因此，三階自相關(guān)函數(shù)在(m，n)平面內(nèi)有六種互換可能，即只要知道m(xù)、n平面上第一象限內(nèi)m＝0、n＝m兩條直線構(gòu)成的區(qū)域內(nèi)的值，便可利用對(duì)稱性質(zhì)得到全面的值。

根據(jù)雙譜的定義，可推出如下性質(zhì):

(1)雙譜為復(fù)數(shù)，可得其幅度譜和相位譜為

Bx(ω1,ω2)=|Bx(ω1,ω2)|exp[jjB(ω1,ω2)]

(10.28)

(2)雙譜是以2π為周期的雙周期函數(shù)，即

Bx(ω1,ω2)=Bx(ω1+2π,ω2+2π)

(10.29)

(3)雙譜具有的對(duì)稱性質(zhì)為

Bx(ω1,ω2)=Bx(ω2,ω1)=B*x(－ω2,－ω1)=B*x(－ω1,－ω2)

=Bx(－ω1,ω2,ω2)=Bx(ω1,－ω1－ω2)

=Bx(－ω1,ω2,ω1)=Bx(ω2,－ω1,ω2)(10.30)

(4)對(duì)于持續(xù)時(shí)間有限的隨機(jī)序列{x(k)}，如果其傅里葉變換X(ω)存在，則雙譜可由下式確定:

Bx(ω1,ω2)=X(ω1)X(ω2)X*(ω1+ω2)

=X(ω1)X(ω2)X(－ω1－ω2)

(10.31)

(5)三階平穩(wěn)0均值非高斯白噪聲序列的功率譜和雙譜均為常數(shù)。

(6)高斯過(guò)程的雙譜恒為0。

10.3雙譜估計(jì)

10.3.1非參數(shù)化雙譜估計(jì)

1.間接法雙譜估計(jì)

設(shè)觀測(cè)數(shù)據(jù){x(1)，x(2)，…，x(N)}為一實(shí)隨機(jī)序列。間接法的核心思想是，由三階累積量定義估計(jì)得到它的三階自相關(guān)函數(shù)，然后對(duì)其進(jìn)行二維DFT變換得到隨機(jī)序列的雙譜估計(jì)。具體算法歸納如下:

(1)有限長(zhǎng)觀測(cè)數(shù)據(jù){x(k)}(k=1，2，…，N)分成K段，每段數(shù)據(jù)有M個(gè)點(diǎn)，即N＝KM;也可以重疊分段，使相鄰數(shù)據(jù)段有一半重疊，即N＝2KM。

(2)除去每段數(shù)據(jù)的均值。

(3)第i段數(shù)據(jù)記為{xi(k)}(k=1，2，…，M;i=1，2，…，K)。計(jì)算每段數(shù)據(jù)的三階累積量:

(10.32)式中，k1=max{0，－m，－n}，k2=min{M，M－m，M－n}。

(4)對(duì)K組進(jìn)行統(tǒng)計(jì)平均，即(10.33)

(5)對(duì)三階累積量進(jìn)行二維傅里葉變換，得雙譜估計(jì)為(10.34)對(duì)于MATLAB6.x和MATLAB7.x以及MATLABR2008版本來(lái)說(shuō)，此函數(shù)不存在。利用MATLAB編寫間接法實(shí)現(xiàn)雙譜估計(jì)的函數(shù)，程序如下:

function[Bspec，waxis]=bispeci(y，nlag，nsamp，overlap，flag，nfft，wind)

%----------parameterchecks-----------------

[ly，nrecs]=size(y);

if(ly==1)y=y(:);ly=nrecs;nrecs=1;end

if(exist(′overlap′)~=1)overlap=0;end

overlap=min(99，max(overlap，0));

if(nrecs>1)overlap=0;end

if(exist(′nsamp′)~=1)nsamp=ly;

end

if(nsamp>ly|nsamp<=0)nsamp=ly;

end

if(exist(′flag′)~=1)flag=′biased′;end

if(flag(1:1)~=′b′)flag=′unbiased′;end

if(exist(′nfft′)~=1)nfft=128;end

if(nfft<=0)nfft=128;end

if(exist(′wind′)~=1)wind=0;end

nlag=min(nlag，nsamp－1);

if(nfft<2*nlag+1)nfft=2^nextpow2(nsamp);end

%------------createthelagwindow------------

Bspec=zeros(nfft，nfft);

if(wind==0)

indx=(1:nlag)′;

window=[1;sin(pi*indx/nlag)./(pi*indx/nlag)];

else

window=ones(nlag+1，1);

end

window=[window;zeros(nlag，1)];

%-----------cumulantsinnon-redundantregion-----------

%definecum(i，j)=Econj(x(n))x(n+i)x(n+j)

%foracomplexprocess，weonlyhavecum(i，j)=cum(j，i)

overlap=fix(nsamp*overlap/100);

nadvance=nsamp-overlap;

nrecord=fix((ly*nrecs-overlap)/nadvance);

c3=zeros(nlag+1，nlag+1);

ind=[1:nsamp]′;

fork=1:nrecord，

x=y(ind);x=x-mean(x);

ind=ind+nadvance;

forj=0:nlag

z=x(1:nsamp-j).*x(j+1:nsamp);

fori=j:nlag

sum=z(1:nsamp－i)′*x(i+1:nsamp);

if(flag(1:1)==′b′)，sum=sum/nsamp;

else

sum=sum/(nsamp－i);

end

c3(i+1，j+1)=c3(i+1，j+1)+sum;

end

end

end

c3=c3/nrecord;

%cumulantselsewherebysymmetry----------------

c3=c3+tril(c3，－1)′;%completeIquadrant

c31=c3(2:nlag+1，2:nlag+1);

c32=zeros(nlag，nlag);c33=c32;c34=c32;

fori=1:nlag，

x=c31(i:nlag，i);

c32(nlag+1－i，1:nlag+1－i)=x′;

c34(1:nlag+1－i，nlag+1－i)=x;

if(i<nlag)

x=flipud(x(2:length(x)));

c33=c33+diag(x，i)+diag(x，-i);

end

end

c33=c33+diag(c3(1，nlag+1:－1:2));

cmat=[[c33，c32，zeros(nlag，1)];[[c34;zeros(1，nlag)]，c3]];

%---------applylag-domainwindow------------------

wcmat=cmat;

if(wind~=－1)

indx=[-nlag:nlag]′;

fork=-nlag:nlag

wcmat(:，k+nlag+1)=cmat(:，k+nlag+1)...

*window(abs(indx－k)+1).*window(abs(indx)+1)...

*window(abs(k)+1);

end

end

%-----compute2d-fft，andshiftandrotateforproperorientation------

Bspec=fft2(wcmat，nfft，nfft);

Bspec=fftshift(Bspec);%axesdandr;origatctr

if(rem(nfft，2)==0)

waxis=[－nfft/2:(nfft/2－1)]/nfft;

else

waxis=[－(nfft－1)/2:(nfft－1)/2]/nfft;

end

contour(waxis，waxis，abs(Bspec)，4)，grid;

title(′Bispectrumestimatedviatheindirectmethod′);

xlabel(′f1′);ylabel(′f2′);

set(gcf，′Name′，′HosaBISPECI′)

return

【例10.2】

有一正弦輸入信號(hào)s(t)=0.8sin(2πft)，f=

68Hz，一個(gè)高斯白噪聲信號(hào)ω(t)疊加到正弦信號(hào)中。采樣頻率fs=5kHz，采樣點(diǎn)數(shù)為5000。試用間接法進(jìn)行雙譜估計(jì)。

MATLAB程序如下:

%MATLABPROGRAM10-2

clc;

fs=5000;

%采樣頻率

t=[0:1/fs:1.0－1/fs];N=1.0*fs;%采樣點(diǎn)數(shù)

Am=0.8;

%輸入信號(hào)幅度

s=Am*sin(2*pi*68*t);

%輸入信號(hào)，頻率68Hz

A=1;

noi=A*randn(1，N);

%高斯白噪聲

d=s+noi;

SNR=10*log10(Am*Am);

%信噪比

SNR

%間接法雙譜估計(jì)

figure(1);

subplot(221);bspec1=bispeci(noi，21，64，0，′unbiased′，128，1);title(′bispectrumofnoi′);

subplot(222);[bspec2,waxis2]=bispeci(d,21,64,0,′unbiased′,128,1);title(′bispectrumofd′);

figure(2);

subplot(121);mesh(abs(bspec1));title(′noi′);%網(wǎng)格圖

subplot(122);mesh(abs(bspec2));title(′d′);

figure(3);

subplot(121);mesh(abs(bspec1));title(′noi′);

view(－37.5，0);

subplot(122);mesh(abs(bspec2));title(′d′);

view(－37.5，0);

程序運(yùn)行結(jié)果如下:

SNR=－1.9382

雙譜估計(jì)曲線如圖10.2所示。圖10.2間接法實(shí)現(xiàn)雙譜估計(jì)

2.直接法雙譜估計(jì)

設(shè)實(shí)平穩(wěn)隨機(jī)過(guò)程x(t)的傅里葉變換X(ω)存在，則有(10.35)因此(10.36)

【例10.3】

設(shè)信號(hào)，其中，f1=0.1，j1=π/6;f2=0.15，j2=π/3;f3=0.25，j3=π/2;f4=0.4，j4=π/4。f3分量是f1分量和f2分量的相位耦合。用直接法雙譜估計(jì)檢測(cè)信號(hào)平方相位耦合信息。

MATLAB程序如下:

%MATLABPROGRAM10-3

N=64;

n=[0:N－1];

f1=0.1;f2=0.15;f3=0.25;f4=0.4;

fai1=pi/6;

fai2=pi/3;

fai3=pi/2;

fai4=pi/4;

s1=cos(2*pi*f1*n+fai1);

s2=cos(2*pi*f2*n+fai2);

s3=cos(2*pi*f3*n+fai3);

s4=cos(2*pi*f4*n+fai4);

x=s1+s2+s3+s4;

bspec=bispecd(x，10);

mesh(abs(bspec));

程序運(yùn)行結(jié)果如圖10.3所示。圖10.3直接法實(shí)現(xiàn)雙譜估計(jì)

3.二維窗函數(shù)

類似于對(duì)有限長(zhǎng)隨機(jī)序列進(jìn)行功率譜估計(jì)，無(wú)論采用間接法還是直接法，為了提高功率譜估計(jì)的質(zhì)量，減少估計(jì)方差，都可以使用適當(dāng)?shù)拇昂瘮?shù)。窗函數(shù)不僅與其寬度有關(guān)，還與其形狀有關(guān)，對(duì)于二維情況，其原理類似，完全可以類推，這里不再贅述。10.3.2參數(shù)化雙譜估計(jì)

1.非高斯MA模型

設(shè)隨機(jī)序列{x(k)}為一q階MA過(guò)程:(10.40)式中，w(k)為一非高斯白噪聲序列，均值為0，其方差為E[w2(k)]=Q，并滿足E[w3(k)]=β≠0，而且x(k)與w(k)統(tǒng)計(jì)獨(dú)立。對(duì)于以上非高斯MA(q)模型，需要根據(jù)有限長(zhǎng)觀測(cè)數(shù)據(jù)確定各模型參數(shù)。一般假定b(0)=1，考慮系統(tǒng)疊加的背景噪聲序列為n(k)，于是有

y(k)=x(k)+n(k)

(10.41)

y(k)的三階累積量為

(10.42)

當(dāng)k=0，m=q，n=i時(shí)，有(10.43)(10.44)

根據(jù)以上分析，可得MA模型的參數(shù)為(10.45)以上稱為非高斯MA模型參數(shù)估計(jì)的閉式解。一旦確定了q個(gè)參數(shù)b(k)，MA(q)模型的雙譜估計(jì)為(10.46)(10.47)

1)函數(shù)maorder

功能:MA模型的定階。

格式:q＝maorder(y，qmin，qman，pfa，flag)

說(shuō)明:maorder返回由序列y估計(jì)的MA階數(shù)q;qmin和qmax分別為預(yù)期的最小和最大階數(shù)，缺省值分別是0和10;pfa為允許的出錯(cuò)概率，缺省值是0.05;flag為非零值時(shí)，函數(shù)顯示算法中相關(guān)統(tǒng)計(jì)量的值，其缺省值為1。

2)函數(shù)maest

功能:MA模型的參數(shù)估計(jì)階。

格式:bvec＝maest(y，q)

bvec＝maest(y,q,norder,samp_seg,overlap,flag)

說(shuō)明:bvec＝maest(y，q)返回由序列y估計(jì)的MA模型參數(shù)bvec，q為指定的模型階數(shù)。norder指定算法中采用何種累積量，可選值為3和4，缺省值為3;samp_seg為每個(gè)數(shù)據(jù)分段中的數(shù)據(jù)點(diǎn)數(shù)，缺省值為序列的長(zhǎng)度;overlap為每個(gè)數(shù)據(jù)段間重疊的點(diǎn)數(shù)(0~99)，缺省值為0;若flag的第一個(gè)字符為′b′，則采用有偏估計(jì)，否則采用無(wú)偏估計(jì)，缺省時(shí)采用有偏估計(jì)。對(duì)于MATLAB6.x和MATLAB7.x以及MATLABR2008版本，函數(shù)maorder和maest均不存在。利用MATLAB編程實(shí)現(xiàn)MA模型的參數(shù)估計(jì)，函數(shù)maest的程序如下:

functionbvec=maest(y，q，norder，samp_seg，overlap，flag)

%------parameterchecks-----------

if(nargin<2)

error(′insufficientnumberofparameters′)

end

[nsamp，nrecs]=size(y);

if(nsamp==1)nsamp=nrecs;nrecs=1;y=y.′;end

if(q<=0)bvec=1;return，end

if(exist(′norder′)~=1)norder=3;end

if(norder~=3&norder~=4)

error(′cumulantordermustbe3or4′)

end

if(exist(′samp_seg′)~=1)samp_seg=nsamp;end

if(exist(′overlap′)~=1)overlap=0;end

overlap=max(0，min(overlap，99));

if(exist(′flag′)~=1)flag=′biased′;end

if(nrecs>1)samp_seg=nsamp;overlap=0;end

%--------estimatecumulantsandcovariances--------

c2=cumest(y，2，q，samp_seg，overlap，flag);

c2=[c2;zeros(q，1)];

cumd=cumest(y，norder，q，samp_seg，overlap，flag，0，0);

cumq=cumest(y，norder，q，samp_seg，overlap，flag，q，0);

cumd=cumd(2*q+1:－1:1);

cumd=[cumd;zeros(q，1)];

cumq(1:q)=zeros(q，1);

%--------TheGM-RCLSalgorithm-----------

cmat=toeplitz(cumd，[cumd(1)，zeros(1，q)]);

rmat=toeplitz(c2，[c2(1)，zeros(1，q)]);

amat0=[cmat，-rmat(:，2:q+1)];

rvec0=c2;

%-------TheTugnaitfix----------------

cumq=[cumq(2*q+1:－1:q+1);zeros(q，1)];

cmat4=toeplitz(cumq，[cumq(1)，zeros(1，q)]);

c3=cumd(1:2*q+1);

amat0=[amat0，zeros(3*q+1，1);zeros(2*q+1，q+1)，cmat4(:，2:q+1)，－c3];

rvec0=[rvec0;－cmat4(:，1)];

%--------GetridofR(0)terms-------------

row_sel=[1:q，2*q+2:3*q+1，3*q+2:4*q+1，4*q+3:5*q+2];

amat0=amat0(row_sel，:);

rvec0=rvec0(row_sel);

%--------SolveforMAparms----------------

bvec=amat0＼rvec0;

b1=bvec(2:q+1)/bvec(1);

b2=bvec(q+2:2*q+1);

if(norder==3)

if(all(b2>0))

b1=sign(b1).*sqrt(0.5*(b1.^2+b2));

else

disp(′MAEST:alternativesolutionb1used′)

end

else

if(sign(b2)==sign(b1))

b1=sign(b1).*(abs(b1)+abs(b2).^(1/3))/2;

else

disp(′MAEST:alternativesolutionb1used′)

end

end

bvec=[1;b1];

return

2.非高斯AR模型

AR模型雙譜估計(jì)方法用于估計(jì)非高斯白噪聲通過(guò)p階AR模型而產(chǎn)生的輸出序列的雙譜?？疾靝階AR過(guò)程{x(k)}，滿足以下關(guān)系:

(10.48)式中，w(k)為一非高斯白噪聲序列，均值為0，其方差E[w2(k)]=0，并滿足E[w3(k)]=β≠0，而且x(k)與w(k)統(tǒng)計(jì)獨(dú)立。如果AR模型為一穩(wěn)定系統(tǒng)，并設(shè)非高斯白噪聲w(k)三階平穩(wěn)，那么x(k)也是三階平穩(wěn)序列。對(duì)式(10.48)求三階自相關(guān)，可得(10.49)式中，δ(n，m)類似于單位沖激函數(shù)，稱為二維單位沖激函數(shù)，即

(10.50)通常稱式(10.49)為三階遞推方程。如果取m=n，則可得{x(k)}的(2p+1)個(gè)三階自相關(guān)R(n，m)的對(duì)角切片值。令n，m=1，2，…，p，可得以下矩陣方程:

R·a=b

(10.51)式中，(10.52)式(10.51)存在的一個(gè)基本條件是AR濾波器的轉(zhuǎn)移函數(shù)滿足穩(wěn)定條件，即(10.53)在滿足條件的情況下，p階穩(wěn)定的AR模型參數(shù)ai(i=1,2,…,p)可由2p+1個(gè)R(n，m)的對(duì)角切片值求得。若給定了2p+1個(gè)對(duì)角切片值R(－p,－p),…,R(0,0),…,R(p，p)，則上述方程便可用來(lái)擬合一個(gè)p階AR模型。這一模型是在AR模型的輸出三階矩序列與其相應(yīng)給定點(diǎn)的采樣值之間的一種良好匹配關(guān)系的意義上被擬合的。若2p+1個(gè)R(n，m)的對(duì)角切片值由一真正的p階過(guò)程的三階矩求得，則模型參數(shù)a1，a2，…，ap隱含了序列{x(k)}的三階矩信息。

1)函數(shù)arorder

功能:AR模型的定階。

格式:p＝arorder(y，norder，pmax，qmax，flag)

說(shuō)明:arorder返回由序列y估計(jì)的AR階數(shù)p;

norder指定算法中采用何種累積量，可選值為2，±3，±4，負(fù)號(hào)意味著同時(shí)使用相關(guān)，缺省值為3;pmax和qmax分別為預(yù)期的最大AR階數(shù)和最大MA階數(shù)，缺省值都是10;flag＝1時(shí)，函數(shù)以波形方式顯示出累積量矩陣的奇異值，其中的峰值對(duì)應(yīng)階數(shù)估計(jì)的結(jié)果由用戶自行確定，否則，函數(shù)直接返回估計(jì)的階數(shù)，flag的缺省值為1。

2)函數(shù)arrcest

功能:AR模型的參數(shù)估計(jì)。

格式:avec＝arrcest(y，p)

avec＝arrcest(y,p,minlag,norder,maxlag,samp_seg,overlap,flag)

說(shuō)明:avec＝arrcest(y，p)由序列y和預(yù)估的模型階數(shù)，估計(jì)AR模型的參數(shù)向量avec;minlag是標(biāo)準(zhǔn)方程中參數(shù)的最小值，如果估計(jì)的模型為ARMA，則minlag應(yīng)大于q;norder指定算法中采用何種累積量，可選值為2，±3，±4，負(fù)號(hào)意味著同時(shí)使用相關(guān)，缺省值為2;maxlag為所使用的最大的累積量延遲，缺省值為p+minlag;samp_seg為各數(shù)據(jù)分段的數(shù)據(jù)點(diǎn)數(shù)，缺省值為整個(gè)時(shí)間序列的長(zhǎng)度;overlap為各數(shù)據(jù)分段間重疊的數(shù)據(jù)點(diǎn)數(shù)，缺省值為0;如果flag的第一個(gè)字符為′b′，則采用有偏估計(jì)，否則采用無(wú)偏估計(jì)，缺省時(shí)采用有偏估計(jì)。對(duì)于MATLAB6.x和MATLAB7.x以及MATLABR2008版本，函數(shù)arorder和arrcest均不存在。利用MATLAB編程可實(shí)現(xiàn)AR模型的參數(shù)估計(jì)，函數(shù)arrcest的程序如下:

functionavec=arrcest(y，p，q，norder，maxlag，samp_seg，overlap，flag)

%---------------parameterchecks--------------

if(nargin<2)

error(′insufficientnumberofparameters′)

end

[nsamp，nrecs]=size(y);

if(nsamp==1)nsamp=nrecs;nrecs=1;y=y.′;end

if(p<=0)avec=1;return，end

if(exist(′q′)~=1)q=0;end

if(q<0)

error(′MAorderqmustbenon-negative′)，

end

if(exist(′norder′)~=1)norder=2;end

if(norder~=2&abs(norder)~=3&abs(norder)~=4)

error(′nordermustbe2，3，4，-3or-4′)

end

if(exist(′maxlag′)~=1)

maxlag=p+q;

end

if(maxlag<p+q)，

disp([′ARRCEST:maxlagchangedfrom′，int2str(maxlag)，...′to′，int2str(p+q)])

maxlag=p+q;

end

if(exist(′samp_seg′)~=1)samp_seg=nsamp;end

if(exist(′overlap′)~=1)overlap=0;end

overlap=max(0，min(overlap，99));

if(exist(′flag′)~=1)flag=′biased′;end

if(nrecs>1)overlap=0;samp_seg=nsamp;end

minlag=－maxlag;

nlags=maxlag-minlag+1;

Amat=[];

rvec=[];

%------------estimatecumulants------------

if(norder~=2)

kslice1=[q－p，q];

kslice2=[0，0];

kslice=(kslice1(2)-kslice1(1)+1)*(kslice2(2)-kslice2(1)+1);

cum_y=zeros(nlags，kslice);

morder=abs(norder);

kloc=0;

fork1=kslice1(1):kslice1(2)

fork2=kslice2(1):kslice2(2)

kloc=kloc+1;

cum_y(:，kloc)=cumest(y，morder，maxlag，

samp_seg，overlap，flag，k1，k2);

end

end

%------setupcumulant-based′normal′equations----

Amat=hankel(cum_y(q－p+1－minlag+1:nlags－p，1)，...cum_y(nlags－p:nlags－1，1));

rvec=cum_y(q+1－minlag+1:nlags，1);

fork=2:kslice

Amat=[Amat;hankel(cum_y(q－p+1－minlag+1:nlags－p，k)，cum_y(nlags－p:nlags－1，k))];

rvec=[rvec;cum_y(q+1－minlag+1:nlags，k)];

end

rvec=－rvec;

end

%-----appendcorrelation-basednormalequations-------

if(norder==2|norder<0)

cor_y=cumest(y，2，maxlag，samp_seg，overlap，flag);

AR=hankel(cor_y(q－p+1－minlag+1:nlags－p，1)，...cor_y(nlags－p:nlags－1，1));

br=－cor_y(q+1－minlag+1:nlags，1);

Amat=[Amat;AR];rvec=[rvec;br];

end

%-------computeLSestimate--------

avec=Amat＼rvec;

avec=[1;avec(p:－1:1)];

return

3.非高斯ARMA模型

設(shè)平穩(wěn)隨機(jī)序列x(k)由以下差分方程確定:

(10.60)式中，a0=1，w(k)為零均值平穩(wěn)非高斯隨機(jī)序列。假定ARMA模型滿足因果性、平穩(wěn)性和非最小相位，則該模型H(z)的部分零點(diǎn)可能位于單位圓外。又因

B(ω1，ω2)=βH(ω1)H(ω2)H*(ω1+ω2)(10.61)式中，(10.62)為了根據(jù)觀測(cè)數(shù)據(jù)確定模型參數(shù)ai和bi，據(jù)前兩式可得

(10.63)故(10.64)(10.65)

對(duì)于m或n>q，則(10.66)因此，只要求得三階累積量，便可確定上式中的各系數(shù)g(i，j)。又因(10.67)可得

(10.68)這樣，利用ak可求出G(ω1，ω2)，進(jìn)一步可得(10.69)式中，(10.70)類似地有(10.71)

1)函數(shù)armarts

功能:采用殘余時(shí)間序列法估計(jì)AR模型參數(shù)。

格式:[avec，bvec]＝armarts(y，p，q)

[avec,bvec]＝armarts(y,p,q,norder,maxlag,samp_seg,overlap,flag)

說(shuō)明:[avec，bvec]＝armarts(y，p，q)采用殘余時(shí)間序列法，由序列y估計(jì)出指定AR階數(shù)p和MA階數(shù)q的ARMA模型參數(shù)，avec為AR參數(shù)向量，bvec為MA參數(shù)向量。norder為算法中采用累積量的階數(shù)，可選值為±3，±4，負(fù)號(hào)意味著同時(shí)使用自相關(guān)函數(shù);maxlag為積累量計(jì)算的最大滯后，缺省值為p+q;samp_seg為每個(gè)數(shù)據(jù)分段的數(shù)據(jù)點(diǎn)數(shù)，缺省值為序列的長(zhǎng)度;overlap為每個(gè)數(shù)據(jù)段間重疊的數(shù)據(jù)點(diǎn)數(shù)(0~99)，缺省值為0;當(dāng)y為矩陣時(shí)，其每一列被視為一段觀測(cè)數(shù)據(jù)，此時(shí)samp_seg等于矩陣的行數(shù)，overlap等于0;當(dāng)flag＝′biased′時(shí)，函數(shù)返回有偏估計(jì)，否則返回?zé)o偏估計(jì);缺省時(shí)為有偏估計(jì)。

2)函數(shù)armaqs

功能:采用q切片算法估計(jì)ARMA模型參數(shù)。

格式:[avec，bvec]=armaqs(y，p，q)

[avec，bvec]=armaqs(y,p,q,norder,maxlag,

samp_seg,overlap,flag)

說(shuō)明:[avec，bvec]=armaqs(y，p，q)采用q切片算法由序列y估計(jì)出指定AR階數(shù)p和MA階數(shù)q的ARMA模型參數(shù)，avec為AR參數(shù)向量，bvec為MA參數(shù)向量。norder為算法中采用累積量的階數(shù)，可選值為±3，±4，負(fù)號(hào)意味著同時(shí)使用自相關(guān)函數(shù)。其他參數(shù)同函數(shù)armarts。

3)函數(shù)bispect

功能:基于ARMA模型計(jì)算雙譜。

格式:[bspec，waxis]＝bispect(ma，ar，nfft)

說(shuō)明:bispect函數(shù)返回ARMA序列的雙譜bspec;waxis為相應(yīng)的頻率點(diǎn)向量;ma為MA參數(shù)向量;ar為AR參數(shù)向量，缺省值為[1.00];nfft為FFT的計(jì)算點(diǎn)數(shù)，缺省值為512。

對(duì)于MATLAB6.x和MATLAB7.x以及MATLABR2008版本，函數(shù)armarts、armaqs和bispect均不存在。利用MATLAB編程可實(shí)現(xiàn)ARMA模型的參數(shù)估計(jì)，函數(shù)bispect的程序如下:

function[bspec，waxis]=bispect(ma，ar，nfft)

if(exist(′ma′)~=1)

error(′insufficientnumberofparameters′)

end

q=length(ma)-1;

if(exist(′ar′)~=1)ar=1;end

p=length(ar)－1;

pq=max(p，q)+1;

if(exist(′nfft′)~=1)nfft=max(2^nextpow2(pq)，512);end

nfft=max(nfft，2^nextpow2(pq));

Xf=freqz(ma，ar，nfft，′whole′);

Xfc=conj(Xf);

bspec=(Xf*Xf′).*hankel(Xfc，Xfc([nfft，1:nfft-1]));

bspec=fftshift(bspec);%centertheorigin;normalaxes

if(rem(nfft，2)==0)

waxis=[－nfft/2:nfft/2－1]′/nfft;

else

waxis=[－(nfft－1)/2:(nfft－1)/2]′/nfft;

end

contour(waxis，waxis，abs(bspec)，6)，gridon

xlabel(′f1′)，ylabel(′f2′)

x=[112/311/30－1－1－2/3－1－1/30];

y=[10－1/3－1－2/3－1－101/312/31];

holdon

fork=1:12

plot([0x(k)]，[0y(k)]，′－－′)

end

holdoff

set(gcf，′Name′，′HosaBISPECT′)

10.4高階譜分析的應(yīng)用

10.4.1利用雙譜進(jìn)行時(shí)延估計(jì)

在雷達(dá)、聲納等應(yīng)用中，時(shí)延估計(jì)是一個(gè)重要課題。設(shè)兩個(gè)不同的接收點(diǎn)位置分別為R1和R2，記錄分別為

x(t)=s(t)+n1(t)

y(t)=s(t－T)+n2(t)

(10.72)

式中，n1(t)和n2(t)均為高斯噪聲，信號(hào)與噪聲統(tǒng)計(jì)獨(dú)立，s(t)為目標(biāo)信號(hào)，T為兩個(gè)不同接收點(diǎn)之間的時(shí)延。采用互相關(guān)法對(duì)時(shí)延參數(shù)T進(jìn)行估計(jì)，可得

Rxy(τ)=E[x(t)y(t+τ)]=Rs(τ－T)+Rn(τ)(10.73)

當(dāng)τ=T時(shí)，Rs(τ)取得最大值，因此通過(guò)確定Rs(τ)的峰值位置便可估計(jì)時(shí)延參數(shù)T。但如果信噪比較小，則估計(jì)質(zhì)量將會(huì)變得很差。為此，采用三階統(tǒng)計(jì)量分析同樣的問(wèn)題，對(duì)接收信號(hào)離散化，可得三階自相關(guān)和三階互相關(guān)函數(shù):

Rxx(m，n)=E[x(k)x(k+n)x(k+m)]

(10.74)

Rxyx(m，n)=E[x(k)y(k+n)x(k+m)]

(10.75)

由此得

Rxx(m，n)=Rss(n，m)

(10.76)

Rxyx(m，n)=Rss(n－T，m)

(10.77)其雙譜分別為

Bxx(ω1,ω2)=Bss(ω1,ω2)

(10.78)

Bxyx(ω1,ω2)=Bss(ω1,ω2)ejωT

(10.79)

于是，(10.80)

參數(shù)T可由下式確定:(10.81)所以，通過(guò)估計(jì)序列a(n)出現(xiàn)沖激的時(shí)間便可估計(jì)時(shí)延參數(shù)。采用雙譜技術(shù)進(jìn)行時(shí)延估計(jì)的明顯優(yōu)點(diǎn)是對(duì)背景噪聲不敏感。對(duì)于MATLAB6.x和MATLAB7.x以及MATLABR2008版本，函數(shù)tdeb不存在。利用MATLAB編程可實(shí)現(xiàn)信號(hào)的時(shí)延估計(jì)，程序如下:

function[