快速傅里叶变换（蝶形算法c++源代码）

 详细理论可以自行百度补充

 图1

图2

为了方便处理数据，程序中会将图1中x[0],x[4],x[2]...x[7]进行倒序，也就是图2中的表。这里就直接上代码的。

    for(j = 0; j< N; j++)
    {
       nInter = 0;
        for(i = 0; i<k;i++)
        {
            if(j&(1<<i))//判断第i位是否为1 &两个结果为真的的才是真的 1&1 = 0 1%0 = 0
            {
                nInter += 1<<(k-i-1); //倒过去
            }
        }
        real[j] = tempreal[nInter];
        imag[j] = tempimag[nInter];
    }

从图1中，我们可以看得出，FFT蝶形算法会使用三重循环，下一层的数据都是由上一层计算得到的。这里直接在代码里面备注一下，比较好解释的。

    for(i = 0; i < k; i++)//第一重循环k=log2N，在这里N = 8，所以k = 2
    {
        step = 1<<(i + 1);
        factor_step = N>>(i + 1);//旋转因数变化速度
        //初始化旋转因子
        factor_real = 1.0;
        factor_imag = 0.0;

        for(j = 0; j < step/2 ; j++)
        {
            for(k1=j;k1<N;k1+=step)
            {
                k2 = k1+step/2;//蝶形运算的两个输入
/*
                temp_real = real[k1] + real[k2]*factor_real-imag[k2]*factor_imag;
                real[k2] = real[k1] - (real[k2]*factor_real-imag[k2]*factor_imag);
                real[k1]= temp_real;

                temp_imag = imag[k1] + real[k2]*factor_imag+imag[k2]*factor_real;
                imag[k2] = imag[k1] - (real[k2]*factor_imag+imag[k2]*factor_real);
                imag[k1] = temp_imag;
               */
               temp_real = real[k2]*factor_real-imag[k2]*factor_imag;
               temp_imag = real[k2]*factor_imag+imag[k2]*factor_real;
               real[k2] = real[k1]-temp_real;
               imag[k2] = imag[k1]-temp_imag;
               real[k1] = real[k1]+temp_real;
               imag[k1] = imag[k1]+temp_imag;
            }
            factor_real = inv*cos(-2*PI*(j+1)*factor_step/N);
            factor_imag = inv*sin(-2*PI*(j+1)*factor_step/N);
        }
    }
    if(inv ==-1)
    {
        for(i = 0;i<=N-1;i++)
        {
            real[i]=real[i]/N;
            imag[i]=imag[i]/N;
        }
    }

完整的代码可以见这里：

#include <iostream>
#include<math.h>
#define PI 3.1415926
using namespace std;
void fit(double real[],double imag[],int N,int k,int inv)
{
    int i,j,k1,k2,m,step,factor_step;
    double temp_real,temp_imag,factor_real,factor_imag;

    if(inv!=1&&inv!=-1)
    {
        cout<<"FFT transformation require:inv=1"<<endl;
        cout<<"FFT inverse transformation require:inv=-1"<<endl;
    }
    //倒序
    double tempreal[8];
    double tempimag[8];
    for(i = 0; i < N; i++)
    {
        tempreal[i] = real[i];
        tempimag[i] = imag[i];
    }
    int nInter = 0;
    for(j = 0; j< N; j++)
    {
       nInter = 0;
        for(i = 0; i<k;i++)
        {
            if(j&(1<<i))//判断第i位是否为1 &两个结果为真的的才是真的 1&1 = 0 1%0 = 0
            {
                nInter += 1<<(k-i-1);
            }
        }
        real[j] = tempreal[nInter];
        imag[j] = tempimag[nInter];
    }
    //蝶形算法
    for(i = 0; i < k; i++)
    {
        step = 1<<(i + 1);
        factor_step = N>>(i + 1);//旋转因数变化速度
        //初始化旋转因子
        factor_real = 1.0;
        factor_imag = 0.0;

        for(j = 0; j < step/2 ; j++)
        {
            for(k1=j;k1<N;k1+=step)
            {
                k2 = k1+step/2;//蝶形运算的两个输入
               temp_real = real[k2]*factor_real-imag[k2]*factor_imag;
               temp_imag = real[k2]*factor_imag+imag[k2]*factor_real;
               real[k2] = real[k1]-temp_real;
               imag[k2] = imag[k1]-temp_imag;
               real[k1] = real[k1]+temp_real;
               imag[k1] = imag[k1]+temp_imag;
            }
            factor_real = inv*cos(-2*PI*(j+1)*factor_step/N);
            factor_imag = inv*sin(-2*PI*(j+1)*factor_step/N);
        }
    }
    if(inv ==-1)
    {
        for(i = 0;i<=N-1;i++)
        {
            real[i]=real[i]/N;
            imag[i]=imag[i]/N;
        }
    }

}
int main(int argc, char *argv[])
{
    cout << "Hello World!" << endl;
    double real[8] = {1,2,3,4,5,6,7,8};
    double imag[8] = {0,0,0,0,0,0,0,0};
    double realTwo1[8][8] = {{1,2,3,4,5,6,7,8},
                        {1,2,3,4,5,6,7,8},
                        {1,2,3,4,5,6,7,8},
                        {1,2,3,4,5,6,7,8},
                        {1,2,3,4,5,6,7,8},
                        {1,2,3,4,5,6,7,8},
                        {1,2,3,4,5,6,7,8},
                        {1,2,3,4,5,6,7,8}};
    double realTwo[8][8]={{1,2,4,1,2,4,2,5},
              {1,2,4,1,4,4,2,0},
              {8,2,4,1,8,4,2,5},
              {1,2,4,1,2,5,2,5},
              {1,8,4,1,2,4,2,5},
              {1,2,8,0,2,4,2,5},
              {1,9,4,1,2,4,8,5},
              {1,9,4,1,2,4,8,5}};
    double imagTwo[8][8] = {{0,0,0,0,0,0,0,0},
                        {0,0,0,0,0,0,0,0},
                        {0,0,0,0,0,0,0,0},
                        {0,0,0,0,0,0,0,0},
                        {0,0,0,0,0,0,0,0},
                        {0,0,0,0,0,0,0,0},
                        {0,0,0,0,0,0,0,0},
                        {0,0,0,0,0,0,0,0},};
    fit(real,imag,8,3,1);
    for (int i = 0; i < 8; i++)
    {
        cout<<real[i]<<" "<<imag[i]<<endl;
    }
    //cout<<imag[1];
    for (int i = 0; i < 8; i++)
    {
        fit(realTwo[i],imagTwo[i],8,3,1);
    }
    //转置数组
    for (int j = 0; j < 8 ; j++)
    {
        for (int i = j; i < 8; i++)
        {
            double temprealTwo = 0.0;
            temprealTwo = realTwo[j][i];
            realTwo[j][i] = realTwo[i][j];
            realTwo[i][j] = temprealTwo;
            double tempimagTwo = 0.0;
            tempimagTwo = imagTwo[j][i];
            imagTwo[j][i] = imagTwo[i][j];
            imagTwo[i][j] = tempimagTwo;
        }
    }
    for (int i = 0; i < 8; i++)
    {
        fit(realTwo[i],imagTwo[i],8,3,1);
    }
    //转置数组
    for (int j = 0; j < 8 ; j++)
    {
        for (int i = j; i < 8; i++)
        {
            double temprealTwo = 0.0;
            temprealTwo = realTwo[j][i];
            realTwo[j][i] = realTwo[i][j];
            realTwo[i][j] = temprealTwo;
            double tempimagTwo = 0.0;
            tempimagTwo = imagTwo[j][i];
            imagTwo[j][i] = imagTwo[i][j];
            imagTwo[i][j] = tempimagTwo;
        }
    }
    for (int j = 0; j < 8 ; j++)
    {
        for (int i = j; i < 8; i++)
        {
        cout<<realTwo[j][i]<<endl;
        }
    }
    return 0;
}