FFT算法實現(fft算法)--快速傅里葉變換算法實現

本文作者：韓申權
作者博客：http://www.cnblogs.com/hsqdboke
轉載請注明出處，侵權必究，保留最終解釋權！

首先先要理解離散傅里葉變換(DAT)，然后再理解其快速計算方法(FFT)的原理，和蝴蝶算法的內涵，否則將寫不出代碼；

蝴蝶算法內涵：

(WN)^n的計算：

本代碼為頻率基2抽取算法，也許只看上面這些內容也不一定能懂得FFT算法，不過多查些資料就可以了！

//快速傅里葉變換(FFT) 只看課本是看不會的,介紹的太籠統...
#include<iostream>
#include<math.h>
#include<string.h>
#define PI 3.14159
using namespace std;
double x1[100000],x2[100000],w1[100000],w2[100000];//角標1代表實數部分,2虛數部分
int visited[100000];
int s[1000];    
void wnp(int M,int L,int N) //求(WN)^p
{
    int i,k;
    int t1=(int)(pow(2.0,M-L))-1;
    double t2=-2*PI/N;
    double a,b;
    memset(w1,0,sizeof(w1));
    memset(w2,0,sizeof(w2));
    k=(int)(pow(2.0,L-1));
    a=cos(t2*k);
    b=sin(t2*k);
    w1[0]=1; //當0的情況
    w2[0]=0;
    for(i=1;i<=t1;i++)
    {
        w1[i]=a*w1[i-1]-b*w2[i-1];
        w2[i]=a*w2[i-1]+b*w1[i-1];
    }
}
void FFT(int N,int M)
{
    int i,j,n=N;
    double a,b;
    for(i=1;i<=M;i++)  //從第1級到M級
    {
        n/=2;
        memset(visited,0,sizeof(visited));//標記數組清零
        wnp(M,i,N);  //調用wnp函數    
        for(j=0;j<N;j++)   //分別求每一級的當前實數部分和虛數部分
        {
            if(!visited[j])
            {
                visited[j]=1;
                visited[j+n]=1;
                a=x1[j]-x1[j+n]; //此處曾出錯 應先計算出其值 因為后面 x1[j]和x2[j]的值會改變
                b=x2[j]-x2[j+n];
                x1[j]=x1[j]+x1[j+n];
                x2[j]=x2[j]+x2[j+n];                
                
                int t=j%(N/(int)(pow(2.0,i*1.0)));
                x1[j+n]=a*w1[t]-b*w2[t];
                x2[j+n]=a*w2[t]+b*w1[t];
            }
        }    
    }
}
void solve(double *x,int N,int M) //數位倒讀
{
    int a,i,j,k;
    double t;
    for(k=0;k<N/2;k++)
    {
        i=k;
        a=0;
        memset(s,0,sizeof(s));    
        for(j=0;j<M;j++)
        {
            s[j]=i%2;
            i/=2;
        }
        for(j=0;j<M;j++)
        {
            a=a+s[j]*(int)(pow(2.0,M-1-j));
        }
        t=x[a];
        x[a]=x[k];
        x[k]=t;
    }    
}
int main()
{
    //freopen("d:\\1.txt","r",stdin);
    int N,i,M;
    cout<<"請輸入區段長度N(N需是2的整數次方): ";
    cin>>N;
    M=floor(log10(N*1.0)/log10(2.0)+0.5);
    cout<<"請分別輸入N個采樣值序列復數的實部和虛部: "<<endl;
    for(i=0;i<N;i++)
    {
        printf("實部x1[%d]=",i);
        cin>>x1[i];
        printf("虛部x2[%d]=",i);
        cin>>x2[i];
    }
    FFT(N,M);
    solve(x1,N,M);
    solve(x2,N,M);
    printf("得到的頻譜值為:\n");
    for(i=0;i<N;i++)
        printf("X[%d]=(%.2lf)+(%.2lf)j\n",i,x1[i],x2[i]);
    return 0;
}

本文作者：韓申權
作者博客：http://www.cnblogs.com/hsqdboke
轉載請注明出處，侵權必究，保留最終解釋權！