经典功率谱估计方法实现问题的研究

1 随机信号的经典谱估计方法估计功率谱密度的平滑周期图是一种计算简单的经典方法。

它的主要特点是与任何模型参数无关，是一类非参数化方法[4]。

它的主要问题是：由于假定信号的自相关函数在数据观测区以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。

在一般情况下，周期图的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。

本章主要介绍了周期图法、相关法谱估计（BT ）、巴特利特(Bartlett)平均周期图的方法和Welch 法这四种方法。

2.1 周期图法周期图法又称直接法。

它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样.周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。

只是1965年FFT 出现后，此法才变成谱估计的一个常用方法。

周期图法[5]包含了下列两条假设：1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。

这当然必然带来误差。

2．由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。

这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。

与相关法相比，相关法在求相关函数)(m R x 时将)(n x N 以外是数据全都看成零，因此相关法认为除)(n x N 外x(n)是全零序列，这种处理方法显然与周期图法不一样。

但是，当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。

通过比较我们发现：如果相关法中M=N ，不加延迟窗，那么就和补充（N-1）个零的周期图法一样了。

简单地可以这样说：周期图法是M=N 时相关法的特例。

因此相关法和周期图法可结合使用。

2.2 相关法谱估计（BT ）法这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。

它是1958年由Blackman 和Tukey 提出。

这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列)(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧序列。

即)()(1)(1m n x n xNm R N n N Nx +=∑-=∧(2-1)这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。

,M-1的傅里叶变换，另一半也就知道了。

第三步：由相关函数的傅式变换求功率谱。

即jwmM M m Xjwx em Re S ----=∧∧∑=)()(1)1( (2-2)以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。

因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。

一般取M<<N ，因为只有当M 较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。

因此，在FFT 问世之前，相关法是最常用的谱估计方法。

当FFT 问世后，情况有所变化。

因为截断后的)(n x N 可视作能量信号，由相关卷积定理可得)(m R x ∧)]()([1m x m x NN N -*=(2-3)这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。

我们可对上式两边取（2N-1）点DFT ，则有2121212)(1)()([1)(K X NK x k x Nm R N N N x ---∧=*=(2-4)于是将时域卷积变为频域乘积，用快速相关求)(m R x ∧的完整方案如下： 1.对N 长)(n x N 的补充（N-1）个零，成为（2N-1）长的。

2.求（2N-1）点的FFT ，得∑-=----=220121212)()(N N mkN N N W n xK X 。

3.求212)(1K X N N -。

由DFT 性质，)(12n x N -是纯实的，)(12k x N -满足共轭偶对称，而212)(1K X NN -一定是实偶的，且以（2N-1）为周期。

4. 求（2N-1）点的IFFT ：mkN N N k N x W K X NN m R -----=-∧∑-=121)1(212)(1121)( (2-5)这里212)(1K X NN -是实偶的，m=-(N-1)...0...N-1。

本来IFFT 求和范围是0至2N-2，由于212)(1K X NN -的实偶性与周期性，求和范围改为-（N-1）至（N-1）不影响计算结果。

同理可将m 的范围改为-（N-1）至（N-1）。

上述的快速相关中，补充零的目的是为了能用圆周卷积代替线性卷积，以便进一步采用快速卷积算法。

快速相关输出是-（N-1）至（N-1）的2N-1点，加)(m W M 窗后截取的是-（M-1）至（M-1）的频段，最后作（2M-1）点FFT ，得)(k S x ∧。

我们注意到：如果数据点数与自相关序列点数相同即M=N ，则（2N-1）点的IFFT 后紧跟一个（2N-1）点的FFT ，利用)(m R x ∧的对称性，FT 运算框的计算式变为=)(K S X ∑---=∧1)1()(N N m Xm R1212---N mkN X W (2-6)由于N=M 并假设窗形状是矩形的，第二次()m W M 的截断就不需要了。

比较式（2-5）和式（2-6）， ](m)R FFT[k S x x∧=）（,])(1[)( 212K X NIFFT m R N x -∧=正反傅氏变换可以抵消，直接得）（k S x=212)(1K X NN - (2-7)为了实行基2FFT ，也可将（2N-1）点换成2N 点，这样做不影响结果的正确性。

2.3 巴特利特(Bartlett)平均周期图法首先让我们来看一下为什么周期图经过某种平均(或平滑)后会使它的方差当∞→N 时趋于零，达到一致估计的目的。

如果L x x x , , ,21 是不相关的随机变量，每一个具有期望值μ，方差2σ，则可以证明它们的数学平均L x x x x l /)...(21+++=的期望值等于μ，数学平均的方差等于L /2σ，即：[][]μμ=⋅=+++=L Lx x x E Lx E L 1121[][][]222)(][))((x E xE x E x E x Var -=-=[]22212)(1μ-+++=L x x x E L()[]2112222121μ-⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧++++=∑∑=≠=L j Lji i j i L x x E x x x E L[][][]∑∑∑∑=≠==≠=⋅=Lj Lji i ijL j ji i jix E x E x x E 1111222)1(μμμμL L L L -=-=所以[][][]⎥⎦⎤⎢⎣⎡-=-⎥⎦⎤⎢⎣⎡-+=∑∑==2122222122211μμμμL x E L L L x E L x Var L i i L i i [][][][][][]{}22222221212])[(])[(])[(1LLxE x E xE x E xE x E L-++-+-=LL L2221σσ==(2-8)由(2-8)可见，L 个平均的方差比每个随机变量的单独方差小L 倍。

当[]0 →∞→x V a r L ，可达到一致谱估计的目的。

因而降低估计量的方差的一种有效方法是将若干个独立估计值进行平均。

把这种方法应用于谱估计通常归功于Bartlett 。

Bartlett 平均周期图的方法是将序列)10( )(-≤≤N n n x 分段求周期图再平均。

设将)(n x 分成L 段，每段有M 个样本，因而LM N =，第i 段样本序列可写成Li M n M iM n x n x i≤≤-≤≤-+=1 ,10 )()(第i 段的周期图为21)(1)(∑-=-=M n nj jiM en xMI ωω如果)( ,m M m xx φ>很小，则可假定各段的周期图)(ωiM I 是互相独立的。

对功率谱密度的概念的讨论，谱估计可定义为L 段周期图的平均，即∑==Li i MxxILP 1)(1)(ˆωω (2-9)于是它的期望值为[][][]∑===Li iM i MxxI E IE LP E 1)()(1)(ˆωωω[]⎰--=ππθωθθπωd eW P P E j B xx xx)()(21)(ˆ)(()[()[⎰-⎥⎦⎤⎢⎣⎡--=ππθθωθωθπd M P Mxx 22/sin 2sin )(21 (2-10)这里L N M /=，因此Bartlett 估计的期望值是真实谱)(ωxx P 与三角窗函数的卷积。

由于三角窗函数不等于δ函数，所以Bartlett 估计也是有偏估计即0≠Bias ，但当∞→N 时，0→Bias。

由于我们假定各段周期图是相互独立的，所以可按式(2-8)得到下式：[][])(1)(ˆωωMxx I Var LP Var =⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛+≈22s i n )s i n (1)(1ωωωM M P L xx (2-11)由此可见，随着L 的增加[])(ˆωxx P Var 是下降的，当∞→L 时，[]0)(ˆ→ωxxP Var 。

因此Bartlett 估计是一致估计。

比较式(2-10)的[])(ˆωxxP E 与式(2-1)的[])(ωN I E 可见在二种情况的估计量的期望值都是真值)(ωxx P 与窗口函数)(ωj B e W 的卷积形式，但后者将前者WB 中N 改为M ，NL N M <=/。

因而使)(ωj B e W 主瓣的宽度增大。

由于主瓣的宽度愈窄愈接近δ函数，偏倚愈小。

今式(3-10)中)(ωj B e W 的主瓣宽度大于后式中的主瓣宽度，因而[][])()(ˆωωNxx I Bias P Bias >，而主瓣愈宽分辨率就愈差。

因此Bias 可用来说明谱的分辨率，Bias 愈大说明谱分辨率愈差。

一个固定的记录长度N ，周期图分段的数目L 愈大将使方差愈小，但M 也愈小，因而使Bias 愈大，谱分辨率变得愈差。

因此Bartlett 方法中Bias 或谱分辨率和估计量的方差间是有互换关系的。

M 和N 的选择一般是由对所研究的信号的预先了解来指导的。

例如，如果我们知道谱有一个窄峰，同时如果分辨出这个峰是重要的，那么我们必须选择M 足够大。

又从方差的表达式我们可以确定谱估计的可接受的方差所要求的记录长度N=(LM)。