今天的低功耗微控制器(μC)也開始集成原先只存在于大型微處理器,、ASIC和DSP中的外設(shè)功能,，使我們有可能以很低的功耗實(shí)現(xiàn)復(fù)雜的算術(shù)運(yùn)算,。本文討論一種快速傅立葉變換(FFT)應(yīng)用,，并在一個(gè)含有單周期硬件乘法器的低功耗μC上實(shí)現(xiàn)該應(yīng)用,。 

       這個(gè)FFT應(yīng)用實(shí)時(shí)計(jì)算一路輸入電壓(圖1中的VIN)的頻譜。為完成該任務(wù),，用一片模數(shù)轉(zhuǎn)換器(ADC)對(duì)VIN進(jìn)行采樣,，獲得的采樣傳送給μC。然后,，μC對(duì)這些采樣執(zhí)行256點(diǎn)FFT運(yùn)算,，獲得輸入電壓的頻譜。為便于檢測(cè),，μC將計(jì)算出的頻譜數(shù)據(jù)傳送給PC,，由PC實(shí)時(shí)顯示出來。

                             
                                        圖1. 利用FFT應(yīng)用計(jì)算輸入電壓的頻譜,。 

       該FFT應(yīng)用的固件針對(duì)MAXQ2000系列中的一款16位,、低功耗μC用C語言編寫。有興趣的讀者可以下載(ZIP,，2.4kb)該項(xiàng)目的固件和電路原理圖,。

       背景知識(shí)

       為確定輸入信號(hào)采樣的頻譜，我們需要對(duì)這些輸入采樣進(jìn)行離散傅立葉變換(DFT),。DFT的定義如下：

       其中N是采樣的數(shù)量,，X(k)是頻譜，x(n)是一組輸入采樣,。利用歐拉等式展開求和符,，并分離輸入采樣和頻譜的實(shí)部和虛部，得到以下等式：

       式2和3中,，求和符中第二項(xiàng)的消失是由于輸入采樣全部為實(shí)數(shù),。假定我們有N個(gè)采樣，直接計(jì)算式2和3需要2N2次乘法和2N(N - 1)次加法,。這樣,，我們的256點(diǎn)輸入采樣DFT將需要進(jìn)行131,072次乘法和130,560次加法運(yùn)算。我們還是將注意力轉(zhuǎn)向FFT吧,！

       有多種FFT算法可供使用,。本應(yīng)用采用普通的radix-2算法，繼續(xù)將DFT分解為兩個(gè)更小的DFT,。為此,，N必須是2的指數(shù)。這種radix-2 FFT算法的步驟可歸納如圖2所示的蝶型運(yùn)算,。觀察這些蝶型運(yùn)算我們可以發(fā)現(xiàn),，radix-2算法僅需(N / 2)log2(N)次乘法和Nlog2(N)次加法。圖2中用到的參數(shù)WN就是通常所謂的“旋轉(zhuǎn)因子”,，可以在執(zhí)行算法前預(yù)先計(jì)算出來,。

                                  圖2. 利用蝶型運(yùn)算實(shí)現(xiàn)N = 8的FFT。

       在圖2中,，F(xiàn)FT的輸入顯示為一種特殊的排列順序,，這種序列是對(duì)原始序列索引號(hào)的二進(jìn)制位反轉(zhuǎn)后得到的,。因此，當(dāng)我們對(duì)N = 8個(gè)采樣執(zhí)行radix-2 FFT算法時(shí),，需要將輸入數(shù)據(jù)的原始序列：

       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還說明，每個(gè)單獨(dú)的蝶型運(yùn)算所得的結(jié)果,，是下一級(jí)FFT運(yùn)算所需的唯一數(shù)據(jù),。由于運(yùn)算過程可“即位”進(jìn)行，新值可替代舊值,，這樣,，計(jì)算N個(gè)采樣的FFT只需要2N個(gè)變量(因?yàn)槊總€(gè)數(shù)據(jù)都包括實(shí)部和虛部?jī)刹糠?。

  & nbsp;    FFT完成后,，結(jié)果為復(fù)數(shù)形式,。式4和5將結(jié)果轉(zhuǎn)換為極坐標(biāo)方式后表示為：

       有關(guān)DSP的文獻(xiàn)中可以找到很多優(yōu)化方法，可使上述DFT/FFT算法更小或更快,。其中最重要的一種優(yōu)化方法(可能也是最容易實(shí)現(xiàn)的)源于這樣一個(gè)事實(shí),，那就是作為一個(gè)實(shí)數(shù)信號(hào)，其DFT幅度是相關(guān)于X(N / 2)對(duì)稱的,，因此：

       編寫FFT代碼絕非易事,。低功耗μC的一些局限又進(jìn)一步使該任務(wù)復(fù)雜化。

       存儲(chǔ)器：我們所選的μC有2kB的RAM,。已經(jīng)知道該算法需要用到2N個(gè)16位變量來存儲(chǔ)FFT數(shù)據(jù),，這樣，我們的μC可以執(zhí)行N最高為512的FFT,。然而,，固件的其他部分也要用到一些RAM。因此,，在此項(xiàng)目中,，我們限制N于256。若采用16位變量來表示每個(gè)值的實(shí)部和虛部,，F(xiàn)FT數(shù)據(jù)總共需要1024字節(jié)的RAM,。

       速度：低功耗μC盡管具有高M(jìn)IPS/mA性能，仍然需要一些優(yōu)化手段來使運(yùn)行FFT的指令數(shù)盡可能少,。好在本應(yīng)用所用的C編譯器(IAR的Embedded Workbench for MAXQ,，見www.iar.com)可提供多種級(jí)別的優(yōu)化和設(shè)置。高效地使用硬件乘法器可使代碼優(yōu)化到可以接受的水平,。

       無浮點(diǎn)能力：所選的μC不具備浮點(diǎn)能力(低功耗產(chǎn)品一般都不具備浮點(diǎn)能力),。因此，所有運(yùn)算都必須采用定點(diǎn)算法。為了表示小數(shù),，固件采用帶符號(hào)的Q8.7表示法,。這樣，在固件中假定：

       第0位至第6位代表小數(shù)部分 
       第7位至第14位代表整數(shù)部分 
       第15位代表符號(hào)位(二的補(bǔ)碼) 

       這樣的安排對(duì)于加法和減法沒有影響,，但在做乘法時(shí)必須注意將數(shù)據(jù)按照Q8.7格式對(duì)齊,。

       所選的數(shù)據(jù)表示法還要適應(yīng)FFT算法可能遇到的最大數(shù)值,，同時(shí)又要提供足夠的精度,。例如，我們的ADC可提供帶符號(hào)的8位采樣,，以二的補(bǔ)碼表示,。如果輸入為最大幅度(對(duì)于帶符號(hào)8位采樣為127)的直流電壓，則其能譜全部包含于X(0)中,，用Q8.7表示為32512,。這個(gè)數(shù)值能夠由單個(gè)帶符號(hào)的16位數(shù)據(jù)表示。

       固件

       以下部分討論在低功耗μC上執(zhí)行radix-2 FFT的固件實(shí)現(xiàn),。信號(hào)采樣由ADC讀出后被存儲(chǔ)在x_n_re數(shù)組中,。這個(gè)數(shù)組代表x(n)的實(shí)部。虛部存儲(chǔ)在x_n_im數(shù)組中,，在開始運(yùn)行FFT前初始化為零,。完成FFT后，計(jì)算結(jié)果取代原始采樣數(shù)據(jù),，被存儲(chǔ)在x_n_re和x_n_im中,。

       獲取采樣 

       FFT算法假定采樣是以固定的取樣頻率獲得的。在為FFT獲取采樣時(shí)如果不加小心將會(huì)產(chǎn)生一些問題,。例如,，采樣間隔的抖動(dòng)就會(huì)給FFT結(jié)果引入誤差，應(yīng)盡力減小之,。
ADC采樣循環(huán)中的判決語句會(huì)造成采樣間隔的抖動(dòng),。例如，我們的系統(tǒng)從ADC讀取帶符號(hào)的8位采樣,，并將其存儲(chǔ)在一組16位變量中,。在下面的程序清單1中給出了兩種偽碼算法，執(zhí)行這種ADC讀取-存儲(chǔ)功能,。算法1給出的方法會(huì)造成采樣間隔的抖動(dòng),，因?yàn)樨?fù)采樣比正采樣需要更多的時(shí)間來讀取并存儲(chǔ)。

       清單1. 兩種ADC采樣偽碼算法,。第二種算法避免了第一種的問題——采樣間隔抖動(dòng),。

// 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

       三角函數(shù)表

       本FFT算法通過查表(LUT)而非計(jì)算得到正弦或余弦函數(shù)值。程序清單2給出了對(duì)于正弦和余弦LUT的申明。實(shí)際固件的注釋中包含了自動(dòng)生成這些LUT的源代碼,，可由程序調(diào)用,。兩個(gè)LUT均含有N / 2分量，因?yàn)樾D(zhuǎn)因子的索引號(hào)變化范圍為從0至N / 2 - 1 (見圖2),。
清單2. 正弦和余弦函數(shù)LUT,。

const int cosLUT[N/2] = {+128,+127,+127, ... ,-127,-127,-127};
const int sinLUT[N/2] = {+0 ,+3 , +6, ... ,+9 , +6, +3};

       這些LUT中的數(shù)組被聲明為const，強(qiáng)制編譯器將它們存儲(chǔ)于代碼空間而非數(shù)據(jù)空間,。由于LUT數(shù)值須采用Q8.7表示法,，它們由正弦和余弦的實(shí)際值乘以27后得到。

       位反轉(zhuǎn)

       位反轉(zhuǎn)排序(N已知)可在運(yùn)行時(shí)通過計(jì)算,、查表或直接利用展開循環(huán)編寫,。所有這些方法都需要在源代碼的尺寸和運(yùn)行速度間進(jìn)行折衷。本FFT應(yīng)用利用展開循環(huán)進(jìn)行位反轉(zhuǎn),，其源代碼較長(zhǎng),，但運(yùn)行速度快。程序清單3顯示了該展開循環(huán)的實(shí)現(xiàn),。本應(yīng)用固件的注釋中包含了用于程序自動(dòng)生成展開循環(huán)的源代碼,。
清單3. 用于實(shí)現(xiàn)N = 256的位反轉(zhuǎn)的展開循環(huán)。

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算法

       采樣按照位反轉(zhuǎn)方式重新排序后就可進(jìn)行FFT運(yùn)算了,。本radix-2 FFT應(yīng)用的固件通過三個(gè)主循環(huán)執(zhí)行圖2所示的蝶型運(yùn)算,。外循環(huán)計(jì)數(shù)log2(N)級(jí)FFT運(yùn)算。內(nèi)循環(huán)執(zhí)行每一級(jí)的蝶型運(yùn)算,。
FFT

算法的核心部分是執(zhí)行蝶型運(yùn)算的一小塊代碼,。程序清單4給出了這一塊代碼，遺憾的是,，它是本應(yīng)用中唯一“不可移植”的固件,。宏MUL_1和MUL_2利用μC的硬件乘法器執(zhí)行單指令周期乘法運(yùn)算。這些宏的內(nèi)容專用于MAXQ2000,，可在實(shí)際固件中全部看到,。

       清單4. 用C編寫的蝶型運(yùn)算。

/* (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;

       復(fù)數(shù)的極坐標(biāo)轉(zhuǎn)換

       為了便于確定VIN頻譜的幅度,，我們須要將復(fù)數(shù)形式的X(k)轉(zhuǎn)換為極坐標(biāo)形式,。實(shí)現(xiàn)該轉(zhuǎn)換的固件示于程序清單5。幅度值取代了原始的FFT結(jié)果,，因?yàn)楣碳辉傩枰@些數(shù)據(jù),。
清單5. FFT結(jié)果從復(fù)數(shù)形式轉(zhuǎn)換為極坐標(biāo)形式。

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; ix_n_re[i] = magnLUT[x_n_re[i]>>11][x_n_im[i]>>11];

x_n_re[N_DIV_2] = magnLUT[x_n_re[N_DIV_2]>>11][0];

       頻譜幅度并非根據(jù)式4計(jì)算得到,，而是通過一個(gè)二維LUT查表得到,。第一索引為頻譜實(shí)部的高4位(MSB),，第二索引為頻譜虛部的高4位。為得到這些數(shù)據(jù),，可將帶符號(hào)的16位數(shù)據(jù)右移11次,。在從頻譜的實(shí)部和虛部取得索引號(hào)前，需首先將它們轉(zhuǎn)換為絕對(duì)值,。因此,，符號(hào)位為零。

       從式6我們已經(jīng)知道,，頻譜的幅度是關(guān)于X(N / 2)對(duì)稱的,，因此我們只需將前(N / 2) + 1個(gè)頻譜數(shù)據(jù)轉(zhuǎn)換為極坐標(biāo)形式。還有,，我們可以看到,，對(duì)于實(shí)數(shù)輸入采樣,，X(0)和X(N / 2)的虛部總為零,。因此這兩條譜線的幅度被單獨(dú)計(jì)算。本項(xiàng)目實(shí)際固件的注釋中包含了用于自動(dòng)生成該LUT的源代碼,，可由程序調(diào)用來計(jì)算X(k)的幅度,。

       Hamming或Hann窗

       此項(xiàng)目固件還包括了對(duì)輸入采樣加Hamming或Hann窗的LUT (Q8.7格式)。加窗函數(shù)可有效降低對(duì)時(shí)域采樣x(n)的舍入操作所引起的頻譜泄漏,。Hamming和Hann窗函數(shù)分別如式7和8所示,。

       程序清單6給出了實(shí)現(xiàn)這些函數(shù)的代碼。同樣,，本項(xiàng)目實(shí)際固件的注釋中包含了用于自動(dòng)生成這些LUT的源代碼,，可由程序調(diào)用來實(shí)現(xiàn)這些窗函數(shù)。

    清單6. 用來實(shí)現(xiàn)Hamming和Hann窗函數(shù)的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
}

       測(cè)試結(jié)果

       為了測(cè)試該FFT應(yīng)用的性能,，固件將X(k)幅度通過μC的UART端口上傳給PC。專門編寫的FFT Graph軟件(隨該項(xiàng)目固件一起提供)用于從PC串口讀取這些幅值,，并以圖形方式實(shí)時(shí)顯示頻譜,。圖3顯示了μC以200ksps采樣四種不同輸入信號(hào)并處理后，由FFT Graph所顯示出來的結(jié)果：

4.3V直流信號(hào)
50kHz正弦信號(hào)
70kHz正弦信號(hào)
6.25kHz方波

                             圖3. FFT Graph軟件顯示的由低功耗μC計(jì)算出的頻譜,。

       接下來干什么,？

       有興趣的讀者還可以花費(fèi)大量的時(shí)間來繼續(xù)優(yōu)化和重新配置該FFT應(yīng)用。盡管在本文中我們選擇了radix-2算法,，還有很多其他算法可以顯著降低加法和乘法運(yùn)算量,。很多本文所未提及的優(yōu)化可以提升FFT的速度。例如,，作為純實(shí)數(shù)的輸入采樣,，其虛部總為零,，頻譜中只有前半部分有實(shí)際意義。利用這一點(diǎn),，第一級(jí)和最后一級(jí)FFT的執(zhí)行速度可進(jìn)一步優(yōu)化,，但需要付出更多的程序空間。

       總之,，本文所討論的算法對(duì)于低功耗μC上的FFT應(yīng)用而言,，提供了一個(gè)很好的出發(fā)點(diǎn)。如果想了解更多信息和具體實(shí)現(xiàn)的細(xì)節(jié),，請(qǐng)查閱我們?yōu)楸緫?yīng)用所提供的,、帶有詳細(xì)注釋的固件信息。