摘 要：本文采用自適應ⅡR陷波器對科里奧利質量流量計信號進行濾波，計算其頻率；采用自適應譜線增強器從含有噪聲的數據中提取出基頻信號。然后采用滑動Goertzel算法實時計算兩路信號之間的相位差，并根據頻率和相位差計算出時間差，測得質量流量。仿真結果表明本文所研究的方法是有效的。

關鍵字：HKC科里奧利質量流量計 自適應ⅡR陷波器 自適應譜線增強 滑動Goertzel算法

    1 引言

    HK-CMF科里奧利質量流量計（以下簡稱為科氏流量計），是一種可以直接實現質量流量測量并可同時獲取流體密度值的流量計，是當前發展迅速的流量計之一，具有廣闊的應用前景?？剖狭髁坑嬕笃湫盘柼幚黼娐纺軌蚓_地測量出來自兩個流量傳感器信號的時間差。目前，國內外絕大部分采用基于模擬和數字電路的信號處理方法，如：放大、濾波、整形和計數等，這種處理方法存在很多局限^[1]，如對噪聲比較敏感、模擬濾波器會改變信號的幅值和相位等。隨著數字信號處理技術和DSP芯片的發展，國內外的相關研究機構開始將數字信號處理技術應用于科氏流量計的信號處理中^[1~6]。

    在文獻[1]中我們提出了一種信號處理方法，從含有各種噪聲干擾的信號中，準確測得質量流量。首先采用兩級多抽一濾波器，對科氏流量計中的兩個磁電式傳感器的輸出信號進行濾波，以減少隨機噪聲的影響；再用自適應陷波濾波方法，抑制確定性噪聲的干擾，提取流量管振動基頻的信號并測得其頻率。然后采用加漢寧窗的DFT頻譜分析方法，得到振動管基頻處的相位差和時間差，從而測得質量流量。但是這種方法存在一些問題，首先兩級多抽一濾波器的計算量比較大，其次所采用的譜分析方法不能做到實時輸出相位差，只能實現每64個點（即一個窗的寬度）計算一次時間差。本文在原有方法的基礎上進行了改進，去除了兩級多抽一濾波器，并采用滑動Goertzel算法^[10]。進行譜分析，從而可以在每個采樣點都測得兩路信號之間的相位差及時間差。

    2 自適應ⅡR陷波器和譜線增強

    本文中我們采用的仍然是零極點約束的ⅡR陷波濾波器^[7]。零極點約束是指陷波器的零極點滿足如下條件：對于每對零極點，零點在單位圓上，并且位于陷波頻率處。極點則在單位圓內，且與零點同一角度，并盡可能地靠近零點，極點與原點的距離為ρ。對于科氏流量計信號可得陷波器的傳遞函數為

         （1）

    假設信號為y（n）=Asin（ωn+φ）+e（n），當α=-2cosω時，

    （1+αz^-1+z^-2）e（n）=（1+αz^-1+z^-2）y（n）    （2）

成立，其中z^-1為單位延遲因子，即y（n）z^-1=y（n-1）。式（2）表明當陷波頻率等于信號頻率且ρ=1時，陷波器的輸出中只含有白噪聲e（n），正弦信號被濾掉，陷波器相當于y（n）的一個白化濾波器。對于0<ρ<1，陷波器的輸出為

        （3）

    當ρ→1時，。通常信號的頻率是未知的，因此需要對α進行估計。陷波器的輸出誤差為，根據遞推預測誤差理論^[8]，取代價函數

        （4）

    α的估計可表示為

        （5）

    由于ρ趨近于1，，因此可由下式遞推計算：

        （6）

    其中

    P（n）可由下式遞推計算：

    λ（n）為遺忘因子：

    λ（n）=λ₀λ（n-1）+（1-λ₀）λ_∞    （10）

    在以上討論中，我們一直假設ρ為定值。但在實際算法中，由于ρ決定每個陷井的帶寬，在輸入信號的先驗知識未知的情況下，如果ρ非常趨近于1，即極點靠近零點，則陷井可能因為過窄而無法落在信號頻率處，也就無法感知信號的存在。因此在陷波的開始階段ρ應取得稍小一點，使算法收斂到期望的傳遞函數上，然后再增大ρ。為此ρ可改寫為ρ（n），

    ρ（n）=ρ₀ρ（n-1）+（1-ρ₀）ρ_∞    （11）

    最后將原先的信號減去陷波器的輸出信號就是我們所希望得到的振動管的信號，同時通過可以計算出信號的頻率。在將計算出的頻率用于相位差的計算之前我們先采用下面的方法對其進行平滑處理，以減少陷波器計算相位差時所引入隨機誤差。

    3 相位差及時間差的計算

    常規Goertzel算法的傳遞函數如下^[9]：

    其頻率為等間隔分布，當信號的頻率恰好對應于某個k值時，可以獲得精確的結果。而當信號頻率落在頻率間隔內時，就會由于泄漏問題而產生較大誤差。泄漏問題對于計算相位的影響要遠遠大于對功率譜的計算^[5]。

    滑動Goertzel算法^[10]克服了常規Goertzel算法的這一缺點，它用實際的信號頻率替換傳遞函數中的2πk/N，因此在信號頻率不變的情況下，經過一定時間的收斂，可以精確地計算出信號的傅立葉系數，并且可以在每個采樣點計算一次傅立葉系數?；瑒覩oertzel算法如圖1所示：

    設固定頻率信號y（n）的傅立葉展開為yn=αsin（nω）+bcos（nω），其z變換為

    V（z^-1）是中間變量v（n）的z變換，對其進行z反變換^[10]，得

    從式（16）-（19）可以看出在n的值比較大時Δ₁和Δ₂可以忽略。經過計算我們得到

    則信號的相位為：

    其中奪和分別為兩路信號的相位估計。

    由于滑動Goertzel算法的相位差計算結果是振蕩收斂，所以對于每一點的計算結果我們也采用了與頻率計算相同的平滑處理方法。

    時間差的計算公式為

    其中為信號頻率的估計，f_s為采樣頻率。

    4 仿真結果

    4.1 自適應陷波器的仿真

    由于陷波器的仿真結果在文獻[1]中已給出，為節約篇幅，我們在這里不再重復，但給出陷波器的一些變量的初始值。

    4.2 整體系統的仿真

    為與原方法對比，我們仍使用文獻[1]中的信號，即信號幅值為10V，諧波頻率為80Hz，幅值為3V，初始相位為零，隨機噪聲的幅值為0.9984V，服從正態分布且去均值，采樣點數為20000。仿真結果如表1所示，表中的數據為Δ_t的計算相對誤差，單位為%。

表1 不同頻率和相位差下的兩種方法計算誤差的比較

    需要說明的是，由于本文所采用的時間差計算方法可以做到實時計算，所以不能與原方法做到真正意義上的一一對比。因此表1所列出的新方法相對誤差是第6001－20000點誤差的均值。圖2給出了頻率為103Hz，相位差為0.1時相位差和時間差的相對計算誤差變化趨勢。

    由表1和圖2可以看出雖然在某些時刻新方法的計算誤差要大于原方法，但整體上改進后的方法比原方法精度更高。重要的是改進后的方法充分利用了采集來的數據，真正實現了在線計算，這是該方法的一大優勢。

    5 結語

    本文對原有的科氏流量計信號處理方法進行了改進，用滑動Goertzel算法替代了原先較為經典的DFF方法以實現實時譜分析和兩信號之間相位差的計算。從仿真結果中可以看出改進后的方法在整體上計算精度要高于原方法。

    參考文獻：

    [1] 徐科軍，倪偉.一種HKC科里奧利質量流量計的信號處理方法，計量學報，2001，22（4）：254-257.
    [2] H.V.Derby，T.Bose，S.Rajan.Method and Apparatus for Adaptive Line Enhancement in Cofiolis Mass Flow Meter Measurement，US Patent No.5555190，Sep.10，1996.
    [3] Yoshimura.Hiroyuki.Phase Diference Measuring Apparatus for Measuring Phase Diference between Input Signals，Euro-pean PatentApplication.EPO791807A2.1997.
    [4] 徐科軍，姜漢科，蘇建徽等.科氏流量計信號處理中頻率跟蹤方法的研究，計量學報，1999，19（4）：304-307.
    [5] 倪偉，用于HK-CMF科里奧利質量流量計的相位差測量方法與裝置，學士學位論文，合肥工業大學自動化研究所，1999.
    [6] 于翠欣，HK-CMF科里奧利質量流量計數字信號處理系統的研制，碩士學位論文，合肥工業大學自動化研究所，2000.
    [7] Nehorai，A Minimal Parameter Adaptive Notch Filter with Comtrained Poles and Zeros，IEEE Traus .Acoust，Speech Signal Processing，1985，ASSP-33：983-996.
    [8] L.Ljung， Identification：Theory for the User（SecondEdition），Prentice Hall FIR，1999，清華大學出版社，2OO2.
    [9] S.K.Mitra，Digital Signal Processing：A Computer-Based Approach（Second Edition），New York，NY：McGraw-Hill，2001，清華大學出版社，2001.
    [10] J.Chicharo and M.Kilani，A Sliding GcertzelAlgorithm，Signal Preessing，1996，52：283-297.