基2 FFT时间抽取和频域抽取算法

目录1、理解离散傅立叶变换2、基2频域算法原理3、基2 FFT的C语言实现(时域)4、基2 FFT的C语言实现(频域)1、理解离散傅立叶变换移步FFT结果的物理意义(编程参考用)2、基2频域算法原理原理请查看“按时间抽取基2 的 FFT算法的实现--杜义君”3、基2 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.imagb1.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(m1;(ff/2)!1;m){;} //求出m为log2 N nm N-2; j N/2; //变址运算对时间进行奇偶分解 for(i1;inm;i)//即xin第一位和最后一位不用操作不用变址其余各位根据码位倒置 { if(i k LH; while(jk){jj-k;kk/2;} j jk; } { for(L1;Lm;L)//运行m级蝶形运算 { lepow(2,L); Ble/2; pi3.14159; for(j0;jB-1;j) { ppow(2,m-L)*j; //k的取值012...(pow(2,L)/2)-1当在第一级时p只为0当在第二 //级时p为0和2当为第三级时p为0123.在递升的级别中不断延续见数字信号处理华中科技大学p83图的系数WN,且相邻不同种基本蝶形的蝶形节系数的增量为2*pi/N*pow(2,m-L) ps2*pi/N*p; w.realcos(ps); w.imag-sin(ps);//求出WNk for(ij;iN-1;iile)//确定与蝶形系数相乘的Xm(q)的下标m此方法为蝶形图的特点来得到 //相邻同种基本蝶形的间距为2的L次方。 { ipiB;//即对频谱进行前后分解 tEE(xin[ip],w);//复数相乘 xin[ip].realxin[i].real-t.real;//基本蝶式运算 xin[ip].imagxin[i].imag-t.imag; xin[i].realxin[i].realt.real; xin[i].imagxin[i].imagt.imag; } } } } return ; }//输入时域数据点为 num则输出频域数据点同为 numnum 是数据长度必须为 2 的整数次幕//其大小由数据采样定理来决定#include math.h #include stdio.h float result[257];//振幅其平方为功率谱 struct compx s[257]; int Num 16;//数据长度必须为2的整数次幕 const float pp 3.14159; void main(void) { int i; for(i0;i16;i) { s[i].realsin(pp*i/32); s[i].imag0; } FFT(s,Num); for(i0;i16;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); } }时域公式频域公式4、基2 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.imagb1.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(m1;(ff/2)!1;m){;} //2^mN { for(Lm;L1;L--) //这里和时域的也有差别 { le pow(2,L); B le/2; //每一级碟形运算间隔的点数 pi 3.14159; for(j0;jB-1;j) { p pow(2,m-L)*j; ps 2*pi/N*p; w.real cos(ps); w.imag -sin(ps); for(ij;iN-1;iile) { ip iB; t xin[i]; xin[i].real xin[i].realxin[ip].real; xin[i].imag xin[i].imagxin[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(i1;inm;i) { if(i k LH; while(jk){jj-k;kk/2;} j jk; } } //main programe #include #include #include float result[257]; struct compx s[257]; int Num 16; const float pp 3.14159; void main(void) { int i; for(i0;i16;i) { s[i].real sin(pp*i/32); s[i].imag 0; } FFT(s,Num); for(i0;i16;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)); } }星光不问赶路人岁月不负有心人。觉得不错动动发财的小手点个赞哦