1 随机信号的经典谱估计方法
估计功率谱密度的平滑周期图是一种计算简单的经典方法。它的主要特点是与任
何模型参数无关,是一类非参数化方法[4]。它的主要问题是:由于假定信号的自相关函数在数据观测区以外等于零,因此估计出来的功率谱很难与信号的真实功率谱相匹配。在一般情况下,周期图的渐进性能无法给出实际功率谱的一个满意的近似,因而是一种低分辨率的谱估计方法。本章主要介绍了周期图法、相关法谱估计(BT)、巴特利特(Bartlett)平均周期图的方法和Welch法这四种方法。 2.1 周期图法
周期图法又称直接法。它是从随机信号x(n)中截取N长的一段,把它视为能量有限x(n)真实功率谱Sx(ejw)的估计Sx(ejw)的抽样.
周期图这一概念早在1899年就提出了,但由于点数N一般比较大,该方法的计算量过大而在当时无法使用。只是1965年FFT出现后,此法才变成谱估计的一个常用方法。周期图法[5]包含了下列两条假设:
1.认为随机序列是广义平稳且各态遍历的,可以用其一个样本x(n)中的一段
xN(n)来估计该随机序列的功率谱。这当然必然带来误差。
2.由于对xN(n)采用DFT,就默认xN(n)在时域是周期的,以及xN(k)在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段xN(n)的周期延拓,这也就是周期图法这个名字的来历。与相关法相比,相关法在求相关函数Rx(m)时将
xN(n)以外是数据全都看成零,因此相关法认为除xN(n)外x(n)是全零序列,这种处理方法显然与周期图法不一样。
但是,当相关法被引入基于FFT的快速相关后,相关法和周期图法开始融合。通过比较我们发现:如果相关法中M=N,不加延迟窗,那么就和补充(N-1)个零的周期图法一样了。简单地可以这样说:周期图法是M=N时相关法的特例。因此相关法和周期图法可结合使用。
2.2 相关法谱估计(BT)法
1
这种方法以相关函数为媒介来计算功率谱,所以又叫间接法。它是1958年由Blackman 和Tukey提出。这种方法的具体步骤是:
第一步:从无限长随机序列x(n)中截取长度N的有限长序列列xN(n) 第二步:由N长序列xN(n)求(2M-1)点的自相关函数Rx(m)序列。即
?1N?1Rx(m)??xN(n)xN(n?m) (2-1)
Nn?0?这里,m=-(M-1)…,-1,0,1…,M-1,MN,Rx(m)是双边序列,但是由自相关函数的偶对称性式,只要求出m=0,。。。,M-1的傅里叶变换,另一半也就知道了。
第三步:由相关函数的傅式变换求功率谱。即
?M?1?Xm??(M?1)Sx(e)?jw?R(m)e?jwm (2-2)
以上过程中经历了两次截断,一次是将x(n)截成N长,称为加数据窗,一次是将x(n)截成(2M-1)长,称为加延迟窗。因此所得的功率谱仅是近似值,也叫谱估计,式中的Sx(ejw)代表估值。一般取M< 当FFT问世后,情况有所变化。因为截断后的xN(n)可视作能量信号,由相关卷积定理可得 Rx(m)??1[xN(m)?xN(?m)] (2-3) N这就将相关化为线性卷积,而线性卷积又可以用快速卷积来实现。我们可对上式两边取(2N-1)点DFT,则有 Rx(m)??112[x2N?1(k)?x2N?1(K)?X2N?1(K) (2-4) NN2 ?于是将时域卷积变为频域乘积,用快速相关求Rx(m)的完整方案如下: 对N长xN(n)的补充(N-1)个零,成为(2N-1)长的。 2N?21. 2. 求(2N-1)点的FFT,得X2N?1(K)?求 N?0?x2N?1(n)W2?Nmk?1。 3. 12X2N?1(K)。由DFT性质,x2N?1(n)是纯实的,x2N?1(k)满足共轭偶对称,N12X2N?1(K)一定是实偶的,且以(2N-1)为周期。 而N4. 求(2N-1)点的IFFT: N?1112?mkRx(m)?X(K)W?2N?12N?1 (2-5) 2N?1k??(N?1)N?12X2N?1(K)是实偶的,m=-(N-1)...0...N-1。本来IFFT求和范围是0至2N-2,N12X2N?1(K)的实偶性与周期性,求和范围改为-(N-1)至(N-1)不影响计由于N这里 算结果。同理可将m的范围改为-(N-1)至(N-1)。 上述的快速相关中,补充零的目的是为了能用圆周卷积代替线性卷积,以便进一步采用快速卷积算法。快速相关输出是-(N-1)至(N-1)的2N-1点,加WM(m)窗后截取的是-(M-1)至(M-1)的频段,最后作(2M-1)点FFT,得Sx(k)。我们注意到:如果数据点数与自相关序列点数相同即M=N,则(2N-1)点的IFFT后紧跟一个(2N-1)点的FFT,利用Rx(m)的对称性,FT运算框的计算式变为 N?1???SX(K)?m??(N?1)?RX(m)W2?Nmk?1X2N?1 (2-6) 由于N=M并假设窗形状是矩形的,第二次WM?m?的截断就不需要了。比较式 ?FFT[Rx(m)] , Rx(m)?IFFT[(2-5)和式(2-6),S(xk)??12X2N?1(K)]正反傅氏N变换可以抵消,直接得 3 ?S(xk)12X2N?1(K) (2-7) N为了实行基2FFT,也可将(2N-1)点换成2N点,这样做不影响结果的正确性。 2.3 巴特利特(Bartlett)平均周期图法 首先让我们来看一下为什么周期图经过某种平均(或平滑)后会使它的方差当 N??时趋于零,达到一致估计的目的。如果x1, x2, ?, xL是不相关的随机变量, 每一个具有期望值?,方差?2,则可以证明它们的数学平均x?(x1?x2?...?xl)/L的期望值等于?,数学平均的方差等于?2/L,即: E?x??11E?x1?x2???xL???L??? LLVar?x??E(x?E(x))2?E[x2]?(E?x?)2 ?122E(x?x???x)?? 12L2L?????LL1?222?2?Ex1?x2???xL???ExixjL?j?1i?1i?j???????2??? ????E?xx???E?x???E?x? ijjij?1i?1i?jj?1i?1i?jLLL?L?(L?1)??L2?2?L?2 所以 1?L1?L?2222?2Var?x??2??Exi?L??L?????2??Exi2?L?2? L?i?1L?i?1?????? ?122Ex12?(E[x1])2?Ex2?(E[x2])2???ExL?(E[xL])2 2L??????????????1?22?2L?? (2-8) LL由(2-8)可见,L个平均的方差比每个随机变量的单独方差小L倍。当 L?? Var?x??0,可达到一致谱估计的目的。因而降低估计量的方差的一种有效 4 方法是将若干个独立估计值进行平均。把这种方法应用于谱估计通常归功于Bartlett。 Bartlett平均周期图的方法是将序列x(n) (0?n?N?1)分段求周期图再平均。设将x(n)分成L段,每段有M个样本,因而N?LM,第i段样本序列可写成 xi(n)?x(n?iM?M) 0?n?M?1, 1?i?L 第i段的周期图为 1iIM(?)?MM?1n?02j?x(n)e?j?n i如果m?M, ?xx(m)很小,则可假定各段的周期图IM对功率(?)是互相独立的。 谱密度的概念的讨论,谱估计可定义为L段周期图的平均,即 1Li?Pxx(?)??IM(?) (2-9) Li?1 于是它的期望值为 1Lii?EPxx(?)??EIM(?)?EIM(?) Li?1???????(?)?1EPxx2???????Pxx(?)WB(ej(???))d? ?21?2?M?sin??????M2?P(?)?d? (2-10) ???xx????sin???/2??这里M?N/L,因此Bartlett估计的期望值是真实谱Pxx(?)与三角窗函数的卷积。由于三角窗函数不等于?函数,所以Bartlett估计也是有偏估计即Bias?0,但当N??时,Bias?0。 由于我们假定各段周期图是相互独立的,所以可按式(2-8)得到下式: ?(?)?1Var?I(?)? VarPxxML5 ??