今天的低功耗微控制器(μC)也開始集成原先只存在于大型微處理器、ASIC和DSP中的外設功能，使我們有可能以很低的功耗實現復雜的算術運算。本文討論一種快速傅立葉變換(FFT)應用，并在一個含有單周期硬件乘法器的低功耗μC上實現該應用。
這個FFT應用實時計算一路輸入電壓(圖1中的VIN)的頻譜。為完成該任務，用一片模數轉換器(ADC)對VIN進行采樣，獲得的采樣傳送給μC。然后，μC對這些采樣執行256點FFT運算，獲得輸入電壓的頻譜。為便于檢測，μC將計算出的頻譜數據傳送給PC，由PC實時顯示出來。
圖1. 利用FFT應用計算輸入電壓的頻譜。
該FFT應用的固件針對MAXQ2000系列中的一款16位、低功耗μC用C語言編寫。有興趣的讀者可以下載(ZIP，2.4kb)該項目的固件和電路原理圖。
背景知識
為確定輸入信號采樣的頻譜，我們需要對這些輸入采樣進行離散傅立葉變換(DFT)。DFT的定義如下：
其中N是采樣的數量，X(k)是頻譜，x(n)是一組輸入采樣。利用歐拉等式展開求和符，并分離輸入采樣和頻譜的實部和虛部，得到以下等式：
式2和3中，求和符中第二項的消失是由于輸入采樣全部為實數。假定我們有N個采樣，直接計算式2和3需要2N2次乘法和2N(N - 1)次加法。這樣，我們的256點輸入采樣DFT將需要進行131,072次乘法和130,560次加法運算。我們還是將注意力轉向FFT吧！
有多種FFT算法可供使用。本應用采用普通的radix-2算法，繼續將DFT分解為兩個更小的DFT。為此，N必須是2的指數。這種radix-2 FFT算法的步驟可歸納如圖2所示的蝶型運算。觀察這些蝶型運算我們可以發現，radix-2算法僅需(N / 2)log2(N)次乘法和Nlog2(N)次加法。圖2中用到的參數WN就是通常所謂的“旋轉因子”，可以在執行算法前預先計算出來。
圖2. 利用蝶型運算實現N = 8的FFT。
在圖2中，FFT的輸入顯示為一種特殊的排列順序，這種序列是對原始序列索引號的二進制位反轉后得到的。因此，當我們對N = 8個采樣執行radix-2 FFT算法時，需要將輸入數據的原始序列：
0 (000b), 1 (001b), 2 (010b), 3 (011b), 4 (100b), 5(101b), 6(110b), 7(111b)
重新排列為：
0 (000b), 4, (100b), 2 (010b), 6 (110b), 1 (001b), 5 (101), 3 (011), 7 (111)
FFT輸出則以正確的順序排列。圖2還說明，每個單獨的蝶型運算所得的結果，是下一級FFT運算所需的唯一數據。由于運算過程可“即位”進行，新值可替代舊值，這樣，計算N個采樣的FFT只需要2N個變量(因為每個數據都包括實部和虛部兩部分)。
& nbsp; FFT完成后，結果為復數形式。式4和5將結果轉換為極坐標方式后表示為：
有關DSP的文獻中可以找到很多優化方法，可使上述DFT/FFT算法更小或更快。其中最重要的一種優化方法(可能也是最容易實現的)源于這樣一個事實，那就是作為一個實數信號，其DFT幅度是相關于X(N / 2)對稱的，因此：
編寫FFT代碼絕非易事。低功耗μC的一些局限又進一步使該任務復雜化。
存儲器：我們所選的μC有2kB的RAM。已經知道該算法需要用到2N個16位變量來存儲FFT數據，這樣，我們的μC可以執行N最高為512的FFT。然而，固件的其他部分也要用到一些RAM。因此，在此項目中，我們限制N于256。若采用16位變量來表示每個值的實部和虛部，FFT數據總共需要1024字節的RAM。
速度：低功耗μC盡管具有高MIPS/mA性能，仍然需要一些優化手段來使運行FFT的指令數盡可能少。好在本應用所用的C編譯器(IAR的Embedded Workbench for MAXQ，見www.iar.com)可提供多種級別的優化和設置。高效地使用硬件乘法器可使代碼優化到可以接受的水平。
無浮點能力：所選的μC不具備浮點能力(低功耗產品一般都不具備浮點能力)。因此，所有運算都必須采用定點算法。為了表示小數，固件采用帶符號的Q8.7表示法。這樣，在固件中假定：
第0位至第6位代表小數部分
第7位至第14位代表整數部分
第15位代表符號位(二的補碼)
這樣的安排對于加法和減法沒有影響，但在做乘法時必須注意將數據按照Q8.7格式對齊。
所選的數據表示法還要適應FFT算法可能遇到的最大數值，同時又要提供足夠的精度。例如，我們的ADC可提供帶符號的8位采樣，以二的補碼表示。如果輸入為最大幅度(對于帶符號8位采樣為127)的直流電壓，則其能譜全部包含于X(0)中，用Q8.7表示為32512。這個數值能夠由單個帶符號的16位數據表示。
固件
以下部分討論在低功耗μC上執行radix-2 FFT的固件實現。信號采樣由ADC讀出后被存儲在x_n_re數組中。這個數組代表x(n)的實部。虛部存儲在x_n_im數組中，在開始運行FFT前初始化為零。完成FFT后，計算結果取代原始采樣數據，被存儲在x_n_re和x_n_im中。
獲取采樣
FFT算法假定采樣是以固定的取樣頻率獲得的。在為FFT獲取采樣時如果不加小心將會產生一些問題。例如，采樣間隔的抖動就會給FFT結果引入誤差，應盡力減小之。
ADC采樣循環中的判決語句會造成采樣間隔的抖動。例如，我們的系統從ADC讀取帶符號的8位采樣，并將其存儲在一組16位變量中。在下面的程序清單1中給出了兩種偽碼算法，執行這種ADC讀取-存儲功能。算法1給出的方法會造成采樣間隔的抖動，因為負采樣比正采樣需要更多的時間來讀取并存儲。
清單1. 兩種ADC采樣偽碼算法。第二種算法避免了第一種的問題——采樣間隔抖動。
// ALGORITHM 1: INCONSISTENT SAMPLING FREQUENCY - BAD!
// sample[] is an array of 16-bit va
riables
for i = 0 to (N-1)
begin
doADCSampleConversion() // Instruct ADC to sample Vin
sample[i] = read8BitSampleFromADC() // Read 8-bit sample from ADC
if (sample[i] & 0x0080) // If the 8-bit sample was negative
sample[i] = sample[i] + 0xFF00 // Make the 16-bit word negative
end
// ALGORITHM 2: FIXED SAMPLING FREQUENCY - GOOD!
// sample[] is an array of 16-bit variables
for i = 0 to (N-1)
begin
doADCSampleConversion() // Instruct ADC to sample Vin
sample[i] = read8BitSampleFromADC() // Read 8-bit sample from ADC
end
for i = 0 to (N-1)
begin
if (sample[i] & 0x0080) // If the 8-bit sample was negative
sample[i] = sample[i] + 0xFF00 // Make the 16-bit word negative
end
三角函數表
本FFT算法通過查表(LUT)而非計算得到正弦或余弦函數值。程序清單2給出了對于正弦和余弦LUT的申明。實際固件的注釋中包含了自動生成這些LUT的源代碼，可由程序調用。兩個LUT均含有N / 2分量，因為旋轉因子的索引號變化范圍為從0至N / 2 - 1 (見圖2)。
清單2. 正弦和余弦函數LUT。
const int cosLUT[N/2] = {+128,+127,+127, ... ,-127,-127,-127};
const int sinLUT[N/2] = {+0 ,+3 , +6, ... ,+9 , +6, +3};
這些LUT中的數組被聲明為const，強制編譯器將它們存儲于代碼空間而非數據空間。由于LUT數值須采用Q8.7表示法，它們由正弦和余弦的實際值乘以27后得到。
位反轉
位反轉排序(N已知)可在運行時通過計算、查表或直接利用展開循環編寫。所有這些方法都需要在源代碼的尺寸和運行速度間進行折衷。本FFT應用利用展開循環進行位反轉，其源代碼較長，但運行速度快。程序清單3顯示了該展開循環的實現。本應用固件的注釋中包含了用于程序自動生成展開循環的源代碼。
清單3. 用于實現N = 256的位反轉的展開循環。
i=x_n_re[ 1]; x_n_re[ 1]=x_n_re[128]; x_n_re[128]=i;
i=x_n_re[ 2]; x_n_re[ 2]=x_n_re[ 64]; x_n_re[ 64]=i;
i=x_n_re[ 3]; x_n_re[ 3]=x_n_re[192]; x_n_re[192]=i;
i=x_n_re[ 4]; x_n_re[ 4]=x_n_re[ 32]; x_n_re[ 32]=i;
...
i=x_n_re[207]; x_n_re[207]=x_n_re[243]; x_n_re[243]=i;
i=x_n_re[215]; x_n_re[215]=x_n_re[235]; x_n_re[235]=i;
i=x_n_re[223]; x_n_re[223]=x_n_re[251]; x_n_re[251]=i;
i=x_n_re[239]; x_n_re[239]=x_n_re[247]; x_n_re[247]=i;
Radix-2 FFT算法
采樣按照位反轉方式重新排序后就可進行FFT運算了。本radix-2 FFT應用的固件通過三個主循環執行圖2所示的蝶型運算。外循環計數log2(N)級FFT運算。內循環執行每一級的蝶型運算。
FFT
算法的核心部分是執行蝶型運算的一小塊代碼。程序清單4給出了這一塊代碼，遺憾的是，它是本應用中唯一“不可移植”的固件。宏MUL_1和MUL_2利用μC的硬件乘法器執行單指令周期乘法運算。這些宏的內容專用于MAXQ2000，可在實際固件中全部看到。
清單4. 用C編寫的蝶型運算。
/* (1) Macro MUL_1(A,B,C): C="A"*B (result in Q8.7)*/
/* (2) Macro MUL_2(A,C) : C="A"*last_B (result in Q8.7)*/
MUL_1(cosLUT[tf],x_n_re[b],resultMulReCos);
MUL_2(sinLUT[tf],resultMulReSin);
MUL_1(cosLUT[tf],x_n_im[b],resultMulImCos);
MUL_2(sinLUT[tf],resultMulImSin);
x_n_re[b] = x_n_re[a]-resultMulReCos+resultMulImSin;
x_n_im[b] = x_n_im[a]-resultMulReSin-resultMulImCos;
x_n_re[a] = x_n_re[a]+resultMulReCos-resultMulImSin;
x_n_im[a] = x_n_im[a]+resultMulReSin+resultMulImCos;
復數的極坐標轉換
為了便于確定VIN頻譜的幅度，我們須要將復數形式的X(k)轉換為極坐標形式。實現該轉換的固件示于程序清單5。幅度值取代了原始的FFT結果，因為固件不再需要這些數據。
清單5. FFT結果從復數形式轉換為極坐標形式。
const unsigned char magnLUT[16][16] =
{
{0x00,0x10,0x20, ... ,0xd0,0xe0,0xf0},
{0x10,0x16,0x23, ... ,0xd0,0xe0,0xf0},
...
{0xe0,0xe0,0xe2, ... ,0xff,0xff,0xff},
{0xf0,0xf0,0xf2, ... ,0xff,0xff,0xff}
};
...
...
/* Compute x_n_re=abs(x_n_re) and x_n_im=abs(x_n_im) */
...
...
x_n_re[0] = magnLUT[x_n_re[0]>>11][0];
for(i=1; i
x_n_re[N_DIV_2] = magnLUT[x_n_re[N_DIV_2]>>11][0];
頻譜幅度并非根據式4計算得到，而是通過一個二維LUT查表得到。第一索引為頻譜實部的高4位(MSB)，第二索引為頻譜虛部的高4位。為得到這些數據，可將帶符號的16位數據右移11次。在從頻譜的實部和虛部取得索引號前，需首先將它們轉換為絕對值。因此，符號位為零。
從式6我們已經知道，頻譜的幅度是關于X(N / 2)對稱的，因此我們只需將前(N / 2) + 1個頻譜數據轉換為極坐標形式。還有，我們可以看到，對于實數輸入采樣，X(0)和X(N / 2)的虛部總為零。因此這兩條譜線的幅度被單獨計算。本項目實際固件的注釋中包含了用于自動生成該LUT的源代碼，可由程序調用來計算X(k)的幅度。
Hamming或Hann窗
此項目固件還包括了對輸入采樣加Hamming或Hann窗的LUT (Q8.7格式)。加窗函數可有效降低對時域采樣x(n)的舍入操作所引起的頻譜泄漏。Hamming和Hann窗函數分別如式7和8所示。
程序清單6給出了實現這些函數的代碼。同樣，本項目實際固件的注釋中包含了用于自動生成這些LUT的源代碼，可由程序調用來實現這些窗函數。
清單6. 用來實現Hamming和Hann窗函數的LUT。
const char hammingLUT[N] = {+10, +10, +10, ... ,+10, +10, +10};
const char hannLUT[N] = { +0, +0, +0, ... , +0, +0, +0};
...
...
for(i=0; i<256; i++)
{
#ifdef WINDOWING_HAMMING
MUL_1(x_n_re[i],hammingLUT[i],x_n_re[i]); // x(n)*=hamming(n);
#endif
#ifdef WINDOWING_HANN
MUL_1(x_n_re[i],hannLUT[i]),x_n_re[i]); // x(n)*=hann(n);
#endif
}
測試結果
為了測試該FFT應用的性能，固件將X(k)幅度通過μC的UART端口上傳給PC。專門編寫的FFT Graph軟件(隨該項目固件一起提供)用于從PC串口讀取這些幅值，并以圖形方式實時顯示頻譜。圖3顯示了μC以200ksps采樣四種不同輸入信號并處理后，由FFT Graph所顯示出來的結果：
4.3V直流信號
50kHz正弦信號
70kHz正弦信號
6.25kHz方波
圖3. FFT Graph軟件顯示的由低功耗μC計算出的頻譜。
接下來干什么？
有興趣的讀者還可以花費大量的時間來繼續優化和重新配置該FFT應用。盡管在本文中我們選擇了radix-2算法，還有很多其他算法可以顯著降低加法和乘法運算量。很多本文所未提及的優化可以提升FFT的速度。例如，作為純實數的輸入采樣，其虛部總為零，頻譜中只有前半部分有實際意義。利用這一點，第一級和最后一級FFT的執行速度可進一步優化，但需要付出更多的程序空間。
總之，本文所討論的算法對于低功耗μC上的FFT應用而言，提供了一個很好的出發點。如果想了解更多信息和具體實現的細節，請查閱我們為本應用所提供的、帶有詳細注釋的固件信息。