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

一、引言

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

變換時(shí)域頻域

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

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

連續(xù)傅里葉變換

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

傅里葉級(jí)數(shù)

連續(xù)形式的傅里葉變換其實(shí)是傅里葉級(jí)數(shù)的推廣，因?yàn)榉e分其實(shí)是一種極限形式的求和算子。對(duì)于周期函數(shù)，其傅里葉級(jí)數(shù)是存在的。

離散時(shí)域傅里葉變換

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

離散傅里葉變換

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

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

二、快速傅里葉變換FFT

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

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

基2 DIT公式推導(dǎo)

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

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

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

將N點(diǎn)的輸入序列按奇偶分為2組分別為N/2點(diǎn)的序列

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

按照蝶形信號(hào)流圖將2的結(jié)果組合為一個(gè)N點(diǎn)的DFT變換結(jié)果

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

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

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

//g++ fft.cpp -o fft -lm #include <complex> #include <cstdio> using namespace std;#define M_PI 3.14159265358979323846 //PI 雙精度//由于蝶形運(yùn)算的需要，根據(jù)奇偶坐標(biāo)將元素分到數(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所有奇下標(biāo)元素b[i] = a[i*2+1];for(int i=0; i<n/2; ++i) //copy所有偶下標(biāo)元素到數(shù)組lower-halfa[i] = a[i*2]; for(int i=0; i<n/2; ++i) //copy所有偶下標(biāo)元素(form heap)到數(shù)組upper-halfa[i+n/2] = b[i]; delete[] b; //delete heap }//N必須是2的整數(shù)次冪 //X[]存儲(chǔ)N個(gè)輸入，FFT后依舊存儲(chǔ)在X中 //由于Nyquit定理，僅僅前N/2 FFT結(jié)果有效(后N/2是鏡射) void fft2(complex<double>* X,int N) {if(N < 2) {//遞歸終止//do nothing,因?yàn)閄[0] = x[0]} else {separate(X,N); //將偶坐標(biāo)元素移至lower half,奇坐標(biāo)元素移至upper halffft2(X, N/2); //遞歸偶坐標(biāo)元素fft2(X+N/2, N/2); //遞歸奇坐標(biāo)元素//合并兩個(gè)遞歸結(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;}}} //測(cè)試 int main() {const int nSamples = 8; complex<double> x[nSamples]; // 存儲(chǔ)采樣數(shù)據(jù)complex<double> X[nSamples]; // 存儲(chǔ)FFT結(jié)果//生成測(cè)試樣本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做準(zhǔn)備}//計(jì)算fftfft2(X,nSamples);for(int i=0; i<nSamples; ++i ) printf("%0.3f \t %0.3f\n",X[i].real(),X[i].imag()); }

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

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

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

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

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
雷德算法就是用于求出這個(gè)倒序的數(shù)列。

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

代碼實(shí)現(xiàn)：

//假設(shè)N為2的整數(shù)次冪 void RaderReverse(int *arr, int N) {int j,k;//第一個(gè)和最后一個(gè)數(shù)位置不變，故不處理for(int i=1,j=N/2; i<N-1; ++i) {//原始坐標(biāo)小于變換坐標(biāo)才交換，防止重復(fù)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)?} }

時(shí)間抽取 DIT Radix－2算法

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

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

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

代碼實(shí)現(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;//第一個(gè)和最后一個(gè)數(shù)位置不變，故不處理for(int i=1,j=N/2; i<N-1; ++i) {//原始坐標(biāo)小于變換坐標(biāo)才交換，防止重復(fù)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) 級(jí)蝶形運(yùn)算 stage */ {l = 1<<i;for(j=0;j<n;j+= 2*l) /*一組蝶形運(yùn)算 group,每組group的蝶形因子乘數(shù)不同*/ {for(k=0;k<l;++k) /*一個(gè)蝶形運(yùn)算 每個(gè)group內(nèi)的蝶形運(yùn)算的蝶形因子乘數(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 } //測(cè)試 int main() {const int nSamples = 8; complex<double> x[nSamples]; // 存儲(chǔ)采樣數(shù)據(jù)complex<double> X[nSamples]; // 存儲(chǔ)FFT結(jié)果//生成測(cè)試樣本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做準(zhǔn)備}//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)容還不錯(cuò)，歡迎將生活随笔推薦給好友。