郭 剛,李 釗,李建軍,胡禮勇,董巍巍
(1.第二炮兵青州士官學(xué)校,山東青州262500;2.第二炮兵駐石家莊地區(qū)軍事代表室,河北石家莊050081)
目前,以單片機(jī)為核心的手持式振動(dòng)信號(hào)采集器配有簡(jiǎn)潔的信號(hào)調(diào)理和顯示電路,能夠完成低頻振動(dòng)信號(hào)的采集、儲(chǔ)存和傳輸,成本低廉,在機(jī)器振動(dòng)狀況檢測(cè)領(lǐng)域應(yīng)用廣泛。但受單片機(jī)硬件資源的制約,其一般不具備現(xiàn)場(chǎng)對(duì)采樣數(shù)據(jù)進(jìn)行FFT的能力[1,2]。然而,FFT是振動(dòng)信號(hào)分析最常用的方法,對(duì)許多現(xiàn)場(chǎng)采集的信號(hào)數(shù)據(jù),只需進(jìn)行少量采樣點(diǎn)的FFT,即可算出機(jī)器振動(dòng)的主頻,進(jìn)而快速判斷運(yùn)轉(zhuǎn)狀態(tài),實(shí)現(xiàn)對(duì)機(jī)器運(yùn)行狀態(tài)的快速檢測(cè)和故障初診,具有很大的實(shí)用價(jià)值。下面將著重闡述在以ATmega64L單片機(jī)為核心的手持式振動(dòng)信號(hào)采集器上實(shí)現(xiàn)128點(diǎn)FFT的過(guò)程和方法。
離散傅里葉變換(DFT)的本質(zhì)就是建立了以時(shí)間為自變量的信號(hào)與以頻率為自變量的頻譜函數(shù)之間的變換關(guān)系。設(shè)X(n)為N點(diǎn)有限長(zhǎng)序列,其DFT公式為:
式中,
按時(shí)間抽取的基2FFT算法,利用旋轉(zhuǎn)因子的對(duì)稱性、周期性和可約性,根據(jù)輸入序列時(shí)間上的奇偶,將一個(gè)N點(diǎn)的DFT分解成2個(gè)N/2點(diǎn)的DFT,再把2個(gè)N/2點(diǎn)的DFT又可按照奇偶分成2個(gè)N/4的DFT,直到最終分解為N/2個(gè)2點(diǎn)的DFT運(yùn)算。本算法的基本運(yùn)算單元是2點(diǎn)蝶形運(yùn)算結(jié)構(gòu)。當(dāng)N=2M時(shí)(這里M為自然數(shù)),整個(gè)運(yùn)算共有M級(jí),每級(jí)又有N/2個(gè)蝶形運(yùn)算組成。蝶形運(yùn)算計(jì)算流程如圖1所示。
圖1 蝶形運(yùn)算單元
蝶形運(yùn)算方程式為:
在處理采集到的實(shí)際數(shù)據(jù)時(shí),輸入信號(hào)序列通常是實(shí)數(shù)。理論上使用N點(diǎn)復(fù)數(shù)FFT完成2N點(diǎn)實(shí)數(shù)序列FFT運(yùn)算,可以提高運(yùn)算速度并且減小存儲(chǔ)空間。但是變換完成后的結(jié)果轉(zhuǎn)化計(jì)算量較大,導(dǎo)致單片機(jī)整體運(yùn)算速度不能顯著提高,又鑒于ATmega 64L單片機(jī)自帶 4 KB的 RAM,完全滿足128個(gè)32位浮點(diǎn)數(shù)進(jìn)行FFT的存儲(chǔ)空間需求,所以此處直接使用128點(diǎn)復(fù)數(shù)FFT算法完成128個(gè)實(shí)數(shù)采樣數(shù)據(jù)的FFT運(yùn)算,在運(yùn)算時(shí)要把原實(shí)數(shù)數(shù)據(jù)作為復(fù)數(shù)的實(shí)部,而對(duì)應(yīng)復(fù)數(shù)虛部需全部置零。
該系統(tǒng)采用ATmega64L單片機(jī)作為控制核心。它是基于RISC結(jié)構(gòu)的8位低功耗CMOS微處理器,大部分指令都可以在一個(gè)時(shí)鐘周期內(nèi)完成,最大可以支持16 MHz的系統(tǒng)時(shí)鐘,片載64 KB的可編程Flash,2B的EEPROM,4KB的SRAM,53個(gè)通用I/O口線,32個(gè)通用工作寄存器和2個(gè)串口、SPI接口和TWI接口等豐富的外圍接口電路。整個(gè)振動(dòng)信號(hào)采集系統(tǒng)的硬件框圖如圖2所示。
圖2 采集系統(tǒng)硬件
使用ATmega64L單片機(jī)進(jìn)行128點(diǎn)采樣數(shù)據(jù)FFT計(jì)算時(shí),應(yīng)首先把輸入的128個(gè)實(shí)信號(hào)虛部置零,使其變?yōu)閺?fù)數(shù),然后將數(shù)據(jù)序列按碼位倒置的方式反序排列,再進(jìn)行按時(shí)間抽取基2FFT運(yùn)算,最后計(jì)算轉(zhuǎn)換后各復(fù)數(shù)的模并比較分析振動(dòng)信號(hào)主頻,具體FFT計(jì)算流程圖如圖3所示。
圖3 128點(diǎn) FFT流程
3.2.1 反序排列算法的實(shí)現(xiàn)
反序排列是由FFT算法的性質(zhì)決定的。首先將128個(gè)數(shù)據(jù)點(diǎn)自然排列的十進(jìn)制數(shù)轉(zhuǎn)換成7位二進(jìn)制數(shù),再將這些二進(jìn)制數(shù)的首位至末尾的順序顛倒過(guò)來(lái)重新?lián)Q算成十進(jìn)制數(shù),即得到反序排列的結(jié)果。為減小單片機(jī)的運(yùn)算負(fù)擔(dān),同時(shí)充分利用ATmega64L程序存儲(chǔ)器大的優(yōu)勢(shì),使用MATLAB軟件算出128點(diǎn)反序排列結(jié)果,并以數(shù)組形式存放在單片機(jī)程序存儲(chǔ)器中,FFT運(yùn)算時(shí),單片機(jī)只需進(jìn)行簡(jiǎn)單的查表,即可快速完成數(shù)據(jù)反序操作。單片機(jī)存儲(chǔ)128個(gè)序列碼為:
const unsigned char inv_tab[128]={0,64,32,96,16,80,……,47,111,31,95,63,127}。
3.2.2 旋轉(zhuǎn)因子的計(jì)算和儲(chǔ)存
進(jìn)行FFT的子函數(shù)運(yùn)算時(shí),實(shí)時(shí)計(jì)算旋轉(zhuǎn)因子的值,對(duì)單片機(jī)而言是十分繁重的任務(wù)。根據(jù)旋轉(zhuǎn)因子的復(fù)數(shù)形式,預(yù)先算出 cos(2πkn/N)和sin(2πkn/N)的值(取小數(shù)點(diǎn)后4位精度),并以數(shù)組的形式存儲(chǔ)在單片機(jī)的程序存儲(chǔ)器內(nèi),供FFT運(yùn)算過(guò)程中隨時(shí)調(diào)用。單片機(jī)存儲(chǔ)旋轉(zhuǎn)因子實(shí)部和虛部的形式為:
const float cos _tab[64]={0,-0.0491,-0.0980,-0.1467,…,-0.9808,-0.9892,-0.9952,-0.9988};
const float sin _tab[64]={1,0.9988,0.9952,…,-0.1467,-0.0980,-0.049}。
3.2.3 FFT各級(jí)蝶形運(yùn)算的實(shí)現(xiàn)
根據(jù)FFT算法的特性,使用3級(jí)循環(huán)完成各級(jí)蝶形運(yùn)算[5]。
for(L=1;L<=7;L++){//L控制蝶形運(yùn)算的級(jí)數(shù)
b=1;i=L-1;while(i>0)
{b=b*2;i--;}//b控制蝶形運(yùn)算2點(diǎn)之間的距離
for(j=0;j<=b-1;j++)
{p=1;i=7-L;
while(i>0)/*p=pow(2,7-L)*j;
{p=p*2;i--;}
p=p*j;//計(jì)算旋轉(zhuǎn)因子虛實(shí)部在數(shù)組中的對(duì)應(yīng)序號(hào)p
for(k=j;k<128;k=k+2*b)//完成本級(jí)N/2個(gè)蝶形運(yùn)算
{TR=dataR[k];TI=dataI[k];temp=dataR[k+b];
dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
dataI[k+b]=TI+temp*sin_ tab[p]-dataI[k+b]*cos_ tab[p];}
}
}
3.2.4 FFT結(jié)果的處理[6]
根據(jù)采樣定理,只要采樣頻率大于信號(hào)頻率的2倍,采樣得到的數(shù)字信號(hào)就可以做FFT變換了。N個(gè)采樣點(diǎn)經(jīng)過(guò)FFT之后就可得到N個(gè)點(diǎn)的FFT結(jié)果。為了方便進(jìn)行FFT運(yùn)算,通常N取2的整數(shù)次方。假設(shè)采樣頻率為fs,信號(hào)頻率為f,采樣點(diǎn)數(shù)為N。那么FFT之后結(jié)果就是一個(gè)為N點(diǎn)的復(fù)數(shù),每個(gè)點(diǎn)就對(duì)應(yīng)有一個(gè)頻率點(diǎn)。這個(gè)點(diǎn)的模值就是該頻率值下的幅度,而每個(gè)點(diǎn)的相位就是該頻率下的信號(hào)的相位。FFT運(yùn)算后第一個(gè)點(diǎn)表示直流分量(即0 Hz),第n個(gè)點(diǎn)所表示的頻率為:fn=(n-1)*fs/N。
為了驗(yàn)證程序的正確性,使用手持式振動(dòng)信號(hào)采集系統(tǒng)對(duì)FG1617函數(shù)波形信號(hào)發(fā)生器產(chǎn)生的標(biāo)準(zhǔn)9 Hz矩形波進(jìn)行數(shù)據(jù)采集,采樣頻率設(shè)為128 Hz(該系統(tǒng)可以在10~1 000 Hz之間靈活設(shè)置采樣頻率)。然后截取采集到的128個(gè)連續(xù)數(shù)據(jù)點(diǎn),使用Lab VIEW編制的專用FFT處理程序和ATmega64L內(nèi)嵌的FFT程序分別進(jìn)行分析,得到的幅頻結(jié)果如表1所示。
表1 FFT結(jié)果比較表
通過(guò)比較分析發(fā)現(xiàn),ATmega64L單片機(jī)對(duì)128點(diǎn)采樣數(shù)據(jù)進(jìn)行FFT,計(jì)算所得到的諧波主頻和幅值與使用專用數(shù)據(jù)處理程序計(jì)算得到的主頻和幅值基本吻合,達(dá)到設(shè)計(jì)要求精度。
FFT是信號(hào)采集領(lǐng)域的重要數(shù)據(jù)處理工具。上述充分利用高性能低成本的ATmega64L單片機(jī)硬件資源和高速運(yùn)算能力,結(jié)合查表手段,可以在220 ms時(shí)間內(nèi)完成128點(diǎn)的32位精度浮點(diǎn)FFT運(yùn)算,而且經(jīng)過(guò)實(shí)驗(yàn)驗(yàn)證轉(zhuǎn)換結(jié)果正確、精度夠用,證明了該型號(hào)單片機(jī)實(shí)現(xiàn)FFT運(yùn)算的可行性,為FFT實(shí)現(xiàn)方法開(kāi)辟了新的應(yīng)用領(lǐng)域。在低成本手持振動(dòng)信號(hào)采集儀器上集成FFT功能,通過(guò)對(duì)現(xiàn)場(chǎng)采集的振動(dòng)信號(hào)進(jìn)行實(shí)時(shí)FFT運(yùn)算,算出被檢測(cè)機(jī)器的振動(dòng)主頻和幅值,實(shí)現(xiàn)了對(duì)采樣數(shù)據(jù)的現(xiàn)場(chǎng)分析,滿足了設(shè)備運(yùn)行狀態(tài)檢測(cè)和故障初判的要求。
[1]肖 飛.手持式振動(dòng)測(cè)試系統(tǒng)的軟件開(kāi)發(fā)[D].南京:東南大學(xué),2005:1-5.
[2]賈民平,劉玉春,鐘秉林,等.便攜式數(shù)據(jù)采集與工況監(jiān)測(cè)分析系統(tǒng)的研制[J].東南大學(xué)學(xué)大報(bào),1997,27(2):99-102.
[3]胡廣書.數(shù)字信號(hào)處理—理論、算法與實(shí)現(xiàn)[M].北京:清華大學(xué)出版社,2003:177-198.
[4]文其林,白曉東,周 洪,等.2048點(diǎn)FFT在TMS320C240x定點(diǎn)DSP上的實(shí)現(xiàn)[J].微計(jì)算機(jī)信息,2006(5-2):159-160.
[5]肖宛昂.嵌人式系統(tǒng)中FFT算法研究[J].單片機(jī)與嵌入式系統(tǒng)應(yīng)用,2003(1):69-71.
[6]李全利,劉長(zhǎng)亮.CCS上FFT運(yùn)算的實(shí)現(xiàn)[J].自動(dòng)化技術(shù)____與應(yīng)用,2009,28(2):61-63.