现代信号处理05a

第五章、谱估计研究二阶平稳随机过程特征－功率谱密度－揭示随机过程中所隐含的周期及相邻的谱峰等有用信息。

则要用有限长的N 个样本数据去估计该平稳随机过程的功率谱密度－谱估计的方法。

此种估计是建立在时间平均的方法之上，并假定具有遍历性。

遍历性：随机平稳过程满足：（1）⎰∞→=TT dt t x Tt x E 0)(1lim )]([（2）⎰+=+∞→TT dt t x t x Tt x t x E 0)()(1lim )]()([ττ和自相关函数Γx (m)是正定的，即∑∑≥-Γijjxi aj i a 0)( 对所有{a j }。

经典谱估计－线性、非参数化方法：周期图法，相关图法等。

采用经典的傅里叶变换及窗口截断。

对长序列有良好估计。

现代谱估计－非线性、参数化方法：最大似然估计，最大熵法，AR 模型法（Auto Regressive method ），预测滤波器法，ARMA 模型（Auto Regressive Move Average ）等。

对短序列的估计精度高，与经典法相互补充。

是融合经典变换理论、统计估计理论、系统辨识、信息论、时间序列分析及计算方法等理论与技术－新学科。

应用广泛，发展迅速。

1807--Fourier Transform1890’s—Schuster {}N n x N ,..2,1,=--Periodogram 周期图22211)(ωωωωφjN N j j e x e x e x N---+++=若0ωω=处有正弦信号与噪声叠加，在0ωω=处周期图出现峰值1930-Wiener-Kinchine 定理，广义谐波分析。

定义自相关函数与功率谱，并证明两者互为傅氏变换。

1958-Blackman-Tukey-BT 法-由通信工程观点测量功率谱 观测数据－自相关估计－乘窗函数－傅氏变换-->功率谱 1965-Cooley-Tukey —FFT —重提周期图法，广泛使用的经典法 1967—现代谱估计－Burg －最大熵谱估计－地震－线性预测 1968-Parzen-自回归AR 谱估计法1971-Van Den Bos-证明最大熵谱估计与AR 谱估计法等价。

自此开始了现代谱估计的深入研究：MA, ARMA 模型等构成现代谱估计的参量法、还有非参量法1969Capon 最大似然谱估计 经典谱估计＋现代谱估计＝完整的谱估计理论 长数据序列－必需应用经典谱估计－性能良好短数据序列－现代谱估计－估计精度分辨度大大高于经典法。

因BT 法及周期图法－无限长序列开窗截断＝有限长序列。

数据开窗、自相关开窗－频域发生“泄漏”，即功率谱的主瓣能量泄漏到旁瓣中，导致弱信号的主瓣被强信号的旁瓣所湮没，造成谱的模糊。

改进窗函数－抑制旁瓣－损失谱的分辨度。

短序列－谱分辨度的极限＝1/序列长度。

现代谱估计适用于短序列的情况，非线性运算，发展很快，新方法层出不穷。

80年代，利用信息论的熵谱估计法，多谱（高阶谱）估计及多维谱估计等。

小波、分数维等。

§5－1、谱密度意义一、 能谱密度设x (t )是确定性的复连续信号，若其绝对可积或其能量有限，即：)15()(2---∞<=⎰∞∞-dt t x E则x (t )的连续傅氏变换存在，由下式给出：)25()2ex p()()(----=⎰∞∞-dt ft j t x f X π根据Parseval 能量定理，有：)35()()(22---==⎰⎰∞∞-∞∞-df f X dt t x E由上式可见，信号能量E 等于信号频谱模值平方2)(f X 在整个频域上的积分，故称)45()()(2--=f X f ε为信号的能谱密度。

当x (t )为广义平稳过程时，其能量通常是无限的，则需研究其功率的频域上的分布，即功率密度。

对于平稳随机过程，谱分析是采用自相关函数：[])55()(*)()(--+=Γτττx t x E xWiener-Kinchine 定理将自相关函数与功率谱密度联系起来：)65()()(---Γ=Φ⎰∞∞--ττωωτd e j x x)75()(21)(---Φ=Γ⎰∞∞-ωωπτωτd e j x x离散形式－{}10,-≤≤N n x N 为平稳、零均值序列；其自相关(协方差)函数为：[])85()()()(--+=Γk n x n x E k 若有：∑∞-∞=--∞<Γk k )95()(则功率谱密度为：∑∞-∞=---Γ=Φk jk e k )105()()(ωω )(ωΦ是以0对称，周期为2π。

反变换为：)115...(2,1,0,)(21)(-±±=Φ=Γ⎰-k d e k jk ωωπππω其自相关函数的估计由时间平均函数给出：)125(,)()(1)(ˆ1--≤+=Γ∑--=N kk j x j x Nk k N j N功率谱的估计为：∑∞-∞=---Γ=Φk jk NNe k )135()(ˆ)(ˆωω 若定义矩形窗⎩⎨⎧----≤≤=)145(,010,1)((else N n n d 即加窗截断为有限长序列，则有：)155()()()()(1)(ˆ--++=Γ∑∞-∞=n N k n x k n d n x n d Nk功率谱的估计可写成：⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛++⎪⎭⎫ ⎝⎛=Φ∑∑∞-∞=+-∞-∞=n k n j n jn N ek n x k n d e n x n d N)()()()()(1)(ˆωωω将m =n +k 代入上式，得：)165(,)(1)()(*1)()(1)(ˆ2---≤≤-==-=ΦπωπωωωωωωN N N N N N X NX X N X X N 式中)175()()(10---=∑-=-N n jn N e n x X ωω为数据序列{}10,-≤≤N n x N 的离散傅氏变换DFT 。

(5-13)式－间接法、相关图法 (5-16)式－直接法、周期图法§5－2、相关图法(Correlogram Method)根据Wiener-Kinchine 定理，先估计出有限长信号x (n ) 的自相关函数，即：)185..(1,)()(1)(ˆ1--≤+=Γ∑-=N m m n x n x Nm N n易见)(ˆm Γ是偶函数，其长度为2N -1.实际计算时，由于x (n )只有N 个观测值，则对于每个延迟m ，可用的数据只有N -1-m 个，故可将上式改写成：)195..(1,)()(1)(ˆ10--≤+=Γ∑--=N m m n x n x Nm mN n第二步，求)(ˆm Γ的DFT,得到x (n )的功率谱估计： ∑-=---≤Γ=ΦMMm jm N M e m )205..(1,)(ˆ)(ˆωω 由于功率谱是经自相关函数间接求出的，故称间接法。

平稳随机过程的自相关函数应为：)(ˆ)(m m Γ≠Γ )215..()()(1lim)(1-+=Γ∑-=∞→N n N m n x n x Nm显见，实际上自相关估计)(ˆm Γ是仅用有限个数据得到的。

现在，来讨论相关图法的性能，即)(ˆm Γ接近)(m Γ的程度。

估计偏差：先求)(ˆm Γ的均值（数学期望） )225)..(()]()([1)()(1)](ˆ[1010-Γ-=+=⎥⎦⎤⎢⎣⎡+=Γ∑∑--=--=m NmN m n x n x E Nm n x n x N E m E mN n mN n 所以估计的偏为：)235)..(()()](ˆ[)](ˆ[-Γ-=Γ-Γ=Γm Nm m m E m Bias 故可得出结论：m ≠ 0，为有偏估计。

当N →∞，m ＝有限，偏→0，故为渐进无偏估计。

估计方差：推导过程略。

(l =n -k )[]{}[])245.()()()(]1[1)()(ˆ)](ˆ[21)1(2--Γ+Γ+Γ+-=Γ-Γ=Γ∑-----=m i m i i N l m Nm m E m Var mN m N i 显然，当N →∞，0)](ˆ[=Γm Var ；又因为： )255..(0)]([lim -=Γ∞→m Var N故)(ˆm Γ为)(m Γ的渐进无偏一致估计。

显然，N 越大，估计精度越高。

当N ,m 均较大时，还可利用FFT 进行计算。

由自相关可通过式(5-20)计算功率谱。

M 越大，分辨率越高，但自相关的偏差及方差也相应增大，通常取M =N /10—N /5，较好。

§5－3、周期图法(Periodogram Method)定义：长度为N 的实平稳随机信号序列{}10,-≤≤N n x N的周期图为：)265(,)(1)(2---=ωωN N X NI 式中DFT en x X N n jn N -=∑-=-1)()(ωω即直接由序列的DFT 计算出来的。

由于x (n )的DFT 有周期性，所以)(ωN I 也有周期性。

周期图法谱估计的性能－估计的偏：)275..()()]([1)1(-Γ-=----=∑ωωjn X N N n N e m NmN I E 所以可见，)(ωN I 不是自相关的DFT ，故为有偏估计。

相当于无限长序列加长度为N 的矩形窗：)285..()()(21)]([--Φ=⎰∞∞θθωθπωd Q I E N N －式中)295()()(---Γ=Φ∑∞-∞=-k jk ek ωω)295..(2/sin 2/sin 1])([1)(22-⎥⎦⎤⎢⎣⎡==∑∞-∞=-ωωωωN N e n d N Q n jn N )(ωN Q 称为F eijer’s Kernel ，其主瓣宽度为2/N .可以证明，当N →∞时，)305..()2(2)(--=∑∞-∞=k N k Q πωπδω)315)...(()]([lim -Φ=∞→ωωN N I E 渐进无偏差估计方差:)325..(sin sin 1)]([24-⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛+=ωωσωN N I Var x N式中：2)]([x N I E σω∆=假定{})(n x 为零均值、正态分布白噪声序列。

可见当N →∞时，方差)]([ωN I Var 并不等于0。

所以，周期图法不是功率谱的一致估计。

这主要是由于时域序列加窗后，频域发生泄漏现象，所以人们作了许多改进，加各种形式的窗。

另外还有分段求估计再求平均的平均周期图法，平滑周期图法等。

比较好的窗有Blackman-Tukey 窗：时域窗：)335(.,,....0cos121)(---≤⎪⎩⎪⎨⎧⎥⎦⎤⎢⎣⎡+=M m elseM m m w π而谱窗：)345).((21)()(41)(-+⎥⎦⎤⎢⎣⎡++-=θπθπθθM M M N D M D M D W式中：2sin)21sin()(θθθ+=M D M 功率谱的Blackman-Tukey 估值为：)355).((ˆ21)(ˆ41)(ˆ41)(000-Φ++Φ+-Φ=ΦωπωπωωM M BT 式中：)365()(ˆ)()(ˆ0---Γ==Φ∑-=-MMm jm NNe m I ωωω 为得到功率谱密度的一致估计，谱窗应满足如下条件： （1）W N (ω)>=0,对所有的N ，ω（2）1)(21=⎰-ππωωπd WN对所有N（3）∞<⎰-ππωωπd W N )(212（4）W N (ω)→0，N →∞，上述条件可以是使W N (ω)作为算子，当N →∞时，为δ函数。