分?jǐn)?shù)階傅立葉變換的最優(yōu)階數(shù)論文

1、現(xiàn)代數(shù)字信號處理學(xué) 號：140808040219學(xué)生所在學(xué)院：測試與光電工程學(xué)院學(xué) 生 姓 名 ：任 課 教 師 ：李 志 農(nóng)教師所在學(xué)院：測試與光電工程學(xué)院2015年1月分?jǐn)?shù)階傅立葉變換的最優(yōu)階數(shù)摘 要：傳統(tǒng)傅立葉變換描述了信號時域或頻域的特性，而不是描述信號時頻特性，于是人們提出了一系列新的時頻分析理論和方法來處理非平穩(wěn)信號，分?jǐn)?shù)階傅立葉變換為其中的一種。本文主要介紹了它的定義、性質(zhì)，還有它的離散算法，介紹了求最優(yōu)階數(shù)的方法，主要是峰值搜索算法。最后進(jìn)行仿真驗證，利用MATLAB對一個已知的chirp信號求解最優(yōu)階數(shù)。關(guān)鍵字：傅立葉變換；分?jǐn)?shù)階傅立葉變換；峰值搜索算法；MATLAB；最優(yōu)階

2、數(shù)1 引言傳統(tǒng)的傅立葉變換在所有的信號處理工具中是應(yīng)用最廣泛、研究最成熟的數(shù)學(xué)工具，作為一種線性算子，傳統(tǒng)傅立葉變換可視為在時頻平面上，信號從時間軸逆時針旋轉(zhuǎn)到頻率軸，而 FRFT 作為 FT 的廣義形式可理解為對信號旋轉(zhuǎn)任意角度的線性算子，從而可以得到信號的任意階次或者任意分?jǐn)?shù)階傅立葉域上的 FRFT 表示，并且在保留了傳統(tǒng)的 FT 所有性質(zhì)和優(yōu)點的基礎(chǔ)之上又增添了新的優(yōu)勢。2 FRFT 的定義及其性質(zhì)1.1 FRFT 的定義如圖1.所示如果把信號的分?jǐn)?shù)階傅立葉變換看作是從時間-頻率平面旋轉(zhuǎn)的話，那么傅立葉變換就相當(dāng)于在時頻平面中逆時針旋轉(zhuǎn)了角度，從時間域變換到頻率域。令，是一個分?jǐn)?shù)，那么就

3、可以在時頻平面內(nèi)以任意角度的旋轉(zhuǎn)定義線性算子，記作，我們就可以把傅立葉變換推廣到任意角度即分?jǐn)?shù)階傅立葉變換。圖1.平面旋轉(zhuǎn)角度變成平面分?jǐn)?shù)階傅立葉變換的定義為： (2-1)其中是核函數(shù)， (2-2)這里，p=1或-1時退化成為常規(guī)的傅立葉變換和逆變換。分?jǐn)?shù)階傅立葉變換和經(jīng)典傅立葉變換具有以下的關(guān)系： 1.分?jǐn)?shù)階傅立葉變換是線性算子 2.周期性 (2-3) (2-4)1.2 分?jǐn)?shù)傅立葉變換的性質(zhì)（1） 線性分?jǐn)?shù)階傅立葉變換為線性變換，滿足疊加原理：若和分別是原函數(shù)和的階分?jǐn)?shù)傅立葉變換，則有 (2-5)（2） 旋轉(zhuǎn)相加行 (2-6)（3） 逆 (2-7)（4） 酋性 (2-8)（5） 交換性 (2

4、-9)（6） 結(jié)合性 (2-10)（7） 周期性 (2-11)（8） 特征函數(shù) (2-12)（9） 卷積、相乘、相關(guān) 函數(shù)、在階次 p 分?jǐn)?shù)階傅立葉域的卷積記作 (2-13)函數(shù)、 在階次 p 分?jǐn)?shù)階傅立葉域的乘積記作 (2-14)函數(shù)、 在階次 p 分?jǐn)?shù)階傅立葉域的相關(guān)定義為 (2-15)（10） 時移特性 (2-16)（11） 頻移特性 (2-17)（12） 尺度特性 (2-18)其中（13） 微分特性 (2-19)（14） 積分特性 (2-20)3分?jǐn)?shù)階傅立葉變換的離散的離散算法分?jǐn)?shù)階傅立葉變換的出現(xiàn)引起了各個領(lǐng)域研究人員的重視，在工程上也有十分廣闊的應(yīng)用前景。在數(shù)字信號處理的應(yīng)用中，必

5、須采用離散形式的分?jǐn)?shù)階傅立葉變換（DFRFT），這使得離散分?jǐn)?shù)階傅立葉變換及其快速算法的研究顯得十分重要。目前，DFRFT 的離散化算法主要有四種：1. 利用來計算離散 FRFT 的核矩陣，從而利用 FFT 來計算DFRFT其中 W 是離散傅立葉變換核矩陣。這種方法實際上缺乏理論基礎(chǔ)，而且其離散 FRFT 矩陣不滿足連續(xù) FRFT 的旋轉(zhuǎn)相加性，因此不能用相同的方法計算逆 FRFT。實際計算所產(chǎn)生的誤差比較大，與連續(xù) FRFT 沒有相似的輸出結(jié)果。其計算復(fù)雜度與傳統(tǒng)傅立葉變換相同，為。2. 分解方法根據(jù) FRFT 的表達(dá)式，將 FRFT 分解為信號的卷積形式，從而利用 FFT 來計算FRFT。

6、這種方法思想比較直觀，計算出的記過與連續(xù) FRFT 的輸出比較接近。但它要經(jīng)過一次 2 倍內(nèi)插和 2 倍抽取，而且還要進(jìn)行坐標(biāo)的無量綱化，實現(xiàn)起來較為煩瑣。其計算復(fù)雜度為。3. 利用矩陣的特征值和特征向量來計算 DFRFT這種方法保持了連續(xù) FRFT 的特征值-特征函數(shù)的關(guān)系，克服了第一種方法中特征值和特征向量不匹配的缺點。采用了兩種正交映射的辦法 DFT 的 Hermite 特征向量，由此開發(fā)出兩種快速方法，即 OPA 方法和 GSA 方法，這兩種方法都有和連續(xù)FRFT 相近的輸出結(jié)果，可逆性好。其計算復(fù)雜度均為。4. 直接對 FRFT 進(jìn)行離散化來計算 DFRFT。 這種方法采用直接將連續(xù)

7、 FRFT 離散化的方法來獲得離散 FRFT 的核矩陣，避開了煩瑣的特征值和特征向量匹配問題以及矩陣的正交歸一化運算，與連續(xù) FRFT有相似的輸出結(jié)果，該算法的計算復(fù)雜度為。在目前的研究中，采用的最多的是分解方法和矩陣特征值和特征向量兩種方法，下面的內(nèi)容重點介紹這兩種方法。3.1 分解方法所謂分解方法是指根據(jù) FRFT 的表達(dá)式，將 FRFT 分解為信號的卷積形式，從而利用 FFT 來計算 FRFT。這種算法由 H.M.Ozaktas 等人提出，其計算速度幾乎與FFT 相當(dāng)，被公認(rèn)為目前計算速度最快的一種 FRFT 數(shù)字計算方法，非常適合于對信號進(jìn)行實時的 FRFT 計算。但這種快速算法的運算

8、機理決定了在進(jìn)行 FRFT 數(shù)值計算之前必須對原始信號進(jìn)行量綱歸一化處理。3.1.1 量綱歸一化原理如果一個信號在時間軸或頻率軸的一個子集取非零值，并且取非零值的條件限定在有限區(qū)間內(nèi)，則稱該信號在時間軸或頻率軸上是緊湊的。從理論上講，任何信號和它的傅立葉變換不可能是同時緊湊的。然而我們實際中需要處理的信號往往是時限的和帶限的。信號的時寬帶寬積可以用來確定信號的采樣頻率和采樣點數(shù)，用于唯一地從離散化的信號中恢復(fù)原始信號。假定原始連續(xù)信號在時間軸和頻率軸上都是緊湊的，其時域表示限定在區(qū)間，而其頻域表示限定在區(qū)間， 和 分別表示信號的時寬和帶寬。信號的時寬帶寬積為，根據(jù)不確定性定理，的值應(yīng)當(dāng)大于 1

9、。由于時域和頻域具有不同的量綱，為了 FRFT 計算處理方便，需要將時域和頻域分別轉(zhuǎn)換成無量綱的域。具體方法是引入一個具有時間量綱的尺度因子，并定義新的尺度化坐標(biāo),新的坐標(biāo)系實現(xiàn)了無量綱化。信號在新坐標(biāo)系中被限定在區(qū)間和內(nèi) 。 為 使 兩 個 區(qū) 間 的 長 度 相 等 ， 選 擇，則兩個區(qū)間長度都等于無量綱量，即兩個區(qū)間歸一化為。歸一化以后信號的 Wigner-Ville 分布限定在以原點為中心，直徑 的圓內(nèi)，如圖 3-2 所示。最后根據(jù)采樣定理對歸一化后的信號進(jìn)行采樣，采樣間隔為 ，采樣點數(shù)為.圖3.1 歸一化后的時頻支撐區(qū)域3.1.2分解算法可以把（2-1）式改寫為 (3-1)式可以具體

10、分解為以下幾個步驟：（1） 用 chirp 信號調(diào)制信號： (3-2)（2） 調(diào)制信號與另一個 chirp 信號卷積： (3-3)（3） 用 chirp 信號調(diào)制卷積后的信號： (3-4)要實現(xiàn) FRFT 的數(shù)值計算必須對以上每個分解步驟都進(jìn)行離散化處理，具體的實現(xiàn)過程如下：（1）對 與 chirp 信號的乘積進(jìn)行采樣 假定分?jǐn)?shù)階次 ，chirp 信號的調(diào)頻率 ，經(jīng) chirp 信號調(diào)制后所得的信號時寬帶寬積可以是原信號的時寬帶寬積的兩倍，所以要求 的采樣間隔為，如果原來的采樣間隔是 ，可通過插值的方法獲得樣本值，然后再 chirp 信號的離散采樣相乘，以得到的采樣值。（2）實現(xiàn) 與 chir

11、p 信號的卷積由于是帶限信號，所以 chirp 信號也可取其帶限形式。所以有： (3-5)其中而則是 chirp 信號傅立葉變換 (3-6)于是，式的離散形式為 (3-7)這一離散卷積可利用 FFT 來計算。（3）計算分?jǐn)?shù)階傅立葉變換的采樣值 由于在第一步操作時對信號作了 2 倍內(nèi)插操作，所以在最后的結(jié)果需要對再進(jìn)行 2 倍抽取，以得到離散采樣?？傊鲜龇椒◤倪B續(xù)信號的 N 個離散采樣開始，最后得到由的 N 個離散樣本值得注意的是上述方法只適用于的情況，對位于該區(qū)間外的情況，可以利用分?jǐn)?shù)階傅立葉變換的周期性將階次變到 后再進(jìn)行計算。3.1.3 特征值和特征向量方法分解方法雖然計算復(fù)雜度較低，但

12、不嚴(yán)格遵守分?jǐn)?shù)階傅立葉變換的旋轉(zhuǎn)特性，因此不能從變換后的信號通過其逆變換精確地恢復(fù)原始信號。為了克服這種方法的缺點和不足，Pei Soo-Chang 等人提出一種新的離散化方法，該方法具有與連續(xù)情況相似的變換性質(zhì)和結(jié)果，并可以通過其逆變換恢復(fù)原始信號。 這種方法就是特征值和特征向量方法，它從連續(xù)傅立葉變換的特征函數(shù)為Hermite 函數(shù)出發(fā)，通過對 Hermite 函數(shù)的離散化近似和正交投影，得到一組與Hermite 函數(shù)形狀相似的 DFT 矩陣的正交化離散 Hermite 特征向量，然后，按照連續(xù)分?jǐn)?shù)階傅立葉變換的核函數(shù)譜分解表達(dá)式，構(gòu)造離散分?jǐn)?shù)階傅立葉變換矩陣。3.1.3.1 DFT 矩陣

13、特征值和特征向量傅立葉變換的特征函數(shù)為Hermite-Gaussian 函數(shù)，其表達(dá)式為： (3-8)其中，為階Hermite多項式。DFT 矩陣的特征值及重復(fù)度如下表：DFT 矩陣 F 的特征值是，共有 1、-j、-1、j 四個值，每個特征值對應(yīng)的特征向量全體組成一個特征子空間，記為，每個特征值的重復(fù)度決定了子空間的秩。矩陣 S 可用于計算 DFT 矩陣 F 的特征向量，S 的表達(dá)式為： (3-9)可以證明矩陣 S 和 F 滿足乘法交換性，即 SF = FS。因此矩陣 S 的特征向量也是矩陣 F 的特征向量，但它們對應(yīng)不同的特征值。3.1.3.2 DFT 矩陣 Hermite 特征向量的計算

14、由于 DFT 矩陣F的特征向量不唯一，我們希望從中找到與Hermite函數(shù)相似的特征向量，這樣的向量稱DFT矩陣Hermite特征向量，為了求Hermite 特征向量，下面給出一系列的重要定理。定理1 對 DFT 矩陣的 Hermite 特征向量而言，它對應(yīng)的連續(xù)函數(shù)的擴展方差應(yīng)當(dāng)為是信號的采樣間隔，連續(xù) Hermite 函數(shù)采樣后得到： (3-10)上式也可以看做是方差為 1 的 Hermite-Gaussian 函數(shù)，以為采樣間隔進(jìn)行采樣得到的序列。定理2 若序列是由單位方差Hermite-Gaussian 函數(shù)經(jīng)過采樣到的，那么，可以證明下列近似等式成立：當(dāng) N 為偶數(shù)時 (3-11)當(dāng)

15、 N 為奇數(shù)時 (3-12)定理3 將序列按照以下方式平移得到當(dāng) N 為偶數(shù)時 (3-13)當(dāng) N 為奇數(shù)時 (3-14)則的離散傅立葉變換近似為，即當(dāng) N 足夠大時 (3-15)從定理 2 和定理 3 可以看出，Hermite 函數(shù)的采樣序列近似為 DFT 矩陣特征向量。對 Hermite 函數(shù)的采樣序列作歸一化，以記為： (3-16)通過 S 矩陣可以得到 DFT 矩陣 F 的一組實正交特征向量。因此可以將這些特征向量作為 DFT 特征子空間的基向量。然后計算向量 在 DFT 特征子空間的投影，從而得到 DFT 矩陣的 Hermite 特征向量。計算方法如下： (3-17)3.1.3.3

16、DFRFT 核矩陣和二維 DFRFT4 最優(yōu)變換階次的選擇本文選擇最優(yōu)變換階次的方法主要是峰值搜索法，目前主要搜索算法有兩種二維搜索算法、擬牛頓搜索算法。盡管擬牛頓算法比二維搜索的計算量有所減但是這兩種算法的計算復(fù)雜度還是不能滿足工程實際的實時處理要求,因此我們一維曲線擬合來代替二維搜索的峰值搜索算法,計算量大大簡化。下面主要介紹幾種算法。4.1 二維搜索算法傳統(tǒng)的二維搜索法若以變換階數(shù)P為變量進(jìn)行掃描搜索的過程中,當(dāng)參數(shù)估計的估計精度要求較高時,為了滿足精度要求,變換階數(shù)P必須選擇小的搜索步長,這樣就會成倍地增加計算的復(fù)雜度,如變換階數(shù)在區(qū)間0,2上,若以步長0.0001進(jìn)行搜索,則需要進(jìn)行

17、20000次的FRFT,計算量很大,通常采用步進(jìn)搜索算法,即先以大步長再以小步長搜索最大值點,詳細(xì)步驟為:先以0.01為步長進(jìn)行階數(shù)從0-2搜索最大值點,再以步長0.0001進(jìn)行階數(shù)從-0.01至+ 0.01搜索最大值點。這樣一共需要進(jìn)行400次FRFT運算,計算量有所減小。4.2擬牛頓搜索算法對于擬牛頓搜索算法,首先介紹一下擬牛頓迭代法的基本思想。采用牛頓法解非線性方程就是把非線性線性化的一種近似方法。把在點展為泰勒級數(shù)得, (4-1)取方程的線性部分作為非線性方程的近似方程,則有 (4-2)假設(shè)則方程的解為 (4-3)再把在點展開為泰勒級數(shù),同樣也取其線性部分作為非線性方程的近似方程,假設(shè)

18、,則有 (4-4)這樣即可得到牛頓法的一個迭代序列： (4-5) (4-6)其中是第和是第n次搜索的結(jié)果,為第n次搜索的步長系數(shù),為函數(shù)在點的尺度矩陣可以通過迭代的方法求得。牛頓迭代法中每次迭代的主要運算為一次一維搜索和函數(shù)的一階偏導(dǎo)數(shù)的計算,這一擬牛頓算法所需要的計算為分?jǐn)?shù)階傅里葉域的一次掃描搜索和一次迭代搜索,因為迭代搜索所需的搜索和函數(shù)的一階偏導(dǎo)數(shù)的計算量遠(yuǎn)遠(yuǎn)小于掃描搜索,若掃描點數(shù)為m,信號的樣本長度為N,則擬牛頓算法的計算復(fù)雜度為O(mNlgN),與其它基于二維時頻分布的算法(運算量一般為)相比,擬牛頓算法計算量較小。5 仿真驗證本文中利用二維搜索法求最優(yōu)階次，該信號是一個chirp

19、信號以相加的形式調(diào)制簡單的高斯信號的混合信號，對該信號進(jìn)行求解它的最優(yōu)階數(shù)。編寫的分?jǐn)?shù)階傅立葉變換函數(shù)見附錄A所示，主函數(shù)見附錄B所示。根據(jù)運行的結(jié)果為P=0.79為所求的最優(yōu)階數(shù)。在該主程序中還用到了階數(shù)P=0.48，P=1.35，和P=1.82和所求的最優(yōu)階數(shù)進(jìn)行對比，很明顯P=0.79是它的幅值最大。初始信號見下圖4.1所示，最優(yōu)階圖見圖4.2所示。圖4.1 初始信號圖圖4.2 最優(yōu)階圖圖4.3 P=0.48，P=1.35，P=1.82和P=0.79分?jǐn)?shù)階傅立葉變換圖參考文獻(xiàn)1 平先軍,陶然,周思永等. 一種新的分?jǐn)?shù)階傅立葉變換快速算法J. 電子學(xué)報,2001,29(3):406-408

20、2 張立浩. 分?jǐn)?shù)階傅立葉變換及其在通信中的應(yīng)用D. 碩士學(xué)位論文,燕山大學(xué),2012.3 張兆祥. 基于分?jǐn)?shù)階傅立葉變換的數(shù)字水印算法研究D. 碩士學(xué)位論文,華北電力大學(xué),2007.4 陳恩慶. 基于分?jǐn)?shù)階fourier變換的ofdm系統(tǒng)關(guān)鍵技術(shù)研究D . 碩士學(xué)位論文,北京理工大學(xué), 2007.5 王仁明, 張春生, 賈玉坤.基于分?jǐn)?shù)階傅里葉域均衡的無線測控系統(tǒng)研究J. 航天電子對抗, 2012, 28(3):47-50.6 程乃平, 席有猷, 郝建華. 基于分?jǐn)?shù)階傅里葉變換的LFM干擾抑制算法J.裝備學(xué)院學(xué)報, 2014, 25(1):74-77.7 陳恩慶, 陶然.一種基于分?jǐn)?shù)階傅里葉

21、變換的OFDM系統(tǒng)及其均衡算法J. 電子學(xué)報, 2007, 35(3):410-4148 殷敬偉,惠俊英,蔡平等. 基于分?jǐn)?shù)階Fourier變換的水聲信道參數(shù)估計J. 系統(tǒng)工程與電子技術(shù), 2007, 29(10):1625-1627.附錄Afunction Faf = frft(f, a)% The fast Fractional Fourier Transform% input: f = samples of the signal% a = fractional power% output: Faf = fast Fractional Fourier transformerror(narg

22、chk(2, 2, nargin);f = f(:);N = length(f);shft = rem(0:N-1)+fix(N/2),N)+1;sN = sqrt(N);a = mod(a,4);% do special casesif (a=0), Faf = f; return; end;if (a=2), Faf = flipud(f); return; end;if (a=1), Faf(shft,1) = fft(f(shft)/sN; return; end if (a=3), Faf(shft,1) = ifft(f(shft)*sN; return; end% reduce

23、to interval 0.5 < a < 1.5if (a>2.0), a = a-2; f = flipud(f); endif (a>1.5), a = a-1; f(shft,1) = fft(f(shft)/sN; endif (a<0.5), a = a+1; f(shft,1) = ifft(f(shft)*sN; end% the general case for 0.5 < a < 1.5alpha = a*pi/2;tana2 = tan(alpha/2);sina = sin(alpha);f = zeros(N-1,1) ; i

24、nterp(f) ; zeros(N-1,1);% chirp premultiplicationchrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.2);f = chrp.*f;% chirp convolutionc = pi/N/sina/4;Faf = fconv(exp(i*c*(-(4*N-4):4*N-4)'.2),f);Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi);% chirp post multiplicationFaf = chrp.*Faf;% normalizing constantFaf = e

25、xp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);function xint=interp(x)% sinc interpolationN = length(x);y = zeros(2*N-1,1);y(1:2:2*N-1) = x;xint = fconv(y(1:2*N-1), sinc(-(2*N-3):(2*N-3)'/2);xint = xint(2*N-2:end-2*N+3);function z = fconv(x,y)% convolution by fftN = length(x(:);y(:)-1;P = 2nextpow2(N);z = ifft( fft(x,P) .* fft(y,P);z = z(1:N);附錄Bclc;clear all;close all;fs=30;