FFT算法实现流程 原理解释:
FFT算法的端口功能可理解为:输入的是ADC采样数据,输出得到对应频谱数据,这里介绍基2的时域抽取实现法。首先我们定义蝶形系数(又叫:旋转因子):WN=cos(2*pi/N)+j*sin(2*pi/N)。图1所示是一个蝶形运算基本单元,接下来就可以看看8点的FFT算法数据流图(图2)。
图1 蝶形运算单元
图2 8点基2时间抽取FFT算法流图
别被吓着了,耐心仔细观察规律,在图2中,显然包括图中三级的相似运算,来一级一级地查看:第一级的蝶形系数均为W0,蝶形节点的距离为1;第二级的蝶形系数为W0,WN/4,蝶形节点的距离为2;第三级的蝶形系数为W0,WN/8,W2N/8,W3N/8,蝶形节点的距离为4;第M级的蝶形系数为W0,W1, … ,W(N/2-1),蝶形节点的距离为N/2。清楚了这个递推关系,再根据图1的小碟形运算定义,整个FFT即可实现。 注意到没有,图2的左侧输入数据有点奇怪,原来这里并非按正常顺序排列,而是预处理为二进制的倒序啦。简单来说,二进制的倒序就是将一个整数按二进制一位一位颠倒过来就是了,图3表示了8点的情况,表1给出了对应关系。
图3 二进制倒序的说明图示
几点说明: 1. 注意该算法是复数的浮点运算,对于输入部分,实部为采样数据,虚部可设为0.0。 2. 输入数据数量要满足N=2 m,如N可取64,128等。当不满足时,可适当做少量的截取或补零处理。 3. 输入要尽量保证输入的是整数个周期波形数据,不能满足时,可考虑加窗处理(如常用的汉宁窗等)。 4. 对应B位的ADC数据,首先要除以2(B-1)化成浮点数再参加运算。 5. 输出数据也是复数,我们往往只关心幅度特性,所以要取复数的模。如:复数Z=a+j*b,则模为|Z|=sqrt(a2+b2)。 6. FFT得到的频谱,具有原理上的对称性,所以只需保存前一半即可。 7. 假设采样率为Fs,对N点数据做FFT,则得到的频谱分辨率Fs/N,每一根谱线依次代表0, Fs/N, 2Fs/N, 3Fs/N… 程序C代码: //FFT的C实现代码 #include <math.h> //包含数学函数库 #define N 240 //定义采集数据量 struct compx {float real,imag;}; //定义复数结构体 struct compx data1[N+1]; //定义变换数组 /********************************** //对采集的数据进行FFT变换 ***********************************/ //完成复数相乘 struct compx EE(struct compx b1,structcompx 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); } //实现FFT变换 /* Input: xin; Output: xin; Count: n */ void FFT(struct compx* xin,int n) { unsigned int f,m,nv2,nm1,i,k,j=1,l; struct compx v,w,t; nv2=n/2; //每一级有n/2个蝶形 f=n; for(m=1;(f=f/2)!=1;m++); //计算共有m级运算 nm1=n-1; //形成逆序数 for(i=1;i<=nm1;i++) { if(i<j) {t=xin[j]; xin[j]=xin; xin=t;} //只对一半进行交换 k=nv2; while(k<j) {j=j-k; k=k/2;} //如果高位是1则清0,同理查看次高位 j=j+k; //一直到次高位为0则置1 } //FFT的程序 { unsigned int le,lei,ip; float pi; for(l=1;l<=m;l++) //控制m级级数循环 { le=pow(2,l); //相同的旋转因子相距le lei=le/2; //共lei个不同的旋转因子 pi=3.14159265; //圆周率的精度位数 v.real=1.0; v.imag=0.0; //旋转因子的初始化 w.real=cos(pi/lei); w.imag=-sin(pi/lei); //每一级的单位旋转因子 for(j=1;j<=lei;j++) //不同的旋转因子进行循环 { for(i=j;i<=n;i=i+le) //相同的旋转因子的蝶形循环 { ip=i+lei; //每一蝶形的两输入相距lei t=EE(xin[ip],v); //旋转因子的形成 xin[ip].real=xin.real-t.real; //一个蝶形运算的完成 xin[ip].imag=xin.imag-t.imag; xin.real=xin.real+t.real; xin.imag=xin.imag+t.imag; } v=EE(v,w); //计算下一个旋转因子 } } } } |