请查看按时抽取基2的原理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 N nm=N-2; j=N/2; ///变址操作，奇偶分解时间
for(i=1;i<=nm;i )//即xin第一个和最后一个不需要操作，不需要改变地址，其余的根据代码位置倒置 { if(i 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 #include #include
floatresult[257] structcompx s[257]; intNum=16.//数据长度必须是2的整数幕 const float pp=3.14159;
main() {
int i; for(i=0;i<16;i ) { s[i].real=sn(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); }}
时域公式
频域公式
----------------------------------------------------------------------------------------------------------------
基于时域和基于频域的算法很相似的，只不过时域里是先乘后加减，对输入序列进行倒序，而频域的算法先加减后乘，输入序列不用倒序，对输出序列进行倒序。调试成功之后，发现两者结果相差很小了。
#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 k=LH; while(j>=k){j=j-k;k=k/2;} j=j+k; }
}
//main programe
#include #include #include
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)); }
}
----------------------------------------------------------------------------------------------------------------