三天的工廠實地監測，在師兄的幫助下，終於理解了原來似懂非懂的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。
