CN102940489B - Filter designing method of weak electrophysiology signal and filtering method - Google Patents

Abstract

The invention discloses a filter designing method of a weak electrophysiology signal and a filtering method. The filter designing method comprises the following steps of supposing that d is equal to (N-1)/2 according to an interpolation factor D and an prolongation period number N; initializing a filter coefficient h(n) to be a matrix with the length of (2d-1) and all the elements being 1, and carrying out upper sampling of a D point on a filter h, namely interpolating (D-1) zero points behind each coefficient h(n) of the filter h; calculating all nonzero coefficients of the filter coefficient h(n) according to the formula; and modifying and calculating a coefficient of a central point so as to design the filter h. The filtering method comprises the following steps of (1), inputting n sample point sets, namely {X<n-(N-1)D/2,?-,X<n-D>,Xn,Xn+D,?-,X<n+(N-1)D/2}, adjacent two sample point of the sample point sets spacing D sampling periods; (2) carrying out accumulation on the n sample points of the sample point sets to acquire the average; and (3) multiplying the central sample point xn by a (1-3/N) value and subtracting an average value calculated by the step (2), namely filter output y(n) at the existing moment, wherein the N is the prolongation period number. Compared with the prior art, the filter designing method of the weak electrophysiology signal and the filtering method, disclosed by the invention, have the advantages of reduced utilization of multiplying units, low calculation complexity, short filtering time and flexible parameter configuration.

Description

A kind of filter design method of faint electricity physiological signal

Technical field

The present invention relates to the filtering technique that medical signals digital processing technology neck is specifically related to remove with China baseline drift interference that faint electricity physiological signal is often subject to and Hz noise.

Background technology

Clinic system carries out long-term electrocardiogram (ECG) monitoring by computer-aided tool to high risk cardiac, requires quick, real-time faint electricity physiological signal analytical method, provides high-precision heart disease automatic diagnosis [1].Like this, pretreatment filtering is carried out to faint electricity physiological signal and just seem most important.

The most commonly Hz noise and base band interference in the multiple external disturbance that faint electricity physiological signal is subject to.Power frequency is fixing (China is 50Hz), but baseline drift changes indefinitely, and baseline drift is common external disturbance, and the breathing of people, human body are slightly mobile etc. all can make faint electricity physiological signal produce drift.The feature of baseline drift is low frequency characteristic, needs to be eliminated by the means of high-pass filtering.

The removal of baseline drift both can be realized by analog filtering, also can be realized by digital filtering.Analog filtering is traditional implementation, and its advantage is that time delay is little, cost is low, but filtering accuracy is low, very flexible; Recently, the removal of baseline drift progressively trend digital filtering realizes, and digital filtering degree of accuracy is high, motility is high, but computation complexity and cost high.Thus the cost (being multiplier and adder resource particularly) and the computation complexity that how to reduce digital filtering are the key problems eliminating electrocardio drift.

And the cut-off frequency of baseline drift filter fixes on academia really and engineering circles exists dispute and contradiction always, if high pass cut off frequency is too high, then mean that baseline drift is eliminated too thorough, this easily also filters out T ripple low-frequency component and affects diagnosis; If cut-off frequency is too low, then means that baseline drift is eliminated very slight, easily the baseline drift that should remove is not filtered out, keep down.

Nineteen ninety U.S. AHA group is pointed out, the longest RR interval is corresponding with the low-limit frequency (typically referring to the heart rate under heartbeat 40 per minute) of faint electricity physiological signal, this just means that high pass cut off frequency is the digital filter of 0.67Hz, can meet most baseline drift to eliminating [2].

There is the filtering method of a lot of scholar to the noise removing faint electricity physiological signal to be studied all the time both at home and abroad, propose a lot of filtering method.Such as foreign scholar H.Govaerts proposes a kind of method that moving average filter of application repeatedly meets the requirement of passband and stopband.Document [3] then proposes FIR filter and the wave filter [4] that combines of rolling average of the design of a kind of Kaiser window; Domestic also have a lot of relevant document [5-9].Although there has been so many wave filter to meet the demands, these methods have had respective shortcoming.Some amounts of calculation are large, are unsuitable for the realization of mini system; Some optimum configurations are dumb, as the selection of the wavelet basis function in document [6], can only choose with experience by experiment, easily affect the effect of filtering; Needing of having carries out pretreatment to faint electricity physiological signal, result in the delay of signal analysis.

The present invention proposes a kind of digital filtering method based on interpolation, and compared to conventional digital filtering method, its computation complexity is very low, and the computational resource expended considerably less (mainly multiplier resources), is suitable for portable electrocardiograph to realize.

List of references:

[1]Kligfield,P.,L.Gettes,J.Bailey,R.Childers,B.Deal,W.Hancock,G.Herpen,J.Kors,P.Macfarlane,D.Mirvis,O.Pahlm,P.Rautaharju,and?G.Wagner.Recommendations?for?the?standardization?and?interpretation?of?the?electrocardiogram,Part?I.The?electrocardiogram?and?its?technology.J.Am.Coll.Cardiol.49:1109–1127,2007.doi:10.1016/j.jacc.2007.01.024.

[2]Serafim?Tabakov,Ivo?Iliev,Vessela?Krasteva.Online?Digital?Filter?and?QRS?Detector?Applicable?in?Low?Resource?ECG?Monitoring?Systems[J].Annals?of?Biomedical?Engineering,2008,36(11):1805-1815.

[3]Faes,Th.,H.Govaerts,B.Tenvoorde,and?O.Rompelman.Frequency?synthesis?of?digital?filters?based?on?repeatedly?applied?unweighed?moving?average?operations.Med.Biol.Eng.Comp.32:698–701,1994.doi:10.1007/BF02524254.

[4]Bai,Y.,W.Chu,Ch.Chen,Y.Lee,Y.Tsai,and?Ch.Tsai.The?combination?of?Kaiser?window?and?moving?average?for?the?low-pass?filtering?of?the?remote?ECG?signals.In:Proc.17th?IEEE?Symposium?on?Computer-Based?Medical?Systems,273,2004.

[5] Zou Bo, Zhou Yuan. the digital filtering algorithm [J] of a kind of filtering baseline drift and Hz noise. Journal of Shenzhen Polytechnic, 2005,4 (4): 3-5

[6] Li An, Li Dong. the application [J] of wavelet analysis method in faint electricity physiological signal digital filtering. Computer Simulation, 2001,18 (6): 70-73

[7] Chen Jing, Qian Lifeng, to thoroughfare, Zhou Linhai. real-time signal-processing method research [J] in portable cardiac monitoring. biomedical engineering is in progress, and 2009,30 (4): 200-204

[8] Chen Tianhua, Chen Qian. a kind of digital filtering method [J] eliminating faint electricity physiological signal noise. World Science technology, 2005,7 (1): 123-126

[9] in tender, Qu Bo. the design [J] of the digital filter in faint electricity physiological signal process. information technology, 2009, the fifth phase: 61-63

Summary of the invention

Based on above-mentioned technical problem, the present invention proposes a kind of filter design method and filtering method of faint electricity physiological signal, propose the method for designing of wave filter in order to the baseline drift of removing faint electricity physiological signal disturbs.

The filter design method of a kind of faint electricity physiological signal of the present invention, comprises the following steps:

Step one, according to interpolation factor D and continuation periodicity N, makes d=(N-1)/2, d represent the group delay of wave filter;

Step 2, is initialized as all 1's matrix that length is 2d-1, wave filter h is carried out to the up-sampling of D point by filter coefficient h (n), namely insert D-1 zero point below each coefficient h (n) of wave filter h; Input n sampling point collection, namely two adjacent sample spacings D sampling periods of this sampling point collection; To the n that sampling point is concentrated individual sampling point carries out cumulative averaging, by center sampling point x _nthe value being multiplied by (1-3/N) deducts meansigma methods, the filtering being current time exports y (n), and N is continuation periodicity,

Step 3, according to above-mentioned parameter, the expression formula of filter coefficient h (n) is

N, m represent integer,

All nonzero coefficients of filter coefficient h (n) are gone out according to above-mentioned formulae discovery;

Step 4 thus design this wave filter h.

Compared with prior art, the efficient filter of the faint electricity physiological signal baseline drift of removal that the present invention proposes, if be applied to practical medical engineering field, can produce following beneficial effect:

1, with little multiplier, the filtering to electrocardiosignal can be realized.

As can be seen from the filtering flow process of the faint electricity physiological signal of Fig. 3 above, no matter periodicity N gets much, all only needs to carry out twice multiplication, this considerably reduces the number of multiplier.

2, computation complexity is low.

The calculating of amount of calculation of the present invention mainly addition.Because multiplier number is few, add calculative N-1 sub-addition, amount of calculation reduces greatly.

3, filtering time is short

Because multiplier used is few, computation complexity is low, therefore comparatively conventional filter, the calculating that filtering spends only has twice multiplying, therefore the filter filtering time that the present invention invents is very short.

4, flexible setting for parameters.

Hz noise frequency f in the present invention _p=50Hz is known constant, generally sample frequency f in practical problem _salso be known, therefore interpolation factor can according to D=f _s/ f _pcalculate.Continuation periodicity N then realizes enumerative technique according to foregoing computer Aided Design and arranges, thus arranges D and N that mate most, make parameter to arrange very flexible.

Stronger to the noise capacity of resisting disturbance of below 0.67Hz.

Example 2 is above verified, this wave filter to the noise of below 0.67Hz also can filtering clean.

Accompanying drawing explanation

Fig. 1 is the method for designing flow chart of wave filter of the present invention;

When Fig. 2 is the sampling point filtering of the present invention's difference, data chooses;

Fig. 3 is the filtering flow process of the faint electricity physiological signal of the present invention;

Fig. 4 is the relation curve of transfer function H of the present invention and periodicity N;

Fig. 5 is the output characteristic curve 5 of wave filter of the present invention: the transfer curve that the coefficient (b) that wherein (a) is wave filter is wave filter;

Fig. 6 is the coefficient of wave filter of the present invention and the partial enlargement of transmission curve;

Fig. 7 is the time domain waveform before and after the faint electricity physiological signal filtering of 900Hz;

Fig. 8 is the time domain waveform added before and after Filtering of ECG Signal of making an uproar of 900Hz;

Fig. 9 is that human body respectively leads the location drawing.

Detailed description of the invention

Below in conjunction with accompanying drawing, the specific embodiment of the present invention is described in detail.

The present invention proposes following wave filter h method for designing: assuming that sampling rate is f _sinterpolation factor is D, continuation periodicity is N, make d=(N-1)/2, then design FIR (the Finite Impulse Response that a length is 2 (d-1) D+1=(N-3) D+1, FIR) wave filter h, its coefficient h (n) is such as formula shown in (1)

Wave filter h coefficient h (n) forming process corresponding with formula (1) as shown in Figure 1.

Design of filter flow process as shown in Figure 1: first according to sample frequency f _swith Hz noise frequency f _pcalculate interpolation factor D, adopt enumerative technique that continuation periodicity N is set, this seasonal d=(N-1)/2; Filter coefficient h (n) is initialized as all 1's matrix that length is 2d-1, then h is carried out to the up-sampling (i.e. interpolation) of D point, D-1 zero point is inserted below like this each coefficient h (n) of h, simultaneously in order to reduce the length of h, the point of last for filter coefficient interpolation is deleted, then all nonzero coefficients h (n) are calculated according to the coefficient expressions of wave filter, finally revise the coefficient of computer center's point, thus design this wave filter.

Such as, make D=3, N=5, filter parameter according to Fig. 1 forms flow process, the coefficient that first can calculate this wave filter is h (n)={-1/5 0 0-1/5 0 0-1/5}, and then the coefficient of central point is determined in amendment, thus obtain h (n)={-0.2000 00 0.4000 0 0-0.2000}.As can be seen from the coefficient h (n) of wave filter, insert D-1 null value between every two nonzero coefficients, and except central point coefficient, other zero coefficient values are all equal.

Filter factor h (n) forming process as can be seen from above-mentioned design of filter flow process, a lot of null value is had in filter coefficient h (n), and these zero valued coefficients are inoperative to filtering, when namely carrying out the filtering of faint electricity physiological signal, only use nonzero coefficient.The memory data output of (D-1)/D can be saved when this means filtering.

Remove the baseline drift in the faint electricity physiological signal of human body and Hz noise, it is good by the data decimation of sampling point filtering first to press shown in Fig. 2, then carries out filtering.As shown in Figure 2, if with x _ncentered by sampling point carry out filtering, take D as data break distance, select first group of sampling point data x respectively symmetrically _n, x _n-D, x _n+D; In like manner, second group with x _n-1centered by sampling point carry out filtering, sampling point data x _n-1, x _n-D-1, x _n+D-1; 3rd group with x _n+1centered by sampling point carry out filtering sampling point data x _n+1, x _n-D+1, x _n+D+1.

Filtering flow process of the present invention comprises the following steps:

1, n sampling point collection is inputted, namely two adjacent sample spacings D sampling periods of this sampling point collection;

2, n the sampling point concentrated sampling point carries out cumulative being averaging;

3, by center sampling point x _nthe value being multiplied by (1-3/N) deducts the meansigma methods that step (ii) calculates, and is current the filtering in moment exports y (n).

As shown in Figure 3, choose 5 sampling points and illustrate the filtering flow process of faint electricity physiological signal:

As can be seen from Figure 3, the filtering of signal very simple, only need carry out n-1 sub-addition and 2 multiplication, the multiplier of consumption is little, therefore computation complexity is very low, and filtration efficiency improves greatly.

The Transmission characteristics of the wave filter one, obtained by above-mentioned method for designing

To filter coefficient h (n) carry out Fourier transformation frequency response function H (e ^{j ω}), this is the important indicator of filter transfer characteristic, namely asks for its Fourier transformation to formula (1):

For probing into the filtering characteristic at power frequency place, specific digital angular frequency might as well be studied, i.e. ω _kabove formula value during=k2 π/D, the analog frequency of its correspondence is f _k=ω _k/ (2 π) f _s=kf _s/ D, then above formula the 2nd becomes the limit item of 0/0 type, according to L' Hospital's ruler, can obtain ultimate value, namely to this differentiate simultaneously

$\begin{array}{c}\underset{\mathrm{\&omega;}\&RightArrow;k2\mathrm{\&pi;}/D}{\lim }H\left(e^{\mathrm{j\&omega;}}\right)=\underset{\mathrm{\&omega;}\&RightArrow;k2\mathrm{\&pi;}/D}{\lim }\frac{(N-2)}{N}-\frac{\cos \left[(N/2-1)\mathrm{D\&omega;}\right]\&CenterDot;(N/2-1)D}{N\cos (\mathrm{D\&omega;}/2)\&CenterDot;D/2}\\ =\frac{(N-2)}{N}-\frac{(N/2-1)D}{N\&CenterDot;D/2}=\frac{(N-2)}{N}-\frac{(N-2)}{N}=0\end{array}---\left(3\right)$

Two, filter parameter method to set up

In order to the interference making wave filter remove the baseline drift composition that the faint electricity physiological signal of below 0.67Hz is subject to, and 50Hz Hz noise is removed clean, reasonably need arrange filter parameter.After sampling rate is determined, optimum configurations here comprises the determination to interpolation factor D, wave filter continuation periodicity N.

1, the setting of interpolation factor D and setting principle thereof

The method of interpolation factor D optimum configurations is: ratio interpolation factor D being set to sampling rate and power frequency.

Its principle is: the installation warrants of interpolation factor D is that Hz noise is removed in requirement.Assuming that Hz noise frequency is f _p(such as China f _p=50Hz), then in order to remove the interference at power frequency place completely, be not difficult to know by inference by formula (3): can k=1 be made, then f _k=f ₁=f _s/ D=f _p, namely satisfy condition D=f _s/ f _ptime, wave filter can remove Hz noise completely.

For example, assuming that f _p=50Hz, f _s=300Hz, then D=f _s/ f _p=6.

2, computer Aided Design realizes enumerative technique parameters N

Before this wave filter of design, be set as removing Hz noise except calculating and except the interpolation factor D that introduces, also need to calculate continuation periodicity N (in order to prevent the phase shift between input signal and output signal, N must be odd number).There is known sample frequency f _swhen with interpolation factor D, (correspond to 3dB cut-off frequency, namely amplitude-frequency response is at 0.67Hz to make cut-off frequency corresponding frequency location), there is strict requirement to periodicity N like this.

In order to calculate the value of N more accurately, the present invention proposes a kind of enumerative technique.According to the transmission characteristic of wave filter:

$H\left(e^{\mathrm{j\&omega;}}\right)=\frac{N-2}{N}-\frac{\sin \left[(\frac{N}{2}-1)\mathrm{D\&omega;}\right]}{N\sin (\mathrm{D\&omega;}/2)}---\left(4\right)$

Thus the method to set up of periodicity N is

(1) ω=2 π f/f is made _s| _f=0.67be updated to $H\left(e^{\mathrm{j\&omega;}}\right)=\frac{N-2}{N}-\frac{\sin \left[(\frac{N}{2}-1)\mathrm{D\&omega;}\right]}{N\sin (\mathrm{D\&omega;}/2)}$ In;

(2) change the N value enumerated, computer Aided Design draws out H (e ^{j ω}) with the change curve of cycle N;

(3) from change curve, H (e is found out ^{j ω}) closest corresponding odd number N value, is the cycle parameter that will arrange.

For example, assuming that f _p=50Hz, f _s=300Hz, then D=6, according to above-mentioned three steps, can draw out H (e ^{j ω}) with enumerating the relation curve of Integer N as shown in Figure 4, for ensureing that frequency response values is 0.707, because N value is set to 59.

Three, filter transfer curve

Assuming that f _p=50Hz, sample frequency f _s=300Hz, in practical application, these two frequencies are normally fixed.Make N=59 again, the output characteristic curve that can obtain this wave filter as shown in Figure 5, can visually see from Fig. 5, and (a) is the coefficient h (n) of wave filter, the length of wave filter is (N-3) D+1, i.e. (59-3) × 6+1=337.And except the coefficients comparison of central point is large, other nonzero coefficient size is all equal, this and previous calculations coefficient are out consistent.B () figure is the transfer curve of wave filter, be accurately 0 at the frequency response H (j2 π f) at 50Hz place, namely trap effect is very good.Although this wave filter has corresponding trap at the integer multiple frequency place of 50Hz backward, but because trap scope is very little, within the frequency being limited in 2Hz, and more than 100Hz faint electricity physiological signal composition is less, therefore little on the impact of faint electricity physiological signal itself.

In order to better analyze and illustrate, now Fig. 5 is carried out partial enlargement as shown in Figure 6.From (a), not only can find out the distribution situation of nonzero coefficient size, and have 5 (i.e. D-1) null values between every two nonzero coefficients.B two red dotted lines in () figure represent that the frequency response at 0.67Hz place is 0.695, error is 0.012, and this illustrates that the wave filter of the present invention's design meets the requirement that cut-off frequency is 0.67Hz.

Five, EXPERIMENTAL EXEMPLIFICATIONThe:

Example 1: with sampling rate f _s=900Hz, primary signal is the human ecg signal that experiment gathers.Sample rate due to this design of filter is 300Hz, therefore first primary signal is carried out to the down-sampling of 3 times, and then carries out filtering with the wave filter of design.Now, each parameter of wave filter is respectively: f _s=300Hz, D=6, N=59, the waveform before and after filtering as shown in Figure 7.As can be seen from the original electro-cardiologic signals waveform of 900Hz, the baseline drift of signal is very large, and the Hz noise of 50Hz is more serious, after filtering after device, substantially eliminate baseline drift and Hz noise, obtain comparatively clean electrocardiosignal, so that the further process of signal.This illustrates that this wave filter can remove the baseline drift of electrocardiosignal really, and has well sunken removing to the Hz noise of 50Hz.

Example 2: with sampling rate f _s=900Hz, primary signal is the human ecg signal that experiment gathers.First in primary signal, adding a frequency is 0.2Hz, and amplitude is 10 ⁴noise signal, then the electrocardiosignal of making an uproar that adds obtained is carried out filtering, the time domain waveform before and after filtering is as shown in Figure 8.Can clearly find out, this wave filter can not only well remove baseline drift and Hz noise, also can well filtering to the noise signal of below 0.67Hz, thus obtains pure electrocardiosignal.

Practical application example of the present invention illustrates: with the collection of the faint electricity physiological signal in the present invention and filtering, and the processing procedure on Matlab7.9 and Visiostudio platform is example, and dominant fraction is according to collection, process, filtering three parts.

1, data acquisition

The collection of the faint electricity physiological signal of human body preserves file on a processor, conducting wire is connected human body and respectively to lead position.Position shown in Fig. 9 connects each conducting wire.

2, date processing

Matlab7.9 platform reads and processes the realization of data.

3, filtering

The each data of leading calculated above are carried out filtering according to the filtering flow process shown in Fig. 3 and can obtain the faint electricity physiological signal of pure human body, faint electricity physiological signal used in the present invention all obtains according to this embodiment.

Claims (5)

1. a filter design method for faint electricity physiological signal, is characterized in that, the method comprises the following steps:

Step one, according to interpolation factor D and continuation periodicity N, makes d=(N-1)/2, d represent the group delay of wave filter;

Step 2, is initialized as all 1's matrix that length is 2d-1, wave filter h is carried out to the up-sampling of D point by filter coefficient h (n), namely insert D-1 zero point below each coefficient h (n) of wave filter h; Input n sampling point collection, namely two adjacent sample spacings D sampling periods of this sampling point collection; Cumulative averaging, by center sampling point x is undertaken to n the sampling point that sampling point is concentrated _nthe value being multiplied by (1-3/N) deducts meansigma methods, and the filtering being current time exports y (n), and N is continuation periodicity;

Step 3, according to above-mentioned parameter, the expression formula of filter coefficient h (n) is

N, m represent integer,

All nonzero coefficients of filter coefficient h (n) are gone out according to above-mentioned formulae discovery;

Step 4 thus design this wave filter h.

2. the filter design method of faint electricity physiological signal as claimed in claim 1, is characterized in that, described step 2 also comprises the step deleted by the point of last for filter coefficient h (n) interpolation.

3. the filter design method of faint electricity physiological signal as claimed in claim 1, is characterized in that, the filtering sampling point of the method choose employing: with x _nand the x be adjacent _n-1and x _n+1centered by sampling point, take D as the sampling period of data decimation, select first group of sampling point data x respectively symmetrically _n, x _n-D, x _n+D; Second group of sampling point data x _n-1, x _n-D-1, x _n+D-1; 3rd group of sampling point data x _n+1, x _n-D+1, x _n+D+1rule choose.

4. the filter design method of faint electricity physiological signal as claimed in claim 1, is characterized in that, described interpolation factor D is set to sample frequency f _swith Hz noise frequency f _pratio D=f _s/ f _p.

5. the filter design method of faint electricity physiological signal as claimed in claim 1, is characterized in that, described continuation periodicity N adopts enumerative technique to try to achieve, and specifically comprises the following steps:

Make ω=2 π f/f _s| _f=0.67be updated to in, w represents frequency;

Change the N value enumerated, computer Aided Design draws out H (e ^{j ω}) with the change curve of cycle N;

From change curve, find out H (e ^{j ω}) closest corresponding odd number N value, is the continuation periodicity N that will arrange.