ECG compression method using Lorentzian functions model

Digital Signal Processing 17 (2007) 319–326 www.elsevier.com/locate/dsp

ECG compression method using Lorentzian functions model Abdelaziz Ouamri a , Amine Naït-Ali b,∗ a Université des Sciences et de la Technologie d’Oran, Département Electronique, Laboratoire LSI, BP 1505, Oran El’Mnouar, 31 Oran, Algeria b Université Paris 12, LISSI EA 3956, 61 avenue du Général de Gaulle, 94010 Créteil, France

Available online 8 August 2006

Abstract An ECG compression algorithm using a combination of Lorentzian functions model is proposed in this paper. In order to estimate the parameters of the Lorentzian functions, the discrete Fourier transform (DFT) is first applied to a mean removed ECG signal from which only the most significant DFT coefficients are retained. The obtained coefficients are, then modeled as the sum of a given number of superimposed exponentially damped sinusoids (EDS), commonly identified by their amplitudes, real damping factors, frequencies and initial phases. Finally, these EDS parameters are estimated, using SVD method, then coded. The algorithm has been tested for its coding efficiency and reconstruction capability by applying it to several popular, benchmark ECG signals. Encouraging results have been obtained. © 2006 Elsevier Inc. All rights reserved. Keywords: ECG compression; SVD; Lorentzian functions; Exponentially damped sinusoids

1. Introduction It is well known that the storage and transmission of digital ECG signals are of great importance for various applications. In these applications, compressing ECG waveforms without any considerable degradation in quality is of a great consideration. The ECG compression techniques can be classified in three broad categories: direct methods, transform-based methods, and parameter extraction methods. In the direct methods category, the original samples are directly compressed [1,2]. In the transformation methods category, the original samples are transformed and the compression is performed in the new domain. Among, the algorithms which employ transform-based techniques, there are several algorithms based on discrete cosine transform [3], and wavelet transform [4–7]. In the category of the methods with extraction of parameters, one extracts some features of the signal that are used later for its reconstruction [8,9]. This paper presents an algorithm based on modeling the ECG signal as a linear combination of Lorentzian functions, to achieve effective ECG data compression. The estimation of unknown parameters of Lorentzian functions from a finite number of noisy data observations has been considered in many applications, particularly in NMR biomedical signals analysis [10]. * Corresponding author. Fax: +33 1 45 17 14 92.

E-mail address: [email protected] (A. Naït-Ali). 1051-2004/$ – see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.dsp.2006.07.003

320

A. Ouamri, A. Naït-Ali / Digital Signal Processing 17 (2007) 319–326

In this paper, we consider the estimation of unknown parameters of Lorentzian functions for the analysis of the data obtained from the ECG signals. The present paper is organized as follows. In Section 2, we briefly describe the SVD ’s algorithm [11] used in estimating the parameters of the model. In Section 3, we describe the proposed algorithm for compressing the ECG signals. Section 4 is devoted to the presentation of experimental results obtained with the MIT-BIH database. Finally, Section 5 concludes this work. 2. SVD algorithm and the Lorentzian functions model 2.1. Data model We consider that the ECG signal x(t) is modeled as a linear combination of M functions fi (t), x(t) =

M 

Ai fi (t) + b(t)

(1)

i=1

with fi (t) is a sinusoidal, power of t, Gaussian or a Lorentzian function, and b(t) is a zero mean white Gaussian with variance σ 2 , representing the instrumental noise. In this paper among the functions fi (t) we choose the Lorentzian functions: αi 1 . (2) fi (t) = π αi2 + (t − ti )2 Since the Lorentzian function has Fourier transform Fi (f ) = F T [fi (t)] = e−j 2πf ti −αi π|f | , the ECG signal has its Fourier transform X(k) as a linear combination of M exponentially damped sinusoids (EDS). Then, the resulting complex-valued X(k) may be expressed as X(k) =

M 

Am exp(j φm ) exp(αm k + j 2πkfm ),

k = 1, N/2,

(3)

m=1

where M is the order of the model, N is the number of samples, Am ∈ ∗+ are the amplitudes, αm are the real damping factors, fm ∈ [−1/2, 1/2] are the frequencies, and φm ∈ [−π, π] denote the initial phases. All these parameters must be estimated. Given N samples, the problem is then to estimate the unknown parameters {Am , fm , αm , φm , }M m=1 . Several algorithms were previously proposed for parameter estimation of exponentially damped sinusoids include those based on linear prediction, those based on systems identification such as KiSS/IQML iterative, and those based on iterative estimation of the parameters of the individual components such as AP and EM algorithms. A direct implementation of an exact maximum likelihood (ML) estimator for EDS model has a prohibitive cost of performing such a large order multidimensional search. Many of these algorithms are either computationally expensive and/or do not have good estimation performance. It is also possible to use a general ARMA model as an alternative to model the EDS signal. In this case we obtain a particular ARMA model with identical autoregression and moving average parameters, so a lot of variants of the SVD may be used. The order of the model M may be assumed known. The first aim of this paper is to verify that, with the model (1) we may achieve a good approximation of signal ECG. The estimation of the model parameters is performed by using the SVD method [11], and then used to reconstruct the ECG signal. One limits the study to the Tufts and Kumaresan (TK) method which allows more efficient estimators of the parameters than others methods like Prony or linear prediction methods. 2.2. Description of the SVD method Consider the following linear system: Aβ = −c,

(4)

A. Ouamri, A. Naït-Ali / Digital Signal Processing 17 (2007) 319–326

321

where A(i, j ) = X(i + j ) and c(i) = X(i), i = 1, N − L, j = 1, L with L such as M  L  N − M is the order of the prediction filter β. The minimum norm solution, proposed in [11] is given by β MN = −

M 

  σm−1 u+ mc vm,

(5)

m=1

where um and v m are the mth left singular vectors, and right singular vectors of the matrix A respectively, and σm is the mth singular value of A. In the case L = N we obtain the well-known Prony method. It can be easily shown that the polynomial equation [11] βMN (z) = zL + β1 zL−1 + · · · + βL = 0 has roots given by zˆ m = exp(αˆ m + j 2π fˆm ), m = 1, M, whenever M  L  N − M. The frequency and the damped factor are then:   zm ) −1 Im(ˆ ˆ αˆ m = ln |ˆzm | and fm = tan . (6) Re(ˆzm ) # # + −1 + ˆ The complexes amplitudes {hm = Am exp(j φm )}M m=1 are determined by h = V X, where V = [V V ] V is the i pseudo-inverse of the van der Monde matrix V defined by V (i, j ) = zˆ j and the vector X is defined by X(j ) with i = 1, M and j = 1, L. The matrix [V ]+ denotes the hermitian transpose of V . The amplitudes Am and the initial phases φm of the various components are computed from the complexes amplitudes:   −1 Im(hm ) ˆ ˆ . (7) Am = |hm | and φm = tan Re(hm )

3. Compression algorithm The compression algorithm consists of three main steps: Fourier transformation, parameters estimation and encoding. The ECG reconstruction is accomplished by inverting the compression operations through the use of model (3) with the estimated parameters, to reconstruct the DFT coefficients, followed by the inverse DFT to get the reconstructed ECG signal. The compression algorithm can be described as follows: • DFT and thresholding The DFT is the first stage of the ECG compression algorithm. It is applied to the preprocessed ECG signal y(n) = x(n) − mx with mx = mean[x(n)]:   Y (k) = DFTN y(n) , k = 0, N − 1,

(8)

where the index N means that N -point direct discrete Fourier transform is used. The DFT coefficients are not all significant and only the first P coefficients greater than a certain threshold level are retained, with P  N/2. In this paper, the value of P is selected by using the following ratio: RP =

EP , Et P

(9)

N/2 where EP = k=1 |Y (k)|2 and Et = k=1 |Y (k)|2 denote the desired retained energy in the thresholded coefficients and the total energy in the DFT coefficients respectively. • Parameters estimation of the model The estimation of the model parameters is achieved by the SVD method presented in Section 5.5.

322

A. Ouamri, A. Naït-Ali / Digital Signal Processing 17 (2007) 319–326

• Coding The model parameters are computed using float point arithmetic. To increase the compression ratio, the parameters should be represented as integers with predefined words lengths, denoted nA , nφ , nf , and nα for amplitude, initial phase, frequency and damping factor, respectively. In the first step of the decompression, the estimated parameters are used to reconstruct the DFT of ECG signal by Yˆ (k) =

M 

    Aˆ m exp j φˆ m exp αˆ m k + j 2πk fˆm ,

k = 1, P .

(10)

m=1

The reconstructed signal is then obtained by adding to the N -point inverse DFT of (10), the mean value mx that was used in the pre-processing step:   (11) x(n) ˆ = IDFTN (Yˆ (k)) + mx , n = 0, N − 1. 4. Performance evaluation The performance of the proposed algorithm is evaluated using the compression ratio (CR) and the percent rootmean–square difference (PRD) in % which is usually used to measure the distortion resulting from the ECG compression. The PRD is defined as  N 2 (x(n) − x(n)) ˆ PRD = n=1 × 100, (12) N 2 n=1 (x(n) − K) where x and xˆ represent, respectively, original and reconstructed signal and K is a constant. Various definitions are used in the literature. If K = 0 the PRD is called PRD1. In the case of MIT-BIH arrhythmia database, the test data included a baseline of 1024 added for storage purposes, hence K is set as 1024 and the corresponding PRD is called PRD2. If K is set as the mean mx of the original signal, the PRD is called PRD3. A low value for these measures does not guarantee a total preservation of the essential features of the original record. The compression ratio is defined as CR = Nx /Nxˆ , where Nx denote the number of bits in the original signal, Nxˆ denote the number of bits used to code the significant parameters. Taking into account the number of the parameters required by this algorithm, we have approximatively: Nx CR = , (13) M(nA + nα + n + nf ) + nN + nmx where nA , nφ , nf , nα , nmx , and nN denote the word length for the amplitude, initial phase, frequency, damping factor, the mean, and N , respectively. 5. Experimental results The MIT-BIH arrhythmia database was used to evaluate the performance of the proposed algorithm. The compression ratio and the PRD will be used for this purpose. The Lead II data set from the MIT-BIH arrhythmia database has been extensively used for a long time as a benchmark for ECG compression algorithms. Each record in the database is acquired with a sampling rate of 360 samples/s and represented with a resolution of 11 bits. For comparison purposes, we ran tests on the same data set S1 which consists of 10 mn of data in record numbers 100, 101, 102, 103, 107, 109, 111, 115, 117, 118, and 119, used by some other ECG compressors from Ref. [7]. In the following, we study the effect of certain parameters of the algorithm on its performance. 5.1. The effect of block size In this experience, we evaluate the effect of signal length N . Various values of N were considered both in the case of a blind segmentation and a supervised segmentation. In the case of a blind segmentation, we do not need any beat’s detection; on the contrary in the case of a supervised segmentation a step of beat’s detection is necessary.

A. Ouamri, A. Naït-Ali / Digital Signal Processing 17 (2007) 319–326

323

Fig. 1. Compressing records using the proposed algorithm in the case of blind segmentation. The first 1024 sample of MIT-BIH record 117, N = 1024, M = 4, CR = 16:1, PRD1 = 0.52%, PRD2 = 2.44%, PRD3 = 8.99% L = 23, P = 120, nA = nf = 11, nal = np = 10.

Fig. 2. Compressing record 117 using the proposed algorithm in the case of supervised segmentation. PRD1 = 0.47%, PRD2 = 2.37%, PRD3 = 9.36%, CR = 15.5:1, M = 5, P = 65, L = 25, nA = nf = 11, nal = np = 10.

Typical reconstructed waveforms of record 117 along with the original ones are shown in Figs. 1 and 2 in the case of blind segmentation and in the case of supervised segmentation, where the diagnostic features in the reconstructed signal are well preserved. However, the performance obtained depends on the record being compressed. We remark that the performance is better in the supervised segmentation case than when any detection is made. This tendency was confirmed by other tests on other records. 5.2. Effect of the order of the model From simulations, we observe that the PRD decreases with M leading to a better reconstruction of the ECG waveforms. However, the CR decreases, this is not surprising because all studies on ECG compression show that there is

324

A. Ouamri, A. Naït-Ali / Digital Signal Processing 17 (2007) 319–326

(a)

(b) Fig. 3. Average performance for record 117 using the proposed algorithm: CR and PRD as function of L obtained with M = 5, P = 65, and N = 512 in the case of blind and supervised segmentation, nA = nf = 11, nal = np = 10.

a compromise between CR and distortion (PRD). The best value of M may be selected as follows: M0 being the best value of M for each segment of any record can be obtained by verifying the following condition: PRD(M0 )  PRDmax , with PRDmax a fixed value of the maximum of the percentage root mean square difference. M0 can be determined by an iterative procedure as follows: • First, define the initial value of M. • Then, check if PRD(M)  PRDmax . If yes stop iterating, otherwise, increment M and repeat the above steps.

A. Ouamri, A. Naït-Ali / Digital Signal Processing 17 (2007) 319–326

325

5.3. Effect of the prediction filter order L It can be shown by simulation that for values of L near the extremes values M or P − M, the proposed algorithm does not perform well. This degradation is due to the performance limitations of the minimum norm method. Whose performances are extensively studied in Ref. [12]. This fact is illustrated by Fig. 3 where the record 117 is compressed by blind and supervised segmentation. We observe that for L lying between 12 and 50, the obtained values of PRD are small compared to those obtained for the values of L close to 5 and 55. The degradation observed near L = M and L = P − M is due to the performance limitations of the minimum-norm method. As consequence, the parameter L can be chosen between P /3 and P /2. 5.4. Effect of the segmentation It is well known that the lossy compression algorithms using the principle of segmentation suffer from the edges effect because adjacent blocks are processed independently. This may result in some discontinuities at the block boundaries. These discontinuities can be reduced by adjusting the parameters M and L in such a way that a low PRD is achieved. It can also be reduced by a simple interpolation in the zone of discontinuities. 5.5. Comparison to other coders In this section, the effects of some parameters of the proposed algorithm have been studied. Some results specific to the record 117 and 119 are shown in Table 1 and the average performance results for compressing records of the first data set using the proposed algorithm are shown in Table 2. The results show that the proposed algorithm has a performance comparable to that of existing wavelet-based or other ECG compression methods. We have summarized the advantages of the proposed method in some few points: • No codebook used for the coding/decoding purpose. Table 1 Performance results for compressing records 117 and 119 using the proposed algorithm: nA = nf = 11, nal = np = 10 Algorithm

Record

Proposed Supervised Segmentation

117 117 119

Blind segmentation

117 117 117

SPIHT [7]

117 119

1.18 5.5

Hilton [6]

117

2.6

8:1

Djohan et al. [3]

117

3.9

8:1

Rajoub [5]

117 119

N

M

400 512 512

P

L

PRD1 (%)

PRD2 (%)

PRD3 (%)

CR

5 5 4

60 50 50

27 27 27

0.39 0.41 1.18

1.94 2.06 5.29

7.68 8.13 12.09

12.68:1 15.25:1 22.46:1

5 6 6

50 50 100

27 22 27

0.64 0.63 0.46

3.21 3.11 2.31

12.6 12.51 9.16

23.16:1 21.06:1 13.56:1

1.06 1.95

8:1 21.6:1

22.19:1 23:1

Table 2 Average performance results for compressing records of the first data set using the proposed algorithm: nA = nf = 11, nal = np = 10, M = 5 First data set Blind Supervised

PRD1 (%)

PRD2 (%)

PRD3 (%)

CR

Parameters

0.67 0.35

5.43 3.0442

9.51 5.6375

11.62:1 9.21:1

N = 2048, P = 300, L = P /3, M = 40 P = 65, L = P /2, M = 5

326

A. Ouamri, A. Naït-Ali / Digital Signal Processing 17 (2007) 319–326

• No learning phase needed as it has been used in our previous work [9], and lower calculation complexity is required. • The technique includes a denoising process. • The CR can be improved when entropy coding is used on the parameters. • Better CR can be obtained if inter-beat correlation is taken into account. The complexity of the proposed method is essentially related to the SVD one. However, for low values of L, this will cause no any major problem for real time processing regarding the performances of the actual signal processing processors. 6. Conclusions In this paper, we propose to use a lossy compression algorithm for ECG data processing. The algorithm is based on modeling ECG signal as the sum of Lorentzian functions. The parameters of the model are estimated by using Minimum-Norm method. The proposed algorithm was tested on various records from MIT-BIH database. The performance obtained depends on the records and on the model parameters. For a suitable choice of the laters, low PRD can be achieved. Greater compression ratios can be achieved if the compressed signal is further encoded using simple encoding scheme such as Huffman encoder and/or run-length encoder and exploiting the strong interbeat correlation among ECG . References [1] S.M. Jalaleddine, C.G. Hutchens, R.D. Strattan, W.A. Coberly, ECG data compression techniques—A unified approach, IEEE Trans. Biomed. Eng. 37 (1990) 329–343. [2] J. Cox, F. Noelle, H. Fozzard, G. Oliver, AZTEC: A preprocessing program for real-time ECG rhythm analysis, IEEE Trans. Biomed. Eng. BME-15 (1968) 128–129. [3] Batista L., Melcher E.U.K., Carvalho L.C., Compression of ECG signals by optimized quantization of discrete cosine transform coefficients, Med. Eng. Phys. 23 (2001) 127–134. [4] S.-G. Miaou, H.-L. Yen, C.-L. Lin, Wavelet-based ECG compression using dynamic vector quantization with tree codevectors in single codebook, IEEE Trans. Biomed. Eng. 49 (2002) 671–680. [5] B.A. Rajoub, An efficient coding algorithm for the compression of ECG signals using the wavelet transform, IEEE Trans. Biomed. Eng. 49 (2002) 355–362. [6] M.L. Hilton, Wavelet and wavelet packet compression of electrocardiograms, IEEE Trans. Biomed. Eng. 44 (1997) 394–402. [7] Z. Lu, D.Y. Kim, W.A. Pearlman, Wavelet compression of ECG signals by the set partitioning in hierarchical trees algorithm, IEEE Trans. Biomed. Eng. 47 (2000) 849–856. [8] G. Nave, A. Cohen, ECG compression using long-term prediction, IEEE Trans. Biomed. Eng. 40 (1993) 877–885. [9] A. Chatterjee, A. Naït-Ali, P. Siarry, An input-delay neural network based approach for piecewise ECG signal compression, IEEE Trans. Biomed. Eng. (2005). [10] S. Umesh, D.W. Tufts, Estimation of parameters of exponentially damped sinusoids using fast maximum likelihood estimation with application to NMR spectroscopy data, IEEE Trans. Signal Process. 44 (9) (1996) 2245–2259. [11] D.W. Tufts, R. Kumaresan, Estimation of frequencies of multiple sinusoids, making linear prediction perform like maximum likelihood, Proc. IEEE 70 (1982) 975–989. [12] H. Clergeot, S. Tressens, A. Ouamri, Performances of high resolution spectral methods compared to the Cramer–Rao bounds, IEEE Trans. Astrophys. Space. Sci. P 37 (1989) 1703–1720.

Abdelaziz Ouamri was born in 1953, he received the BSc degree in electrical engineering from ENSI (CAEN), the DEA degree in automatic and signal processing from University Paris XI in 1979, and the PhD degree in signal processing from University Paris XI in 1986. He is currently a Professor at University of Sciences and Technology of Oran, Algeria. His research interests are focused on high resolution spectral methods, deconvolution, compression, and source separation. Amine Naït-Ali was born in 1972 in Oran, Algeria, he received the BSc degree in electrical engineering from University of Sciences and Technology of Oran, the DEA degree in automatic and signal processing from University Paris XI, and the PhD degree in biosignal processing from University Paris XII in 1998. He is presently an Associate Professor at the same university. His research interests are focused on physiological signal processing, optimisation, modeling, and medical signal and image compression.