ECG Signal Denoising and Features Extraction Using Unbiased FIR Smoothing

Methods of the electrocardiography (ECG) signal features extraction are required to detect heart abnormalities and different kinds of diseases. However, different artefacts and measurement noise often hinder providing accurate features extraction. One of the standard techniques developed for ECG signals employs linear prediction. Referring to the fact that prediction is not required for ECG signal processing, smoothing can be more efficient. In this paper, we employ the p-shift unbiased finite impulse response (UFIR) filter, which becomes smooth by p < 0. We develop this filter to have an adaptive averaging horizon: optimal for slow ECG behaviours and minimal for fast excursions. It is shown that the adaptive UFIR algorithm developed in such a way provides better denoising and suboptimal features extraction in terms of the output signal-noise ratio (SNR). The algorithm is developed to detect durations and amplitudes of the P-wave, QRS-complex, and T-wave in the standard ECG signal map. Better performance of the algorithm designed is demonstrated in a comparison with the standard linear predictor, UFIR filter, and UFIR predictive filter based on real ECG data associated with normal heartbeats.


Introduction
The electrocardiography (ECG) signals play a key role in diagnosing diverse kinds of heart diseases. Because the pulses produced by heart may have subtle differences from each other and noise affects the decision accuracy, the ECG is commonly organized using precise electronic equipment [1]. Accurate measurements are especially required when data are used to extract features of ECG signals and make decisions about different kinds of heart diseases employing special software. However, even very precise measurements are typically contaminated by artefacts and noise. Artefacts may result from a variety of internal and external causes, such as the Parkinsonian muscle tremors drying electrode gel. Different kinds of noises may contaminate the ECG signal during its acquisition and transmission, such as the high frequency noise (electromyogram noise, additive white Gaussian noise, and power line interference) and low frequency noise (baseline wandering). Because noise may lead to wrong interpretation, ECG signal denoising is required. Therefore, significant attention has been paid during the last decades to develop mathematical methods and computation algorithms to extract the ECG features from regular (noisy) data with an accuracy sufficient for medical needs [2][3][4][5][6][7][8][9][10][11][12].
The Fourier transform-based approach has been developed in [13] to extract ECG signal features in the frequency domain. But, this method omits the time resolution, which affects the estimation accuracy. This issue has been circumvented in some other works by providing the time-frequency analysis without significantly affecting the resolution. In [14][15][16][17], the wavelet transform-based algorithms were developed to find applications in some medical areas. In the wavelet domain, a compromise between the frequency and time resolutions is achieved easier and one can select a proper wavelet to provide a reasonable accuracy. However, a choice of an optimal wavelet is still challenging [18] and the approach has low efficiency in smoothing ECG signals. Other algorithms tested for such needs include the principal component 2 BioMed Research International analysis (PCA) [19], linear discriminant analysis (LDA) [20], independent component analysis (ICA) [21], support vector machine [22], and neural networks [23].
One of the widely recognized approaches proposed in [24] provides noise reduction and features extraction from ECG data by employing linear prediction based on the theory developed in [25]. The approach suggests that all main features of ECG signals can be saved and gained using a onestep linear predictor. Accordingly, features extraction in the QRS complex (region of fast ECG excursions) is provided from an analysis of residual errors between the data and estimates. The approach has manifested itself as useful in the detection of arrhythmias. In other works employing onestep prediction [26,27], automatic classification of the ECG cardiac abnormalities is provided using Gaussian mixtures. Later, the prediction-based approach has been recognized as one of the standard techniques suitable for ECG signals [28].
It has to be remarked now that, from the standpoint of optimal filtering, prediction is less accurate in noise reduction than filtering and much less accurate than smoothing. On the other hand, the ECG signal processing problems do not imply predicting future values and smoothing with some time-lag may be a better choice for cardiac analysis. A classic example is the Savitzky-Golay filter (smoother) [29], which has found wide applications in diverse areas [30][31][32][33][34][35].
An optimal approach to provide smoothing and state estimation in linear models has been proposed in [36] to minimize the mean square error (MSE). A solution was found on a horizon [ − , − ] of data points, where corresponds to a fixed discrete point of the ECG signal, = − + 1, and is a discrete shift. The derived optimal FIR (OFIR) filter becomes smoothing with lag = − by < 0, provides filtering with = 0, and becomes -step predictive when > 0. However, the -shift OFIR filter requires information about noise, which is not completely available for ECG signals.
A special case of the -shift OFIR filter is the -shift unbiased FIR (UFIR) filter [36][37][38][39], which completely ignores zero mean noise and is thus more suitable for ECG signals. As being more general, the -shift UFIR filter generalizes the Savitzky-Golay filter by = −( − 1)/2 and linear predictor with > 0. Although such a filter does not require the noise statistics except for the zero mean assumptions, it provides nice near optimal estimates if is set optimally as opt by minimizing the MSE [36].
In this paper, we develop an adaptive-horizon UFIR smoothing filtering algorithm for denoising ECG signals and features extraction. We also investigate the trade-off between the UFIR smoothing filter, UFIR filter, and UFIR predictive filter and compare them to the standard linear predictor suggested in [24]. We base our investigations on the MIT-BIH Arrhythmia Database available for free use from [40,41]. Focused on the design of efficient algorithms, in this paper we limit our investigations by data associated with normal heartbeats and postpone an analysis of different kinds of heart diseases to future investigations. The rest of the paper is organized as follows. In Section 2, we describe the database, theory of the algorithms proposed, and validation. The experimental results are showed in Section 3, where we provide a comparison between the UFIR, UFIR smoothing, and UFIR predictive filtering algorithms. A discussion of the results is provided in Section 4 and generalizations with concluding remarks are given in Section 5.

Material and Methods
. . Materials. This work employs the MIT-BIH Arrhythmia Database as a benchmark. This database contains 48 ECG recordings applying two leads (e.g., MLII, V1), obtained from 47 subjects studied. The recordings have a sampled frequency 360Hz per channel with 11-bit resolution over a 10 mV range. In general, the lead most common in this database is MLII where the morphology of signal ECG is seen clearly. For testing process the MLII lead and sinusal normal rhythm are considered. Moreover, an additional test with ECG synthetic data is provided by specialized software of MATLAB designed by Karthik Raviprakash. The simulated ECG signals are based on principle of Fourier series. Here, the signal is corrupted by different levels of Gaussian white noise. Below is a brief description of the characteristics of the ECG signal.
. . . ECG Signal. The morphology of heartbeat is fundamental for extracting features of ECG signals, which are quasiperiodic as sketched in Figure 1. The heartbeat pulse can be represented with four fundamental features: P-wave (left slow excursion), QRS-complex (central fast excursion), T-wave (first right slow excursion), and U-wave (second right slow excursion).
Several problems arise while processing ECG signals shown in Figure 1: (i) Measurement data are commonly contaminated by noise, which may not be Gaussian and white.
(ii) Standard features depicted in Figure 1 must be estimated with highest accuracy to avoid medical mistakes.
(iii) The ground truth (reference model) is not available to tune an estimator optimally.
Under such conditions, two approaches relying on accurate identification of heartbeat pulses are commonly considered to extract ECG signal features: fiducial and nonfiducial. The fiducial approach refers to the characteristics such as amplitude and heart rate, which are related to the duration, amplitude, and wave shape [44][45][46][47][48][49]. The nonfiducial approach refers to quasiperiodicity of ECG signals [28] and all features are separated into three main categories based on autocorrelation, phase-space, and frequency-domain analysis.
In view of the fact that noise V in (1) may not be white Gaussian and its statistics are commonly not well-known, the best way to avoid large estimation errors is using filters that do not require information about the statistics of noise. The -shift UFIR filter, which completely ignores noise and the initial conditions, can thus be considered as a good candidate.
On a finite horizon [ − , − ] of points, the ECG signal can be represented with a degree polynomial and theshift UFIR filter [37] applied to remove noise. In accordance with [37], the UFIR estimatê| − of via data taken from [ − , − ] can be found in the convolution-based form of̂| where ℎ ( ) ≜ ℎ ( , ) is the { , }-variant impulse response of the -degree UFIR filter, the extended measurement vector S is and the filter gain matrix is given by To satisfy the unbiasedness condition where { } means an average of , then ℎ ( ) can be represented as [37,50] where ∈ [ , − 1 + ] and the coefficients ( ) are defined by [37] ( ) = (−1) Here, and ( +1)1 ( ) is the minor of D( ). Function ℎ ( , ) has the following fundamental properties [37,50]: For low-degrees, = 1 and = 2, one can find ℎ ( , ) in Appendix A. For higher degrees, ℎ ( , 0) can be computed using a recurrence relation [51,52] and then ℎ ( , ) is obtained by a projection. Of importance is that the UFIR estimate (2b) does not require the noise statistics and initial values. The zero mean noise V is allowed to have any distribution and covariance [53,54] that is a fundamental difference with optimal estimates.
. . . ECG Signal Denoising on Adaptive Horizons. The determination of optimal horizon opt is critical in UFIR filtering and smoothing [42]. Because a reference signal is unavailable for ECG data, opt can be found following [55] via the mean square value (MSV) ( , ) = { ( , ) 2 } of the measurement residual ( , ) = −̂| − ( ). It has been shown in [55] that opt ( ) can be estimated by minimizing the derivative ( , )/ aŝ opt ( ) = arg min ( , ) + 1. To optimize the horizon, let us consider a single ECG pulse shown in Figure 2. As can be seen, the ECG pulse is slowly changing, except for a fast excursion in the QRS region. The slow background requires an optimal horizon opt ≜ opt ( ) in order to provide best denoising with no essential bias. On the contrary, the QRS region requires a minimum horizon min ≜ min ( ) of + 1 points to track the behaviour exactly. The horizon must thus be adaptive.
. . . General UFIR Smoothing Algorithm. The general UFIR smoothing algorithm is represented with a pseudocode listed as Algorithm 1. It requires values of S , , and described in Section 2.1.1. Function CalculateG provides a vector G = [1 0 ⋅ ⋅ ⋅ 0] and CalculateV calculates matrix V given by (8). Vector B contains the UFIR filter coefficients (6). Provided there are V and B, the -degree matrix W is computed and estimatê| − ( ) is provided by (2b). We will use this algorithm at different horizons as smoothingUFIR function.
. . . Computing for ECG Data. Optimal horizon opt is provided by the algorithm designed, with a pseudocode listed as Algorithm 2. This algorithm requires the following variables: heartbeats data , filter degree , a set of heartbeats , the number of heartbeats , and the window width V . By defining (8) and (9) and then analysing (12), the filter coefficients specified by (13) are obtained for given , , and . Next, coefficients are computed for (11) and estimate (1) is provided as̃| − ( ). Function IntervalQRS is introduced to detect Q and S via data and a value called V , which is related to the window width.  The window covers a region including Q and S and is used as min . Because opt will produce highly biased estimates around Q and S, the window is split into three parts: where points Q int and S int determinate the window width for V . The horizon opt is applied to the first part (13) and third part (15). In the second part (14), estimation is provided with min .
Function Cat is used to concatenate estimates (13)-(15) and compute the final estimatê| − ( ). Provided we havê| − ( ), function ( ) is calculated for in the scale. This variable is saved as ( ) to represent a whole set of data ( ) for different heartbeats. Provided there are Algorithm 3: Algorithm for estimatinĝ| − in ECG signals.
various values of MSV for each , an average of ( ) is computed as avg . Because avg is accompanied with ripples causing ambiguities, it is further approximated with a cubic polynomial using function CubicFit. The derivative applied to smoothed avg while solving the optimization problem (12) yields opt .
. . . Denoising Algorithm for ECG Signals. Provided there are min and opt , the UFIR smoothing algorithm can be designed for ECG signals with a pseudocode listed as Algorithm 3. In this algorithm, function smoothingUFIR is applied to different ECG signals parts with different horizons.
Five parts of the ECG signal are recognized by function smoothingUFIR over points Q int , S int , S, and Q: In Figure 2, the first and fifth parts are defined by (16) and (20), respectively, to apply opt . The third part represents an estimate, which is equal to the original ECG signal without noise reduction on [Q, S]. The adaptive horizon apt is applied to (17) and (19). Here, is decreased from opt to with a one-time step in the QRS complex region. Beyond the QRS complex, is gradually increased from min to opt with a one-time step. Finally, function Cat provides the ECG signal estimate at the last fifth part.
. . . UFIR-Based Algorithm for Features Extraction. Provided there is denoising by Algorithm 3, in this section we develop an efficient computation algorithm for ECG signal features extraction. To this end, we first localize special points on the ECG heartbeat pulse and then compute relevant amplitudes, durations, and an angle. Unlike the approaches developed in [15,56,57], this algorithm is based on theshift and -order UFIR smoothing filter exploited with = 2 and < 0. It was found out for data used that opt = 21 suites smooth parts of the discrete ECG signal and min = 3 fits the QRS complex. Note that opt and min must be specified for each of the measured ECG signals.
Step-by-step events representing the strategy of ECG signal denoising and features extraction are shown in Figure 3. The original discrete-time ECG signal (a) is smoothed as (b) using Algorithm 3. Then the ECG signal features are extracted as in the following: is estimated asR and a window is introduced with two points, Q and S . The estimateQ of Q is found as the least in the interval between Q andR. The estimateŜ of S is found as the least betweenR and S .
(ii) Figure 3(d): provided there areQ,R, andŜ, the QRS complex is suppressed to save only P and T waves. Then the estimatesP of P andT of T are obtained similarly by suppressing one of the waves.
(iii) Figure 3(e): provided there isP, the P wave is split into two segments, P 1 and P 2 , where P 1 is extended from the initial point to P 2 . In segment P 1 , we apply the derivative. Next, we consider a small section of the resulting signal and find a global maximum. We consider it as a start point of P wave and call it P onset . In segment P 2 , we also apply the derivative, consider a small portion of the resulting signal, and find a global minimum. This minimum, which corresponds to the end of P wave, is called P offset . Values of P onset and P offset are located at points (Although P on p and P off p are omitted in Figure 3, their values represent the temporal line in the ECG signal. These points can be used to compute features of the duration and applied to R p , Q p S p , P p , T p , T on p , T off p , and S * p , which are described in Algorithm 4.) P on p and P off p , respectively. Then, the duration of P wave is computed as P dur =

BioMed Research International
Original ECG signal where a and b are vectors created fromŜ and S * . These values are localized in S p and S * p . Vectors a and b have two components dependent onŜ and S * . We consider a flat part, where S p andŜ represent the origin zero point. We sum a temporal unity from the origin, obtain S * and S * p , and rename S * as S * y and S * p as S * x from xy plane. We then compute a = S * x +S * y and b = 0 + S * y and estimate via (18).

. . . Algorithm Design for Features Extraction of ECG Signals.
A pseudocode of the algorithm designed to extract features of ECG signal is shown as Algorithm 4. Here, is the smoothed ECG signal represented asx in Figure 3; b is the number of heartbeats; is a variable, which represents the reference line; is the data sample frequency; V is a value, which determines the window width to cover Q and S points ( Figure 3). The algorithm output consists of estimates of the ECG signal features such asP of P, P amp of the P amplitude, P dur of the P duration, QRS of the QRS amplitude, QRS dur of the QRS duration,T of T, T amp of the T amplitude, T dur of the T duration, and̂of the ST angle . All these features are extracted from the smoothed signal .
The algorithm starts by computingR as the ECG signal maximum, using function max. Function IntervalQRS is applied to compute points Q and S . The V variable determines the window width to cover the QRS complex and obtainQ andŜ as two minima between points Q and S . Function min is used to find the above-mentioned points. The supress function is used to suppress the QRS complex. Function max is used to estimate P and T. Function diff is introduced to compute the derivatives in the P 1 , P 2 , T 1 , and T 2 intervals. Functions max and min with function diff are used to find P onset , P on p , P offset , P off p , T onset T on p , and T offset T off p . Provided the above-mentioned values are considered, the duration is estimated of P and T features. Function length is introduced to compute the signal length. The variable determines the reference line for computing the amplitude features. This variable is equal to P offset . Function vector is used to provide vectors a and b based on S p ,Ŝ, S * p , and S * . Finally, function arcos is used to compute an angle between vectors a and b. Note that all the above introduced functions are available from the authors by request. . . Validation. Several methods have been proposed during decades for ECG signal features extraction. Among these methods, the Linear Predict approach proposed in [25] and developed by Martis [28] has been recognized as one of most efficient. The method employs the following model: in which ( ) is the original ECG pulse, is the estimator order, and ( ) is the linear prediction coefficient. The estimatê( ) is provided as a linear weighted combination of ( − ), = [1, ]. The residual error is considered as the ECG signal fraction, which cannot be predicted. To compare with the UFIR filter, we will assign = 2 as suggested by Lin et al. [24]. The UFIR filter predicts estimates with > 0 and both the prediction estimator (17) and the UFIR predictive filter (1) employ discrete linear prediction of the undergoing process via its noisy data. Even so, there are some zones in the ECG picture where linear predictors are unsuccessful in extracting ECG features. Therefore, a comparative analysis of different methods developed in [36][37][38][39] is required.
The real ECG data has unknown model and noise. A suitable metric is the concentration of error which is the difference between the estimate and measurement for different parts of ECG signals. The box plot allows giving indices related to the error dispersion and concentration. Moreover, a critical measure of denoising efficiency in any estimator is the MSE at its output. We provide the relevant study based on synthetic ECG signals generated using MATLAB. The ECG signal is contaminated by zero mean additive white Gaussian noise (AWGN) providing different SNR values.
The assessing performance for the features extraction is to analyse the concentration of the features seeing the effect of noise in the estimated features.
As can be seen, behaves similarly for different degrees . It can also be observed that opt generally grows with and elevates to opt = 27 when = 4. Particularly in Algorithm 3, an analysis of estimation errors produced by the 2-degree and 3-degree UFIR filters reveals no significant differences, except for the horizon length, which inherently grows with . This is explained by the fact that = −( − 1)/2 makes the noise power gain (NPG) of both filters equal [37]. The role of on the smoothing filter NPG has been studied by Shmaliy et al. in [37]. However, choosing = 2 reduces the computational complexity, while saving the estimation accuracy, and we accept = 2 as near optimal. Effect of on the estimation accuracy is illustrated in Figure 5.
. . Critical Evaluation of Denoising Algorithms. In Figure 6(a), we illustrate typical denoising errors produced by the predictive filter, filter, and smoothing filter, all having batch structures. A part of the ECG signal taken from [120:200] is zoomed in Figure 6(b). The denoising errors are sketched in Figure 7.
As can be seen, all UFIR filters are successful in denoising with consistent errors. Even so, the UFIR smoothing filter does it more precisely while the predictive filter produces more errors. The medians of errors produced by the algorithms and represented with the dispersion are listed in  . This figure suggests that the UFIR smoothing filter outperforms both the UFIR filter and the standard linear predictor developed in [28] for ECG signals. An analysis of the signal-to-noise rations (SNRs) at the filters outputs will be provided next.
. . Effect of SNR on the Estimator MSE. The root MSEs (RMSEs) are shown in Figure 8 as functions of the SNR depicted in decibels (dB) at 18 discrete points with a step of 5dB. It follows that the UFIR smoothing filter outperforms other solutions in a wide range of SNR values. For 0 ⩽ SNR < 15 dB, higher accuracy is achieved with a constant and, for SNR > 15 dB, with an adaptive .
. . Applications to ECG Signals. Based upon the above developed UFIR-based approach, we now apply Algorithm 3 to the ECG signal database and extract special features depicted in Figure 1. The results obtained using the designed UFIR smoothing Algorithm 2 (UfirSmooth), UFIR predictive algorithm (Predictor UFIR), and basic linear predictor (Linear Predict) [25] are sketched in Figures 9 and 10. In these figures, 100 synthetic heartbeats are processed at each time index. This synthetic ECG signal is contaminated by AWGN at 35 dB with properties similar to the original data.
In Figure 11, we show dispersions and concentrations of the estimated features about their means. Shadowed areas represent features extracted by smoothing and it follows that the outputs of the filter and linear predictor are more vulnerable. Furthermore, noise dominates in the predictive filters outputs. This experiment was based on healthy records of MIT-Arrhythmia database (lead MLII) analysing 1000 heartbeats. Overall, the UFIR smoothing approach developed in this paper always produced better estimates than other linear methods considered.

Discussion
The purpose of this study is denoising the attached noise in ECG signals using a UFIR smoothing filter for features extraction. This work is focused on the morphological features extraction individual ECG signal processing with normal rhythm. A principal finding in applying the proposed method is the considerable reduction of noise with an optimum and adaptive horizon for real ECG data. This reduction contributes determining with better precision the features associated with the heartbeat.
From analysis of errors variability in real ECG signals and SNRs based on ECG synthetic data in different estimators has shown that the UFIR smoothing filter with adaptive horizon outperforms the linear predictor [25][26][27][28] and other UFIR solutions such as the UFIR filter and UFIR predictive filter on MIT-BIH arrhythmia dataset. Let us notice again that the approaches based on linear prediction were recognized   as standard for the ECG signal features extraction [28]. In this regard, better performance of the smoothing algorithm developed in this paper opens new horizons in achieving higher accuracy and reliability in detecting different kinds of heart diseases.
The UFIR smoothing filter performance was optimized by making the averaging horizon adaptive. Note that such an opportunity has not been used in the design of known linear predictors for ECG data. As a result, we have achieved the following improvements: (1) Suboptimal denoising of ECG signals with no requirements to noise, except for the zero mean assumption.
(2) Unbiased filtering in the QRS region, in which the ECG signal demonstrates rapid excursions.
Such abilities of the UFIR smoothing filter have resulted in higher estimation accuracy, namely, in smaller variability of the estimated features around their mean values. In this regard, let us notice that larger variability in the standard linear predictor is due to larger errors and instability caused by unknown future data and errors in the determination of the predictor coefficients determined by the correlation method. Accordingly, errors in the determination of the prediction function lead to larger prediction errors (random and regular). This has appeared to be particularly true for the P amp and T amp values, which are estimated by other methods with much larger errors. Estimates of QRS e and QRS dur by different methods have appeared to be consistent, because these values are not affected by noise as much as other features. Nevertheless, the UFIR smoothing filter has demonstrated smaller errors even for QRS amp . In the cases of both T dur and angle , one watches for highly unstable estimates provided by the prediction-based filters, while the proposed UFIR smoothing filter has produced acceptable estimates. Also, it is important to clarify that the evaluation of features is analysed from the consistence of data near the average of the measurement. This is shown analysing the number of outliers. However, in this scenario, the quality features are not strictly analysed because the ECG signal used is just under normal conditions.

Conclusions
The UFIR smoothing filtering approach developed in this paper for ECG signals denoising and features extraction has demonstrated an ability to outperform the linear predictorbased one [25], which is recognized as one of the standard techniques for ECG signals. That has become possible by optimizing the order and averaging horizon for the UFIR smoothing filter in a way such that the horizon has become adaptive to different parts of ECG signals. A comparison of the UFIR predictive, filtering, and smoothing estimates has revealed a considerable difference in denoising in favor of the smoothing one. The results have also indicated that features extracted using the smoothing filter are more reliable and less prone to large deviations from average values. This is definitely an important advantage for medical needs. As a future work, we consider extracting features of ECG signals in discrete-time state-space by developing the fast iterative UFIR smoothing filtering algorithm and optimize it for different orders and kinds of heart diseases.