基于FFT與自相關(guān)函數(shù)的快速功率譜估計方法＊

李春林 伍 勇

（92755部隊67分隊1） 臨高 571820）（92098部隊電子對抗科2） 陵水 572425）

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

根據(jù)以上式子，可以得到有限數(shù)據(jù)的功率譜為

1 引言

在無線電技術(shù)（比如通信、雷達(dá)、電子戰(zhàn)、遙控遙測等）領(lǐng)域，無線電系統(tǒng)所處理的主要對象是電磁波［1～2］。無線電系統(tǒng)中的接收天線首先把電磁波轉(zhuǎn)化為電信號，然后饋入到接收系統(tǒng)進行必要的分析和處理。在20世紀(jì)70年代以前，對電磁信號的處理主要以模擬處理方法為主［3～5］，包含有濾波、放大、混頻、檢波等等各種模擬處理環(huán)節(jié)。模擬處理方法不僅缺乏靈活性、可擴展性，而且每一環(huán)節(jié)都會不同程度地引入各種非線性失真，導(dǎo)致接收系統(tǒng)性能下降，功能降低。

自從1965年J.W.Tukey和T.W.Coody在《計算數(shù)學(xué)》雜志上發(fā)表了著名的《機器計算傅立葉級數(shù)的一種算法》［6］論文并幾經(jīng)人們改進之后，很快形成了一套高效的算法FFT。80年代以后，隨著微電子技術(shù)的迅猛發(fā)展，基于FFT的數(shù)字信號處理技術(shù)開始獲得廣泛應(yīng)用，并逐步顯現(xiàn)出模擬處理方式所無法達(dá)到的優(yōu)越性［7～8］。

本文將引入FFT用于基于自相關(guān)函數(shù)的功率譜估計中，加快估計的計算速度。在信號功率譜估計中有許多的高分辨率譜估計算法，之所以討論FFT在關(guān)于譜估計中的應(yīng)用，是因為FFT已經(jīng)廣泛應(yīng)用到數(shù)字接收機的設(shè)計中。文中首先簡要介紹了有偏和無偏自相關(guān)函數(shù)的定義，然后給出了基于自相關(guān)函數(shù)的譜估計算法，最后就如何借助FFT來實現(xiàn)信號的功率譜估計給出了詳細(xì)過程，并給出了計算實例。

2 自相關(guān)函數(shù)的定義及其計算過程

假設(shè)有N點輸入數(shù)據(jù)x（n），n＝0，…，N－1，其自相關(guān)定義為［1，9］：

式中m稱之為自相關(guān)的延遲變量。嚴(yán)格意義下，上述自相關(guān)函數(shù)應(yīng)該稱為取樣自相關(guān)，它逼近于E［x（n）x（n＋m）］，式中E［·］表示期望值。m值既可以正也可以負(fù)。如果變量m為負(fù)數(shù)，那么其自相關(guān)函數(shù)與正m的自相關(guān)函數(shù)的關(guān)系為R（－m）＝R（m）＊，其中R（k）是復(fù)數(shù)，如果R（k）是實數(shù)，則R（－m）＝R（m）。如果延遲變量是m，那么把0～N－1的輸入數(shù)據(jù)分成長度相同的兩組，一組從0～N－M－1，另一組從m～N－1。這兩組中數(shù)據(jù)項相乘的情況如0所示。所有乘積項的和等于自相關(guān)函數(shù)R（m）的N倍。

圖1 R（m）的計算過程描述

自相關(guān)可以視為兩組數(shù)據(jù)之間相似性的度量。如果兩組數(shù)據(jù)很相似，那么其自相關(guān)函數(shù)就大，相反就小。當(dāng)m＝0時，兩組數(shù)據(jù)是相同的，所以在所有自相關(guān)函數(shù)值中R（0）最大。當(dāng)m值比較大時，就會只有很少幾項參加求和運算，但式（1）中的分母N是一個固定常數(shù)，所以求和式被N除之后，通常使R（m）的幅度變得很小，而理論上該值可能是比較大的。因此，式（1）所定義的自相關(guān)函數(shù)稱為自相關(guān)函數(shù)的有偏形式。

而無偏自相關(guān)函數(shù)的定義為

在這一定義中，當(dāng)m變大時，分母變小，而且等于求和的項數(shù)。由于這種形式會產(chǎn)生負(fù)的功率譜，所以很少用在譜估計中。如果R（m）比較小，那么對所計算的功率譜的影響就會小一些。由于R（m）只有很少幾個點，所以預(yù)計R（m）值會比較小，這樣對功率譜的影響也就小了。

3 基于自相關(guān)函數(shù)的功率譜估計的結(jié)構(gòu)

可以采用兩種方法來估計功率譜。第一種方法就是計算輸入信號的FFT，然后對其結(jié)果求平方，這一方法叫周期圖法。這一結(jié)果可以直接從FFT運算得到。第二種方法，功率譜是從輸入數(shù)據(jù)的自相關(guān)函數(shù)得到的，可以表示為［9］

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

根據(jù)以上式子，可以得到有限數(shù)據(jù)的功率譜為

式中，m取值區(qū)間為從－M～M，k是頻率分量，ts是采樣間隔。該方法通常稱為Blackman－Tukey方法［11～12］。求和項的總數(shù)為2M＋1，因為求和式中包含了m＝0的項。

為了獲得正的功率譜，在上述等式中通常采用有偏自相關(guān)R（m），可以把有偏自相關(guān)R（m）作為一個加窗函數(shù)。當(dāng)m比較大時，R（m）通常就小。有時把求和項大約限制在m＝－N／10～N／10之間，推薦的最大值為－N／5～N／5之間。另外還可以對有偏自相關(guān)加一個窗函數(shù)，以進一步減小旁瓣。

加窗后的功率譜可以寫為

式中，W（m）為窗函數(shù)。該窗函數(shù)必須是對稱的，保證P（k）為偶函數(shù)。

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

R（m）的期望值為

這是由于R（m）與n無關(guān)。求和式等于N－｜m｜，所以式（7）可以寫為

式中，（N－｜m｜）／N表示一個三角窗。由此可見，該方法受固有加窗效應(yīng)的內(nèi)在限制，而m＞N時，R（m）的值假定為0。

4 基于FFT的功率譜估計步驟

用FFT來計算式（4）的結(jié)果。式（4）與DFT非常類似，但k值是任意的。比如，從以上方程可以計算出任意給定k值所對應(yīng)的結(jié)果。如果對k嚴(yán)格加以限制，使其與DFT完全一致，那么就可以采用FFT算法來計算，從而節(jié)省計算時間。為此，必須把式（4）變?yōu)楹线m的形式。DFT可寫為

式中，X（K）為頻域響應(yīng)，x（n）為時域取樣數(shù)據(jù)點。在上述方程中，n和k是離散的，N通常取為2的冪次方。

對式（4）進行適當(dāng)?shù)淖冃危璨扇∫韵聨讉€步驟：

第一步：改變標(biāo)記符t。通過比較式（4）和式（9）可以看出，令式中N是2的冪次方的數(shù)。這里假定信號總的持續(xù)時間T為單位1。

第二步：方程式（4）有2M＋1項參加求和，并不滿足2的冪次方的要求。為了利用FFT，就必須把求和項數(shù)變換為2的冪次方。為此，就需要對式（4）進行補零。

第三步：方程式（4）的求和是從－M到M，而在式（9）中求和是從0開始的。所以必須對式（4）的求和進行重排，使其從0開始。

為了做到這一點，須滿足：

但N必須是2的冪次方。

求和式（4）可以分解為兩項：

對第二個求和式作如下變量代換：

那么，式（12）可以改寫為

由于m和m′都為子變量，可以用n來替代，則

在該方程中，要注意到兩個問題：

1）從n＝M＋1到N－M－1添加了很多零點，右邊的總項數(shù)為N，且為2的冪次方；

2）必須對R（n）進行重新排列。對于n＝0到M，R（n）保持不變；對于第二個求和項，就需要把R（n）的變量從式（12）的n＝－M到－1范圍變?yōu)閚＝N－M到N－1。這樣上述方程可以寫為

比較式（16）與式（9），可以發(fā)現(xiàn)他們是完全一樣的，因此可以用FFT來完成功率譜計算。

5 計算實例

設(shè)R（n）的變量從－4～4（M＝4）取值，如圖2（a）所示，一共有9項。根據(jù)式（11），N的選取范圍為9～18，且要求N必須是2的冪次方，所以取N＝16，重新排列后的R（n）如圖2（b）所示，經(jīng)過這樣的重新排列后，就可以用FFT來計算P（k）。在上面的討論中把式（4）直接變換為式（9）的形式，所以保持了正確的相位關(guān)系。對于P（k）所得到的結(jié)果通常是一個復(fù)數(shù)值，如果需要計算功率譜，就可以取絕對值。

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

如果對R（n）稍作不同的排列，也可以得到同樣的結(jié)果。這一方法就是把整個R（n）從n＝－M到M移到n＝0～2M，并在數(shù)據(jù)串的尾部補零，如圖2（c）所示。這樣排列的計算結(jié)果與圖2（b）計算結(jié)果是不一樣的，有一個相位差。但是P（k）的絕對值是一樣的，絕對值對功率譜估計非常重要。這種對R（n）進行移位的方法實現(xiàn)起來要容易一些。

至于FFT計算后實際頻率軸的分量的確定可以參考文獻(xiàn)［10，13］給出的方法。

6 結(jié)語

由于基于自相關(guān)函數(shù)的譜估計的結(jié)構(gòu)與離散傅立葉變換（DFT）的計算過程非常類似，差別僅在于每次的計算結(jié)果所對應(yīng)的頻譜分量是任意的。因此本文考慮利用FFT算法來實現(xiàn)功率譜估計的快速計算。文章首先對頻譜分量所對應(yīng)的參數(shù)進行嚴(yán)格限制，使其與DFT完全一致，然后給出了如何借助FFT來實現(xiàn)信號的功率譜估計給出了詳細(xì)過程，最后給出了計算實例，論證了該算法的有效性，從而節(jié)省了功率譜計算的時間。

［1］林茂庸，柯有安.雷達(dá)信號理論［M］.北京：國防工業(yè)出版社，1984：1～2

［2］M.I.Skolnik.雷達(dá)系統(tǒng)導(dǎo)論［M］.第三版.左群聲，徐國良，馬琳，等譯.北京：電子工業(yè)出版社，2006：35～91

［3］M.A.Richards.雷達(dá)信號處理基礎(chǔ)［M］.邢孟道，王彤，李真芳，等譯.北京：電子工業(yè)出版社，2008：40～71

［4］張明友，呂明.信號檢測與估計［M］.北京：電子工業(yè)出版社，2005：105～111

［5］丁玉美，高西全.數(shù)字信號處理［M］.西安：西安電子科技大學(xué)出版社，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］薛定宇，陳陽泉.高等應(yīng)用數(shù)學(xué)問題的 Matlab求解［M］.北京：清華大學(xué)出版社，2004：146～151

［9］陸大纟金.隨機過程及其應(yīng)用［M］.北京：清華大學(xué)出版社，2007：335～337

［10］Cleve B.Molter.Matlab數(shù)值計算［M］.喻文建，譯.北京：機械工業(yè)出版社，2007：201～209

［11］李仕專，李維濤，姜全賢，等.一種基于并行計算的快速FFT IP核設(shè)計［J］.計算機與數(shù)字工程，2010，38（4）

［12］張驥.一種基于FFT計算離散小波變換的方法［J］.計算機與數(shù)字工程，2009，37（10）

［13］James B.Y.Tsui.寬帶數(shù)字接收機［M］.楊小牛，陸安南，譯.北京：電子工業(yè)出版社，2002：77～103