电信往年卷-现代数字信号处理spectrum

五、功率谱估 确定性信号：Fourier变换－频谱特周期信随机信号：不满足绝对可积条 TXT()T

x(t)ejwtS()limE

()2 T NN (ej)x(n)eNnXN2(XN2(ejS(ej)lim N （Wiener-Khinchin）定理Sx()

m

R(m)e 谱分析：利用有限从观察样本估计信号的功率分布情寻找数据的“隐周期性语音处 生物工无源声纳信号处信号处理：回波信号的功率谱密度分 应用例1，脑电波信号1（正常的脑电波脑电波信号2（癫 信号1

9信号2 不同方法得到的功率谱估 3,序 4,

ˆ(z)101

0dB

6

0

10dB 19世纪末（1899），舒斯特（Schuster）－周期1930年，维纳：维纳－ 20世纪50年代，BT1967年，Burg－最大熵谱1968年，Parzen－AR谱估计1971年，VanDenBos证明最大熵谱法与AR谱估1969年，Capon 经典谱估计法、现代谱估计算非参量方法、参量方

AR,ARMA模 Pisarenko谐波分 偏方谱分辨 估计方信号模外推方 概S()1x(0)x(1)ejx(N1)ej(N1)NxS(ejx

XNXN

N (ejw)

(ej)

N x(n)ejnNN lim NNN

limN

n ˆx() N

(ej) NNN

N x(n)en

()

N

ˆ(m)

jm(N1)

ˆ()

X(ej)

N N Nx(i)e N

N

kNN

(k)x

(ik))ej(ik i0kNN

N N

(k (km)e

(mik m(N1)k(N

ˆ

em(Nˆ(m)

Nx*(n)x (n其中

N

ˆx(m)R(m);ˆx()Sx() 功率谱估Sˆx()

N

ˆ

(m)ejmm(N1)E[ˆx()]

N

E[ˆ

(m)]ejmm(N1)E[ˆx()]F[E(ˆx(m))] 由E{ˆ(m)}N|m|R(m)w(m)R

N|m

|m|

(m10.80.60.40.20(N1)

N

E[ˆx()]F[w(m)Rxˆ()]1S()W()1S()W(

2 NW()

1sin(2) N

N )

m(N MagnitudeMagnitude--

Normalized

N

E[ˆx(m)]R(m);E[ˆx()]ˆx(实际中数据长度N有限，偏差总是存1xbia(ˆ())1x

Sx()W()dSxx(n)sin(2f*n)u(n) Cv(ˆˆxE(ˆxE(ˆxE(ˆx11 S()D0(1)D0(2)d11

S()D0(1)D0(2)d

|n|N d0(n)1 11

D02 S()D0()D0(

E(SN

ˆx{E[ˆx

ˆx

x|S( x 111 Var(Sx

))

S()D0()D0(

E(Sx假设窗函数旁瓣为0，主瓣宽度为

B(B D0()D0() 2Var(ˆxE(ˆx2 211ˆ Cov(S

1),S(2))

11

S()D0

)D0

2

|12|有：D01)D02

),

)) D0((a

(b

B2

2

D0( D0( D0(B22(( D0( 1202(c (d

1 1

D0(D0(1D0( 2B B2 随着N的增大，周期图的谱估计起伏增 真实谱中两个靠得很近的谱峰能被分辨的能ˆ()]S()W()1 1958年由Blackman与Tukey提1N11R(m) x(n)x(n|m|) n S( R(m)e jm 估计自相关函 1N

|M|<=N-Sˆ(

ˆ(m)v(m)ejMm M ˆBTˆPER BTBT法的方差小于周期图法的方BT法的分辨率比周期图 方差性能不好，不是一致估分辨率 假定x(n)各态历 R(m) 1 x*(n)x(nN 2N1N 2NS()2NN

x(n)ejnnN

2 Sˆ() N

N x(n)ejn 选取主瓣窄、旁瓣幅值小且衰减快的窗函 5.2.3方差性能不好，起伏大：缺少统计平分辨率平平 5.2.3.2Bartlett平均，改善方差性牺牲分辨基本原理：将长度为N的数据分为LM，对每段数据用周期法进行谱估计MxM1xM

i(n)x[n(i1)MS1

Si()M

M x M x MMLSˆ() 1Si( i 5.2.3.2Bartlett 设新的

X L L 方差为2 5.2.3.2Bartlett

Sˆ()L

Si([ˆ()]1S()W Var(Sˆ())1S()L 5.2.3.3Welch分每段数据用平滑窗平滑处M

E[ˆ 2

1 |(n)| nVar(Sˆ())9LS noverlap：数据分段 的样本 缺省值：nfft： 50-- - - - 50-- - - - 50-- - - - 例 3 21

0

0

用法：pwelch(xnfft：FFT长noverlap：数据分段 的样本 5.3对随机信号建根据观测值估计模型参求功率x(n)A1exp(j2nf1)A2exp(j2nf2)v(n)S(f)A(ff)A(ff) 5.3H(zuH(z

X(n 5.3离散随机信号的线性模型- kx(nk)ku(nk k

00qz qz

i p i 5.3 1

i

jSx()H(

jw)

*(

jw)

i p

1选择合适的模型描述观测数选择合适的模型阶

i

i

j选择合适的方法估计模型的参 5.3ARMA 估计模型H(z)的参由H(z)的参数估计x(n)功率选择合适的模选择合适的模型阶选择合适的方法估计模型的参 p阶AR模px(n) x(ni)u(npi iRx( i),

(m ip p i

iR

( i)2

矩阵形Rx

x) Rx(p)1

R

Rx(p1)1 0x x

R R

R(0) 0

p 2R

估计观察序列x(n)自相关函选择适当的AR模型阶数解线性方程组，求出参计算功率谱，公式高斯消元法，Cholesky分解 AR参数谱估计与最佳线性预测的关计算量大－>快速算Toeplitz矩阵的性利用Toeplitz矩阵性 Rx(p)R R R(pR R( R(p R(0) 已知kk阶的参数求出k+数?阶数递ˆ0 ˆ 1利用x(n-p),x(n-p+1x(n-1)的线性组合预测ppxˆf(n k

kx( k

e(n)x(n)ˆf(nE[e2(n)]

x(n)

ppk

kx(nk

2 E[x(nm)e(n)]

ˆ m1,2,,yxx0eRx(m)kRx(mkk最小预测均方误差 R(0)p ˆ

kRx

投影 线性预测Wiener－Hopf方程：ak和Rx Rx

Rx(p) 1

R

R

R(p1)

0

1 R(

R(p

R

0 p 与AR模型的Yule－Walker方程相 kk k1,2,, 2ppe(n)x(n)ˆf(n)x(n)x(nk)kk (u(n(u(n)1A(z)11az1az (b)A(x(n)ˆ(b)A(x(n)ˆ(ˆ(e(n)1A(11az1az1paz1az 前向预测：利用p个n时刻之前的数据预测 Forward

Backward xˆf(n)af(k)x(nk

xˆ(np)a(k)x(npk k

k (n)x(n)

e(n)x(np)ˆ(n

Eef(n)2

bEeb(n)2

利用正交性原理，可以得到后向预测的Weiner－Hopf方程0 0bbbmi

Rx(0)ab(k)k

x(kpRx(m)ab(k)Rx(mk),k

m1,2, b

f

af(k)ab(k 低阶AR参数＝高阶AR模型参AR模型与线性预测模型的参数一递推线性预测系数＝>AR模型参am1,am2,…,

mi

mi 00

5.3.3Levinson-Durbin0ef(n) x(n0 Eef(n)

E

2(n)R(0

R

R(0) R(1)1

10

R(0)a

0

R(0)1a2a2 E{ef(n)eb(n mm1R(m1)amiR(m1 证明：mm+1

1 m

R(m

0

R(m

0

mm

M+1M+1

R(m

m1 R(m

0

m1,1

R(m) R(m

a0 m1,m a0 R(m)

0 M+1M+1 R(m) R(m 1

R(m

R(m)

m1 R(m) R(m

a

mm

R(m)

0

mm1R(m1)amiR(m1m R(m

R(m)

1

m1 R(m) R(m

m1 R(m

a

mm

R(m 0 m1R(mRR(mR(m)R(m)R(m R(m) mm

m1

1

M+1M+1

R(m

m1 R(m

0

R(m) R(m

0

m1,m R(m)

0 Km1

R(m1)amkR(m1 m

R(m

0 R(m1 m1 Km1

mm

m1

0

1

m10 0

0 0 Km1

R(m

m1 R(m

0

m1,1

R(m

0

m1,m

R(m)

0

R(m1)1

0am1,m1 amkkm1am,m1k k1,, (1k

R(m1)amkR(m1kkm1

m

k m

R(m

m1 R(m

0

m1,1

R(m

0

m1,m

R(m)

0 5.4.1AR0ef(n)x(n0K0 Eef(n)2

R(0 5.4.1AR11阶递ef(n)x(n) x(n1) a11K1R(1)(1K2)(1K2 22阶递ef(n)x(n) x(n1) x(n a21(n)a11K2a11(1K2a22K2R(2)a11R(1)(1K2 p0 R(0) Ex2(n)0K 1,2,...,am1, am Km1am,m1am1,m1 Km (1K )m m

1k Km

m R(m1)amkR(m1k) SˆA

(

k 2p e ek

jp 递推过程的每一步运算量为所有递推过程的计算量为11-|km|当|km|=0时，递推算法停 m2m-12…2 随着模型阶数的增从最小均方误差准则的意义，算法收 ef(n)e (n)k m (n (n

(n1)ke

mE (nmE

(n)eb

(nKm

EE

,m1,2,...,ef(n)eb(n)x(n

ef

ef ef

ef0

ˆ(n)

kp eb

eb Levinson-Durbin递归算 Burg自相关函数估 确－>不估分辨率较自相关函数法线 一谱的置分裂为条谱线），模可产虚值 E[e2(n)]

x(n)

ppk

kx(nk

2 ˆ(m NnNn

emf(nem

Nn

2bem(nb 1 b ef(n)e

(n)k

(n

eb(n)eb (n1)ke

(n)

N2e

(n

ˆm ,

m1,2,, Ne

N

(n

|ef(m1)|2|

mˆm

1,2,,m, m,

1kmm,

m (1| m

m110ˆ和1阶AR1计算m阶时的AR模型系重复上述过程，直到m=p，求出所有阶次的AR参 不需估计自相关函N较小时，性能不亚于Levinson算满足|ki|1，AR模型总是稳定Bur位的敏感问题 fb

1( 2

p p

N

x(n)

f(k)x(nk pn

k

N

N1n

x(n)apkp

(k)x(nk f

(

f ( aˆ(i

ˆi)

N1 ˆNpx(n)a(k)x(nk) (ni)np

k N1p (n)a(k) (nk)x(ni)n

k

N

N1

x(nk)x(ni)a(k)x( k) (n i) k

n N1x(n)x(ni)

N1n

x(n)x(ni) N

N1cx(i,k)

2(N

n

x(ni)x(n

k)

n0

x(nk)x(ni) cx

cx ˆ()

c

c

ˆ(2)

c(2,0)

c( c(

c(p,p)aˆ(

c(x不是矩

x(n)ˆ(k)x(nk)x2(Np)np

k N1p

(n)a(k (nk)

k

cx(0,0)ˆ(k)cx(0,kk

ARAR谱比经典谱平AR谱分辨率比经典谱要 AR号初相位的影响，并且可能出现“谱线” ARFPE(k)Nk

N AIC(k)Nlnˆx

CAT(k)

N 例f1f2x(n)0.5exp(j2nf1)exp(j2nf2) 50

0

50000

00

qx(n)u(n)b(k)u(nkkrx(m)

(ej)21b(k)ejkq2 b(k)b(km)2

q2b(k)b(k

m0,1,,

k kS(ej) (ej)MA

非线性方程非线性方程很难直接求qB(ej)2q

xr(m)exmqSˆBT(ej) rˆx(m)ejm 谱分解高阶的AR模型近似MA模最大似然估计法（最小二乘法 5.7.2MA H(z)

1A(z

1kq

1a(k)zH(z)Hq(z)1b(k)zkA(z)B(z)

qa(m)b(k)a(mk)k 5.7.2MA取有限的pqˆb(k)ˆkb(k)可以看作ˆp)b(k)可以看作实际上，e(m)不为ˆM

e(m |m S(ej)21b(k)ekAIC(q)Nln(ˆMA)2q 5.7.3ARMA谱估 2a(k)rx(mk) G(k)b(m k0,1,,kp

ka(k)rx(m

m k rx(q rx(q1)

rx(q1)rx(q

rx(qp1)a(1) (qp2)a(2)

(qp1) (qp2)

(q a(p