C语言实现fft理论基础与工程应用的实例分析
三天的工厂实地监测，在师兄的帮助下，终于理解了原来似懂非懂的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。
这个频率点的幅值为：该点复数的模值除以N/2（n=1时是直流分量，其幅值是该点的模值除以N）。
下面先来简要分析下封装好的FFT的C程序包
1 /*********************************************************************
2 快速福利叶变换C程序包
3 函数简介：此程序包是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依
4 赖硬件。此程序包采用联合体的形式表示一个复数，输入为自然顺序的复
5 数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的
6 复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表，
7 以后的可采用查表法计算耗时较多的sin和cos运算，加快可计算速度.与
8 Ver1.1版相比较，Ver1.2版在创建正弦表时只建立了1/4个正弦波的采样值，
9 相比之下节省了FFT_N/4个存储空间
10 使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的
11 应该为2的N次方，不满足此条件时应在后面补0。若使用查表法计算sin值和
12 cos值，应在调用FFT函数前调用create_sin_tab()函数创建正弦表
13 函数调用：FFT(s);
14 作 者：吉帅虎
15 时 间：2010-2-20
16 版 本：Ver1.2
17 参考文献：
18 **********************************************************************/
19 #include <math.h>
20 #include "fft.h"
21
22 float *SIN_TAB;//定义正弦表的存放空间
23 int FFT_N = 1024;//定义采样点大小
24 /*******************************************************************
25 函数原型：struct compx EE(struct compx b1,struct compx b2)
26 函数功能：对两个复数进行乘法运算
27 输入参数：两个以联合体定义的复数a,b
28 输出参数：a和b的乘积，以联合体的形式输出
29 *******************************************************************/
30 struct compx EE(struct compx a,struct compx b)
31 {
32 struct compx c;
33 c.real=a.real*b.real-a.imag*b.imag;
34 c.imag=a.real*b.imag+a.imag*b.real;
35 return(c);
36 }
37
38 /******************************************************************
39 函数原型：void create_sin_tab(float *sin_t，int PointNum)
40 函数功能：创建一个正弦采样表，采样点数与福利叶变换点数相同
41 输入参数：*sin_t存放正弦表的数组指针,PointNum采样点数
42 输出参数：无
43 ******************************************************************/
44 void create_sin_tab(float *sin_t,int PointNum)
45 {
46 int i;
47 SIN_TAB=sin_t;
48 FFT_N=PointNum;
49 for(i=0;i<=FFT_N/4;i++)
50 SIN_TAB[i]=sin(2*PI*i/FFT_N);
51 }
52 /******************************************************************
53 函数原型：void sin_tab(float pi)
54 函数功能：采用查表的方法计算一个数的正弦值
55 输入参数：pi 所要计算正弦值弧度值，范围0--2*PI，不满足时需要转换
56 输出参数：输入值pi的正弦值
57 ******************************************************************/
58 float sin_tab(float pi)
59 {
60 int n=0;
61 float a=0;
62 n=(int)(pi*FFT_N/2/PI);
63
64 if(n>=0&&n<=FFT_N/4)
65 a=SIN_TAB[n];
66 else if(n>FFT_N/4&&n<FFT_N/2)
67 {
68 n-=FFT_N/4;
69 a=SIN_TAB[FFT_N/4-n];
70 }
71 else if(n>=FFT_N/2&&n<3*FFT_N/4)
72 {
73 n-=FFT_N/2;
74 a=-SIN_TAB[n];
75 }
76 else if(n>=3*FFT_N/4&&n<3*FFT_N)
77 {
78 n=FFT_N-n;
79 a=-SIN_TAB[n];
80 }
81
82 return a;
83 }
84 /******************************************************************
85 函数原型：void cos_tab(float pi)
86 函数功能：采用查表的方法计算一个数的余弦值
87 输入参数：pi 所要计算余弦值弧度值，范围0--2*PI，不满足时需要转换
88 输出参数：输入值pi的余弦值
89 ******************************************************************/
90 float cos_tab(float pi)
91 {
92 float a,pi2;
93 pi2=pi+PI/2;
94 if(pi2>2*PI)
95 pi2-=2*PI;
96 a=sin_tab(pi2);
97 return a;
98 }
99 /*****************************************************************
100 函数原型：void FFT(struct compx *xin)
101 函数功能：对输入的复数组进行快速傅里叶变换（FFT）
102 输入参数：*xin复数结构体组的首地址指针，struct型
103 输出参数：无
104 *****************************************************************/
105 void FFT(struct compx *xin)
106 {
107 int f,m,i,k,l,j=0;
108 register int nv2,nm1;
109 struct compx u,w,t;
110
111 nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
112 nm1=FFT_N-1;
113 for(i=0;i<nm1;i++)
114 {
115 if(i<j) //如果i<j,即进行变址
116 {
117 t=xin[j];
118 xin[j]=xin[i];
119 xin[i]=t;
120 }
121 k=nv2; //求j的下一个倒位序
122 while(k<=j) //如果k<=j,表示j的最高位为1
123 {
124 j=j-k; //把最高位变成0
125 k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
126 }
127 j=j+k; //把0改为1
128 }
129
130 {
131 int le,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算
132 f=FFT_N;
133 for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数
134 ;
135 for(m=1;m<=l;m++) // 控制蝶形结级数
136 { //m表示第m级蝶形，l为蝶形级总数l=log（2）N
137 le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点
138 lei=le/2; //同一蝶形结中参加运算的两点的距离
139 u.real=1.0; //u为蝶形结运算系数，初始值为1
140 u.imag=0.0;
141 //w.real=cos(PI/lei); //不适用查表法计算sin值和cos值
142 // w.imag=-sin(PI/lei);
143 w.real=cos_tab(PI/lei); //w为系数商，即当前系数与前一个系数的商
144 w.imag=-sin_tab(PI/lei);
145 for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结
146 {
147 for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结
148 {
149 ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点
150 t=EE(xin[ip],u); //蝶形运算，详见公式
151 xin[ip].real=xin[i].real-t.real;
152 xin[ip].imag=xin[i].imag-t.imag;
153 xin[i].real=xin[i].real+t.real;
154 xin[i].imag=xin[i].imag+t.imag;
155 }
156 u=EE(u,w); //改变系数，进行下一个蝶形运算
157 }
158 }
159 }
160 }
1 #ifndef FFT_H 2 #define FFT_H 3 //定义福利叶变换的点数 4 #define PI 3.1415926535897932384626433832795028841971 //定义圆周率值 5 6 struct compx {float real,imag;}; //定义一个复数结构 7 8 extern void create_sin_tab(float *sin_t,int PointNum); 9 extern void FFT(struct compx *xin); 10 11 #endif // FFT_H
程序主要分为两个部分，第一部分主要采用雷德算法来实现倒位序，第二部分主要是利用已经生成好的旋转矩阵做蝶形变换。（程序来源于《数字信号处理》清华大学出版社，程佩青，按时间抽选的基-2 FFT蝶形图）
1.第一个要解决的问题就是码位倒序，以雷德算法为例。
下面假如使用A[I]存的是顺序位序，而B[J]存的是倒位序。I<J的时候需要变序，I>J的时候就不用，不然就白忙活了。
倒位序 顺序 二进制表示 倒位序 顺序
①第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。下面对具体程序和实际情况分析。
1 void FFT_handle(u16 Input_voltage) 2 { 3 4 int i = 0; 5 6 s[t].real = Input_voltage; 7 s[t].imag = 0; 8 t++; 9 if(t >= NPT) 10 { 11 t = 0; 12 value_ACsignalMAX = 0; 13 FFT(s); 14 for(i=0;i<NPT / 2;i++) 15 { //求变换后结果的模值，存入复数的实部部分 16 // s[i].real=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2)); 17 table[i]=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2)); 18 19 if(i == 0) 20 { 21 value_DCsignal = table[i]; 22 } 23 else if(0 < i && i < NPT / 2) 24 { 25 if(table[i] > value_ACsignalMAX) 26 { 27 if(i < 15) // 滤除高频干扰 28 { 29 value_ACsignalMAX = table[i]; 30 MAX_palce = i; 31 } 32 else{} 33 } 34 } 35 else{} 36 37 jww=value_DCsignal; 38 jwww=value_ACsignalMAX; 39 }
首先要求出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。
