Wavelet transform in electrocardiography—data compression

lntam~aMl

Jaumal Cal

Medical lnformatics ELSEVIER

International

Journal

of Medical

Informatics

45 (1997)

11 l- 128

Wavelet transform in electrocardiography-data I. Provaznik Department

of‘ Biomedical

Engineering,

compression

*, J. Kozumplik

Faculty of’ Elecrrical Engineering and Computer Purkynova 9la, 612 00 Bmo, Czech Republic

Science,

Technical

University

Brno,

Abstract An application of the wavelet transform to electrocardiography is described in the paper. The transform is used as a first stage of a lossy compression algorithm for efficient coding of rest ECG signals. The proposed technique is based on the decomposition of the ECG signal into a set of basic functions covering the time-frequency domain. Thus, non-stationary character of ECG data is considered. Some of the time-frequency signal components are removed because of their low influence to signal characteristics. Resulting components are efficiently coded by quantization, composition into a sequence of coefficients and compression by a run-length coder and a entropic Huffman coder. The proposed wavelet-based compression algorithm can compress data to average code length about 1 bit/sample. The algorithm can be also implemented to a real-time processing system when wavelet transform is computed by fast linear filters described in the paper. 0 1997 Elsevier Science Ireland Ltd. Keywords:

ECG signals; Discrete wavelet transform;

Data compression;

1. Introduction

Compression algorithms are proposed for real-time data communication and/or for low-cost data storing. Electrocardiographical (ECG) signals are usually sampled at frequency 500 Hz and digitized to 12 bits in l-8 channels that represents data flow of 600048 000 bps. To transfer an one-lead ECG *Corresponding author. Tel.: + 42 5 759310; fax: f42 5 746757; e-mail: [email protected] alternatively [email protected] web: http://www.fee.vutbr/cz 1386.5056:97/$17.00

Q 1997 Elsevier

PII SOO20-7101(97)00040-8

Science

Ireland

Ltd.

All rights

Huffman

coding

signal through common communication telephone line of 4800 bps, the data must be reduced. Also, memory requirements for storing one 10-s 12-lead ECG record digitized as described above are 90 kb. The ECG signals digitized as described above are redundant. The redundancy follows from the fact that distribution of signal values is not uniform. The optimum method to suppress redundancy is entropic coding if the samples are statistically independent. The most familiar method is Huffman code that is used in many versions in almost all published reserved

112

I. Prol;aznik,

J. Kozumplik

/International

x(n)

Journal

of Medical

Informatics

45 (1997)

111~ 128

x’(d)

YIO

.

Hh -‘12

‘12

---)

Ydd

Hd ‘12

Hh -‘12

...

Fig. 1. Block diagram of dyadic wavelet transform (left) and its associated inverse transform (right). high-pass filters, Hd, Fd are low-pass filters, 12 denotes decimation by a factor 2, T2 denotes interpolation 2.

compression systems. Nevertheless, entropic coding is not sufficient compression scheme because of high entropy of ECG signals. The entropy can be a posteriori decreased in several ways. A deterministic part of signals can be removed by predictive, interpolate and differentiate estimation which are optimum in minimization of residual error signal entropy [l]. A more successful method to decrease the entropy of ECG signal is using extensive Markov process which considers both the coding redundancy and the intersample redundancy [2]. The average code length for the extensive Markov system on the 2”d difference signals is about 2.5 bit/s while average Huffman code length for 2”d difference signals is about 3.3 bit/s [2]. The described entropic coding of residuals can be without loss, i.e. the compressed signal can be fully recovered. An interesting approach is subband coding that was inspired by speech processing. The first articles describe a division into six octave bands [3], or into four bands of the same bandwidth [4]. The sampling frequency of signal of particular frequency bands is converted and the signals are coded. The total number of signal samples is the same and dependency of samples in decomposed signals

HA, Fh are by a factor

is relatively high. The following suppression of deterministic parts produces error signals with low entropy. The result efficiency is then significantly increased. The dependency of corresponding samples in ECG signal cycles also increases the redundancy. A prediction of the sample from corresponding samples of previous cycles can be done before entropic coding [5]. This approach requires synchronization of heart cycles and can be extended for arrhythmia processing. This can be also modified by subtraction of averaged cycle. The last possible way of increasing the redundancy of multi-lead ECG is a decorrelation of leads. Signals can be decorrelated by Karhunen-Loeve transform and then compressed by subband coding [6]. A comparison of the methods described above is difficult because their efficiency strongly depends on test conditions. Results are influenced by chosen sampling frequency, quantization step, nature of test signal and the error of the reconstructed signal. Lossy compression techniques are successful when the clinical quality of the signals is kept. The aim of this paper is to propose an efficient method for compression of ECG signals. Since our primary field of interest

I. Prouaznik,

J. Kozumplik

/International

Journal

concerns time-frequency digital signal processing, we tried to develop a compression algorithm that would exploit non-stationary characteristics of ECG signals to obtain max-

x(n)

threshold d

,quantization --)

--)

threshold Y

quantization +

Informaficr

45 (1997)

1 ll-

128

DTWT

packing +

threshold d

quantization b

i

Fig. 3. Block diagram of wavelet-based lossy compression.

compressed Huffman ’ decoding --) data

113

imum compression ratio (CR, see Appendix A). The method proposed here is based on the decomposition of the ECG signal into a set of basic functions covering the time-frequency domain. This is also known under the general term of wavelet transform. Through the definition of a direct and an inverse transform both discrete and continuous, the full formulation will be developed for computer implementation. The wavelet transform will be finally applied to ECG data compression and its properties will be illustrated on data from CSE library. We should note that the wavelet transform has been proposed for many application of ECG processing [12] e.g. detection of significant points in ECG signals [13], recognizing cardiac patterns [14], and ventricular late potential analysis [ 151.

Fig. 2. Two-channel analysis (left) and synthesis (right) filter banks as two basic blocks of dyadic discrete-time wavelet transform. The high-pass filter Hh and the low-pass filter Hd decompose the signal x(n) into two subbands that are consequently decimated. The outputs of the analysis block represent wavelet coefficients y,(k) and y,(k + 1). To reconstruct the original signal, the outputs are interpolated and filtered by the highpass filter Fh and the low-pass filter F&

-+

of Medical

run-length decoding -

_)

DTWT-’

unpacking _)

Fig. 4. Block diagram of wavelet-based lossy decompression.

x(n)

114

I. Provaznik,

J. Kozumplik

j International

magnitude I

I

frequency I

Journal

responses 1

of Medical

Injhnatics

IFOI=IHO( and IF1 =IHlI, wavelets I I I

45 (1997)

I I I ~ 128

Meyer I

I

Fig. 5. Magnitude frequency responses of QMF lHhl = /Fhl and lHdl = 1~~1derived from Meyer’s wavelets.

2. Wavelet realization

transform

and its computer

To perform a time-frequency decomposition of a signal, we need to choose a set of basic time functions depending on two parameters: time shift and frequency. A classical approach is a short-time Fourier transform consisting of a Fourier transform computed on a windowed part of the signal. The window length is fixed and it leads to poor time resolution of either high or low frequency components. To overcome this drawback, basic functions called wavelets are used. The functions are defined again by two but different parameters: time shift and time dilation of a unique mother wavelet [7-lo].

To decompose a signal x(t) on the base of functions, the projection of the signal on each wavelet g(e) can be computed. Thus, the wavelet transform (WT) of the signal x(t) is defined as

Eq. (1) represents cross-correlation of the signal to each wavelet, where 7 is time shift, R is time dilation, and star denotes the complex conjugate. J2 is a normalization coefficient which allows each wavelet to keep the same energy. The squared values of y(i,z) represent the energy of a part of the signal around time z in a frequency band related to A.

I. Proaaznik,

J. Kozumplik

magnitude

Fig. 6. Magnitude

/International

frequency

Journal

responses I

1

Injbrmatics

IF01 and (Fll, wavelets I I

blorthogonal I

45 (1997)

11 l-128

a

= x(t)g(2 r 2” s - zc

-“t - k) dt.

(2)

With these notations, the Fourier transform of m-th wavelet is given by

f*(q)]

115

3.5 I

frequency responses of filters (a) IHhl, 117~1,and (b) (FAl, lFdl related to biorthogonal

For practical applications, we are faced with the problem of properly sampling in the time-frequency plane, i.e. selecting the discrete set of wavelet coefficients y(/i,z) to be used. One possibility leads to discrete wavelet transform (DWT) of a continuous-time signal x(t) when /z= ZJ’ and r = $kr, , where A0> 1, r, > 0, m, k are integers. DWT is often proposed with a discretization of the form 3,, = 2, r,, = 1, i.e. /z= 2”, r = 2”k [7]. Then, the DWT is called dyadic wavelet transform: y(m,k)

of Medical

1

wavelets 3/S.

=JFG(2mco)e-jor2mk, (3)

where F[g(t)] = G(o) denotes the Fourier transform of the mother wavelet. DWT leads to an infinite discrete set of coefficients y(m,k) regularly spaced in cycle-octave domain. In other words, the number of time positions of y(m,k) doubles from m to m + 1. From Eq. (3), it follows that time expansion of the mother wavelet corresponds to its spectrum compression in a dyadic way. The technique represented by DWT in Eq. (2) actually consists of applying cross-correlation to wavelets organized in octave bands having a constant relative bandwidth.

I. Prooaznik,

116

Table

J. Kozumplik

/International

Journal

of Medical

Informatics

45 (1997)

I1 l-128

1 Orthogonal

Wordlength [bit] Average PRD [“/;I] Average CR [-] Average BR ACL [bit/s] [bit/sample]

Meyer’s wavelet

Biorthogonal

wavelet 3/5

-I 5.43 18.5

8 4.35 13.0

7 5.25 16.8

8 4.85 (4.43)” 14.0 (14.3)”

3240.65

4620.92

3510.71

4290.86 (420)(0.84>a

Compression efficiency and reconstruction error. Thresholding 30 ,uV in three highest frequency a Denotes results for data set after powerline noise removal (see the text).

bands.

Table 2 Orthogonal Wordlength [bit] Average PRD [‘XI] Average CR [-] Average BR [bit/s] Average ACL [bit/sample] Compression

Meyer’s wavelet

I 3.75 13.1 458 0.92

efficiency and reconstruction

Biorthogonal

8 3.00 9.2 652 1.30

error. Thresholding

wavelet 315

7 4.50 14.0 429 0.86

8 3.94 10.3 583 1.17

30 ,uV in two highest frequency

bands.

To keep possibility of perfect recover of the original signal x(t) from wavelet coefficients, DWT should possesses orthogonality property. Then:

From a computational point of view, dyadic discrete-time wavelet transform (DTWT) must be defined as:

x(t)

y,(k) =

= F- ’ c fi m

c y(m,k)e k

pj’02’“kG(2Mw)

1 .

(4)

We can also write Eq. (2) in term of convolution the signal x(t) with time reverse functions: m x(t)h,(2”k - t) dt, Yhk) = (5) s -cc where 1 h,(2”k - t) = -g(2 - “t - k). J- 2” From Eq. (5) and Eq. (6), DWT consists of applying the signal to a bank of octave band filters with impulse responses h,(O).

=

f x(n)h,(2”k n= --52 f h,(n)x(2”k-n). n= --?i

- n) (7)

Eq. (7) represents an analysis of a discretetime signal x(n) by a bank of digital octave band filters. Coefficients yn?(k) represent the undersampled output of m-th filter when each 2”-th sample is taken, i.e. sampling frequency is 2”-times lower than sampling frequency of x(n). The principal block scheme of dyadic DTWT is shown in Fig. 1. In Fig. 1, one observes a repeating tree structure of two basic blocks of DTWT: two-channel analysis and two-channel synthesis filter banks. The blocks are depicted in Fig. 2.

I. Prooaznik,

J. Kozumplik

/International

Journal

CSE Ilbmy. threshold

of Me&d

30. bns 8. wavelets:

Informatics

45 (1997)

11 I - 128

117

Meyer

0

Fig. 7. Distribution of CR and PRD for CSE library of rest ECG signals- (signals No. 1-12.5, leads I, II, Vl, V2, V3, V4, V5, V6). Compression algorithm used Meyer’s wavelets, threshold 30 PV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Average PRD = 4.350/o, average CR = 13.

For the perfect reconstruction of the original signal, i.e. x(n) = x’(n - r), two conditions have to be fulfilled: F&)&(z) q=Kd

+ Fh(Z)H/JZ) = 22 - ?)

(8)

- z) + f-/z(Z)ffh( -z) = 0,

(9)

where r is time delay of a cascade Hd(z)Fd(z) or IY~(z)F,~(z). Eq. (9) is met if: Fd(z) = Hh( - Z) and Fh(z) = - Hd( - z). (10)

where PJz) is a low-pass filter and Pd( - z) is a high-pass filter with symmetric frequency response. PJz) is a halfband filter and is characterized by /Ham= ,+ = 0.5pqo)(, = 0. To design the two-channel analysis and synthesis filter banks from Fig. 2, a suitable Pd(z) must be selected such as Pd(z) = Hd(z)Fd(-7). Among the various systems Pd(z) proposed in the wavelet literature [8,9], the following halfband low-pass filter may be used for its relative simplicity:

We can also write Eq. (8) as: Fd(Z)&(Z)

Pd(Z) =(- 1 +9zP2+

- Fd( - z)&( - Z)Pd(Z)

-P(,(-Z)=2zrr,

(11)

16~-~+9z-~-~--)/32. (12)

118

I. Prooaznik,

J. Kozumplik

/international

Journal

threshol& I

I

of Medical

30. bits: 8. wavelets. I

Informatics

4.5 (I 997) Ill-

128

b135 I

1

0 0 8

08 0

0

0

0

0

0

0 0

I

I

5

10

I

I

1

20

25

Fig. 8. Distribution of CR and PRD for CSE library of rest ECG signals (signals No. 1-125, leads I, II, Vl, V2, V3, V4, V5, V6). Compression algorithm used biorthogonal wavelets 3/5, threshold 30 ,uV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Average PRD = 4.85%, average CR = 14.

This integer-coefficient filters can be divided into a cascade of filters: fW) =(-

Hh(z) = Hd( - z).

(14)

Inserting Eq. (14) into Eq. (10) we have F&l = H&X

1 +2~-‘+6~-~+2-7-~-~--4)/8,

F,Jz) = - Hh(z).

(15)

Due to the orthogonality property, impulse responses of the filters are reverse. Thus, for the impulse response ,{h,(n), n = 0, 1, . . ., N l> of the filter Hd(z), we have

Fh(z)=(1+2z-1-6z-2+2z-3+z-4)/8. (13) As a special case, quadrature mirror filters (QMF) can be used. Then Hd(z) and H,,(z) have the same magnitude frequency responses but folded at o = n/2:

FJZ) = z -(A’- l)&(zFh(z) = zrcN-‘)Hh(zrl),

I), (16)

and Eq. (8) reduces after some straightforward calculations to

I. Prouaznik,

J. Kozumplik

/international filtered I

I 5

I 10

Journal

CSE Ikbrary, threshold. ,

of Medical

informatics

30. bits 8. wevelets: I

I 15

, 20

45 (1997)

Ill-

128

119

bi35 I

I

I 25

I 30

I 35

PRD P4

Fig. 9. Distribution of CR and PRD for filtered CSE library of rest ECG signals (signals No. l- 125, leads I, II, Vl, V2, V3, V4, V5, V6). Powerline noise suppressed by a notch filter at 60 and 180 Hz. Compression algorithm used biorthogonal wavelets 315, threshold 30 PV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Average PRD = 4.43%, average CR = 14.3

-z)fqj-zp’)=2 (17)

The simplest QMFs are defined as: f&(z) = (1 + 2 - 9/z, Fd(Z) = (1 + Z- ‘w, l&(z) = (1 -z-

‘)/2,

Fh(Z) = - (1 - z ~ ‘)/2.

(18) The bank of filters Eq. (18) represents dyadic DTWT with Haar wavelets. Practical problems emerge when using the wavelet approach to analyze a discrete-time signal of finite duration. When computing the convolution formula Eq. (7) over such sig-

nals, distortions occur at the beginning and at the end of the signals, mainly due to the low-frequency wavelets. To overcome this problem, periodic orthogonal wavelets are employed [lo]. The basic concept relies on the use of DTWT obtained e.g. by cyclic convolution in frequency domain:

~(4 = DFT,&mMWL,(~)1,

(19)

where X(u) = DFT, [x(n)], DFT is discrete Fourier transform, m = 1,2, . . .. A4, k= 0, 1, . ..) N/2”I,u=O 713 *.., N- 1, M is a number of wavelets, N = 2M is duration of the signal and also a number of generated wavelet coefficients; and H,(u)

= J2mG*(2”‘v)

(20)

120

I. Provaznik,

J. Kozumplik

/International

Journal

of‘ Medical

Infbrmatics

4.5 (1997)

I I I 128

1 16

I 1.8

origmalsrgnalwO06.vZ

I 0

1 02

I

I 04

decompressedsl I

I 06

I

I 08

gnal,threshold:30,~te:8, I

I 1 t [set]

PRD=2

‘n

I 12

I 14

643,CR=12 I

89,waveiets I

'

A'

/

c

:; 08

I

1

t [set]

I

h~eyer

500 -

I

I 2

I

I

I

12

1.4

16

I

I -

18

Fig. 10. Original signal No. 6, lead V2 (a) and related reconstructed signal (b). The compression algorithm used Meyer’s wavelets, threshold 30 PV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Resulted PRD = 2.64X, CR = 12.9.

is a frequency response of m-th filter related to m-th wavelet. Subscript N,N/2” in Eq. (19) represents N-point inverse discrete Fourier transform (DFT- ‘) when each 2”-th sample is taken. Bandwidth of m-th filter is 2 - “fs 12, where Ji is sampling frequency, and its output signal is fV, = 2 - “N samples long. Inverse DTWT is defined by: x(n) = DFT-

‘[DFTN,,n,,~~,~(k>]Hm(v)l

+ Y(O)/N. (21)

A set of orthogonal wavelets can also be constructed in frequency domain by sampling their spectra. Detailed description of periodical orthogonal Meyer wavelets is in [lo].

3. Wavelet-based signals

compression

3.1. The compression

of ECG

concept

The wavelet transform is a valuable technique for lossy compression of signals having high-frequency components of relatively short duration. Such signals are, for instance, ECG records. Their QRS complexes contain the highest frequency components and cover about 10% of each heart cycle. This strictly determines the lowest sampling frequencyusually 500 Hz. Parts of the ECG signal between QRS complexes are five or even more times oversampled [ 111. Using dyadic wavelet approach to analyze the signal allows

I. Prot’aznik,

J. Kozumplik

/International

Journal

of Medical

Informatics

121

45 (1997)

1 I I - 128

16

1.8

2

1.8

2

ongmalsignalwO06.v2

0

02

I

I

04

06

I

0

1 02

1 t [set]

12

1.4

decompressedslgnal,threshold:3O,bits:8,PRD=2694,CR=1457,wa~lets 1

'

-2ooot

08

I 04

I 06

bi35

';.:I:+

08

1

t [sac]

12

Fig. 11. Original signal No. 6, lead V2 (a) and related reconstructed biorthogonal wavelets 315, threshold 30 PV in three highest frequency 8 bits. Resulted PRD = 2.69X, CR = 14.6.

decomposition into frequency bands with perfect time localization of high-frequency components. Ideally, the highest bands (represented by wavelet coefficients y,(k) for the lowest m) consist of long series of zeroes alternated with short sequences of samples related to QRS complexes. In practice, the long series consist of ‘almost’ zero-valued samples. From the above, the following compression scheme is proposed. Firstly, the original signal x(n) is decomposed by dyadic DTWT. The wavelet coefficients yrn(k) are thresholded, i.e. y,(k) of value lower than certain threshold are set to zero. This later leads to efficient coding. Wavelets coefficients are gen-

14

16

signal (b). The compression algorithm used bands, quantization of wavelet coefficients to

erally real-valued numbers. The next step is then quantization. At this point, the data are prepared for the following lossless coding that increase compression efficiency. All wavelet coefficients are packed into one sequence. Due to described properties of am, the zero-valued samples can be efficiently stored by a run-length coder followed by a standard entropic coder. To recover the original data, the discussed steps of compression scheme are applied in reverse order: Huffman decoding, run-length decoding, data unpacking, and inverse wavelet transform. The block diagram of the compression and decompression schemes are in Figs. 3 and 4.

I. Provaznik,

122

J. Kozumplik

/International

Journal

WOO6 v2 slgnel decomposition.

-4OC

-50 -

I 50

I 100

I 150

I 200

of Medical lhreshold:30.

I 250

Informatics wavelets-

I 300

45 (1997)

1 I1 ~ 128

Meyer

I 350

I 400

1 450

1 500

"

-100 I 50

I 100

I 150

1 200

I 250

200 100 0 -100 -200

Y

I 20

I 40

60

I 80

I 100

I 120

Fig. 12. Three highest frequency bands: (a) 125-250 Hz, (b) 62.5-125 Hz, (c) 31.25-62.5 Hz of the decomposed original signal No. 6, lead V2. Horizontal lines represent used thresholds 30 ,uV. The compression algorithm used Meyer’s wavelets, quantization of wavelet coefficients to 8 bits. Resulted PRD = 2.64%, CR = 12.9.

3.2. The test ECG signals

3.3. The wavelet transform

Compression efficiency of the proposed wavelet-based algorithm has been tested on common standard in electrocardiology (CSE) library of rest ECG signals. The library consists of 125 10-s records. Eight leads (I, II, Vl, V2, V3, V4, V5, V6) from each record have been chosen according to usual standards for ECG data archiving. The signals have been sampled at frequency of 500 Hz with quantization step of 5 PV and represented by integer numbers. 12-bit wordlength has been assumed as suitable for data storing.

For our first experiments, Meyer’s periodic wavelets have been used. The generating Meyer wavelet g(t) is defined [lo] by its Fourier transform G(f) as follows: IG(f)l = 0 for (fl~[O; 4 - e] (G(f)1 = SW--4) for ]fl~G - E;$ + E] lG(f)l = 1 for lfl~[$-IG(f)l=S(

-g&i)

E; 1 - 2~1 for If]~[l-22~;1+2~]

I. Provaznik,

J. Kozumplik

/International

wUO6v2

Journal

of Medical

Informatics

s~gnaldecomposi~on,tkeshold:30,wavelets:

45 (1997)

Ill-

128

123

bt35 I

-100 I 50

-500

I 100

I 20

I 40

I 150

, 60

I 250

I 200

1 80

I 100

I 120

Fig. 13. Three highest frequency bands: (a) 125-250 Hz, (b) 62.5-125 Hz, (c) 31.25-62.5 Hz of the decomposed original signal No. 6, lead V2. Horizontal lines represent used thresholds 30 pV. The compression algorithm used biorthogonal wavelets 3/5, quantization of wavelet coefficients to 8 bits. Resulted PRD = 2.69%, CR = 14.6.

arg G(f) = - rcf

onal wavelets 3/5. The discrete-time generating wavelet g(n) is defined in time domain [9] according to Eq. (13). Magnitude frequency responses related to mentioned biorthogonal 3/5 wavelets are in Fig. 6.

where

W> = df) for fW 4 Kf) = $733 forfd S(f) = l/G

5 01

3.4. Thresholding

for f= 0

and

r(f) = 1 - expb*/(f-

~1~1,

a = log[l - l/,:2], 0 < EI l/6. Magnitude frequency responses of QMFs derived from Meyer’s wavelets are in Fig. 5. Also, the tests have been done with biorthog-

Of the various methods to modify the wavelet coefficients (such as hard thresholding, soft thresholding, quantile thresholding, universal thresholding etc.), soft thresholding has been chosen. Using this method, values, of the wavelet coefficients are shifted to zero level in selected high frequency bands. We thresholded the coefficients in two and three frequency highest bands (125-250, 62.5-125,

124

I. Provaznik,

J. Koxmplik

/International

Journal

of Medical

origmal slgnal

Informatics

45 (1997)

11 I ~ 128

WOO6 i

250

decompressed I

slgnel,~reshold:30,b~ts I I

8, PRD=14 I

17,CR=ll I

45,wevelets I

Meyer I

Fig. 14. Original signal No. 6, lead I, distorted by powerline noise (a) and related reconstructed signal (b). The compression algorithm used Meyer’s wavelets, threshold 30 ,uV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Resulted PRD = 14.17%, CR = 11.5.

and 31.25-62.5 Hz for sampling frequency 500 Hz). All bands have been thresholded at the same level of 30 ,uV. This value has been chosen as a compromise according to compression efficiency and error of reconstruction. 3.5. Quantization

and run-length

coding

The modified wavelet coefficients have been stored with wordlength of 7 and 8 bits. For lower values, the error of compression dramatically increased having a negative effect on compression efficiency. Zero-valued coefficients are stored by a run-length coder. Relatively long sequences of coefficients related to parts of the original signal between

QRS complexes are then represented by a two-number code and thus efficiently stored. 3.6. En tropic coding

A Huffman entropic coder has been used for further compression. This lossless coder increased compression ratio 1.7 times in average for the test data. The best results have been achieved by the Huffman coder using no predictor. 4. Results and discussion The proposed wavelet-based compression algorithm has been tested for the following

I. Prouaznik,

J. Kozumplik

/International

Journal

of Medical

Informatics

45 (1997)

111~ I28

16

18

125

ongmalslgnalw006.1

300 250 200 7 E 150 f > 100 50

0

02

I

I

04

06

08

1 t [set]

12

14

decompresseds~gnal,threshold30,b~ts8,PRD=1512,CR=1283,wavelets I I I I I

I

b135 I

2

I

250

I

I

I

/

t

0

02

04

06

08

I

1 t [set]

I

I

I

12

14

16

I

1.8

2

Fig. 15. Original signal No. 6, lead I, distorted by powerline noise (a) and related reconstructed signal (b). The compression algorithm used biorthogonal wavelets 3/5, threshold 30 ,uV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Resulted PRD = 15.12%, CR = 12.8.

parameters: Meyer’s and biorthogonal 3/5 generating wavelets, thresholding in two and three highest frequency bands with threshold of 30 pV, quantization with 7 and 8 bit resolution. Compression efficiency has been evaluated by compression ratio (CR), bit rate (BR), and average code length (ACL), (see Appendix A). Error of reconstruction has been evaluated by percent root-mean-square difference (PRD), (see Appendix B). For example, one can choose biorthogonal 3/5 wavelets, with 8 bit resolution, and thresholding in three bands. Then, resulted compression ratio CR for tests on all 1000 signals CR = 14.0 (BR = 429 bit/s, ACL = 0.86 bit/ sample) when reconstruction error PRD = 4.85%. Results of other tests are in Table 1

and Table 2. One should note PRD and CR in both tables are average results. Due to distortion of some records by powerline noise and/or myopotentials, PRD can significantly arise as the noise is suppressed by thresholding. The remaining noise shorten sequences of zerovalued samples in particular bands and CR is then typically low. This is well visible in Figs. 7 and 8 in which distribution of resulting CR and PRD for all tested signals is shown. The results for the compression algorithm using Meyer’s wavelets are in Fig. 7, for the algorithm using biorthogonal wavelets 3/5 in Fig. 8. The distribution of CR and PRD for a set of test signals with removed 60 Hz powerline noise is in Fig. 9.

126

I. Pmvaznik,

J. Kozumplik

/ Imternational

Journal ongmal

decompressed

0

02

04

signal, threshoM30.

0.6

08

of Medical

Informatics

45 (I 997) I1 1- 128

signal wUO6f I

b&8.

PRD=B

1 t [set]

849, CR=16

55. wavelets

1.2

1.4

b135

16

18

2

Fig. 16. Original signal No.6, lead I, filtered by a notch filter at 60 and 180 Hz to suppress powerline noise (a) and related reconstructed signal (b). The compression algorithm used biorthogonal wavelets 3/5, threshold 30 PV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Resulted PRD = 8.85%. CR = 16.6.

In Figs. 10 and 11, examples of the original signal No. 6, lead V2, from CSE library and the related reconstructed signal for the compression algorithm using Meyer’s wavelets (Fig. 10) and for the algorithm using biorthogonal wavelets 3/5 (Fig. 11) are shown. Fig. 12 and Fig. 13 corresponding to Figs. 10 and 11, respectively, show three highest frequency bands of the decomposed original signal where horizontal lines represent the used thresholds. An example of the original signal distorted by powerline noise (CSE signal No. 6, lead I) and the related reconstructed signal is in Fig. 14 and Fig. 15. Fig. 14 represents coding by the compression algorithm using Meyer’s wavelets. Fig. 15 represents coding by the

algorithm using biorthogonal wavelets 3/5. Original signal No. 6, lead I, filtered by a notch filter at 60 and 180 Hz to suppress powerline noise and related reconstructed signal are shown in Fig. 16. The compression algorithm used biorthogonal wavelets 3/5, threshold 30 ,uV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Thresholding in three highest frequency bands can sporadically disturb shape of Pwaves when their spectral contents exceeds 30 Hz. In order to avoid this shortcoming, we propose thresholding in two highest frequency bands resulted in lower CR (see Table 2). Another solution consists of subdividing the third frequency band (31.25-62.5 Hz)

I. Provaznik,

J. Kozumplik

/International

Journal

into two bands of same bandwidth and thresholding in a higher band (46.875-62.5 Hz). Then, CR would vary between values from Tables 1 and 2. For modification of wavelet coefficients, we have proposed soft thresholding using an empirically derived threshold. It should be noted that the use of optimum thresholding where threshold value is derived from noise characteristics can increase compression efficiency. Thus, the algorithm using such thresholding can minimize distortion of QRS complexes. More sophisticated algorithms could employ a QRS complex detector and modify the wavelet coefficients between QRS complexes. For the best results, we have proposed to remove powerline noise before compression. Further PRD reduction could be reached by using different quantization in particular frequency bands. Also, the results are influenced by the selection of the used wavelets. 5. Conclusions In this paper, we propose to use a lossy compression algorithm for ECG data processing. The algorithm is based on orthogonal or biorthogonal wavelet transforms. Unimportant wavelet coefficients are set to zero by soft thresholding (30 pV) in two or three highest frequency bands. Thresholding and quantization to 7 or 8 bits combined with Huffman coding leads to average code length from 0.65 to 1.3 bit per sample with corresponding relative root-mean-square error (PRD) between 5.43 and 3.00%. This compression is about 3-5 times more efficient than direct lossless entropic coding. The proposed method was tested on the CSE ECG signal library.

of Medical

Informaties

45 (1997)

I II-

128

127

Acknowledgements This project was partly supported by the grant No. 102/94/0374 of the Czech State Grant Agency.

Appendix A. Compression efficiency Compression ratio (CR) is defined by the formula: CR =

amount of original data [bit] amount of compressed data [bit] L-1 (23)

Bit rate (BR) is defined by the formula: BR = amount of compressed data [bit] @its - ‘1 original signal duration [s] (24) Average code length (ACL) is defined by the formula: ACL = amount of compressed data [bit] amount of original data [sample] [bit.sample - ‘1

(25)

Appendix B. Error of compression The error of reconstruction is characterized difference by percent root-mean-square (PRD) by the formula:

ccx