CN1194210C

CN1194210C - Ke's mass flowmeter digital signal processing system based on AFF and SGA

Info

Publication number: CN1194210C
Application number: CNB031089445A
Authority: CN
Inventors: 徐科军; 徐文福
Original assignee: Hefei University of Technology
Current assignee: Hefei University of Technology
Priority date: 2003-04-04
Filing date: 2003-04-04
Publication date: 2005-03-23
Anticipated expiration: 2023-04-04
Also published as: CN1467485A

Abstract

The present invention relates to a digital signal processing system of a Coriolis mass flow meter based on AFF and SGA. the digital signal processing system is composed of a signal conditioning circuit, an analog-digital converter (A/D), a digital signal processor (DSP), a liquid crystal display (LCD) circuit, a keyboard input circuit, a digital-to-analog converter (D/A), a power amplification circuit and software. The digital signal processing system uses multi-to-one filtration, self-adaptive funnel type filtration (AFF) and sliding Goertzel algorithm (SGA) to process sampled data and calculate the vibration fundamental frequency and phase of a flow tube. The time difference of two ways of signals is figured out by the phase difference of the frequency and the two ways of singles; consequently, a mass flow value and a density value are obtained. A frequency value obtained by the calculation of DSP generates a driving signal which is sent to D/A to be amplified in power for making the flow tube vibrate.

Description

Coriolis mass flowmeter digital signal processing system based on AFF and SGA

技术领域technical field

本发明涉及用于确定科里奥利(简称科氏)质量流量传感器的基频和相位差一种方法和装置，特别是以数字信号处理器(Digital Signal Processor，缩写为DSP)为核心，采用自适应漏斗形滤波器(Adaptive Funnel Shaped Filter，缩写为AFF)，以自适应线性增强信号、跟踪和测量信号频率；采用滑动Goertzel算法(SlidingGoertzel Algorithm，缩写为SGA)，以测量两路信号的相位差。The present invention relates to a method and device for determining the fundamental frequency and phase difference of a Coriolis (referred to as Coriolis) mass flow sensor, especially with a digital signal processor (Digital Signal Processor, abbreviated as DSP) as the core, using Adaptive Funnel Shaped Filter (AFF for short) to adaptively enhance the signal linearly, track and measure the frequency of the signal; use the sliding Goertzel algorithm (SlidingGoertzel Algorithm, short for SGA) to measure the phase of the two signals Difference.

背景技术Background technique

科里奥利质量流量计(以下简称为科氏质量流量计)可以直接测量质量流量，是当前发展最为迅速的流量计之一，具有广阔的应用前景。科氏质量流量计要求其信号处理电路精确地测量来自两个流量传感器信号的相位差，并跟踪其频率的变化。基于模拟和数字电路的科氏质量流量传感器的信号处理方式和系统存在一些局限。为此，人们将数字信号处理方法和DSP应用于科氏质量流量传感器的信号处理，目前有以下三种方案。The Coriolis mass flowmeter (hereinafter referred to as the Coriolis mass flowmeter) can directly measure the mass flow rate. It is one of the most rapidly developing flowmeters and has broad application prospects. A Coriolis mass flow meter requires its signal processing circuitry to accurately measure the phase difference of the signals from the two flow sensors and track changes in their frequency. There are some limitations in the signal processing methods and systems of Coriolis mass flow sensors based on analog and digital circuits. To this end, people apply digital signal processing methods and DSP to the signal processing of Coriolis mass flow sensors. At present, there are the following three schemes.

(1)基于离散傅里叶变换的方法(1) Method based on discrete Fourier transform

美国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的信号处理系统，解决了美国专利中没有说明的技术难点，并在细测和频率跟踪方面做了改进。Paul Romano of the American Micro Motion company uses discrete Fourier transform (DFT) to process the output signal of the Coriolis mass flowmeter, and uses the DSP of the TMS series as the processing core of the secondary instrument ("Coriolicsmass flow rate meter having a substantially increased noise immunity" , US Patent No.4934196, Jun.19, 1990). When non-full-period sampling is performed, the calculation error of DFT cannot meet the requirements of instrument accuracy. For this reason, the ideas of rough measurement, fine measurement and frequency tracking are proposed. However, some key technologies are not disclosed. For example, when the frequency changes, how to collect zero-crossing points, etc. Referring to his ideas, Hefei University of Technology developed a signal processing system based on ADSP series DSP using DFT, which solved the technical difficulties not described in the US patent, and made improvements in fine measurement and frequency tracking.

此类方法存在的问题是：(a)实时性较差。在频率跟踪时，要不断地变化采样频率进行采样和计算，再比较功率谱值的大小，以确定实现整周期采样的频率，其时间长达10秒以上；(b)当输入信号中含有噪声时，频率计算的精度很高，但是，相位差的精度受到影响。其原因在于计算相位差时要用到幅度，而DFT计算出的幅度受噪声影响大。The problems of this type of method are: (a) poor real-time performance. In frequency tracking, it is necessary to continuously change the sampling frequency for sampling and calculation, and then compare the size of the power spectrum value to determine the frequency for realizing the whole cycle sampling, which lasts for more than 10 seconds; (b) when the input signal contains noise When , the accuracy of frequency calculation is high, but the accuracy of phase difference is affected. The reason is that the magnitude is used when calculating the phase difference, and the magnitude calculated by DFT is greatly affected by noise.

(2)基于信号幅值的方法(2) Method based on signal amplitude

日本富士公司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 0702212A2，20.03.1996)。此方法的问题是：两个传感器信号的相位差很小(一般小于4度)，故差动放大器的幅值很小，极易受到噪声干扰；补偿方法耗时太多，因为计算一次相位差需要采集6路信号。Yoshimura Hiroyuki of Fuji Corporation of Japan amplifies the signals from the two sensors in the Coriolis mass flowmeter, and sends the two signals to the differential amplifier at the same time to obtain the difference between the signals of the two sensors. The multiplexer sends these three signals to the analog/digital converter in sequence, and then sends them to the DSP for DFT to calculate the phase difference of the two sensor signals, and select one signal as a reference to compensate for the difference in the characteristics of each transmission channel.造成的误差(“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 0702212A2，20.03.1996)。 The problem with this method is: the phase difference of the two sensor signals is very small (generally less than 4 degrees), so the amplitude of the differential amplifier is very small, and it is very susceptible to noise interference; the compensation method takes too much time, because the calculation of the phase difference once Need to collect 6 signals.

(3)基于自适应线性增强的方法(3) Method based on adaptive linear enhancement

美国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值(窄陷波带)滤波器以产生比以往的技术或者上述第一实施例中更强的抑制噪声和谐波的作用。American Micro Motion Company Howard V.Derby et al. designed a DSP-based digital signal processing system, using Adaptive Linear Enhancement (ALE) technology to determine the frequency of the vibrating tube and the phase difference between the two signals, so as to more accurately measure Mass flow ("Method and apparatus for adaptive line enhancement in Coriolis mass flow meter measurement", US Patent No.5555190, Sep.10, 1996; "Adaptive linear enhancement method and device for Coriolis mass flow meter measurement" , Chinese Invention Patent Application Publication, CN 1190461, on August 12th, 1998. Because these two patents are exactly the same in content, hereinafter referred to as " the patent of Derby "). There are two embodiments in Derby's patent, and the signal processing link is composed of three parts: multi-decimation and one-filtering, frequency estimation/linear enhancement and phase difference (time difference) calculation. Multi-decimation and one-filtering adopt 8:1 and 6:1 two-stage decimation, and phase difference calculation adopts improved Goertzel algorithm. The frequency estimation/linear enhancement part is different, and the ALE of the first embodiment is realized by an adaptive notch filter (Adaptive Notch Filter, abbreviated as ANF) to two parameter estimates simultaneously, the ANF that this embodiment adopts ANF with minimal parameters (Arye Nehoral, "a minimal parameter adaptive notch filter with constrained poles and zeros", IEEE Transactions on Acoutics, Speech. and Signal Processing, Vol.ASSP-33, No.4, August 1985, pp.987-996 ) is not exactly the same, it relaxes the restriction on zero-pole pairs, and does not fix the zero-point radius to 1, but is automatically adjusted by the algorithm, which can take into account both the accuracy of frequency tracking and the convergence speed of the algorithm (that is, the speed of frequency tracking). Similarly, the second embodiment uses four adaptive notch filters, that is, two notch filters are connected in series on the left channel and the right channel respectively. The two ANFs placed in the left and right channels respectively are cascaded, where the first filter uses a low-Q (wide notch band) filter to produce limited signal enhancement, however, it converges quickly To the variation range of the fundamental frequency of the vibrating flow tube, the signal output from the first cascaded ANF is then transmitted to the second cascaded ANF; the second cascaded ANF uses a high-Q (narrow notch band) filter to produce a stronger effect of suppressing noise and harmonics than in the prior art or in the above-mentioned first embodiment.

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's patent also has limitations. In the first embodiment, two parameters are estimated, increasing the complexity of the algorithm. In the second embodiment, two cascaded ANFs are used, and the complexity of the algorithm is higher than that of the first embodiment; moreover, since the second-level ANF is based on the first-level ANF, only the first The second stage will converge after the first stage converges. Therefore, in the case of a large frequency change, the convergence speed will decrease instead. In addition, the Goertzel algorithm used to calculate the phase difference may overflow when implemented with a fixed point (see J.A.Beraldin and W.Steenaart, "Overflow analysis of a fixed-point implementation of the Goertzel algorithm, IEEE Trans. Circuits and Systems, Vol.36 , No.2, February 1989, pp.322-324).

因此，需要一种既能兼顾频率跟踪精度和算法收敛速度，又不明显增加算法复杂程度的信号处理方法，而且该方法用定点实现时不易溢出。Therefore, there is a need for a signal processing method that can balance frequency tracking accuracy and algorithm convergence speed without significantly increasing the complexity of the algorithm, and the method is not easy to overflow when implemented with fixed points.

发明内容Contents of the invention

本发明的目的在于解决Derby的专利中存在的缺点，又不失去其优点，即设计了一种信号处理方法和装置，它既能兼顾频率跟踪精度和算法的收敛速度，又不明显增加算法的复杂程度，而且用定点实现时不易溢出。The purpose of the present invention is to solve the shortcomings of Derby's patent without losing its advantages, that is, to design a signal processing method and device, which can not only take into account the frequency tracking accuracy and the convergence speed of the algorithm, but also not significantly increase the algorithm. Complexity, and it is not easy to overflow when implemented with fixed point.

本发明为了实现发明目的，采用了如下技术方案。信号处理环节由三部分组成：多抽一滤波、频率估计/线性增强和相位差(时间差)计算。In order to realize the object of the invention, the present invention adopts the following technical solutions. The signal processing link consists of three parts: multi-decimation and one-filtering, frequency estimation/linear enhancement and phase difference (time difference) calculation.

本发明的频率估计/线性增强部分与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对大的、突然的频率变化无法跟踪的缺点。The frequency estimation/linear enhancement part of the present invention is completely different from Derby's patent, and what adopts is a novel filter-Adaptive Funnel Shaped Filter (Adaptive Funnel Shaped Filter, abbreviated as AFF) (Sergio M.Savaresi, "Funnel filter: a New Class of Filters for Frequency Estimation of Harmonic Signals”, Automatic, Vol.33, No.9, September, 1997, pp.1711-1718). This filter is different from the quadratic adaptive notch filter, which is often used for frequency estimation and linear enhancement, and is characterized by only a single design parameter (Pole limiting factor or debiasing factor, represented by α in Derby's patent), the selection of this parameter usually needs to take into account the tracking ability and tracking accuracy. For ANF with only a single design parameter, the two are contradictory : (1) ρ is close to 1, the biggest advantage is that the variance and deviation of frequency estimation are very small, that is, the tracking accuracy is high; but the tracking ability has declined, it cannot track relatively large frequency changes, and it is sensitive to initial conditions; (2) When ρ is far from 1, the tracking ability is strong, it is not sensitive to initial conditions, and it can track relatively large frequency changes; but the variance and deviation of frequency estimation in steady state are large (that is, the tracking accuracy is not high). Therefore, when designing an ANF, a trade-off can only be made between tracking capability and tracking accuracy. However, AFF has two design parameters, which have stronger restrictions on zero and pole points, and its order is higher than that of ANF at the same level. For example, in the case of estimating a single parameter, ANF only needs 2 times, while AFF needs 4 times, but the algorithm The complexity does not increase much, because the parameter that needs to be estimated is only one. The characteristics of the AFF cost function are: wide and smooth mouth, steep and narrow bottom; compared with ANF, the wide place is wider than it, and the narrow place is narrower than it, so it can take into account both tracking ability and tracking accuracy. For large sudden frequency changes, its tracking ability is better than ANF, but for small slow frequency changes, the tracking performance is about the same. It has the characteristics of high tracking accuracy of ANF, and at the same time overcomes the disadvantage that ANF cannot track large and sudden frequency changes.

本发明在相位差计算部分采用更易于定点实现的、且非常适合于时变正弦信号的滑动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)非常适合于时变正弦信号。The present invention adopts the Sliding Goertzel Algorithm (Sliding Goertzel Algorithm, abbreviated as SGA) (Sliding Goertzel Algorithm, abbreviated as SGA) (Sliding Goertzel Algorithm, abbreviated as SGA) (Joe F.Chicharo, Mehdi T.Kilani, "A Sliding Goertzel algorithm", Signal Processing, Vol.52, No.3, August, 1996, pp.283-297), this algorithm has many advantages compared with traditional ones (including the improved Goertzel algorithm): (1) It can The Fourier coefficients are calculated within the number of points of a signal cycle, which has a faster acquisition time; (2) when implemented with a fixed-point algorithm, numerical overflow is not easy to occur; (3) it is very suitable for time-varying sinusoidal signals.

附图说明Description of drawings

图1是本发明的系统图示。Figure 1 is a system diagram of the present invention.

图2是DSP完成的主要功能。Figure 2 is the main function completed by DSP.

图3是本发明采用的信号处理方法Fig. 3 is the signal processing method that the present invention adopts

图4是ANF的代价函数的频率特性。Figure 4 is the frequency characteristics of the cost function of ANF.

图5是ANF与AFF代价函数频率特性比较。Figure 5 is a comparison of the frequency characteristics of ANF and AFF cost functions.

图6是数字信号处理部分的流程。Figure 6 is a flow chart of the digital signal processing section.

图7是频率估计/线性增强的自适应算法。Figure 7 is an adaptive algorithm for frequency estimation/linear enhancement.

图8是SGA流程图。Fig. 8 is a flowchart of SGA.

图9是SGA的实现框图。Figure 9 is a block diagram of the implementation of the SGA.

图10是整个跟踪过程图示。Figure 10 is a diagram of the entire tracking process.

图11是频率突变前的稳态情况。Figure 11 is the steady state situation before the frequency mutation.

图12是频率突变后的稳态情况。Figure 12 is the steady state situation after frequency mutation.

图13是自适应算法的两个设计参数的调整。Figure 13 shows the adjustment of two design parameters of the adaptive algorithm.

具体实施方式Detailed ways

下面结合附图和实施例对本发明作详细说明。The present invention will be described in detail below in conjunction with the accompanying drawings and embodiments.

科氏流量计是基于科里奥利力的原理而设计的。流体流过测量管时，如果测量管以某一频率振动，则振动的测量管相当于一个匀速转动的参考系。由于流体与测量管具有相对运动，所以会受到科里奥利力的作用。其反作用力作用在测量管的两边上，使测量管发生扭曲。测量管两边装着磁电传感器，输出信号是正比于测量管振动速度的电压信号。当振动管是以一定频率振动时，其角速度按正弦规律变化。故由磁电传感器输出的信号是一正弦信号，其频率为角速度的变化频率，大小正比于角速度。测出其两路信号的相位差和信号频率或者时间差就可以求出流体的质量流量。Coriolis flow meters are designed based on the principle of Coriolis force. When the fluid flows through the measuring tube, if the measuring tube vibrates at a certain frequency, the vibrating measuring tube is equivalent to a reference frame rotating at a constant speed. Due to the relative motion of the fluid and the measuring tube, it will be affected by Coriolis force. Its reaction force acts on both sides of the measuring tube, causing the measuring tube to twist. The two sides of the measuring tube are equipped with magnetoelectric sensors, and the output signal is a voltage signal proportional to the vibration speed of the measuring tube. When the vibrating tube vibrates at a certain frequency, its angular velocity changes according to the sinusoidal law. Therefore, the signal output by the magnetoelectric sensor is a sinusoidal signal whose frequency is the change frequency of the angular velocity, and its magnitude is proportional to the angular velocity. The mass flow rate of the fluid can be obtained by measuring the phase difference and signal frequency or time difference of the two signals.

图1是本发明的系统图示。两路传感器信号经过信号调理和A/D转换后，进入DSP。信号调理部分包括低通滤波和放大。由于从传感器出来的信号幅值比较小，而且有用的正弦信号落入许多工业噪声范围内，为了充分利用A/D转换器的量程，需要对从传感器出来的信号进行放大；同时为了尽可能地削弱工业噪声的干扰，还采用了低通滤器器。低通滤波和放大环节的设计是本领域的技术人员所熟知的。本发明选用巴特沃斯二阶低通滤波器。Figure 1 is a system diagram of the present invention. The two sensor signals enter the DSP after signal conditioning and A/D conversion. The signal conditioning section includes low-pass filtering and amplification. Since the signal amplitude from the sensor is relatively small, and the useful sinusoidal signal falls into the range of many industrial noises, in order to make full use of the range of the A/D converter, it is necessary to amplify the signal from the sensor; To reduce the interference of industrial noise, a low-pass filter is also used. The design of low-pass filtering and amplification links is well known to those skilled in the art. The present invention selects Butterworth second-order low-pass filter.

A/D的采样频率固定为38.4KHz。多抽一部分采用Derby的专利所描述的抽选率和滤波器结构。另外还可以采用48KHz的采样频率以及12∶1和6∶1两级多抽一滤波。由于不只一路模拟输入信号，可以采用每路信号一个A/D；或者多路信号经过多路选通后再经过一个A/D。The sampling frequency of A/D is fixed at 38.4KHz. The decimation rate and filter structure described in Derby's patent are used for extra decimation. In addition, 48KHz sampling frequency and 12:1 and 6:1 two-stage multi-decimation and one-filtering can also be used. Since there are more than one analog input signal, one A/D for each signal can be used; or multiple signals can pass through an A/D after multi-channel selection.

DSP是信号处理系统的核心，它对两路传感器信号进行数字滤波、频率估计/线性增强和相位差(时间差)计算。对DSP的选择也没有特别限制，几乎可以用任何市面上可以买到的DSP芯片。本发明所叙述的软件由于采用了比Derby的专利更适合定点实现的算法，因此采用TI公司的定点DSP芯片。DSP is the core of the signal processing system, which performs digital filtering, frequency estimation/linear enhancement and phase difference (time difference) calculation on the two sensor signals. There is no special restriction on the choice of DSP, and almost any DSP chip available on the market can be used. The software described in the present invention adopts the fixed-point DSP chip of TI Company because it adopts an algorithm that is more suitable for fixed-point realization than Derby's patent.

DSP还实现数字驱动的功能。DSP also realizes the function of digital drive.

本系统还包括LCD显示、键盘输入电路。The system also includes LCD display, keyboard input circuit.

以下的描述是针对典型的科氏流量计应用进行的，在这种流量计中振动流管的基频大约是100Hz。很容易理解，本发明的装置和方法可以应用于处理任何常规的流量计振动基频。The following description is for a typical Coriolis flowmeter application where the fundamental frequency of the vibrating flow tube is around 100 Hz. It will be readily understood that the apparatus and method of the present invention can be applied to address any conventional flow meter vibration fundamental frequency.

图2是本系统中DSP完成的主要功能图示。DSP主要完成信号采集、信号处理和数字驱动。系统上电后，DSP首先在一段时间内(定时时间)使振动管振动起来，然后进入正常的工作状态。当定时时间没到时，执行振动管起振算法：DSP先施加一带限随机信号，经过D/A转换、功率放大后施加给振动管(可以连续施加一小段时间)，振动系统就能够振动起来(虽然振动不是很理想)。之所以施加随机噪声，是因为随机噪声的频谱很宽，具有很多种频率分量；又由于振动系统具有选频特性，其中，噪声中频率与振动管固有频率相同的分量就会起主要作用，振动管开始振动起来。定时时间到后，振动已经趋于稳定，即振动管以固有频率及稳定的幅值振动，如果振动还没有稳定(这可以通过比较DSP连续的几个信号处理周期计算的信号幅值来判断)，重新设定定时时间，再执行起振算法，直到振动稳定为止。这时DSP进入正常的工作状态：读入传感器信号、对信号进行处理(去噪、求频率和相位)、产生驱动信号(用求出来的频率、相位构造新的正弦驱动信号)。关于定时时间，可以设定为1～5s，可以根据实际情况具体选定。Fig. 2 is a schematic diagram of the main functions completed by DSP in this system. DSP mainly completes signal acquisition, signal processing and digital drive. After the system is powered on, the DSP first makes the vibrating tube vibrate for a period of time (timing time), and then enters the normal working state. When the timing is not up, execute the vibrating tube start-up algorithm: DSP first applies a band-limited random signal, after D/A conversion and power amplification, it is applied to the vibrating tube (it can be applied continuously for a short period of time), and the vibrating system can vibrate (Although the vibration is not ideal). The reason why random noise is applied is that the frequency spectrum of random noise is very wide and has many kinds of frequency components; and because the vibration system has a frequency-selective characteristic, the component in the noise with the same frequency as the natural frequency of the vibrating tube will play a major role, and the vibration The tube began to vibrate. After the timing time is up, the vibration has tended to be stable, that is, the vibrating tube vibrates at a natural frequency and a stable amplitude. If the vibration is not stable (this can be judged by comparing the signal amplitude calculated by the DSP for several consecutive signal processing cycles) , reset the timing time, and then execute the vibration algorithm until the vibration is stable. At this time, the DSP enters the normal working state: read in the sensor signal, process the signal (denoise, find the frequency and phase), generate the drive signal (construct a new sinusoidal drive signal with the found frequency and phase). Regarding the timing time, it can be set to 1-5s, which can be selected according to the actual situation.

图3是本发明所采用的信号处理方法。经过转换的数字信号通过路径300(301)传输到抽选单元，抽选单元分为两级。第一级抽选为8∶1，使采样频率从38.4KHz(记为f_s1)，降低到4.8KHz(记为f_s2)，采用的滤波器的传函为：Fig. 3 is a signal processing method adopted by the present invention. The converted digital signal is transmitted to the sampling unit through the path 300 (301), and the sampling unit is divided into two stages. The first-stage decimation is 8:1, so that the sampling frequency is reduced from 38.4KHz (denoted as f _s1 ) to 4.8KHz (denoted as f _s2 ), and the transfer function of the filter used is:

$G G (({z z}^{- - 11})) = = \frac{{(({11 - - z z}^{- - 88}))}^{55}}{{(({11 - - z z}^{- - 11}))}^{55}} - - - - - - - - - - - - ((11))$

零极点对消后得到一个36抽头的FIR滤波器，这种滤波器在二次采样频率的各个倍数点具有5个零点，这可以极大的消除混叠在第二级滤波器通带中的那些频率；该滤波器具有最小整数系数，可以用一种精确的计算机运算来表示，从而简化卷积运算的复杂性和提高计算速度。After zero-pole cancellation, a 36-tap FIR filter is obtained. This filter has 5 zero points at each multiple of the sub-sampling frequency, which can greatly eliminate aliasing in the passband of the second-stage filter. those frequencies; the filter has the smallest integer coefficients and can be represented by an exact computer operation, thereby simplifying the complexity of the convolution operation and increasing the calculation speed.

第二级抽选为6∶1，使采样频率从4.8KHz降低到800Hz(记为f_s)，采用的滤波器是Remez转换算法设计的一个131抽头的FIR滤波器。通带为直流到250Hz，阻带起点为400Hz；通带加权为10^-5，阻带加权为1。The second stage of decimation is 6:1, which reduces the sampling frequency from 4.8KHz to 800Hz (marked as f _s ). The filter used is a 131-tap FIR filter designed by Remez conversion algorithm. The passband is from DC to 250Hz, and the start of the stopband is 400Hz; the weight of the passband is 10 ^-5 , and the weight of the stopband is 1.

利用两级抽选滤波器具有较高的抗混迭性能。信号经过两级多抽一滤波后通过路径304、305进入自适应漏斗型滤波器(Adaptive Funnel filter，缩写为AFF)环节，该环节完成频率估计/线性增强的功能，它通过自适应算法对输入信号的频率进行估计，同时滤去混在正弦信号中的宽带噪声，输出增强信号；增强后的信号进入滑动Goertzel算法(Sliding Goertzel Algorithm，缩写为SGA)环节，该环节计算增强信号的离散傅立叶系数，从而求出信号的相位(差)；相位差、频率都求出后，就可以求出两路信号的时间差。The use of two-stage decimation filter has high anti-aliasing performance. The signal enters the adaptive funnel filter (Adaptive Funnel filter, abbreviated as AFF) link through the paths 304 and 305 after being filtered by two stages, and this link completes the function of frequency estimation/linear enhancement. The frequency of the signal is estimated, and the broadband noise mixed in the sinusoidal signal is filtered out, and the enhanced signal is output; the enhanced signal enters the sliding Goertzel algorithm (Sliding Goertzel Algorithm, abbreviated as SGA) link, which calculates the discrete Fourier coefficient of the enhanced signal, Thus, the phase (difference) of the signal is obtained; after the phase difference and the frequency are obtained, the time difference of the two signals can be obtained.

下面描述频率计算/线性增强部分。The frequency calculation/linear enhancement section is described below.

设两路传感器信号经过两级多抽一后可以表示为：Assuming that the two-way sensor signals can be expressed as:

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

y_R(n)＝Asin(2πf₀nT+θ₂)+ξ₂(n) (3)y _R (n)＝Asin(2πf ₀ nT+θ ₂ )+ξ ₂ (n) (3)

$T T = = \frac{11}{{f f}_{s the s}} . . - - - - - - ((44))$

式中，y_L，y_R分别表示左、右两路信号，A表示信号幅值，f₀是信号的真实频率，θ₁，θ₂分别为两路信号的相位(单位：弧度)，ξ₁，ξ₂分别为两路信号经过两级多抽一滤波后残留在有用信号中的噪声，f_s是经过两级多抽一后的采样频率，T是采样间隔。In the formula, y _L and y _R represent the left and right signals respectively, A represents the signal amplitude, f ₀ is the real frequency of the signal, θ ₁ and θ ₂ are the phases of the two signals respectively (unit: radian), ξ ₁ and _ξ2 are respectively the noise remaining in the useful signal after two-stage decimation filtering of the two signals, f _s is the sampling frequency after two-stage decimation, and T is the sampling interval.

本发明采用的自适应算法就是要估计出信号的频率f₀，并且当信号频率发生变化时要能够跟踪信号频率的变化，同时对信号进行信号增强(即有效地滤出噪声分量ξ₁和ξ₂)。The adaptive algorithm adopted in the present invention is to estimate the frequency f ₀ of the signal, and when the signal frequency changes, it will be able to track the change of the signal frequency, and simultaneously carry out signal enhancement to the signal (i.e. effectively filter out the noise components ξ ₁ and ξ ₂ ).

本发明通过AFF来实现频率计算/线性增强功能。AFF的一般形式为The present invention realizes the function of frequency calculation/linear enhancement through AFF. The general form of AFF is

$H h (({z z}^{- - 11})) = = \frac{11}{22} [[\frac{C C (({z z}^{- - 11}))}{{D D.}_{11} (({z z}^{- - 11}))} + + \frac{C C (({z z}^{- - 11}))}{{D D.}_{22} (({z z}^{- - 11}))}]] - - - - - - - - - - - - - - - - - - ((55))$

其中 $\frac{C (z^{- 1})}{D_{1} (z^{- 1})}, 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)决定。in $\frac{C (z^{- 1})}{{D.}_{1} (z^{- 1})}, i = 1,2$ are the transfer functions of the two notch filters, so it can be considered that AFF is derived from ANF. Since there are many forms of notch filters, correspondingly, there are many forms of AFF. The present invention adopts the AFF derived from the theoretically unbiased form of ANF (called Type II ANF) (called Type II AFF, see Sergio M. Savaresi, Funnel filter: a New Class of Filters for Frequency Estimation of Harmonic Signals, Automatic , Vol.33, No.9, September, 1997, the structure given in pp.1711-1718). It should be pointed out that other forms of AFF, such as AFF derived from the direct form of ANF (called Type I ANF) (called Type I AFF), can also achieve this function, and different forms have different advantages. shortcoming. The AFF used in the present invention is a Type II AFF, and its transfer function is determined by formula (6) to formula (8).

C(z^-1)＝1+2az^-1+z^-2 (6)C(z ^-1 )＝1+2az ^-1 +z ^-2 (6)

${D D.}_{11} (({z z}^{- - 11})) = = 11 + + (({11 + + ρ ρ}_{11}^{22})) {az az}^{- - 11} + + {ρ ρ}_{11}^{22} {z z}^{- - 22} - - - - - - - - - - - - - - - - ((77))$

${D D.}_{22} (({z z}^{- - 11})) = = 11 + + ((11 + + {ρ ρ}_{22}^{22})) a a {z z}^{- - 11} + + {ρ ρ}_{22}^{22} {z z}^{- - 22} - - - - - - - - - - - - - - ((88))$

${\overset{^^}{Ω Ω}}_{00} = = {cos cos}^{- - 11} ((- - a a)) - - - - - - - - - - - - - - - - ((99))$

${\overset{^^}{f f}}_{00} = = \frac{\overset{^^}{Ω Ω}}{22 π π} {f f}_{s the s} - - - - - - - - - - - - - - - - ((1010))$

式中，a是待估计的参数，是待估计的信号频率(数字频率，单位：弧度)，a与之间有固定的函数关系，即只要估计出参数a，就可以估计出信号的频率。

是信号的频率的估计值(单位：Hz)。自适应算法以使式(11)最小为准则调整系数a。In the formula, a is the parameter to be estimated, is the signal frequency to be estimated (digital frequency, unit: radian), a and There is a fixed functional relationship between them, that is, as long as the parameter a is estimated, the frequency of the signal can be estimated.

is an estimate of the frequency of the signal (unit: Hz). The self-adaptive algorithm adjusts the coefficient a based on the minimum of formula (11).

$((Ω Ω)) = = {Σ Σ}_{s the s = = 11}^{n no} {λ λ}^{n no - - s the s} {ϵ ϵ ((s the s,, Ω Ω))}^{22} - - - - - - - - - - - - - - ((1111))$

式(11)称为代价函数(Cost Function)。其中Equation (11) is called the cost function (Cost Function). in

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

是两个极点限制因子，是我们可以利用的设计参数。

is the two-pole limiting factor and is a design parameter we can exploit.

当ρ₁＝ρ₂时AFF就完全等同于ANF，从这个角度讲，ANF是AFF的特例。由于AFF有两个设计参数，因而算法上就可以采取比ANF更为灵活的方法，既考虑跟踪能力又兼顾跟踪精度。When ρ ₁ =ρ ₂ , AFF is completely equivalent to ANF. From this perspective, ANF is a special case of AFF. Since AFF has two design parameters, the algorithm can adopt a more flexible method than ANF, considering both tracking ability and tracking accuracy.

图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时所对应的频率点)来描述吸引域，根据这样的定义，有下面的结论(考虑相同的估计方差的情况下， $ρ^{2} = 0.98, ρ_{1}^{2} = 0.84, ρ_{2}^{2} = 0.988) :$ Figure 4 is the frequency characteristics of the ANF cost function. The frequency characteristics of the cost function of the notch filter are divided into five regions: a1, a small flat region near the notch frequency; a2, a large flat region to the left of the notch frequency; a3, a large flat region to the right of the notch frequency ; b1, the steep slope area on the left side of the notch frequency; b2, the steep slope area on the right side of the notch frequency. There will be a lot of ripples in the a2 and a3 areas, so in the area far from the notch frequency, the algorithm for finding the minimum value of the cost function will easily converge to a local minimum, which is not desirable. It is natural to take the width of the notch as a measure of the domain of attraction. Literature (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) uses Ω ⁺ -Ω ^- (Ω ⁺ is The corresponding frequency point when the slope of the cost function is +1, Ω ^- is the corresponding frequency point when the slope of the cost function is -1) to describe the attraction domain, according to this definition, the following conclusions (considering the same estimate In the case of variance, $ρ^{2} = 0.98, ρ_{1}^{2} = 0.84, ρ_{2}^{2} = 0.988) :$

可见AFF的吸引域比ANF大很多，因此其跟踪能力也就比ANF强很多，能跟踪上比较大的频率变化。It can be seen that the attraction field of AFF is much larger than that of ANF, so its tracking ability is much stronger than that of ANF, and it can track relatively large frequency changes.

图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算法的有关公式，同时，经过大量的仿真确定了两个设计参数ρ₁ ²，ρ₂ ²的选择。Figure 5 is a comparison of the frequency characteristics of AFF and ANF cost functions. From the comparison of their cost functions, it can be seen that the tracking ability of AFF is stronger than that of ANF (the mouth is wider than it), and the tracking accuracy is similar to it (the bottom is almost as narrow) ; It can also be seen from the figure that the adaptive algorithm of ANF may converge to a local minimum, while AFF can avoid it. It should be pointed out that the present invention has only adopted literature (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 ) proposed the AFF structure, and derived the relevant formulas suitable for the RML algorithm. At the same time, after a large number of simulations, the selection of the two design parameters ρ ₁ ² and ρ ₂ ² was determined.

借鉴Derby的专利，采用SNR故障检测。自适应算法中的输入信号是经过两级多抽一后的左、右路信号，用符号y_L和y_R表示。Using Derby's patent for reference, SNR fault detection is adopted. The input signal in the self-adaptive algorithm is the left and right channel signals after two stages of multi-pumping, represented by symbols y _L and y _R.

左信道和右信道的输出信号都可以被用作加权自适应单元的反馈信号。虽然在加权自适应单元中同时使用两信道的输出信号作反馈信号是可能的，但是这样所产生的益处与增加的计算复杂程度相比并不占优，因此只用其中一个信道(比如右信道)作反馈，而计算出的加权自适应参数被同时用在左右信道的ANF，从而使两个传感器输出信道经过同样的处理。在本发明所采用的算法中，用右信道的信号作为加权自适应的反馈信号，由该反馈信号确定好滤波器的系数后，两路信号又通过相同的滤波器(AFF)，得到增强后的信号，用x_L和x_R表示。在下面的公式中，没有严格区分y_L和y_R，而是用y表示；也没有严格区分x_L和x_R，而是用x表示。Both the output signals of the left channel and the right channel can be used as the feedback signal of the weight adaptation unit. Although it is possible to use the output signals of two channels simultaneously as feedback signals in the weighted adaptive unit, the benefits produced by this are not dominant compared with the increased computational complexity, so only one of the channels (such as the right channel ) as feedback, and the calculated weighted adaptive parameters are used in the ANF of the left and right channels at the same time, so that the two sensor output channels undergo the same processing. In the algorithm adopted in the present invention, the signal of the right channel is used as a weighted adaptive feedback signal. After the coefficient of the filter is determined by the feedback signal, the two signals pass through the same filter (AFF) again to obtain an enhanced The signal, denoted by x _L and x _R. In the following formula, y _L and y _R are not strictly distinguished, but are represented by y; and x _L and x _R are not strictly distinguished, but are represented by x.

图7是频率估计/线性增强部分的自适应算法(AFF环节)的流程图，本发明采用的AFF的自适应算法的主要公式为：Fig. 7 is the flow chart of the adaptive algorithm (AFF link) of frequency estimation/linear enhancement part, and the main formula of the adaptive algorithm of AFF that the present invention adopts is:

(1)有关向量、矩阵及符号的意义：(1) The meanings of vectors, matrices and symbols:

输入信号：y(包括y_L，y_R)Input signal: y (including y _L , y _R )

先验预测误差：ε(n)Prior prediction error: ε(n)

后验预测误差： ε(n)Posterior prediction error: ε(n)

设计参数(极点限制因子和去偏置因子)：ρ₁ ²，ρ₂ ² Design parameters (pole limiting factor and debiasing factor): ρ ₁ ² , ρ ₂ ²

遗忘因子：λForgetting factor: λ

中间状态变量：F(n)Intermediate state variable: F(n)

自回归向量：Φ(n)＝[φ₁(n)，φ₂(n)]^T；Autoregressive vector: Φ(n)=[φ ₁ (n), φ ₂ (n)] ^T ;

负梯度向量：Ψ(n)＝[₁(n)，₂(n)]^T；Negative gradient vector: Ψ(n)=[ ₁ (n),  ₂ (n)] ^T ;

待估计参数a(t)构成的向量：θ(n)＝[a(n)，a²(n)]^T；A vector composed of parameters to be estimated a(t): θ(n)=[a(n), a ² (n)] ^T ;

协方差矩阵： $P (n) = [\begin{matrix} p_{1} (n), p_{2} (n) \\ p_{3} (n), p_{4} (n) \end{matrix}]$ Covariance matrix: $P (no) = [\begin{matrix} p_{1} (no), p_{2} (no) \\ p_{3} (no), p_{4} (no) \end{matrix}]$

中间状态矩阵：M(n)Intermediate state matrix: M(n)

增强后的信号：x(包括x_L，x_R)Enhanced signal: x (including x _L , x _R )

数字角频率：Ω、

Ω_k等Digital angular frequency: Ω,

Ω _k etc.

正弦信号频率估计值：

Estimated frequency of a sinusoidal signal:

(2)初始化：(2) Initialization:

θ(0)＝0，P(0)＝10^-2I，P_SNR＝10^-4I，Ψ(0)＝Φ(0)＝0；y(-i)＝0，i＝1，2θ(0)=0, P(0)=10 ^-2 I, P _SNR =10 ^-4 I, Ψ(0)=Φ(0)=0; y(-i)=0, i=1, 2

${ρ ρ}_{1010}^{22} = = {0.7 0.7}^{22},, {ρ ρ}_{11 \infty \infty}^{22} = = {0.995 0.995}^{22},, {ρ ρ}_{11 decay decay} = = 0.998 0.998,, {ρ ρ}_{11 SNR SNR} = = 0.80 0.80$

${ρ ρ}_{2020}^{22} = = {0.8 0.8}^{22},, {ρ ρ}_{22 \infty \infty}^{22} = = {0.995 0.995}^{22},, {ρ ρ}_{22 decay decay} = = 0.99 0.99,, {ρ ρ}_{22 SNR SNR} = = 0.95 0.95$

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

(3)主循环：(3) Main loop:

$F f ((n no)) = = y the y ((n no)) + + \frac{22 + + {ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22}}{22} y the y ((n no - - 22)) + + \frac{{ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22}}{22} + + y the y ((n no - - 44)) - - (({ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22})) ϵ ϵ ((n no - - 22)) - - {ρ ρ}_{11}^{22} {ρ ρ}_{22}^{22} ϵ ϵ ((n no - - 44)) - - - - - - - - - - - - - - - - ((1414))$

ε(n)＝F(n)-Ψ^T(n)θ(n-1) (15)ε(n)＝F(n) ^-ΨT (n)θ(n-1) (15)

$P P ((n no)) = = \frac{11}{λ λ ((n no))} [[P P ((n no - - 11)) - - \frac{P P ((n no - - 11)) Ψ Ψ ((n no)) {Ψ Ψ}^{T T} ((n no)) P P ((n no - - 11))}{λ λ ((n no)) + + {Ψ Ψ}^{T T} ((n no)) P P ((n no - - 11)) Ψ Ψ ((n no))}]] - - - - - - - - - - ((1616))$

θ(n)＝θ(n-1)+P(n)Ψ(n)ε(n) (16)θ(n)=θ(n-1)+P(n)Ψ(n)ε(n) (16)

ε(n)＝F(n)-Ψ^T(n)θ(n) (17)ε(n)=F(n) ^-ΨT (n)θ(n) (17)

x(n)＝y(n)- ε(n)x(n)=y(n)- ε(n)

${φ φ}_{11} ((n no)) = = - - \frac{66 + + {ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22}}{22} y the y ((n no - - 11)) - - \frac{22 + + 33 {ρ ρ}_{11}^{22} + + 33 {ρ ρ}_{22}^{22}}{22} y the y ((n no - - 33)) + + ((22 + + {ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22})) ϵ ϵ ((n no - - 11)) - - - - - - - - ((1818))$

$+ + (({ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22} + + 22 {ρ ρ}_{11}^{22} {ρ ρ}_{22}^{22})) ϵ ϵ ((n no - - 33))$

${φ φ}_{22} ((n no)) = = - - ((22 + + {ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22})) y the y ((n no - - 22)) + + ((11 + + {ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22} + + {ρ ρ}_{11}^{22} {ρ ρ}_{22}^{22})) ϵ ϵ ((n no - - 22)) - - - - - - - - ((1919))$

Φ(n)＝[φ₁(n)，φ₂(n)]^T (20)Φ(n)＝[φ ₁ (n), φ ₂ (n)] ^T (20)

$M m ((n no)) = = [[((22 + + {ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22})) Ψ Ψ ((n no - - 11)) + + (({ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22} + + 22 {ρ ρ}_{11}^{22} {ρ ρ}_{22}^{22})) Ψ Ψ ((n no - - 33)),, ((11 + + {ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22} + + {ρ ρ}_{11}^{22} {ρ ρ}_{22}^{22})) Ψ Ψ ((n no - - 22))]] - - - - - - - - - - - - ((21 twenty one))$

$Ψ Ψ ((n no)) = = Φ Φ ((n no)) - - M m ((n no)) θ θ ((n no)) - - (({ρ ρ}_{11}^{22} + + {ρ ρ}_{22}^{22})) Ψ Ψ ((n no - - 22)) - - {ρ ρ}_{11}^{22} {ρ ρ}_{22}^{22} Ψ Ψ ((n no - - 44)) - - - - - - - - - - - - ((22 twenty two))$

${ρ ρ}_{11}^{22} ((n no)) = = {ρ ρ}_{11}^{22} ((n no - - 11)) {ρ ρ}_{11 decay decay} + + ((11 - - {ρ ρ}_{11 decay decay})) {ρ ρ}_{11}^{22} ((\infty \infty)) - - - - - - - - ((23 twenty three))$

${ρ ρ}_{22}^{22} ((n no)) = = {ρ ρ}_{22}^{22} ((n no - - 11)) {ρ ρ}_{22 decay decay} + + ((11 - - {ρ ρ}_{22 decay decay})) {ρ ρ}_{22}^{22} ((\infty \infty)) - - - - - - - - - - - - ((24 twenty four))$

λ(n)＝λ(n-1)λ_decay+(1-λ_decay)λ(∞) (25)λ(n)=λ(n-1)λ _decay +(1-λ _decay )λ(∞) (25)

频率的估计采用的公式为The frequency is estimated using the formula

$\overset{^^}{Ω Ω} ((n no)) = = {cos cos}^{- - 11} ((- - a a ((n no)))) - - - - - - - - - - - - ((2626))$

$\overset{^^}{f f} ((n no)) = = \frac{\overset{^^}{Ω Ω} ((n no))}{22 π π} {f f}_{s the s} - - - - - - - - - - - - - - - - - - ((2727))$

加权自适应单元计算所得值(即a(t)的值)用于计算相位和在SGA环节中计算相位差、时间差Δt。频率计算单元计算频率并产生Goertzel滤波器的加权信息。相位计算单元从频率计算单元接受Goertzel加权信息和频率信息。相位计算单元利用具有两个Hanning窗口的傅立叶分析技术确定滤波信号的相位选择8个流管周期作为优选的窗口长度。假设一个给定的频率期望值为f_q，则优选的窗口长度(L＝2N)由下式确定：The value calculated by the weighted adaptive unit (that is, the value of a(t)) is used to calculate the phase and calculate the phase difference and time difference Δt in the SGA link. The frequency calculation unit calculates frequencies and generates weighting information of the Goertzel filter. The phase calculation unit receives Goertzel weighting information and frequency information from the frequency calculation unit. The phase calculation unit uses the Fourier analysis technique with two Hanning windows to determine the phase of the filtered signal and selects 8 flow tube periods as the optimal window length. Assuming a given frequency expectation value f _q , the preferred window length (L=2N) is determined by the following formula:

$L L = = 22 N N = = 88 floor floor ((\frac{{f f}_{s the s}}{{f f}_{q q}})) - - - - - - - - - - - - ((2828))$

式中，floor(x)表示不大于x的整数。In the formula, floor(x) represents an integer not greater than x.

f_s为两级多抽一后的频率，比如f_s＝800，f_q＝100，则2N＝64。Hanning窗表示为：f _s is the frequency after two stages of extra sampling, for example, f _s =800, f _q =100, then 2N=64. The Hanning window is expressed as:

$h h ((n no)) = = \frac{11}{22} ((11 - - cos cos \frac{22 πn πn}{22 N N - - 11})),, n no = = 0,1 0,1,, \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot,, 22 N N - - 11 - - - - - - - - - - - - - - ((2929))$

为了尽可能最大限度的利用时域数据，本发明借鉴Derby的专利采用的重叠Hanning窗口的平行计算。In order to utilize the time domain data as much as possible, the present invention refers to the parallel calculation of overlapping Hanning windows adopted in the patent of Derby.

增强后的信号经过Hanning窗后，采用Goertzel算法计算其离散傅立叶系数。采用重叠Hanning窗口的平行计算后每个半窗长(每N个经过两级多抽一后的时域点)进行一次相位、频率和Δt的计算。此时频率的计算采用N个估计值的平均值求得：After the enhanced signal passes through the Hanning window, its discrete Fourier coefficients are calculated using the Goertzel algorithm. The phase, frequency and Δt are calculated once for each half-window length (every N time-domain points after two-level decimation) after the parallel calculation of overlapping Hanning windows. At this time, the calculation of the frequency is obtained by using the average value of N estimated values:

${\overset{^^}{Ω Ω}}_{00} ((n no)) = = {cos cos}^{- - 11} ((- - \frac{11}{N N} {Σ Σ}_{i i = = 11}^{N N} a a ((i i)))) - - - - - - - - - - - - - - ((3030))$

${\overset{^^}{f f}}_{00} ((n no)) = = \frac{{\overset{^^}{Ω Ω}}_{00} ((n no))}{22 π π} {f f}_{s the s} - - - - - - - - - - - - - - - - - - - - ((3131))$

下面描述相位差计算部分。The phase difference calculation section is described below.

本发明采用滑动Goertzel算法(Joe F.Chicharo，Mehdi T.Kilani，A Sliding Goertzelalgorithm，signal processing，Vol.52，No.3，August，1996，pp.283-297)，用该算法分别计算两路增强后的信号的离散傅立叶系数，根据离散傅立叶系数求出两路信号的相位，再求出两路信号之间的相位差，最后计算出两路信号之间的时间差。由于AFF环节已经估计出了信号的频率并同时对信号进行了线性增强，采用Goertzel类的算法可以直接计算给定频率点的离散傅立叶系数，从而大大减少了计算量(与传统的计算离散傅立叶系数的DFT算法、FFT算法相比较)。The present invention adopts the sliding Goertzel algorithm (Joe F.Chicharo, Mehdi T.Kilani, A Sliding Goertzel algorithm, signal processing, Vol.52, No.3, August, 1996, pp.283-297), calculates two paths respectively with this algorithm The discrete Fourier coefficient of the enhanced signal is used to obtain the phase of the two signals according to the discrete Fourier coefficient, and then obtain the phase difference between the two signals, and finally calculate the time difference between the two signals. Since the AFF link has estimated the frequency of the signal and linearly enhanced the signal at the same time, using the Goertzel algorithm can directly calculate the discrete Fourier coefficient of a given frequency point, thereby greatly reducing the amount of calculation (compared with the traditional calculation of the discrete Fourier coefficient Compared with the DFT algorithm and FFT algorithm).

图8是SGA的流程图，与此相对应，下面给出本发明采用的SGA的有关步骤和公式：Fig. 8 is the flowchart of SGA, and corresponding to this, provide the relevant steps and the formula of the SGA that the present invention adopts below:

(1)确定需要计算的频率点 (1) Determine the frequency points that need to be calculated

(2)状态条件初始化为零，即v_k(0)＝0；v_k(-i)＝0，k＝1，2(2) The state condition is initialized to zero, that is, v _k (0)=0; v _k (-i)=0, k=1, 2

(3)计算共振滤波器的的输出(3) Calculate the output of the resonance filter

${v v}_{k k} ((n no)) = = x x ((n no)) + + 22 cos cos {\overset{^^}{Ω Ω}}_{00} \cdot &Center Dot; {v v}_{k k} ((n no - - 11)) - - {v v}_{k k} ((n no - - 22)) - - - - - - - - - - ((3232))$

${\overset{^^}{A A}}_{k k} = = \frac{22 (({v v}_{k k} ((N N)) - - {v v}_{k k} ((N N - - 11)) cos cos {\overset{^^}{Ω Ω}}_{00}}{N N},, k k = = 1,2 1,2 - - - - - - - - - - - - - - ((3333))$

${\overset{^^}{B B}}_{k k} = = - - \frac{22 {v v}_{k k} ((N N - - 11)) sin sin {\overset{^^}{Ω Ω}}_{00}}{N N},, k k = = 1,2 1,2 - - - - - - - - - - - - ((3434))$

$[\begin{matrix} {a a}_{k k} \\ {b b}_{k k} \end{matrix}] = = [\begin{matrix} cos cos {\overset{^^}{Ω Ω}}_{00} n no,, - - sin sin {\overset{^^}{Ω Ω}}_{00} n no \\ sin sin {\overset{^^}{Ω Ω}}_{00} n no,, cos cos {\overset{^^}{Ω Ω}}_{00} n no \end{matrix}] [\begin{matrix} {\overset{^^}{A A}}_{k k} ((n no)) \\ {\overset{^^}{B B}}_{k k} ((n no)) \end{matrix}],, k k = = 1,2 1,2 - - - - - - - - - - ((3535))$

(4)计算相位(4) Calculate the phase

${\overset{^^}{θ θ}}_{k k} = = {tan the tan}^{- - 11} \frac{{a a}_{k k}}{{b b}_{k k}} - - - - - - - - - - - - - - ((3636))$

(5)两路信号的相位差(5) The phase difference of the two signals

$\overset{^^}{φ φ} = = {\overset{^^}{θ θ}}_{11} - - {\overset{^^}{θ θ}}_{22} - - - - - - - - - - - - - - ((3737))$

(6)计算时间差(6) Calculate the time difference

$Δt Δt = = \frac{\overset{^^}{φ φ}}{{\overset{^^}{Ω Ω}}_{00}} \cdot &Center Dot; \frac{11}{{f f}_{s the s}} - - - - - - - - - - - - ((3838))$

图9是实现SGA算法的结构图。Fig. 9 is a structural diagram for realizing the SGA algorithm.

需要指出的是，由于自适应算法估计出参数α，因而在上面的式(32)～式(34)中实际上并不需要计算余弦值和正弦值而可以分别采用式(39)、式(40)代替：It should be pointed out that since the parameter α is estimated by the adaptive algorithm, it is not actually necessary to calculate the cosine value in the above equations (32) to (34) and the sine Instead, formula (39) and formula (40) can be used to replace:

$cos cos {\overset{^^}{Ω Ω}}_{00} = = {W W}_{11} - - - - - - - - - - - - ((3939))$

$sin sin {\overset{^^}{Ω Ω}}_{00} = = {W W}_{22} - - - - - - - - - - - - - - ((4040))$

其中，in,

${W W}_{11} = = - - \frac{11}{N N} {Σ Σ}_{i i = = 11}^{N N} a a ((i i)) - - - - - - - - - - - - - - ((4141))$

${W W}_{22} = = \sqrt{11 - - {W W}_{11}^{22}} - - - - - - - - - - - - ((4242))$

图10～图13是部分结果。下面作简要说明。Figures 10 to 13 are some results. A brief description is given below.

输入信号的幅值为10mV，频率为100Hz，相位差为4°，所加噪声为正态分布噪声，服从N(0，0.4)分布，噪声幅值为0.89mV(按式(43)的定义信噪比为34)。所用的数据在两级多抽一后为4000个，即进入DSP是有48×4000个。数据窗长度为64(即L＝64，N＝32)，采用重叠窗的并行处理，每个半窗长(称为一个信号处理周期)计算一次频率、相位差、时间差。考虑这样的情况：信号的频率在第2000点时突然发生了变化，从100Hz～110Hz，这样的变化幅度对科氏质量流量计来说是比较大了，因而足以说明问题。The amplitude of the input signal is 10mV, the frequency is 100Hz, and the phase difference is 4°. The added noise is normally distributed noise, obeying the N(0, 0.4) distribution, and the noise amplitude is 0.89mV (according to the definition of formula (43) The signal-to-noise ratio is 34). The number of data used is 4000 after one extra draw in two stages, that is, there are 48×4000 data entering DSP. The length of the data window is 64 (that is, L=64, N=32), and the parallel processing of overlapping windows is adopted, and the frequency, phase difference, and time difference are calculated once for each half window length (called a signal processing cycle). Consider such a situation: the frequency of the signal changes suddenly at the 2000th point, from 100Hz to 110Hz, such a change range is relatively large for the Coriolis mass flowmeter, so it is enough to explain the problem.

$SNR SNR = = 101101 g g \frac{{A A}_{s the s}^{22}}{{σ σ}^{22}} - - - - - - - - - - - - - - ((4343))$

式中，A_s为正弦量的幅值，σ²为正态分布噪声的方差。In the formula, A _s is the magnitude of the sine quantity, and σ ² is the variance of the normal distribution noise.

图10是整个跟踪过程。图10a)表示频率估计的相对误差；图10b)表示相位差估计的相对误差；图10c)表示频率估计值。图10～13中，横坐标k表示第k个信号处理周期，而n表示第n个数据点 Figure 10 is the whole tracking process. Figure 10a) shows the relative error of the frequency estimate; Figure 10b) shows the relative error of the phase difference estimate; Figure 10c) shows the frequency estimate. In Figures 10 to 13, the abscissa k represents the kth signal processing cycle, and n represents the nth data point

图11表示频率突变前的稳态情况，即自适应算法从初始状态(任意)到收敛到突变前的频率100Hz时的稳态情况。Fig. 11 shows the steady-state situation before the frequency mutation, that is, the steady-state situation when the adaptive algorithm converges from the initial state (arbitrary) to the frequency 100 Hz before the mutation.

图12表示频率突变后的稳态情况，从图中可以看出，当信号频率从100Hz变化到110Hz时，自适应算法能够在几个周期内跟踪上信号频率的变化。下面是一些数据：Figure 12 shows the steady-state situation after a sudden frequency change. It can be seen from the figure that when the signal frequency changes from 100Hz to 110Hz, the adaptive algorithm can track the change of the signal frequency within several cycles. Here is some data:

突变前，稳态时的频率估计误差≤5.8679×10^-5 Before the sudden change, the frequency estimation error at steady state is ≤5.8679×10 ^-5

突变前，稳态时的相位差估计误差≤6.6973×10^-4 Before the sudden change, the phase difference estimation error at steady state ≤ 6.6973×10 ^-4

突变后，稳态时的频率估计误差≤1.4124×10^-5 After the sudden change, the frequency estimation error at steady state is ≤1.4124×10 ^-5

突变后，稳态时的相位差估计误差≤8.1189×10^-4 After the sudden change, the phase difference estimation error at steady state is ≤8.1189×10 ^-4

图13是本发明采用的算法中，设计参数的调整曲线。由于采用了SNR故障检测，当信号频率发生变化时，设计参数先保持在某一数值上，待算法接近收敛时，再按指数函数变化。Fig. 13 is an adjustment curve of design parameters in the algorithm adopted by the present invention. Due to the adoption of SNR fault detection, when the signal frequency changes, the design parameters are kept at a certain value first, and then changed according to the exponential function when the algorithm is close to convergence.

Claims

1. A Coriolis mass flowmeter digital signal processing system based on adaptive funnel filter and sliding Goertzel algorithm is composed of a signal conditioning circuit, an analog-to-digital converter (A/D), a Digital Signal Processor (DSP), a Liquid Crystal Display (LCD) circuit, a keyboard input circuit, a digital-to-analog converter (D/A), a power amplification circuit and software; two magnetoelectric sensors in the flowmeter convert sensed flow signals into electric signals; the electric signal output by the magnetoelectric sensor is sent to the DSP through the signal conditioning circuit and the A/D; the DSP processes the sampling data, calculates the frequency and the phase position of the vibration fundamental frequency of the flow tube, and calculates the time difference of the two paths of signals according to the phase difference of the frequency and the two paths of signals, thereby obtaining a mass flow value and a density value for LCD display; a driving signal is generated by a frequency value obtained by calculation of the DSP and is sent to the D/A, and is sent to a coil of an electromagnetic vibration exciter after power amplification, so that the flow tube vibrates; the DSP adopts a self-adaptive funnel filter to carry out frequency estimation and tracking on a sampling signal subjected to two-stage multi-pumping and one-filtering, can consider both tracking capability and tracking precision, simultaneously filters broadband noise mixed in a sinusoidal signal, and outputs an enhanced signal; and calculating the discrete Fourier coefficient of the enhanced signal by adopting a sliding Goertzel algorithm which is easier to realize at fixed points and is very suitable for the time-varying sinusoidal signal so as to solve the phase difference of the signal.

2. The coriolis mass flowmeter digital signal processing system of claim 1 characterized in that said digital signal processing tasks are performed in real time using a DSP; the DSP is used as the core control system of the system to coordinate with other parts.

3. The coriolis mass flowmeter digital signal processing system of claim 1 characterized in that an adaptive funnel filter is used to perform adaptive notch filtering on the signal after passing through a plurality of filters, calculate the fundamental frequency of the flow tube vibration, and track the change in the fundamental frequency; and subtracting the signal passing through the self-adaptive funnel type filter from the signal after the one-more-pumping and one-filtering so as to enhance the flow signal.

4. The coriolis mass flowmeter digital signal processing system of claim 1 wherein said DSP applies a band-limited random signal to said flow tube after D/a conversion and power amplification to vibrate said flow tube; under normal working condition, DSP uses the base frequency and phase position to construct new sine driving signal to drive flow tube.