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

一、引言

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

變換時域頻域

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

通過表格后兩列可以發(fā)現(xiàn)：
時域的周期對應(yīng)頻域的離散、時域的連續(xù)對應(yīng)頻域的非周期，反過來也是如此。
一個域的周期對應(yīng)另一個域離散、一個域的連續(xù)對應(yīng)另一個域的非周期

連續(xù)傅里葉變換

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

傅里葉級數(shù)

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

離散時域傅里葉變換

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

離散傅里葉變換

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

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

二、快速傅里葉變換FFT

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

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

基2 DIT公式推導

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

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

那么在實現(xiàn)時步驟如下：

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

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

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

三、快速傅里葉變換遞歸實現(xiàn)

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

參考了官網(wǎng)的代碼，進行8點FFT的實現(xiàn)：

//g++ fft.cpp -o fft -lm #include <complex> #include <cstdio> using namespace std;#define M_PI 3.14159265358979323846 //PI 雙精度//由于蝶形運算的需要，根據(jù)奇偶坐標將元素分到數(shù)組的前后各半部分 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所有偶下標元素到數(shù)組lower-halfa[i] = a[i*2]; for(int i=0; i<n/2; ++i) //copy所有偶下標元素(form heap)到數(shù)組upper-halfa[i+n/2] = b[i]; delete[] b; //delete heap }//N必須是2的整數(shù)次冪 //X[]存儲N個輸入，FFT后依舊存儲在X中 //由于Nyquit定理，僅僅前N/2 FFT結(jié)果有效(后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); //遞歸奇坐標元素//合并兩個遞歸結(jié)果for(int k=0; k<N/2; ++k) {complex<double> e = X[k]; //偶complex<double> o = X[k+N/2]; //奇//w是蝶形系數(shù)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]; // 存儲采樣數(shù)據(jù)complex<double> X[nSamples]; // 存儲FFT結(jié)果//生成測試樣本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()); }

四、快速傅里葉變換迭代實現(xiàn)

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

位反轉(zhuǎn)算法- 雷德算法

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

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

發(fā)現(xiàn)最終的序列變?yōu)榱?/p>

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

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

代碼實現(xiàn)：

//假設(shè)N為2的整數(shù)次冪 void RaderReverse(int *arr, int N) {int j,k;//第一個和最后一個數(shù)位置不變，故不處理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;// 該位變?yōu)?k = k/2;// 用于比較下一高位}j = j+k;// 判斷為0的位變?yōu)?} }

時間抽取 DIT Radix－2算法

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

快速傅里葉變換迭代算法實現(xiàn)

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

代碼實現(xiàn)：

//g++ fft2.cpp -o fft2 -lm #include <complex> #include <cstdio> using namespace std;#define M_PI 3.14159265358979323846 //PI 雙精度 //假設(shè)N為2的整數(shù)次冪 void RaderReverse(complex<double> *arr, int N) {int j,k;//第一個和最后一個數(shù)位置不變，故不處理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;// 該位變?yōu)?k = k/2;// 用于比較下一高位}j = j+k;// 判斷為0的位變?yōu)?} } 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是蝶形系數(shù)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的蝶形因子乘數(shù)不同*/ {for(k=0;k<l;++k) /*一個蝶形運算 每個group內(nèi)的蝶形運算的蝶形因子乘數(shù)成規(guī)律變化*/ {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]; // 存儲采樣數(shù)據(jù)complex<double> X[nSamples]; // 存儲FFT結(jié)果//生成測試樣本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()); }

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。