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, 116, 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, 851866, 1999.
[6] NAGY J. G., Iterative technics for the solution of Toeplitz systems, Siam News, 28, 116, 1995.
[7] MATHWORKS, Wavelet toolbox: Users 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, 10621070, 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, 8598, 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, 200210, 1997.