MSP430F1611周期圖譜校正FFT

引 言

基于FFT的頻譜分析方法可以從含有噪聲的信號中提取有用的信息，在儀器儀表的數據處理中具有重要的應用價(jià)值。為了保證頻譜分析的精度，往往進(jìn)行多點(diǎn)的FFT運算，例如，1024點(diǎn)、2048點(diǎn)等，這樣運算量大、所占內存也大，只有采用DSP(數字信號處理器)才能實(shí)現實(shí)時(shí)的處理。目前，在工業(yè)現場(chǎng)普遍使用的兩線(xiàn)制、低功耗自動(dòng)化儀表，由于儀表本身消耗電流必須控制在4 mA之內，所以無(wú)法采用DSP等運算能力強的芯片，只能采用低功耗單片機;而低功耗單片機的運算速度和內存容量都很有限，所以，至今未見(jiàn)用其進(jìn)行多點(diǎn)數FFT運算的報道。為了能夠用低功耗單片機實(shí)時(shí)做FFT運算，以提高自動(dòng)化儀表信息處理的能力，我們用匯編語(yǔ)言編制FFT程序，在程序中用定點(diǎn)數運算(以下簡(jiǎn)稱(chēng)定點(diǎn)FFT)，采取措施防止數據溢出，保證計算精度，合理分配內存。測試結果表明，我們編制的程序在MSP430F、1611單片機上，完成一次2048點(diǎn)的基于FFT的頻譜分析和校正只需要500 ms，精度也達到要求，可以用于以低功耗單片機為核心的儀表中，實(shí)時(shí)完成信號處理任務(wù)。

1 定點(diǎn)運算

1.1 數據表示

在MSP430中使用C語(yǔ)言實(shí)現FFT運算，其乘法和加法運算都是默認使用浮點(diǎn)實(shí)現的。于MSP430屬于定點(diǎn)單片機，因此浮點(diǎn)運算必須由大量的定點(diǎn)指令模擬，這將耗費大量的時(shí)間。因此我們針對MSP430的特點(diǎn)，使用匯編語(yǔ)言編制FFT程序，在程序中用定點(diǎn)數運算，并將數據統一使用16位定點(diǎn)數表示。16位定點(diǎn)數中最高位(左邊的第1位)作為符號位，剩下的15位用于存放數值。數據格式如圖1所示。

1.2 數據定標

定點(diǎn)單片機參與數值運算的數都是16位的整型數，但是運算過(guò)程中的數不一定都是整數。那么，定點(diǎn)計算過(guò)程中如何處理小數呢?這其中的關(guān)鍵就是由程序員來(lái)確定一個(gè)數的小數點(diǎn)處于16位中的哪一位。這就是數的定標。

通過(guò)設定小數點(diǎn)在16位數中的不同位置，就可以表示不同大小和不同精度的小數了。數的定標有Q表示法和S表示法兩種。表1列出了一個(gè)16位數的16種Q表示、S表示及它們所能表示的十進(jìn)制數值范圍。

從表1中可以看出，同樣一個(gè)16位數，若小數點(diǎn)設定的位置不同，它所表示的數也就不同。例如，十六進(jìn)制數2000H=8192，用Q0表示;十六進(jìn)制數2000H=O.25，用Q15表示;但對于定點(diǎn)運算來(lái)說(shuō)，處理方法是完全相同的。下面簡(jiǎn)要介紹如何使用定點(diǎn)數乘法運算模擬浮點(diǎn)數乘法。

設浮點(diǎn)乘法運算的表達式為：float x，y，z;z=xy。假設經(jīng)過(guò)統計(這里“統計”的意思是所有計算中數據范圍都在定標范圍內)后x的定標值為Qx，y的定標值為Qy，乘積z的定標值為Qz，則z=xy;zq×2-Qz=xq×yq×2-(Qx+Qy);zq=(xqyq)2Qz-(Qx+Qy)。所以，定點(diǎn)表示的乘法為：

1.3 FFT計算過(guò)程中的數據定標

為了在以MSP43F1611為處理器的儀表系統上進(jìn)行基于FFT的功率譜估計，必須先由MSP430F1611的ADC進(jìn)行采樣，而ADC采樣得到的數據需要經(jīng)過(guò)定標后才能進(jìn)行定點(diǎn)計算。定標過(guò)程為：ADC的采樣電壓范圍為0～2.5 V，因此，采樣過(guò)程實(shí)際上就是將信號電壓除以2.5進(jìn)行歸一化，使得采樣得到的數據范圍為O～1 V，此時(shí)數據就可用Q15表示，即將ADC的12位采樣結果寄存器中的數據右移4位保存起來(lái)，維持12位精度，轉換為Q15定點(diǎn)數表示。由于FFT運算過(guò)程中，蝶形輸出相對蝶形輸入數據被放大了3倍，因此蝶形輸出數據范圍為一3～+3 V。此時(shí)數據如果仍然使用Q15表示，就會(huì )發(fā)生溢出，故改用Q13表示數據，即將12位ADC數據右移1位。實(shí)際上經(jīng)過(guò)處理后，ADC數據精度沒(méi)有變化，但使用Q13表示數據比用Q12表示數據，其蝶形輸出的數據精度高。這是由于定點(diǎn)計算時(shí)需要對蝶形輸出右移以防止溢出，而使用Q13表示數據比使用Q12表示數據少右移了1位，因此多了1位有效數據。FFT運算過(guò)程中使用Q13表示數據，就使得加法乘法運算都可以直接使用定點(diǎn)指令實(shí)現，減少了很多判斷處理，提高了運算速度。使用Q13表示數據，即最高位(左邊的第1位)是符號位，剩下的15位表示數據的大小。表示數據大小的15位中，高2位(左邊的第2位和第3位)用來(lái)表示數據中的整數部分，在計算中作為保護位;最低的13位(右邊的13位)表示數據中的小數部分，如果經(jīng)過(guò)某次蝶形單元運算后，最大值正好被放大3倍，此時(shí)數據就由13位擴大到15位，保證數據不會(huì )增大到16位，沖走符號位，發(fā)生溢出。運算完成后將FFT計算過(guò)程中的這一級所有結果都右移2位，就能夠使得這一級計算結果的最大值仍然可用13位表示，同時(shí)也可將這一級所有蝶形運算輸出的數據同時(shí)縮小，保證下級計算。表示數值大小的15位數據的數據格式如圖2所示。

1.4 旋轉因子數據定標

FFT運算過(guò)程使用定點(diǎn)計算，且使用有符號乘法，必須始終保留1位作為符號位;而旋轉因子范圍為-1～1，因此可定標為Q14，轉換為16位定點(diǎn)數。其轉換過(guò)程為：根據參與FFT運算的數據點(diǎn)數計算出旋轉因子的正余弦表，然后將正余弦表乘以16384，即左移14位，最后四舍五人取整。如果使用Q15表示數據，即需要左移15位，那么正余弦表中最大值1，經(jīng)過(guò)上述處理后成為-1，發(fā)生溢出。

2 防止溢出，保證精度

FFT中的蝶形運算如圖3所示。設蝶形輸入為：X1(k)，實(shí)部為X1(k)r，虛部為X1(k)i;X2(k)，實(shí)部為X2(k)r，虛部為X2(k)i。設蝶形輸出為：X(k)，實(shí)部為X(k)，，虛部為X(k)i;X(k+N/2)，實(shí)部為X(k+N/2)r，虛部為X(k+N/2)i。則有：

由式(1)和式(2)可以看出，蝶形輸出的實(shí)部和虛部是由3個(gè)數相加得到的，因此數據可能會(huì )放大3倍。如果計算過(guò)程中的數據始終使用定點(diǎn)數表示，隨著(zhù)級數的增加，就會(huì )發(fā)生溢出。例如，使用16位定點(diǎn)數表示，其最高位(從左數第1位)為符號位，其表示的數據范圍為-32 768～32 767。如果采樣得到的數據最大值為4 096，經(jīng)過(guò)兩級計算后蝶形最大輸出就可能為4096×3×3=36 864，超出了16位定點(diǎn)數的表示范圍。

下面介紹保證數據計算精度的方法。

為了提高計算速度，系統中使用定點(diǎn)數法運算FFT，且使用Q13表示數據。蝶形運算中，其蝶形輸出的數據的實(shí)部和虛部都使用3次加法運算，即每級蝶形運算都可能使數據擴大3倍，因此，蝶形輸出的實(shí)部和虛部結果都需要右移2位(縮小4倍)以防止溢出。但隨著(zhù)計算級數的增加，移位將會(huì )使數據變得越來(lái)越小。例如，128點(diǎn)FFT，總共需要7級運算，數據最終將移位2×7=14位(縮小47=16 384倍)，因此當信號幅值不夠大時(shí)，經(jīng)過(guò)多級運算可能會(huì )無(wú)法分辨出主信號頻率。

設FFT運算結果的主信號頻率點(diǎn)的對應實(shí)部為r，虛部為i，其幅值為A(ADC的量化值)，參與運算的數據點(diǎn)數為N，由FFT功率譜計算的性質(zhì)可得：

設經(jīng)過(guò)定點(diǎn)FFT運算，也就是運算過(guò)程中有移位，則該主信號頻率點(diǎn)的模為K，即：

聯(lián)立式(3)和式(4)，得

由于功率譜估計是找出功率譜中的最大值，確定主信號的頻率，根據經(jīng)驗，使用定點(diǎn)數運算FFT，當實(shí)部和虛部的模的平方K2為2時(shí)，就無(wú)法由功率譜分辨出主信號頻率。由式(5)可得：

因此，當K2為2，N為128時(shí)，A=128×1.414=180.992=181，即當信號的幅值為18l/4 096×2.538=112 mV，就分辨不出主信號頻率?？紤]K2為2的極限情況，當A為724，N為512時(shí)，即給定信號幅值為724/4 096×2.538=449 mV時(shí)，就分辨不出主信號頻率。

為了防止計算結果經(jīng)過(guò)多次移位后，數據太小無(wú)法分辨主信號，系統針對定點(diǎn)FFT運算采取如下處理：由于FFT定點(diǎn)運算中，一般情況下，為了處理方便，每級蝶形運算中乘法結果都限制在-1～1范圍內，即乘法運算的結果始終為小數(只有經(jīng)過(guò)加法運算，數據才有可能超出-1～1范圍)，因此，通過(guò)判斷蝶形輸出的結果，決定是否移位。當發(fā)現超出-1～1范圍，就將本級的所有蝶形運算的輸出結果右移2位，沒(méi)有超出就不進(jìn)行移位。

3 內存分配

由式(3)可知，功率譜估算時(shí)需要另外開(kāi)辟一段內存空間存儲功率譜結果。例如，當進(jìn)行2048點(diǎn)基于FFT的功率譜分析時(shí)，需用1024個(gè)浮點(diǎn)數存放功率譜計算結果，這將占有很大一段內存。但實(shí)際運算中，每個(gè)頻率點(diǎn)功率，只與其FFT運算結果中的對應頻率點(diǎn)的實(shí)部、虛部有關(guān)，而與其他頻率點(diǎn)無(wú)關(guān)。因此功率譜運算中，可采取以下步驟將存放實(shí)部的空間存放功率譜：

①實(shí)部、虛部數據平方計算。由于MSP430F1611內部集成了硬件乘法器，因此可將乘法器的第一操作數寄存器(OP1)、第二操作數寄存器(OP2)寫(xiě)入相同的數據實(shí)現平方運算。

②平方結果移位。平方結果需要右移13位，使用Q13表示，同時(shí)使用16位的臨時(shí)變量將平方結果保存。

③功率譜計算結果保存。實(shí)部平方結果、虛部平方結果相加后再存人原來(lái)的實(shí)部單元。

經(jīng)過(guò)上述步驟后，就可將原來(lái)存放實(shí)部、虛部數據的內存單元再次利用。

定點(diǎn)FFT運算過(guò)程中，還可將用來(lái)存放采集數據的內存空間，再次用作存放FFT運算過(guò)程中的實(shí)部數據，另外再開(kāi)辟同等大小的內存空間，存放虛部數據。例如，對于RAM空間為10 KB的MSP430F16ll來(lái)說(shuō)，使用16位定點(diǎn)數運算FFT，最多能夠運算2 048點(diǎn)。因為實(shí)部、虛部結果都需4 096 KB，故共需8.192 KB，正好小于10KB;而運算4 096點(diǎn)FFT時(shí)，共需16.384 KB，超出10 KB。

4 程序實(shí)現

算法實(shí)現時(shí)使用如下方法簡(jiǎn)化了程序運算過(guò)程：

①C程序調用匯編FFT程序，同時(shí)為了處理方便將功率譜運算過(guò)程也用C語(yǔ)言實(shí)現。為了使匯編程序中使用的內存空間與C程序中的內存空間地址不發(fā)生沖突，匯編程序中所需的變量都在C文件中定義。

②由于實(shí)部、虛部都使用C語(yǔ)音數組來(lái)存儲，當計算點(diǎn)數很多時(shí)，數組將很大。例如，當運算2 048點(diǎn)FFT時(shí)，就需定義兩個(gè)長(cháng)度為2 048的整形數組，這兩段數組不能用堆棧局部空間存儲，只能用全局數組，由于C語(yǔ)言規定全局變量默認初始化為0，MSP430的IAR編譯環(huán)境，進(jìn)入main函數之前的cstart函數中就用cstar_inh_zero函數對全局變量進(jìn)行初始化，由于定義的數組太長(cháng)，初始化需要很長(cháng)時(shí)間，導致程序還沒(méi)有進(jìn)入main函數，看門(mén)狗就已經(jīng)復位。因此定義全局數組時(shí)，加上_no_init關(guān)鍵字。例如，定義一個(gè)數據長(cháng)度為2 048的不需要初始化的整型數組，使用語(yǔ)句no_init int fft[2048]。

③為節約RAM內存空間，將旋轉因子對應的正余弦表制作成表格存放在ROM空間中，而蝶形運算時(shí)可通過(guò)查表得到旋轉因子后再進(jìn)行FFT運算。

④由于FFT運算過(guò)程中的旋轉因子是通過(guò)左移14位取整得到的，因此蝶形運算過(guò)程中需要乘法運算結果右移14位。MSP430F1611單片機乘法器的結果寄存器，由高16位乘法結果寄存器RESFII、低16位乘法結果寄存器RESLO組成，右移14位操作就是保留高位結果寄存器所有內容和低位結果寄存器中的高兩位，因此RESHI、RESLO一起向左移2位，然后保留高位結果寄存器作為乘法結果就可實(shí)現右移14位過(guò)程。

5 算法測試

為了驗證算法的實(shí)時(shí)性和正確性，通過(guò)信號發(fā)生器產(chǎn)生標準信號送入所研制的基于MSP430f1611為核心的處理系統，對算法進(jìn)行了全面的測試。

(1)測試算法運行時(shí)間

測試對2 048點(diǎn)數據進(jìn)行功率譜估計所需要的總時(shí)間。預先設置MSP430F1611單片機的P5.4引腳為普通I/O，且為輸出方式，接著(zhù)，循環(huán)執行FFT運算和功率譜估計程序，且每次開(kāi)始FFT計算和功率譜估算前，將P5.4輸出電平翻轉，因此P5.4輸出的相鄰兩次翻轉電平的時(shí)間間隔就是FFT運算和功率譜估計的總時(shí)間。通過(guò)數字存儲示波器觀(guān)測P5.4引腳輸出的信號波形，如圖4所示。

每次高低電平翻轉的時(shí)間間隔約為500 ms，即對2 048點(diǎn)數據進(jìn)行功率譜估計總共需要500 ms。

(2)測試算法計算精度

由于FFT運算的最大誤差發(fā)生在非整周期采樣時(shí)，所以，選擇這些最大誤差點(diǎn)來(lái)進(jìn)行測試。由于泄漏誤差，信號基頻表示為

f0=(k+d)fs/N (7)

式(7)中，k為整數，d為小數(定義d為泄漏誤差系數)。由于泄漏誤差不超過(guò)頻率分辨率的二分之一，所以|d|≤O.5。當d=O時(shí)，即為整周期采樣情況;當d=O.5時(shí)，就是最大非整周期采樣的地方。因為所研制的基于FFT的頻譜分析方法將應用于數字渦街流量計，在此，針對氣體40口徑頻率范圍為69～1 380 Hz，設定采樣頻率為3 717.472 199 Hz，數據點(diǎn)數為2 048，選擇不同的k值得到不同的頻率信號，由信號發(fā)生器產(chǎn)生幅值為60 mV的這些信號，送入兩線(xiàn)制渦街流量計信號處理系統低通濾波器前端，然后經(jīng)過(guò)預放大電路和低通濾波電路后，送入MSP430F1611進(jìn)行頻率估計。