《数字信号处理》课程教学资源（教案讲义）ch6 功率谱估计

第6章功率谱估计-离散随机序列的特征描述■平稳随机序列通过LTI系统-经典功率谱估计国现代功率谱估计6.1离散随机序列的特征描述国随机过程的分布函数厦随机信号的数字特征平稳各态遍历随机信号的时域描述晨平稳各态遍历随机信号的频域描述（功率谱密度）随机过程的分布函数(X[K],keZ)表示一个随机过程一维分布函数F(x,k)= P(X[kl≤x)二维分布函数F(x,x,;,k,k,)=P(Xkl<≤x,X[k,l≤x)N维分布函数F(xj,x2,"",xn;kj,k2,"-.kn)= P(X[k,]≤x,"..X[kn]≤xn)二、随机信号的数字特征均值m.[k]= E(X[k]方差α"[K]= E((X[K]- m,[k])?} = E(X?[k]} -m[k]自相关函数R[k,k]= E(X[k]X[k,]互相关函数R,[k],k,]= E(X[k,]Y[k,]]三、平稳各态遍历随机信号的时域描述1平稳随机序列指统计特性不随时间的平移而变化的那一类随机序列严平稳随机序列：F(xj,2x,2,..)=F(,x2,xn,,+n, +n,..+n)宽平稳随机序列：E(X[K]) =mE(X[K]X[k+n]} = R,[n]平稳随机信号自相关函数特性(1)对称性R.[n]= R,[-n]R [n]=R[-n]

第 6 章 功率谱估计 ◼ 离散随机序列的特征描述 ◼ 平稳随机序列通过 LTI 系统 ◼ 经典功率谱估计 ◼ 现代功率谱估计 6.1 离散随机序列的特征描述 ◼ 随机过程的分布函数 ◼ 随机信号的数字特征 ◼ 平稳各态遍历随机信号的时域描述 ◼ 平稳各态遍历随机信号的频域描述（功率谱密度） 一、 随机过程的分布函数 {X[k], kZ}表示一个随机过程 一维分布函数 二维分布函数 N 维分布函数 二、随机信号的数字特征 均值 方差 自相关函数 互相关函数 三、平稳各态遍历随机信号的时域描述 1 平稳随机序列 指统计特性不随时间的平移而变化的那一类随机序列 严平稳随机序列： 宽平稳随机序列： mx E{X[k]} = E{X[k]X[k n]} R [n] + = x 平稳随机信号自相关函数特性 （1） 对称性 R [n] R [ n] x = x − [ ] [ ] * Rx n = Rx −n ( , , , ; , , ) ( , , , ; , , ) 1 2 1 2 1 2 1 2 F x x x k k k F x x x k n k n k n N N N N   =  + +  + F(x, k) = P(X[k]  x) ( , ; , ) ( [ ] , [ ] ) F x1 x2 k1 k2 = P X k1  x1 X k2  x2 ( , , , ; , , ) F x1 x2  xN k1 k2 kN ( [ ] , [ ] ) = P X k1  x1 X kN  xN m [k] E{X[k]} x = [ ] {( [ ] [ ]) } { [ ]} [ ] 2 2 2 2 k E X k m k E X k m k  x = − x = − x [ , ] { [ ] [ ]} 1 2 1 2 R k k E X k X k x = [ , ] { [ ] [ ]} 1 2 1 2 R k k E X k Y k xy =

（2）极限值n=0R,[O]=E(X-[k]]R,[]=m?n-→（3）不等式R,[0]≥ R.[n]2.各态遍历随机信号集平均等于时间平均1m, = E(X[k]) = limx[k]m2N+1kN1? = E([X[k]-m,]} = lim[x[k]-m]Nm2N+1k=NN12x[kx[k+n]R.[n]==E(Xk]X[k+n=lim2N+1kN四、平稳各态遍历随机信号的频域描述功率谱密度1[2N +1/Fx(2, N))P(2) = lim El-->a维纳一一辛钦公式P(Q)= ZR(n)e-mn[P(2)ej"dR(n)=2元当自相关函数绝对可积时，平稳随机信号的自相关函数和功率谱密度是一对傅里叶变换对

（2）极限值 n = 0 [0] { [ ]} 2 R E X k x = n → 2 [ ] Rx  = mx （3）不等式 R [0] R [n] x  x 2. 各态遍历随机信号 集平均等于时间平均 =− → + = = N k N N x x k N m E X k [ ] 2 1 1 { [ ]} lim 四、平稳各态遍历随机信号的频域描述 功率谱密度 维纳——辛钦公式   =− −   = n j n x x P ( ) R (n)e =    − R n P e d j n x x ( ) 2 1 ( )    当自相关函数绝对可积时，平稳随机信号的自相关函数和 功率谱密度是一对傅里叶变换对。 ( , ) ] 2 1 1 ( ) lim [ 2 F N N P E X N x  +  = → [ ] [ ] 2 1 1 [ ] { [ ] [ ]} lim x k x k n N R n E X k X k n N k N N x + + == + = =− →  =− → − + = − = N k N x N x x x k m N E X k m 2 2 2 [ [ ] ] 2 1 1  {[ [ ] ] } lim

6.2平稳随机序列通过LTI离散时间系统输出序列的均值输出序列的自相关函数输出序列的功率谱心输入/输出序列的互相关函数及互功率谱平稳随机序列通过LTI系统Y[Kk]h[k]X[K]Y[k]= h[n]X[k -n]n、输出序列的均值m,[k]= E([k]) = h[n]E(x[k-n]]=m,h[n]H(ejo)m,[K] = m,H(ej0)二、输出序列的自相关函数R,[n]= R,[n]* R.[n]Ry[n]是系统单位脉冲响应h[k]的自相关函数Rh[n]与输入随机序列X[K]的自相关函数 Rx[n]的卷积.系统单位脉冲响应h[K]是确定信号，其自相关函数定义为= DTFT(R,[n]}-DTFT(R,[n])明证R,[n] = E(y[k]y[k +n])= E( h[]x[k -1](Zh[m]x[k + n - m])m-ZZ h[]h[m]E(x[k-]x[k + n -m]R,[n-m+]] [][n+]* R,[n] = R,[n]* R,[n]R,[n]

6.2 平稳随机序列通过 LTI 离散时间系统 ◼ 输出序列的均值 ◼ 输出序列的自相关函数 ◼ 输出序列的功率谱 ◼ 输入/输出序列的互相关函数及互功率谱 平稳随机序列通过 LTI 系统 X[k] h[k] Y[k] Y[k] h[n]X[k n] n =  − 一、输出序列的均值 H(ej0) 二、输出序列的自相关函数 R [n] R [n] R [n] y = h  x Ry[n]是系统单位脉冲响应 h[k]的自相关函数 Rh[n]与输入随机序列 X[k] 的自相关函数 Rx[n]的卷积. 系统单位脉冲响应 h[k]是确定信号，其自相关函数定义为 DTFT{R [n]} DTFT{R [n]} h x =  Ry [n] = E{y[k]y[k + n]}       = ( [ ] [ − ])( [ ] [ + − ]) l m E h l x k l h m x k n m h[l]h[m]Ex[k l]x[k n m] l m =  − + − R [n m l] x − + h[l]h[n l]* R [n] x l =  + R [n] R [n] = h  x Rh [n] 证 明 m [k] E{y[k]} h[n]E{x[k n]} n y = =  − =  n mx h[n] [ ] ( ) j0 m k m H e y = x

三、输出序列的功率谱P,(2)= DTFT(R,[n]} = DTFT(R,[n]* R,[n])DFTRDIFTRDTFT(h[n)* h- n]) =|H(ej)2P(2)P,(2)=H(e) P(2)四、输入/输出序列的互相关函数及互功率谱互相关R,[n]= E(Y[k]X[k +n]) = h[-n]* R,[n]R,[n]= E(x[k]y[n+ k]] = h[n]* R[n]互功率谱P(2)= DIFT(R,[nl) = H(ej)P(2)Pr(2)=DIFT(R,[nl) = H*(e/°)P(2)[例]一离散时间平稳白噪声通过一阶IIR数字滤波器[α[alH(2)=1 az-lh[k]=α*u[k]1H(ejn)ae-oH(en1-2αcosQ+α2（1）计算输出的自相关函数

三、输出序列的功率谱 P ( ) DTFT{R [n]} DTFT{R [n] R [n]} y  = y = h  x ( ) ( ) ( ) 2  =   x j Py H e P DTFT{R [n]} DTFT{R [n]} h x =  2 DTFT{ [ ]* ]} ( )  − = j h n h n H e Px() 四、输入/输出序列的互相关函数及互功率谱互相关 R [n] E{Y[k]X[k n]} h[ n] R [n] yx = + = −  x R [n] E{x[k]y[n k]} h[n]*R [n] xy = + = x 互功率谱 () = DTFT{ [ ]} = ( ) ()  x j Pxy Rxy n H e P () = DTFT{ [ ]} = *( ) ()  x j Pyx Ryx n H e P [例]一离散时间平稳白噪声通过一阶 IIR 数字滤波器 求输出的自相关函数、平均功率和功率谱。 零均值白噪声的特征 E{X[k]} = 0 2 P() = d [ ] 2π 1 [ ] 2 j 2 π π R n e Ω n n x =  =   − 解： z a az H z  − = −1 1 1 ( ) h[k] u[k] k = −   − = j j e H e 1  1 ( ) 2 2 1 2 cos 1 ( ) −   + = j H e （1）计算输出的自相关函数 y[k] = x[k]+y[k −1]  1

α"SWαta"n≥01-αR,[n]=h[k]h[k +n]-Wα""k=0α*α+kn<o1-α2Qll1-α2ao0lR,[n]= R,[n]* R,[n] =α?R,[n]1-α2（2）输出平均功率1个ay2[K]= R,[0] :lim1-α3N2N+1k-N（3）输出功率谱92P,()=H(en) P(2)=1-2αcos +α经典功率谱估计6.3谱估计的质量--相关法(间接法)-周期图法(直接法)-周期图法的改进-利用MATLAB实现功率谱估计一、谱估计的质量1.估计量的偏差bia-E-02．估计量的方差var(0) = E((0- E(0)))若bia(@)=0，limvar(①)=0，则称为的一致估计3.估计量的均方差MSE(0)=E((0-0)1 =var(0)+bia(0)二、相关法（间接法)进行功率谱估计-相关法的理论基础一自相关函数估计的计算-相关法进行功率谱估计

=  + =  =− R [n] h[k]h[k n] k h         − =  − = − +  =− +  =   0 1 0 1 2 2 0 n n n k n k k n n k n k k         2 1   − = n R [n] R [n] R [n] y = h  x [ ] 2 = Rh n 2 2 1    − = n （2）输出平均功率 [ ] [0] 2 1 1 lim 2 y N k N N y k R N = + =− → 2 2 1   − = （3）输出功率谱 ( ) ( ) ( ) 2  =   x j Py H e P 2 2 1 2 cos   −  + = 6.3 经典功率谱估计 ◼ 谱估计的质量 ◼ 相关法(间接法) ◼ 周期图法(直接法) ◼ 周期图法的改进 ◼ 利用 MATLAB 实现功率谱估计 一、谱估计的质量 1．估计量的偏差  = }− ˆ } { ˆ bia{ E 2．估计量的方差 }) } ˆ { ˆ } {( ˆ var{ 2  = E  − E  若bia{ ˆ } = 0, lim var{ ˆ } = 0, 则称 ˆ为的一致估计 N→ 3．估计量的均方差 } ˆ } bia{ ˆ ) } var{ ˆ } {( ˆ MSE{ 2  = E  − =  +  二、相关法(间接法)进行功率谱估计 ◼ 相关法的理论基础 ◼ 自相关函数估计的计算 ◼ 相关法进行功率谱估计

M功率谱估计的质量1.维纳一辛钦定理R,[n]<FP(2)计算方法：（1）由随机序列一个样本的N个观测值计算自相关函数的估计（2）对R[n|进行DTFT即得该随机序列的功率谱估计P(Q)2.自相关函数估计的计算X[K]是宽平稳各态遍历随机信号，x[K]是其一个样本1R,[n]= lim x[k]x[k+n]N-→2N+1k=-N已知x[K]的N个观测值x[0],x[1].x[N-1]，则自相关函数的估计为1HR,[n]=x[k]x[k+n]N1x[-n]*x[n] -(N-1)≤n≤N-1N[例]已知平稳各态遍历的实随机序列XIK]的单一样本的N个观测值为x[K]=(1,0，-1]，试计算该随机序列的自相关函数估计。解：R,[n] =[-1)*x[ =(-1,0, 2,0, -1)c利用MATLAB计算相关函数的估计1.利用conv函数计算2.利用数字处理工具箱中提供的函数xcorrxcorr(xy);%随机序列X和Y的互相关xcorr(x)：%随机序列X的自互相关利用DFT计算自相关函数的估计1)对x[]补零形成L点序列x,[k](L≥2N-1)2) X,[m]= DFT(x,[K])

◼ 功率谱估计的质量 1. 维纳—辛钦定理 [ ]⎯→ () x F Rx n P 计算方法： （1） 由随机序列一个样本的 N 个观测值计算自相关函数的估计 （2） 对 [ ] ˆ Rx n 进行 DTFT 即得该随机序列的功率谱估计 ( ) P ˆ x  2. 自相关函数估计的计算 X[k]是宽平稳各态遍历随机信号，x[k]是其一个样本 =− → + + = N k N N x x k x k n N R n [ ] [ ] 2 1 1 [ ] lim 已知 x[k]的 N 个观测值 x[0],x[1],,x[N-1]，则自相关函数 的估计为 [ ] [ ] 1 [ ] ˆ 1 0 x k x k n N R n N k x =  + − = [ ]* [ ] 1 x n x n N = − − (N −1)  n  N −1 [例]已知平稳各态遍历的实随机序列 X[k]的单一样本的 N 个观 测值为 x[k]={1, 0，-1}，试计算该随机序列的自相关函数估计。 解: [ ]* [ ] 1 [ ] ˆ x n x n N Rx n = − { 1, 0, 2, 0, 1} 3 1 = − −  利用 MATLAB 计算相关函数的估计 1. 利用 conv 函数计算 2. 利用数字处理工具箱中提供的函数 xcorr xcorr(x,y); %随机序列 X 和 Y 的互相关 xcorr(x); %随机序列 X 的自互相关 利用 DFT 计算自相关函数的估计 1) 对 x[k]补零形成 L 点序列 x [k](L  2N −1) L 2) X [m] DFT{x [k]} L = L 1 0 −1 −1 0 1 −1 0 1 1 −1 0 0 0 0

IDFT(X,[m]3) R[n]=-3.相关法进行功率谱估计估计=R,[n]DTFT,DFT=P(O),P,[m)x[k]-1 N-1-/alR,[n =ktNZR,[ne-mn,P(2)=[≤N-1n=-LN1Wix[K]x[k+] ≤N-1R,[n]=Nx[K]N100≤n≤N-1x[k+n], n≥0N-1-n-nR.[n]=x[k]x[k +n]N=0x[k+n], n≤0N-1-n-n-(N-1)≤n≤01 N!N-1+1R,[n]=2x[]x/-n]Zx[k]x[k+n]=NNks-n1=0[例]已知实平稳随机序列XIK]单一样本的N个观测值为x[k]=(1, 0,-1,试利用相关法估计其功率谱。解：X[K]的自相关函数估计值为R.[n]=m*[=(-1,0,2,0. -1)

3) IDFT{ [ ] } 1 [ ] ˆ 2 X m N Rx n = L 3. 相关法进行功率谱估计 [ ] ˆ ( ), DTFT,DFT ˆ [ ] ˆ x[k]⎯估计⎯⎯→Rx n ⎯⎯⎯⎯⎯⎯→Px  Px m  − − =  = + N n k x x k x k n N R n 1 0 [ ] [ ] 1 [ ] [ ] , 1 ˆ ( ) ˆ  =  − −  =− P R n e L N jn L n L x x [ ] [ ] 1 1 [ ] 1 0 =  +  − − − = x k x k n n N N R n N n k x 0  n  N −1 [ ] [ ] 1 [ ] ˆ 1 0 x k x k n N R n N n k x =  + − − = − (N −1)  n  0   − + = − =− = + = − N n l N k n x x l x l n N x k x k n N R n 1 0 1 [ ] [ ] 1 [ ] [ ] `1 [ ] ˆ [例]已知实平稳随机序列 X[k]单一样本的 N 个观测值为 x[k]={1, 0，-1}， 试利用相关法估计其功率谱。 解：X[k]的自相关函数估计值为 { 1, 0, 2, 0, 1} 3 1 [ ]* [ ] 1 [ ] ˆ = − = − −  x n x n N Rx n 0 N−1 x[k] −n N−1−n x[k+n], n  0 N−1−n x[k+n], n  0 −n

对R[n]进行傅里叶变换得XIk]的功率谱估计4Px(Q)= DTFT(R,[n]) = (-e/20 +2-e-j20 ) -(1-cos22)224.相关法功率谱估计的质量功率谱估计的质量与自相关函数估计的质量密切相关R,mbia(R,[n]) =L1E (R'[]+ R[r+n]R[r-n])var(R[n] ~ -NN→o0，偏差、方差趋于零，是一致估计。N固定时，n→N，偏差、方差较大三、周期图法(直接法)进行功率谱估计周期图法功率谱估计的计算-周期图法功率谱估计的质量1.周期图法功率谱估计的计算[x[k]k=0,1,...,N-1已知x[K]=lo其它方法基础：1R,[n]=x[-n*[n]由维纳一辛钦定理IDTFT (x[n)*×[-n])P(Q) =N1X(e/)Xx(e-1n)=x(en)NN周期图法功率谱估计的步骤功率谱估计P:(2)=区(ein)DIFT→X(ej))x~[k]-NIx[K] DET→X[m]功率进注 P-[m=X[m入!Exx[k]e-ink其中 X(ej)=DTFT(x[K])=k=0N!-j2元mkZxw[k]e""X[m]=DFT(x[K]} = >k=0

对 [ ] ˆ Rx n 进行傅里叶变换得 X[k]的功率谱估计 [ ]} ˆ Px () = DTFT{Rx n  { 2 } 3 1 2 − 2 = − + − j j e e (1 cos 2 ) 3 2 = −  4. 相关法功率谱估计的质量 功率谱估计的质量与自相关函数估计的 质量密切相关 bia{ [ ]} R [n] N n Rx n = − x  ( [ ] [ ] [ ]) 1 [ ]} ˆ var{ 2 R r R r n R r n N R n r   + + −  =− N→，偏差、方差趋于零，是一致估计。 N 固定时，n→ N，偏差、方差较大 三、周期图法(直接法)进行功率谱估计 ◼ 周期图法功率谱估计的计算 ◼ 周期图法功率谱估计的质量 1. 周期图法功率谱估计的计算 已知:    = − = 0 其它 [ ] 0,1, , 1 [ ] x k k N x k N  方法基础: [ ]* [ ] 1 [ ] ˆ x n x n N Rx n = N − N 由维纳—辛钦定理 DTFT{ [ ]* [ ]} 1 ( ) ˆ x n x n N Px  = N N − ( ) ( ) 1  −  = j N j N X e X e N 2 ( ) 1  = j N X e N 周期图法功率谱估计的步骤 2 DTFT ( ) 1 [ ] ( ) ( )    ⎯⎯⎯→ ⎯⎯⎯⎯→  = j N x j N N X e N x k X e P 功率谱估计 DFT 2 [ ] 1 [ ] [ ] [ ] X m N x k X m P m N x N ⎯ ⎯→ N ⎯⎯⎯⎯→ =  功率谱估计 其中  − =  −  = = 1 0 ( ) DTFT{ [ ]} [ ] N k j k N N j N X e x k x k e  − = − = = 1 0 2 [ ] DFT{ [ ]} [ ] N k mk N j N N N X m x k x k e 

[例]已知实平稳随机序列X[]单一样本的N个观测值为x[K]=(1,0，-1]，试利用周期图法估计其功率谱。解：对x[K]进行离散时间傅里叶变换N-X~(e/)-Zx[k]e-10k =1-e-/2k=0功率谱估计为：X(ejn)?1 Xn(en)xx(ein)I N(2) =NIN-(1-e-j29)(1-e/22)3平稳高斯白噪声功率谱估计结果（周期图法）N=64N=12820(apea(ap)leaaC-2020asamod40-403FrequencyFrequencyN=256N=5122020(ep)(ap)le)1eerw005-2020amodamoe1040LV3L4FrequencyFrequency2.周期图法功率谱估计的质量 R[n]e-jon均值EI(Q)=n=(N-I)NNo0-→，E(IN(W))=Px(W))，渐进无偏估计(sin( N2))var(1(2)) =0*(1+方差Nsin QN增加，方差不减小，不是一致估计四、周期图法的改进-问题的提出平滑周期图(Blackman-Tukey法)-平均周期图法(Bartlett法)■-重叠平均周期图法(Welch法)1.问题的提出周期图法进行功率谱估计，方差不随N的增加减小

[例]已知实平稳随机序列 X[k]单一样本的 N 个观测值为 x[k]={1, 0，-1}，试利用周期图法估计其功率谱。 解：对 x[k]进行离散时间傅里叶变换 −  − =  −  =  = − 2 1 0 ( ) [ ] 1 j N k j j k N X e x k e e 功率谱估计为： ( ) ( ) 1 ( ) 1 ( ) j * j 2 j    = = X e X e N X e N I N N N N (1 )(1 ) 3 1 − 2 2 = − − j j e e 平稳高斯白噪声功率谱估计结果(周期图法) 2. 周期图法功率谱估计的质量 均值 j n N n N N R n e N N n E I −  − =− − − { ()} =  [ ] 1 ( 1) N→，E{IN(W)}= Px(W)}，渐进无偏估计 方差 } sin sin( ) var{ ( )} {1 2 4          = + N N I N  N 增加， 方差不减小，不是一致估计 四、周期图法的改进 ◼ 问题的提出 ◼ 平滑周期图(Blackman-Tukey 法) ◼ 平均周期图法(Bartlett 法) ◼ 重叠平均周期图法(Welch 法) 1. 问题的提出 周期图法进行功率谱估计, 方差不随 N 的增加减小

如何提高谱估计质量？减小方差方法：1.对自相关函数估计值加窗2.将N个观测值分段，计算各段的周期图，再取平均2.平滑周期图(Blackman-Tukey法)对自相关函数估计值加窗，将误差较大的估计值截去Wn]R,[n]e-jna?PM(2)=n=(N1)窗函数w[nl(MKNM)满足下述条件[0≤w[n]≤w[0]=1Jn bia[I ~(2]平均周期图方差减小的代价之一是偏差增大

如何提高谱估计质量? 减小方差方法: 1. 对自相关函数估计值加窗 2. 将 N 个观测值分段, 计算各段的周期图, 再取平均 2. 平滑周期图(Blackman-Tukey 法) 对自相关函数估计值加窗, 将误差较大的估计值截去 −  − =− −  =  jn x N n N MP w n R [n]e ˆ ( ) [ ] 1 ( 1) 窗函数 w[n] (M<N)满足下述条件      =  = −    = w n n M w n w n n M w n w [ ] 0 [ ] [ ] 0 [ ] [0] 1 B-T 法进行功率谱估计的主要步骤 ◼ (1) 利用观测数据估计自相关序列。 ◼ (2) 对自相关函数估计值加窗。 ◼ (3) 计算加窗后自相关函数的 DTFT。 优点：PM(W)波动比 IN (W)小，可证是一致估计 缺点：降低了频率分辨率 3. 平均周期图法(Welch-Bartlett 法) 将随机序列 X[k]的 N 个观测值分成 A 段 x [k] = x[iM + k];i = 0,1, A−1;k = 0,1, M −1 i   第 i 段序列的周期图为 2 ( ) 1 ( )   = j M i M X e M I 平均周期图 ( ) 1 ( ) 1 0  =   − = i M A i A M I A P 平均周期图法估计质量 var{ ( )} 1 var{ ()} =  i M A M I A P A→, 方差为零，是一致估计 因为  − =− − −   = − 1 ( 1) bia[ ( )] [ ] N n N j n N R n e N n I 所以 bia[ ()]  bia[ ()] N i M I I 平均周期图方差减小的代价之一是偏差增大