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