1 引言
科里奧利質(zhì)量流量計（以下簡稱科氏流量計）可直接高精度地測量流體的質(zhì)量流量且可同時獲取流體密度值，是當(dāng)前發(fā)展最為迅速的流量計之一。
科氏流量計通過測量2個磁電傳感器輸出信號之間的相位差（時間差）來得到流體的質(zhì)量流量。為了能夠從含有噪聲的傳感器輸出信號中準(zhǔn)確地提取出流量信息，國內(nèi)外研究者采用數(shù)字信號處理技術(shù)來計算相位差/時間差。但均是以時不變信號為研究對象的，即認(rèn)為在一批處理數(shù)據(jù)（例如2048點）中，信號的頻率、相位和幅值是不變的。實際上，由于受到管內(nèi)流體的流速、密度和流體脈動等因素變化的影響，科氏流量計信號會隨著時間發(fā)生變化。為此，我們率先建立了時變信號模型，并提出了基于陷波濾波器和滑動Goertzel算法的科氏流量計信號處理算法，克服了基于DFT頻譜分析整周期采樣的影響，在每個采樣點都可輸出計算結(jié)果，并且能跟蹤變化的流量信號[1-2]。但是，算法的收斂過程較長，跟蹤精度不高，且跟蹤不到信號的微小變化，同時用于實際系統(tǒng)時計算量大且會發(fā)生數(shù)值溢出。在此基礎(chǔ)上，文獻(xiàn)[3-4]提出了計及負(fù)頻率影響的DTFT遞推算法，即在計算正弦信號的傅里葉系數(shù)時不忽略負(fù)頻率成分的影響，從而提高計算精度，極大地縮短了收斂過程。但該算法只適合一段時間內(nèi)信號恒定的情況，對于流量變化如批料流或有微小波動的時變信號時就失去了作用。同時，文獻(xiàn)[2-3]沒有給出算法的具體實施方案、DSP芯片的型號、信號處理硬件系統(tǒng)和采樣頻率等信息。
為此，針對時變信號，提出將多抽一濾波、自適應(yīng)格型陷波濾波和負(fù)頻率修正的滑動DTFT遞推算法（以下簡稱SDTFT）有機(jī)地組合起來形成一套完整的科氏流量計信號處理算法，并在DSP系統(tǒng)上實現(xiàn)，實時處理流量計信號，既能用于非時變信號又能跟蹤頻率、相位變化的信號，更加接近實際應(yīng)用情況。且算法計算量較小，不會發(fā)生數(shù)值溢出，便于實時實現(xiàn)。
2 時變信號
根據(jù)科氏流量計的工作原理，在理想情況下，傳感器輸出的信號為標(biāo)準(zhǔn)的正弦信號。但在實際應(yīng)用中，由于受到管內(nèi)流體特性如流速、密度、流體脈動及流場等因素變化的影響，信號的頻率和相位并不是恒定不變的。因此，我們把這種頻率、相位等隨時間發(fā)生變化的信號稱為時變信號。我們通過觀察MicroMotion公司生產(chǎn)的型號為CNG050微彎型科氏質(zhì)量流量傳感器的實際工作過程發(fā)現(xiàn)：
1）當(dāng)流量較穩(wěn)定時，兩路信號的相位隨時間會發(fā)生緩慢而微小的變化，這種變化是隨機(jī)的沒有規(guī)律的，但變化的幅度并不像文獻(xiàn)[1]中建立的信號模型那樣大，一般不超過給定相位的1%。在小流量測量中，如1kg/min的流量對應(yīng)的相位差為0.020左右，此時相位波動不超過0.00020，變化范圍很小，但要提高小流量測量精度這種波動是不能忽略的。
2）當(dāng)流體密度確定后，傳感器信號的頻率也會發(fā)生變化，但變化范圍比相位要小很多，最大不超過流量管振動頻率的萬分之一。同時信號頻率反映流體的密度，因此算法必須適用于測試不同的流體，也需要跟蹤信號的頻率。
為此，改進(jìn)時變信號模型如式（1）和（2），更加逼真地模擬實際科氏流量計信號，便于算法仿真。
式（1）和（2）表示信號的相位、頻率在變化，每一時刻的值是前一時刻的值加一個隨機(jī)數(shù)，其中eΦ（n）和eω（n）為零均值、正態(tài)分布、方差為1且互不相關(guān)的白噪聲，σΦ和σω分別控制eφ（n）和eω（n）的變化幅度，當(dāng)信號變化緩慢時，兩者減小，當(dāng)信號突變時，兩者增大。λΦ和λΦ分別控制Φ（n）和ω（n）的變化幅度，相位變化幅度應(yīng)低于給定相位的1%，頻率變化低于振動頻率的0.01%，這樣比較符合實際情況。
3 算法原理及推導(dǎo)
3.1 多抽一濾波
為了增強(qiáng)對噪聲的抑制，先用16kHz較高的采樣頻率對科氏流量計的輸出信號進(jìn)行采樣，然后用多抽一濾波器進(jìn)行抗混疊濾波和抽取。多抽一濾波器分為兩級[4]，第一級為2抽1，使實際的采樣頻率從16kHz降低到8kHz，主要目的是減少數(shù)據(jù)量。第二級為4抽1，采樣頻率降低為2kHz。同時采用30階FIR低通濾波器，不僅保證線性相位，而且在實際的實現(xiàn)中，可以只對抽取的點進(jìn)行濾波，然后再抽取，這樣可以減少計算量節(jié)省時間。多抽一濾波器的系數(shù)在確定截止頻率后通過計算機(jī)輔助設(shè)計的方法得到。仿真結(jié)果表明該方法由于盡可能多地獲取了原信號的信息，所以比單純用2kHz采樣、濾波所得的效果要好。
3.2 自適應(yīng)格型陷波濾波器
自適應(yīng)陷波濾波器參數(shù)可以根據(jù)信號特征收斂并可估算信號的頻率。采用的格型IIR陷波器[1]是由全極點和全零點兩個格型濾波器級聯(lián)而成，傳遞函數(shù)為：
（3）
為了減少計算負(fù)擔(dān)，通過將零點固定在單位圓上，使得只調(diào)整一個參數(shù)就能達(dá)到自適應(yīng)陷波的目的。將零點固定在單位圓上，即令k1=1，k0在經(jīng)過一段時間自適應(yīng)后收斂到−cosω，ω是信號的歸一化頻率，α決定陷阱的寬度，k0使用Burg算法[1]進(jìn)行自適應(yīng)調(diào)整。
由于科氏流量計流體的密度反映為頻率的變化，需要及時跟蹤流體信號的頻率變化。通過大量仿真研究發(fā)現(xiàn)，通過調(diào)整ρ和λ的終值，適當(dāng)?shù)丶哟笙莶ㄆ飨葳宓膶挾?，便能在保證精度的同時實現(xiàn)對頻率變化的跟蹤，ρ和λ計算公式如式（4）和（5）所示。
（4）
（5）
格型自適應(yīng)陷波濾波器的計算量大大降低，且參數(shù)少調(diào)整方便。調(diào)整ρ和λ的終值以及變化的步長就能方便的跟蹤頻率的變化，同時亦能達(dá)到很高的精度。
3.3 SDTFT遞推算法及測量相位差原理
3.3.1 SDTFT遞推算法
離散時間序列的傅里葉變換（DTFT）為：
（6）
DTFT是從第一個采樣點開始通過不斷增加計算的序列長度來實現(xiàn)指定頻率處傅里葉系數(shù)的計算，如果信號在一段時間內(nèi)恒定不變，這種算法是可行的。但是，無法用于時變信號。時變信號的每個采樣點都包含著相位變化的新信息，DTFT將相位變化的新舊信息全部混淆疊加在一起，對相位的變化根本無法靈敏的反映。因此，我們提出滑動DTFT來處理時變信號。
給所觀測的信號加一個N點的時間窗，矩形窗是最簡單的時間窗，并讓這個時間窗隨著采樣點數(shù)的增加不斷向前滑動，如圖1所示。隨著窗函數(shù)的滑動，在每個采樣點計算N點有限長序列的傅里葉變換即為滑動的或滑動窗的DTFT（SDTFT）。面向時變信號時，計算新采樣點的傅里葉系數(shù)時僅利用的是當(dāng)前采樣點之前的N點（N是可以改變的），更新新的相位信息并摒除舊的相位信息，這樣時間窗隨著新采樣點不斷向前滑動，計算的相位差才能跟蹤上實際相位差的變化。
圖1 N點滑動時間窗
如圖1（a）所示，對于觀察信號x（t），設(shè)在m時刻采樣得到N個采樣數(shù)據(jù)x（0），x（1），…，x（N–1），首次構(gòu)成N點有限長序列，其離散時間傅里葉變換為：
（7）
式中：ω為數(shù)字角頻率，單位為rad，t表示采樣點的序號。
圖1（b）所示在m+1時刻，得到新的采樣點x（N），則該點與之前的N–1點重新構(gòu)成一個N點有限長序列，該序列在處的離散時間傅里葉變換為：
以此遞推，當(dāng)新的采樣點與其之前的N–1個采樣點組成第k個N點時間窗即采樣點序號為（N+k–1）時，該序列在ω處的傅里葉變換如式（9）所示。
式（9）即為SDTFT的遞推算法的遞推公式?？梢?，每采入一點新的信號，雖然需計算N點傅里葉變換，但通過遞推公式并沒有增大計算量，比滑動Goertzel算法計算量小，可以滿足科氏流量計實時性要求。同時，每計算一個采樣點的傅里葉系數(shù)始終是N點疊加的計算結(jié)果，不存在序列不斷疊加溢出的問題，非常利于實際系統(tǒng)的實現(xiàn)。
3.3.2 窗長度N的選取
科氏流量計傳感器信號是正弦信號，具有周期性性質(zhì)。同時，格型濾波器估計頻率精度雖然較高，但仍然與真實頻率有偏差。這樣，應(yīng)用DTFT及SDTFT仿真計算周期信號在有偏頻率下的傅里葉系數(shù)時會出現(xiàn)如圖2的現(xiàn)象。
圖2 計算有偏頻率下的相位差
圖2中點劃線為真實相位差值，實線為兩種方法的估計值。在圖2的上圖中，可以看出由DTFT計算的每個采樣點輸出的相位差會發(fā)生波動，但波動是偏離真實值的，且真實值偏上的部分比偏下的部分幅度小，這樣平均處理后得到的結(jié)果會比真實值偏?。欢鴪D2中下圖所示的SDTFT計算結(jié)果中，每個采樣點輸出的相位差是在真實值上下波動的，雖然比上圖中波動的幅度要大，但是真實值上下波動幅度相當(dāng)，這樣平均處理后會比DTFT更接近真實值，精度更高。出現(xiàn)圖2仿真結(jié)果是因為鑒于處理信號的周期性，我們根據(jù)格型濾波器估計的頻率值選擇SDTFT的窗函數(shù)長度N盡量接近信號周期的整數(shù)倍，同時SDTFT對整周期的要求要遠(yuǎn)遠(yuǎn)低于SDFT[6-7]的要求。因此SDTFT的時間窗函數(shù)長度的選擇要針對處理信號特點選擇。
同時，窗長度N的選取還要根據(jù)實際信號的具體變化靈活選取，如果信號的相位變化比較緩慢，可以將N增長，不僅能跟蹤上變化，同時點數(shù)多可以提高計算精度；如果信號的相位差變化快速，可將N縮短，增加跟蹤速度，但不可避免地犧牲了精度，此時可以通過改變窗函數(shù)的形狀如漢寧窗等來提高計算精度。
3.3.3 SDTFT遞推算法計算相位差
由于科氏質(zhì)量流量計傳感器信號為正弦信號，因此可進(jìn)行計及負(fù)頻率的修正[2]，減小頻譜中負(fù)頻率成分的影響，增加計算相位差的精度，縮短收斂過程。具體推到公式如下：
根據(jù)式（12）得到相位差和時間差后，即可根據(jù)儀表系數(shù)得到瞬時流量和累積流量。本文測試相位差的方法用于時變信號時不僅能跟蹤微小緩慢的變化，而且跟蹤速度和精度均優(yōu)于滑動Goertzal算法。
4 MATLAB仿真結(jié)果
型號為CNG050的科氏質(zhì)量流量傳感器，滿管振動時傳感器信號基頻為188.64Hz，因此仿真時信號頻率采用188.64Hz，著重仿真小流量對應(yīng)的相位差。
根據(jù)時變信號模型產(chǎn)生的相位隨采樣點數(shù)變化的曲線以及生成的時變信號波形如圖3所示。
圖（3）、圖（4）仿真參數(shù)為：
（a）相位變化
（b）時變信號波形
圖3 按照時變信號模型產(chǎn)生的時變信號
圖4 本文方法跟蹤相位差變化效果
圖4所示即相位在[0.009940，0.010140]內(nèi)變化時，本文算法的跟蹤效果圖，可見波動幅度為0.01的1%時，本文方法仍具有很好的跟蹤速度和跟蹤精度，而SGA算法已經(jīng)無法跟蹤此時相位的微小變化了。
圖5所示的是相位在[0.0080，0.0150]內(nèi)變化時，本文算法與SGA算法跟蹤相位變化的效果對比，可以看出，相位變化超過0.010的70%時，SGA才能跟蹤上變化，但本文算法的跟蹤速度和精度都優(yōu)于SGA。
圖5 SDTFT與SGA跟蹤相位比較
[0.010，0.20]內(nèi)幾個相位的仿真結(jié)果數(shù)據(jù)如表1所示，仿真參數(shù)同圖6。從表1中可以看出，本文方法具較高的精度。
5 系統(tǒng)研制
本課題組研制了相應(yīng)的數(shù)字式科氏質(zhì)量流量變送器系統(tǒng)，一方面為激振器提供激振信號，另一方面對兩路傳感器信號進(jìn)行處理，得到流量信息。
5.1 系統(tǒng)硬件組成
系統(tǒng)硬件基于TI公司TMS320F28335浮點DSP，由驅(qū)動模塊、信號調(diào)理采集系模塊、脈沖和4～20mA電流輸出及人機(jī)接口組成。
5.2 系統(tǒng)軟件研制
系統(tǒng)軟件設(shè)計采取模塊化設(shè)計方案，將完成特定功能或者類似功能的子程序組合成功能模塊，主要功能模塊如圖6軟件總體框圖所示。
圖6 系統(tǒng)軟件總體框圖
5.2.1 主要模塊
初始化模塊負(fù)責(zé)系統(tǒng)內(nèi)可編程器件和全局變量等的初始化，包括TMS320F28335初始化、全局變量初始化、算法初始化和LCD初始化等。F28335的初始化[8]主要包括系統(tǒng)配置、用到的28335內(nèi)部集成功能模塊的初始化，如看門狗、多通道緩沖串口McBSP、DMA、系統(tǒng)外部接口XINTF和數(shù)字復(fù)用口等以及外部可編程器件如CodecCS4271芯片。全局變量以及算法初始化主要包括反映系統(tǒng)工作狀態(tài)的狀態(tài)變量和算法計算的初始變量等。LCD初始化則包括控制引腳初始化、寫命令字和初始化顯示內(nèi)容。
中斷模塊由外設(shè)中斷擴(kuò)展控制器來處理所有片內(nèi)外設(shè)和外部引腳中斷的優(yōu)先級以及中斷的響應(yīng)。本系統(tǒng)采用了三個中斷，DMACH1接收中斷、DMACH1發(fā)送中斷、定時器0的周期中斷觸發(fā)。Codec啟動后開始采傳感器的信號并進(jìn)行AD轉(zhuǎn)換，之后送到DSP的McBSP的接收引腳上，通過DMA傳輸?shù)絻?nèi)部RAM中，當(dāng)RAM放滿數(shù)據(jù)后產(chǎn)生DMACH1接收中斷，中斷響應(yīng)后將數(shù)據(jù)放于外部RAM一個長的循環(huán)數(shù)組中供算法調(diào)用，并修改DMA的目的地址，為下次傳輸做準(zhǔn)備；定時器0的周期中斷觸發(fā)用于定時刷新LCD的顯示值。
算法模塊內(nèi)包含了系統(tǒng)所進(jìn)行的大部分?jǐn)?shù)學(xué)運算子程序，主要有信號處理算法、瞬時流量計算和累計流量計算以及溫度補(bǔ)償算法。信號處理算法在主監(jiān)控程序中調(diào)用，先調(diào)用多抽一算法程序進(jìn)行抽取，隨后調(diào)用格型濾波估算頻率，之后再利用負(fù)頻率修正的SDTFT遞推算法計算信號的相位差，平均處理后，再結(jié)合儀表系數(shù)計算當(dāng)前的瞬時流量。還需根據(jù)特定公式對結(jié)果進(jìn)行溫度補(bǔ)償。測量結(jié)果包括瞬時質(zhì)量流量、累積質(zhì)量流量、密度和溫度，由LCD顯示。
主監(jiān)控程序[9]是整個系統(tǒng)的總調(diào)度程序，實現(xiàn)儀表所要求的功能。系統(tǒng)上電復(fù)位后，立即進(jìn)行初始化，完成對系統(tǒng)及各模塊的初始化工作，然后監(jiān)控程序開始依次查詢標(biāo)志寄存器的有關(guān)位，如采樣數(shù)據(jù)標(biāo)志位、頻率跟蹤標(biāo)志位、流量計算標(biāo)志位、通信標(biāo)志位等，如果有標(biāo)志位被置位，則調(diào)用相應(yīng)的子程序，并清除標(biāo)志位為下一次操作做好準(zhǔn)備。一次查詢完后，監(jiān)控程序進(jìn)入等待狀態(tài)，直到有關(guān)中斷發(fā)生時開始進(jìn)入下一次查詢，周期循環(huán)，實現(xiàn)儀表流量的實時檢測。
5.2.2 關(guān)鍵實現(xiàn)技術(shù)
1）TMS320F28335是TI公司推出的一款新的浮點運算DSP芯片，由于負(fù)頻率修正的SDTFT算法計算小相位差的關(guān)鍵在于保證公式（20）中分子得到的數(shù)據(jù)的有效位數(shù)，經(jīng)過研究發(fā)現(xiàn)用32位float型變量類型在小流量時只能保證3位有效數(shù)字，使得計算精度降低，因此我們使用64位longdouble型變量類型；
2）本文提出的一套數(shù)字信號處理算法有大量的變量，且變量采用64位的數(shù)據(jù)類型，這就需要較大的存儲空間，因此外擴(kuò)了64k的SARAM，為了提高計算速率，在存儲器映射設(shè)置時將全局變量全部映射到外部存儲器中，將使用頻率高的局部變量和程序映射到內(nèi)部RAM中，使得算法成功實現(xiàn)；
3）經(jīng)過測試F28335定點CPU計算外部RAM中一點數(shù)據(jù)的兩級多抽一算法時間為116μs，計算格型濾波算法的時間為83μs，計算SDTFT遞推算法時間為221μs，再加上對結(jié)果的補(bǔ)償及平均處理，得到一點流量的時間不超過500μs，因此CodecAD采樣頻率設(shè)置為16kHz，多抽一后為2kHz，在保證精度的情況下完全實現(xiàn)了連續(xù)采樣數(shù)據(jù)的實時計算；
4）AD設(shè)置較高的采樣頻率，這就要求傳送大量的數(shù)據(jù)，一點一點傳送顯然不滿足實時性的要求，為此我們利用了TMS320F28335的多通道緩沖串口McBSP與CodecAD通訊，運用DMA傳輸大量的數(shù)據(jù)而不降低CPU的使用效率；
5）在系統(tǒng)軟件的初始化模塊中對各個部分的初始化順序有嚴(yán)格的要求，否則系統(tǒng)無法正常運行或者不穩(wěn)定。例如在F28335系統(tǒng)配置之后，要先初始化通用輸入輸出口；又例如要先初始化Codec芯片再初始化McBSP和DMA模塊等。
6 系統(tǒng)測試
考慮到在實際流量標(biāo)定中，很難產(chǎn)生變化規(guī)律已知的時變信號，所以，我們采用以下方式對本文提出的時變信號處理方法進(jìn)行測試。首先根據(jù)時變信號模型，利用MATLAB產(chǎn)生已知參數(shù)的時變信號數(shù)據(jù)；其次，將信號數(shù)據(jù)導(dǎo)入Fluke282任意波形信號發(fā)生器，由信號發(fā)生器產(chǎn)生兩路存在相位差且時變的信號；再將信號發(fā)生器的兩路輸出分別接至信號處理系統(tǒng)的兩個輸入通道，運行DSP程序，對信號進(jìn)行處理。相位差在[0.010，20]范圍、變化為1%的測試結(jié)果如表2所示。
由表2可以看出小相位差的相對誤差大于0.1%，結(jié)果并沒有仿真時精度高，其誤差來源為：1）科氏流量計變送器系統(tǒng)采用Codec芯片采集數(shù)據(jù)，其AD位數(shù)為13位，從而限制了計算精度；2）MATLAB產(chǎn)生的信號是64位、雙精度浮點數(shù)，但Fluke282信號發(fā)生器的數(shù)據(jù)是12位，這造成了截斷誤差。
7 結(jié)論
1）面向時變信號，提出由多抽一濾波、自適應(yīng)格型陷波濾波、負(fù)頻率修正的SDTFT遞推算法組合而成的一套數(shù)字信號處理算法，由仿真結(jié)果可以看出該算法能在每個采樣點上輸出傅里葉系數(shù)，實時性好，能跟蹤信號微小的變化，在小流量時仍具有較高的精度，性能優(yōu)越，且算法計算量小，不存在數(shù)值溢出；
2）將整套算法在TMS320F28335DSP搭建的數(shù)字式科氏質(zhì)量流量變送器系統(tǒng)上實時實現(xiàn)，得到了較好的效果?？梢娺@種處理方法能夠?qū)崟r獲得信號間的相位差和時間差，可以提高科氏流量計的動態(tài)響應(yīng)速度，滿足一些測量場合的需要，具有較強(qiáng)的實用性。
參考文獻(xiàn)
[1]徐科軍，倪偉，陳智淵.基于時變信號模型和格型陷波器的科里奧利質(zhì)量流量計信號處理方法[J].儀器儀表學(xué)報，2006，27（6）：596-601.
[2]李祥剛，徐科軍.科氏質(zhì)量流量管非線性幅值控制方法研究[J].電子測量與儀器學(xué)報，2009，6：85-99.
[3]張海濤，涂亞慶.計及負(fù)頻率影響的科里奧利質(zhì)量流量計信號處理方法[J].儀器儀表學(xué)報，2007，3（3）：540-541.ZHANGHT，
[4]TUYQ，ZHANGHT.MethodforCMFsignalprocess-ingbasedontherecursiveDTFTalgorithmwithnegativefrequencycontribution[J].IEEETransactionsOnInstru-mentationAndMeasurement，2008，57（11）：2647-2653.
[5]DERBYHV，BOSET，RAJANS.Methodandapparatusforadaptivelineenhancementincoriolismassflowmetermeasurement[P].USPatentNo.5555190，Sep.10，1996
[6]JACOBSENE，LYONSR.TheslidingDFT[J].IeeeSig-nalProcessingMagazine，vol.20，pp.74-80，Mar2003.
[7]JACOBSENE，LYONSR.AnupdatetotheslidingDFTieeesignalprocessingmagazine，2004.
[8]TexasInstruments.TMS320F28335，TMS320F28334，TMS320F28332，TMS320F28235，DigitalSignalCon-trollers（DSCs）DataManual.February，2008.
[9]徐科軍，于翠欣，蘇建徽，等.基于DSP的科式質(zhì)量流量計信號處理系統(tǒng)[J].儀器儀表學(xué)報，2002，23（2）：
客服1