通俗易懂的講解FFT的讓你快速了解FFT

相信網(wǎng)上現(xiàn)在有很多關(guān)于FFT的教程，我曾經(jīng)也參閱了很多網(wǎng)上的教程，感覺都不怎么通俗易懂。在基本上的研究FFT，并且通過編程的形式實(shí)現(xiàn)之后。我決定寫一篇通俗易懂的關(guān)于FFT的講解。因此我在接下來的敘述中盡量非常通俗細(xì)致的講解。

本人最早知道傅里葉變換的時(shí)候是沉迷于音樂的頻譜跳動(dòng)無法自拔，當(dāng)時(shí)就很想做一個(gè)音樂頻譜顯示器。搜閱了很多資料之后，才了解到傅里葉變換，和FFT。當(dāng)然這都是以前的事情了，經(jīng)過了系統(tǒng)的學(xué)習(xí)+2個(gè)星期的研究，自制了一個(gè)FFT的算法，不敢說速度上的優(yōu)勢(shì)，但是個(gè)人認(rèn)為是一種通俗易懂的實(shí)現(xiàn)方法。經(jīng)過實(shí)際的VC++模擬實(shí)驗(yàn)、和STM32跑的也很成功。

首先，要會(huì)FFT，就必須要對(duì)DFT有所了解，因?yàn)閮烧咧g本質(zhì)上是一樣的。在此之前，先列出離散傅里葉變換對(duì)(DFT):

,k=0,1,…N-1

n=0,1…N-1

其中：

但是FFT之所以稱之為快速傅里葉變換，就利用了以下的幾個(gè)性質(zhì)（重中之重?。?/p>

周期性：

對(duì)稱性：

可約性：

先把這仨公式放到這，接下來會(huì)用到。

根據(jù)這幾個(gè)特性，就可以將一個(gè)長的DFT運(yùn)算分解為若干短序列的DFT運(yùn)算的組合，從而減少運(yùn)算量。在這里，為了方便理解，我就用了一個(gè)按時(shí)間抽取的快速傅里葉變換（DIT-FFT）的方法。

首先，將一個(gè)序列x(n)一分為二

對(duì)于,k=0,1,…N-1? ?設(shè)N是2的整次冪 也就是N=2^M

將x(n)按照奇偶分組

x(2r)=x1(r)

x(2r+1)=x2(r) r=0,1,…..(N/2)-1

那么將DFT也分為兩組來預(yù)算

(第一項(xiàng)是偶第二項(xiàng)是奇)

這個(gè)時(shí)候，我們上邊提的三個(gè)性質(zhì)中的可約性就起到作用了：

回顧一下：

那么這個(gè)式子就可以化為：

(式1-1)

則X1(k)和X2(k)都是長度為N/2的序列 x1(k)和x2(k)的N/2點(diǎn)的離散傅里葉變換

其中：

K=0,1,2…N/2-1

至此，一個(gè)N點(diǎn)的DFT就被分解為2個(gè)N/2的DFT。但是X1(k),和X2(k)只有N/2個(gè)點(diǎn)，也就是說式子（1-1）只是DFT前半部分。要求DFT的后半部分可以利用其對(duì)稱性求出后半部分為：

(式1-2)

那么式（1-1）和（1-2）就可以專用一個(gè)蝶形信號(hào)流圖符號(hào)來表示。如圖：

以N=8為例，可以用下圖表示：

通過這樣的分解，每一個(gè)N/2點(diǎn)DFT只需(N^2)/4次復(fù)數(shù)相乘計(jì)算，明顯的節(jié)省了運(yùn)算量。

以此類推，繼續(xù)將已經(jīng)得出的X1(k)和X2（k）按照奇偶分解，還是和上邊一樣的方法。那么就得出了百度上都可以找到的一大推的這個(gè)圖片了。（笑）

對(duì)于這張圖片我想強(qiáng)調(diào)的一點(diǎn)就是第二階蝶形運(yùn)算的時(shí)候，WN0之后不應(yīng)該是WN1嗎，為什么是W2N了，這個(gè)問題之前困擾了我一段時(shí)間，但是不要忘了，前者的 W0N的展開是W0N/2因?yàn)榇藭r(shí)N已經(jīng)按照奇偶分開了，所以是N/2而W2N/2也就是W2N是根據(jù)可約性得出的，在這里不能混淆.。

對(duì)于運(yùn)算效率就不用多提了

以上就是FFT算法的理論內(nèi)容了，接下來就是用C語言對(duì)這個(gè)算法的實(shí)現(xiàn)了，對(duì)于FFT算法C語言的實(shí)現(xiàn)，網(wǎng)上的方法層出不窮，介于本人比較懶（懶得看別人的程序），再加上自給自足豐衣足食的原則，我自己也寫了一個(gè)個(gè)人認(rèn)為比較通俗易懂的程序，并且為了幫助讀者理解，我特意盡量減少了庫函數(shù)的使用，一些基本的函數(shù)都是自己寫的（難免有很多BUG），但是作為FFT算法已經(jīng)夠用了，目前這個(gè)程序只能處理2^N的數(shù)據(jù)，理論上來講如果不夠2^N的話可以在后面數(shù)列補(bǔ)0這種操作為了簡(jiǎn)約我也就沒有實(shí)現(xiàn)，但是主要的功能是具備的，下面是代碼：

/*

時(shí)間：2018年11月24日

功能：將input里的數(shù)據(jù)進(jìn)行快速傅里葉變換 并且輸出

*/

#include

#include

#define PI 3.1415926

#define FFT_LENGTH 8

double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};

struct complex1{ //定義一個(gè)復(fù)數(shù)結(jié)構(gòu)體

double real; //實(shí)部

double image; //虛部

};

//將input的實(shí)數(shù)結(jié)果存放為復(fù)數(shù)

struct complex1 result_dat[8];

/* 虛數(shù)的乘法

*/

struct complex1 con_complex(struct complex1 a,struct complex1 b){

struct complex1 temp;

temp.real=(a.real*b.real)-(a.image*b.image);

temp.image=(a.image*b.real)+(a.real*b.image);

return temp;

}

/*

簡(jiǎn)單的a的b次方

*/

int mypow(int a,int b){

int i,sum=a;

if(b==0)return 1;

for(i=1;i

sum*=a;

}

return sum;

}

/*

簡(jiǎn)單的求以2為底的正整數(shù)

*/

int log2(int n){

unsigned i=1;

int sum=1;

for(i;;i++){

sum*=2;

if(sum>=n)break;

}

return i;

}

/*

簡(jiǎn)單的交換數(shù)據(jù)的函數(shù)

*/

void swap(struct complex1 *a,struct complex1 *b){

struct complex1 temp;

temp=*a;

*a=*b;

*b=temp;

}

/*

dat為輸入數(shù)據(jù)的數(shù)組

N為抽樣次數(shù) 也代表周期 必須是2^N次方

*/

void fft(struct complex1 dat[],unsigned char N){

/*最終 dat_buf計(jì)算出 當(dāng)前蝶形運(yùn)算奇數(shù)項(xiàng)與W 乘積

dat_org存放上一個(gè)偶數(shù)項(xiàng)的值

*/

struct complex1 dat_buf,dat_org;

/* L為幾級(jí)蝶形運(yùn)算 也代表了2進(jìn)制的位數(shù)

n為當(dāng)前級(jí)蝶形的需要次數(shù) n最初為N/2 每級(jí)蝶形運(yùn)算后都要/2

i j為倒位時(shí)要用到的自增符號(hào) 同時(shí) i也用到了L碟級(jí)數(shù) j是計(jì)算當(dāng)前碟級(jí)的計(jì)算次數(shù)

re_i i_copy均是倒位時(shí)用到的變量

k為當(dāng)前碟級(jí) cos(2*pi/N*k)的 k 也是e^(-j2*pi/N)*k 的 k

*/

unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;

//經(jīng)過觀察，發(fā)現(xiàn)每級(jí)蝶形運(yùn)算需要N/2次運(yùn)算，共運(yùn)算N/2*log2N 次

unsigned char fft_counter=0;

//在此要進(jìn)行補(bǔ)2 N必須是2^n 在此略

//蝶形級(jí)數(shù) (L級(jí))

L=log2(N);

//計(jì)算每級(jí)蝶形計(jì)算的次數(shù)(這里只是一個(gè)初始值) 之后每次要/2

//n=N/2;

//對(duì)dat的順序進(jìn)行倒位

for(i=1;i

i_copy=i;

re_i=0;

for(j=L-1;j>0;j--){

//判斷i的副本最低位的數(shù)字 并且移動(dòng)到最高位 次高位 ..

//re_i為交換的數(shù) 每次它的數(shù)字是不能移動(dòng)的 并且循環(huán)之后要清0

re_i|=((i_copy&0x01)<

i_copy>>=1;

}

swap(&dat[i],&dat[re_i]);

}

//進(jìn)行fft計(jì)算

for(i=0;i

fft_flag=1;

fft_counter=0;

for(j=0;j

if(fft_counter==mypow(2,i)){ //控制隔幾次，運(yùn)算幾次,

fft_flag=0;

}else if(fft_counter==0){ //休止結(jié)束，繼續(xù)運(yùn)算

fft_flag=1;

}

//當(dāng)不判斷這個(gè)語句的時(shí)候 fft_flag保持 這樣就可以持續(xù)運(yùn)算了

if(fft_flag){

dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));

dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));

dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);

//計(jì)算 當(dāng)前蝶形運(yùn)算奇數(shù)項(xiàng)與W 乘積

dat_org.real=dat[j].real;

dat_org.image=dat[j].image; //暫存

dat[j].real=dat_org.real+dat_buf.real;

dat[j].image=dat_org.image+dat_buf.image;

//實(shí)部加實(shí)部 虛部加虛部

dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;

dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;

//實(shí)部減實(shí)部 虛部減虛部

k++;

fft_counter++;

}else{

fft_counter--; //運(yùn)算幾次，就休止幾次

k=0;

}

}

}

}

void main(){

int i;

//先將輸入信號(hào)轉(zhuǎn)換成復(fù)數(shù)

for(i=0;i

result_dat[i].image=0;

//輸入信號(hào)是二維的，暫時(shí)不存在復(fù)數(shù)

result_dat[i].real=input[i];

//result_dat[i].real=10;

//輸入信號(hào)都為實(shí)數(shù)

}

fft(result_dat,FFT_LENGTH);

for(i=0;i

input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);

//取模

printf("%lf\n",input[i]);

}

while(1);

}

這個(gè)程序中input這個(gè)數(shù)組是輸入信號(hào)，在這里只模擬抽樣了8次，輸出的數(shù)據(jù)也是input，如果想看其它序列的話，可以改變FFT_LENGTH的值以及input里的內(nèi)容，程序輸出的是實(shí)部和虛部的模，如果單純想看實(shí)部或者虛部的幅度的話，請(qǐng)自行修改程序~這就是本次淺談FFT以及FFT算法的基本實(shí)現(xiàn)的全部?jī)?nèi)容了參考書籍：數(shù)字信號(hào)處理