（1.海軍航空工程學(xué)院電子信息工程系 煙臺(tái) 264001）（2.中國(guó)人民解放軍91960部隊(duì) 汕頭 515074）

1 引言

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

2 分?jǐn)?shù)階傅里葉變換定義

分?jǐn)?shù)傅里葉變換算子Fa通過(guò)實(shí)變量a將函數(shù)x變換為Xa＝Fa（x），定義可以表述為整體積分變換：

其中，a為FrFT 的階數(shù)，u為分?jǐn)?shù)階域（FrFD），旋轉(zhuǎn)角α＝a，則FrFT 變換核Ka（u，t）為

式中

Xa（u）的逆變換為

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

3 分?jǐn)?shù)階傅里葉變換數(shù)值計(jì)算

3.1 離散分?jǐn)?shù)階傅里葉變換定義

仿照離散傅里葉變換（DFT）的定義，x＝［x（0），…，x（N-1）］T的離散分?jǐn)?shù)階傅里葉變換（DFrFT）可以理解為x經(jīng)過(guò)DFrFT算子Fa作用的結(jié)果［5］，即Xa＝Fax：

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

3.2 主要的數(shù)值計(jì)算方法

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

1）特征分解型

特征分解算法通過(guò)求DFT 的核矩陣F（k，n）的特征值和特征向量構(gòu)造DFT 核矩陣的分?jǐn)?shù)冪，以此作為Fa（k，n）來(lái)計(jì)算DFrFT。Candan等在文獻(xiàn)［7］中對(duì)DFrFT 進(jìn)行了詳細(xì)描述。S.C.Pei等在文獻(xiàn)［8］提出了一種基于正交投影的特征分解型算法。此類(lèi)方法的優(yōu)點(diǎn)是保持了連續(xù)Fr-FT 的許多性質(zhì)，但該法計(jì)算復(fù)雜度高，不能用來(lái)實(shí)時(shí)處理。

2）線性組合型

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

3）離散采樣型

離散采樣型算法實(shí)質(zhì)上是，對(duì)連續(xù)FrFT 的輸入輸出信號(hào)進(jìn)行采樣獲得DFrFT 核矩陣。由Ozaktas等學(xué)者提出的計(jì)算量與FFT 相當(dāng)?shù)目焖俳品謹(jǐn)?shù)階傅里葉變換算法（FAFrFT），其本質(zhì)是實(shí)現(xiàn)信號(hào)的離散取樣值與連續(xù)FrFT的離散取樣值之間的映射，計(jì)算復(fù)雜度僅為O（NlogN）［11］。該算法引起了國(guó)內(nèi)外學(xué)者的廣泛關(guān)注，使FrFT 從理論逐步走向工程實(shí)踐，從而在信號(hào)處理領(lǐng)域得到廣泛應(yīng)用。具體步驟是：

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

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

圖1 FrFT 的1相實(shí)現(xiàn)原理圖

這類(lèi)算法將FrFT 分解為三個(gè)步驟來(lái)實(shí)現(xiàn)，即步驟一信號(hào)先與Chirp信號(hào)乘積，步驟二與Chirp信號(hào)卷積運(yùn)算，步驟三與Chirp 信號(hào)乘積。將此類(lèi)方法實(shí)現(xiàn)原理理解為1相實(shí)現(xiàn)，如圖1所示。

這類(lèi)算法的機(jī)理決定了運(yùn)算前，假設(shè)信號(hào)在時(shí)域和頻域上是緊支撐的，然后對(duì)時(shí)頻域進(jìn)行量綱歸一化，時(shí)域和頻域都被限定在區(qū)間［-Δ／2，Δ／2］，寬度為Δ，采樣點(diǎn)數(shù)為N＝Δ2，采樣間隔［11］。而且每個(gè)步驟必須考慮采樣間隔，以符合香農(nóng)采樣定理。考慮到信號(hào)與Chirp信號(hào)乘積、卷積運(yùn)算的Wigner分布區(qū)域，需要將支撐區(qū)域限定在區(qū)間［-Δ，Δ］。為了恢復(fù)原來(lái)的采樣樣本值，我們以采樣間隔1／2Δ獲取樣本值，從而得到2N個(gè)樣本值。并得到離散形式如下：

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

文獻(xiàn)［5］將Ozaktas算法［11］與Pei算法［8］進(jìn)行了比較，前者是一種實(shí)現(xiàn)近似計(jì)算連續(xù)FrFT 更直接和快速的方法，而后者僅當(dāng)采樣點(diǎn)數(shù)N比較大時(shí)，才可以用于近似計(jì)算連續(xù)FrFT。

3.3 二相數(shù)值算法

Ozaktas等學(xué)者提出的FAFrFT 算法，得到了廣泛應(yīng)用與推廣。但是，該算法在計(jì)算的時(shí)候需要上采樣，在計(jì)算結(jié)束時(shí)再下采樣，使得計(jì)算前后的樣本數(shù)N一致。從而導(dǎo)致實(shí)際計(jì)算過(guò)程中，有些數(shù)據(jù)沒(méi)有必要計(jì)算，延長(zhǎng)了算法運(yùn)算時(shí)間。為避免此情況，文獻(xiàn)［12］提出了一種高效的多相算法，文獻(xiàn)［13］提出了一種簡(jiǎn)化的2相算法。

2相算法在計(jì)算過(guò)程中，將信號(hào)分段成偶數(shù)樣本和奇數(shù)樣本部分，分別對(duì)這兩個(gè)部分進(jìn)行計(jì)算，避免沒(méi)有必要的計(jì)算，縮短了算法運(yùn)算時(shí)間。假設(shè)信號(hào)x（t） 有N個(gè)等間隔的樣本點(diǎn)數(shù)，k的取值可以是整數(shù)或者半整數(shù)。在計(jì)算時(shí)，N是偶數(shù)或者奇數(shù)導(dǎo)致離散樣本中的索引值取得整數(shù)值或者半整數(shù)值。不妨設(shè)定將值存放在偶數(shù)索引項(xiàng)部分，稱(chēng)之為偶數(shù)樣本。在這些偶數(shù)樣本之間定義奇數(shù)樣本，而這些值是需要通過(guò)sinc插值得到的，共有N-1個(gè)樣本。

偶數(shù)樣本：

奇數(shù)樣本：

上式含有卷積運(yùn)算，可以等價(jià)于兩個(gè)FFT 運(yùn)算乘積的逆FFT 運(yùn)算（IFFT）。

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

圖2 FrFT 的2相實(shí)現(xiàn)原理圖

4 數(shù)值算法性能分析

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

圖3 a階誤差系數(shù)圖

圖4 兩種算法運(yùn)算時(shí)間圖

定義誤差系數(shù)

式中abs（·）為取模計(jì)算。經(jīng)計(jì)算，rmax＝0.0059。

綜合仿真結(jié)果可以看出，算法2與算法1的精度相當(dāng)，而計(jì)算速度有較大的提高。

5 結(jié)語(yǔ)

FrFT 是對(duì)傳統(tǒng)FT 的推廣，已經(jīng)廣泛應(yīng)用于雷達(dá)、通信、聲納、水印等許多實(shí)際工程領(lǐng)域。由于DFrFT 的計(jì)算比DFT 的計(jì)算要復(fù)雜的多，盡管?chē)?guó)內(nèi)外學(xué)者提出了很多種FrFT 的快速數(shù)值計(jì)算方法，但兼顧FrFT 的各種性質(zhì)和實(shí)際計(jì)算復(fù)雜度，目前還沒(méi)有一種FrFT 的快速計(jì)算方法滿足以上所有要求。該文首先給出了分?jǐn)?shù)階傅里葉變換的整體積分定義、離散分?jǐn)?shù)階傅里葉變換定義，對(duì)現(xiàn)有的主要幾種數(shù)值算法進(jìn)行了研究和比較，最后基于Matlab分析了一種二相FrFT 快速算法，經(jīng)比較其計(jì)算速度更快、有更好的實(shí)用價(jià)值，可以進(jìn)一步貼近工程領(lǐng)域?qū)崟r(shí)性需求。

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

［2］齊林，陶然，周思永，等.基于分?jǐn)?shù)階fourier變換的多分量LFM信號(hào)的檢測(cè)和參數(shù)估計(jì)［J］.中國(guó)科學(xué)，2003，33（8）：749-759.

［3］鄧兵，陶然，齊林，等.分?jǐn)?shù)階Fourier變換與時(shí)頻濾波［J］.系統(tǒng)工程與電子技術(shù)，2004，26（10）：1357-1359.

［4］謝德光，張賢達(dá).基于分?jǐn)?shù)階Fourier變換的雷達(dá)目標(biāo)識(shí)別［J］.清華大學(xué)學(xué)報(bào)，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］陶然，鄧兵，王越.分?jǐn)?shù)階Fourier變換在信號(hào)處理領(lǐng)域的研究進(jìn)展［J］.中國(guó)科學(xué)，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］張明，柳超，張志剛.鞭天線線性陣列方向圖的數(shù)值計(jì)算［J］.計(jì)算機(jī)與數(shù)字工程，2009（7）.