基于FFT与自相关函数的快速功率谱估计方法＊

李春林 伍 勇

（92755部队67分队1） 临高 571820）（92098部队电子对抗科2） 陵水 572425）

式中，E［·］表示求期望值。

根据以上式子，可以得到有限数据的功率谱为

1 引言

在无线电技术（比如通信、雷达、电子战、遥控遥测等）领域，无线电系统所处理的主要对象是电磁波［1～2］。无线电系统中的接收天线首先把电磁波转化为电信号，然后馈入到接收系统进行必要的分析和处理。在20世纪70年代以前，对电磁信号的处理主要以模拟处理方法为主［3～5］，包含有滤波、放大、混频、检波等等各种模拟处理环节。模拟处理方法不仅缺乏灵活性、可扩展性，而且每一环节都会不同程度地引入各种非线性失真，导致接收系统性能下降，功能降低。

自从1965年J.W.Tukey和T.W.Coody在《计算数学》杂志上发表了著名的《机器计算傅立叶级数的一种算法》［6］论文并几经人们改进之后，很快形成了一套高效的算法FFT。80年代以后，随着微电子技术的迅猛发展，基于FFT的数字信号处理技术开始获得广泛应用，并逐步显现出模拟处理方式所无法达到的优越性［7～8］。

本文将引入FFT用于基于自相关函数的功率谱估计中，加快估计的计算速度。在信号功率谱估计中有许多的高分辨率谱估计算法，之所以讨论FFT在关于谱估计中的应用，是因为FFT已经广泛应用到数字接收机的设计中。文中首先简要介绍了有偏和无偏自相关函数的定义，然后给出了基于自相关函数的谱估计算法，最后就如何借助FFT来实现信号的功率谱估计给出了详细过程，并给出了计算实例。

2 自相关函数的定义及其计算过程

假设有N点输入数据x（n），n＝0，…，N－1，其自相关定义为［1，9］：

式中m称之为自相关的延迟变量。严格意义下，上述自相关函数应该称为取样自相关，它逼近于E［x（n）x（n＋m）］，式中E［·］表示期望值。m值既可以正也可以负。如果变量m为负数，那么其自相关函数与正m的自相关函数的关系为R（－m）＝R（m）＊，其中R（k）是复数，如果R（k）是实数，则R（－m）＝R（m）。如果延迟变量是m，那么把0～N－1的输入数据分成长度相同的两组，一组从0～N－M－1，另一组从m～N－1。这两组中数据项相乘的情况如0所示。所有乘积项的和等于自相关函数R（m）的N倍。

图1 R（m）的计算过程描述

自相关可以视为两组数据之间相似性的度量。如果两组数据很相似，那么其自相关函数就大，相反就小。当m＝0时，两组数据是相同的，所以在所有自相关函数值中R（0）最大。当m值比较大时，就会只有很少几项参加求和运算，但式（1）中的分母N是一个固定常数，所以求和式被N除之后，通常使R（m）的幅度变得很小，而理论上该值可能是比较大的。因此，式（1）所定义的自相关函数称为自相关函数的有偏形式。

而无偏自相关函数的定义为

在这一定义中，当m变大时，分母变小，而且等于求和的项数。由于这种形式会产生负的功率谱，所以很少用在谱估计中。如果R（m）比较小，那么对所计算的功率谱的影响就会小一些。由于R（m）只有很少几个点，所以预计R（m）值会比较小，这样对功率谱的影响也就小了。

3 基于自相关函数的功率谱估计的结构

可以采用两种方法来估计功率谱。第一种方法就是计算输入信号的FFT，然后对其结果求平方，这一方法叫周期图法。这一结果可以直接从FFT运算得到。第二种方法，功率谱是从输入数据的自相关函数得到的，可以表示为［9］

式中，E［·］表示求期望值。

根据以上式子，可以得到有限数据的功率谱为

式中，m取值区间为从－M～M，k是频率分量，ts是采样间隔。该方法通常称为Blackman－Tukey方法［11～12］。求和项的总数为2M＋1，因为求和式中包含了m＝0的项。

为了获得正的功率谱，在上述等式中通常采用有偏自相关R（m），可以把有偏自相关R（m）作为一个加窗函数。当m比较大时，R（m）通常就小。有时把求和项大约限制在m＝－N／10～N／10之间，推荐的最大值为－N／5～N／5之间。另外还可以对有偏自相关加一个窗函数，以进一步减小旁瓣。

加窗后的功率谱可以写为

式中，W（m）为窗函数。该窗函数必须是对称的，保证P（k）为偶函数。

可按如下方程求解R（m）的期望值。因为R

R（m）的期望值为

这是由于R（m）与n无关。求和式等于N－｜m｜，所以式（7）可以写为

式中，（N－｜m｜）／N表示一个三角窗。由此可见，该方法受固有加窗效应的内在限制，而m＞N时，R（m）的值假定为0。

4 基于FFT的功率谱估计步骤

用FFT来计算式（4）的结果。式（4）与DFT非常类似，但k值是任意的。比如，从以上方程可以计算出任意给定k值所对应的结果。如果对k严格加以限制，使其与DFT完全一致，那么就可以采用FFT算法来计算，从而节省计算时间。为此，必须把式（4）变为合适的形式。DFT可写为

式中，X（K）为频域响应，x（n）为时域取样数据点。在上述方程中，n和k是离散的，N通常取为2的幂次方。

对式（4）进行适当的变形，需采取以下几个步骤：

第一步：改变标记符t。通过比较式（4）和式（9）可以看出，令式中N是2的幂次方的数。这里假定信号总的持续时间T为单位1。

第二步：方程式（4）有2M＋1项参加求和，并不满足2的幂次方的要求。为了利用FFT，就必须把求和项数变换为2的幂次方。为此，就需要对式（4）进行补零。

第三步：方程式（4）的求和是从－M到M，而在式（9）中求和是从0开始的。所以必须对式（4）的求和进行重排，使其从0开始。

为了做到这一点，须满足：

但N必须是2的幂次方。

求和式（4）可以分解为两项：

对第二个求和式作如下变量代换：

那么，式（12）可以改写为

由于m和m′都为子变量，可以用n来替代，则

在该方程中，要注意到两个问题：

1）从n＝M＋1到N－M－1添加了很多零点，右边的总项数为N，且为2的幂次方；

2）必须对R（n）进行重新排列。对于n＝0到M，R（n）保持不变；对于第二个求和项，就需要把R（n）的变量从式（12）的n＝－M到－1范围变为n＝N－M到N－1。这样上述方程可以写为

比较式（16）与式（9），可以发现他们是完全一样的，因此可以用FFT来完成功率谱计算。

5 计算实例

设R（n）的变量从－4～4（M＝4）取值，如图2（a）所示，一共有9项。根据式（11），N的选取范围为9～18，且要求N必须是2的幂次方，所以取N＝16，重新排列后的R（n）如图2（b）所示，经过这样的重新排列后，就可以用FFT来计算P（k）。在上面的讨论中把式（4）直接变换为式（9）的形式，所以保持了正确的相位关系。对于P（k）所得到的结果通常是一个复数值，如果需要计算功率谱，就可以取绝对值。

图2 重新排列后的R（n）

如果对R（n）稍作不同的排列，也可以得到同样的结果。这一方法就是把整个R（n）从n＝－M到M移到n＝0～2M，并在数据串的尾部补零，如图2（c）所示。这样排列的计算结果与图2（b）计算结果是不一样的，有一个相位差。但是P（k）的绝对值是一样的，绝对值对功率谱估计非常重要。这种对R（n）进行移位的方法实现起来要容易一些。

至于FFT计算后实际频率轴的分量的确定可以参考文献［10，13］给出的方法。

6 结语

由于基于自相关函数的谱估计的结构与离散傅立叶变换（DFT）的计算过程非常类似，差别仅在于每次的计算结果所对应的频谱分量是任意的。因此本文考虑利用FFT算法来实现功率谱估计的快速计算。文章首先对频谱分量所对应的参数进行严格限制，使其与DFT完全一致，然后给出了如何借助FFT来实现信号的功率谱估计给出了详细过程，最后给出了计算实例，论证了该算法的有效性，从而节省了功率谱计算的时间。

［1］林茂庸，柯有安.雷达信号理论［M］.北京：国防工业出版社，1984：1～2

［2］M.I.Skolnik.雷达系统导论［M］.第三版.左群声，徐国良，马琳，等译.北京：电子工业出版社，2006：35～91

［3］M.A.Richards.雷达信号处理基础［M］.邢孟道，王彤，李真芳，等译.北京：电子工业出版社，2008：40～71

［4］张明友，吕明.信号检测与估计［M］.北京：电子工业出版社，2005：105～111

［5］丁玉美，高西全.数字信号处理［M］.西安：西安电子科技大学出版社，2000：97～101

［6］J.W.Tukey，T.W.Coody.A Fast Computation Algorithm of Fourier Series by Machine［J］.Math.Computation，1965，19

［7］徐萃微，孙绳武.计算方法引论［M］.北京：高等教育出版社，2008：66～75

［8］薛定宇，陈阳泉.高等应用数学问题的 Matlab求解［M］.北京：清华大学出版社，2004：146～151

［9］陆大纟金.随机过程及其应用［M］.北京：清华大学出版社，2007：335～337

［10］Cleve B.Molter.Matlab数值计算［M］.喻文建，译.北京：机械工业出版社，2007：201～209

［11］李仕专，李维涛，姜全贤，等.一种基于并行计算的快速FFT IP核设计［J］.计算机与数字工程，2010，38（4）

［12］张骥.一种基于FFT计算离散小波变换的方法［J］.计算机与数字工程，2009，37（10）

［13］James B.Y.Tsui.宽带数字接收机［M］.杨小牛，陆安南，译.北京：电子工业出版社，2002：77～103