《2022年C语言实现FFT .pdf》由会员分享，可在线阅读，更多相关《2022年C语言实现FFT .pdf（18页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、#include #include /* 快速福利叶变换C 函数函数简介：此函数是通用的快速傅里叶变换C 语言函数，移植性强，以下部分不依赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序的复数（输入实数是可令复数虚部为0） ，输出为经过FFT 变换的自然顺序的复数使用说明：使用此函数只需更改宏定义FFT_N 的值即可实现点数的改变，FFT_N 的应该为 2 的 N 次方，不满足此条件时应在后面补0 函数调用： FFT(s); 时间： 2010-2-20 版本： Ver1.0 参考文献：*/ #include #define PI 3.1415926535897932384626433
2、832795028841971 /定义圆周率值#define FFT_N 128 / 定义福利叶变换的点数struct compx float real,imag; /定义一个复数结构struct compx sFFT_N; /FFT 输入和输出：从 S1 开始存放，根据大小自己定义/* 函数原型： struct compx EE(struct compx b1,struct compx b2) 函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数a,b 输出参数： a 和 b 的乘积，以联合体的形式输出*/ struct compx EE(struct compx a,struc
3、t compx b) struct compx c; c.real=a.real*b.real-a.imag*b.imag; c.imag=a.real*b.imag+a.imag*b.real; return(c); /* 函数原型： void FFT(struct compx *xin,int N) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 18 页 - - - - - - - - - 函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数： *xin 复
4、数结构体组的首地址指针，struct 型*/ void FFT(struct compx *xin) int f,m,nv2,nm1,i,k,l,j=0; struct compx u,w,t; nv2=FFT_N/2; /变址运算，即把自然顺序变成倒位序，采用雷德算法nm1=FFT_N-1; for(i=0;inm1;i+) if(ij) /如果 ij, 即进行变址 t=xinj; xinj=xini; xini=t; k=nv2; /求 j 的下一个倒位序while(k=j) /如果 k=j,表示 j 的最高位为1 j=j-k; /把最高位变成0 k=k/2; /k/2，比较次高位，依次类
5、推，逐个比较，直到某个位为0 j=j+k; /把 0 改为 1 int le,lei,ip; /FFT 运算核，使用蝶形运算完成FFT 运算f=FFT_N; for(l=1;(f=f/2)!=1;l+) /计算 l 的值，即计算蝶形级数; for(m=1;m=l;m+) / 控制蝶形结级数 /m 表示第m 级蝶形， l 为蝶形级总数l=log（2） N le=2(m-1); /le 蝶形结距离，即第m 级蝶形的蝶形结相距 le 点lei=le/2; /同一蝶形结中参加运算的两点的距离u.real=1.0; /u 为蝶形结运算系数，初始值为1 u.imag=0.0; w.real=cos(PI/
6、lei); /w 为系数商，即当前系数与前一个系数的商w.imag=-sin(PI/lei); for(j=0;j=lei-1;j+) /控制计算不同种蝶形结，即计算系数不同的名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 18 页 - - - - - - - - - 蝶形结 for(i=j;i=FFT_N-1;i=i+le) /控制同一蝶形结运算，即计算系数相同蝶形结 ip=i+lei; /i，ip 分别表示参加蝶形运算的两个节点t=EE(xinip,u); /蝶形运算
7、，详见公式xinip.real=xini.real-t.real; xinip.imag=xini.imag-t.imag; xini.real=xini.real+t.real; xini.imag=xini.imag+t.imag; u=EE(u,w); /改变系数，进行下一个蝶形运算 /* 函数原型： void main() 函数功能：测试FFT 变换，演示函数使用方法输入参数：无输出参数：无*/ void main() int i; for(i=0;iFFT_N;i+) /给结构体赋值 si.real=sin(2*3.141592653589793*i/FFT_N); /实部为正弦波F
8、FT_N 点采样，赋值为1 si.imag=0; /虚部为 0 FFT(s); /进行快速福利叶变换for(i=0;iFFT_N;i+) /求变换后结果的模值，存入复数的实部部分si.real=sqrt(si.real*si.real+si.imag*si.imag); while(1); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 18 页 - - - - - - - - - #include #include /* 快速福利叶变换C 程序包函数简介：此程序包是通用
9、的快速傅里叶变换C 语言函数，移植性强，以下部分不依赖硬件。此程序包采用联合体的形式表示一个复数，输入为自然顺序的复数（输入实数是可令复数虚部为0） ，输出为经过FFT 变换的自然顺序的复数 .此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表，以后的可采用查表法计算耗时较多的sin 和 cos 运算，加快可计算速度使用说明：使用此函数只需更改宏定义FFT_N 的值即可实现点数的改变，FFT_N 的应该为 2 的 N 次方，不满足此条件时应在后面补0。若使用查表法计算sin 值和cos 值，应在调用FFT 函数前调用create_sin_tab()函数创建正弦表函数调用
10、： FFT(s); 时间： 2010-2-20 版本： Ver1.1 参考文献：*/ #include #define FFT_N 128 / 定义福利叶变换的点数#define PI 3.1415926535897932384626433832795028841971 /定义圆周率值struct compx float real,imag; /定义一个复数结构struct compx sFFT_N; /FFT 输入和输出：从S0开始存放，根据大小自己定义float SIN_TABFFT_N/2; /定义正弦表的存放空间/* 函数原型： struct compx EE(struct compx
11、 b1,struct compx b2) 函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数a,b 输出参数： a 和 b 的乘积，以联合体的形式输出*/ struct compx EE(struct compx a,struct compx b) struct compx c; c.real=a.real*b.real-a.imag*b.imag; c.imag=a.real*b.imag+a.imag*b.real; return(c); /* 函数原型： void create_sin_tab(float *sin_t) 名师资料总结 - - -精品资料欢迎下载 - - -
12、 - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 18 页 - - - - - - - - - 函数功能：创建一个正弦采样表，采样点数与福利叶变换点数相同输入参数： *sin_t 存放正弦表的数组指针输出参数：无*/ void create_sin_tab(float *sin_t) int i; for(i=0;i=0&n=FFT_N/2&n2*PI) pi2-=2*PI; a=sin_tab(pi2); return a; /* 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - -
13、 - - - - - 名师精心整理 - - - - - - - 第 5 页，共 18 页 - - - - - - - - - 函数原型： void FFT(struct compx *xin,int N) 函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数： *xin 复数结构体组的首地址指针，struct 型输出参数：无*/ void FFT(struct compx *xin) int f,m,nv2,nm1,i,k,l,j=0; struct compx u,w,t; nv2=FFT_N/2; /变址运算，即把自然顺序变成倒位序，采用雷德算法nm1=FFT_N-1; for(i
14、=0;inm1;i+) if(ij) /如果 ij, 即进行变址 t=xinj; xinj=xini; xini=t; k=nv2; /求 j 的下一个倒位序while(k=j) /如果 k=j,表示 j 的最高位为1 j=j-k; /把最高位变成0 k=k/2; /k/2，比较次高位，依次类推，逐个比较，直到某个位为0 j=j+k; /把 0 改为 1 int le,lei,ip; /FFT 运算核，使用蝶形运算完成FFT 运算f=FFT_N; for(l=1;(f=f/2)!=1;l+) /计算 l 的值，即计算蝶形级数; for(m=1;m=l;m+) / 控制蝶形结级数 /m 表示第
15、m 级蝶形， l 为蝶形级总数l=log（2）N le=2(m-1); /le 蝶形结距离，即第m 级蝶形的蝶形结相距le 点lei=le/2; /同一蝶形结中参加运算的两点的距离u.real=1.0; /u 为蝶形结运算系数，初始值为1 u.imag=0.0; /w.real=cos(PI/lei); /不适用查表法计算sin 值和 cos 值/ w.imag=-sin(PI/lei); w.real=cos_tab(PI/lei); /w 为系数商，即当前系数与前一个系数的商w.imag=-sin_tab(PI/lei); 名师资料总结 - - -精品资料欢迎下载 - - - - - -
16、- - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 18 页 - - - - - - - - - for(j=0;j=lei-1;j+) /控制计算不同种蝶形结，即计算系数不同的蝶形结 for(i=j;i=FFT_N-1;i=i+le) /控制同一蝶形结运算，即计算系数相同蝶形结 ip=i+lei; /i，ip 分别表示参加蝶形运算的两个节点t=EE(xinip,u); /蝶形运算，详见公式xinip.real=xini.real-t.real; xinip.imag=xini.imag-t.imag; xini.real=xini.rea
17、l+t.real; xini.imag=xini.imag+t.imag; u=EE(u,w); /改变系数，进行下一个蝶形运算 /* 函数原型： void main() 函数功能：测试FFT 变换，演示函数使用方法输入参数：无输出参数：无*/ void main() int i; create_sin_tab(SIN_TAB); for(i=0;iFFT_N;i+) /给结构体赋值 si.real=sin(2*3.141592653589793*i/FFT_N); /实部为正弦波FFT_N 点采样，赋值为1 si.imag=0; /虚部为 0 FFT(s); /进行快速福利叶变换for(i=
18、0;iFFT_N;i+) /求变换后结果的模值，存入复数的实部部分si.real=sqrt(si.real*si.real+si.imag*si.imag); while(1); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 18 页 - - - - - - - - - #include #include /* 快速福利叶变换C 程序包函数简介：此程序包是通用的快速傅里叶变换C 语言函数，移植性强，以下部分不依赖硬件。此程序包采用联合体的形式表示一个复数，输入为自然顺序
19、的复数（输入实数是可令复数虚部为0） ，输出为经过FFT 变换的自然顺序的复数 .此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表，以后的可采用查表法计算耗时较多的sin 和 cos 运算，加快可计算速度.与Ver1.1 版相比较， Ver1.2 版在创建正弦表时只建立了1/4 个正弦波的采样值，相比之下节省了FFT_N/4 个存储空间使用说明：使用此函数只需更改宏定义FFT_N 的值即可实现点数的改变，FFT_N 的应该为 2 的 N 次方，不满足此条件时应在后面补0。若使用查表法计算sin 值和cos 值，应在调用FFT 函数前调用create_sin_tab()
20、函数创建正弦表函数调用： FFT(s); 时间： 2010-2-20 版本： Ver1.2 参考文献：*/ #include #define FFT_N 128 /定义福利叶变换的点数#define PI 3.1415926535897932384626433832795028841971 / 定义圆周率值struct compx float real,imag; /定义一个复数结构struct compx sFFT_N; /FFT 输入和输出：从S0开始存放，根据大小自己定义float SIN_TABFFT_N/4+1; /定义正弦表的存放空间/* 函数原型： struct compx EE
21、(struct compx b1,struct compx b2) 函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数a,b 输出参数： a 和 b 的乘积，以联合体的形式输出*/ struct compx EE(struct compx a,struct compx b) struct compx c; c.real=a.real*b.real-a.imag*b.imag; c.imag=a.real*b.imag+a.imag*b.real; return(c); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - -
22、 名师精心整理 - - - - - - - 第 8 页，共 18 页 - - - - - - - - - /* 函数原型： void create_sin_tab(float *sin_t) 函数功能：创建一个正弦采样表，采样点数与福利叶变换点数相同输入参数： *sin_t 存放正弦表的数组指针输出参数：无*/ void create_sin_tab(float *sin_t) int i; for(i=0;i=0&nFFT_N/4&n=FFT_N/2&n=3*FFT_N/4&n2*PI) pi2-=2*PI; a=sin_tab(pi2); return a; /* 函数原型： void F
23、FT(struct compx *xin,int N) 函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数： *xin 复数结构体组的首地址指针，struct 型输出参数：无*/ void FFT(struct compx *xin) int f,m,nv2,nm1,i,k,l,j=0; struct compx u,w,t; nv2=FFT_N/2; /变址运算，即把自然顺序变成倒位序，采用雷德算法nm1=FFT_N-1; for(i=0;inm1;i+) if(ij) /如果 ij, 即进行变址 t=xinj; xinj=xini; xini=t; k=nv2; /求 j 的下一
24、个倒位序while(k=j) /如果 k=j,表示 j 的最高位为1 j=j-k; /把最高位变成0 k=k/2; /k/2，比较次高位，依次类推，逐个比较，直到某个位为0 j=j+k; /把 0 改为 1 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 18 页 - - - - - - - - - int le,lei,ip; /FFT 运算核，使用蝶形运算完成FFT 运算f=FFT_N; for(l=1;(f=f/2)!=1;l+) /计算 l 的值，即计算蝶形级数
25、; for(m=1;m=l;m+) / 控制蝶形结级数 /m 表示第m 级蝶形， l 为蝶形级总数l=log（2） N le=2(m-1); /le 蝶形结距离，即第m 级蝶形的蝶形结相距 le 点lei=le/2; /同一蝶形结中参加运算的两点的距离u.real=1.0; /u 为蝶形结运算系数，初始值为1 u.imag=0.0; /w.real=cos(PI/lei); /不适用查表法计算sin 值和 cos 值/ w.imag=-sin(PI/lei); w.real=cos_tab(PI/lei); /w 为系数商，即当前系数与前一个系数的商w.imag=-sin_tab(PI/lei
26、); for(j=0;j=lei-1;j+) /控制计算不同种蝶形结，即计算系数不同的蝶形结 for(i=j;i=FFT_N-1;i=i+le) /控制同一蝶形结运算，即计算系数相同蝶形结 ip=i+lei; /i，ip 分别表示参加蝶形运算的两个节点t=EE(xinip,u); /蝶形运算，详见公式xinip.real=xini.real-t.real; xinip.imag=xini.imag-t.imag; xini.real=xini.real+t.real; xini.imag=xini.imag+t.imag; u=EE(u,w); /改变系数，进行下一个蝶形运算 /* 函数原型：
27、 void main() 函数功能：测试FFT 变换，演示函数使用方法输入参数：无输出参数：无名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 18 页 - - - - - - - - - */ void main() int i; create_sin_tab(SIN_TAB); for(i=0;iFFT_N;i+) /给结构体赋值 si.real=sin(2*3.141592653589793*i/FFT_N); /实部为正弦波FFT_N 点采样，赋值为1 si.im
28、ag=0; /虚部为 0 FFT(s); /进行快速福利叶变换for(i=0;iFFT_N;i+) /求变换后结果的模值，存入复数的实部部分si.real=sqrt(si.real*si.real+si.imag*si.imag); while(1); #include struct compx float real ; float imag ; /* 定义一个复数结构struct compx s257; /*FFT 输入，输出均从 s1开始1.struct compx EE( struct compx, struct compx); 2.void FFT(struct compx * ,in
29、t) 3.4.struct compx EE( struct compx b1, struct compx b2) 5. struct compx b3; 6.b3.real =b1.real*b2.real-b1.imag*b2.imag; 7. b3.imag=b1.real*b2.imag+b1.imag*b2.real; 8. return( b3); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 12 页，共 18 页 - - - - - - - - - 9. 10./*
30、输入 xin 的实部，虚部，输出：xin（实部，虚部）N 为 FFT 点数11.12.void FFT( struct compx *xin ,N) 13. 14.int f,m,nv2,nm1,i,k,j=1,l; 15.int le,lei,ip; 16.float pi,x,y; 17.struct compx v,w,t; 18.nv2=N/2; 19.f=N; 20.for(m=1;(f=f/2)!=1;m+) ; /*求 N 为 2 的多少次幂放入m 21.nm1=N-1; 22. /* 变变址运算，实现数据运算前的位倒序*/ 23.for(i=1;i=nm1;i+)24. if(
31、ij) t=xinj;xinj=xini;xini=t; 25. k=nv2; 26. while(kj) j=j-k;k=k/2; 27. j=j+k; 28. 29./* FFT*/ 30.for(l=1;lm;l+) 31. 32.le=pow(2,l); 33.lei=le/2; 34.pi=3.14159265; 35.v.real=1.0;v.imag=0.0; 36.w.real=cos(pi/lei);w.imag=-sin(pi/lei); 37. for(j=1;j=lei;i+) 38. 39. for(i=j;i=N;i=i+le) ip=i+lei; 40. t=EE
32、(xinip,v); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 13 页，共 18 页 - - - - - - - - - 41. xinip.real=xini.real-t.real; 42. xinip.imag=xini.imag-t.imag; 43. xini.real=xini.real+t.real; 44. xini.imag=xini.imag+t.imag; 45. v=EE(v,w); 46. 47.48. 49. 注意（用的是基二时域算法）且点数是的n
33、 次方#include math.h #include stdio.hstruct compx double real; double imag; ;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; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - -
34、 第 14 页，共 18 页 - - - - - - - - - return(b3); FFT 运算函数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+); nm=N-2; j=N/2; for(i=1;i=nm;i+) if(i=k)j=j-k;k=k/2; j=j+k; for(L=1;L=m;L+) le=pow(2,L); B=le/2
35、; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 15 页，共 18 页 - - - - - - - - - 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=EE(xinip,w); xinip.real=xini.real-t.real; xinip.imag=xini.imag-t
36、.imag; xini.real=xini.real+t.real; xini.imag=xini.imag+t.imag; return ; 主函数#include #include #include float result257; struct compx s257; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 16 页，共 18 页 - - - - - - - - - int Num= ; const float pp=3.14159; main() int i; for
37、(i=0;i ;i+) si.real=sin(pp*i/16); si.imag=0; FFT(s,Num); for(i=0;i ;i+) printf(%.4f,si.real); printf(+%.4fjn,si.imag); resulti=sqrt(pow(si.real,2)+pow(si.imag,2); 其函数中取的是个点for(i=0;i ;i+)如果要改变点数或者函数只需在主函数中改变for(i=0;i ;i+)si.real=sin(pp*i/16)名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 17 页，共 18 页 - - - - - - - - - 便可名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 18 页，共 18 页 - - - - - - - - -