摘 要：采用比Ⅱ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）

    k₀在经过一段时间自适应后应该收敛到-cosω，ω是信号y（n）的归一化频率。α可决定着陷波的带宽。

    自适应算法如下^[10]：

    

    

    其中是k₀（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值较大时Δ₁和Δ₂可以忽略。经计算得到

    

    则信号的相位为φ（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算法实现了实时谱分析和两信号之间相位差的计算。仿真结果表明，本文的方法计算精度较高，取得了令人满意的效果。