引 言

　　基于FFT的頻譜分析方法可以從含有噪聲的信號中提取有用的信息,，在儀器儀表的數(shù)據(jù)處理中具有重要的應用價值,。為了保證頻譜分析的精度，往往進行多點的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)存,。測試結果表明,，我們編制的程序在MSP430F、1611單片機上,，完成一次2048點的基于FFT的頻譜分析和校正只需要500 ms,，精度也達到要求，可以用于以低功耗單片機為核心的儀表中,，實時完成信號處理任務,。

　　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ù)呢?這其中的關鍵就是由程序員來確定一個數(shù)的小數(shù)點處于16位中的哪一位,。這就是數(shù)的定標。

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

　　從表1中可以看出，同樣一個16位數(shù),，若小數(shù)點設定的位置不同,，它所表示的數(shù)也就不同。例如,，十六進制數(shù)2000H=8192，用Q0表示,；十六進制數(shù)2000H=O．25,，用Q15表示；但對于定點運算來說,，處理方法是完全相同的,。下面簡要介紹如何使用定點數(shù)乘法運算模擬浮點數(shù)乘法。

　　設浮點乘法運算的表達式為：float x,，y,，z；z=xy,。假設經(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位采樣結果寄存器中的數(shù)據(jù)右移4位保存起來,，維持12位精度,，轉換為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計算過程中的這一級所有結果都右移2位,，就能夠使得這一級計算結果的最大值仍然可用13位表示，同時也可將這一級所有蝶形運算輸出的數(shù)據(jù)同時縮小,，保證下級計算,。表示數(shù)值大小的15位數(shù)據(jù)的數(shù)據(jù)格式如圖2所示,。

　　1．4 旋轉因子數(shù)據(jù)定標

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

　　2 防止溢出，保證精度

　　FFT中的蝶形運算如圖3所示,。設蝶形輸入為：X1(k),，實部為X1(k)r，虛部為X1(k)i,；X2(k)，實部為X2(k)r,，虛部為X2(k)i,。設蝶形輸出為：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倍,，因此,，蝶形輸出的實部和虛部結果都需要右移2位(縮小4倍)以防止溢出。但隨著計算級數(shù)的增加,，移位將會使數(shù)據(jù)變得越來越小,。例如，128點FFT,，總共需要7級運算,，數(shù)據(jù)最終將移位2×7=14位(縮小47=16 384倍)，因此當信號幅值不夠大時,，經(jīng)過多級運算可能會無法分辨出主信號頻率,。

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

　　設經(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時,，就分辨不出主信號頻率。

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

　　3 內(nèi)存分配

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

　?、賹嵅?、虛部數(shù)據(jù)平方計算。由于MSP430F1611內(nèi)部集成了硬件乘法器,，因此可將乘法器的第一操作數(shù)寄存器(OP1),、第二操作數(shù)寄存器(OP2)寫入相同的數(shù)據(jù)實現(xiàn)平方運算,。

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

　?、酃β首V計算結果保存。實部平方結果,、虛部平方結果相加后再存人原來的實部單元,。

　　經(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點,。因為實部、虛部結果都需4 096 KB,，故共需8．192 KB,，正好小于10KB；而運算4 096點FFT時,，共需16．384 KB,，超出10 KB。

　　4 程序實現(xiàn)

　　算法實現(xiàn)時使用如下方法簡化了程序運算過程：

　?、貱程序調用匯編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ù)組太長，初始化需要很長時間,，導致程序還沒有進入main函數(shù),，看門狗就已經(jīng)復位。因此定義全局數(shù)組時,，加上_no_init關鍵字,。例如，定義一個數(shù)據(jù)長度為2 048的不需要初始化的整型數(shù)組,，使用語句no_init int fft[2048],。

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

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

　　5 算法測試

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

　　(1)測試算法運行時間

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

　　每次高低電平翻轉的時間間隔約為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的頻譜分析方法將應用于數(shù)字渦街流量計，在此,，針對氣體40口徑頻率范圍為69～1 380 Hz,，設定采樣頻率為3 717．472 199 Hz，數(shù)據(jù)點數(shù)為2 048,，選擇不同的k值得到不同的頻率信號,，由信號發(fā)生器產(chǎn)生幅值為60 mV的這些信號，送入兩線制渦街流量計信號處理系統(tǒng)低通濾波器前端,，然后經(jīng)過預放大電路和低通濾波電路后,，送入MSP430F1611進行頻率估計。