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

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

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

Technical field

The present invention relates to be used for determining fundamental frequency and a kind of method and apparatus of phase differential of Coriolis (being called for short Ke Shi) mass flow sensor, particularly with digital signal processor (Digital Signal Processor, be abbreviated as DSP) be core, adopt self-adaptation infundibulate wave filter (Adaptive Funnel Shaped Filter, be abbreviated as AFF), with adaptive line enhancing signal, tracking and measuring-signal frequency; Adopt slip Goertzel algorithm (SlidingGoertzel Algorithm is abbreviated as SGA), to measure the phase differential of two paths of signals.

Background technology

Coriolis mass flowmeters (being designated hereinafter simply as Coriolis mass flowmeter) can directly be measured mass rate, is one of current flowmeter with the fastest developing speed, has broad application prospects.Coriolis mass flowmeter requires its signal processing circuit accurately to measure phase differential from two flow sensor signals, and follows the tracks of the variation of its frequency.There are some limitations in signal processing mode and system based on the coriolis mass flow sensor of analog-and digital-circuit.For this reason, people have following three kinds of schemes at present with the signal Processing that digital signal processing method and DSP are applied to coriolis mass flow sensor.

(1) based on the method for discrete Fourier transformation

The Paul Romano of U.S. Micro Motion company handles the output signal of Coriolis mass flowmeter with discrete Fourier transformation (DFT), with the DSP of TMS series processing core (" Coriolicsmass flow rate meter having a substantially increased noise immunity " as secondary instrument, US Patent No.4934196, Jun.19,1990).When non-integer-period sampled, the error of calculation of DFT can not satisfy the requirement of accuracy of instrument.For this reason, proposed bigness scale, carefully survey and the thinking of frequency-tracking.But, some gordian techniquies are not wherein disclosed.For example, when frequency change, how to gather zero crossing, etc.HeFei University of Technology is with reference to its thinking, developed adopt DFT, based on the signal processing system of the DSP of ADSP series, solved the technological difficulties that not have explanation in the United States Patent (USP), and improved aspect thin survey and the frequency-tracking.

The problem that these class methods exist is: (a) real-time is relatively poor.When frequency-tracking, constantly change sample frequency and sample and calculate, the size of power spectral value relatively again, to determine realizing integer-period sampled frequency, its time reached more than 10 seconds; (b) when containing noise in the input signal, the precision of frequency computation part is very high, and still, the precision of phase differential is affected.Will use amplitude when its reason is to calculate phase differential, and the amplitude that DFT calculates is affected by noise big.

(2) based on the method for signal amplitude

The Yoshimura Hiroyuki of Japanese fuji company amplifies the signal from two sensors in the coriolis mass flowmeters, simultaneously these two signals is sent into differential amplifier, obtains signal poor of two sensors.Traffic pilot is sent these three signal sequences into A/D converter, send into DSP again and be DFT, calculate the phase differential of two sensor signals, and therefrom select a signal as a reference, go to compensate the error (" Phase difference measuring apparatus for measuring phasedifference between input signals " that each transmission channel characteristics difference is caused, European patent application, EP0791807A2.27.08.1997; " Phase difference measuring apparatus and flowmeterthereof ", European patent application, EP 0702212A2,20.03.1996).The problem of the method is: the phase differential of two sensor signals very little (being generally less than 4 degree) so the amplitude of differential amplifier is very little, very easily is subjected to noise; Compensation method is consuming time too many, need gather 6 road signals because calculate a phase differential.

(3) method that strengthens based on adaptive line

People such as the Howard V.Derby of U.S. Micro Motion company have designed a digital information processing system based on DSP, adopt adaptive line to strengthen frequency and two phase difference between signals that (ALE) technology is determined vibrating tube, thereby measure mass rate (" Method and apparatus for adaptive lineenhancement in Coriolis mass flow meter measurement " more accurately, US Patent No.5555190, Sep.10,1996; " the adaptive line Enhancement Method and the device that are used for coriolis mass flow meter measurement ", Chinese invention patent ublic specification of application, CN on August 12nd, 1190461,1998.Because these two patents are just the same in terms of content, abbreviate " patent of Derby " below as).Two embodiment are arranged in the patent of Derby, and its signal Processing link all is made up of three parts: take out a filtering, Frequency Estimation/linearity enhancing and phase differential (mistiming) more and calculate.Take out a filtering more and adopted 8: 1 and two-stage extraction in 6: 1, phase difference calculating adopts improved Goertzel algorithm.It is different that Frequency Estimation/linearity strengthens part, the ALE of first embodiment is by simultaneously to adaptive notch filter (the Adaptive Notch Filter of two parameter estimation, be abbreviated as ANF) realize, the ANF that this embodiment adopts and ANF (the Arye Nehoral of minimum parameter, " 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) incomplete same, it has relaxed the zero right restriction of limit, datum radius is not fixed as 1, but regulate automatically by algorithm, can take into account the precision and the algorithm the convergence speed (being the speed of frequency-tracking) of frequency-tracking like this.Same, second embodiment used four adaptive notch filters, promptly at left channel and right channel two notch filters of connecting respectively.Two ANF that are separately positioned on left channel and right channel are cascades, wherein, first wave filter adopts a low reactance-resistance ratio (wide trap band) wave filter, strengthen to produce limited signal, but, can converge to rapidly on the variation range of vibrating flow tube fundamental frequency, will be transferred to the ANF of second cascade from the signal of first cascade ANF output then; The ANF of second cascade adopts a high Q value (narrow trap band) wave filter to produce the inhibition noise stronger among the technology that compares in the past or above-mentioned first embodiment and the effect of harmonic wave.

The patent of Derby also has limitation.In first embodiment, two parameters are estimated, increased the complicacy of algorithm.In a second embodiment, adopt the ANF of two cascades, the complexity of algorithm is higher than first embodiment; And, because partial ANF is based on the ANF of the first order, only afterwards just can restrain the second level in first order convergence, therefore, under the bigger situation of frequency change, speed of convergence reduces on the contrary.In addition, the Goertzel algorithm computation phase differential that adopts, might overflow (referring to J.A.Beraldin and W.Steenaart when realizing with fixed point, " Overflow analysis of a fixed-point implementationof the Goertzel algorithm; IEEE Trans.Circuits and Systems, Vol.36, No.2; February1989, pp.322-324).

Therefore, need a kind ofly can to take into account frequency-tracking precision and algorithm the convergence speed, the signal processing method of not obvious again increase algorithm complexity, and also this method is difficult for overflowing during with the fixed point realization.

Summary of the invention

The objective of the invention is to solve the shortcoming that exists in the patent of Derby, do not lose its advantage again, a kind of signal processing method and device have promptly been designed, it can take into account frequency-tracking precision and convergence of algorithm speed, the complexity of not obvious again increase algorithm, and be difficult for overflowing when realizing with fixed point.

The present invention has adopted following technical scheme in order to realize goal of the invention.The signal Processing link is made up of three parts: take out a filtering, Frequency Estimation/linearity enhancing and phase differential (mistiming) more and calculate.

It is different fully with the patent of Derby that Frequency Estimation of the present invention/linearity strengthens part, that adopt is a kind of novel wave filter one self-adaptation funnel type wave filter (Adaptive Funnel Shaped Filter, be abbreviated as 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).This wave filter is different with the secondary adaptive notch filter, and the secondary adaptive notch is characterized in having only single design parameter through being usually used in Frequency Estimation, the linear enhancing (the limit restriction factor or the title debiasing factor, in the patent of Derby, represent) with α, the selection of this parameter need take into account tracking power and tracking accuracy usually, for the ANF that has only single design parameter, the two is conflicting: (1) ρ is near 1, maximum benefit is that the variance and the deviation of Frequency Estimation is very little, i.e. the tracking accuracy height; But tracking power but descends to some extent, can not follow the tracks of bigger frequency change, and to the starting condition sensitivity; (2) ρ is away from 1, and tracking power is strong, and is insensitive to starting condition, can follow the tracks of bigger frequency change; But the variance of Frequency Estimation and deviation big (tracking accuracy is not high) during stable state.Therefore, when design ANF, can only between tracking power and tracking accuracy, get compromise.And AFF has two design parameters, zero limit has been carried out stronger restriction, order is also than the ANF height of peer, such as to estimating the situation of single parameter, ANF only needs 2 times, and AFF needs 4 times, but how many complexities of algorithm do not increase, because the parameter of required estimation yet is one.The characteristics of AFF cost function are: wide and level and smooth mouth, the steep and narrow end; Compare with ANF, wide place is wideer than it, and narrow place is narrower than it, thereby can take into account tracking power and two aspects of tracking accuracy.For big unexpected frequency change, its tracking power is stronger than ANF, and for little frequency change slowly, tracking performance is similar with it.It has the high characteristics of ANF tracking accuracy, overcomes the shortcoming that ANF can't follow the tracks of big, unexpected frequency change simultaneously.

The present invention adopts in phase difference calculating section and is easier to fix a point to realize, and become slip Goertzel algorithm (the Sliding Goertzel Algorithm of sinusoidal signal when being very suitable for, be 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 is compared a lot of benefits with traditional (comprising improved Goertzel algorithm): it can calculate fourier coefficient (1) in being less than the counting an of signal period, have the time that obtains faster; (2) when realizing, be not easy to take place numerical value and overflow with fixed-point algorithm; Become sinusoidal signal when (3) being very suitable for.

Description of drawings

Fig. 1 is a system diagram of the present invention.

Fig. 2 is the major function that DSP finishes.

Fig. 3 is the signal processing method that the present invention adopts

Fig. 4 is the frequency characteristic of the cost function of ANF.

Fig. 5 is that ANF and AFF cost function frequency characteristic compare.

Fig. 6 is the flow process of digital signal processing part.

Fig. 7 is the adaptive algorithm that Frequency Estimation/linearity strengthens.

Fig. 8 is the SGA process flow diagram.

Fig. 9 is the realization block diagram of SGA.

Figure 10 is whole tracing process diagram.

Figure 11 is the stable situation before the frequency discontinuity.

Figure 12 is the stable situation after the frequency discontinuity.

Figure 13 is the adjustment of two design parameters of adaptive algorithm.

Embodiment

Below in conjunction with drawings and Examples the present invention is elaborated.

The principle that Coriolis flowmeter is based on Coriolis force designs.When fluid flows through measuring tube, if measuring tube with a certain frequency vibration, then Zhen Dong measuring tube is equivalent to the reference frame of a uniform rotation.Because fluid and measuring tube have relative motion, so can be subjected to the effect of Coriolis force.Its reaction force acts is distorted measuring tube on the both sides of measuring tube.Magnetoelectric transducer is being adorned on the measuring tube both sides, and output signal is the voltage signal that is proportional to the measuring tube vibration velocity.When vibrating tube is when vibrating with certain frequency, its angular velocity changes by sinusoidal rule.So the signal by magnetoelectric transducer output is a sinusoidal signal, its frequency is the change frequency of angular velocity, and size is proportional to angular velocity.Measure the mass rate that the phase differential of its two paths of signals and signal frequency or mistiming just can obtain fluid.

Fig. 1 is a system diagram of the present invention.After two-way sensor signal process signal condition and the A/D conversion, enter DSP.Signal condition partly comprises low-pass filtering and amplification.Because the signal amplitude that comes out from sensor is smaller, and useful sinusoidal signal falls in many man-made noise scopes, in order to make full use of the range of A/D converter, need the signal that come out from sensor be amplified; In order to weaken the interference of man-made noise as much as possible, also adopted the low pass filter device simultaneously.The design of low-pass filtering and amplifying element is that those skilled in the art is known.The present invention selects the Butterworth second-order low-pass filter for use.

The sample frequency of A/D is fixed as 38.4KHz.Take out described rate and the filter construction selected of patent that a part adopts Derby more.Can also adopt the sample frequency of 48KHz and 12: 1 and 6: 1 two-stages to take out a filtering in addition more.Because not only one road analog input signal can adopt A/D of every road signal; Perhaps multiple signals through behind the multi channel selecting again through an A/D.

DSP is the core of signal processing system, and it carries out digital filtering, Frequency Estimation/linearity enhancing and phase differential (mistiming) to the two-way sensor signal and calculates.Selection to DSP also is not particularly limited, almost can be with any dsp chip that can buy on the market.Therefore the software that the present invention narrated adopt the fixed-point DSP chip of TI company owing to adopted patent than Derby to be more suitable for the algorithm of fixing a point to realize.

DSP also realizes the function of digital drive.

Native system comprises that also LCD shows, the keyboard input circuit.

Following description is used at typical Coriolis flowmeter and is carried out, and the fundamental frequency of vibrating flow tube approximately is 100Hz in this flowmeter.Be readily appreciated that apparatus and method of the present invention can be applied to handle the flowmeter fundamental vibration frequency of any routine.

Fig. 2 is the major function diagram that DSP finishes in the native system.DSP mainly finishes signals collecting, signal Processing and digital drive.After system powers on, DSP at first in a period of time (timing) vibrating tube is vibrated, enter normal duty then.When timing is less than, carry out vibrating tube starting of oscillation algorithm: DSP and apply a band limit random signal earlier, through D/A conversion, the after-applied vibrating tube (can apply a bit of time continuously) of giving of power amplification, vibrational system just can be vibrated (though vibration is not very desirable).Why applying random noise, is because the frequency spectrum of random noise is very wide, has a variety of frequency components; Because vibrational system has selecting frequency characteristic, and wherein, the component that the noise medium frequency is identical with the vibrating tube natural frequency will play a major role, vibrating tube begins vibration again.After timing arrives, vibration tends towards stability, be that vibrating tube is with natural frequency and stable amplitude vibration, if vibration also unstable (this can judge by the signal amplitude that compares the continuous several signal Processing computation of Period of DSP), reset timing, carry out the starting of oscillation algorithm again, till vibration is stable.At this moment DSP enters normal duty: read in sensor signal, signal is handled (denoising, ask frequency and phase place), produced drive signal (with obtaining the frequency of coming, the sinusoidal drive signals that the phase place structure is new).About timing, can be set at 1～5s, can be specifically selected according to actual conditions.

Fig. 3 is a signal processing method of the present invention.The digital signal of process conversion is transferred to by path 300 (301) selects the unit, selects the unit and is divided into two-stage.It is 8: 1 that the first order is selected, and makes sample frequency (be designated as f from 38.4KHz _S1), be reduced to 4.8KHz and (be designated as f _S2), the biography letter of the wave filter of employing is:

$G\left(z^{-1}\right)=\frac{{\left({1-z}^{-8}\right)}^5}{{\left({1-z}^{-1}\right)}^5}------\left(1\right)$

Obtain the FIR wave filter of one 36 tap behind the pole zero cancellation, this wave filter has 5 zero points at each multiple point of double sampling frequency, and this can eliminate those frequencies that are aliasing in the filter transmission band of the second level greatly; This wave filter has the smallest positive integral coefficient, can represent with a kind of accurate Computing, thereby simplify the complicacy of convolution algorithm and improve computing velocity.

It is 6: 1 that the second level is selected, and makes sample frequency be reduced to 800Hz from 4.8KHz and (is designated as f _s), the wave filter of employing is the FIR wave filter of one 131 tap of Remez transfer algorithm design.Passband be direct current to 250Hz, the stopband starting point is 400Hz; Passband is weighted to 10 ^-5, stopband is weighted to 1.

Utilize the two-stage decimation filter to have the higher anti-repeatedly performance of mixing.Signal enters self-adaptation funnel type wave filter (Adaptive Funnel filter by path 304,305 after taking out a filtering through two-stage more, be abbreviated as AFF) link, this link is finished the function that Frequency Estimation/linearity strengthens, it is estimated the frequency of input signal by adaptive algorithm, elimination simultaneously is mixed in the broadband noise in the sinusoidal signal, the output enhancing signal; Signal after the enhancing enters slip Goertzel algorithm (Sliding Goertzel Algorithm is abbreviated as SGA) link, and this link is calculated the discrete Fourier coefficient of enhancing signal, thereby obtains the phase place (poor) of signal; After phase differential, frequency are all obtained, just can obtain the mistiming of two paths of signals.

Frequency computation part/linearity is described below strengthens part.

If the two-way sensor signal can be expressed as after taking out one through two-stage more:

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)$

In the formula, y _L, y _RRepresent left and right two paths of signals respectively, A represents signal amplitude, f ₀Be the actual frequency of signal, θ ₁, θ ₂Be respectively two paths of signals phase place (unit: radian), ξ ₁, ξ ₂Be respectively two paths of signals and take out the noise that remains in after the filtering in the useful signal, f more through two-stage _sBe that the process two-stage is taken out the sample frequency after more, T is a sampling interval.

The adaptive algorithm that the present invention adopts is exactly the frequency f that will estimate signal ₀, and when signal frequency changes, want can the tracking signal frequency variation, simultaneously signal is carried out signal and strengthens and (promptly leach noise component ξ effectively ₁And ξ ₂).

The present invention realizes frequency computation part/linear enhancement function by AFF.The general type of AFF is

$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)$

Wherein $\frac{C\left(z^{-1}\right)}{D_1\left(z^{-1}\right)},i=\mathrm{1,2}$ Be respectively the transport function of two notch filters, therefore can think that AFF derives out from ANF.Because notch filter has various ways, corresponding, AFF also has various ways.The present invention adopts not to be had the AFF that the ANF (being called II type ANF) of inclined to one side form derives out from theory and (is called II type AFF, referring to document Sergio M.Savaresi, Funnel filter:a New Class of Filters for FrequencyEstimation of Harmonic Signals, Automatic, Vol.33, No.9, September, 1997, the structure that pp.1711-1718 provides).It is to be noted the AFF that adopts other form, the AFF (AFF that is called the I type) such as the ANF (being called I type ANF) from direct form derives out also can realize this function, and different forms has different relative merits.The AFF that the present invention adopts is the AFF of II type, and its transport function is by formula (6)～formula (8) decision.

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

$D_1\left(z^{-1}\right)=1+\left({1+\mathrm{\ρ}}_1^2\right){\mathrm{az}}^{-1}+{\mathrm{\ρ}}_1^2{z}^{-2}--------\left(7\right)$

$D_2\left(z^{-1}\right)=1+(1+{\mathrm{\ρ}}_2^2)a{z}^{-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)$

In the formula, a is a parameter to be estimated, Be treat the estimated signals frequency (numerical frequency, unit: radian), a with Between fixing funtcional relationship is arranged, as long as promptly estimate parameter a, just can estimate the frequency of signal.

It is the estimated value (unit: Hz) of the frequency of signal.Adaptive algorithm is adjusted coefficient a so that formula (11) minimum is a criterion.

$\left(\mathrm{\Ω}\right)=\underset{s=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}^{n-s}{\mathrm{\ϵ}(s,\mathrm{\Ω})}^2-------\left(11\right)$

Formula (11) is called cost function (Cost Function).Wherein

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

Being two limit restriction factors, is our utilizable design parameter.

Work as ρ ₁=ρ ₂The time AFF just be equal to ANF fully, from this angle, ANF is the special case of AFF.Because AFF has two design parameters, thereby just can take method more flexibly on the algorithm, not only consider tracking power but also take into account tracking accuracy than ANF.

Fig. 4 is the frequency characteristic of ANF cost function.The frequency characteristic of the cost function of notch filter is divided into 5 district: a1, near the little flat region the trap frequency; A2, the big flat region in trap frequency left side; A3, the big flat region on trap frequency right side; B1, the escarpment district in trap frequency left side; B2, the escarpment district on trap frequency right side.Have a large amount of ripples in a2, a3 district, therefore, in the zone away from trap frequency, the algorithm of seeking the cost function minimum value is easy to converge to local minimum, and this does not wish to see.Very naturally, can be with the width of trap tolerance as basin of attraction.Document (Sergio M.Savaresi, Funnelfilter:a New Class of Filters for Frequency Estimation of HarmonicSignals, Automatic, Vol.33, No.9, September, 1997, pp.1711-1718) adopt Ω ⁺-Ω ^-(Ω ⁺The slope that is cost function is+1 o'clock pairing Frequency point, Ω ^-The slope that is cost function is-1 o'clock pairing Frequency point) basin of attraction is described, according to such definition, have following conclusion (consider under the situation of identical estimation variance, ${\mathrm{\ρ}}^2=0.98,{\mathrm{\ρ}}_1^2=0.84,{\mathrm{\ρ}}_2^2=0.988):$

As seen the basin of attraction of AFF is big more a lot of than ANF, so its tracking power is also just strong a lot of than ANF, can follow the tracks of bigger frequency change.

Fig. 5 is the comparison of the frequency characteristic of AFF and ANF cost function, and from their relatively the learning of cost function, the tracking power of AFF is than ANF strong (mouthful wide than it), and tracking accuracy is with its almost (end is about the same narrow); The adaptive algorithm that it can also be seen that ANF from figure might converge to local minimum, and AFF can avoid.It is to be noted, the present invention has only adopted document (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) structure of the AFF of Ti Chuing, be suitable for the relevant formula of RML algorithm and derived voluntarily, simultaneously, determined two design parameter ρ through a large amount of emulation ₁ ², ρ ₂ ²Selection.

Fig. 6 is the flow process of digital signal processing part.

Use for reference the patent of Derby, adopt the SNR fault detect.Input signal in the adaptive algorithm is that the process two-stage is taken out the left and right road signal after more, uses symbol y _LAnd y _RExpression.

The output signal of left channel and right channel can be used as the feedback signal of Weighted adaptive unit.Though using the output signal of two channels in the Weighted adaptive unit simultaneously is possible as feedback signal, but such benefit that produces is compared with the computational complexity of increase and is not dominant, therefore only feed back with one of them channel (such as right channel), and the Weighted adaptive parameter that calculates is used in the ANF of left and right sides channel simultaneously, thereby makes two sensor delivery channels through same processing.In algorithm of the present invention, with the signal of right channel feedback signal as Weighted adaptive, determined the coefficient of wave filter by this feedback signal after, by identical wave filter (AFF), the signal after being enhanced is used x to two paths of signals again _LAnd x _RExpression.In the formula below, there is not the strict y of differentiation _LAnd y _R, but represent with y; There is not the strict x of differentiation yet _LAnd x _R, but represent with x.

Fig. 7 is the process flow diagram that Frequency Estimation/linearity strengthens the adaptive algorithm (AFF link) of part, and the main formulas of the adaptive algorithm of the AFF that the present invention adopts is:

(1) meaning of relevant vector, matrix and symbol:

Input signal: y (comprises y _L, y _R)

Priori predicated error: ε (n)

Posteriority predicated error: ε (n)

Design parameter (the limit restriction factor and the debiasing factor): ρ ₁ ², ρ ₂ ²

Forgetting factor: λ

Intermediateness variable: F (n)

Autoregression vector: Φ (n)=[φ ₁(n), φ ₂(n)] ^T

Negative gradient vector: Ψ (n)=[ ₁(n), ₂(n)] ^T

Treat the vector that estimated parameter a (t) constitutes: θ (n)=[a (n), a ²(n)] ^T

Covariance matrix: $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]$

Intermediateness matrix: M (n)

Signal after the enhancing: x (comprises x _L, x _R)

Numeral angular frequency: Ω,

Ω _kDeng

The sinusoidal signal frequency estimated value:

(2) initialization:

θ(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,{\mathrm{\ρ}}_{1\mathrm{decay}}=0.998,{\mathrm{\ρ}}_{1\mathrm{SNR}}=0.80$

${\mathrm{\ρ}}_{20}^2={0.8}^2,{\mathrm{\ρ}}_{2\∞}^2={0.995}^2,{\mathrm{\ρ}}_{2\mathrm{decay}}=0.99,{\mathrm{\ρ}}_{2\mathrm{SNR}}=0.95$

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

(3) major cycle:

$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)

The formula that the estimation of frequency is adopted is

$\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)$

The Weighted adaptive unit calculates income value (being the value of a (t)) and is used for calculating phase place and calculates phase differential, mistiming Δ t in the SGA link.Frequency computation part unit calculated rate also produces the weighted information of Goertzel wave filter.Phase calculation unit is accepted Goertzel weighted information and frequency information from the frequency computation part unit.The phase calculation unit utilization has the Fourier analysis technology of two Hanning windows and determines that 8 stream pipe cycles of phase place selection of filtering signal are as preferred length of window.Suppose that a given frequency expectation value is f _q, then preferred length of window (L=2N) is determined by following formula:

$L=2N=8\mathrm{floor}\left(\frac{f_s}{f_q}\right)------\left(28\right)$

In the formula, floor (x) expression is not more than the integer of x.

f _sFor two-stage is taken out frequency after one more, such as f _s=800, f _q=100,2N=64 then.The Hanning window is expressed as:

$h\left(n\right)=\frac{1}{2}(1-\cos \frac{2\mathrm{\πn}}{2N-1}),n=\mathrm{0,1},\·\·\·,2N-1-------\left(29\right)$

The parallel computing of the overlapping Hanning window that adopts for the patent of utilizing time domain data, the present invention to use for reference Derby as far as possible to greatest extent.

Behind the signal process Hanning window after the enhancing, adopt its discrete Fourier coefficient of Goertzel algorithm computation.Adopt each half-window after the parallel computing of overlapping Hanning window long (every N is passed through two-stage and take out time domain point after one more) to carry out the calculating of phase place, frequency and a Δ t.This moment, the calculating of frequency adopted the mean value of N estimated value to try to achieve:

${\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)$

Phase difference calculating section is described below.

The present invention adopts slip Goertzel algorithm (Joe F.Chicharo, Mehdi T.Kilani, A Sliding Goertzelalgorithm, signal processing, Vol.52, No.3, August, 1996, pp.283-297), the discrete Fourier coefficient of the signal after calculating two-way respectively and strengthen with this algorithm is obtained the phase place of two paths of signals according to the discrete Fourier coefficient, obtain the phase differential between the two paths of signals again, calculate the mistiming between the two paths of signals at last.Because the AFF link has estimated the frequency of signal and simultaneously signal has been carried out linear enhancing, adopt the algorithm of Goertzel class can directly to calculate the discrete Fourier coefficient of given Frequency point, thereby significantly reduced calculated amount (comparing) with DFT algorithm, the fft algorithm of traditional calculating discrete Fourier coefficient.

Fig. 8 is the process flow diagram of SGA, and is corresponding, provides relevant step and the formula of the SGA of the present invention's employing below:

(1) determines calculative Frequency point

(2) status condition is initialized as zero, i.e. v _k(0)=0; v _k(-i)=0, k=1,2

(3) calculate resonance filter output

$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) calculate phase place

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

(5) phase differential of two paths of signals

$\widehat{\mathrm{\φ}}={\widehat{\mathrm{\θ}}}_1-{\widehat{\mathrm{\θ}}}_2-------\left(37\right)$

(6) computing time is poor

$\mathrm{\Δt}=\frac{\widehat{\mathrm{\φ}}}{{\widehat{\mathrm{\Ω}}}_0}\·\frac{1}{f_s}------\left(38\right)$

Fig. 9 is a structural drawing of realizing the SGA algorithm.

It is pointed out that because adaptive algorithm estimates parameter alpha, thereby in fact do not need to calculate cosine value in superincumbent formula (32)～formula (34) And sine value And can adopt formula (39), formula (40) to replace respectively:

$\cos {\widehat{\mathrm{\Ω}}}_0=W_1------\left(39\right)$

$\sin {\widehat{\mathrm{\Ω}}}_0=W_2-------\left(40\right)$

Wherein,

${\mathrm{<figure-callout\; id="1"\; label="W"\; filenames="C0310894400151.png,C0310894400161.png"\; state="\{\{state\}\}">W</figure-callout>}}_1=-\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}a\left(i\right)-------\left(41\right)$

$W_2=\sqrt{1-{{\mathrm{<figure-callout\; id="1"\; label="W"\; filenames="C0310894400151.png,C0310894400161.png"\; state="\{\{state\}\}">W</figure-callout>}}_1}^2}------\left(42\right)$

Figure 10～Figure 13 is a partial results.Make brief description below.

The amplitude of input signal is 10mV, and frequency is 100Hz, and phase differential is 4 °, and institute's plus noise is the normal distribution noise, obeys N (0,0.4) and distributes, and the noise amplitude is 0.89mV (the definition signal to noise ratio (S/N ratio) by formula (43) is 34).Used data are 4000 after two-stage takes out one more, promptly entering DSP is to have 48 * 4000.Data window length is 64, and (be L=64, N=32), adopt the parallel processing of overlapping window, each half-window long (being called a signal Processing cycle) calculates a secondary frequencies, phase differential, mistiming.Consider such situation: variation has taken place suddenly in the frequency of signal in the time of the 2000th, and from 100Hz～110Hz, such amplitude of variation is bigger to Coriolis mass flowmeter, thereby is enough to say something.

$\mathrm{SNR}=101g\frac{A_s^2}{{\mathrm{\&sigma;}}^2}-------\left(43\right)$

In the formula, A _sBe the amplitude of sinusoidal quantity, σ ²Variance for the normal distribution noise.

Figure 10 is whole tracing process.Figure 10 a) represents the relative error of Frequency Estimation; Figure 10 b) relative error of expression phase difference estimation; Figure 10 c) expression frequency estimation.In Figure 10～13, horizontal ordinate k represents k signal Processing cycle, and n represents n data point

Figure 11 represents the stable situation before the frequency discontinuity, and promptly adaptive algorithm is from the stable situation of original state (arbitrarily) when converging to the frequency 100Hz that suddenlys change preceding.

Figure 12 represents the stable situation after the frequency discontinuity, as can be seen from the figure, when signal frequency when 100Hz changes to 110Hz, adaptive algorithm can be followed the tracks of the variation of signal frequency in several cycles.Be some data below:

Before the sudden change, the estimated frequency error during stable state≤5.8679 * 10 ^-5

Before the sudden change, phase difference estimation error≤6.6973 * 10 during stable state ^-4

After the sudden change, the estimated frequency error during stable state≤1.4124 * 10 ^-5

After the sudden change, phase difference estimation error≤8.1189 * 10 during stable state ^-4

Figure 13 is in the algorithm that adopts of the present invention, the adjustment curve of design parameter.Owing to adopted the SNR fault detect, when signal frequency changed, design parameter remained on earlier on a certain numerical value, when treating algorithm near convergence, pressed exponential function again and changed.

Claims (4)

1. the Ke's mass flowmeter digital signal processing system based on self-adaptation infundibulate wave filter and slip Goertzel algorithm is made up of signal conditioning circuit, A/D converter (A/D), digital signal processor (DSP), liquid crystal display (LCD) circuit, keyboard input circuit, D/A (D/A), power amplification circuit and software; Two magnetoelectric sensors in the described flowmeter convert experiencing flow signal to electric signal; The electric signal of described magnetoelectric sensor output is delivered to DSP through described signal conditioning circuit and A/D; Described DSP handles sampled data, and the frequency and the phase place at calculated flow rate tube vibration fundamental frequency place are obtained the mistiming of two paths of signals by the phase differential of frequency and two paths of signals, thereby obtain mass flow value and density values, shows for LCD; The frequency values that is calculated by DSP produces drive signal, delivers to D/A, and through power amplification, delivers to the coil of electromagnetic exciter, makes the flowtube vibration; It is characterized in that described DSP adopts self-adaptation funnel type wave filter to carry out Frequency Estimation and tracking to taking out a filtered sampled signal through two-stage more, can take into account two aspects of tracking power and tracking accuracy, elimination simultaneously is mixed in the broadband noise in the sinusoidal signal, the output enhancing signal; That employing is easier to fix a point to realize and become the discrete Fourier coefficient of the slip Goertzel algorithm computation enhancing signal of sinusoidal signal when being very suitable for, thus obtain the phase differential of signal.

2. Ke's mass flowmeter digital signal processing system according to claim 1 is characterized in that adopting DSP to finish the digital signal processing task in real time; DSP is as other part co-ordinations of system core control system.

3. Ke's mass flowmeter digital signal processing system according to claim 1, it is characterized in that adopting self-adaptation funnel type wave filter that the signal after too much taking out a wave filter is carried out adaptive trap filter, the fundamental frequency of calculated flow rate tube vibration, and the variation of tracking fundamental frequency; With taking out the signal that a filtered signal deducts process self-adaptation funnel type wave filter, flow signal is enhanced more.

4. Ke's mass flowmeter digital signal processing system according to claim 1 is characterized in that DSP applies a band limit random signal earlier, through D/A conversion, the after-applied flowtube of giving of power amplification, flowtube is vibrated; In normal operation, DSP constructs new sinusoidal drive signals with obtaining the fundamental frequency and the phase place of coming, and goes to drive flowtube.

Images