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

相信網上現在有很多關於FFT的教程，我曾經也參閱了很多網上的教程，感覺都不怎麼通俗易懂。在基本上的研究FFT，並且通過編程的形式實現之後。我決定寫一篇通俗易懂的關於FFT的講解。因此我在接下來的敘述中儘量非常通俗細緻的講解。

本人最早知道傅立葉變換的時候是沉迷於音樂的頻譜跳動無法自拔，當時就很想做一個音樂頻譜顯示器。搜閱了很多資料之後，才了解到傅立葉變換，和FFT。當然這都是以前的事情了，經過了系統的學習+2個星期的研究，自製了一個FFT的算法，不敢說速度上的優勢，但是個人認為是一種通俗易懂的實現方法。經過實際的VC++模擬實驗、和STM32跑的也很成功。

首先，要會FFT，就必須要對DFT有所了解，因為兩者之間本質上是一樣的。在此之前，先列出離散傅立葉變換對(DFT):

,k=0,1,…N-1

n=0,1…N-1

其中： 

但是FFT之所以稱之為快速傅立葉變換，就利用了以下的幾個性質（重中之重！）

周期性： 

對稱性： 

可約性： 

先把這仨公式放到這，接下來會用到。

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

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

對於 ,k=0,1,…N-1   設N是2的整次冪 也就是N=2^M

將x(n)按照奇偶分組

x(2r)=x1(r)

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

那麼將DFT也分為兩組來預算 

  (第一項是偶  第二項是奇)

這個時候，我們上邊提的三個性質中的可約性就起到作用了：

回顧一下： 

那麼這個式子就可以化為：

(式1-1)

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

其中：

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

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

(式1-2)

那麼式（1-1）和（1-2）就可以專用一個蝶形信號流圖符號來表示。如圖：

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

通過這樣的分解，每一個N/2點DFT只需(N^2)/4次複數相乘計算，明顯的節省了運算量。

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

對於這張圖片我想強調的一點就是第二階蝶形運算的時候，WN0之後不應該是WN1嗎，為什麼是W2N了，這個問題之前困擾了我一段時間，但是不要忘了，前者的    W0N的展開是W0N/2  因為此時N已經按照奇偶分開了，所以是N/2  而W2N/2也就是  W2N是根據可約性得出的，在這裡不能混淆.。

對於運算效率就不用多提了

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

/*

時間：2018年11月24日

功能：將input裡的數據進行快速傅立葉變換  並且輸出

*/

#include

#include

#define PI 3.1415926

#define FFT_LENGTH                8

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

struct complex1{                //定義一個複數結構體

double real;        //實部

double image;        //虛部

};        

//將input的實數結果存放為複數

struct complex1 result_dat[8];

/*        虛數的乘法

*/

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;

}

/*

簡單的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;

}

/*

簡單的求以2為底的正整數

*/

int log2(int n){

unsigned i=1;

int sum=1;

for(i;;i++){

sum*=2;

if(sum>=n)break;

}

return i;

}

/*

簡單的交換數據的函數

*/

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

struct complex1 temp;

temp=*a;

*a=*b;

*b=temp;

}

/*

dat為輸入數據的數組

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

*/

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

/*最終  dat_buf計算出 當前蝶形運算奇數項與W  乘積

dat_org存放上一個偶數項的值

*/

struct complex1 dat_buf,dat_org;

/*        L為幾級蝶形運算    也代表了2進位的位數

n為當前級蝶形的需要次數  n最初為N/2 每級蝶形運算後都要/2

i j為倒位時要用到的自增符號  同時  i也用到了L碟級數   j是計算當前碟級的計算次數 

re_i i_copy均是倒位時用到的變量

k為當前碟級  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;

//經過觀察，發現每級蝶形運算需要N/2次運算，共運算N/2*log2N  次  

unsigned char fft_counter=0;

//在此要進行補2   N必須是2^n   在此略

//蝶形級數  (L級)

L=log2(N);        

//計算每級蝶形計算的次數(這裡只是一個初始值)  之後每次要/2

//n=N/2;

//對dat的順序進行倒位

for(i=1;i

i_copy=i;

re_i=0;

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

//判斷i的副本最低位的數字  並且移動到最高位  次高位  ..

//re_i為交換的數   每次它的數字是不能移動的 並且循環之後要清0

re_i|=((i_copy&0x01)<

i_copy>>=1;

}

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

}

//進行fft計算

for(i=0;i

fft_flag=1;

fft_counter=0;

for(j=0;j

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

fft_flag=0;

}else if(fft_counter==0){                //休止結束，繼續運算

fft_flag=1;

}

//當不判斷這個語句的時候  fft_flag保持  這樣就可以持續運算了

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);

//計算 當前蝶形運算奇數項與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;                

//實部加實部   虛部加虛部

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

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

//實部減實部        虛部減虛部

k++;

fft_counter++;

}else{

fft_counter--;                                //運算幾次，就休止幾次

k=0;

}

}

}

}

void main(){

int i;

//先將輸入信號轉換成複數

for(i=0;i

result_dat[i].image=0;        

//輸入信號是二維的，暫時不存在複數

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

//result_dat[i].real=10;

//輸入信號都為實數

}

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);

}

這個程序中input這個數組是輸入信號，在這裡只模擬抽樣了8次，輸出的數據也是input，如果想看其它序列的話，可以改變FFT_LENGTH的值以及 input裡的內容，程序輸出的是實部和虛部的模，如果單純想看實部或者虛部的幅度的話，請自行修改程序~這就是本次淺談FFT以及FFT算法的基本實現的全部內容了參考書籍：數位訊號處理 

打開APP閱讀更多精彩內容

聲明：本文內容及配圖由入駐作者撰寫或者入駐合作網站授權轉載。文章觀點僅代表作者本人，不代表電子發燒友網立場。文章及其配圖僅供工程師學習之用，如有內容圖片侵權或者其他問題，請聯繫本站作侵刪。 侵權投訴