C语言编写FFT程序

DSP课程作业用C语言编写FFT程序1，快速傅里叶变换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的优越性。2，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其中再利用周期性求X(k)的后半部分：n为奇数n为偶数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）、算法特点原位计算1111()()()()()()rmmmNrmmmNXkXkXjWXjXkXjW2102()()xnnnnn蝶形运算两节点的第一个节点为k值，表示成L位二进制数，左移L–m位，把右边空出的位置补零，结果为r的二进制数。倒位序对N=2L点FFT，输入倒位序，输出自然序，第m级运算每个蝶形的两节点距离为2m–1第m级运算：蝶形运算输入序列x(n):N个存储单元系数：N/2个存储单元3，快速傅立叶变换的C语言实现方法我们要衡量一个处理器有没有足够的能力来运行FFT算法，根据以上的简单介绍可以得出以下两点：1.处理器要在一个指令周期能完成乘和累加的工作，因为复数运算要多次查表相乘才能实现。2.间接寻址，可以实现增/减1个变址量，方便各种查表方法。FFT要对原始序列进行反序排列，处理器要有反序间接寻址的能力。#includestdio.h#includemath.h#includestdlib.h#defineN1000typedefstruct{doublereal;doubleimg;}complex;1111111()()(2)(2)()(2)mrmmmNmmrmmmNXkXkXkWXkXkXkWrNW的确定rNW蝶形运算两节点的第一个节点为k值，表示成L位二进制数，左移L–m位，把右边空出的位置补零，结果为r的二进制数。存储单元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,complexb,complex*c){c-real=a.real-b.real;c-img=a.img-b.img;}voiddivi(complexa,complexb,complex*c){c-real=(a.real*b.real+a.img*b.img)/(b.real*b.real+b.img*b.img);c-img=(a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);}