Fast Time-Varying Linear Filters for Suppression of Baseline Drift in Electrocardiographic Signals

Fast time‑varying linear filters forsuppression ofbaseline drift inelectrocardiographic signals Jiří Kozumplík and Ivo Provazník* Background Heart frequency in humans can vary between around 0.67 to 3Hz (40–180 beats/min) depending on age, sex, stress, health state and a number of other factors. The lower limit of the range can be found in only a small number of physically trained persons in rest, usually in supine position. The upper limit is usually reached only in extreme physical stress. Heart frequency is usually denoted as heart rate (HR) measured by the number of contractions of the heart/min. Baseline wander is a noise with slow and usually large changes of the signal offset. Its frequency spectrum interferes with the frequency spectrum of the useful part of the signal—the ECG including its main waves and intervals: PR, ST, TP intervals, PQ segment, Abstract Background: The paper presents a method of linear time-varying filtering, with extremely low computational costs, for the suppression of baseline drift in electrocardiographic (ECG) signals. An ECG signal is not periodic as the length of its heart cycles vary. In order to optimally suppress baseline drift by the use of a linear filter, we need a high-pass filter with time-varying cut-off frequency controlled by instant heart rate. Methods: Realization of the high-pass (HP) filter is based on a narrow-band low-pass (LP) filter of which output is subtracted from the delayed input. The base of an LP filter is an extremely low computational cost Lynn’s filter with rectangular impulse response. The optimal cut-off frequency of an HP filter for baseline wander suppression is identical to an instantaneous heart rate. Instantaneous length of heart cycles (e.g. RR intervals) are interpolated between QRS complexes to smoothly control cut-off frequency of the HP filter that has been used. Results and conclusions: We proved that a 0.5 dB decrease in transfer function, at a time-varying cut-off frequency of HP filter controlled by an instant heart rate, is acceptable when related to maximum error due to filtering. Presented in the article are the algorithms that enable the realization of time-variable filters with very low computational costs. We propose fast linear HP filters for the suppression of baseline wander with time-varying cut-off frequencies controlled by instant heart rate. The filters fulfil accepted professional standards and increase the efficiency of the noise suppression. Keywords: Baseline drift, ECG signal, Time-varying linear filter Open Access © The Author(s) 2017. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. RESEARCH Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 DOI 10.1186/s12938‑017‑0316‑0 BioMedical Engineering OnLine *Correspondence: [email protected] Department of Biomedical Engineering, Brno University of Technology, Technická 12, 61200 Brno, Czech Republic Page 2 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 ST segment, and QRS complex (see Fig.1). The main goal of filtering is to suppress the noise, while the useful signal cannot be distorted more than specified in a standard recommendation. If the ECG signal is (hypothetically) periodic, its first harmonic frequency would be identical with the heart frequency. Lower frequency components would only be composed of noise. Removing these components would not distort the shape of the ECG signal. However, the ECG signal is not periodical but quasiperiodic (repetitive). Its heart frequency varies due to physiological or pathological reasons, thus it does not allow for the use of ideally set filters. Van Alsté etal. recommend attenuation of −0.5dB at heart rate. In the case of on-line processing of longer signals, they recommend −0.5dB at a fixed cut-off frequency 0.8Hz [1]. The used filter may not introduce phase distortion. Cardiac electrophysiology societies recommend the use of a linear HP filter with cut-off frequency of 0.67Hz and 3dB attenuation. The AHA reports [2] and [3] recommend an amplitude response flat within<−0.5, 0.5>dB, within the range of 1.0–30Hz. The reports recommend that low-frequency cut-off be 0.05Hz to avoid possible distortion of ST segments, but this frequency can be relaxed up to 0.67Hz (−3dB) for linear digital filters with zero phase distortion. Abacherli etal. refers in [4] to standards which recommend an HP filter without phase distortion with −3dB at 0.67Hz to suppress baseline drift during monitoring. In diagnostic devices, standards recommend attenuation of −0.9dB, at the same cut-off frequency of 0.67Hz. Luo etal. refers in [5] to the same values and recommends attenuation not more than 0.5dB at 1Hz for stress-test ECG. All mentioned recommendations and standards only deal with baseline wander suppression by linear filters with the fixed cut-off frequency. However, the main disadvantage of such filtering is that it sets a universal cut-off frequency which causes a lower efficacy in filtering ECG signals with a higher HR. It is generally known that baseline drift spectrum can significantly overlay spectrum of the useful part of ECG signals. Thus, it is desirable to use the highest possible cut-off frequency of the high-pass filter but acceptable regarding distortion of the useful part of ECG signals. This has been the reason for development of a number of alternative (non-linear) filtering methods. Meyer etal. approximated baseline drift by generating cubic splines from knots in PR intervals where we expect zero line of the ECG signal [6]. The main disadvantage of this method was the necessity of PR interval detection. The method became more efficient RR interval PR interval ST interval PQ segment QRS complex TP interval ST segment Fig. 1 Main peaks (Q, R, S), waves (T, P), time intervals (PR, ST, RR) and segments (PQ, ST) in an ECG signal Page 3 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 with increasing HRs when we obtained higher density of knots, while useful parts of the signal remained uncorrupted. Thakor etal. used a simple adaptive filter with a constant reference signal and a single weight [7]. However, this filtering method was a source of certain ST segment distortion. Jane etal. [8] described a method based on a cascade of two adaptive filters. The first, simple, adaptive filter with a constant reference input and a single weight represented a simple HP filter with cut-off frequency of about 0.3Hz. Its output fed a QRS complex detector that produced impulses derived from a rhythm of detected QRS complexes. The impulses entered the reference input of the second adaptive filter with a number of weights equal to a number of samples of the ECG cycle. The filter suppressed signals not correlated with the useful part of the ECG signal. ST segments were not distorted thanks to their direct relation to QRS complexes. A cascade adaptive filter was also used by Laguna etal. [9]. Blanco-Velasco et al. exploited methods based on empirical mode decomposition (EMD) [10]. EMD decomposed the signal on a sum of intrinsic mode functions. These were derived directly from an analysed signal and represented a simple oscillatory mode as a counterpart to the simple harmonic function used in Fourier analysis. Shusterman etal. developed a two-step procedure to correct baseline drift [11]. Firstly, two infinite impulse response filters were applied in a backward and forward direction to avoid phase distortion and obtained ECG signals free of large baseline wander. Secondly, QRS complexes were detected and the rest of the baseline drift was interpolated from determined PQ and TP intervals. Shin etal. used modified non-linear methods originally designed for the detrendization of heart rate variability signals to suppress baseline drift [12]. The resulting trend was derived from an estimation of overlapping short-time trends and was based on a smoothness prior approach. Fasano etal. applied an approach of baseline wander estimation and its removal in ECG signals based on the approximation of quadratic variation (measure of variability for discrete signals) reduction. Baseline wander was estimated by solving a constrained convex optimization problem where quadratic variation entered as a constraint [13]. Sharma etal. [14] described a method based on Hilbert vibration decomposition. The method considered the first component of the decomposition when applied to an ECG signal that corresponds to baseline wander of the signal. Zivanovic etal. introduced a baseline wander modelling using low-order polynomials [15]. Hao etal. designed in [16] filtering based on an estimation of baseline wander using the mean–median filter and discrete wavelet transform. This paper presents an application of a linear filter with a time-varying impulse response. This allows us to fulfil accepted professional standards and to increase the efficiency of the noise suppression. The main aim is to reach a maximum possible attenuation based on an instant HR. Linear filters provide the correct filtering and it is widely accepted by the biomedical engineering community. At the same time, this filter cannot be considered as optimal due to its variable heart frequency. For more effective suppression of baseline drift, an Page 4 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 HP filter with time-varying cut-off frequency related to instant heart frequency should be used. Sörnmo proposed in [17] and [18] a time-varying filter. In [17], he used a bank of low pass filters with cut-off frequencies 0.5, 0.75, 1.0, 1.25 a 1.5Hz (at −6dB), the output of the filters were subtracted from the delayed input signal. Selection of a filter from the bank was based on the length of RR interval, or estimation of drift. Sampling frequency was decimated from 500to 12.5Hz to decrease computational cost of the filtering. However, decimation and interpolation caused a higher phase delay of the filter. We propose a time-varying linear HP filter which does not introduce any phase distortion and excels with an extremely low computational load. The frequency response of the filter is adapted to an instant (interpolated) HR in each signal sample. Methods Filter design Linear phase frequency characteristics beginning at the origin of axes of the phase frequency response are a strict requirement to prevent phase distortion that could decline the ST segment. This requirement can be fulfilled by using a finite impulse response (FIR) linear filter with symmetric impulse response. The considered filters are a relatively narrow-band; thus their impulse responses are relatively long (up to hundreds samples). Direct realization of classical FIR filters leads to a high load of signal response computation which is not mainly suitable in real time applications incorporating signal processors. Low computational costs can be achieved by an elegant solution employing Lynn’s LP filters. These are called simple moving-average filters with a rectangular impulse response [19]. Realization of the required HP filter HHP is based on a narrow-band LP filter HLP of which output is subtracted from the delayed input Lynn’s LP filter is a comb filter with N zeroes uniformly positioned on the unit circle in z-plain. The first zero is at z=1. The LP filter is constructed by inserting a single pole to z=1. It results in a recursive FIR filter G with rectangular impulse response. Its transfer function is The filter may be described in its non-recursive form with the transfer function H Lynn’s LP filter as defined by (2) has a high stop-band ripple. Thus, it is recommended to use a cascade of two identical filters with transfer function GLP (see Fig.2). (1) HHP (z)=z−τ−HLP (z). (2) G (z)= z N −1 NzN−1(z−1) = 1−z −N N  1−z−1 . (3) H (z)=  1+z−1+z−2+···+z−(N−1)  /N . (4) G LP (z)=G(z)G(z)=  1−z−N N  1−z−1 2 . Page 5 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 Module of the transfer function GHP has an acceptable passband ripple from 0.0 to −0.4dB according to [2]. Module of transfer function GHP reaches 1 at fs/N, where fs is the sampling frequency. The cascade GLP can be realized in a non-recursive form with transfer function HLP. Both the recursive and non-recursive realizations of the cascade of two identical filters GLP, or HLP respectively, have a triangular impulse response. The fundamental frequency of an idealized periodic ECG signal is where NRR is a number of samples of an ECG cycle that ideally has a constant length, and TS is a sampling period. When module frequency response of an HP filter is expected to be 1 at frequency fECG, then where fs is a sampling frequency. If fS>> fECG, then (5) H LP (z)=H(z)H(z)=  1+2z−1+···+Nz−(N−1)+···+2z−2(N−1)−1+z−2(N−1)  /N2 . (6) fECG = 1 (NRR −1)TS , (7) N RR = f S fECG + 1, (8) N =round f S fECG  ≈NRR . 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 0.2 0.4 0.6 0.8 1 f [Hz] |H LP (f)| f c Fig. 2 Example of a cascade of two identical Lynn’s LP filters. The amplitude frequency response GLP (HLP) for fs = 500 Hz, N = 500, and fc = fs/N = 1 Hz Page 6 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 Thus, N can be directly derived from a number of samples of a RR interval provided that the RR interval represents the ECG cycle. A number of samples of the symmetric impulse response of the HP filter realized using a cascade of two identical LP filters and subtraction are always odd and the phase delay of the HP filter is an integer In this case, the module frequency response value will be 1 at frequency fC ≈ fECG . If we require the filter gain to be equal to −0.5dB at the frequency fC (transfer 0.9441), we need to decrease the value of N that leads to widening the stop-band of the HP filter. Considering that N corresponds to the frequency fC=fECG for zero gain decrease, the required value of NC at frequency fC for 0.5dB gain decrease is computed by multiplication or division by an appropriate constant. As we can consider the ratio of two frequencies with transfers 1 and 0.9441 (−0.5 dB) constant, we can write according to Fig.3 The constant c can be evaluated as follows. The high-pass filter HLP is derived from a low-pass filter with recursive realization described by (4). Its amplitude frequency response GLP is (9) NHP =2N−1, (10) τ HP = N HP −1 2 =N− 1. (11) c = f 1 f0 = f 2 f1 →f2=cf1 . (12)    GLPejωTs   =     1−e−jωTsN N  1−e−jωTs       2 =     e−jωTsN/2ejωTsN/2−e−jωTsN/2 Ne−jωTs/2  ejωTs/2−e−jωTs/2       2 =    sin(ωTsN/2) Nsin(ωTs/2)    2 . Fig. 3 Graphical representation of the ratio of a couple of frequencies with transfers 1 and 0.9441 (−0.5 dB). The amplitude frequency response GHP (HHP) of the derived HP filter GHP(z) = z−τ − GLP(z) for fs = 500 Hz and fc ≈ 1 Hz Page 7 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 For ω=ωc Then where 0.0559 is transfer of a low-pass filter GLP (HLP) at fc and corresponds to transfer 0.9441 of a high-pass filter HHP at fc=fECG (gain equals to −0.5dB). As fc≪fs, we can write We can easily derive that f c f0 =c= 1.253 . As the cut-off frequency and the length of the impulse response are inversely related, we can write Fixed filter realization Presented above was the idea of an optimal HP filter with its impulse response length controlled by the instant length of an ECG cycle. Such a filter has a maximum possible attenuation in a frequency band below fECG that can be reached by a linear system of this type. Further, the proposed filter is linear and it has linear phase frequency characteristics that are required for the processing of ECG signals. Recursive realization of the Lynn’s filter is not an appropriate solution. Although the single pole on a unit circle counteracts with a zero at the same position, there are rounding errors due to division by a large number N2. This negatively influences filtration. Non-recursive realization of the convolution leads to large impulse responses, thus it can be computationally expensive and slow. However, non-recursive realization can be represented by a cascade of two non-recursive (moving-average) filters with a low number of necessary operations per sample interval. The idea is based on the use of a filter H with a rectangular impulse response where we add a new input sample to a sum, then we subtract the oldest input sample and finally divide by a constant N in each sampling interval. Two such filters in a series represent an LP filter with triangular impulse response. The needed HP filter requires one more subtraction. The realized filter represents a fixed system based on Lynn’s filter with a low number of required operations. Its cut-off frequency can be chosen in advance. However, such a solution is the appropriate basis to design an elegant filter with a time-varying impulse response (and thus time-varying cut-off frequency). (13) ω TSN=2πN f c fs =2π f c f0 . (14)       sinπfc f0 Nsin  πfc f s       2 = 0.0559, (15)       sinπfc f0 πfc f0       2 ∼ = 0.0559. (16) N c= N c ≈round N 1.253 . Page 8 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 Time‑varying impulse response filter realization An ECG signal is not periodic—the length of its heart cycle(s) vary. To suppress baseline drift optimally, we need an HP filter with time-varying cut-off frequency controlled by an instant HR. The heart frequency in each time instant can only be estimated as we usually measure heart cycles from detected QRS complexes. However, the instant length of heart cycles (e.g. RR intervals) can be interpolated to obtain a signal NRR(n) to smoothly control the cut-off frequency of the HP filter being used. We use simple 1st order interpolation (by a line). Fundamental frequency of the ECG signal is then varying When the module frequency response of an HP filter is expected to be equal to 1 at frequency fECG(n), then the number of samples of the rectangular impulse response in n-th cycle is Thus, we can compute N(n) for each n directly from interpolated values of RR intervals. In other words, we design a new LP filter that always has an odd number of impulse response samples NLP(n) for each n by the above simple procedure The impulse response is triangular; its values can be easily derived. Direct realization ofan LP filter withminimum delay The designed HP filter must possess a constant phase delay despite the time-varying length of its impulse response. Therefore, the phase delay τ of the final HP filter is adapted to the maximum desirable delay that corresponds to the longest expected RR interval. The longest expected RR interval is derived from the lowest expected heart rate 40 beats/min (i.e. 0.67Hz) [2, 3]. Interpolated instant values of RR intervals are stored in a circular buffer that contains Nmax samples corresponding to the longest possible impulse response of the Lynn’s filter. The transfer function of the LP filter for current N in each n It is obvious from (17) that the LP filter impulse response has always an odd number of samples. The corresponding difference equation in non-casual form for l=n−τ is (17) fECG(n)= 1 (NRR(n)−1)TS . (18) N (n)=round f S fECG(n). (19) NLP (n)=2N(n)−1. (20) τ = N HPmax −1 2 =Nmax − 1. (21) H LP (z)=z −τ H(z)H(z) =z−(Nmax−1)  z−1+2zN−2+···+N+···+2z−(N−2)+z−(N−1)  /N2 . Page 9 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 where we used N=N(l)=N(n−τ) for simplicity of equational notation. The principle of computation of the output sample is presented in Fig.4. We should note that if N(n) varies with time, the impulse response can be gradually extended or shortened with a minimum step of two samples to keep its symmetry along the middle sample. Direct realization of the LP filter with the triangular impulse response with 2N−1 samples (see Fig.4) has no advantage of low computational complexity due to constantly changing all weights of the filter in time. Realization ofan LP filter bya cascade oftwo Lynn’s filters (knot insideQRS complexes) Using a cascade of two LP filters is more beneficial because both filters in a series have the same rectangular impulse responses (see Fig.5). A new sample is added if we consider a fixed length of the impulse response and the oldest sample is subtracted from a sum in each cycle. Under the condition that both impulse responses must be symmetrical along their middle sample (as required for integer delay of the final filter), i.e. N must be odd, the impulse response of each filter will vary with a minimum step of two samples. This results in a minimum step of four samples for two filters in a series. (22) yLP ( l )=[ x ( l + N −1)+2 x ( l + N −2)+···+ Nx ( l )+···+2 x ( l − N +2)+ x ( l − N +1)]/ N2, Fig. 4 Schematic representation of direct realization of the LP filter with minimum delay. Buffer A buffer of RR intervals (Nmax length), buffer B buffer of the input signal samples (2Nmax − 1 length), filter a filter with impulse response h(n) = {1, 2, 3,…, N,…, 3, 2, 1}, NRR number of sampling intervals, NRRmax number of samples of the longest expected RR interval, x(n) current input sample Page 16 of 16 Kozumplík and Provazník BioMed Eng OnLine (2017) 16:24 • We accept pre-submission inquiries • Our selector tool helps you to find the most relevant journal • We provide round the clock customer support • Convenient online submission • Thorough peer review • Inclusion in PubMed and all major indexing services • Maximum visibility for your research Submit your manuscript at www.biomedcentral.com/submit Submit your next manuscript to BioMed Central and we will help you at every step: 4. Abacherli R, Schmid HJ. Meet the challenge of high-pass filter and ST-segment requirement with a DC-coupled digital electrocardiogram amplifier. J Electrocardiol. 2009;42:574–9. 5. Luo S, Johnston P. A review of electrocardiogram filtering. J Electrocardiol. 2010;43:486–96. 6. Meyer CR, Keiser HN. Electrocardiogram baseline noise estimation and removal using cubic-splines and State-space computation techniques. Comput Biomed Res. 1977;10:459–70. 7. Thakor NV, Zhu Y-S. Applications of adaptive filtering to ECG analysis: noise cancellation and arrhythmia detection. IEEE T Bio-Med Eng. 1991;38(8):785–94. 8. Jane R, Laguna P, Thakor NV, Caminal P. Adaptive baseline wander removal in the ECG: comparative analysis with cubic spline technique. Comput Cardiol. 1992:143–6. 9. Laguna P, Jané R, Caminal P. Adaptive filtering of ECG baseline wander. In: 14th Annual International Conference of the IEEE EMBS, conference proceedings. Engineering in Medicine and Biology Society; 1992, p. 509–10. 10. Blanco-Velasco M, Weng BW, Barner KE. ECG signal denoising and baseline wander correction based on the empirical mode decomposition. Comput Biol Med. 2008;38:1–13. 11. Shusterman V, Shah SI, Beigel A, Anderson KP. Enhancing the precision of ECG baseline correction: selective filtering and removal of residual error. Comput Biomed Res. 2000;33:144–60. 12. Shin SW, Kim KS, Song CG, Lee JW, Kim JH, Jeung GW. Removal of baseline wandering in ECG signal by improved detrending method. Bio-Med Mater Eng. 2015;26:S1087–93. 13. Fasano A, Villani V. Baseline wander removal for biological signals by qudratic variation reduction. Signal Process. 2014;99:48–57. 14. Sharma H, Sharma KK. Baseline wander removal of ECG signals using Hilbert vibration decomposition. Electron Lett. 2015;51(6):447–9. 15. Zivanovic M, González-Izal M. Simultaneous powerline interference and baseline wander removal from ECG and EMG signals by sinusoidal modeling. Med Eng Phys. 2013;5:1431–41. 16. Hao W, Chen Y, Xin Y. ECG Baseline wander correction by mean-median filter and discrete wavelet transform. In: 33rd Annual International Conference of the IEEE EMBS, conference proceedings. Engineering in Medicine and Biology Society, Boston, MA, USA; 2011. p. 2712–15. 17. Sörnmo L. Time-varying digital filtering of ECG baseline wander. Med Biol Eng Comput. 1993;31:603–8. 18. Sörnmo L, Laguna P. Bioelectrical signal processing in cardiac and neurological applications. Cambridge: Elsevier Academic Press; 2005. 19. Lynn PA, Fuerst W. Introductory digital signal processing with computer applications. Hoboken: Wiley; 1992. 20. Willems J, Arnaud P, van Bemmel JH, Bourdillon PJ, Degani R, Denis B, Harms FMA, Macfarlane PW, Mazzocca G, Meyer J, van Eck HJR, de Medina EOR, Zywietz C. Establishment of a reference library for evaluating computer ECG measurement programs. Comput Biomed Res. 1985;18(5):439–57.