fft c代码以及工程应用实例
来源:互联网 发布:内网渗透 端口转发 编辑:程序博客网 时间:2024/06/06 01:48
转自:http://www.cnblogs.com/guluxuanyuan/p/4047771.html
三天的工厂实地监测,在师兄的帮助下,终于理解了原来似懂非懂的FFT变换的工程意义,废话少说,直入正题。
一、理论分析
快速傅里叶变换(Fast Fourier Transform)是离散傅里叶变换的一种快速算法,简称FFT,通过FFT可以将一个信号从时域变换到频域。
模拟信号经过A/D转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真,采样频率必须大于信号中最高频率成分的2倍,这称之为采样定理。
假设采样频率为fs,采样点数为N,那么FFT结果就是一个N点的复数,每一个点就对应着一个频率点,某一点n(n从1开始)表示的频率为:fn=(n-1)*fs/N。
举例说明:用1kHz的采样频率采样128点,则FFT结果的128个数据即对应的频率点分别是0,1k/128,2k/128,3k/128,…,127k/128 Hz。
fft.c
/********************************************************************* 快速福利叶变换C程序包函数简介:此程序包是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依 赖硬件。此程序包采用联合体的形式表示一个复数,输入为自然顺序的复 数(输入实数是可令复数虚部为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()函数创建正弦表函数调用:FFT(s);作 者:吉帅虎时 间:2010-2-20版 本:Ver1.2参考文献:**********************************************************************/#include <math.h>#include "fft.h"float *SIN_TAB;//定义正弦表的存放空间int FFT_N = 1024;//定义采样点大小 /*******************************************************************函数原型:struct compx EE(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);}/******************************************************************函数原型:void create_sin_tab(float *sin_t,int PointNum)函数功能:创建一个正弦采样表,采样点数与福利叶变换点数相同输入参数:*sin_t存放正弦表的数组指针,PointNum采样点数输出参数:无******************************************************************/void create_sin_tab(float *sin_t,int PointNum){ int i; SIN_TAB=sin_t; FFT_N=PointNum; for(i=0;i<=FFT_N/4;i++) SIN_TAB[i]=sin(2*PI*i/FFT_N);}/******************************************************************函数原型:void sin_tab(float pi)函数功能:采用查表的方法计算一个数的正弦值输入参数:pi 所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换输出参数:输入值pi的正弦值******************************************************************/float sin_tab(float pi){ int n=0; float a=0; n=(int)(pi*FFT_N/2/PI); if(n>=0&&n<=FFT_N/4) a=SIN_TAB[n]; else if(n>FFT_N/4&&n<FFT_N/2) { n-=FFT_N/4; a=SIN_TAB[FFT_N/4-n]; } else if(n>=FFT_N/2&&n<3*FFT_N/4) { n-=FFT_N/2; a=-SIN_TAB[n]; } else if(n>=3*FFT_N/4&&n<3*FFT_N) { n=FFT_N-n; a=-SIN_TAB[n]; } return a;}/******************************************************************函数原型:void cos_tab(float pi)函数功能:采用查表的方法计算一个数的余弦值输入参数:pi 所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换输出参数:输入值pi的余弦值******************************************************************/float cos_tab(float pi){ float a,pi2; pi2=pi+PI/2; if(pi2>2*PI) pi2-=2*PI; a=sin_tab(pi2); return a;}/*****************************************************************函数原型:void FFT(struct compx *xin)函数功能:对输入的复数组进行快速傅里叶变换(FFT)输入参数:*xin复数结构体组的首地址指针,struct型输出参数:无*****************************************************************/void FFT(struct compx *xin){ int f,m,i,k,l,j=0; register int nv2,nm1; struct compx u,w,t; nv2=FFT_N/2; //变址运算,即把自然顺序变成倒位序,采用雷德算法 nm1=FFT_N-1; for(i=0;i<nm1;i++) { if(i<j) //如果i<j,即进行变址 { t=xin[j]; xin[j]=xin[i]; xin[i]=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表示第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); for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结,即计算系数不同的蝶形结 { for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算,即计算系数相同蝶形结 { ip=i+lei; //i,ip分别表示参加蝶形运算的两个节点 t=EE(xin[ip],u); //蝶形运算,详见公式 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; } u=EE(u,w); //改变系数,进行下一个蝶形运算 } } }}
fft.h
#ifndef FFT_H#define FFT_H //定义福利叶变换的点数#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值struct compx {float real,imag;}; //定义一个复数结构extern void create_sin_tab(float *sin_t,int PointNum);extern void FFT(struct compx *xin);#endif // FFT_H
程序主要分为两个部分,第一部分主要采用雷德算法来实现倒位序,第二部分主要是利用已经生成好的旋转矩阵做蝶形变换。(程序来源于《数字信号处理》清华大学出版社,程佩青,按时间抽选的基-2 FFT蝶形图)
1.第一个要解决的问题就是码位倒序,以雷德算法为例。
下面假如使用A[I]存的是顺序位序,而B[J]存的是倒位序。I小于J的时候需要变序,I>J的时候就不用,不然就白忙活了。
例如 N = 8 的时候,
倒位序 顺序 二进制表示 倒位序 顺序
0 0 000 000
4 1 100 001
2 2 010 010
6 3 110 011
1 4 001 100
5 5 101 101
3 6 011 110
7 7 111 111
由上面的表可以看出,按自然顺序排列的二进制数,其下面一个数总是比其上面一个数大1,即下面一个数是上面一个数在最低位加1并向高位进位而得到的。而倒位序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位进位而得到。
I、J都是从0开始,若已知某个倒位序J,要求下一个倒位序数,则应先判断J的最高位是否为0,这可与k=N/2相比较,因为N/2总是等于100..的。如果k>J,则J的最高位为0,只要把该位变为1(J与k=N/2相加即可),就得到下一个倒位序数;如果K<=J,则J的最高位为1,可将最高位变为0(J与k=N/2相减即可)。然后还需判断次高位,这可与k=N\4相比较,若次高位为0,则需将它变为1(加N\4即可)其他位不变,既得到下一个倒位序数;若次高位是1,则需将它也变为0。然后再判断下一位。具体代码实现参见代码fft.c。
2.第二个要解决的问题就是蝶形运算
①第1级(第1列)每个蝶形的两节点“距离”为1,第2级每个蝶形的两节点“距离”为2,第3级每个蝶形的两节点“距离”为4,第4级每个蝶形的两节点“距离”为8。由此推得,第m级蝶形运算,每个蝶形的两节点“距离”L=2m-1。
②对于16点的FFT,第1级有16组蝶形,每组有1个蝶形;第2级有4组蝶形,每组有2个蝶形;第3级有2组蝶形,每组有4个蝶形;第4级有1组蝶形,每组有8个蝶形。由此可推出,对于N点的FFT,第m级有N/2L组蝶形,每组有L=2m-1个蝶形。
③旋转因子的确定
为提高FFT的运算速度,我们可以事先建立一个旋转因子数组,然后通过查表法来实现
complex code WN[N](旋转因子数组)
为节省CPU计算时间,旋转因子采用查表处理
根据实际FFT的点数N,该表数据需自行修改
以下结果通过Excel自动生成
WN[k].real=cos(2*PI/N*k)
WN[k].img=-sin(2*PI/N*k)
二、工程应用
为了找到基波位置,我们对第一次平均值做fft变换。第一次平均值即反映了ADC采样的值得整体趋势,只是在幅值上有所降低。ADC的定时器触发的采样频率为20KHz,40个采样点一次平均后的采样频率为500Hz,故fft的采样频率为500Hz(也即fft变换后最大频率不会达到500Hz)。
fft的采样点为1024点。故fft的分辨率约为0.49Hz。下面对具体程序和实际情况分析。
fft handle.c
void FFT_handle(u16 Input_voltage){ int i = 0; s[t].real = Input_voltage; s[t].imag = 0; t++; if(t >= NPT) { t = 0; value_ACsignalMAX = 0; FFT(s); for(i=0;i<NPT / 2;i++) { //求变换后结果的模值,存入复数的实部部分 // s[i].real=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2)); table[i]=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2)); if(i == 0) { value_DCsignal = table[i]; } else if(0 < i && i < NPT / 2) { if(table[i] > value_ACsignalMAX) { if(i < 15) // 滤除高频干扰 { value_ACsignalMAX = table[i]; MAX_palce = i; } else{} } } else{} jww=value_DCsignal; jwww=value_ACsignalMAX; }
首先要求出fft变化后的各频率点的模值,i=0对应直流分量,其次找到最大频率点(即基波),但程序界定i<15,故滤除高于15*1.49=7.4Hz的分量。
在debug全速运行时,程序稳定状态下测得Max i=9,故对应4.4Hz。
用信号发生器测得被ADC采样波形和将其做fft变换。可以看出基波频率为4.148Hz,而示波器fft变换后最大模值点为4.19Hz。
由此得出结论:模拟测试4.14Hz与程序结果4.4Hz(i=9)基本一致,误差来源于fft采样分辨率0.49Hz。
- fft c代码以及工程应用实例
- FFT算法应用实例
- FFT变换c代码
- FFT c代码的收集
- FFT和IFFT的C代码实现
- JS中封装解析xml的代码,以及实例应用
- fft代码
- FFT.c
- C实例---Makefile工程管理
- S3C2416裸机开发系列八_MDK启动代码工程应用实例
- S3C2416裸机开发系列九_GCC启动代码工程应用实例
- c++:sprintf应用实例
- C#----笔画应用实例
- opencl fft实例整理
- 队列代码C实例
- C语言实例代码
- C语言代码实例
- Redis常用命令以及代码实例
- Math.Round()方法可能与我们想象的大不一样(JS、C#)
- C#下载报错:An exception occurred during a WebClient request.,不知道如何处理。
- 关于做支付的笔记
- sqlite 获取当天当周当月等数据
- delphi webbrowser 去掉边框
- fft c代码以及工程应用实例
- 关于阅读体验
- SpringMVC+JQuery Ajax简单实例
- 模板_数组实现邻接表
- 前端css及js处理手记
- 解决AndroidStudio的logcat显示超长字符串的问题
- 一般处理程序里使用session对象为null,未将对象引用到实例化
- bootstrap的学习
- 使用ztree的心得