分數階傅里葉變換的數值計算方法研究＊

（1.海軍航空工程學院電子信息工程系 煙臺 264001）（2.中國人民解放軍91960部隊 汕頭 515074）

1 引言

對信號進行分數階傅里葉變換（FrFT）可以解釋為將信號坐標軸在時頻平面繞原點逆時針旋轉［1］。隨著階數a從0連續變化到1，表現出信號從時域到頻域變化的全部特征，因此，可以看做是一種統一的時頻變換。作為一維線性變換，可以有效地避免信號的交叉項干擾，具有較好的時頻聚集性。因此，特別適合于線性調頻信號（Chirp）的檢測與參數估計［2］、時頻濾波［3］、目標識別［4］等用途。

2 分數階傅里葉變換定義

分數傅里葉變換算子Fa通過實變量a將函數x變換為Xa＝Fa（x），定義可以表述為整體積分變換：

其中，a為FrFT 的階數，u為分數階域（FrFD），旋轉角α＝a，則FrFT 變換核Ka（u，t）為

式中

Xa（u）的逆變換為

由上式可知，信號x（t）可以分解為一組系數為Xa（u） 的正交Chirp基K-a（u，t）的線性組合。

3 分數階傅里葉變換數值計算

3.1 離散分數階傅里葉變換定義

仿照離散傅里葉變換（DFT）的定義，x＝［x（0），…，x（N-1）］T的離散分數階傅里葉變換（DFrFT）可以理解為x經過DFrFT算子Fa作用的結果［5］，即Xa＝Fax：

式中Fa（k，n）為DFrFT 核矩陣，Fa＝ΕΛaΕT，F＝ΕΛΕT為DFT 矩陣的特征值分解。

3.2 主要的數值計算方法

按照Fa（k，n）核矩陣的求法差異，我們可以將近年來主要的一些分數階傅里葉變換的數值算法分為以下三類［6］：

1）特征分解型

特征分解算法通過求DFT 的核矩陣F（k，n）的特征值和特征向量構造DFT 核矩陣的分數冪，以此作為Fa（k，n）來計算DFrFT。Candan等在文獻［7］中對DFrFT 進行了詳細描述。S.C.Pei等在文獻［8］提出了一種基于正交投影的特征分解型算法。此類方法的優點是保持了連續Fr-FT 的許多性質，但該法計算復雜度高，不能用來實時處理。

2）線性組合型

線性組合方法是對特定階數的一些DFrFT 的線性組合來表示任意階的DFrFT。文獻［9］利用Taylor級數展開以及Caylay-Hamilton定理，將其表示為階數為0，1，2，3的線性組合。雖然該方法計算量最小，但計算精度太低。文獻［10］提出了另外一種DFrFT 線性組合型算法，將上述特殊階數0，1，2，3變為從0階開始間隔N／4取得N個階數，其中N為樣本個數。利用FrFT 的旋轉相加性，可將一次變換的結果輸入到后一次變換，因此算法可以采用串行的方法實現，可用超大規模集成電路（VLSI）實現。但缺點是，必須知道至少1個特定階數變換的結果，同時這些階數還與輸入信號的樣本個數N有關。

3）離散采樣型

離散采樣型算法實質上是，對連續FrFT 的輸入輸出信號進行采樣獲得DFrFT 核矩陣。由Ozaktas等學者提出的計算量與FFT 相當的快速近似分數階傅里葉變換算法（FAFrFT），其本質是實現信號的離散取樣值與連續FrFT的離散取樣值之間的映射，計算復雜度僅為O（NlogN）［11］。該算法引起了國內外學者的廣泛關注，使FrFT 從理論逐步走向工程實踐，從而在信號處理領域得到廣泛應用。具體步驟是：

代入式（1），利用 （cotα-cscα） ＝-tan（α／ 2）化簡得到：

令Tr（x）＝exp（jπrx2），這里r可以取兩個值v＝-tan（α／2），c＝cscα，定義算子g＝Tv·x，算子h＝Tc＊g。式（7）可以進一步簡化為

圖1 FrFT 的1相實現原理圖

這類算法將FrFT 分解為三個步驟來實現，即步驟一信號先與Chirp信號乘積，步驟二與Chirp信號卷積運算，步驟三與Chirp 信號乘積。將此類方法實現原理理解為1相實現，如圖1所示。

這類算法的機理決定了運算前，假設信號在時域和頻域上是緊支撐的，然后對時頻域進行量綱歸一化，時域和頻域都被限定在區間［-Δ／2，Δ／2］，寬度為Δ，采樣點數為N＝Δ2，采樣間隔［11］。而且每個步驟必須考慮采樣間隔，以符合香農采樣定理?？紤]到信號與Chirp信號乘積、卷積運算的Wigner分布區域，需要將支撐區域限定在區間［-Δ，Δ］。為了恢復原來的采樣樣本值，我們以采樣間隔1／2Δ獲取樣本值，從而得到2N個樣本值。并得到離散形式如下：

這就是對連續FrFT 的近似值，直接對其計算得到的結果精度低。于是，學者對這類算法進行了改進。目前主要有兩種實現途徑：一種是由A.M.Kutay實現，另一種是由J.O’Neill實現。前者只允許采樣點數N為偶數，且a的取值限定在［-2，2］；后者要求采樣點數N為奇數。文獻［5］對兩者進行了比較，結果顯示后者在實際運算中能得到更精確的結果，并對后者進行了改進，消除了限制條件，提出了“最好”的FAFrFT 算法，可實現任意的階數a和采樣點數N的計算。

文獻［5］將Ozaktas算法［11］與Pei算法［8］進行了比較，前者是一種實現近似計算連續FrFT 更直接和快速的方法，而后者僅當采樣點數N比較大時，才可以用于近似計算連續FrFT。

3.3 二相數值算法

Ozaktas等學者提出的FAFrFT 算法，得到了廣泛應用與推廣。但是，該算法在計算的時候需要上采樣，在計算結束時再下采樣，使得計算前后的樣本數N一致。從而導致實際計算過程中，有些數據沒有必要計算，延長了算法運算時間。為避免此情況，文獻［12］提出了一種高效的多相算法，文獻［13］提出了一種簡化的2相算法。

2相算法在計算過程中，將信號分段成偶數樣本和奇數樣本部分，分別對這兩個部分進行計算，避免沒有必要的計算，縮短了算法運算時間。假設信號x（t） 有N個等間隔的樣本點數，k的取值可以是整數或者半整數。在計算時，N是偶數或者奇數導致離散樣本中的索引值取得整數值或者半整數值。不妨設定將值存放在偶數索引項部分，稱之為偶數樣本。在這些偶數樣本之間定義奇數樣本，而這些值是需要通過sinc插值得到的，共有N-1個樣本。

偶數樣本：

奇數樣本：

上式含有卷積運算，可以等價于兩個FFT 運算乘積的逆FFT 運算（IFFT）。

依照此思路，按照FrFT 的1相實現的三個步驟，可以得 到FrFT 的2相實現［13］，如圖2所示。其中As＝Aα／

圖2 FrFT 的2相實現原理圖

4 數值算法性能分析

通過計算機仿真對文獻［5］中的算法1與文獻［13］算法2進行比較。仿真平臺：cpu Pentium（R）Dual-Core 2GHz，cache1MB，系統windows xp sp3，MATLAB R2009a。設定信號為Chirp信號，解析式exp（jπkt2），離散樣本個數N＝2n，n＝6，…15，調頻率k＝1Hz／s，采樣頻率fs＝1000Hz／s。FrFT階數a等間隔從區間［0，2］取100個點。誤差系數如圖3所示，算法運算時間為Matlab的cputime如圖4所示。

圖3 a階誤差系數圖

圖4 兩種算法運算時間圖

定義誤差系數

式中abs（·）為取模計算。經計算，rmax＝0.0059。

綜合仿真結果可以看出，算法2與算法1的精度相當，而計算速度有較大的提高。

5 結語

FrFT 是對傳統FT 的推廣，已經廣泛應用于雷達、通信、聲納、水印等許多實際工程領域。由于DFrFT 的計算比DFT 的計算要復雜的多，盡管國內外學者提出了很多種FrFT 的快速數值計算方法，但兼顧FrFT 的各種性質和實際計算復雜度，目前還沒有一種FrFT 的快速計算方法滿足以上所有要求。該文首先給出了分數階傅里葉變換的整體積分定義、離散分數階傅里葉變換定義，對現有的主要幾種數值算法進行了研究和比較，最后基于Matlab分析了一種二相FrFT 快速算法，經比較其計算速度更快、有更好的實用價值，可以進一步貼近工程領域實時性需求。

［1］Almeida L B.The fractional Fourier transform and time-frequency representations［J］.IEEE Trans on SP，1994，42（11）：3084-3091.

［2］齊林，陶然，周思永，等.基于分數階fourier變換的多分量LFM信號的檢測和參數估計［J］.中國科學，2003，33（8）：749-759.

［3］鄧兵，陶然，齊林，等.分數階Fourier變換與時頻濾波［J］.系統工程與電子技術，2004，26（10）：1357-1359.

［4］謝德光，張賢達.基于分數階Fourier變換的雷達目標識別［J］.清華大學學報，2010，50（4）：485-488.

［5］Bultheel A and Sulbaran H E M.Computation of the Fractional Fourier Trasform［J］.Appl.Comput.Harmonic Anal，2004，16（3）：182-202.

［6］陶然，鄧兵，王越.分數階Fourier變換在信號處理領域的研究進展［J］.中國科學，2006，36（2）：113-116.

［7］Candan C.The discrete Fractional Fourier Transform［D］.MS Thesis，Bilkent University，Ankara，1998.

［8］Pei S C，Yeh M H，Tseng C C.Digital fractional Fourier transform based on orthogonal projections［J］.IEEE Trans on SP，1999，47（5）：1335-1348.

［9］Santhanam B，McClellan J H.The discrEte rotational Fourier transform［J］.IEEE Trans on SP，1996，44（4）：994-998.

［10］Yeh M H，Pei S C.A method for the discrete fractional Fourier transform computation［J］.IEEE Trans on SP，2003，51（3）：889-891.

［11］Ozaktas H M，Ankan Orhan.Digital Computation of the Fractional Fourier Transform［J］.IEEE Trans on SP，1996，44（9）：2141-2150.

［12］TaoR，Liang G，Zhao X H.An Efficient FPGA-based Implementation of Fractional Fourier Transform Algorithm［J］.Sign Process Syst，2010，60（1）：47-58.

［13］Bultheel Adhermar.A two-phase implementation of the fractional Fourier transform［R］.Report TW 588，Dept.Computer Science，K.U.Leuven，March，2011.

［14］張明，柳超，張志剛.鞭天線線性陣列方向圖的數值計算［J］.計算機與數字工程，2009（7）.