|
引言 基于FFT的頻譜分析方法可以從含有噪聲的信號中提取有用的信息，在儀器儀表的數(shù)據(jù)處理中具有重要的應(yīng)用價值。為了保證頻譜分析的精度，往往進行多點的 FFT運算，例如，1024點、2048點等，這樣運算量大、所占內(nèi)存也大，只有采用DSP(數(shù)字信號處理器)才能實現(xiàn)實時的處理。目前，在工業(yè)現(xiàn)場普遍使用的兩線制、低功耗自動化儀表，由于儀表本身消耗電流必須控制在4 mA之內(nèi)，所以無法采用DSP等運算能力強的芯片，只能采用低功耗單片機；而低功耗單片機的運算速度和內(nèi)存容量都很有限，所以，至今未見用其進行多點數(shù) FFT運算的報道。為了能夠用低功耗單片機實時做FFT運算，以提高自動化儀表信息處理的能力，我們用匯編語言編制FFT程序，在程序中用定點數(shù)運算(以下簡稱定點FFT)，采取措施防止數(shù)據(jù)溢出，保證計算精度，合理分配內(nèi)存。測試結(jié)果表明，我們編制的程序在MSP430F、1611單片機上，完成一次 2048點的基于FFT的頻譜分析和校正只需要500 ms，精度也達到要求，可以用于以低功耗單片機為核心的儀表中，實時完成信號處理任務(wù)。 1 定點運算 1．1 數(shù)據(jù)表示 在MSP430中使用C語言實現(xiàn)FFT運算，其乘法和加法運算都是默認使用浮點實現(xiàn)的。于MSP430屬于定點單片機，因此浮點運算必須由大量的定點指令模擬，這將耗費大量的時間。因此我們針對MSP430的特點，使用匯編語言編制FFT程序，在程序中用定點數(shù)運算，并將數(shù)據(jù)統(tǒng)一使用16位定點數(shù)表示。 16位定點數(shù)中最高位(左邊的第1位)作為符號位，剩下的15位用于存放數(shù)值。數(shù)據(jù)格式如圖1所示. 
1．2 數(shù)據(jù)定標 定點單片機參與數(shù)值運算的數(shù)都是16位的整型數(shù)，但是運算過程中的數(shù)不一定都是整數(shù)。那么，定點計算過程中如何處理小數(shù)呢?這其中的關(guān)鍵就是由程序員來確定一個數(shù)的小數(shù)點處于16位中的哪一位。這就是數(shù)的定標。 通過設(shè)定小數(shù)點在16位數(shù)中的不同位置，就可以表示不同大小和不同精度的小數(shù)了。數(shù)的定標有Q表示法和S表示法兩種。表1列出了一個16位數(shù)的16種Q表示、S表示及它們所能表示的十進制數(shù)值范圍。 從表1中可以看出，同樣一個16位數(shù)，若小數(shù)點設(shè)定的位置不同，它所表示的數(shù)也就不同。例如，十六進制數(shù)2000H=8192，用Q0表示；十六進制數(shù) 2000H=O．25，用Q15表示；但對于定點運算來說，處理方法是完全相同的。下面簡要介紹如何使用定點數(shù)乘法運算模擬浮點數(shù)乘法。
設(shè)浮點乘法運算的表達式為：float x，y，z；z=xy。假設(shè)經(jīng)過統(tǒng)計(這里“統(tǒng)計”的意思是所有計算中數(shù)據(jù)范圍都在定標范圍內(nèi))后x的定標值為Qx，y的定標值為Qy，乘積z的定標值為 Qz，則z=xy；zq×2-Qz=xq×yq×2-(Qx+Qy)；zq=(xqyq)2Qz-(Qx+Qy)。所以，定點表示的乘法為：
1．3 FFT計算過程中的數(shù)據(jù)定標 為了在以MSP43F1611為處理器的儀表系統(tǒng)上進行基于FFT的功率譜估計，必須先由MSP430F1611的ADC進行采樣，而ADC采樣得到的數(shù)據(jù)需要經(jīng)過定標后才能進行定點計算。定標過程為：ADC的采樣電壓范圍為0～2．5 V，因此，采樣過程實際上就是將信號電壓除以2．5進行歸一化，使得采樣得到的數(shù)據(jù)范圍為O～1 V，此時數(shù)據(jù)就可用Q15表示，即將ADC的12位采樣結(jié)果寄存器中的數(shù)據(jù)右移4位保存起來，維持12位精度，轉(zhuǎn)換為Q15定點數(shù)表示。由于FFT運算過程中，蝶形輸出相對蝶形輸入數(shù)據(jù)被放大了3倍，因此蝶形輸出數(shù)據(jù)范圍為一3～+3 V。此時數(shù)據(jù)如果仍然使用Q15表示，就會發(fā)生溢出，故改用Q13表示數(shù)據(jù)，即將12位ADC數(shù)據(jù)右移1位。實際上經(jīng)過處理后，ADC數(shù)據(jù)精度沒有變化，但使用Q13表示數(shù)據(jù)比用Q12表示數(shù)據(jù)，其蝶形輸出的數(shù)據(jù)精度高。這是由于定點計算時需要對蝶形輸出右移以防止溢出，而使用Q13表示數(shù)據(jù)比使用Q12 表示數(shù)據(jù)少右移了1位，因此多了1位有效數(shù)據(jù)。FFT運算過程中使用Q13表示數(shù)據(jù)，就使得加法乘法運算都可以直接使用定點指令實現(xiàn)，減少了很多判斷處理，提高了運算速度。使用Q13表示數(shù)據(jù)，即最高位(左邊的第1位)是符號位，剩下的15位表示數(shù)據(jù)的大小。表示數(shù)據(jù)大小的15位中，高2位(左邊的第2 位和第3位)用來表示數(shù)據(jù)中的整數(shù)部分，在計算中作為保護位；最低的13位(右邊的13位)表示數(shù)據(jù)中的小數(shù)部分，如果經(jīng)過某次蝶形單元運算后，最大值正好被放大3倍，此時數(shù)據(jù)就由13位擴大到15位，保證數(shù)據(jù)不會增大到16位，沖走符號位，發(fā)生溢出。運算完成后將FFT計算過程中的這一級所有結(jié)果都右移 2位，就能夠使得這一級計算結(jié)果的最大值仍然可用13位表示，同時也可將這一級所有蝶形運算輸出的數(shù)據(jù)同時縮小，保證下級計算。表示數(shù)值大小的15位數(shù)據(jù)的數(shù)據(jù)格式如圖2所示。
1．4 旋轉(zhuǎn)因子數(shù)據(jù)定標 FFT運算過程使用定點計算，且使用有符號乘法，必須始終保留1位作為符號位；而旋轉(zhuǎn)因子范圍為-1～1，因此可定標為Q14，轉(zhuǎn)換為16位定點數(shù)。其轉(zhuǎn)換過程為：根據(jù)參與FFT運算的數(shù)據(jù)點數(shù)計算出旋轉(zhuǎn)因子的正余弦表，然后將正余弦表乘以16384，即左移14位，最后四舍五人取整。如果使用Q15表示數(shù)據(jù)，即需要左移15位，那么正余弦表中最大值1，經(jīng)過上述處理后成為-1，發(fā)生溢出。 2 防止溢出，保證精度 FFT中的蝶形運算如圖3所示。設(shè)蝶形輸入為：X1(k)，實部為X1(k)r，虛部為X1(k)i；X2(k)，實部為X2(k)r，虛部為 X2(k)i。設(shè)蝶形輸出為：X(k)，實部為X(k)，，虛部為X(k)i；X(k+N／2)，實部為X(k+N／2)r，虛部為X(k+N／2)i。則有：
由式(1)和式(2)可以看出，蝶形輸出的實部和虛部是由3個數(shù)相加得到的，因此數(shù)據(jù)可能會放大3倍。如果計算過程中的數(shù)據(jù)始終使用定點數(shù)表示，隨著級數(shù)的增加，就會發(fā)生溢出。例如，使用16位定點數(shù)表示，其最高位(從左數(shù)第1位)為符號位，其表示的數(shù)據(jù)范圍為-32 768～32 767。如果采樣得到的數(shù)據(jù)最大值為4 096，經(jīng)過兩級計算后蝶形最大輸出就可能為4096×3×3=36 864，超出了16位定點數(shù)的表示范圍。 下面介紹保證數(shù)據(jù)計算精度的方法。 為了提高計算速度，系統(tǒng)中使用定點數(shù)法運算FFT，且使用Q13表示數(shù)據(jù)。蝶形運算中，其蝶形輸出的數(shù)據(jù)的實部和虛部都使用3次加法運算，即每級蝶形運算都可能使數(shù)據(jù)擴大3倍，因此，蝶形輸出的實部和虛部結(jié)果都需要右移2位(縮小4倍)以防止溢出。但隨著計算級數(shù)的增加，移位將會使數(shù)據(jù)變得越來越小。例如，128點FFT，總共需要7級運算，數(shù)據(jù)最終將移位2×7=14位(縮小47=16 384倍)，因此當信號幅值不夠大時，經(jīng)過多級運算可能會無法分辨出主信號頻率。 設(shè)FFT運算結(jié)果的主信號頻率點的對應(yīng)實部為r，虛部為i，其幅值為A(ADC的量化值)，參與運算的數(shù)據(jù)點數(shù)為N，由FFT功率譜計算的性質(zhì)可得：
設(shè)經(jīng)過定點FFT運算，也就是運算過程中有移位，則該主信號頻率點的模為K，即：
聯(lián)立式(3)和式(4)，得
由于功率譜估計是找出功率譜中的最大值，確定主信號的頻率，根據(jù)經(jīng)驗，使用定點數(shù)運算FFT，當實部和虛部的模的平方K2為2時，就無法由功率譜分辨出主信號頻率。由式(5)可得：
因此，當K2為2，N為128時，A=128×1．414=180．992=181，即當信號的幅值為18l／4 096×2．538=112 mV，就分辨不出主信號頻率。考慮K2為2的極限情況，當A為724，N為512時，即給定信號幅值為724／4 096×2．538=449 mV時，就分辨不出主信號頻率。 為了防止計算結(jié)果經(jīng)過多次移位后，數(shù)據(jù)太小無法分辨主信號，系統(tǒng)針對定點FFT運算采取如下處理：由于FFT定點運算中，一般情況下，為了處理方便，每級蝶形運算中乘法結(jié)果都限制在-1～1范圍內(nèi)，即乘法運算的結(jié)果始終為小數(shù)(只有經(jīng)過加法運算，數(shù)據(jù)才有可能超出-1～1范圍)，因此，通過判斷蝶形輸出的結(jié)果，決定是否移位。當發(fā)現(xiàn)超出-1～1范圍，就將本級的所有蝶形運算的輸出結(jié)果右移2位，沒有超出就不進行移位。 3 內(nèi)存分配 由式(3)可知，功率譜估算時需要另外開辟一段內(nèi)存空間存儲功率譜結(jié)果。例如，當進行2048點基于FFT的功率譜分析時，需用1024個浮點數(shù)存放功率譜計算結(jié)果，這將占有很大一段內(nèi)存。但實際運算中，每個頻率點功率，只與其FFT運算結(jié)果中的對應(yīng)頻率點的實部、虛部有關(guān)，而與其他頻率點無關(guān)。因此功率譜運算中，可采取以下步驟將存放實部的空間存放功率譜： ①實部、虛部數(shù)據(jù)平方計算。由于MSP430F1611內(nèi)部集成了硬件乘法器，因此可將乘法器的第一操作數(shù)寄存器(OP1)、第二操作數(shù)寄存器(OP2) 寫入相同的數(shù)據(jù)實現(xiàn)平方運算。 ②平方結(jié)果移位。平方結(jié)果需要右移13位，使用Q13表示，同時使用16位的臨時變量將平方結(jié)果保存。 ③功率譜計算結(jié)果保存。實部平方結(jié)果、虛部平方結(jié)果相加后再存人原來的實部單元。 經(jīng)過上述步驟后，就可將原來存放實部、虛部數(shù)據(jù)的內(nèi)存單元再次利用。 定點FFT運算過程中，還可將用來存放采集數(shù)據(jù)的內(nèi)存空間，再次用作存放FFT運算過程中的實部數(shù)據(jù)，另外再開辟同等大小的內(nèi)存空間，存放虛部數(shù)據(jù)。例如，對于RAM空間為10 KB的MSP430F16ll來說，使用16位定點數(shù)運算FFT，最多能夠運算2 048點。因為實部、虛部結(jié)果都需4 096 KB，故共需8．192 KB，正好小于10KB；而運算4 096點FFT時，共需16．384 KB，超出10 KB。 4 程序?qū)崿F(xiàn) 算法實現(xiàn)時使用如下方法簡化了程序運算過程： ①C程序調(diào)用匯編FFT程序，同時為了處理方便將功率譜運算過程也用C語言實現(xiàn)。為了使匯編程序中使用的內(nèi)存空間與C程序中的內(nèi)存空間地址不發(fā)生沖突，匯編程序中所需的變量都在C文件中定義。 ②由于實部、虛部都使用C語音數(shù)組來存儲，當計算點數(shù)很多時，數(shù)組將很大。例如，當運算2 048點FFT時，就需定義兩個長度為2 048的整形數(shù)組，這兩段數(shù)組不能用堆棧局部空間存儲，只能用全局數(shù)組，由于C語言規(guī)定全局變量默認初始化為0，MSP430的IAR編譯環(huán)境，進入 main函數(shù)之前的cstart函數(shù)中就用cstar_inh_zero函數(shù)對全局變量進行初始化，由于定義的數(shù)組太長，初始化需要很長時間，導(dǎo)致程序還沒有進入main函數(shù)，看門狗就已經(jīng)復(fù)位。因此定義全局數(shù)組時，加上_no_init關(guān)鍵字。例如，定義一個數(shù)據(jù)長度為2 048的不需要初始化的整型數(shù)組，使用語句no_init int fft[2048]。 ③為節(jié)約RAM內(nèi)存空間，將旋轉(zhuǎn)因子對應(yīng)的正余弦表制作成表格存放在ROM空間中，而蝶形運算時可通過查表得到旋轉(zhuǎn)因子后再進行FFT運算。 ④由于FFT運算過程中的旋轉(zhuǎn)因子是通過左移14位取整得到的，因此蝶形運算過程中需要乘法運算結(jié)果右移14位。MSP430F1611單片機乘法器的結(jié)果寄存器，由高16位乘法結(jié)果寄存器RESFII、低16位乘法結(jié)果寄存器RESLO組成，右移14位操作就是保留高位結(jié)果寄存器所有內(nèi)容和低位結(jié)果寄存器中的高兩位，因此RESHI、RESLO一起向左移2位，然后保留高位結(jié)果寄存器作為乘法結(jié)果就可實現(xiàn)右移14位過程。 5 算法測試 為了驗證算法的實時性和正確性，通過信號發(fā)生器產(chǎn)生標準信號送入所研制的基于MSP430f1611為核心的處理系統(tǒng)，對算法進行了全面的測試。 (1)測試算法運行時間 測試對2 048點數(shù)據(jù)進行功率譜估計所需要的總時間。預(yù)先設(shè)置MSP430F1611單片機的P5．4引腳為普通I／O，且為輸出方式，接著，循環(huán)執(zhí)行FFT運算和功率譜估計程序，且每次開始FFT計算和功率譜估算前，將P5．4輸出電平翻轉(zhuǎn)，因此P5．4輸出的相鄰兩次翻轉(zhuǎn)電平的時間間隔就是FFT運算和功率譜估計的總時間。通過數(shù)字存儲示波器觀測P5．4引腳輸出的信號波形，如圖4所示。
每次高低電平翻轉(zhuǎn)的時間間隔約為500 ms，即對2 048點數(shù)據(jù)進行功率譜估計總共需要500 ms。 (2)測試算法計算精度 由于FFT運算的最大誤差發(fā)生在非整周期采樣時，所以，選擇這些最大誤差點來進行測試。由于泄漏誤差，信號基頻表示為 f0=(k+d)fs／N (7) 式(7)中，k為整數(shù)，d為小數(shù)(定義d為泄漏誤差系數(shù))。由于泄漏誤差不超過頻率分辨率的二分之一，所以|d|≤O．5。當d=O時，即為整周期采樣情況；當d=O．5時，就是最大非整周期采樣的地方。因為所研制的基于FFT的頻譜分析方法將應(yīng)用于數(shù)字渦街流量計，在此，針對氣體40口徑頻率范圍為 69～1 380 Hz，設(shè)定采樣頻率為3 717．472 199 Hz，數(shù)據(jù)點數(shù)為2 048，選擇不同的k值得到不同的頻率信號，由信號發(fā)生器產(chǎn)生幅值為60 mV的這些信號，送入兩線制渦街流量計信號處理系統(tǒng)低通濾波器前端，然后經(jīng)過預(yù)放大電路和低通濾波電路后，送入MSP430F1611進行頻率估計。 參考文獻 1. 徐科軍 信號處理技術(shù) 2001 2. 徐科軍.黃云志 定點DSP原理、開發(fā)與應(yīng)用 2002 3. 吳湘淇 信號、系統(tǒng)與信號處理 1996 4. 周利清.蘇菲 數(shù)字信號處理基礎(chǔ) 2005 5. 胡廣書 數(shù)字信號處理 1997 6. 沈建華.楊艷琴.翟驍曙 MSP430 系列16位超低功耗單片機原理與應(yīng)用 2004 7. Texas Instruments Inc MSP430x15x,MSP430x16x,MSP430x-161x Mixed Signal Microcontroller Datasheet 2002 8. Texas Instruments Inc MSP430x1xx Family User's Guide 2006 作者：劉家祥 劉三山 朱志海(合肥工業(yè)大學(xué)) 徐科軍(合肥工業(yè)大學(xué)工業(yè)自動化安徽省工程技術(shù)研究中心) 來源：單片機與嵌入式系統(tǒng)應(yīng)用 2009（5） |
