第五章功率譜估計(jì)(Power

Spectrum

Estimation)1參考章節(jié)1：第11章奧本海姆：第11章離散時(shí)間信號(hào)處理：第11章232§5.1

估值與一致估計(jì)條件1.算術(shù)平均及與概率平均的關(guān)系算術(shù)平均：n

x

lim

nN

1Nn0m

E[x

]x

nN

1

平穩(wěn)且x

各態(tài)歷經(jīng)(1a)自相關(guān)算術(shù)平均：

13nx

x

*(n

m)

x

x

limN

1n0N

Nn

nm

xx

(m)

E[xn

x

*(n

m)]平穩(wěn)且各態(tài)歷經(jīng)(2)§5.1

估值與一致估計(jì)條件由3-(26b)，mx＝0時(shí)：xxxxm(m)e

jmP

()

以<xnxn+m>取代xx(m)，則功率譜算術(shù)平均：24(2)1lim

lim

1N

1N

1

nxx

xxN

Nx

x

*(n

m)

em

N

N

n0x(n)e

jnn0jm

P

()

P

()

見(jiàn)附頁(yè)(3)（這也稱為周期圖，其中X

N

()是有限長(zhǎng)xN

(n)的付氏變換）5§5.1

估值與一致估計(jì)條件值平均的估值：2.從算術(shù)平均到有限長(zhǎng)段的估N

1

x(n)n0

1Nm?

x

(1b)自相關(guān)估值：N

1n0n0xxxN

(n)xN

(n

m)N

1x(n)x(n

m)

1NN(m)

1?(4)由式(3)，功率譜估值：211NxxX

()

2N

NP?

(m)

N

1n0

x(n)e

jn

(5)§5.1

估值與一致估計(jì)條件3.一致估計(jì)條件用于衡量某估值是否接近真值（概率平均）的準(zhǔn)則。（1）偏畸（Bias）定義：

B[?]

E[?]式中，

為真值，

?

為估值。(6)(7)當(dāng)B=0時(shí)，

E[?]，稱為無(wú)偏（畸）估計(jì)。（2）估值方差（variance）定義：?6??var[22?

E[

]]

}

E{[]

7§5.1

估值與一致估計(jì)條件(8a)由于：及得：即，估值的均方誤差為偏畸平方與估值方差之和。2∴E[e

]取決于B與var，只需考慮B、（3）估值的均方誤差定義：E[e2

]

E[(?

)2

](7)2?

E[?

2

]

E

2

[?](8b)(6)B2

2

2E[?]

E

2[?]

B2

E[?2

]

2E[?]

2

E[e2

]2?(9)2?8§5.1

估值與一致估計(jì)條件?xx

xx（4）一致估計(jì)條件用于：N->∞時(shí)，估值是否接近真值，若是，為正確估計(jì)。比較：(1a)與(1b)、(2)與(4)，當(dāng)N->∞，m?

x

mx

,估值即真值。估值即真值，稱一致估計(jì)條件。(10)var[?]

0即：當(dāng)N->∞時(shí)滿足：

B[?]

0§5.1

估值與一致估計(jì)條件(11)4.Davenport定理（隨當(dāng)隨機(jī)過(guò)程{xn}為白色、實(shí)量xn的四階矩）分布時(shí)，有：E[x1

x2

x3

x4

]

E[x1

x2

]E[x3

x4

]

E[x1

x3

]E[x2

x4

]

E[x1

x4

]E[x2

x3

]9§5.2

自相關(guān)(自協(xié)方差)估值及其表示式1.?xx(m)的一般表示式由(4)式求?xx(m)，僅給出N個(gè)x(n)相繼值，故也可僅給出(N-m)個(gè)x(n+m)，因此(4)式可修改為：(12)'xx10CN

|m|1(m)

x(n)x(n

m), |

m

|

N1N

|

m

|n0用|m|可使m也適用于負(fù)值。上式是否滿足一致估計(jì)條件？§5.2

自相關(guān)(自協(xié)方差)估值及其表示式（1）是否N->∞時(shí)B->0（或B=0）(13)xxN

|

m

|1

E[x(n)x(n

m)]1N

|

m

|E[C

(m)]

N

|m|1N

|m|1n0'xx

(m)

xx

(m),

|

m

|

Nn0(6)'(m)

E[C

(m)]

0xx

xx

B

是無(wú)偏估計(jì)。（2）是否N->∞時(shí)，Var[C’xx(m)]->0xx11xxxxxxxxxx

E[C

(m)]

E

[C

(m)]2'

2(13)2

'2'(8b)(m)由(7)式，var[C

'C

'

E[C

(m)]

(m),|

m

|

N(m)]

2(14)§5.1

估值與一致估計(jì)條件由(1)及(2)的結(jié)果，以(12)式表示的C’xx(m)是滿足一致估計(jì)條件的。但顯然只有|m|<<N時(shí)才滿足，當(dāng)(對(duì)于一個(gè)固定的延遲)|m|越接近N時(shí)，Var[C’xx(m)](的偏差)就越大。其中利用(11)式求得E[C’xx

(m)]代入(14)式，再變量置換，2可得下式：212'

limvar[C'

(m)]

0N

1(N

|

m

|)2

N(r

m) (r

m)](N

|

m

|)2NxxN

N

|m|1r

(

N

|m|1)xxxx

xxxxvar[C

(m)]

N

|

m

|時(shí)其中的[

(r)

(15)§5.2

自相關(guān)(自協(xié)方差)估值及其表示式2.自相關(guān)估值的改進(jìn)表示式若將(12)式修改為：(17)NxxC

(m),

|

m

|

N1Cxx

(m)

N

x(n)x(n

m)N

|m|1n0'(12)N

|

m

|(16b)'xx13(6)(17)

偏畸B[Cxx

(m)]

xx

(m)

E[Cxx

(m)]滿足N

時(shí)，B

0(m)

N

|

m

|

(m)

|

m

|

(m)NNxx

(m)N

|

m

|(13)N

|

m

|xx

E[C

(m)]

E[C (m)]

Nxx

xxNxx(16a)(18)§5.2

自相關(guān)(自協(xié)方差)估值及其表示式(19)'142222

22

((14

)

()

Var

[C

(m

)]

Var

[C

(m

)]NN

|

m

|)

{E[CNN

|

m

|NN

|

m

|NN

|

m

|(m

)]Var

[C

(m

)]

E[C

2

(m

)]

E

2

[C2

'xx

xxxx'xxxx2

')

C

xx

(m

)]

((16

b

),(17

)

E[(xx(8

b

)xx

xx(m

)]

(m

)})

(m

)當(dāng)然滿足N

時(shí)，Var

[C

xx

(m

)]

0∴Cxx(m)滿足一致估計(jì)條件，且其估值方差比C’xx(m)的小。∴選擇Cxx(m)作為自相關(guān)估計(jì)。又，§5.3

直接用周期圖作功率譜估計(jì)1.周期圖的一般定義1NNNNNNNxxxxxxNxN()

|2

I()

15

1

XN(e

j

)

X

*

(e

j

)

1

|

XN

(16

a

)(m)e

jm

,

|

m

|

Nl

n

(n)e

jnl

n

mm

(n

m)e

j

(

n

m

)n

1

(n)e

jnm

n

m

N

|m|1n

0m

x

(l

)e

jl

x

N

x

1

N

1

N

x

(n)

x

(n

m)e

jmNN

x(n)

x(n

m)e

jmP?()

P?

(

)

C(20)XN(ω)為有限長(zhǎng)序列的付氏變換，具有周期性；I

(ω)具有功率譜N的量綱，作為功率譜估值，稱為周期圖。16§5.3

直接用周期圖作功率譜估計(jì)2.用周期圖直接作功率譜估值是否滿足一致估計(jì)條件xN

(n)x

N

(n

m)1Nx(n)x(n

m)

N

|m|1n0

n(16a)Cxx

(m)

1N(20)(1)周期圖的偏畸BNNNNxxNxx

1

[

R

(

n

)

RNxx(

m

)]

(

m

)

w

(

m

)

(

m

)BN(

n

)

RN(

n

m

)

E

[

x

(

n

)

x

(

n

m

)

]NN

1

Rn

1

E

[

R

(

n

)

x

(

n

)

R

(

n

m

)

x

(

n

m

)

]n

E

[

C

(

m

)]

1

E

[

x

N

(

n

)

x

N

(

n

m

)

]n

(21)§5.3

直接用周期圖作功率譜估計(jì)xxBN

B

xxW

()

P

(

)d2E[I

()]

W

()

P

()

1

(23)其中，三角窗序列如圖，可表示為WB(m)-N

0

N

m1

0,

elseN

N

|

m

|

,

|

m

|

Nw

(m)

B(22a)sin

/

21

sin(N

/

2)

2WB

()

N

(22b)可見(jiàn)，除非WB(θ)=2πδ(θ)，否則E[IN(ω)]≠Pxx(ω),即B[IN

()]

Pxx

()

E[IN

()]0,

∴IN(ω)為有偏估計(jì)。17§5.3

直接用周期圖作功率譜估計(jì)但當(dāng)N很大時(shí)，由(22a)可見(jiàn)，N→∞時(shí)wB(m)→1∴當(dāng)N→∞時(shí)，WB(m)→2πδ(ω)則：(24a)N

lim

E[IN

()]

Pxx

()于是，N

，B[IN

()]

0(2)周期圖的方差var[IN(ω)](25)22NN(8b)Nvar[I

()]

E

[I)]

E[I

(其中，可求得：N18()]4xE

2[I()]

(26)(24b)這表明直接用周期圖作功率譜估計(jì)不滿足一致估計(jì)條件，需修19改§5.3

直接用周期圖作功率譜估計(jì)并利用(11)式展開(kāi)、化簡(jiǎn)可得：

21

1

2

(27)21

1

2

41

N

sin(

)

/

2

N

sin(

)

/

2

sin

N

(

)

/

2

2

sin

N

(

)

/

2

2

E[I

(

)I

(

)]

x

N

1

N

2

1421

2

1

N

sin

sin

N

2)]

E[I

(

(

)]

)I

E[I

(xN

1

N

2

N(28)將(28)式、(26)式代入到(25)式得：

44

N

sin

sin

N

2

N

sin

sin

N

2

x

1

4xvar[I

N

()]

x

2

(29)

04xN)]

N

時(shí)，var[

I

(§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)(30a)5.4.1

特（Bar

ett）平均法1.

不相關(guān)隨

量間的算術(shù)平均方差設(shè)x1,x2,…xk為互不相關(guān)的隨

量，每個(gè)期望均為：kx

ii1i

x

x

1

2

kkE[x

]

m

,

m?

(x

x

x

)

/

k

1xxkx

1

2

kE

[

m?

]

E[x

x

x

]

/

k

1

km

m(30b)kx

x

xvar[xi

]??2?

2(8b)var[m

]

E[m

]

E

[m

]

(31)k個(gè)不相關(guān)隨

量的平均方差比單個(gè)隨

量方差小k倍。20§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)為：2.特法－平均周期圖法將序列x(n)(0≤n≤N-1)分成前后不相接的k段，求各段的周期圖平均。設(shè)每段有M個(gè)樣本，則N=kM，第i段樣本序列可寫(xiě)xi

(n)

x[n

(i

1)M

]

x(n

iM

M

)每段長(zhǎng)度0

n

M

1,

1

i

kMI

iM,

1

i

k12M

1

xi

(n)e

jnn0()

(32)xx

(m)隨m的增大而減少（即相關(guān)性減少）。第i段的周期圖，據(jù)(20)為：21比較(34a)與(29)式可見(jiàn)：平均周期圖法的方差僅是周期圖法方差的1/k?！喽螖?shù)k→∞時(shí)，var[Bxx(ω)]→0,為一致估計(jì)。22§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)若m>M（段數(shù)分得越多，每段長(zhǎng)度M很小時(shí)），可認(rèn)為各段的周期圖IMi(ω)相互度獨(dú)立，不相關(guān)，k段的周期圖的平均（譜估計(jì)）為kiMxxi11kI

()B

()

2(29)

1)

1

M

sin

sin

M

2

or

:

P

(3-(38)k

M

sin

sin

M

2

k1kvar[B

()]

xxxvar[I

i

()]

4

1Mxx白色時(shí)1(33)上式即(30a)式的形式。故利用(31)式可得：(34a)(34b)§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)5.4.2用窗口處理法平滑周期圖1.周期圖的協(xié)方差（covariance）cov[定IN

(義1)：,IN

(2

)]

E{(IN

(1)

E[IN

(1)])(IN

(2

)

E[IN

(2

)])}

E{IN

(1)IN

(2

)

IN

(2

)E[IN

(1)]

IN

(1)E[IN

(2

)]

E[IN

(1)]E[IN

(2

)]}

E[IN

(1)IN

(2

)]

E[IN

(1)]E[IN

(2

)]注意，當(dāng)ω1=ω2=ω時(shí)，(35a)(35b)(25)232241212N

1

N

2NN

NN

1

N

2cov[I

()])]

E[(I

()]

E

[(I

(

),

I

(

)]

var[I

((27)

sin

N(

)

/

22sin

N(

)

/

22E[I

(

)I

(

)]

11212

x

N

sin(

)

/

2

N

sin(

)

/

2

(35c)§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)2xN

1

N

2將以上兩式代入到(35b)，得：由(26)式的推導(dǎo)過(guò)程知：

E[I

(

)]

E[I

(

)]

(36)124

12

1

2

1

)

/

2sin

N(

)

/

2

2

sin

N(

2

N

sin(

)

/

2N

sin(

2)

/

2

2

),

I

(

)]

cov[I

(4xN

1

N

2§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)M

1m(

M

1)C

xx

(m)w(m)e

jmSxx

()

2.用窗口法平滑周期圖(38)其中，w(m)是長(zhǎng)度為2M-1的窗序列。Nxx2E[I

(

)]W(e

j

(

)

)dE[S

()]

1(39)xx25N定義：Sd2j

(

)I

()W

e()

1

(37)即用窗函數(shù)與周期圖卷積后作為譜估計(jì)，以降低方差§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)由(22b)：1

sin

N

/

2

23.用窗口法平滑后的譜方差26

lim

2

(m),

1

sin

N

/

2

2N

N

sin

/

2

WB

()

N

sin

/

2

由(24a)：N

時(shí)，WB

()

2

()(40)27§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)

4xj

(

)j

(

)

ddN

sin(

)

/

2

N

sin(

)

/

2

1

sin(

)N

/

2

2

1

sin(

)N

/

2

2

]W

[e]W[e4

2

N

1cov[Sxx

(1

),

Sxx

(2

)]

21lim

1

sin(

)N

/

2

2

(40)N

N

sin(

)

/

2

2

(

)

得：先求出Sxx(ω)在ω

1、ω2間的協(xié)方差：(41)式中，12

124114xxj

(d2

N

[W

(e2

NN

limN

j

(

)j

(

)j

(

)

))]dlim

cov[Sxx(1

),

Sxx

(2

)]

W[e

j

(1

)

]W[e

j

(2

)

][

(

)

(

)]d

limN

)W

(e

)

W

(e

)W

(e§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)124j

(

)j

(

)2j

(

)w2

(m)M

1m(

M

1)1)W

(e

)

W(e

)]d12

N

4

x

2

NN

大時(shí)，W(e

j

)主瓣窄前一項(xiàng)互偏離大，可略W

2

()為偶函數(shù)，窄帶時(shí)xxvar[S

()]

xx4x28N

W(e

)]d4

2j

(

)

[W

(e2

N

cov[Sxx

(1

),

Sxx

(2

)]

var[Sxx

()](35c)∴N很大時(shí)：W

2

(e

j

)]d移位積分不影響數(shù)值利用帕氏定理且M足夠大時(shí)(42)§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)2914xxNw2

(m)var[S

()]var[Isin

Nm(

M

1)R

(42)()]

(29)

sin

N

21

N

sin

0, 0

R

N

1w2

(m)

2M

1

(

2M)N

NM

1m(

M

1)M

1N

sin

N

較大

1,

0,1其中，

與周期圖法的方差比：∴當(dāng)M比N小很多時(shí)，R<1∴用窗口法可降低譜估計(jì)的方差即調(diào)整窗的總長(zhǎng)度2M-1及作譜估計(jì)的記錄長(zhǎng)度N,總可使R<1§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)......03M/2-

2M-11M/2 M-1i=2Mi=L-1Mi=1Mi=3Mi=L5.4.3

（Welch）修定義：

結(jié)合

特法與窗口法進(jìn)行了修正，并證明了對(duì)長(zhǎng)度N的實(shí)平穩(wěn)隨機(jī)序列x(n)以2:1覆蓋分為L(zhǎng)段時(shí)，不僅滿足一致估計(jì)條件，且譜估計(jì)的方差最小。2:1覆蓋分段：MMn3031§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)分段后的每個(gè)段序列乘上窗序列w(n)，求付氏變換：M

1X

i

(e

j

)

xi

(n)w(n)e

jn

,

i

1,2,...,

Ln0其中，w(n)的長(zhǎng)度為M，可為矩形、三角、升余弦窗。如三角窗時(shí)：,2

2

n

M

1

M-1

22n M

1, 0

n

M

1w(n)

M

12n或：2M

1w(n)

1

2

, 0

n

M

1n

M

132§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)升余弦窗（

M

1

）時(shí)：w(n)

0.54

0.46

cos

2n

,0

n

M

1每段的修正周期圖：MUJ

iM12,

i

1,2,...LM

1ix

(n)w(n)e

jnn0()

(43)與(32)式相比，多乘一個(gè)w(n)，且分母多了U，U代表窗函數(shù)能量：21M

1n0MU

w

(n)(44)33§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)整個(gè)長(zhǎng)度為N的實(shí)平穩(wěn)隨機(jī)序列的功率譜估值是上述L個(gè)修正周期圖的平均：LxxWxxMULiMi1(43)1

1LB

(L M

1

2

xi

(n)w(n)e

jni1

n0J()

)

J

()

(45)這稱welch分段修正周期圖平均法（可以證明此種以2：1覆蓋分段法時(shí)的方差最小）。Welch已證明：xx

xxWxxvar[B

()]

var[S

()]

var[B

()]21N

1n0其中不分段時(shí)簡(jiǎn)化為：WxxN

n0N

1U

1

w2

(n)NUB

()

x(n)w(n)e

jn(46)34§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)在計(jì)算機(jī)上實(shí)現(xiàn)時(shí)，對(duì)(45)式進(jìn)行頻域取樣：,1MULB

(NL M

1

xi

(n)w(n)ei1

n0xxNWxx2

j

2

kn

2

k)

|

BW

(k)

k

0,1,...,

M

1可采用FFT降低計(jì)算量。(47)35§5.4

對(duì)周期圖作平滑后的功率譜估計(jì)功率譜估計(jì)小結(jié)：1.直接周期圖估計(jì)IN(ω)，∵N→∞時(shí)var[IN(ω)]≠0,件，一般不用（對(duì)譜估計(jì)精度要求高時(shí)）。不滿足一致估計(jì)條分段周期圖平滑法Bxx(ω)，滿足一致估計(jì)條件，但必須每段長(zhǎng)度M<<N時(shí)才適應(yīng)。窗口平滑法Sxx(ω)，也滿足一致估計(jì)條件（但仍需窗口長(zhǎng)度

M<<N），優(yōu)點(diǎn)是var[Sxx(ω)]≤var[IN(ω)]且可按要求的方差大小，調(diào)節(jié)窗口長(zhǎng)度M。Welch修改法BWxx(ω)，滿足一致估計(jì)條件，2:1覆蓋分段時(shí)可得到比以上方法更小的方差，缺點(diǎn)是增加計(jì)算復(fù)雜度，但便于用FFT快速運(yùn)算。參數(shù)化譜估計(jì)（略），具有不降低譜分辨率的特點(diǎn)，用于

、目標(biāo)分類、醫(yī)學(xué)圖像處理，而經(jīng)典譜估計(jì)用于一般的工程上。§5.5

功率譜的測(cè)量方法1.偏置（奇型）離散變換（ODFT）N

1ODFT[x(n)]

D(k)

x(n)en0

j

2

(k

1

)nN

2(48)x(n)為實(shí)序列，0≤k≤N-1。N

2

D*

(k)D(N

1

k)

x(n)e

x(n)en0N

1

x(n)en0N

1

x(n)en0N

1N

1j

2

(k

1

)n

(48)2

1j

(k

)nN2

j

2nen0

j

2

(

N

k

1

)nN

2

j

2

(

N

1k

1

)nN

2(49)36§5.5

功率譜的測(cè)量方法2.由ODFT取代DFT運(yùn)算ODFT的偶數(shù)譜線：1)2,

k

0,1,2,...,(D(2k

)

x(n)e(48)

N

1

j

2

(2k

1

)nN

2n0將x(n)分成前后兩半：N37x(n

x(n)

2

n

N

1）（1相當(dāng)于x(n),),

0

n

0

n

N

1222x(n),NNN38§5.5

功率譜的測(cè)量方法]4N,

k

0,1,...,(

1)2)]e

eN2N[x(n)

jx(n

)]2

eN[x(n)

x(n

)e2)eN2

x(n

D(2k

)

x(n)eN

/

21

u(n)en0N

/

21

[x(n)