6.2離散序列的相關(guān)函數(shù)

第6章功率譜估計(jì)

離散隨機(jī)序列的特征描述

平穩(wěn)隨機(jī)序列通過LTI系統(tǒng)

經(jīng)典功率譜估計(jì)

■現(xiàn)代功率譜估計(jì)

6.1離散隨機(jī)序列的特征描述

隨機(jī)過程的分布函數(shù)

隨機(jī)信號的數(shù)字特征

平穩(wěn)各態(tài)遍歷隨機(jī)信號的時(shí)域描述

平穩(wěn)各態(tài)遍歷隨機(jī)信號的頻域描述(功率譜密度)

一、隨機(jī)過程的分布函數(shù)

{X[k],keZ}表示一個(gè)隨機(jī)過程一維分布函數(shù)F(x,k)=P(X[k]<x)

二維分布函數(shù)

F(x,x;k,k)=P(X[k]<x,X[k]<x)

N維分布函1數(shù)2121122

自相關(guān)函數(shù)

F(xi,x2,…,x；ki，k2，"n)=P(X[k]<x，…X[k]<x二、隨機(jī)信號的數(shù)字1特征1NN

m[k]二E{X[k]}

方差

x

均值

b2[k]二E{(X[k]-m[k])2}二E{X2[k]}-m2[k]

x

x

R[k,k]二E{X[k]X[k]}

互相關(guān)函數(shù)x1212

R[k,k]二E{X[k]Y[k]}

三、平穩(wěn)各態(tài)遍歷隨機(jī)信號的時(shí)域描述

1平穩(wěn)隨機(jī)序列

指統(tǒng)計(jì)特性不隨時(shí)間的平移而變化的那一類隨機(jī)序列

嚴(yán)平穩(wěn)隨機(jī)序列：

+n)

F(X1,3…XN；k1,k2，?…J=F(X1,3…XN；k1+仏k2+冷…N

寬平穩(wěn)隨機(jī)序列：

E{X[k]}二m

x

E{X[k]X[k+n]}二R[n]

x

平穩(wěn)隨機(jī)信號自相關(guān)函數(shù)特性

對稱性

R[n]二R[-n]R*[n]=R[-n]

xxxx

極限值

n二0R[0]=E{X2[k]}

x

nTaR[g]=m2

xx

(3)不等式

R[0]>R[n]

xx

各態(tài)遍歷隨機(jī)信號

集平均等于時(shí)間平均

二E{X[k]}二lim

NTa

￡x[k]

k=-N

b2二E{[X[k]-m]2}=lim1￡[x[k]-m]2

xxnTa2N+1稈x

k=-N

R[n]==E{X[k]X[k+n]}=lim￡x[k]x[k+n]

XNTa2N+1

k=-N

四、平穩(wěn)各態(tài)遍歷隨機(jī)信號的頻域描述

功率譜密度

lim

NTa

E[

'Fx

維納——辛欽公式

P(0)=￡R(n)e-j0n

n=—a

1.

R(n)=FP(0)ej0nd0

2兀-兀

當(dāng)自相關(guān)函數(shù)絕對可積時(shí)，平穩(wěn)隨機(jī)信號的自相關(guān)函數(shù)和

功率譜密度是一對傅里葉變換對。

6.2平穩(wěn)隨機(jī)序列通過LTI離散時(shí)間系統(tǒng)

輸出序列的均值

輸出序列的自相關(guān)函數(shù)

輸出序列的功率譜

輸入/輸出序列的互相關(guān)函數(shù)及互功率譜

平穩(wěn)隨機(jī)序列通過LTI系統(tǒng)

X[k\

h[k]

*mi

Y比|=Yh[n]X[k-n\

、輸出序列的均值

m[k]二E{y[k]}=工h[n]E{x[k-n]}

=m工h[n]

x

nH(ej0)

m[k]=mH(ej0)

二、輸出序列的自相關(guān)函數(shù)x

R[n]=R[n]*R[n]

yhx

Ry[n]是系統(tǒng)單位脈沖響應(yīng)h[k]的自相關(guān)函數(shù)Rh[n]與輸入隨機(jī)序列X[k]的自相關(guān)函數(shù)Rx[n]的卷積.

系統(tǒng)單位脈沖響應(yīng)h[k]是確定信號，其自相關(guān)函數(shù)定義為

二DTFT{R[n]}-DTFT{R[n]}

hx

uE

^yW-E{y[k]y\kn]}

[Fm

二工工h\l]h\mtK{x\k-/\x\k

三、輸出序列的功率譜

P(0)=DTFT{R[n]}=DTFT{R[n]*R[n]}

yyhx

DTFT{R[n]}DTFT{R[n]}

Px(0)

H(ej0)=

1

1-ae-j0

|h(ej0)|2=

1

1—2acos0+a2

DTFT{h[n]*h-n]}=H(ej0)

P(0)=H(ej0)2P(0)

y

四、輸入/輸出序列的互相關(guān)函數(shù)及互功率譜互相關(guān)

R[n]=E{Y[k]X[k+n]}=h[-n]*R[n]

TOC\o"1-5"\h\z

yxx

R[n]=E{x[k]y[n+k]}=h[n]*R[n]

xyx

互功率譜

P(0)=DTFT{R[n]}=H(ej0)P(0)

xyxyx

P(0)=DTFT{R[n]}=H*(ej0)P(0)

yxyxx

[例]一離散時(shí)間平穩(wěn)白噪聲通過一階IIR數(shù)字濾波器

y[k]=x[k]+ay[k-1]p|<1

求輸出的自相關(guān)函數(shù)、平均功率和功率譜。零均值白噪聲的特征

E{X[k]}=0P(0)=b2

R[n]=J"c2ejn0dQ=Q28[n]

x2n-n

解：

H(z)=-z>|a|

1一az-1

h[k]=aku[k]

1）計(jì)算輸出的自相關(guān)函數(shù)

R[n]=蘭

h

k=一8

h[k]h[k+n]=

k=-n

an

akan+k=

1一a2

a一n

akan+k=

1一a2

anl

1-a2

c2anl

1-a2

R[n]=R[n]*R[n]=c2R[n]yhxh

2）輸出平均功率

c2

1-a2

1

lim—蘭y2[k]=R⑹=nT82N+1y

Nk=-N

3）輸出功率譜

P(⑵=H(ea)2P(0)=

yx

c2

1—2acos0+a2

6.3經(jīng)典功率譜估計(jì)

■譜估計(jì)的質(zhì)量

■相關(guān)法（間接法）

周期圖法（直接法）

周期圖法的改進(jìn)

利用MATLAB實(shí)現(xiàn)功率譜估計(jì)

一、譜估計(jì)的質(zhì)量

1．估計(jì)量的偏差

bia{0}=E{0}—0

2．估計(jì)量的方差

var{0}=E{（0—E{0}）2}

若bia{0}=0,limvar{0}=0,則稱0為0的一致估計(jì)

NT8

3．估計(jì)量的均方差

MSE{0}=E{（0—0）2}=var{0}+bia{0}

二、相關(guān)法（間接法）進(jìn)行功率譜估計(jì)

相關(guān)法的理論基礎(chǔ)

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

相關(guān)法進(jìn)行功率譜估計(jì)

功率譜估計(jì)的質(zhì)量

1.維納—辛欽定理

R[n]P(Q)

xx

計(jì)算方法：

(1)由隨機(jī)序列一個(gè)樣本的N個(gè)觀測值計(jì)算自相關(guān)函數(shù)的估計(jì)

(2)對R[n]進(jìn)行DTFT即得該隨機(jī)序列的功率譜估計(jì)P(Q)

xx

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

X[k]是寬平穩(wěn)各態(tài)遍歷隨機(jī)信號，x[k]是其一個(gè)樣本

R[n]二

x

lim

NT8

迓x[k]x[k+n]

k=-N

已知x[k]的N個(gè)觀測值x[O],x[l],...,x[N-l],則自相關(guān)函數(shù)的估計(jì)為

R[n]=x[k]x[k+n]

xN

k=0

=x[-n]*x[n]-(N一1)<n<N一1

N

[例]已知平穩(wěn)各態(tài)遍歷的實(shí)隨機(jī)序列X[k]的單一樣本的N個(gè)觀測值為x[k]={1,0,-1}，試計(jì)算該隨機(jī)序列的自相關(guān)函數(shù)估計(jì)。解:

R[n]=丄x[-n]*x[n]=丄{-1,0,2,0,-1}

xN3

-1

-1

/■

0.

0.-

0.-

-1

1.--

利用MATLAB計(jì)算相關(guān)函數(shù)的估計(jì)

利用conv函數(shù)計(jì)算

2.利用數(shù)字處理工具箱中提供的函數(shù)xcorrxcorr(x,y);%隨機(jī)序列X和Y的互相關(guān)xcorr(x);%隨機(jī)序列X的自互相關(guān)

利用DFT計(jì)算自相關(guān)函數(shù)的估計(jì)

對x[k]補(bǔ)零形成L點(diǎn)序列

x[k](L>2N-1)

L

X[m]=DFT{x[k]}

LL

R[n]二；IDFTx厶[m]|2}

相關(guān)法進(jìn)行功率譜估計(jì)

x[k]估計(jì)>R[n]DTFT,DFT>p(⑵,p[m]xxx

.in乎nii

&[n]二乙x[k]x[k+|n|]

k=0

L"'

P(⑵=乙R[n]e-g,|L<N-1

xx

n=-L

1N-1-n

R[n]=乙x[k]x[k+|n]|n|<N—1

xN

k=0

x[k]

0N-1

0<n<N-1

-n

x[k+n],n>0

N-1-n

R[n]=f"x[k]x[k+n]xN

k=0

x[k+n],n<0

-nN-1-n

-(N-1)<n<0

R

k=-n

x[k]x[k+n]=-

N

fx[l]x[l-n]

l=0

[例]已知實(shí)平穩(wěn)隨機(jī)序列X[k]單一樣本的N個(gè)觀測值為x[k]={1,0,-1}，

試?yán)孟嚓P(guān)法估計(jì)其功率譜。

解：X[k]的自相關(guān)函數(shù)估計(jì)值為

2,0,-1}

R[n]=—x[-n]*x[n]=]{-1,0,

N3

對RR[n]進(jìn)行傅里葉變換得X[k]的功率譜估計(jì)

x

八12

P(0)二DTFRR[n]｝二一｛-ej2□+2-e-j2q｝二一(1-cos20)xx33

相關(guān)法功率譜估計(jì)的質(zhì)量

功率譜估計(jì)的質(zhì)量與自相關(guān)函數(shù)估計(jì)的

質(zhì)量密切相關(guān)

bia{R[n]}=—呱[n]

xNx

var｛R[n]｝沁+n]R[r—n])

N

r=—g

NT8,偏差、方差趨于零，是一致估計(jì)。

N固定時(shí)，nTN，偏差、方差較大

三、周期圖法(直接法)進(jìn)行功率譜估計(jì)

周期圖法功率譜估計(jì)的計(jì)算

周期圖法功率譜估計(jì)的質(zhì)量

周期圖法功率譜估計(jì)的計(jì)算

已知：xjk]=；；伙]

k二0,1,,N-1

其它

方法基礎(chǔ):

[-n]*x[n]

N

1R[n]=xxN由維納—辛欽定理

1

P(Q)=DTFT{x[n]*x[-n]}

TOC\o"1-5"\h\z

xNNN

=—X(ejQ)X(e-jQ)=—X(ejQ)NNNN'

周期圖法功率譜估計(jì)的步驟

x[k]-dtft>X(ejQ)—功率譜估計(jì)>P(Q)=—|x(ejQ)

NNxN''

x[k]—dft>X[m]—功率譜估計(jì)~tP[m]=—X[m]l2

NNxNN

其中X(ejQ)二DTFT{x[k]}二藝x[k]e-jQk

NNN

k=0

X[m]=DFT{x[k]}=匸x[k]e-jNmk

NNN

k=0

[例]已知實(shí)平穩(wěn)隨機(jī)序列X[k]單一樣本的N個(gè)觀測值為

x[k]={1,0，-1}，試?yán)弥芷趫D法估計(jì)其功率譜。

解：對x[k]進(jìn)行離散時(shí)間傅里葉變換

X(ejo)=￡x[k]e-jok=1—e-j20

N

k=0

功率譜估計(jì)為：

NN

I(0)=—X(ejo)2=—X(ejo)X*(ejo)

NN'■“'

=3(1—e—j2o)(1—ej20)

M=t29

平穩(wěn)高斯白噪聲功率譜估計(jì)結(jié)果(周期圖法)

ffiE-eJoadul-la

01?30

Frequency

■2

遠(yuǎn)pjra二肚alw-li

Frequency

N=512

I2

Finetjuenc/

o

■30

周期圖法功率譜估計(jì)的質(zhì)量

均值E{I(o)}=七-nR[n]e-㈣

NN

n=—(N—1)

E{IN(W)}=Px(W)},漸進(jìn)無偏估計(jì)

、■、八

方差

var{I(o)}=g

N

4{1+

'sinN。)丫}

.Nsi￡2丿

N增加,方差不減小,不是一致估計(jì)

四、周期圖法的改進(jìn)

問題的提出

平滑周期圖(Blackman-Tukey法)

平均周期圖法(Bartlett法)

重疊平均周期圖法(Welch法)

1.問題的提出

周期圖法進(jìn)行功率譜估計(jì)，方差不隨N的增加減小.

如何提高譜估計(jì)質(zhì)量?減小方差方法:

對自相關(guān)函數(shù)估計(jì)值加窗

將N個(gè)觀測值分段,計(jì)算各段的周期圖,再取平均

平滑周期圖(Blackman-Tukey法)

對自相關(guān)函數(shù)估計(jì)值加窗,將誤差較大的估計(jì)值截去

P(Q)二龍1w[n]R[n]e-加。

Mx

n=-(N-1)

窗函數(shù)w[n](M<N)滿足下述條件

0<w[n]<w[0]二1

<w[n]二w[-n]|n<M

w[n]二0|n|>M

B-T法進(jìn)行功率譜估計(jì)的主要步驟

(1)利用觀測數(shù)據(jù)估計(jì)自相關(guān)序列。

(2)對自相關(guān)函數(shù)估計(jì)值加窗。

(3)計(jì)算加窗后自相關(guān)函數(shù)的DTFT。

優(yōu)點(diǎn)：PM(W)波動(dòng)比IN(W)小，可證是一致估計(jì)

缺點(diǎn)：降低了頻率分辨率

平均周期圖法(Welch-Bartlett法)

將隨機(jī)序列X[k]的N個(gè)觀測值分成A段

xi[k]二x[iM+k];i=0,1,…A-1;k二0,1,…M-1

第i段序列的周期圖為

Ii(⑵二丄|X(e；n)|2

MMIM

平均周期圖

Pa(⑵二丄t1Ii(⑵

MAM

i=0

平均周期圖法估計(jì)質(zhì)量

var{PA(0)}=—var{Ii(0)}

MAM

A8f方差為零，是一致估計(jì)

因?yàn)閎ia[I(0)]=-匸旭R[n]e-㈣

NN

n=-(N-1)

所以bia[Ii(0)]>bia[I(0)]

MN

平均周期圖方差減小的代價(jià)之一是偏差增大

4.重疊平均周期圖法(Welch法)

平均周期圖法優(yōu)點(diǎn)：減小方差

缺點(diǎn)：增加估計(jì)的偏差，降低了譜的分辨率原因：分段即加窗，段越多，窗越短，主瓣寬度越大解決方法：將各段數(shù)據(jù)有一定程度的重疊

i=1i=3i=A

NT

.MXM■??LM

i=2

平穩(wěn)高斯白噪聲功率譜估計(jì)結(jié)果(Welch法)

0123

Frequency

)Bd(lartcepSrewoP

0123

Frequency

0123

Frequency

)Bd(lartcepSrewo

20

1

0

-20

M=32

-40

0123

Frequency

五、利用MATLAB進(jìn)行非參數(shù)功率譜估計(jì)

周期圖法

[Pxx,F]=PERIODOGRAM(X,WINDOW,NFFT,Fs)

X：進(jìn)行功率譜估計(jì)的輸入有限長序列；

WINDOW:指定窗函數(shù)，默認(rèn)值為矩形窗(boxcar)；NFFT：DFT的點(diǎn)數(shù)，NFFT>X,默認(rèn)值為256；

Fs：繪制功率譜曲線的抽樣頻率，默認(rèn)值為1；Pxx：功率譜估計(jì)值；

F：Pxx值所對應(yīng)的頻率點(diǎn)

Welch-Bartlett平均周期圖法的MATLAB實(shí)現(xiàn)

[Pxx,F]=PSD(X,NFFT,Fs,WINDOW,NOVERLAP)

X,NFFT,Fs用法同periodogram函數(shù)；

WINDOW：指定窗函數(shù)，默認(rèn)值為hanning窗；

NOVERLAP指定分段重疊的樣本數(shù)。

如果使用boxcar窗，且NOVERLAP=O，則可得到

Bartlett法的平均周期圖。

如果N0VERLAP=L/2，則可得到重疊50%的Welch法平均周期圖。

[Pxx,F]=PWELCH(X,WINDOW,NOVERLAP,NFFT,Fs)[例]一序列含有白噪聲和兩個(gè)頻率間隔很近的余弦信號，設(shè)

x[k]=cos(0.3兀k)+cos(0.32兀k)+w[k]

分別采用周期圖法和Welch法估計(jì)該序列的功率譜。

%PowerSpectralEstimation:Periodogram

N=512;Nfft=1024;Fs=2*pi;

n=0:N-1;

xn=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n));

XF=fft(xn,Nfft);

Pxx=abs(XF)?八2/length(n);

index=0:round(Nfft/2-1);

f=index*Fs/Nfft;

plot(f,10*log(Pxx(index+1))),grid

或直接采用periodogram函數(shù)

%PowerSpectralDensityUsingWelchAlgorithmN=512;Nfft=1024;Fs=2*pi;

n=0:N-1;xn=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n));

L=input('L=')

window=boxcar(L);

noverlap=L/2;

[Pxx2f]=psd(xn,Nfft,Fs,window,noverlap);plot(f,10*log(Pxx2)),grid周期圖法進(jìn)行功率譜估計(jì)的結(jié)果

PSDEstimation:Periodogram

000

24

)zH/Bd(ytisneDlartcepSrewoP

-60

0

0.5

1.5

2.5

Frequency(Hz)

Frequency

Welch法，M=128的譜估計(jì)結(jié)果

6.4現(xiàn)代譜估計(jì)簡介

問題提出

平穩(wěn)隨機(jī)信號的參數(shù)模型

AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系

AR模型參數(shù)與線性預(yù)測濾波器的關(guān)系

Y-W方程的L-D遞推算法

伯格(Burg)遞推算法

利用MATLAB進(jìn)行AR模型功率譜估計(jì)

一、問題提出經(jīng)典法存在問題:

方差性能不好，不是Px(W)的一致估計(jì)

平滑周期圖和平均周期圖改善了周期圖的方差性能，但卻降低了譜分辨率和增大了偏差。

可能使短序列的功率譜估計(jì)出現(xiàn)錯(cuò)誤的結(jié)果出現(xiàn)問題的原因:

將觀測數(shù)據(jù)以外的數(shù)據(jù)一律視為零，與實(shí)際不符。

參數(shù)模型法的基本思想根據(jù)所研究信號的先驗(yàn)知識，對觀測數(shù)據(jù)以外的數(shù)據(jù)作出某種比較合理的假設(shè)。

方法：選擇一個(gè)好的模型，在輸入是沖激函數(shù)或白噪聲的情況下，使其輸出等于所研究的信號，至少也是對該信號的一個(gè)良好近似。利用已知的自相關(guān)函數(shù)或數(shù)據(jù)求模型的參數(shù)。)利用求出的模型參數(shù)或數(shù)據(jù)估計(jì)該信號的功率譜。

二、平穩(wěn)隨機(jī)信號的參數(shù)模型

AR模型

1

AZ)

n=0

MA模型

H(z)=1+為bz-

l

l=1

B(z)

麗

工bz-i

ARMA模型

l

—1=0

1+az-n

n

n=0

若輸入白噪聲的功率譜

P(Q)=G2e

則輸出序列的功率譜為

P(Q)=H(ejQ)|2P(Q)=b2H(ejQ)|2

若能確定模型中各參數(shù)an和b就可以求得功率譜Px(Q)

三、AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系

y[k]+Xay[k—n]=w[k]

n

n=1

Yule-Walker(Y-W)方程

_R⑹

R[1]

A?

R[2]

A?

…R[p]

A?

_1-

b2

y

R[1]

y

R[0]

A?

y

R[1]

A?

7

?…R[p—1]

A?

a

1

a

2

0

y

R[2]y.

y

R[1]

y

R[0]

y

7

…?R[p—2].y.

=

0

R[p]

y

R[p-1]

y

R[p-2]

y

…R[0]

y

a

p

0

若已知Ry[n]，由Y-W方程解出各參數(shù)al,a2,…,ap，則可由AR模型參數(shù)獲得功率譜Py(Q)的估計(jì)值。

四、AR模型參數(shù)與線性預(yù)測濾波器的關(guān)系

前向線性預(yù)測濾波器

y[k]的預(yù)測值y[k]由其過去值y[k-1],y[k-2],...,y[k-p]的線性加權(quán)得到。

前向預(yù)測誤差

ef[k]二y[k]-y[k]二y[k]+》a(n)y[k-n]

pp

n=1

前向預(yù)測誤差濾波器系統(tǒng)函數(shù)

y[k]——?A(z)——?ef[k]

p

后向線性預(yù)測濾波器

由y[k],y[k-1],…,y[k-p+1]p個(gè)數(shù)據(jù)預(yù)測數(shù)據(jù)y[k-p]

后向預(yù)測誤差

eb[k]=y[k-p]-y[k-p]

p

=y[k—p]+另a(n)y[k+n—p]

p

n=1

前向預(yù)測誤差濾波器系統(tǒng)函數(shù)

y[k-p]Ab(z)——?ep[k]

Ab(z)=z-p[1+》a(n)z“]=z-pA(z-1)

p

n=1

五、Y-W方程的L-D遞推算法

(1)計(jì)算自相關(guān)函數(shù)的估計(jì)值

(2)由自相關(guān)函數(shù)的估計(jì)值，遞推a,a，…,a,c2

12pp

a1(1)=

R[1]

—1

R⑹

y

c2

=R[O](l-|a|)2

y11

Ry[p]+#

a(p)=—n=1

p

a(n)R[p-n]

p-1y

c2

p-1

a(n)=a(n)+a(p)a(p-n)pp-1pp-1

(n=1,2,…,p一1)

c2

p

=[1—a

p

(3)求出功率譜估計(jì)

P(Q)=

AR2

1ae-jQn

n

n=1

［例］利用L-D算法對x［k］=cos(0.3p)+cos(0.32p)+耳［k］進(jìn)行譜估計(jì)

cniIiiil

00511.522S3

Fneqiuency(Hi)

-匸

頻率分辨率比N=128點(diǎn)高

L-D1p^B01N=512

■日口IIiiii

00.511.522.53

equency(Hz)

六、伯格(Burg)遞推算法

L-D算法缺點(diǎn)：在計(jì)算相關(guān)函數(shù)估計(jì)時(shí)，對N個(gè)觀測數(shù)據(jù)以外的數(shù)據(jù)作零的假設(shè)，故譜估計(jì)誤差較大。

伯格(Burg)遞推算法基本思想：

直接從觀測的數(shù)據(jù)利用線性預(yù)測器的前向和后向預(yù)測的總均方誤差之和為最小的準(zhǔn)則來估計(jì)反射系數(shù)，進(jìn)而通過

L-D算法的遞推公式求出AR模型優(yōu)化的參數(shù)。伯格(Burg)遞推算法步驟

(1)確定初始條件

ef[k]=eb[k]=y[k]

00

(2)從p=1開始迭代計(jì)算:

-2藝

b2=*y[k]

0N

k=0

計(jì)算AR模型參數(shù)

ef[k]eb[k-1]

p-1p-1

K=k=p

P藝{ef[k]2+eb[k-1]2}