**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**

_{}

_{}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:

- unify existing methods,
- give guidelines on how to design new methods and improve existing ones,
- provide a mathematical theory for stability and accuracy of these methods."

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

1) _{} is a set of subspaces of _{}

2) _{}

_{}

_{}

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)

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

- The magnitude of the noise is small compared with the signal
_{}(SNR>10dB [10]). - The approximation procedure is good enough, so the residuals in the details are expected to be zero or close to zero at most of the points.

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:

- Efficiency
- No information about the magnitude or variance of the noise is required
- Asymptotic optimality.

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**

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. [ Links ]

**[2]** MURIO D. A., The mollification method and the numerical solution of ill-posed problems, John Wiley, New York, 1993. [ Links ]

**[3]** MURIO D. A., C. E. MEJÍA, AND S. ZHAN, Discrete mollification and automatic numerical differentiation, Computers Math. Applic., 35, 1–16, 1998. [ Links ]

**[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. [ Links ]

**[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. [ Links ]

**[6]** NAGY J. G., Iterative technics for the solution of Toeplitz systems, Siam News, 28, 1–16, 1995. [ Links ]

**[7]** MATHWORKS, Wavelet toolbox: User’s guide - version 2.1 for use with Matlab, The MathWorks Inc., Natick, 2001. [ Links ]

**[8]** CHUI C., An introduction to wavelets, Academic Press, 1992. [ Links ]

**[9]** E W. AND ENGQUIST B., Multiscale modeling and computation, Notices Amer. Math. Soc., 50, 1062–1070, 2003. [ Links ]

**[10]** JANSEN M., Noise reduction by wavelet thresholding, Springer, 2001. [ Links ]

**[11]** HANSEN P. C., Rank-deficient and discrete ill-posed problems: Numerical aspects of linear inversion, SIAM , 1997. [ Links ]

**[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. [ Links ]

**[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. [ Links ] ]]>