FFT算法理解与c语言的实现

来源：互联网发布：淘宝卖写真集违规吗编辑：程序博客网时间：2024/05/17 23:18

学习DIP第5天
为什么需要FFT

第一个问题是为什么要创造FFT，简单的说，为了速度。我们承认DFT很有用，但是我们发现他的速度不是很快，1D的DFT原始算法的时间复杂度是O(n^2)，这个可以通过公式观察出来，对于2D的DFT其时间复杂度是O（n^4），这个速度真的很难接受，也就是说，你计算一幅1024x768的图像时，你将等大概。。。大概。。。我也没试过，反正很久。不信的自己去试试。所以找到一种快速方法的方法计算FFT势在必行。
以下为DFT公式
计算一个4点DFT。计算量如下：

如何得到FFT

通过观察DFT公式，我们发现DFT计算每一项X[u]的时候都遍历了完整的x[n]所以，我们的想法就是能不能通过其他的X[u+m]（m为一个整数，可正可负）得到X[u]，也就是充分利用之间的计算结构来构建现在的结果，这种方法就很容易表现成迭代的形式。本文介绍基2的FFT，及离散信号x[n]的个数为2的k次方，即如果完整的离散信号中有N个信分量{x1，x2，x3....xN}其中N=1<<k。

数学基础

FFT的数学基础其实就是：

旋转因子具有以下性质：

这些性质使得我们可以利用前面的计算结果来完成出后续的计算。

2-FFT

观察：一个只有两个值的离散信号，假设N=2，利用性质2-对称性可以得到

4-FFT

与上面2点的一样，推广到4点，将会是这样，其中方框内的操作为上述2-FFT：

最终算法：

8-FFT

废话不多说，和上面一样：

计算过程为：

结果为：

完整的计算过程见下面的动图。CSDN不会动，网盘自己下来看：

4-FFT.gif

2^n-FFT

同理，推广到2^n，可以得到类似的结果，不过我们发现，为了使输出结果为顺序结果，输入的顺序经过了一系列的调整，而且每一步WN的幂次参数也是变化，我们必须得到其变化规律才能更好的编写程序：

观察WN的变化规律为：

节点距离（设为z）就是从WN的0次幂开始连续加到WN的z次幂。然后间隔为z个式子，再次从0次幂加到z次。重复以上过程：

例如第二级，间隔z=2，节点为实心黑色点，节点0,1，不做操作，节点2*W8^0,节点3*W8^2，节点4，5不做操作。。

其内在规律就是节点i是否乘以WN：

if（i%z==奇数）

节点i*WN^(step)；

step每次增加的数量由当前的所在的计算级决定

step=1<<（k（总级数）-i（当前级数））

输入参数重新排序

i行表示数组下标，蓝色数字表示实际数据：其内在规律就是，下标为偶数的将被映射到自己本身下标的1/2处，下表为奇数时被映射到数组后半段（size_n/2）的（下标-1）1/2处，将排列后的数据分为前后两个部分，递归次过程，直到只有两个元素。则停止。过程如下：

=================================================================

==================================================================================================================================

==================================================================================================================================

观察算法

观察至此，我们已经基本把FFT蝶形算法的所有特征都搞定了，接下来就是使用代码来实现了。

实现代码

////  main.c//  Fourer1D////  Created by Tony on 14/11/16.//  Copyright (c) 2014年 Tony. All rights reserved.//#include <stdio.h>#include <stdlib.h>#include <time.h>#include <math.h>//mac下M_PI在math.h中有宏定义，所以这里我们选择行的宏定义#ifndef M_PI#define M_PI 3.14159265358979323846#endif#define SIZE 1024*16#define VALUE_MAX 1000//////////////////////////////////////////////////////////////////////定义一个复数结构体///////////////////////////////////////////////////////////////////struct Complex_{    double real;    double imagin;};typedef struct Complex_ Complex;//////////////////////////////////////////////////////////////////////定义一个复数计算，包括乘法，加法，减法///////////////////////////////////////////////////////////////////void Add_Complex(Complex * src1,Complex *src2,Complex *dst){    dst->imagin=src1->imagin+src2->imagin;    dst->real=src1->real+src2->real;}void Sub_Complex(Complex * src1,Complex *src2,Complex *dst){    dst->imagin=src1->imagin-src2->imagin;    dst->real=src1->real-src2->real;}void Multy_Complex(Complex * src1,Complex *src2,Complex *dst){    double r1=0.0,r2=0.0;    double i1=0.0,i2=0.0;    r1=src1->real;    r2=src2->real;    i1=src1->imagin;    i2=src2->imagin;    dst->imagin=r1*i2+r2*i1;    dst->real=r1*r2-i1*i2;}//////////////////////////////////////////////////////////////////////在FFT中有一个WN的n次方项，在迭代中会不断用到，具体见算法说明///////////////////////////////////////////////////////////////////void getWN(double n,double size_n,Complex * dst){    double x=2.0*M_PI*n/size_n;    dst->imagin=-sin(x);    dst->real=cos(x);}//////////////////////////////////////////////////////////////////////随机生成一个输入，显示数据部分已经注释掉了//注释掉的显示部分为数据显示，可以观察结果///////////////////////////////////////////////////////////////////void setInput(double * data,int  n){    //printf("Setinput signal:\n");    srand((int)time(0));    for(int i=0;i<SIZE;i++){        data[i]=rand()%VALUE_MAX;        //printf("%lf\n",data[i]);    }}//////////////////////////////////////////////////////////////////////定义DFT函数，其原理为简单的DFT定义，时间复杂度O（n^2）,//下面函数中有两层循环，每层循环的step为1，size为n，故为O（n*n）,//注释掉的显示部分为数据显示，可以观察结果///////////////////////////////////////////////////////////////////void DFT(double * src,Complex * dst,int size){    clock_t start,end;    start=clock();        for(int m=0;m<size;m++){        double real=0.0;        double imagin=0.0;        for(int n=0;n<size;n++){            double x=M_PI*2*m*n;            real+=src[n]*cos(x/size);            imagin+=src[n]*(-sin(x/size));                }        dst[m].imagin=imagin;        dst[m].real=real;       /* if(imagin>=0.0)            printf("%lf+%lfj\n",real,imagin);        else            printf("%lf%lfj\n",real,imagin);*/    }    end=clock();    printf("DFT use time :%lf for Datasize of:%d\n",(double)(end-start)/CLOCKS_PER_SEC,size);}//////////////////////////////////////////////////////////////////////定义IDFT函数，其原理为简单的IDFT定义，时间复杂度O（n^2）,//下面函数中有两层循环，每层循环的step为1，size为n，故为O（n*n）,///////////////////////////////////////////////////////////////////void IDFT(Complex *src,Complex *dst,int size){    clock_t start,end;    start=clock();    for(int m=0;m<size;m++){        double real=0.0;        double imagin=0.0;        for(int n=0;n<size;n++){            double x=M_PI*2*m*n/size;            real+=src[n].real*cos(x)-src[n].imagin*sin(x);            imagin+=src[n].real*sin(x)+src[n].imagin*cos(x);                    }        real/=SIZE;        imagin/=SIZE;        if(dst!=NULL){            dst[m].real=real;            dst[m].imagin=imagin;        }        if(imagin>=0.0)            printf("%lf+%lfj\n",real,imagin);        else            printf("%lf%lfj\n",real,imagin);    }    end=clock();    printf("IDFT use time :%lfs for Datasize of:%d\n",(double)(end-start)/CLOCKS_PER_SEC,size);    }//////////////////////////////////////////////////////////////////////定义FFT的初始化数据，因为FFT的数据经过重新映射，递归结构///////////////////////////////////////////////////////////////////int FFT_remap(double * src,int size_n){    if(size_n==1)        return 0;    double * temp=(double *)malloc(sizeof(double)*size_n);    for(int i=0;i<size_n;i++)        if(i%2==0)            temp[i/2]=src[i];        else            temp[(size_n+i)/2]=src[i];    for(int i=0;i<size_n;i++)        src[i]=temp[i];    free(temp);    FFT_remap(src, size_n/2);    FFT_remap(src+size_n/2, size_n/2);    return 1;}//////////////////////////////////////////////////////////////////////定义FFT，具体见算法说明，注释掉的显示部分为数据显示，可以观察结果///////////////////////////////////////////////////////////////////void FFT(double * src,Complex * dst,int size_n){    FFT_remap(src, size_n);   // for(int i=0;i<size_n;i++)    //    printf("%lf\n",src[i]);clock_t start,end;    start=clock();    int k=size_n;    int z=0;    while (k/=2) {        z++;    }    k=z;    if(size_n!=(1<<k))        exit(0);    Complex * src_com=(Complex*)malloc(sizeof(Complex)*size_n);    if(src_com==NULL)        exit(0);    for(int i=0;i<size_n;i++){        src_com[i].real=src[i];        src_com[i].imagin=0;    }    for(int i=0;i<k;i++){        z=0;        for(int j=0;j<size_n;j++){            if((j/(1<<i))%2==1){                Complex wn;                getWN(z, size_n, &wn);                Multy_Complex(&src_com[j], &wn,&src_com[j]);                z+=1<<(k-i-1);                Complex temp;                int neighbour=j-(1<<(i));                temp.real=src_com[neighbour].real;                temp.imagin=src_com[neighbour].imagin;                Add_Complex(&temp, &src_com[j], &src_com[neighbour]);                Sub_Complex(&temp, &src_com[j], &src_com[j]);            }            else                z=0;        }        }       /* for(int i=0;i<size_n;i++)        if(src_com[i].imagin>=0.0){            printf("%lf+%lfj\n",src_com[i].real,src_com[i].imagin);        }        else            printf("%lf%lfj\n",src_com[i].real,src_com[i].imagin);*/for(int i=0;i<size_n;i++){dst[i].imagin=src_com[i].imagin;dst[i].real=src_com[i].real;}end=clock();    printf("FFT use time :%lfs for Datasize of:%d\n",(double)(end-start)/CLOCKS_PER_SEC,size_n);    }////////////////////////////////////////////////////////////////////int main(int argc, const char * argv[]) {    double input[SIZE];    Complex dst[SIZE];    setInput(input,SIZE);    printf("\n\n");    DFT(input, dst, SIZE);    printf("\n\n");    FFT(input, dst, SIZE);    //IDFT(dst, NULL, SIZE);    getchar();}