將之前學過的知識重新拾起來，仔細理解并實現。
參考：《算法導論》第30章
從頭到尾徹底理解傅里葉變換算法、上
Cooley–Tukey FFT algorithm
FFT(快速傅里葉) c語言版
數字信號處理–FFT與蝶形算法
在線MATLAB

一、引言

首先回顧信號與系統的知識，傅里葉變換是一種從時間域轉換到頻率域的變換，下面列出的幾種變體。

變換時域頻域

連續傅里葉變換(FT)	連續、非周期	非周期、離散
傅里葉級數	連續、周期	非周期、離散
離散時間傅里葉變換(DTFT)	離散、非周期	周期、連續
離散傅里葉變換(DFT)	離散、周期	周期、離散

通過表格后兩列可以發現：
時域的周期對應頻域的離散、時域的連續對應頻域的非周期，反過來也是如此。
一個域的周期對應另一個域離散、一個域的連續對應另一個域的非周期

連續傅里葉變換

連續傅里葉變換將平方可積的函數f(t)表示成復指函數的積分或者級數形式。

傅里葉級數

連續形式的傅里葉變換其實是傅里葉級數的推廣，因為積分其實是一種極限形式的求和算子。對于周期函數，其傅里葉級數是存在的。

離散時域傅里葉變換

離散時域傅里葉變換DTFT在時域是離散的，在頻域是周期的。DTFT可以看做是傅里葉級數的逆變換。

離散傅里葉變換

離散傅里葉變換DFT是連續傅里葉變換在時域和頻域都離散的。同時在時域和頻譜序列通常是有限長的，實際上將他們認為是離散周期信號的主值序列，對其進行周期延拓即可得到周期信號。
為了在科學計算和數字信號處理等領域使用計算機進行傅里葉變換，必須將輸入定義在離散點而非連續域內，且須滿足有限性或周期性條件。這是使用離散傅里葉變換DFT的原因。

離散傅里葉變換可以將連續的頻譜轉化成離散的頻譜去計算，這樣就易于計算機編程實現。時間復雜度O(n^2)。
進而引出下文，快速傅里葉變換FFT的出現，使得DFT的計算速度更快。時間復雜度O(nlgn)。

二、快速傅里葉變換FFT

快速傅立葉變換（Fast Fourier Transform，FFT）是離散傅立葉變換（Discrete Fourier transform，DFT）的快速算法，它是根據離散傅立葉變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。它對傅立葉變換的理論并沒有新的發現，但是對于在計算機系統或者說數字系統中應用離散傅立葉變換，可以說是進了一大步。
設 Xn 為 N 項的復數序列，由 DFT 變換，任一 Xi 的計算都需要 N 次復數乘法和 N -1 次復數加法，而一次復數乘法等于四次實數乘法和兩次實數加法，一次復數加法等于兩次實數加法，即使把一次復數乘法和一次復數加法定義成一次“運算”（四次實數乘法和四次實數加法），那么求出 N 項復數序列的 Xi ，即 N 點 DFT 變換大約就需要 N^2 次運算。

舉例：當 N =1024 點的時候，
使用DFT，需要 N^2 = 1048576 次運算。
使用 FFT ，利用 ωn 的周期性和對稱性，把一個 N 項序列（設 N 為偶數），分為兩個 N / 2 項的子序列，每個 N / 2點 DFT 變換需要 (N / 2)^2 次運算，再用 N 次運算把兩個 N / 2點的 DFT 變換組合成一個 N 點的 DFT 變換。這樣變換以后，總的運算次數就變成 N + 2 * (N / 2)^2 = N + N^2 / 2。當N =1024 時，總的運算次數就變成了525312 次。
二者對比可以看到，節省了大約 50% 的運算量。
而如果我們將這種“一分為二”的思想不斷進行下去，直到分成兩兩一組的 DFT 運算單元，那么N 點的 DFT 變換就只需要 N * log2N 次的運算，N = 1024 點時，運算量僅有 10240 次，是先前的直接算法的1% ，點數越多，運算量的節約就越大，這就是 FFT 的優越性。

基2 DIT公式推導

核心：FFT算法是把長序列的DFT逐次分解為較短序列的DFT。

綜合以上推導我們可以得到如下結論：一個N點的DFT變換過程可以用兩個N/2點的DFT變換過程來表示。
上式中Ek為偶數項分支的離散傅立葉變換，Ok為奇數項分支的離散傅立葉變換。其中的一個計算單元可以使用蝶形算法流圖直觀地表示出來。

那么在實現時步驟如下：

將N點的輸入序列按奇偶分為2組分別為N/2點的序列

分別對每組序列進行DFT變換得到兩組點數為N/2的DFT變換值X1和X2

按照蝶形信號流圖將2的結果組合為一個N點的DFT變換結果

三、快速傅里葉變換遞歸實現

類似于歸并排序，采用分治算法自定向下進行將問題劃分為同等規模的小問題。
FFT 的實現也可以自頂而下，采用遞歸實現。

參考了官網的代碼，進行8點FFT的實現：

//g++ fft.cpp -o fft -lm #include <complex> #include <cstdio> using namespace std;#define M_PI 3.14159265358979323846 //PI 雙精度//由于蝶形運算的需要，根據奇偶坐標將元素分到數組的前后各半部分 void separate(complex<double>* a,int n) {complex<double>* b = new complex<double>[n/2]; // get temp heapfor(int i=0; i<n/2; ++i) //copy所有奇下標元素b[i] = a[i*2+1];for(int i=0; i<n/2; ++i) //copy所有偶下標元素到數組lower-halfa[i] = a[i*2]; for(int i=0; i<n/2; ++i) //copy所有偶下標元素(form heap)到數組upper-halfa[i+n/2] = b[i]; delete[] b; //delete heap }//N必須是2的整數次冪 //X[]存儲N個輸入，FFT后依舊存儲在X中 //由于Nyquit定理，僅僅前N/2 FFT結果有效(后N/2是鏡射) void fft2(complex<double>* X,int N) {if(N < 2) {//遞歸終止//do nothing,因為X[0] = x[0]} else {separate(X,N); //將偶坐標元素移至lower half,奇坐標元素移至upper halffft2(X, N/2); //遞歸偶坐標元素fft2(X+N/2, N/2); //遞歸奇坐標元素//合并兩個遞歸結果for(int k=0; k<N/2; ++k) {complex<double> e = X[k]; //偶complex<double> o = X[k+N/2]; //奇//w是蝶形系數complex<double> w = exp( complex<double>(0,-2.0*M_PI*k/N) );X[k ] = e + w * o;X[k+N/2] = e - w * o;}}} //測試 int main() {const int nSamples = 8; complex<double> x[nSamples]; // 存儲采樣數據complex<double> X[nSamples]; // 存儲FFT結果//生成測試樣本for(int i=0; i<nSamples; ++i) {x[i] = complex<double> (0.0,0.0);x[i].real() = (double)i;x[i].imag() = 0.0;X[i] = x[i]; //拷貝至X，為FFT做準備}//計算fftfft2(X,nSamples);for(int i=0; i<nSamples; ++i ) printf("%0.3f \t %0.3f\n",X[i].real(),X[i].imag()); }

四、快速傅里葉變換迭代實現

FFT 的實現可以自頂而下，采用遞歸，但是對于硬件實現成本高，對于軟件實現都不夠高效，改用迭代較好，自底而上地解決問題。感覺和歸并排序的迭代版很類似，不過先要采用“位反轉置換”的方法把 Xi 放到合適的位置。

位反轉算法- 雷德算法

一個小算法的感覺。
拿一個0到2^n-1的自然數序列。
比方說
0 1 2 3 4 5 6 7
我們轉換為二進制狀態，那么這個序列就是

000 001 010 011 100 101 110 111

接下來我們模擬FFT的位置交換，即：

0 1 2 3 4 5 6 7
0 2 4 6 1 3 5 7
0 4 2 6 1 5 3 7

發現最終的序列變為了

000 100 010 110 001 101 011 111
雷德算法就是用于求出這個倒序的數列。

由上面的表可以看出，按自然順序排列的二進制數，其后面一個數總是比其前面一個數大1，即后面一個數是前面一個數在最低位加1并向高位進位而得到的。
而倒位序二進制數的后面一個數是前面一個數在最高位加1并由高位向低位進位而得到。
I、J都是從0開始，若已知某個倒位序J，要求下一個倒位序數：
應先判斷J的最高位是否為0。
這可與k=N/2相比較，因為N/2總是等于100的。
如果k>J，則J的最高位為0，只要把該位變為1（J與k=N/2相加即可），就得到下一個倒位序數；
如果K<=J，則J的最高位為1，可將最高位變為0（J與k=N/2相減即可）。然后還需判斷次高位，這可與k=N\4相比較，若次高位為0，
則需將它變為1（加N\4即可）其他位不變，既得到下一個倒位序數；若次高位是1，則需將它也變為0。然后再判斷下一位……

代碼實現：

//假設N為2的整數次冪 void RaderReverse(int *arr, int N) {int j,k;//第一個和最后一個數位置不變，故不處理for(int i=1,j=N/2; i<N-1; ++i) {//原始坐標小于變換坐標才交換，防止重復if(i<j) {int temp = arr[j];arr[j] = arr[i];arr[i] = temp;}k = N/2; // 用于比較最高位while(k <= j) { // 位判斷為1j = j-k;// 該位變為0k = k/2;// 用于比較下一高位}j = j+k;// 判斷為0的位變為1} }

時間抽取 DIT Radix－2算法

該算法的特征是：
1）輸入序列順序位反轉，輸出所有頻率值都是按順序出現的。
2）計算可以“就地”完成，也就是蝶形所使用的存儲位置可以被重寫。
3）從圖中我們可以看到對于點數為N = 2^L的FFT運算，可以分解為L階蝶形圖級聯，每一階蝶形圖內又分為M個蝶形組，每個蝶形組內包含K個蝶形。（迭代算法實現：根據這一點我們就可以構造三階循環來實現蝶形運算。編程過程需要注意旋轉因子與蝶形階數和蝶形分組內的蝶形個數存在關聯。）

快速傅里葉變換迭代算法實現

該代碼在上文的基礎上，參考了FFT(快速傅里葉) c語言版的思想。

代碼實現：

//g++ fft2.cpp -o fft2 -lm #include <complex> #include <cstdio> using namespace std;#define M_PI 3.14159265358979323846 //PI 雙精度 //假設N為2的整數次冪 void RaderReverse(complex<double> *arr, int N) {int j,k;//第一個和最后一個數位置不變，故不處理for(int i=1,j=N/2; i<N-1; ++i) {//原始坐標小于變換坐標才交換，防止重復if(i<j) {complex<double> temp = arr[j];arr[j] = arr[i];arr[i] = temp;}k = N/2; // 用于比較最高位while(k <= j) { // 位判斷為1j = j-k;// 該位變為0k = k/2;// 用于比較下一高位}j = j+k;// 判斷為0的位變為1} } void fft(complex<double> *x,int n) {int i=0,j=0,k=0,l=0;complex<double> up,down,product;RaderReverse(x,n);//w是蝶形系數complex<double>* W = new complex<double>[n]; for(int i=0;i<n;++i) {W[i] = exp( complex<double>(0,-2.0*M_PI*i/n) );//printf("%0.3f \t %0.3f\n",W[i].real(),W[i].imag());}for(i=0; i<log(n)/log(2); ++i) /*log(n)/log(2) 級蝶形運算 stage */ {l = 1<<i;for(j=0;j<n;j+= 2*l) /*一組蝶形運算 group,每組group的蝶形因子乘數不同*/ {for(k=0;k<l;++k) /*一個蝶形運算 每個group內的蝶形運算的蝶形因子乘數成規律變化*/ {product = x[j+k+l]*W[n*k/2/l];up = x[j+k] + product;down = x[j+k] - product;x[j+k] = up;x[j+k+l] = down;}}}delete[] W; //delete W } //測試 int main() {const int nSamples = 8; complex<double> x[nSamples]; // 存儲采樣數據complex<double> X[nSamples]; // 存儲FFT結果//生成測試樣本for(int i=0; i<nSamples; ++i) {x[i] = complex<double> (0.0,0.0);x[i].real() = (double)i;x[i].imag() = 0.0;X[i] = x[i]; //拷貝至X，為FFT做準備}//RaderReverse(X,nSamples);//for(int i=0; i<nSamples; ++i ) // printf("%0.3f \t %0.3f\n",X[i].real(),X[i].imag());fft(X,nSamples);for(int i=0; i<nSamples; ++i ) printf("%0.3f \t %0.3f\n",X[i].real(),X[i].imag()); }

總結

以上是生活随笔為你收集整理的快速傅里叶变换学习及C语言实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。