淺談FFT以及FFT算法的基本實(shí)現

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

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

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

其中：

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

周期性：

對稱(chēng)性：

可約性：

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

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

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

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

將x(n)按照奇偶分組

x(2r)=X₁(r)

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

那么將DFT也分為兩組來(lái)預算

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

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

回顧一下：

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

則X₁和X₂都是長(cháng)度為N/2的序列x1和x₂的N/2點(diǎn)的離散傅里葉變換

                   其中：

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

(式1-2)

那么式（1-1）和（1-2）就可以專(zhuān)用一個(gè)蝶形信號流圖符號來(lái)表示。如圖：

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

那么，通過(guò)這樣的分解，每一個(gè)N/2點(diǎn)DFT只需(N^2)/4次復數相乘計算，明顯的節省了運算量。

那么以此類(lèi)推，繼續將已經(jīng)得出的X_1(K)和X₂(k)按奇偶繼續分解，還是上邊一樣的方法。那么就得出了百度上都可以找到的一大堆的這個(gè)圖片了。

對于這張圖片我想強調的一點(diǎn)就是第二階蝶形運算的時(shí)候，W⁰_N之后不應該是W¹_N嗎，為什么是W²_N了，這個(gè)問(wèn)題之前困擾了我一段時(shí)間，但是不要忘了，前者的W⁰_N的展開(kāi)是 W⁰_N/2因為此時(shí)N已經(jīng)按照奇偶分開(kāi)了，所以是N/2  而W²_N也就是 W¹_N/2是根據可約性得出的，在這里不能混淆.。

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

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

/*
    時(shí)間：2018年11月24日
    功能：將input里的數據進(jìn)行快速傅里葉變換  并且輸出
*/
 
#include <stdio.h>
#include <math.h>
#define PI 3.1415926
#define FFT_LENGTH      8
double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};
struct complex1{        //定義一個(gè)復數結構體
    double real;    //實(shí)部
    double image;   //虛部
};  
//將input的實(shí)數結果存放為復數
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;
}
 
/*
    簡(jiǎn)單的a的b次方
*/
int mypow(int a,int b){
    int i,sum=a;
    if(b==0)return 1;
    for(i=1;i<b;i++){
        sum*=a; 
    }
    return sum;
}
/*
    簡(jiǎn)單的求以2為底的正整數
*/
int log2(int n){
    unsigned i=1;
    int sum=1;
    for(i;;i++){
        sum*=2;
        if(sum>=n)break;
    }
    return i;
}
/*
    簡(jiǎn)單的交換數據的函數
*/
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存放上一個(gè)偶數項的值
    */
    struct complex1 dat_buf,dat_org;
    /*  L為幾級蝶形運算    也代表了2進(jìn)制的位數
        n為當前級蝶形的需要次數  n最初為N/2 每級蝶形運算后都要/2
        i j為倒位時(shí)要用到的自增符號  同時(shí)  i也用到了L碟級數   j是計算當前碟級的計算次數 
        re_i i_copy均是倒位時(shí)用到的變量
        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;
    //經(jīng)過(guò)觀(guān)察，發(fā)現每級蝶形運算需要N/2次運算，共運算N/2*log2N  次  
    unsigned char fft_counter=0;
    //在此要進(jìn)行補2   N必須是2^n   在此略
    //蝶形級數  (L級)
    L=log2(N);  
    //計算每級蝶形計算的次數(這里只是一個(gè)初始值)  之后每次要/2
    //n=N/2;
 
    //對dat的順序進(jìn)行倒位
    for(i=1;i<N/2;i++){
        i_copy=i;
        re_i=0;
        for(j=L-1;j>0;j--){
            //判斷i的副本最低位的數字  并且移動(dòng)到最高位  次高位  ..
            //re_i為交換的數   每次它的數字是不能移動(dòng)的 并且循環(huán)之后要清0
            re_i|=((i_copy&0x01)<<j);       
            i_copy>>=1;
        }
        swap(&dat[i],&dat[re_i]);
    }
    //進(jìn)行fft計算
    for(i=0;i<L;i++){
        
        fft_flag=1;
        fft_counter=0;
        for(j=0;j<N;j++){
            if(fft_counter==mypow(2,i)){        //控制隔幾次，運算幾次,
                fft_flag=0;
            }else if(fft_counter==0){       //休止結束，繼續運算
                fft_flag=1;
            }
            //當不判斷這個(gè)語(yǔ)句的時(shí)候  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;       
                                                    //實(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--;              //運算幾次，就休止幾次
                k=0;
            }
        }
    }
 
    
}
void main(){
    int i;
    //先將輸入信號轉換成復數
    for(i=0;i<FFT_LENGTH;i++){
        result_dat[i].image=0;  
                                //輸入信號是二維的，暫時(shí)不存在復數
        result_dat[i].real=input[i];
        //result_dat[i].real=10;
                                //輸入信號都為實(shí)數
    }
    fft(result_dat,FFT_LENGTH);
    for(i=0;i<FFT_LENGTH;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);
}

這就是本次淺談FFT以及FFT算法的基本實(shí)現的全部?jì)热萘?/span>

原創(chuàng )內容轉載其注明原作者.  

參考書(shū)籍：數字信號處理