稀疏快速傅里葉變換算法理論基礎(chǔ)概述目錄TOC\o"1-3"\h\u234741.1理論基礎(chǔ) 1241281.1.1傅里葉變換 1234271.1.2離散傅里葉變換 246431.1.3快速傅里葉變換 3275161.1.4時間復(fù)雜度 336991.2稀疏快速傅里葉變換 4257191.1.1誤差約束準則 44801.1.2稀疏快速傅里葉變換理論框架 5303931.1.3頻譜重排 6283231.1.4加窗 826761.1.5降采樣 9223991.1.6定位 10198011.1.7估值 11229251.1.8外循環(huán) 111.1理論基礎(chǔ)1.1.1傅里葉變換傅里葉變換是由法國數(shù)學(xué)家、物理學(xué)家傅里葉最早在1807年發(fā)表的論文提出的。傅里葉認為任何一個周期函數(shù)都可以由一組適當?shù)恼仪€組合而成。傅里葉所作出的早期研究為傅里葉級數(shù)的理論開辟了研究方向。傅里葉級數(shù)表示如下：(2-1)上式中n取正整數(shù)，直流分量、余弦系數(shù)和正弦系數(shù)由以下公式計算：(2-2)周期函數(shù)能夠用傅里葉級數(shù)進行分解的前提是滿足Dirichlet條件，另外一種表示方式如下：(2-3)其中：，。將作為橫軸，作為縱軸畫出頻譜圖得到的結(jié)果稱為幅度譜，而作為縱軸畫出的頻譜圖稱為相位譜。指數(shù)形式的傅里葉級數(shù)可表示為：(2-4)其中為指數(shù)形式的傅里葉系數(shù)。對于非周期信號，人們?nèi)匀幌Ｍ軌蛲ㄟ^傅里葉級數(shù)分解的方式來方便對信號的分析。因此對傅里葉級數(shù)進行推廣，令上述式子中的信號的周期趨向于無窮大，那么從系數(shù)的表達式可以看出頻譜的間隙就會趨向于無窮小，最終得到一個連續(xù)的頻譜圖像。此時，傅里葉系數(shù)經(jīng)極限運算后也變成一個連續(xù)的函數(shù)：(2-5)傅里葉變換可表示為：(2-6)逆變換：(2-7)用“”來表示傅里葉變換運算符號，令為的傅里葉變換，則有以下性質(zhì)：對稱性：；線性：；尺度變換：；位移性質(zhì)：；。1.1.2離散傅里葉變換在使用計算機處理信號數(shù)據(jù)時，則有另一套運算規(guī)則——離散傅里葉變換(discretefouriertransform,DFT)。它可以對離散信號進行變換，也可以對連續(xù)信號經(jīng)過采樣后得到的采樣信號進行變換。其定義為：，(2-8)在式中，的絕對值是對應(yīng)頻譜中相應(yīng)頻率的諧波振幅，對應(yīng)的頻率為，是采樣間隔。1.1.3快速傅里葉變換一般來說，對于能夠大大減少離散傅里葉變換的運算時間的算法，我們統(tǒng)稱為快速傅里葉變換。其中最經(jīng)典的也是被提及到最多的是1965年Cooley等人提出的變換算法。這里要介紹的就是Cooley等人提出的快速傅里葉變換。在處理信號時，對信號長度為n點的時間序列直接應(yīng)用離散傅里葉變換，通常會因為對每個點都直接運算而導(dǎo)致運算量的大幅增加。但是考慮到信號具有周期性的特點，我們可以對信號進行分離。在快速傅里葉變換中，將N點的序列分離成兩個部分。改變公式中的加法順序，讓下標為奇數(shù)和偶數(shù)分別寫在兩個求和符號中，最終得到兩個序列。顯然，得到的兩個序列長度為，根據(jù)離散傅里葉變換的線性性質(zhì)，可以了解到如果對這兩個序列分別進行離散傅里葉變換，得到兩個結(jié)果之和就是原來信號序列的離散傅里葉變換。繼續(xù)對分解之后的序列進行分解，分別得到兩個序列，按照這個規(guī)則繼續(xù)分解下去，最終可以得到一個N行1列的序列（信號長度n為2的整數(shù)次冪）。在離散傅里葉變換表達式中的等效表現(xiàn)為：(2-9)出于簡化公式方便推導(dǎo)的目的，令前后兩個序列和分別為和，那么上式變?yōu)椋?2-10)根據(jù)離散傅里葉變換的周期運算又可以得到：(2-11)通過這樣改變運算方式和順序，可以看到只要求出在[0，N/2-1]區(qū)間內(nèi)各個整數(shù)對應(yīng)的離散傅里葉變換系數(shù)，那么就可以將整個區(qū)間的系數(shù)算出來。得到同樣的結(jié)果，但是在快速傅里葉變換中只需要計算一半的次數(shù)。1.1.4時間復(fù)雜度時間復(fù)雜度的定義為運行一次算法所需要花費時間的量級。在實際評測算法的性能過程中，算法中計算的次數(shù)、整體運算規(guī)模、運行環(huán)境的配置等因素都會影響時間的消耗，所以需要借助一個量級來評估，這就是時間復(fù)雜度的意義。實驗時，控制其他的參數(shù)變量，最主要的能夠明顯影響時間復(fù)雜度的因素是輸入信號的長度。由于算法具有一定的復(fù)雜性，另外信號長度也會因為其較大的數(shù)量級而導(dǎo)致我們無法準確地計算出確切的運行時間與長度的關(guān)系，只能得到一個模糊的與長度相關(guān)的趨向?！癘(·)”描述的是函數(shù)數(shù)量級逼近的上界，是1894年由PaulBachmann引入的，通常用它來作為評估函數(shù)的增長速率，在時間復(fù)雜度方面我們也采取這樣的評估方式。信號長度通常是一個十分巨大的數(shù)，因此在比較不同的時間復(fù)雜度就可以按照幾項時間復(fù)雜度中的最大的一項來比較得到結(jié)果。目前經(jīng)常見到的時間復(fù)雜度的比較順序如下：借助這些可以了解到快速傅里葉變換的時間復(fù)雜度為，遠遠小于離散傅里葉變換由于沒有借助信號的性質(zhì)運算而得到的時間復(fù)雜度。1.2稀疏快速傅里葉變換相比離散傅里葉變換的O()，快速傅里葉變換的確大大加快了運算的速度。但是隨著需要處理的信號與數(shù)據(jù)量越來越大，而技術(shù)實現(xiàn)對于實時性的要求也與日俱增，O()的時間復(fù)雜度也逐漸不能滿足對算速的要求。因此本文探索時間復(fù)雜度更低的稀疏快速傅里葉變換，通過分析原理得到算法的實現(xiàn)關(guān)鍵技術(shù)和在時間復(fù)雜度上稀疏快速傅里葉變換更優(yōu)的原因。1.1.1誤差約束準則由于稀疏快速傅里葉變換是受到各個領(lǐng)域需要處理的信號具有稀疏這樣的共同特點啟發(fā)而逐步發(fā)展起來的算法優(yōu)化的成果，它并不適用于對所有信號的處理。它所處理的信號必須滿足在頻譜中是稀疏的特點。因此我們首先定義稀疏信號：當一個信號進行離散傅里葉變換后得到的頻譜中只有少數(shù)非零的傅里葉系數(shù)，或者其他傅里葉系數(shù)與這幾個少數(shù)的傅里葉系數(shù)相比可忽略，那么就稱信號是稀疏的。不可忽略的傅里葉系數(shù)數(shù)量為，定義為稀疏度，是遠小于信號長度的。這些傅里葉系數(shù)的值為，。要判斷需要處理的信號是否為稀疏信號，在運算之前必須進行預(yù)評估才能實現(xiàn)良好的效果。一般情況下稀疏信號需要滿足一定的準則才能令稀疏傅里葉變換有效，這里設(shè)為變換得到的近似結(jié)果，為傅里葉變換，設(shè)定一個大于1的數(shù)，將它作為近似的誤差因子，是一個稀疏度為，長度為的稀疏序列。要求令范數(shù)取最小值。那么信號需要滿足范數(shù)準則，即：(2-12)在實際情況中，如果我們把那些不包含重要信息的小值頻點忽略為0，僅僅保留不可被忽略的點，那么就可以被稀疏表示，而稀疏度遠小于信號的長度，因此計算結(jié)果就可以得到非常清晰的頻譜。而這樣的計算能夠使得時間復(fù)雜度近似與信號長度呈現(xiàn)線性關(guān)系。這類算法也被稱為亞線性算法。對于稀疏快速傅里葉變換算法來說，要得到更良好的頻點的幅度估計，需要滿足更加嚴格的約束條件，對此我們引進一個可以對每一個幅值都能進行控制，令輸出結(jié)果更加合理的范數(shù)準則：(2-13)式中其他參數(shù)向量同上，不同的是增加了一個大于零的常數(shù)和精度參數(shù)。對于我們要求它要等于 。而這項準則之所以更嚴格，是因為它是用頻譜中各項幅值誤差的平方和來約束整體，因此能夠達到稀疏快速傅里葉變換后每個值都能得到良好估計的效果。1.1.2稀疏快速傅里葉變換理論框架稀疏快速傅里葉變換算法進行運算時，核心思想在于快速傅里葉變換所用到的點數(shù)減少，保留重要信息而忽略近零的頻點，降低信號的復(fù)雜程度。而通過給信號“分筐”的操作，將大值點隨機分列到“筐”中，再將每個“筐”視為采樣點進行傅里葉變換，就可以達到這樣的效果，如圖2-1。對頻譜分筐的操作是通過降采樣進行的，而在運算過程中需要構(gòu)建一個在原信號N維度到分筐后信號B維度的映射。本文采取的是構(gòu)建一種名為“哈?！钡挠成渌惴ǎ瑢⒃盘栃蛄凶鴺藢?yīng)到“筐”中。由于本身是概率性算法，在映射過程中有概率會發(fā)生不同大值點被分到同一個“筐”中，進而導(dǎo)致變換失敗。而且在實際應(yīng)用中，大值頻點在頻譜中密集存在的情況非常多，所以需要在“分筐”前對原信號的這些大值點進行隨機重新排列，來避免上述情況帶來的算法失效問題。重排后的信號仍然不能直接用于降采樣，在此之前需要加上窗函數(shù)來截斷信號，避免發(fā)生由于直接截斷而產(chǎn)生頻譜泄露問題，造成很大的誤差。圖2-1分筐示意圖降采樣FFT后完成了“分筐”操作，這時我們需要將散列到“筐”中的大值頻點對應(yīng)到原函數(shù)頻譜中的坐標和幅值，方便后續(xù)進行的恢復(fù)頻譜操作?；謴?fù)坐標需要進行定位操作，主要是通過哈希逆映射對其進行定位。而估值則是通過一定的函數(shù)對每個“筐”中的頻點進行對應(yīng)的幅值估計。得到了定位操作返回的坐標結(jié)果和估值操作返回的幅值估計結(jié)果，再進行相應(yīng)的頻譜重構(gòu)就可以得到稀疏快速傅里葉變換得到的結(jié)果。而由于重排的概率性和單次定位估值的低置信度，我們需要把從重排到估值這些操作按照信號長度來進行相應(yīng)次數(shù)的循環(huán)，最終得到的統(tǒng)計學(xué)結(jié)果才能效果更好。稀疏快速傅里葉變換雖然有概率重構(gòu)失敗，但是在對信號稀疏度進行誤差約束的預(yù)評估之后，進行足夠的整體循環(huán)，能夠讓重構(gòu)成功的概率盡量接近1。綜上所述，對于稀疏快速傅里葉變換算法理論的整體框架可以用圖2-2來說明：圖2-2稀疏快速傅里葉變換理論框架1.1.3頻譜重排在稀疏快速傅里葉變換的過程中，稀疏信號的頻譜中大值頻點距離過近對于“分筐”操作影響較大，會導(dǎo)致很大的沖突問題影響后續(xù)重構(gòu)信號頻譜的準確度。然而在研究過程中發(fā)現(xiàn)大部分信號都會有大值頻點聚集的現(xiàn)象。因此需要對信號進行頻譜重排。但是算法的輸入信號得到的是時域的信息，不能直接進行頻譜上的操作，所以重排的工作必須在時域上進行。根據(jù)傅里葉變換的性質(zhì)，我們可以通過時域上的位移和縮放來間接地對頻譜進行重排。傅里葉變換有如下性質(zhì)可供利用：位移性質(zhì)：若，則；縮放性質(zhì)：若，則；其中，，坐標溢出時對取模，在隨機變換后要對輸入信號進行時域上的完全遍歷，因此設(shè)定兩個隨機整數(shù)作為參數(shù)，分別是和，是一個與互質(zhì)的正整數(shù)，。設(shè)為原信號時域坐標，則重排后時域坐標對應(yīng)為。根據(jù)完全剩余理論，這里重排的信號對原信號時域坐標進行了遍歷。設(shè)重排后的信號為，則原信號與重排信號之間的關(guān)系可以表示為：(2-14)通過推導(dǎo)可以看到頻譜的關(guān)系為：(2-15)在式(2-15)中可以得出結(jié)論，對輸入信號進行時域重新排列后頻譜的位置滿足的對應(yīng)關(guān)系，另外是在相位方面影響信號頻譜的，并不會對信號的幅度譜產(chǎn)生影響。所以可以認為成功地對信號的幅度譜進行了重新隨機的排列。下面通過一個簡單信號的例子來驗證重排的實際效果：設(shè)定一個長度為256，由頻率為15和65的兩個正弦信號合成的數(shù)字信號作為重排的輸入信號，其中通過在MATLAB中構(gòu)建函數(shù)隨機生成一個符合條件的整數(shù)，為了方便起見取1。如圖2-3上側(cè)為信號時域圖像，下側(cè)為信號頻域圖像：圖2-3示例信號的時域(上)和頻域(下)進行頻譜重排的處理之后可以得到如圖2-4所示的時域效果和圖2-5所示的頻域效果：圖2-4原信號與重排信號的時域?qū)Ρ葓D2-5原信號與重排信號的頻域?qū)Ρ仍趫D2-4中可以大致看出在頻譜重排后信號時域的幅度沒有發(fā)生改變，而圖2-5中則能夠更加直觀地看出頻譜重排產(chǎn)生的效果。1.1.4加窗計算機在處理信號的時候只能處理有限長數(shù)字信號，而實際應(yīng)用當中很難有理想的有限長數(shù)字信號。因此面對無限長連續(xù)信號時需要進行一定的技術(shù)處理。處理的方式主要包括采樣和截斷。截斷是指只分析無限長信號的有限部分，在信號處理過程中是通過用窗函數(shù)序列與信號序列對應(yīng)相乘。然而直接截斷會對處理后的信號產(chǎn)生很嚴重的影響。具體影響表現(xiàn)在導(dǎo)致原來單一窄帶的譜改變?yōu)閷拵У倪B續(xù)頻譜，使得頻譜的分辨能力大大降低。造成這些影響的原因主要是因為加窗的操作是時域相乘，根據(jù)傅里葉變換的性質(zhì)，在頻譜中是進行卷積運算的。在卷積運算的過程中，由于窗函數(shù)自身不能做到理想化，在主瓣邊緣會產(chǎn)生一些旁瓣，這些旁瓣與信號進行卷積后導(dǎo)致結(jié)果變得失真。某些不理想的窗自身旁瓣幅度大，對信號的破壞能力更強，甚至會導(dǎo)致真實譜中的某些低值被覆蓋消失，另外還會在頻譜中創(chuàng)造出一些預(yù)期之外的譜峰。為了提高頻譜的分辨能力，需要選擇主瓣更窄的窗函數(shù)進行加窗。如果主瓣達不到相應(yīng)參數(shù)的要求，那么就會在一定程度上降低分辨能力。對于窗函數(shù)旁瓣也有一定的要求，比如旁瓣衰減速率和旁瓣峰值電平。前者要求速率盡量更快以達到旁瓣變窄，降低影響；后者則要求峰值更低，削弱旁瓣對頻譜產(chǎn)生的影響。以上的指標一般作為對窗函數(shù)性能評估的重要指標。而對于優(yōu)化算法計算量來說，形式是否簡單也應(yīng)當作為考慮的因素。稀疏快速傅里葉算法對信號也是采用部分點進行變換，因此對信號加窗的環(huán)節(jié)十分關(guān)鍵。為了避免對頻譜造成能量的泄露，選擇一個合適的窗函數(shù)是十分必要的。所以對窗函數(shù)的參數(shù)研究和形式選擇上需要更深入的討論。為了使截斷頻譜更加平滑，在稀疏快速傅里葉變換中選擇平坦窗函數(shù)對信號進行處理。為了達到上文中性能指標的標準，構(gòu)建的平坦窗函數(shù)在時域中相比長度來說只有一小部分點幅度非零。而在頻域中則通帶處有較大值，其余近乎為零。平坦窗函數(shù)的構(gòu)建通常是由高斯窗函數(shù)與sinc函數(shù)相乘構(gòu)建，也可以通過標準窗函數(shù)與高斯函數(shù)通過卷積運算獲得。窗函數(shù)長度設(shè)為，滿足(2-16)其中為窗函數(shù)頻譜幅值，為大于的整數(shù)，稱為阻帶截斷因子，為振蕩紋波。則認為是一個標準窗函數(shù)，參數(shù)為。在實際的生成中，只需要對滿足條件的高斯窗進行截斷即可。截斷位置在。將得到的標準窗與Box-car函數(shù)進行卷積，得到平坦窗函數(shù)。平坦窗函數(shù)參數(shù)為。圖2-6平坦窗函數(shù)時域圖2-7平坦窗函數(shù)頻域如圖2-6和圖2-7給出了在示例信號長度為2048時平坦窗函數(shù)的時域和頻域的圖像。平坦窗函數(shù)的參數(shù)不受輸入信號的影響。對于輸入信號，在進行變換之前已經(jīng)進行了預(yù)評估，先行設(shè)置好信號的稀疏度、需要分筐數(shù)量。而窗函數(shù)的參數(shù)只與這些預(yù)定系數(shù)有關(guān)，因此輸入信號的總體性質(zhì)不變的話不需要更改窗函數(shù)參數(shù)的，總體的性質(zhì)改變后窗函數(shù)才需要再另行改變參數(shù)。1.1.5降采樣降采樣技術(shù)是稀疏快速傅里葉變換對計算次數(shù)突破的關(guān)鍵。但是由于獲取的輸入信息只有信號的時域，不能對它的頻域進行直接降采樣。因此需要對信號進行時域上的混疊?；殳B依據(jù)的是信號傅里葉變換的性質(zhì)求：(2-17)根據(jù)推導(dǎo)得到時域信號關(guān)系：(2-18)其中的和都在0到范圍內(nèi)。由式(2-18)可以看出，對于信號時域的混疊能夠達到在頻域中采樣間隔的降低。后續(xù)的處理中，算法不再是對原來信號的頻譜進行操作。而是對于降采樣的頻譜，借助哈希算法進行處理，也就是前文中所說的“分筐”。通過混疊成功避免了對信號直接的快速傅里葉變換，雖然后面也需要對降采樣信號進行變換，但是由于(一般情況下)是遠遠小于信號的長度的，在這里應(yīng)用了快速傅里葉變換，將時間復(fù)雜度與信號長度的關(guān)聯(lián)改變?yōu)榕c設(shè)定的“筐”數(shù)的關(guān)聯(lián)，將快速傅里葉變換的時間復(fù)雜度通過降采樣運算降低到了達到了算速的優(yōu)化目的。1.1.6定位定位的作用在于確定非零頻點的位置。定位嚴格來說是在稀疏快速傅里葉變換中進行定位循環(huán)。它不是單純的對信號頻域中的大值頻點進行定位，而是包括前文中所描述的頻譜隨機重排、加窗截斷和降采樣FFT的模塊構(gòu)成了哈希映射，最后再通過哈希逆映射進行對大值頻點坐標的處理。哈希映射表示的是通過一個函數(shù)(哈希函數(shù))對原信號頻點的坐標向降采樣譜的坐標映射，其表達式為：(2-19)從而實現(xiàn)將其散列至在降采樣譜中預(yù)先分好的“筐”中，而這時降采樣譜中每個“筐”的幅值取散列至筐中的點的幅值之和。接下來是對信號大值頻點的操作：先對降采樣頻譜中的頻點進行排序，生成一個個大值頻點集中在一側(cè)的向量集，同時獲取排序前的坐標。收集被集中到一側(cè)的這個大值頻點的坐標放入數(shù)集中。經(jīng)過哈希逆映射將降采樣譜中的頻點向原譜映射，得到在數(shù)集中大值頻點坐標在原譜中對應(yīng)的坐標的集合。按照出現(xiàn)次數(shù)對這些坐標進行排序，其中出現(xiàn)次數(shù)最多的個坐標便可認為是頻點在原譜中的位置。在逆映射的過程中往往會出現(xiàn)多個相近的原譜坐標，因此需要設(shè)置一個范圍，該范圍可以通過以下函數(shù)實現(xiàn)：(2-20)其中為降采樣譜的坐標。通過逆映射得到的坐標在這樣的范圍內(nèi)就可以認定為大值頻點，取其中出現(xiàn)次數(shù)最多的個坐標得到輸出數(shù)集，等待后續(xù)的估值處理。如下圖2-8展示了定位循環(huán)的結(jié)構(gòu)：圖2-8定位循環(huán)結(jié)構(gòu)1.1.7估值經(jīng)定位獲得了大值頻點對應(yīng)原譜