優雅的FFT演算法

2022-08-09|由 APPLE頻道 發表于 母嬰

優雅的FFT演算法

簡 介：

利用FFT演算法實現快速傅立葉變換， 在理論、工程中具有非常廣泛的應用。 除了能夠在合適的計算平臺完成FFT演算法，同時還需要注意到它在頻譜分析中可能帶來的

頻率混疊

以及

頻率洩露

等問題。

關鍵詞： FFT，演算法實現

01 Python演算法

今天下午的訊號與系統， 給同學們介紹了

離散傅立葉

變換的基本應用， 並且介紹了快速傅立葉變換（FFT）的主要思想與演算法。 FFT演算法因其優異的效能和廣泛的應用， 堪稱資訊處理領域的

原子武器

。 實現FFT程式語言很多， 比較來比較去， 利用Python語音所描述的該演算法最為簡明和優雅。

1。1 FFT演算法程式碼

下面的程式碼是在

The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?

影片中給出的 FFT 遞迴演算法形式， 最大精度反映了FFT演算法核心。

這個程式碼實現了DIF（時域抽取快速傅立葉變換）， 利用遞迴定義，將FFT核心演算法中的

分而治之

體現的淋漓盡致， 突出了遞迴核心中的核心思想。

def FFT（P）： n = len（P） if n * 1： return P ye = FFT（P［0：：2］） yo = FFT（P［1：：2］） y = ［0］*n w = exp（-1j*2*pi/n） for j in range（n//2）： yow = w**j * array（yo） y［j］ = ye［j］ + yow［j］ y［j+n//2］ = ye［j］ - yow［j］ return y

利用Python語音中對於陣列切片操作語法， 還可以將上面FFT演算法中的迴圈部分都替換成關於陣列的操作， 使得實際運算速度得到提高。

def FFT1（P）： n = len（P） if n * 1： return P ye = FFT（P［0：：2］） yo = FFT（P［1：：2］） w = exp（-1j*2*pi/n）**array（list（range（n//2））） yow = w*yo y = ［0］*n y［：n//2］ = ye + yow y［n//2：］ = ye - yow return y

1。2 FFT 演算法測試

為了測試演算法的有效性， 下面對於一個方波訊號計算對應的FFT結果。

測試演算法程式碼如下：

LEN = 1024oneLEN=10p1 = ［1］*oneLEN+［0］*（LEN-oneLEN）y = FFT（p1）plt。plot（abs（array（y））， label=‘abs（FFT）’）plt。plot（p1， label=‘Data’）plt。xlabel（“y”）plt。ylabel（“abs（FFT（y））”）plt。grid（True）plt。legend（loc=‘upper right’）plt。tight_layout（）plt。show（）

下面是測試利用Python語言實現的FFT演算法計算結果。

優雅的FFT演算法

▲ 圖1。2。1 利用Python語音實現的FFT演算法測試結果

02 其它語言FFT

FFT演算法貴在計算效率，前面使用Python實現FFT，雖然形式上優雅，但實際執行效率不高。 提高執行效率，還是需要使用編譯語言。

2。1 Fortran FFT演算法

在我上大學期間所學的程式語言為Fortran， 估計現在沒有多少同學學習這個演算法語言。 下面給出了利用Fortran語言實現的FFT演算法程式。

演算法整體上包括有兩個階段：

第一個階段實現了對輸入資料進行

倒讀順序

排列；

第二階段利用三重迴圈實現了分組蝶形運算。

當然了，時過三十年再看Fortran感覺十分酸爽， 但它簡練語言和執行高效還是讓我們回憶起當年程式設計時所感覺到的快樂。

優雅的FFT演算法

▲ 圖 Fortran 語言實現的FFT演算法

2。2 C語言FFT演算法

下面是在網路上博文

C++ Program to Compute Discrete Fourier Transform using Fast Fourier Transform Approach[1]

給出的FFT演算法， 沒有對其功能進行測試。 相比於前面利用Python，Fortran來看， C語言實現FFT就顯得非常囉嗦了。

#include #include #include #include using namespace std；unsigned int bitReverse（unsigned int x， int log2n） { int n = 0； int mask = 0x1； for （int i = 0； i < log2n； i++） { n <<= 1； n |= （x & 1）； x >>= 1； } return n；}const double PI = 3。1415926536；templatevoid fft（Iter_T a， Iter_T b， int log2n） { typedef typename iterator_traits：：value_type complex； const complex J（0， 1）； int n = 1 << log2n； for （unsigned int i = 0； i < n； ++i） { b［bitReverse（i， log2n）］ = a［i］； } for （int s = 1； s <= log2n； ++s） { int m = 1 << s； int m2 = m >> 1； complex w（1， 0）； complex wm = exp（-J * （PI / m2））； for （int j = 0； j < m2； ++j） { for （int k = j； k < n； k += m） { complex t = w * b［k + m2］； complex u = b［k］； b［k］ = u + t； b［k + m2］ = u - t； } w *= wm； } }}

※ 總 結 ※

利用FFT演算法實現快速傅立葉變換， 在理論、工程中具有非常廣泛的應用。 除了能夠在合適的計算平臺完成FFT演算法，同時還需要注意到它在頻譜分析中可能帶來的

頻率混疊

以及

頻率洩露

等問題。

參考資料

［1］

C++ Program to Compute Discrete Fourier Transform using Fast Fourier Transform Approach

：

https：//www。sanfoundry。com/cpp-program-compute-discrete-fourier-transform-using-fast-fourier-transform-approach/

TAG： FFT演算法intpltpython