C语言编写FFT程序

徐士良老师编写的c语言算法程序下载链接：https://pan.baidu.com/s/1zDV6iLeYeXmZaoZlP4yRAA
提取码：8opo

一、什么是FFT？

FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

二、FFT的作用

FFT可以用来加速多项式乘法。假设有两个n−1次多项式A(x)和B(x)，我们的目标是——把它们乘起来。
普通的多项式乘法的复杂度是O(n2)的，我们要枚举A(x)中的每一项，分别与B(x)中的每一项相乘，来得到一个新的多项式C(x)。
但是，如果A(x)，B(x)两个多项式用点值表示的方法进行相乘，复杂度是O(n)的。具体方法：C(xi)=A(xi)×B(xi)，所以枚举xi即可。
要是我们把两个多项式转换成点值表示，再相乘，再把新的点值表示转换成多项式岂不就可以O(n)的复杂度来解决多项式乘法了！
显然，把多项式转换成点值表示的朴素算法是O(n2) O(n^2)O(n2 )的。难道大整数乘法就只能是O(n2) O(n^2)O(n2)吗？不甘心的同学可以发现，大整数乘法复杂度的瓶颈可能在“多项式转换成点值表示”这一步做改进，只要完成这一步就可以O(n)的复杂度求答案了。

三、n次单位根

四、FFT的核心

五，徐士良老师编写的FFT程序

1.FFT源程序【FFT.c】
抽样点为：64
输入信号：0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)

#include "stdio.h"
#include "math.h"
//#include "kfft.c"
#include <D:\temp\kfft.c>
#include <time.h>#define PI 3.1415926535main()
{ int i,j;double pr[64],pi[64],fr[64],fi[64],t[64];clock_t begin, end;double cost1 ,cost2, persent;for (i=0; i<=63; i++)  //生成输入信号{ t[i] = i*0.001;pr[i]=0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i); pi[i]=0.0; //0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i)}begin = clock();                                  //开始记录时间kfft(pr,pi,64,6,fr,fi);  //调用FFT函数end = clock();                                       //结束记录时间cost1 = (double)(end - begin) / CLOCKS_PER_SEC;for (i=0; i<64; i++){ printf("%d\t%lf\n",i,pr[i]); //输出结果}
}

2.源程序【kfft.c】

  #include "math.h"void kfft(pr,pi,n,k,fr,fi)int n,k;double pr[],pi[],fr[],fi[];{ int it,m,is,i,j,nv,l0;double p,q,s,vr,vi,poddr,poddi;for (it=0; it<=n-1; it++)  //将pr的实部和虚部循环赋值给fr[]和fi[]{ m=it; is=0;for(i=0; i<=k-1; i++){ j=m/2; is=2*is+(m-2*j); m=j;}fr[it]=pr[is]; fi[it]=pi[is];}pr[0]=1.0; pi[0]=0.0;p=6.283185306/(1.0*n);pr[1]=cos(p); //将w=e^-j2pi/n用欧拉公式表示pi[1]=-sin(p);for (i=2; i<=n-1; i++)  //计算pr[]{ p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);pr[i]=p-q; pi[i]=s-p-q;}for (it=0; it<=n-2; it=it+2)  { vr=fr[it]; vi=fi[it];fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];}m=n/2; nv=2;for (l0=k-2; l0>=0; l0--) //蝴蝶操作{ m=m/2; nv=2*nv;for (it=0; it<=(m-1)*nv; it=it+nv)for (j=0; j<=(nv/2)-1; j++){ p=pr[m*j]*fr[it+j+nv/2];q=pi[m*j]*fi[it+j+nv/2];s=pr[m*j]+pi[m*j];s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);poddr=p-q; poddi=s-p-q;fr[it+j+nv/2]=fr[it+j]-poddr;fi[it+j+nv/2]=fi[it+j]-poddi;fr[it+j]=fr[it+j]+poddr;fi[it+j]=fi[it+j]+poddi;}}for (i=0; i<n; i++){ pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);  //幅度的计算}return;}

3、TCC编译

4、图像显示
跟之前DFT画的图像一致

六、自己编写的FFT程序

（不足之处还望大家指出，谢谢亲们了！）
1、程序源码


#include <stdio.h>
#include <math.h>
#include <time.h>#define N 64                                         //设置抽样点数
#define PI 3.1415   //定义圆周率typedef struct                                           //定义复数结构体变量
{double real;double imag;
}complex;void c_jiafa(complex, complex, complex *);             //复数加运算
void c_jianfa(complex, complex, complex *);             //复数减运算
void c_chengfa(complex, complex, complex *);                //复数乘运算
void Wn_i(int, int, complex *);                         //FFT旋转因子
void Wn_ik(int, int, int, complex *);                   //DFT旋转因子int main()
{complex f[N];                          //f[N]为输出的FFT序列int LH, K, J, B, L, k, N1, P, M, K1;             //L表示第L级蝶形  p旋转因子指数 B两序列间隔点数k   第k个序列double T,y[N];double t_r,t_i;M =6;LH = N / 2;J = LH;N1 = N - 1;for (int i = 0; i < N; i++)                           //为FFT运算提供初始序列{f[i].real =0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i) ;f[i].imag = 0;}//printf("*******************************  级数为：%d  *******************************\n", M);//.......................................................................................................................for (int I = 0; I < N1; I++)                      //定义倒序序列函数{if (I < J){t_r=f[I].real;t_i=f[I].imag;f[I].real=f[J].real;f[I].imag=f[J].imag;f[J].real=t_r;f[J].imag=t_i;}K = N>>1;while(K<=J){J = J - K;K>>1;}J=J+K;}//FFT运算for (L = 1; L <= M; L++)                          {B = (int)(pow(2, L - 1));//第L级，每个蝶形的两个数据有B=2^(L-1)个点for (J = 0; J < B; J++){P = (int)(J*pow(2, M - L));//每级有B个旋转因子，一旋转因子对应着2^(M-1)个蝶形for (k = J; k < N; k = (int)(k + pow(2, L))){K1 = k + B;complex wn, t;Wn_i(N, P, &wn);c_chengfa(f[K1], wn, &t);                    //。。。。。。。。。。。。c_jianfa(f[k], t, &(f[K1]));              //蝶形运算c_jiafa(f[k], t, &(f[k]));                //。。。。。。。。。。。。}}}for (int i = 0; i < N; i++)                          //快速傅里叶变换输出{y[i]=sqrt(f[i].real*f[i].real+f[i].imag*f[i].imag);printf("%d %lf\n", i, y[i]);}return 0;
}void c_jiafa(complex a, complex b, complex *c)         //复数加法
{c->real = a.real + b.real;c->imag = a.imag + b.imag;
}
void c_jianfa(complex a, complex b, complex *c)         //复数减法
{c->real = a.real - b.real;c->imag = a.imag - b.imag;
}
void c_chengfa(complex a, complex b, complex *c)            //复数乘法
{c->real = a.real*b.real - a.imag*b.imag;c->imag = a.real*b.imag + a.imag*b.real;
}void Wn_i(int n1, int i, complex *Wn)                  //定义FFT旋转因子
{Wn->real = cos(2 * PI*i / n1);Wn->imag = -sin(2 * PI*i / n1);
}

2、TCC编译

3、图像显示

gnuplot> plot [0:64] [0:12.0] "<7.exe" u 1:2 w l