CN1467485A - 基于aff和sga的科氏质量流量计数字信号处理系统 - Google Patents

Abstract

基于AFF和SGA的科氏质量流量计数字信号处理系统，由信号调理电路、模/数转换器(A/D)、数字信号处理器(DSP)、液晶显示(LCD)电路、键盘输入电路、数/模转换器(D/A)、功率放大电路以及软件组成。它采用多抽一滤波、自适应漏斗型滤波(AFF)和滑动Goertzel算法(SGA)处理采样数据，计算其流量管振动的基频和相位，由频率和两路信号的相位差求出两路信号的时间差，从而得到质量流量值和密度值；由DSP计算得到的频率值产生驱动信号，送到D/A，并经过功率放大，使流量管振动。

Description

基于AFF和SGA的科氏质量流量计数字信号处理系统

技术领域

本发明涉及用于确定科里奥利(简称科氏)质量流量传感器的基频和相位差一种方法和装置，特别是以数字信号处理器(Digital Signal Processor，缩写为DSP)为核心，采用自适应漏斗形滤波器(Adaptive Funnel Shaped Filter，缩写为AFF)，以自适应线性增强信号、跟踪和测量信号频率；采用滑动Goertzel算法(SlidingGoertzel Algorithm，缩写为SGA)，以测量两路信号的相位差。

背景技术

科里奥利质量流量计(以下简称为科氏质量流量计)可以直接测量质量流量，是当前发展最为迅速的流量计之一，具有广阔的应用前景。科氏质量流量计要求其信号处理电路精确地测量来自两个流量传感器信号的相位差，并跟踪其频率的变化。基于模拟和数字电路的科氏质量流量传感器的信号处理方式和系统存在一些局限。为此，人们将数字信号处理方法和DSP应用于科氏质量流量传感器的信号处理，目前有以下三种方案。

(1)基于离散傅里叶变换的方法

美国Micro Motion公司Paul Romano用离散傅里叶变换(DFT)处理科氏质量流量计的输出信号，用TMS系列的DSP作为二次仪表的处理核心(“Coriolicsmass flow rate meter having a substantially increased noise immunity”，US Patent No.4934196，Jun.19，1990)。当非整周期采样时，DFT的计算误差不能满足仪表精度的要求。为此，提出了粗测、细测和频率跟踪的思路。但是，对其中的一些关键技术没有披露。例如，当频率变化时，如何采集过零点，等。合肥工业大学参考其思路，研制了采用DFT的、基于ADSP系列的DSP的信号处理系统，解决了美国专利中没有说明的技术难点，并在细测和频率跟踪方面做了改进。

此类方法存在的问题是：(a)实时性较差。在频率跟踪时，要不断地变化采样频率进行采样和计算，再比较功率谱值的大小，以确定实现整周期采样的频率，其时间长达10秒以上；(b)当输入信号中含有噪声时，频率计算的精度很高，但是，相位差的精度受到影响。其原因在于计算相位差时要用到幅度，而DFT计算出的幅度受噪声影响大。

(2)基于信号幅值的方法

日本富士公司Yoshimura Hiroyuki对来自科里奥利质量流量计中两个传感器的信号进行放大，同时将这两个信号送入差动放大器，得到两个传感器的信号之差。多路转换器将这三个信号顺序送入模/数转换器，再送入DSP做DFT，计算出两个传感器信号的相位差，并从中选一个信号作为参考，去补偿各个传输通道特性差异所造成的误差(“Phase difference measuring apparatus for measuring phasedifference between input signals”，European patent application，EP0791807A2，27.08.1997；“Phase difference measuring apparatus and flowmeterthereof’，European patent application，EP 07022l2A2，20.03.1996)。此方法的问题是：两个传感器信号的相位差很小(一般小于4度)，故差动放大器的幅值很小，极易受到噪声干扰；补偿方法耗时太多，因为计算一次相位差需要采集6路信号。

(3)基于自适应线性增强的方法

美国Micro Motion公司Howard V.Derby等人设计了一个基于DSP的数字信号处理系统，采用自适应线性增强(ALE)技术确定振动管的频率和两个信号之间的相位差，从而更精确地测量质量流量(“Method and apparatus for adaptive lineenhancement in Coriolis mass flow meter measurement”，US Patent No.5555190，Sep.10，1996；“用于科里奥利质量流量计测量的自适应线性增强方法和装置”，中国发明专利申请公开说明书，CN 1190461，1998年8月12日。由于这两项专利在内容上完全一样，下面简称为“Derby的专利”)。Derby的专利中有两个实施例，其信号处理环节都由三部分组成：多抽一滤波、频率估计/线性增强和相位差(时间差)计算。多抽一滤波采用8:1和6:1两级抽取，相位差计算采用改进的Goertzel算法。频率估计/线性增强部分有所不同，第一个实施例的ALE是通过同时对两参数估计的自适应陷波滤波器(Adaptive Notch Filter，缩写为ANF)来实现的，该实施例采用的ANF与最小参数的ANF(Arye Nehoral，“a minimalparameter adaptive notch filter with constrained poles and zeros”，IEEE Transactionson Acoutics，Speech.and Signal Processing，Vol.ASSP-33，No.4，August 1985，pp.987-996)不完全相同，它放宽了对零极点对的限制，没有把零点半径固定为1，而是由算法自动调节，这样可以兼顾频率跟踪的精度和算法收敛速度(即频率跟踪的速度)。同样的，第二个实施例使用了四个自适应陷波滤波器，即在左信道和右信道分别串联两个陷波滤波器。分别设置在左信道和右信道的两个ANF是级联的，其中，第一个滤波器采用一个低Q值(宽陷波带)滤波器，以产生有限的信号增强，但是，能迅速收敛到振动流管基频的变化范围上，然后将从第一个级联ANF输出的信号传输到第二级联的ANF；第二级联的ANF采用一个高Q值(窄陷波带)滤波器以产生比以往的技术或者上述第一实施例中更强的抑制噪声和谐波的作用。

Derby的专利也有局限。在第一实施例中，对两个参数进行估计，增加了算法的复杂性。在第二实施例中，采用两个级联的ANF，算法的复杂程度比第一实施例更高；而且，由于第二级的ANF是以第一级的ANF为基础的，只有在第一级收敛后第二级才会收敛，因此，对于频率变化较大的情况下，收敛速度反而降低。此外，采用的Goertzel算法计算相位差，用定点实现时有可能发生溢出(参见J.A.Beraldin and W.Steenaart，“Overflow analysis of a fixed-point implementationof the Goertzel algorithm，IEEE Trans.Circuits and Systems，Vol.36，No.2，February1989，pp.322-324)。

因此，需要一种既能兼顾频率跟踪精度和算法收敛速度，又不明显增加算法复杂程度的信号处理方法，而且该方法用定点实现时不易溢出。

发明内容

本发明的目的在于解决Derby的专利中存在的缺点，又不失去其优点，即设计了一种信号处理方法和装置，它既能兼顾频率跟踪精度和算法的收敛速度，又不明显增加算法的复杂程度，而且用定点实现时不易溢出。

本发明为了实现发明目的，采用了如下技术方案。信号处理环节由三部分组成：多抽一滤波、频率估计/线性增强和相位差(时间差)计算。

本发明的频率估计/线性增强部分与Derby的专利完全不同，采用的是一种新型的滤波器—自适应漏斗型滤波器(Adaptive Funnel Shaped Filter，缩写为AFF)(Sergio M.Savaresi，“Funnel filter：a New Class of Filters for FrequencyEstimation of Harmonic Signals”，Automatic，Vol.33,No.9，September，1997，pp.1711-1718)。该滤波器与二次自适应陷波滤波器不同，二次自适应陷波经常用于频率估计、线性增强，其特点是只有单一的设计参数

(极点限制因子或称去偏置因子，在Derby的专利中用α表示)，该参数的选择通常需要兼顾跟踪能力和跟踪精度，对于只有单一设计参数的ANF来说，这两者是相互矛盾的：(1)ρ接近1，最大的好处是频率估计的方差和偏差非常小，即跟踪精度高；但是跟踪能力却有所下降，不能跟踪上比较大的频率变化，而且对初始条件敏感；(2)ρ远离1，跟踪能力强，对初始条件不敏感，能够跟踪比较大的频率变化；但稳态时频率估计的方差和偏差大(即跟踪精度不高)。因此，在设计ANF时，只能在跟踪能力和跟踪精度之间取折衷。而AFF有两个设计参数，对零极点进行了更强的限制，阶次也比同级的ANF高，比如对估计单个参数的情况，ANF只需要2次，而AFF需要4次，但算法的复杂程度并没有增加多少，因为所需要估计的参数也还只是一个。AFF代价函数的特点是：宽而平滑的口、陡而窄的底；与ANF相比较，宽的地方比它宽，窄的地方比它窄，因而能够兼顾跟踪能力和跟踪精度两个方面。对于大的突然的频率变化，其跟踪能力比ANF强，而对于小的缓慢的频率变化，跟踪性能与之差不多。它具有ANF跟踪精度高的特点，同时克服ANF对大的、突然的频率变化无法跟踪的缺点。

本发明在相位差计算部分采用更易于定点实现的、且非常适合于时变正弦信号的滑动Goertzel算法(Sliding Goertzel Algorithm，缩写为SGA)(Joe F.Chicharo，Mehdi T.Kilani，“A Sliding Goertzel algorithm”，Signal Processing，Vol.52，No.3，August，1996，pp.283-297)，该算法与传统的(包括改进的Goertzel算法)相比有很多好处：(1)它能够在少于一个信号周期的点数内计算出傅立叶系数，具有较快的获得时间；(2)当用定点算法来实现时，不容易发生数值溢出；(3)非常适合于时变正弦信号。

附图说明

图1是本发明的系统图示。

图2是DSP完成的主要功能。

图3是本发明采用的信号处理方法

图4是ANF的代价函数的频率特性。

图5是ANF与AFF代价函数频率特性比较。

图6是数字信号处理部分的流程。

图7是频率估计/线性增强的自适应算法。

图8是SGA流程图。

图9是SGA的实现框图。

图10是整个跟踪过程图示。

图11是频率突变前的稳态情况。

图12是频率突变后的稳态情况。

图13是自适应算法的两个设计参数的调整。

具体实施方式

下面结合附图和实施例对本发明作详细说明。

科氏流量计是基于科里奥利力的原理而设计的。流体流过测量管时，如果测量管以某一频率振动，则振动的测量管相当于一个匀速转动的参考系。由于流体与测量管具有相对运动，所以会受到科里奥利力的作用。其反作用力作用在测量管的两边上，使测量管发生扭曲。测量管两边装着磁电传感器，输出信号是正比于测量管振动速度的电压信号。当振动管是以一定频率振动时，其角速度按正弦规律变化。故由磁电传感器输出的信号是一正弦信号，其频率为角速度的变化频率，大小正比于角速度。测出其两路信号的相位差和信号频率或者时间差就可以求出流体的质量流量。

图1是本发明的系统图示。两路传感器信号经过信号调理和A/D转换后，进入DSP。信号调理部分包括低通滤波和放大。由于从传感器出来的信号幅值比较小，而且有用的正弦信号落入许多工业噪声范围内，为了充分利用A/D转换器的量程，需要对从传感器出来的信号进行放大；同时为了尽可能地削弱工业噪声的干扰，还采用了低通滤器器。低通滤波和放大环节的设计是本领域的技术人员所熟知的。本发明选用巴特沃斯二阶低通滤波器。

A/D的采样频率固定为38.4KHz。多抽一部分采用Derby的专利所描述的抽选率和滤波器结构。另外还可以采用48KHz的采样频率以及12∶1和6∶1两级多抽一滤波。由于不只一路模拟输入信号，可以采用每路信号一个A/D；或者多路信号经过多路选通后再经过一个A/D。

DSP是信号处理系统的核心，它对两路传感器信号进行数字滤波、频率估计/线性增强和相位差(时间差)计算。对DSP的选择也没有特别限制，几乎可以用任何市面上可以买到的DSP芯片。本发明所叙述的软件由于采用了比Derby的专利更适合定点实现的算法，因此采用TI公司的定点DSP芯片。

DSP还实现数字驱动的功能。

本系统还包括LCD显示、键盘输入电路。

以下的描述是针对典型的科氏流量计应用进行的，在这种流量计中振动流管的基频大约是100Hz。很容易理解，本发明的装置和方法可以应用于处理任何常规的流量计振动基频。

图2是本系统中DSP完成的主要功能图示。DSP主要完成信号采集、信号处理和数字驱动。系统上电后，DSP首先在一段时间内(定时时间)使振动管振动起来，然后进入正常的工作状态。当定时时间没到时，执行振动管起振算法：DSP先施加一带限随机信号，经过D/A转换、功率放大后施加给振动管(可以连续施加一小段时间)，振动系统就能够振动起来(虽然振动不是很理想)。之所以施加随机噪声，是因为随机噪声的频谱很宽，具有很多种频率分量；又由于振动系统具有选频特性，其中，噪声中频率与振动管固有频率相同的分量就会起主要作用，振动管开始振动起来。定时时间到后，振动已经趋于稳定，即振动管以固有频率及稳定的幅值振动，如果振动还没有稳定(这可以通过比较DSP连续的几个信号处理周期计算的信号幅值来判断)，重新设定定时时间，再执行起振算法，直到振动稳定为止。这时DSP进入正常的工作状态：读入传感器信号、对信号进行处理(去噪、求频率和相位)、产生驱动信号(用求出来的频率、相位构造新的正弦驱动信号)。关于定时时间，可以设定为1～5s，可以根据实际情况具体选定。

图3是本发明所采用的信号处理方法。经过转换的数字信号通过路径300(301)传输到抽选单元，抽选单元分为两级。第一级抽选为8∶1，使采样频率从38.4KHz(记为f_s1)，降低到4.8KHz(记为f_s2)，采用的滤波器的传函为： $G\left(z^{-1}\right)=\frac{{(1-z^{-8})}^5}{{(1-z^{-1})}^5}----\left(1\right)$

零极点对消后得到一个36抽头的FIR滤波器，这种滤波器在二次采样频率的各个倍数点具有5个零点，这可以极大的消除混叠在第二级滤波器通带中的那些频率；该滤波器具有最小整数系数，可以用一种精确的计算机运算来表示，从而简化卷积运算的复杂性和提高计算速度。

第二级抽选为6∶1，使采样频率从4.8KHz降低到800Hz(记为f_s)，采用的滤波器是Remez转换算法设计的一个131抽头的FIR滤波器。通带为直流到250Hz，阻带起点为400Hz；通带加权为10^-5，阻带加权为1。

利用两级抽选滤波器具有较高的抗混迭性能。信号经过两级多抽一滤波后通过路径304、305进入自适应漏斗型滤波器(Adaptive Funnel filter，缩写为AFF)环节，该环节完成频率估计/线性增强的功能，它通过自适应算法对输入信号的频率进行估计，同时滤去混在正弦信号中的宽带噪声，输出增强信号；增强后的信号进入滑动Goertzel算法(Sliding Goertzel Algorithm，缩写为SGA)环节，该环节计算增强信号的离散傅立叶系数，从而求出信号的相位(差)；相位差、频率都求出后，就可以求出两路信号的时间差。

下面描述频率计算/线性增强部分。

设两路传感器信号经过两级多抽一后可以表示为：

y_L(n)＝Asin(2πf₀nT+θ₁)+ξ₁(n)                  (2)

y_R(n)＝Asin(2πf₀nT+θ₂)+ξ₂(n)                  (3) $T=\frac{1}{f_s}.----\left(4\right)$ 式中，y_L，y_R分别表示左、右两路信号，A表示信号幅值，f₀是信号的真实频率，θ₁，θ₂分别为两路信号的相位(单位：弧度)，ξ₁，ξ₂分别为两路信号经过两级多抽一滤波后残留在有用信号中的噪声，f_s是经过两级多抽一后的采样频率，T是采样间隔。

本发明采用的自适应算法就是要估计出信号的频率f₀，并且当信号频率发生变化时要能够跟踪信号频率的变化，同时对信号进行信号增强(即有效地滤出噪声分量ξ₁和ξ₂)。

本发明通过AFF来实现频率计算/线性增强功能。AFF的一般形式为 $H\left(z^{-1}\right)=\frac{1}{2}[\frac{C\left(z^{-1}\right)}{D_1\left(z^{-1}\right)}+\frac{C\left(z^{-1}\right)}{D_2\left(z^{-1}\right)}]----\left(5\right)$ 其中

i＝1，2分别为两个陷波滤波器的传递函数，因此可以认为AFF是从ANF衍生出来的。由于陷波滤波器有多种形式，相应的，AFF也有多种形式。本发明采用从理论无偏形式的ANF(称为II型ANF)衍生出来的AFF(称为II型AFF，参见文献Sergio M.Savaresi，Funnel filter：a New Class of Filters for FrequencyEstimation of Harmonic Signals，Automatic，Vol.33，No.9，September，1997，pp.1711-1718给出的结构)。需要指出的是采用其它形式的AFF，比如从直接形式的ANF(称为I型ANF)衍生出来的AFF(称为I型的AFF)，也可以实现这一功能，不同的形式有不同的优缺点。本发明采用的AFF是II型的AFF，其传递函数由式(6)～式(8)决定。

C(z^-1)＝1+2az^-1+z^-2                       (6) $D_1\left(z^{-1}\right)=1+(1+{\mathrm{\ρ}}_1^2)a{z}^{-1}+{\mathrm{\ρ}}_1^2{z}^{-2}----\left(7\right)$ $D_2\left(z^{-1}\right)=1+(1+{\mathrm{\ρ}}_2^2){\mathrm{az}}^{-1}+{\mathrm{\ρ}}_2^2{z}^{-2}----\left(8\right)$ ${\widehat{\mathrm{\Ω}}}_0={\cos }^{-1}(-a)----\left(9\right)$ ${\hat{f}}_0=\frac{\widehat{\mathrm{\Ω}}}{2\mathrm{\π}}{f}_s----\left(10\right)$ 式中，a是待估计的参数，

是待估计的信号频率(数字频率，单位：弧度)，a与

之间有固定的函数关系，即只要估计出参数a，就可以估计出信号的频率。是信号的频率的估计值(单位：Hz)。自适应算法以使式(11)最小为准则调整系数a。 $J\left(\mathrm{\Ω}\right)=\underset{s=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}^{n-s}\mathrm{\ϵ}(s,\mathrm{\Ω}{)}^2----\left(11\right)$ 式(11)称为代价函数(Cost Function)。其中

ε(n，Ω)＝H(z^-1)y(n)                    (12)

是两个极点限制因子，是我们可以利用的设计参数。当ρ₁＝ρ₂时AFF就完全等同于ANF，从这个角度讲，ANF是AFF的特例。由于AFF有两个设计参数，因而算法上就可以采取比ANF更为灵活的方法，既考虑跟踪能力又兼顾跟踪精度。

图4是ANF代价函数的频率特性。陷波滤波器的代价函数的频率特性分为5个区：a1，陷波频率附近的小平坦区；a2，陷波频率左侧的大平坦区；a3，陷波频率右侧的大平坦区；b1，陷波频率左侧的陡斜坡区；b2，陷波频率右侧的陡斜坡区。在a2、a3区会有大量的波纹，因此，在远离陷波频率的区域，寻找代价函数最小值的算法很容易收敛到局部最小值，这是不希望看到的。很自然地，可以将陷波的宽度作为吸引域的度量。文献(Sergio M.Savaresi，Funnelfilter：a New Class of Filters for Frequency Estimation of HarmonicSignals，Automatic，Vol.33，No.9，September，1997，pp.1711-1718)采用Ω⁺-Ω^-(Ω⁺是代价函数的斜率为+1时所对应的频率点，Ω^-是代价函数的斜率为-1时所对应的频率点)来描述吸引域，根据这样的定义，有下面的结论(考虑相同的估计方差的情况下，ρ²＝0.98， ${\mathrm{\ρ}}_1^2=0.84,$ ${\mathrm{\ρ}}_2^2=0.988$ )：

可见AFF的吸引域比ANF大很多，因此其跟踪能力也就比ANF强很多，能跟踪上比较大的频率变化。

图5是AFF与ANF代价函数的频率特性的比较，从它们的代价函数的比较得知，AFF的跟踪能力比ANF的强(口比它宽)，跟踪精度跟它差不多(底差不多一样窄)；从图中还可以看出ANF的自适应算法有可能收敛到局部最小值，而AFF可以避免。需要指出的是，本发明仅仅采用了文献(Sergio M.Savaresi，Funnelfilter：a New Class of Filters for Frequency Estimation of Harmonic Signals，Automatic，Vol.33，No.9，September，1997，pp.1711-1718)提出的AFF的结构，而自行推导了适合于RML算法的有关公式，同时，经过大量的仿真确定了两个设计参数ρ₁ ²，ρ₂ ²的选择。

图6是数字信号处理部分的流程。

借鉴Derby的专利，采用SNR故障检测。自适应算法中的输入信号是经过两级多抽一后的左、右路信号，用符号y_L和y_R表示。

左信道和右信道的输出信号都可以被用作加权自适应单元的反馈信号。虽然在加权自适应单元中同时使用两信道的输出信号作反馈信号是可能的，但是这样所产生的益处与增加的计算复杂程度相比并不占优，因此只用其中一个信道(比如右信道)作反馈，而计算出的加权自适应参数被同时用在左右信道的ANF，从而使两个传感器输出信道经过同样的处理。在本发明所采用的算法中，用右信道的信号作为加权自适应的反馈信号，由该反馈信号确定好滤波器的系数后，两路信号又通过相同的滤波器(AFF)，得到增强后的信号，用x_L和x_R表示。在下面的公式中，没有严格区分y_L和y_R，而是用y表示；也没有严格区分x_L和x_R，而是用x表示。

图7是频率估计/线性增强部分的自适应算法(AFF环节)的流程图，本发明采用的AFF的自适应算法的主要公式为：

(1)有关向量、矩阵及符号的意义：

输入信号：y(包括y_L，y_R)

先验预测误差：ε(n)

后验预测误差： ε(n)

设计参数(极点限制因子和去偏置因子)：ρ₁ ²，ρ₂ ²

遗忘因子：λ

中间状态变量：F(n)

自回归向量：Φ(n)＝[φ₁(n)，φ₂(n)]^T；

负梯度向量：Ψ(n)＝[₁(n)，₂(n)]^T；

待估计参数a(t)构成的向量：θ(n)＝[a(n)，a²(n)]^T；

协方差矩阵： $P\left(n\right)=\left[\begin{array}{c}{p}_1\left(n\right),p_2\left(n\right)\\ p_3\left(n\right),p_4\left(n\right)\end{array}\right]$

中间状态矩阵：M(n)

增强后的信号：x(包括x_L，x_R)

数字角频率：Ω、

Ω_k等

正弦信号频率估计值：

(2)初始化：

θ(0)＝0，P(0)＝10^-2I，P_SNR＝10^-4I，Ψ(0)＝Φ(0)＝0；y(-i)＝0，i＝1，2 ${\mathrm{\ρ}}_{10}^2={0.7}^2,$ ${\mathrm{\ρ}}_{1\∞}^2={0.995}^2,$ ρ_1decay＝0.998，ρ_1SNR＝0.80 ${\mathrm{\ρ}}_{20}^2={0.8}^2,$ ${\mathrm{\ρ}}_{2\∞}^2={0.995}^2,$ ρ_2decay＝0.99，ρ_2SNR＝0.95

λ₀＝0.95，λ_∞＝1，λ_decay＝0.995，λ_SNR＝0.97

(3)主循环： $F\left(n\right)=y\left(n\right)+\frac{2+{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2}{2}y(n-2)+\frac{{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2}{2}y(n-4)-({\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2)\mathrm{\ϵ}(n-2)-{\mathrm{\ρ}}_1^2{\mathrm{\ρ}}_2^2\mathrm{\ϵ}(n-4)----\left(14\right)$ ε(n)＝F(n)-Ψ^T(n)θ(n-1)                                    (15) $P\left(n\right)=\frac{1}{\mathrm{\λ}\left(n\right)}[P(n-1)-\frac{P(n-1)\mathrm{\Ψ}\left(n\right){\mathrm{\Ψ}}^T\left(n\right)P(n-1)}{\mathrm{\λ}\left(n\right)+{\mathrm{\Ψ}}^T\left(n\right)P(n-1)\mathrm{\Ψ}\left(n\right)}]----\left(16\right)$ θ(n)＝θ(n-1)+P(n)Ψ(n)ε(n)                                (16)ε(n)＝F(n)-Ψ^T(n)θ(n)                                     (17)x(n)＝y(n)- ε(n) ${\mathrm{\φ}}_1\left(n\right)=-\frac{6+{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2}{2}y(n-1)-\frac{2+3{\mathrm{\ρ}}_1^2+3{\mathrm{\ρ}}_2^2}{2}y(n-3)+(2+{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2)\mathrm{\ϵ}(n-1)----\left(18\right)$ $+({\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2+2{\mathrm{\ρ}}_1^2{\mathrm{\ρ}}_2^2)\mathrm{\ϵ}(n-3)$ ${\mathrm{\φ}}_2\left(n\right)=-(2+{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2)y(n-2)+(1+{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2+{\mathrm{\ρ}}_1^2{\mathrm{\ρ}}_2^2)\mathrm{\ϵ}(n-2)----\left(19\right)$ Φ(n)＝[φ₁(n)，φ₂(n)]^T                                  (20) $M\left(n\right)=[(2+{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2)\mathrm{\Ψ}(n-1)+({\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2+2{\mathrm{\ρ}}_1^2{\mathrm{\ρ}}_2^2)\mathrm{\Ψ}(n-3),(1+{\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2+{\mathrm{\ρ}}_1^2{\mathrm{\ρ}}_2^2)\mathrm{\ψ}(n-2)]----\left(21\right)$ $\mathrm{\ψ}\left(n\right)=\mathrm{\Φ}\left(n\right)-M\left(n\right)\mathrm{\θ}\left(n\right)-({\mathrm{\ρ}}_1^2+{\mathrm{\ρ}}_2^2)\mathrm{\Ψ}(n-2)-{\mathrm{\ρ}}_1^2{\mathrm{\ρ}}_2^2\mathrm{\Ψ}(n-4)----\left(22\right)$ ${\mathrm{\ρ}}_1^2\left(n\right)={\mathrm{\ρ}}_1^2(n-1){\mathrm{\ρ}}_{1\mathrm{decay}}+(1-{\mathrm{\ρ}}_{1\mathrm{decay}}){\mathrm{\ρ}}_1^2(\∞)----\left(23\right)$ ${\mathrm{\ρ}}_2^2\left(n\right)={\mathrm{\ρ}}_2^2(n-1){\mathrm{\ρ}}_{2\mathrm{decay}}+(1-{\mathrm{\ρ}}_{2\mathrm{decay}}){\mathrm{\ρ}}_2^2(\∞)----\left(24\right)$ λ(n)＝λ(n-1)λ_decay+(1-λ_decay)λ(∞)                   (25)

频率的估计采用的公式为 $\widehat{\mathrm{\Ω}}\left(n\right)={\cos }^{-1}(-a\left(n\right))----\left(26\right)$ $\hat{f}\left(n\right)=\frac{\widehat{\mathrm{\Ω}}\left(n\right)}{2\mathrm{\π}}{f}_s----\left(27\right)$

加权自适应单元计算所得值(即a(t)的值)用于计算相位和在SGA环节中计算相位差、时间差Δt。频率计算单元计算频率并产生Goertzel滤波器的加权信息。相位计算单元从频率计算单元接受Goertzel加权信息和频率信息。相位计算单元利用具有两个Hanning窗口的傅立叶分析技术确定滤波信号的相位选择8个流管周期作为优选的窗口长度。假设一个给定的频率期望值为f_q，则优选的窗口长度(L＝2N)由下式确定： $L=2N=8\mathrm{floor}\left(\frac{f_s}{f_q}\right)----\left(28\right)$ 式中，floor(x)表示不大于x的整数。

f_s为两级多抽一后的频率，比如f_s＝800，f_q＝100，则2N＝64。Hanning窗表示为： $h\left(n\right)=\frac{1}{2}(1-\cos \frac{2\mathrm{\πn}}{2N-1}),n=\mathrm{0,1},\·\·\·2N-1----\left(29\right)$

为了尽可能最大限度的利用时域数据，本发明借鉴Derby的专利采用的重叠Hanning窗口的平行计算。

增强后的信号经过Hanning窗后，采用Goertzel算法计算其离散傅立叶系数。采用重叠Hanning窗口的平行计算后每个半窗长(每N个经过两级多抽一后的时域点)进行一次相位、频率和Δt的计算。此时频率的计算采用N个估计值的平均值求得： ${\widehat{\mathrm{\Ω}}}_0\left(n\right)={\cos }^{-1}(-\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}a\left(i\right))----\left(30\right)$ ${\hat{f}}_0\left(n\right)=\frac{{\widehat{\mathrm{\Ω}}}_0\left(n\right)}{2\mathrm{\π}}{f}_s----\left(31\right)$

下面描述相位差计算部分。

本发明采用滑动Goertzel算法(Joe F.Chicharo，Mehdi T.Kilani，A Sliding Goertzelalgorithm，signal processing，Voi.52，No.3，August，1996，pp.283-297)，用该算法分别计算两路增强后的信号的离散傅立叶系数，根据离散傅立叶系数求出两路信号的相位，再求出两路信号之间的相位差，最后计算出两路信号之间的时间差。由于AFF环节已经估计出了信号的频率并同时对信号进行了线性增强，采用Goertzel类的算法可以直接计算给定频率点的离散傅立叶系数，从而大大减少了计算量(与传统的计算离散傅立叶系数的DFT算法、FFT算法相比较)。

图8是SGA的流程图，与此相对应，下面给出本发明采用的SGA的有关步骤和公式：

(1)确定需要计算的频率点

(2)状态条件初始化为零，即ν_k(0)＝0；ν_k(-i)＝0，k＝1，2

(3)计算共振滤波器的的输出 $v_k\left(n\right)=x\left(n\right)+2\cos {\widehat{\mathrm{\Ω}}}_0\·v_k(n-1)-v_k(n-2)----\left(32\right)$ ${\hat{A}}_k=\frac{2{(v}_k\left(N\right)-v_k(N-1)\cos {\widehat{\mathrm{\Ω}}}_0}{N},k=\mathrm{1,2}----\left(33\right)$ ${\hat{B}}_k=-\frac{2v_k(N-1)\sin {\widehat{\mathrm{\Ω}}}_0}{N},k=\mathrm{1,2}----\left(34\right)$ $\left[\begin{array}{c}{a}_k\\ b_k\end{array}\right]=\left[\begin{array}{c}\cos {\widehat{\mathrm{\Ω}}}_{0}n,-\sin {\widehat{\mathrm{\Ω}}}_{0}n\\ \sin {\widehat{\mathrm{\Ω}}}_{0}n,\cos {\widehat{\mathrm{\Ω}}}_{0}n\end{array}\right]\left[\begin{array}{c}{\hat{A}}_k\left(n\right)\\ {\hat{B}}_k\left(n\right)\end{array}\right],k=\mathrm{1,2}----\left(35\right)$

(4)计算相位 ${\widehat{\mathrm{\θ}}}_k={\tan }^{-1}\frac{a_k}{b_k}----\left(36\right)$

(5)两路信号的相位差 $\widehat{\mathrm{\φ}}={\widehat{\mathrm{\θ}}}_1-{\widehat{\mathrm{\θ}}}_2----\left(37\right)$

(6)计算时间差 $\mathrm{\Δt}=\frac{\widehat{\mathrm{\φ}}}{{\widehat{\mathrm{\Ω}}}_0}\·\frac{1}{f_s}----\left(38\right)$

图9是实现SGA算法的结构图。

需要指出的是，由于自适应算法估计出参数a，因而在上面的式(32)～式(34)中实际上并不需要计算余弦值 和正弦值

而可以分别采用式(39)、式(40)代替： $\cos {\widehat{\mathrm{\Ω}}}_0=W_1----\left(39\right)$ $\sin {\widehat{\mathrm{\Ω}}}_0=W_2----\left(40\right)$ 其中， $W_1=-\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}a\left(i\right)----\left(41\right)$ $W_2=\sqrt{1-{W_1}^2}----\left(42\right)$

图10～图13是部分结果。下面作简要说明。

输入信号的幅值为10mV，频率为100Hz，相位差为4°，所加噪声为正态分布噪声，服从N(0，0.4)分布，噪声幅值为0.89mV(按式(43)的定义信噪比为34)。所用的数据在两级多抽一后为4000个，即进入DSP是有48×4000个。数据窗长度为64(即L＝64，N＝32)，采用重叠窗的并行处理，每个半窗长(称为一个信号处理周期)计算一次频率、相位差、时间差。考虑这样的情况：信号的频率在第2000点时突然发生了变化，从100Hz～110Hz，这样的变化幅度对科氏质量流量计来说是比较大了，因而足以说明问题。 $\mathrm{SNR}=101g\frac{A_s^2}{{\mathrm{\σ}}^2}----\left(43\right)$ 式中，A_s为正弦量的幅值，σ²为正态分布噪声的方差。

图10是整个跟踪过程。图10a)表示频率估计的相对误差；图10b)表示相位差估计的相对误差；图10c)表示频率估计值。图10～13中，横坐标k表示第k个信号处理周期，而n表示第n个数据点(

)。

图11表示频率突变前的稳态情况，即自适应算法从初始状态(任意)到收敛到突变前的频率100Hz时的稳态情况。

图12表示频率突变后的稳态情况，从图中可以看出，当信号频率从100Hz变化到110Hz时，自适应算法能够在几个周期内跟踪上信号频率的变化。下面是一些数据：

突变前，稳态时的频率估计误差≤5.8679×10^-5

突变前，稳态时的相位差估计误差≤6.6973×10^-4

突变后，稳态时的频率估计误差≤1.4124×10^-5

突变后，稳态时的相位差估计误差≤8.1189×10^-4

图13是本发明采用的算法中，设计参数的调整曲线。由于采用了SNR故障检测，当信号频率发生变化时，设计参数先保持在某一数值上，待算法接近收敛时，再按指数函数变化。

Claims (6)

1.一种基于AFF和SGA的科氏质量流量计数字信号处理系统，由信号调理电路、模/数转换器(A/D)、数字信号处理器(DSP)、液晶显示(LCD)电路、键盘输入电路、数/模转换器(D/A)、功率放大电路以及软件组成；所述的流量计中的两个磁电式传感器将所感受到流量信号转换成电信号；所述的磁电式传感器输出的电信号经过所述的信号调理电路和A/D送到DSP；所述的DSP处理采样数据，计算流量管振动基频处的频率和相位，由频率和两路信号的相位差求出两路信号的时间差，从而得到质量流量值和密度数值，供LCD显示；由DSP计算得到的频率值产生驱动信号，送到D/A，并经过功率放大，送到电磁激振器的线圈，使流量管振动；其特征在于所述的DSP采用自适应漏斗型滤波器(AFF)对经过两级多抽一滤波后的采样信号进行频率估计和跟踪，能兼顾跟踪能力和跟踪精度两个方面，同时滤去混在正弦信号中的宽带噪声，输出增强信号；采用更易于定点实现的、且非常适合于时变正弦信号的滑动Goertzel算法(SGA)计算增强信号的离散傅立叶系数，从而求出信号的相位差。

2.根据权利要求1所述的一种基于AFF和SGA的科氏质量流量计数字信号处理系统，其特征在于采用DSP实时完成数字信号处理任务；DSP作为系统核心控制系统其他部分协调工作。

3.根据权利要求1所述的一种基于AFF和SGA的科氏质量流量计数字信号处理系统，其特征在于采用两级多抽一滤波器对两路传感器信号进行滤波，以消除随机噪声的影响。

4.根据权利要求1所述的一种基于AFF和SGA的科氏质量流量计数字信号处理系统，其特征在于采用自适应漏斗型滤波器对经过多抽一滤波器后的信号进行自适应陷波滤波，计算流量管振动的基频，并跟踪基频的变化；将多抽一滤波后的信号减去经过自适应漏斗型滤波器的信号，使流量信号得到增强。

5.根据权利要求1所述的一种基于AFF和SGA的科氏质量流量计数字信号处理系统，其特征在于采用滑动Goertzel算法计算基频处的离散傅立叶系数，从而求出信号的相位差；当相位差和频率求出后，求出两路信号的时间差。

6.根据权利要求1所述的一种基于AFF和SGA的科氏质量流量计数字信号处理系统，其特征在于DSP先施加一带限随机信号，经过D/A转换、功率放大后施加给流量管，使流量管振动起来；在正常工作状态下，DSP用求出来的基频和相位构造新的正弦驱动信号，去驱动流量管。

Images