摘 要：采用比ⅡR型計算更為簡單的自適應格型陷波濾波器對科里奧利質量流量計的流量傳感器的輸出信號進行濾波并求得其頻率；并應用自適應譜線增強器從含有噪聲的信號中提取出流量管振動的基頻信號；然后采用滑動Goertzel算法計算兩路信號之間的實時相位差，再根據頻率和相位差計算出時間差最終求得質量流量。
關鍵字：計量學 科里奧利質量流量計 自適應格型陷波器 自適應譜線增強 滑動Goertzel算法
1 引言
科里奧利質量流量計（以下簡稱為質量流量計）可以直接測量質量流量，并可同時獲取流體密度值，應用廣泛，有很好的發展前景。但是，目前的國內外產品采用基于模擬和數字電路的信號處理方法存在很多局限。例如，對噪聲比較敏感、模擬濾波器會改變信號的幅值和相位、測量的是合成波的相位差等，使其在現場的測量精度低于實驗室標定的精度[1]。為了解決上述問題，也為提升流量測量儀表的性能，國內外研究機構和公司將數字信號處理方法和數字信號處理器（DSP）應用于HKC質量流量計的信號處理。文獻[2]采用基于DFT（離散傅里葉變換）的方法計算科氏質量流量傳感器的頻率和相位差，提出了粗測、細測和頻率跟蹤的思路。文獻[3，4]采用DFT和基于DSP的信號處理系統，并在細測和頻率跟蹤方面做了改進。文獻[5，6]研究了基于信號幅值的計算相位差的方法來處理科氏質量流量傳感器的信號，并進行補償。這些方法均存在一些不足之處。文獻[7，8]提出了基于陷波濾波器的信號處理方法，首先采用兩級多抽一濾波器，對HKC質量流量計中的兩個磁電式傳感器的輸出信號進行濾波，以減少隨機噪聲的影響；再采用自適應陷波濾波方法，抑制確定性噪聲的干擾，提取流量管振動基頻的信號并測得其頻率。然后，采用加漢寧窗的DFT或常規Goertzel算法進行譜分析，得到振動管基頻處的相位差和時間差，從而測得質量流量。但是，這種方法存在一些問題，如兩級多抽一濾波器的計算量比較大，所采用的直接型ⅡR自適應陷波器計算復雜且對初始值的選取較為敏感，譜分析方法不能做到實時輸出相位差，只能每64個點（即一個漢寧窗的寬度）計算一次時間差。
本文研究了一種新的信號處理方法，該方法中沒有使用多抽一濾波器，從而節約了計算時間，并采用一種計算更為簡單的自適應格型陷波器[9，10]來對信號進行陷波濾波。最后采用滑動Goertzel算法[12]進行實時譜分析，從而可以在每個采樣點都能測得兩路信號之間的相位差及時間差，實現在線計算。
2 自適應格型陷波器
圖1所示的為一種級聯的全極點和全零點格型陷波器，其傳遞函數為
（1）
k0在經過一段時間自適應后應該收斂到-cosω，ω是信號y（n）的歸一化頻率。α可決定著陷波的帶寬。
自適應算法如下[10]：
其中是k0（n）的估計值，是經過平滑處理后的值，γ=0.5。y（n）-x（n）即為去除噪聲后的增強信號。
在將計算出的頻率用于相位差的計算之前我們先采用下面的方法對其進行平滑處理，以減少陷波器計算相位差時所引入的隨機誤差。
3 相位差及時間差的計算
常規Goertzel算法的傳遞函數如下：
其頻率為等間隔分布，當信號的頻率恰好對應于某個k值時，可以獲得精確的結果。而當信號頻率落在頻率間隔內時，就會由于泄漏問題而產生較大誤差。泄漏問題對于計算相位的影響要遠遠大于對計算功率譜的影響[7]。
滑動Goertzel算法[12]克服了常規Goertzel算法的這一缺點，它用實際的信號頻率替換傳遞函數中的2πk/N，因此在信號頻率不變的情況下，經過一定時間的收斂，可以精確地計算出信號的傅里葉系數，并且可以在每個采樣點計算一次傅里葉系數。滑動Goertzel算法如圖2所示。
設固定頻率信號Y（n）的傅里葉展開為y（n）=asinnω+bcosnω，其z變換為
V（z-1）是中間變量v（n）的z變換，對其進行z反變換[12]，得
從式（15）~（18）可以看出在n值較大時Δ1和Δ2可以忽略。經計算得到
則信號的相位為φ（n）=arctan[b（n）/a（n）] （21）
相位差的計算公式為
其中和為兩路信號的相位估計值。
由于滑動Goertzel算法的相位差計算結果是振蕩收斂，所以，對于每一點的計算結果也采用了與頻率計算相同的平滑處理方法。時間差的計算公式為
其中為信號頻率的估計值，fs為采樣頻率。
4 仿真結果
仿真實踐中，使用的信號幅值為10V，隨機噪聲的幅值為1.0547V，服從正態分布且去均值，采樣點數為20000。仿真結果如表1所示，表中的數據為第6001~20000點Δt的計算相對誤差的均值，單位為% 。
表1 不同頻率和相位差下的計算誤差
相位差 | 頻率/Hz | ||
97.0 | 100.0 | 103.0 | |
0.0001 | 0.017357 | 0.034520 | 0.031693 |
0.0010 | 0.017356 | 0.034521 | 0.013702 |
0.0030 | 0.017356 | 0.034521 | 0.013701 |
0.0050 | 0.017356 | 0.034521 | 0.013701 |
0.0100 | 0.017357 | 0.034520 | 0.013701 |
0.0500 | 0.017362 | 0.034516 | 0.013701 |
0.1000 | 0.017368 | 0.034511 | 0.013702 |
0.5000 | 0.017418 | 0.034472 | 0.013706 |
1.0000 | 0.017479 | 0.034423 | 0.013711 |
1.5000 | 0.017540 | 0.034373 | 0.013717 |
2.0000 | 0.017601 | 0.034323 | 0.013724 |
4.0000 | 0.017845 | 0.034122 | 0.013760 |
圖3給出了頻率為103Hz、相位差為0.1時相位差和時間差的相對計算誤差變化趨勢。
值得注意的是將本文所研究的方法應用于實際的系統時要注意到DSP有限精度的影響。從式（15）和式（16）中可以看出隨著采樣點n的增大而增大。當v增大到一定程度時會超出DSP所能表示的數的范圍，造成溢出，這一點在16位定點DSP上尤其明顯。對于這種情況的解決辦法是一旦檢測到v溢出，就將滑動Goertzel算法的所有變量清零，重新開始計算。這種方法相當于用一個矩形窗對數據進行截取，而這個窗的寬度是不固定的，取決于DSP的位數，所以，本文方法比較適合于以32位浮點DSP為處理核心的系統。
5 結束語
格型陷波器比原先的直接型ⅡR陷波器計算更為簡單并且初始值易于選擇。滑動Goertzel算法實現了實時譜分析和兩信號之間相位差的計算。仿真結果表明，本文的方法計算精度較高，取得了令人滿意的效果。
