MULTISCALE ANALYSIS BY MEANS OF DISCRETE MOLLIFICATION FOR ECG NOISE REDUCTION

ANÁLISIS MULTIESCALA POR MOLIFICACIÓN DISCRETA PARA LA REDUCCIÓN DE RUIDO EN SEÑALES ECG

 

JUAN PULGARÍN-GIRALDO
Grupo de Investigación en Ingeniería Biomédica G-BIO, Universidad Autónoma de Occidente, Cali, jdpulgarin@uao.edu.co

CARLOS ACOSTA-MEDINA
Departamento de Matemáticas y Estadística, Universidad Nacional de Colombia, Manizales, cdacostam@unal.edu.co

GERMÁN CASTELLANOS-DOMÍNGUEZ
Grupo de Control y Procesamiento Digital de Señales, Universidad Nacional de Colombia, Manizales, gcastell9@gmail.com

 

Recibido para revisar abril 11 de 2008, aceptado diciembre 16 de 2008, versión final enero 29 de 2009

 

ABSTRACT: Multiscale analysis and computation is a rapidly evolving area of research that have had a fundamental impact on computational science and applied mathematics and have influenced the way we view the relation between mathematics and science. Even though multiscale problems have been longly studied in mathematics, such techniques suffer of the ill-posedness nature of the problem. In the solution of several ill-posed problems, discrete mollification has been used for regularization. In this paper, we propose a new technique (procedure) for multiscale analysis by using discrete mollification. The multiscale scheme is based on numerical linear algebra results combined with the mollification method applied to the Mallat algorithm. The new technique has a simple theory, an efficient implementation and compares fairly well with classical wavelet transform procedures. An application on electrocardiographic signals contaminated with typical non-white noise is considered.

KEYWORDS: ECG, GCV, multiscale analysis, mollification, non-white noise, regularization, thresholding.

RESUMEN: El análisis multiescala es un área de gran actividad investigativa con fuerte impacto en computación científica y matemática aplicada, ocupando un lugar de privilegio en la forma como se entiende la relación entre la matemática y las demás ciencias. Aunque el estudio matemático de problemas multiescala está bastante documentado, estas técnicas heredan en el entorno discreto la naturaleza mal condicionada del problema. La molificación discreta ha sido empleada con éxito en la solución numérica de diversos problemas mal condicionados. Este árticulo propone una nueva técnica de análisis multiescala, basada en molificación discreta. El procedimiento aquí propuesto usa resultados de algebra lineal numérica para implementar molificación discreta en el algoritmo de Mallat. La nueva técnica tiene una teoría simple, una implementación eficiente y proporciona resultados de calidad comparable a técnicas multiescala clásicas tipo onditas (wavelet). El proceso es aplicado en señales electrocardiográficas contaminadas con ruido no blanco, para fines de comparación.

PALABRAS CLAVE: ECG, GCV, análisis multiescala, molificación, ruido no blanco, regularización, umbralización.

 

1. MOLLIFICATION

The mollification method is a filtering procedure based on convolution, which is appropriate for the regularization of ill-posed problems and for the stabilization of explicit schemes in the solution of partial differential equations. As a regularization method for ill-posed problems, the method is well documented, for instance in [1], [2] and [3]. For the definition of mollification, the implementation of numerical boundary conditions for discrete mollification and the automatic selection of mollification parameters, we recommend [1] and [4].

1.1 Abstract Setting

Let and

(1)

We work with the following truncated Gaussian kernel:

(2)

This kernel satisfies: is zero outside and

Given locally integrable, we define its , denoted as the convolution of with the kernel

That is,

(3)

1.2 Unbounded Discrete Domain

Definition 1: Let be a discrete domain with and given real numbers and Let be a function defined by Set

(4)

with the characteristic function of Then for and a given non-negative integer, we define the of as the of with

(5)

that is,

(6)

We are particularly interested in the value of at the points in Let

(7)

Then we can write

(8)

Furthermore,

Thus,

(9)

where

(10)

Notice that satisfies

(11)

Theorem 2: Let be a function defined on with fourth derivative continuous and bounded in Let be its discrete version defined on If is another discrete function defined on and such that

(12)

then, there exists a constant such that

(13)

Furthermore, if is smooth enough, there exists a constant such that

(14)

where and are the forward, backward and central finite difference operators respectively.

For a proof, see [4].

1.3 Discrete Mollification
For working with finite length data some kind of boundary condition has to be assumed, [5], [6]. The most common and documented are Zero Padding (Dirichlet), Scaled version (Non-homogeneous Dirichlet), Even Reflection (Neumann) and Periodic. All these options lead to linear operators with Toeplitz, Scaled Toeplitz, Toeplitz plus Hankel and Circulant matrix representation, respectively.

Under periodic boundary conditions the FFT allows us to know the actual spectral action of the mollification kernel on the data, [6]. In the case of even reflection the Discrete Cosine Transform (DCT) results useful, [5].

 

2. MULTISCALE ANALYSIS (MSA)

2.1 Multiscale Decomposition Tree
The clue in the multiscale analysis is the Decomposition Tree, [7], [8] (Fig. 1).


Figure 1. Multiscale decomposition tree

The are filtered versions of at different scales or resolutions. The are the complementary details for obtaining the original source. The idea is that the approximations give us an overview of the signal and the details give us its specials features. The approximations are usually obtained from successive applications of lowpass filters, and the details from associated highpass filters.

2.2 The Idea
A decomposition tree could be set up by using discrete mollification at different resolutions as approximations and residuals as details (Fig. 2).


Figure 2. Multiscale decomposition tree using discrete mollification

2.3 Discrete MSA
Weinan E and Bjorn Engquist in [9] state:

"Given such a variety of multiscale methods in many different applications, it is natural to ask whether a general framework can be constructed.

The general framework should ideally:

Definition 3: A 3-upla of collections is called a DMSA if satisfy

1) is a set of subspaces of

2)

3)

4)

5) is an upsampling operator and is a downsampling operator.

6) if and then and

Example 4: The decomposition generated by orthogonal DWTs is a DMSA.

Definition 5: A vector is said to be decomposed to levels if it has been written in the form

(15)

where and

 

3. DMSA BY MOLLIFICATION

3.1 Kernel Adjusting

For , consider the Gaussian kernel

(16)

Its Fourier transform is

(17)

Then

(18)

So, a dyadic dilation in the kernel's parameter generates a dyadic contraction in its spectrum.

Furthermore, under periodic boundary conditions, by taking

(19)

the discrete mollification process cuts the upper half of the frequency components off.

Then the following decomposition scheme will produce a dyadic spectral decomposition of a vector in (Fig. 3).


Figure 3. Dyadic spectral decomposition

To this procedure we could associate a sub-band DMSA in the frequency domain. The problem is that, because of the shape of the Gaussian kernel's spectrum, the residuals (i.e. details) contain low-frequency information of the vector and consequently So, our scheme is not a perfect DMSA.

3.2 Mallat's Algorithm
In a perfect DMSA we can apply the Mallat's Algorithm for obtaining the MSA decomposition of a vector, [7], [8]. However, if this algorithm is directly applied to our scheme we will not obtain a perfect recovery structure.

In our case, we propose the following algorithm (Fig. 4), whose action takes place in the frequency domain and consists on applying the downsampling operator to suppress the entries associated to the upper half of the spectrum.


Figure 4. Decomposition tree (dyadic case) using discrete mollification

3.3 Non-dyadic case
For non dyadic vectors the DMSA by Discrete Mollification is implemented by increasing the values (Fig. 2).

 

4. REGULARIZATION BY THRESHOLDING

4.1 Introduction
A measured vector contaminated with additive noise, can be modeled as

(20)

where represents the exact data function and the corresponding added noise, usually a random variable. Let

(21)

Now we apply levels of DMSA to to obtain a decomposition of the form

(22)

We expect to be somewhat similar to and the details to depend in some way on the noise. At this point we do some assumptions

Under these hypotheses we make zero the detail components with magnitude under a certain threshold value, because they are suspicious of being noise. The entries with magnitude above this value are reduced in magnitude, for continuity purposes (soft-thresholding). With this procedure we obtain a modified version of the measured data which is expected to be smooth and to preserve the most important features of the signal, [7], [10].

A different threshold value is selected for each level, because each level of detail deals with different frequency bands and probably different noise behavior.

Theorem 6 (Numerical Convergence): Let such that

(23)

If and denote the results of reconstructing after levels of DMSA by mollification and soft-thresholding with threshold values applied to and respectively, then

(24)

Consequently, it holds the convergence

(25)

4.2 Threshold value selection
As usual in regularization methods, the most difficult and relevant task is the automatic selection of regularization parameters that allow an optimal balance between the errors of regularization and perturbation, [10], [11]. The procedure selected for this task is Generalized Cross Validation (GCV), which offers:

The GCV function in this case is

(26)

where denotes de number of components in with magnitude under the threshold value

 

5. EXPERIMENTAL SETUP

5.1 ECG Database
The experiments were performed on 55 registers extracted from the MIT-BIH arrhythmia database and corrupted, at a 6dB signal-to-noise ratio, with three different types of synthesized noise: electromyographic interference, 60Hz powerline interference and electrosurgical noise [12],[13].

This test shows how DMSA by discrete mollification works in noise reduction in comparison with another popular multiscale procedure: wavelet. Four levels of decomposition were used for both methods. The wavelet of choice was “db4”.

5.2 Scoring criterion
For each register, the scouring routine obtains absolute error between original and filtered signal for both methods. The result on the i-th register is denoted by.

Afterwards, mean value, variance, maximum and minimum values are computed as follows:

(27)

 

6. RESULTS

Both methods are aceptable for non-white noise reduction in ECG signals, but a discriminat factor could be the maximus error. Because of that, DMSA by mollification still shows a little pikes reduction on ECG signals.

Table 1. Absolute error in ECG signals after regularization by thresholding. Four levels, wavelet of choice: ‘db4’


Figure 5. Result of applying DMSA by Discrete Mollification for filtering ECG signals

 

7. CONCLUSION

The results show how the new procedure is capable of making an automatic filtering preserving not only the low frecuencies but also the main high frecuency features of the signal.

Additionally, this work offers a fully discrete mathematical analysis of the tool. This could be useful in the convergence analysis of algorithms using DMSA for the solution of multiscale problems.

All this imply that DMSA by mollification could be applied to a wide range of problems with strong multiscale interaction. Additional work on this matter has to be done.

 

8. ACKNOWLEDGMENT

This work has been partially supported by COLCIENCIAS, project "Centro ARTICA" and DIMA, project “Esquemas Molificados para Problemas no Lineales de Convección-Difusión” number 20201005215.

 

REFERENCES

[1] MURIO D. A. Mollification and space marching, In: Inverse Engineering Handbook (ed. K. Woodbury), CRC Press, 2002.
[2] MURIO D. A., The mollification method and the numerical solution of ill-posed problems, John Wiley, New York, 1993.
[3] MURIO D. A., C. E. MEJÍA, AND S. ZHAN, Discrete mollification and automatic numerical differentiation, Computers Math. Applic., 35, 1–16, 1998.
[4] ACOSTA C. D. AND C.E. MEJÍA, Stabilization of explicit methods for convection diffusion equations by discrete mollification, Computers Math. Applic., 55, 368 – 380, 2008.
[5] NG M. K., CHAN R. H., AND TANG W., A fast algorithm for deblurring models with Neumann boundary conditions, SIAM J. Sci. Comput., 21, 851–866, 1999.
[6] NAGY J. G., Iterative technics for the solution of Toeplitz systems, Siam News, 28, 1–16, 1995.
[7] MATHWORKS, Wavelet toolbox: User’s guide - version 2.1 for use with Matlab, The MathWorks Inc., Natick, 2001.
[8] CHUI C., An introduction to wavelets, Academic Press, 1992.
[9] E W. AND ENGQUIST B., Multiscale modeling and computation, Notices Amer. Math. Soc., 50, 1062–1070, 2003.
[10] JANSEN M., Noise reduction by wavelet thresholding, Springer, 2001.
[11] HANSEN P. C., Rank-deficient and discrete ill-posed problems: Numerical aspects of linear inversion, SIAM , 1997.
[12] FRIESEN G. M., JANNETT T. C., JADALLAH M. A., YATES S. L., QUINT S. R. AND NAGLE H. T., A comparison of the noise sentitivity of nine qrs detection algorithms, IEEE Transactions on Biomedical Engineering, 37, 85–98, 1990.
[13] OROZCO R., PÉREZ M., LORENZO J. V., GRAU R. AND RAMOS R., Evaluation of qrs morphological classifiers in the presence of noise, Computers and Biomedical Research, 30, 200–210, 1997.