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