在MSP430F1611上 實(shí)現(xiàn)周期圖譜分析及校正

摘要 基于FFT的周期圖譜分析方法可以有效地從含有噪聲的信號(hào)中提取有用信息；但是，由于低功耗單片機(jī)的速度和內(nèi)存有限，所以無(wú)法實(shí)時(shí)地完成FFT運(yùn)算。為此，我們采用匯編語(yǔ)言編制FFT程序。采用定點(diǎn)運(yùn)算來(lái)減少運(yùn)算量和存儲(chǔ)量；采用先判斷再移位的方式，既防止了數(shù)據(jù)溢出，又保證了足夠的計(jì)算精度，從而在MSP430F1611上實(shí)現(xiàn)了2048點(diǎn)FFT。其精度和實(shí)時(shí)性滿足儀表的要求，并成功應(yīng)用于低功耗、兩線制數(shù)字渦街流量計(jì)。
關(guān)鍵詞 功率譜估計(jì) 定點(diǎn)計(jì)算 MSP430F1611

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

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

1．2 數(shù)據(jù)定標(biāo)
    定點(diǎn)單片機(jī)參與數(shù)值運(yùn)算的數(shù)都是16位的整型數(shù)，但是運(yùn)算過(guò)程中的數(shù)不一定都是整數(shù)。那么，定點(diǎn)計(jì)算過(guò)程中如何處理小數(shù)呢?這其中的關(guān)鍵就是由程序員來(lái)確定一個(gè)數(shù)的小數(shù)點(diǎn)處于16位中的哪一位。這就是數(shù)的定標(biāo)。
    通過(guò)設(shè)定小數(shù)點(diǎn)在16位數(shù)中的不同位置，就可以表示不同大小和不同精度的小數(shù)了。數(shù)的定標(biāo)有Q表示法和S表示法兩種。表1列出了一個(gè)16位數(shù)的16種Q表示、S表示及它們所能表示的十進(jìn)制數(shù)值范圍。
    從表1中可以看出，同樣一個(gè)16位數(shù)，若小數(shù)點(diǎn)設(shè)定的位置不同，它所表示的數(shù)也就不同。例如，十六進(jìn)制數(shù)2000H=8192，用Q0表示；十六進(jìn)制數(shù)2000H=O．25，用Q15表示；但對(duì)于定點(diǎn)運(yùn)算來(lái)說(shuō)，處理方法是完全相同的。下面簡(jiǎn)要介紹如何使用定點(diǎn)數(shù)乘法運(yùn)算模擬浮點(diǎn)數(shù)乘法。

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

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

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

2 防止溢出，保證精度
     FFT中的蝶形運(yùn)算如圖3所示。設(shè)蝶形輸入為：X1(k)，實(shí)部為X1(k)r，虛部為X1(k)i；X2(k)，實(shí)部為X2(k)r，虛部為X2(k)i。設(shè)蝶形輸出為：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è)數(shù)相加得到的，因此數(shù)據(jù)可能會(huì)放大3倍。如果計(jì)算過(guò)程中的數(shù)據(jù)始終使用定點(diǎn)數(shù)表示，隨著級(jí)數(shù)的增加，就會(huì)發(fā)生溢出。例如，使用16位定點(diǎn)數(shù)表示，其最高位(從左數(shù)第1位)為符號(hào)位，其表示的數(shù)據(jù)范圍為-32 768～32 767。如果采樣得到的數(shù)據(jù)最大值為4 096，經(jīng)過(guò)兩級(jí)計(jì)算后蝶形最大輸出就可能為4096×3×3=36 864，超出了16位定點(diǎn)數(shù)的表示范圍。
    下面介紹保證數(shù)據(jù)計(jì)算精度的方法。
    為了提高計(jì)算速度，系統(tǒng)中使用定點(diǎn)數(shù)法運(yùn)算FFT，且使用Q13表示數(shù)據(jù)。蝶形運(yùn)算中，其蝶形輸出的數(shù)據(jù)的實(shí)部和虛部都使用3次加法運(yùn)算，即每級(jí)蝶形運(yùn)算都可能使數(shù)據(jù)擴(kuò)大3倍，因此，蝶形輸出的實(shí)部和虛部結(jié)果都需要右移2位(縮小4倍)以防止溢出。但隨著計(jì)算級(jí)數(shù)的增加，移位將會(huì)使數(shù)據(jù)變得越來(lái)越小。例如，128點(diǎn)FFT，總共需要7級(jí)運(yùn)算，數(shù)據(jù)最終將移位2×7=14位(縮小47=16 384倍)，因此當(dāng)信號(hào)幅值不夠大時(shí)，經(jīng)過(guò)多級(jí)運(yùn)算可能會(huì)無(wú)法分辨出主信號(hào)頻率。
    設(shè)FFT運(yùn)算結(jié)果的主信號(hào)頻率點(diǎn)的對(duì)應(yīng)實(shí)部為r，虛部為i，其幅值為A(ADC的量化值)，參與運(yùn)算的數(shù)據(jù)點(diǎn)數(shù)為N，由FFT功率譜計(jì)算的性質(zhì)可得：

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

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

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

   
    因此，當(dāng)K2為2，N為128時(shí)，A=128×1．414=180．992=181，即當(dāng)信號(hào)的幅值為18l／4 096×2．538=112 mV，就分辨不出主信號(hào)頻率?？紤]K2為2的極限情況，當(dāng)A為724，N為512時(shí)，即給定信號(hào)幅值為724／4 096×2．538=449 mV時(shí)，就分辨不出主信號(hào)頻率。
    為了防止計(jì)算結(jié)果經(jīng)過(guò)多次移位后，數(shù)據(jù)太小無(wú)法分辨主信號(hào)，系統(tǒng)針對(duì)定點(diǎn)FFT運(yùn)算采取如下處理：由于FFT定點(diǎn)運(yùn)算中，一般情況下，為了處理方便，每級(jí)蝶形運(yùn)算中乘法結(jié)果都限制在-1～1范圍內(nèi)，即乘法運(yùn)算的結(jié)果始終為小數(shù)(只有經(jīng)過(guò)加法運(yùn)算，數(shù)據(jù)才有可能超出-1～1范圍)，因此，通過(guò)判斷蝶形輸出的結(jié)果，決定是否移位。當(dāng)發(fā)現(xiàn)超出-1～1范圍，就將本級(jí)的所有蝶形運(yùn)算的輸出結(jié)果右移2位，沒(méi)有超出就不進(jìn)行移位。

3 內(nèi)存分配
    由式(3)可知，功率譜估算時(shí)需要另外開辟一段內(nèi)存空間存儲(chǔ)功率譜結(jié)果。例如，當(dāng)進(jìn)行2048點(diǎn)基于FFT的功率譜分析時(shí)，需用1024個(gè)浮點(diǎn)數(shù)存放功率譜計(jì)算結(jié)果，這將占有很大一段內(nèi)存。但實(shí)際運(yùn)算中，每個(gè)頻率點(diǎn)功率，只與其FFT運(yùn)算結(jié)果中的對(duì)應(yīng)頻率點(diǎn)的實(shí)部、虛部有關(guān)，而與其他頻率點(diǎn)無(wú)關(guān)。因此功率譜運(yùn)算中，可采取以下步驟將存放實(shí)部的空間存放功率譜：
    ①實(shí)部、虛部數(shù)據(jù)平方計(jì)算。由于MSP430F1611內(nèi)部集成了硬件乘法器，因此可將乘法器的第一操作數(shù)寄存器(OP1)、第二操作數(shù)寄存器(OP2)寫入相同的數(shù)據(jù)實(shí)現(xiàn)平方運(yùn)算。
    ②平方結(jié)果移位。平方結(jié)果需要右移13位，使用Q13表示，同時(shí)使用16位的臨時(shí)變量將平方結(jié)果保存。
    ③功率譜計(jì)算結(jié)果保存。實(shí)部平方結(jié)果、虛部平方結(jié)果相加后再存人原來(lái)的實(shí)部單元。
    經(jīng)過(guò)上述步驟后，就可將原來(lái)存放實(shí)部、虛部數(shù)據(jù)的內(nèi)存單元再次利用。
    定點(diǎn)FFT運(yùn)算過(guò)程中，還可將用來(lái)存放采集數(shù)據(jù)的內(nèi)存空間，再次用作存放FFT運(yùn)算過(guò)程中的實(shí)部數(shù)據(jù)，另外再開辟同等大小的內(nèi)存空間，存放虛部數(shù)據(jù)。例如，對(duì)于RAM空間為10 KB的MSP430F16ll來(lái)說(shuō)，使用16位定點(diǎn)數(shù)運(yùn)算FFT，最多能夠運(yùn)算2 048點(diǎn)。因?yàn)閷?shí)部、虛部結(jié)果都需4 096 KB，故共需8．192 KB，正好小于10KB；而運(yùn)算4 096點(diǎn)FFT時(shí)，共需16．384 KB，超出10 KB。

4 程序?qū)崿F(xiàn)
    算法實(shí)現(xiàn)時(shí)使用如下方法簡(jiǎn)化了程序運(yùn)算過(guò)程：
    ①C程序調(diào)用匯編FFT程序，同時(shí)為了處理方便將功率譜運(yùn)算過(guò)程也用C語(yǔ)言實(shí)現(xiàn)。為了使匯編程序中使用的內(nèi)存空間與C程序中的內(nèi)存空間地址不發(fā)生沖突，匯編程序中所需的變量都在C文件中定義。
    ②由于實(shí)部、虛部都使用C語(yǔ)音數(shù)組來(lái)存儲(chǔ)，當(dāng)計(jì)算點(diǎn)數(shù)很多時(shí)，數(shù)組將很大。例如，當(dāng)運(yùn)算2 048點(diǎn)FFT時(shí)，就需定義兩個(gè)長(zhǎng)度為2 048的整形數(shù)組，這兩段數(shù)組不能用堆棧局部空間存儲(chǔ)，只能用全局?jǐn)?shù)組，由于C語(yǔ)言規(guī)定全局變量默認(rèn)初始化為0，MSP430的IAR編譯環(huán)境，進(jìn)入main函數(shù)之前的cstart函數(shù)中就用cstar_inh_zero函數(shù)對(duì)全局變量進(jìn)行初始化，由于定義的數(shù)組太長(zhǎng)，初始化需要很長(zhǎng)時(shí)間，導(dǎo)致程序還沒(méi)有進(jìn)入main函數(shù)，看門狗就已經(jīng)復(fù)位。因此定義全局?jǐn)?shù)組時(shí)，加上_no_init關(guān)鍵字。例如，定義一個(gè)數(shù)據(jù)長(zhǎng)度為2 048的不需要初始化的整型數(shù)組，使用語(yǔ)句no_init int fft[2048]。
    ③為節(jié)約RAM內(nèi)存空間，將旋轉(zhuǎn)因子對(duì)應(yīng)的正余弦表制作成表格存放在ROM空間中，而蝶形運(yùn)算時(shí)可通過(guò)查表得到旋轉(zhuǎn)因子后再進(jìn)行FFT運(yùn)算。
    ④由于FFT運(yùn)算過(guò)程中的旋轉(zhuǎn)因子是通過(guò)左移14位取整得到的，因此蝶形運(yùn)算過(guò)程中需要乘法運(yùn)算結(jié)果右移14位。MSP430F1611單片機(jī)乘法器的結(jié)果寄存器，由高16位乘法結(jié)果寄存器RESFII、低16位乘法結(jié)果寄存器RESLO組成，右移14位操作就是保留高位結(jié)果寄存器所有內(nèi)容和低位結(jié)果寄存器中的高兩位，因此RESHI、RESLO一起向左移2位，然后保留高位結(jié)果寄存器作為乘法結(jié)果就可實(shí)現(xiàn)右移14位過(guò)程。

5 算法測(cè)試
    為了驗(yàn)證算法的實(shí)時(shí)性和正確性，通過(guò)信號(hào)發(fā)生器產(chǎn)生標(biāo)準(zhǔn)信號(hào)送入所研制的基于MSP430f1611為核心的處理系統(tǒng)，對(duì)算法進(jìn)行了全面的測(cè)試。
    (1)測(cè)試算法運(yùn)行時(shí)間
    測(cè)試對(duì)2 048點(diǎn)數(shù)據(jù)進(jìn)行功率譜估計(jì)所需要的總時(shí)間。預(yù)先設(shè)置MSP430F1611單片機(jī)的P5．4引腳為普通I／O，且為輸出方式，接著，循環(huán)執(zhí)行FFT運(yùn)算和功率譜估計(jì)程序，且每次開始FFT計(jì)算和功率譜估算前，將P5．4輸出電平翻轉(zhuǎn)，因此P5．4輸出的相鄰兩次翻轉(zhuǎn)電平的時(shí)間間隔就是FFT運(yùn)算和功率譜估計(jì)的總時(shí)間。通過(guò)數(shù)字存儲(chǔ)示波器觀測(cè)P5．4引腳輸出的信號(hào)波形，如圖4所示。

    每次高低電平翻轉(zhuǎn)的時(shí)間間隔約為500 ms，即對(duì)2 048點(diǎn)數(shù)據(jù)進(jìn)行功率譜估計(jì)總共需要500 ms。
    (2)測(cè)試算法計(jì)算精度
    由于FFT運(yùn)算的最大誤差發(fā)生在非整周期采樣時(shí)，所以，選擇這些最大誤差點(diǎn)來(lái)進(jìn)行測(cè)試。由于泄漏誤差，信號(hào)基頻表示為
    f0=(k+d)fs／N (7)
    式(7)中，k為整數(shù)，d為小數(shù)(定義d為泄漏誤差系數(shù))。由于泄漏誤差不超過(guò)頻率分辨率的二分之一，所以|d|≤O．5。當(dāng)d=O時(shí)，即為整周期采樣情況；當(dāng)d=O．5時(shí)，就是最大非整周期采樣的地方。因?yàn)樗兄频幕?strong>FFT的頻譜分析方法將應(yīng)用于數(shù)字渦街流量計(jì)，在此，針對(duì)氣體40口徑頻率范圍為69～1 380 Hz，設(shè)定采樣頻率為3 717．472 199 Hz，數(shù)據(jù)點(diǎn)數(shù)為2 048，選擇不同的k值得到不同的頻率信號(hào)，由信號(hào)發(fā)生器產(chǎn)生幅值為60 mV的這些信號(hào)，送入兩線制渦街流量計(jì)信號(hào)處理系統(tǒng)低通濾波器前端，然后經(jīng)過(guò)預(yù)放大電路和低通濾波電路后，送入MSP430F1611進(jìn)行頻率估計(jì)。