用c语言实现的FFT

1、一、对FFT的介绍1. FFT （ Fast Fourier Transformation ），即为快速傅里叶变换，是离散 傅里叶变换的快速算法，它是根据离散傅里叶变换的奇、偶、虚、实等 特性，对离散傅里叶变换的算法进行改进获得的。2. FFT算法的基本原理FFT算法是把长序列的DFT逐次分解为较短序列的 DFT。按照抽取方式的不同可分为 DIT-FFT（按时间抽取）和DIF-FFT（按 频率抽取）算法。按蝶形运算的构成不同可分为基 2，基4，基8，以及 任意因子的类型。3. 迭代关系冷 -1e2 r ? ! NfjJ =0y / 2 -1N J 2 - 1二 Ze2刃比(2 J)/Nfl j

2、 +y長此（2十：八0J = o.V 2 - 1N2-1二 2e(2 j ) yf 2 j +W KY总j = Q/ = 0二 F /+W疋 jza4、本次程序的基本过程我们这次所研究的是数字信号处理中的FFT算法，我们这次所用的数字信号是复数类型的。(1) 所以首先，我们先定义了一个复数结构体，因为是进行复数的运算，我们又相继定义复数的加减乘运算的函数。(2) 紧接着，我们定义了进行 FFT计算的fft()快速傅里叶变换函数ini tW()初始化变换核函数即旋转因子的计算，cha nge()变址函数，output()输出傅里叶变换的结果的函数。(3) 定义主函数，并调用定义好的相关子函数，利

3、用fft()中的蝶形运算以及change()函数来完成从时间域上选取的DIT-FFT。二、FFT中码位倒置排序1、码位倒置的实现方法：(1) 简单的利用按位与、或循环实现(2) 利用公式推导的迭代方法2、为什么要进行码位倒置因为由于FFT的计算特性，如果按照正常顺序输入，而没有进行码位倒置的话，就会以乱序输出，就不便于我们后续对信号的相关性质进行研究了，所以DIT-FFT算法就是在进行FFT计算之前，进行分奇偶后的码位倒置运算，即二进制数的倒位。3、倒位序由奇偶分组造成，以 N=8为例，说明如下：Op心乐15卩三、蝶形运算由心（力）：心（疋除示x （k助运算是一种特殊运算1x (疋)=兀 1(

4、£ )+ W x2(k )前一半x (Ar ) = Xj (k )- fF x2(k )后一半(k =(k = 1按照上述公式的规律进行逐级分解，直到2点DFT，如下是N=8时的蝶形算法分析图：四、FFT算法中蝶形算法的基本思想分析（1）我们知道N点FFT运算可以分成Iog2 （ N ）级，每一级都有N/2 个碟形，FFT的基本思想是用3层循环完成全部运算（N点FFT）。（2）第一层循环：由于 N=2m 需要m级计算，第一层循环对运算的 级数进行控制。（stages）（3）第二层循环：由于第L级有2YL-1）个蝶形因子（乘数），第二层循环根据乘数进行控制，保证对于每一个蝶形因子第三层

5、循环要执行 一次，这样，第三层循环在第二层循环控制下，每一级要进行2A（L-1）次循环计算.（选择W）（4）第三层循环：由于第L级共有N/2AL即2A （n-L）个群，并且同一级内不同群的乘数分布相同，当第二层循环确定某一乘数后，第三层循环要将本级中每个群中具有这一乘数的蝶形计算一次，即第三层循环每执行完一次要进行 N/2AL个碟形计算。（执行不同group中具有 相同W的蝶形运算）(5 )可以得出结论：在每一级中，第三层循环完成N/2U 个碟形计算;第二层循环使第三层循环进行2YL-1)次，因此，第二层循环完成时，个碟形计算。实质是：第二、第三层循共进行 2A(L-1) *N/2AL=N/2

6、环完成了第L级的计算。五、用c语言实现的FFT算法如下:<span style="font-size:18px;">#include <stdio.h>#in clude <math.h>#i nclude <stdlib.h>#defi ne N 1000/*定义复数类型*/typedef structdouble real;double img;complex;complex xN, *W; /*输入序列，变换核*/int size_x=0; /*输入序列的大小，在本程序中仅限2的次幕*/double PI;/*圆周率*/

7、void fft(); /*快速傅里叶变换*/void ini tW(); /*初始化变换核*/void cha nge(); /*变址*/void add(complex ,complex ,complex *); /*复数加法*/void mul(complex ,complex ,complex *); /*复数乘法*/void sub(complex ,complex ,complex *); /*复数减法*/void output();/*输出快速傅里叶变换的结果*/int mai n()int i;system("cls");PI=ata n*4; printf

8、("果n");/*输出结果*/输岀DIT方法实现的FFT结输入序列的大小输入序列的实部和虚部prin tf("Please in put the size of x: n");/ sea nf("%d", &size_x);printf("Please in put the data in xN:n");/ for(i=0;i<size_x;i+)printf(" 请输入第%d个序列：",i); sea nf("%lf%lf", &xi.real, &a

9、mp;xi.img);printf(“输岀倒序后的序列n");ini tW();/ 调用变换核fft();/调用快速傅里叶变换printf("输岀FFT后的结果n");output();/调用输岀傅里叶变换结果函数return 0;/*快速傅里叶变换*/void fft()int i=O,j=O,k=O,l=O;complex up,dow n,product;eha nge(); / 调用变址函数 for(i=0;i< log(size_x)/log(2) ;i+)/*l=1<<i;for(j=0;j<size_x;j+= 2*l )/*

10、p的蝶形因子乘数不同*/for(k=0;k<l;k+)/*蝶形运算*/一级蝶形运算stage */一组蝶形运算group,每个grou一个蝶形运算每个group内的mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product，&up);sub(xj+k,product，& dow n);xj+k=up; xj+k+l=dow n;/*初始化变换核，定义一个变换核，相当于旋转因子WAP*/void ini tW() int i;生成变换核用欧拉公式计算旋转因子W=(complex *)malloc(sizeof(compl

11、ex) * size_x); / for(i=0;i<size_x;i+)Wi.real=cos(2*PI/size_x*i); / Wi.img=-1*si n(2*PI/size_x*i);/*变址计算，将 x(n)码位倒置*/void cha nge()complex temp;un sig ned short i=O,j=O,k=O;double t;for(i=0;i<size_x;i+)k=i;j=0;t=(log(size_x)/log (2);while( (t-)>0 ) /利用按位与以及循环实现码位颠倒j=j<<1;j|=(k & 1)

12、;k=k>>1;if(j>i) /将x(n)的码位互换temp=xi;xi=xjxj=temp;output();/*输出傅里叶变换的结果*/void output()int i;printf("The result are as follows: n");for(i=0;i<size_x;i+)prin tf("%.4f',xi.real);if(xi.img>=0.0001)pri ntf("+%.4fjn",xi.img);else if(fabs(xi.img)<0.0001)pr in tf

13、("n"); else prin tf("%.4fjn",xi.img);void add(complex a,complex b,complex *c) / c->real=a.real+b.real; c->img=a.img+b.img;void mul(complex a,complex b,complex *c) / c->real=a.real*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real;复数加法的定义复数乘法的定义void sub(complex a,complex b,complex *c) /复数减法的定义c->real=a.real-b.real; c->img=a.img-b.img;</spa n>六、FFT原理的理解和程序设计中遇到的相关问题及解决方法1、遇到的相关问题:(1)首先一开始不知道什么是 FFT,以及FFT原理是什么不理解FFT中迭代关系的推导以及缘由(3)不理解变址计算的原理(4) 对蝶形运算的推导与原理的理解不透彻(5) 编程过程中对变址计算即对按位与的变换形式的不理解2、解决的方法：(1) 到图书馆借相关的书籍理解相关的原理和过程(2) 到百度收索相关的资料促进理解相