基于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