*基二FFT算法*/#include "math.h"#include "stdio.h"struct compx{ double real; double imag;} compx ;struct compx EE(struct compx b1,struct compx b2)//复数相乘{struct compx b3;b3.real=b1.real*b2.real-b1.imag*b2.imag;b3.imag=b1.real*b2.imag+b1.imag*b2.real;return(b3);}void FFT(struct compx *xin,int N){int f,m,LH,nm,i,k,j,L;double p , ps ;int le,B,ip;float pi;struct compx v,w,t;LH=N/2;f=N;for(m=1;(f=f/2)!=1;m++){;}//求出m为log2 Nnm=N-2; j=N/2;//变址运算，对时间进行奇偶分解for(i=1;i<=nm;i++)//即xin第一位和最后一位不用操作，不用变址，其余各位根据码位倒置{if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}k=LH;while(j>=k){j=j-k;k=k/2;}j=j+k;}{for(L=1;L<=m;L++)//运行m级蝶形运算{ le=pow(2,L);B=le/2;pi=3.14159; for(j=0;j<=B-1;j++) { p=pow(2,m-L)*j;//k的取值（0，1，2，...，(pow(2,L)/2)-1），当在第一级时，p只为0，当在第二//级时，p为0和2，当为第三级时，p为0，1，2，3.在递升的级别中不断延续，见数字信号处理（华中//科技大学p83图的系数WN）,且相邻不同种基本蝶形的蝶形节系数的增量为（2*pi/N）*pow(2,m-L) ps=（2*pi/N）*p; w.real=cos(ps); w.imag=-sin(ps);//求出WNk for(i=j;i<=N-1;i=i+le)//确定与蝶形系数相乘的Xm(q)的下标m，此方法为蝶形图的特点来得到，//相邻同种基本蝶形的间距为2的L次方。 { ip=i+B;//即对频谱进行前后分解 t=EE(xin[ip],w);//复数相乘 xin[ip].real=xin[i].real-t.real;//基本蝶式运算 xin[ip].imag=xin[i].imag-t.imag; xin[i].real=xin[i].real+t.real; xin[i].imag=xin[i].imag+t.imag; } }}} return ;}//输入时域数据点为num,则输出频域数据点同为num,num是数据长度，必须为2的整数次幕，//其大小由数据采样定理来决定#include <math.h>#include <stdio.h>#include <stdlib.h>float result[257];//振幅，其平方为功率谱struct compx s[257];int Num=16;//数据长度，必须为2的整数次幕const float pp=3.14159;main(){int i;for(i=0;i<16;i++){s[i].real=sin(pp*i/32);s[i].imag=0;}FFT(s,Num);for(i=0;i<16;i++){printf("%.4f",s[i].real);printf("+%.4fj\n",s[i].imag);result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));//pow功 能: 指数函数(x的y次方) 用 法: double pow(double x, double y);}} 基二FFT的C语言实现 频域基于时域和基于频域的算法很相似的，只不过时域里是先乘后加减，对输入序列进行倒序，而频域的算法是先加减后乘，输入序列不用倒序，对输出序列进行倒序。调试成功之后，发现两者结果相差很小了。 #include "math.h"#include "stdio.h"struct compx{ double real; double imag;} compx ;struct compx EE(struct compx b1,struct compx b2){struct compx b3;b3.real=b1.real*b2.real-b1.imag*b2.imag;b3.imag=b1.real*b2.imag+b1.imag*b2.real;return(b3);}void FFT(struct compx *xin,int N){int f,m,LH,nm,i,k,j,L;double p , ps ;int le,B,ip;float pi;struct compx v,w,t;LH=N/2; f=N;for(m=1;(f=f/2)!=1;m++){;} //2^m=N{for(L=m;L>=1;L--) //这里和时域的也有差别{ le=pow(2,L);B=le/2; //每一级碟形运算间隔的点数pi=3.14159; for(j=0;j<=B-1;j++) { p=pow(2,m-L)*j; ps=2*pi/N*p; w.real=cos(ps); w.imag=-sin(ps); for(i=j;i<=N-1;i=i+le) { ip=i+B; t=xin[i]; xin[i].real=xin[i].real+xin[ip].real; xin[i].imag=xin[i].imag+xin[ip].imag; xin[ip].real=xin[ip].real-t.real; xin[ip].imag=xin[ip].imag-t.imag; xin[ip]=EE(xin[ip],w); } }}}//变址运算nm=N-2; j=N/2;for(i=1;i<=nm;i++){if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}k=LH;while(j>=k){j=j-k;k=k/2;}j=j+k;}}//main programe#include <math.h>#include <stdio.h>#include <stdlib.h>float result[257];struct compx s[257];int Num=16;const float pp=3.14159;main(){int i;for(i=0;i<16;i++){s[i].real=sin(pp*i/32);s[i].imag=0;}FFT(s,Num);for(i=0;i<16;i++){printf("%.4f",s[i].real);printf("+%.4fj\n",s[i].imag);result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));}}
