快速傅里叶变换FFT的C语言实现及应用

快速傅里叶变换FFT的C语言实现及应用孙娟0909083014，蔡云0909083001，周硕举09090828301.中南大学升华公寓27栋,湖南410032E-mail:[e-mail]摘要:快速傅里叶变换简介，FFT算法的基本原理，快速傅里叶变换的C语言实现方法，快速傅里叶变换的发展前景，快速傅里叶变换的应用领域，感想。关键词：快速傅里叶变换C语言FFT应用快速傅里叶变换简介计算离散傅里叶变换的一种快速算法,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化快速傅里叶变换成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT).1965年，Cooley和Tukey提出了计算离散傅里叶变换（DFT）的快速算法，将DFT的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基２DIT和基２DIF。FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设快速傅里叶变换x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实快速傅里叶变换数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。FFT算法的基本原理FFT算法的基本思想：利用DFT系数的特性，合并DFT运算中的某些项，吧长序列的DFT—短序列的DFT，从而减少其运算量。FFT算法分类:时间抽选法DIT:Decimation-In-Time；频率抽选法DIF:Decimation-In-Frequency按时间抽选的基-2FFT算法1、算法原理设序列点数N=2L，L为整数。若不满足，则补零。N为2的整数幂的FFT算法称基-2FFT算法。将序列x(n)按n的奇偶分成两组：则x(n)的DFT:12221xrxrxrxr0,1,...,12Nr111000NNNnknknkNNNnnnXkxnWxnWxnWn为奇数n为偶数其中再利用周期性求X(k)的后半部分：221121200221NNrkrkNNrrxrWxrW2211221200NNrkrkkNNNrrxrWWxrW22111/22/200NNrkkrkNNNrrxrWWxrW(,0,1,...1)2Nrk12kNXkWXk2222111100()()(2)NNNNrkrkrrXkxrWxrW2222111100()()(2)NNNNrkrkrrXkxrWxrW(0,1,...1)2Nk121122,222NXkXkNNXkXkXkXk是以为周期的1212()()()()()()2kNkNXkXkWXkNXkXkWXk22NNkkkNNNN又2、运算量当N=2L时，共有L级蝶形，每级N/2个蝶形，每个蝶形有1次复数乘法2次复数加法。2、运算量222()2()loglog2FFmDFTNNNmFFTNN比较DFT2logFaNLNN复数加法：2log22FNNmLN复数乘法：3、算法特点1）原位计算1111()()()()()()rmmmNrmmmNXkXkXjWXjXkXjW2102()()xnnnnn蝶形运算两节点的第一个节点为k值，表示成L位二进制数，左移L–m位，把右边空出的位置补零，结果为r的二进制数。2）倒位序对N=2L点FFT，输入倒位序，输出自然序，第m级运算每个蝶形的两节点距离为2m–1第m级运算：3）蝶形运算输入序列x(n):N个存储单元系数：N/2个存储单元快速傅立叶变换的C语言实现方法有了傅立叶变换，我们可以从信号的频域特征去分析信号。尤其在无线通信系统中，傅里叶变换的重要性就更加明显了，无论是设计者还是测试工程师，在工作中都会和傅立叶变换打交道。我们要衡量一个处理器有没有足够的能力来运行FFT算法，根据以上的简单介绍可以得出以下两点：1.处理器要在一个指令周期能完成乘和累加的工作，因为复数运算要多次查表相乘才能实现。2.间接寻址，可以实现增/减1个变址量，方便各种查表方法。FFT要对原始序列进行反序排列，处理器要有反序间接寻址的能力。下面为一份FFT（快速傅立叶变换）的源码（基于C）/************FFT***********/#includestdio.h#includemath.h#includestdlib.h1111111()()(2)(2)()(2)mrmmmNmmrmmmNXkXkXkWXkXkXkWrNW的确定rNW蝶形运算两节点的第一个节点为k值，表示成L位二进制数，左移L–m位，把右边空出的位置补零，结果为r的二进制数。4）存储单元#defineN1000typedefstruct{doublereal;/*实部*/doubleimg;/*虚部*/}complex;voidfft();/*快速傅里叶变换*/voidifft();/*快速傅里叶逆变换*/voidinitW();/*初始化变化核*/voidchange();/*变址*/voidadd(complex,complex,complex*);/*复数加法*/voidmul(complex,complex,complex*);/*复数乘法*/voidsub(complex,complex,complex*);/*复数减法*/voiddivi(complex,complex,complex*);/*复数除法*/voidoutput();/*输出结果*/complexx[N],*W;/*输出序列的值*/intsize_x=0;/*输入序列的长度，只限2的N次方*/doublePI;intmain(){inti,method;system(cls);PI=atan(1)*4;/*pi等于4乘以1.0的正切值*/printf(Pleaseinputthesizeofx:\n);/*输入序列的长度*/scanf(%d,&size_x);printf(Pleaseinputthedatainx[N]:(suchas:56)\n);/*输入序列对应的值*/for(i=0;isize_x;i++)scanf(%lf%lf,&x[i].real,&x[i].img);initW();/*选择FFT或逆FFT运算*/printf(UseFFT(0)orIFFT(1)?\n);scanf(%d,&method);if(method==0)fft();elseifft();output();return0;}/*进行基-2FFT运算*/voidfft(){inti=0,j=0,k=0,l=0;complexup,down,product;change();for(i=0;ilog(size_x)/log(2);i++)/*一级蝶形运算*/{l=1i;for(j=0;jsize_x;j+=2*l)/*一组蝶形运算*/{for(k=0;kl;k++)/*一个蝶形运算*/{mul(x[j+k+l],W[size_x*k/2/l],&product);add(x[j+k],product,&up);sub(x[j+k],product,&down);x[j+k]=up;x[j+k+l]=down;}}}}voidifft(){inti=0,j=0,k=0,l=size_x;complexup,down;for(i=0;i(int)(log(size_x)/log(2));i++)/*一级蝶形运算*/{l/=2;for(j=0;jsize_x;j+=2*l)/*一组蝶形运算*/{for(k=0;kl;k++)/*一个蝶形运算*/{add(x[j+k],x[j+k+l],&up);up.real/=2;up.img/=2;sub(x[j+k],x[j+k+l],&down);down.real/=2;down.img/=2;divi(down,W[size_x*k/2/l],&down);x[j+k]=up;x[j+k+l]=down;}}}change();}/*初始化变化核*/voidinitW(){inti;W=(complex*)malloc(sizeof(complex)*size_x);for(i=0;isize_x;i++){W[i].real=cos(2*PI/size_x*i);W[i].img=-1*sin(2*PI/size_x*i);}}/*变址计算，将x(n)码位倒置*/voidchange(){complextemp;unsignedshorti=0,j=0,k=0;doublet;for(i=0;isize_x;i++){k=i;j=0;t=(log(size_x)/log(2));while((t--)0){j=j1;j|=(k&1);k=k1;}if(ji){temp=x[i];x[i]=x[j];x[j]=temp;}}}voidoutput()/*输出结果*/{inti;printf(Theresultareasfollows\n);for(i=0;isize_x;i++){printf(%.4f,x[i].real);if(x[i].img=0.0001)printf(+%.4fj\n,x[i].img);elseif(fabs(x[i].img)0.0001)printf(\n);elseprintf(%.4fj\n,x[i].img);}}voidadd(complexa,complexb,complex*c){c-real=a.real+b.real;c-img=a.img+b.img;}voidmul(complexa,complexb,complex*c){c-real=a.real*b.real-a.img*b.img;c-img=a.real*b.img+a.img*b.real;}voidsub(complexa,comple