VIEWS: 77 PAGES: 18 CATEGORY: Alternative Medicine POSTED ON: 6/1/2010 Public Domain
lntam~aMl Jaumal Cal Medical lnformatics ELSEVIER International Journal of Medical Informatics 45 (1997) 11 l- 128 Wavelet transform in electrocardiography-data compression I. Provaznik *, J. Kozumplik Department of‘ Biomedical Engineering, Faculty of’ Elecrrical Engineering and Computer Science, Technical University Brno, Purkynova 9la, 612 00 Bmo, Czech Republic 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; Huffman coding 1. Introduction signal through common communication tele- phone line of 4800 bps, the data must be Compression algorithms are proposed for reduced. Also, memory requirements for stor- real-time data communication and/or for ing one 10-s 12-lead ECG record digitized as low-cost data storing. Electrocardiographical described above are 90 kb. (ECG) signals are usually sampled at fre- The ECG signals digitized as described quency 500 Hz and digitized to 12 bits in l-8 above are redundant. The redundancy fol- channels that represents data flow of 6000- lows from the fact that distribution of signal 48 000 bps. To transfer an one-lead ECG values is not uniform. The optimum method to suppress redundancy is entropic coding if the samples are statistically independent. The *Corresponding author. Tel.: + 42 5 759310; fax: f42 5 746757; e-mail: provazni@dbme.fee.vutbr.cz alternatively most familiar method is Huffman code that is kozumpli@dbme.fee.vutbr.cz web: http://www.fee.vutbr/cz used in many versions in almost all published 1386.5056:97/$17.00 Q 1997 Elsevier Science Ireland Ltd. All rights reserved PII SOO20-7101(97)00040-8 112 I. Prol;aznik, J. Kozumplik /International Journal of Medical Informatics 45 (1997) 111~ 128 x(n) YIO x’(d) Hh -‘12 . ‘12 ---) Ydd Hd ‘12 Hh -‘12 ... Fig. 1. Block diagram of dyadic wavelet transform (left) and its associated inverse transform (right). HA, Fh are high-pass filters, Hd, Fd are low-pass filters, 12 denotes decimation by a factor 2, T2 denotes interpolation by a factor 2. compression systems. Nevertheless, entropic is relatively high. The following suppression coding is not sufficient compression scheme of deterministic parts produces error signals because of high entropy of ECG signals. The with low entropy. The result efficiency is then entropy can be a posteriori decreased in sev- significantly increased. eral ways. A deterministic part of signals can The dependency of corresponding samples be removed by predictive, interpolate and in ECG signal cycles also increases the redun- differentiate estimation which are optimum in dancy. A prediction of the sample from cor- minimization of residual error signal entropy responding samples of previous cycles can be [l]. A more successful method to decrease the done before entropic coding [5]. This ap- entropy of ECG signal is using extensive proach requires synchronization of heart cy- Markov process which considers both the cles and can be extended for arrhythmia coding redundancy and the intersample re- processing. This can be also modified by sub- dundancy [2]. The average code length for the traction of averaged cycle. extensive Markov system on the 2”d differ- The last possible way of increasing the ence signals is about 2.5 bit/s while average redundancy of multi-lead ECG is a decorrela- Huffman code length for 2”d difference sig- tion of leads. Signals can be decorrelated by nals is about 3.3 bit/s [2]. The described Karhunen-Loeve transform and then com- entropic coding of residuals can be without pressed by subband coding [6]. loss, i.e. the compressed signal can be fully A comparison of the methods described recovered. above is difficult because their efficiency An interesting approach is subband coding strongly depends on test conditions. Results that was inspired by speech processing. The are influenced by chosen sampling frequency, first articles describe a division into six octave quantization step, nature of test signal and bands [3], or into four bands of the same the error of the reconstructed signal. Lossy bandwidth [4]. The sampling frequency of compression techniques are successful when signal of particular frequency bands is con- the clinical quality of the signals is kept. verted and the signals are coded. The total The aim of this paper is to propose an number of signal samples is the same and efficient method for compression of ECG dependency of samples in decomposed signals signals. Since our primary field of interest I. Prouaznik, J. Kozumplik /International Journal of Medical Informaficr 45 (1997) 1 ll- 128 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-fre- quency domain. This is also known under the Fig. 2. Two-channel analysis (left) and synthesis (right) general term of wavelet transform. Through filter banks as two basic blocks of dyadic discrete-time the definition of a direct and an inverse trans- wavelet transform. The high-pass filter Hh and the form both discrete and continuous, the full low-pass filter Hd decompose the signal x(n) into two formulation will be developed for computer subbands that are consequently decimated. The outputs of the analysis block represent wavelet coefficients implementation. The wavelet transform will y,(k) and y,(k + 1). To reconstruct the original signal, be finally applied to ECG data compression the outputs are interpolated and filtered by the high- and its properties will be illustrated on data pass filter Fh and the low-pass filter F& from CSE library. We should note that the wavelet transform has been proposed for concerns time-frequency digital signal pro- many application of ECG processing [12] e.g. cessing, we tried to develop a compression detection of significant points in ECG signals algorithm that would exploit non-stationary [13], recognizing cardiac patterns [14], and characteristics of ECG signals to obtain max- ventricular late potential analysis [ 151. -+ threshold d ,quantization --) x(n) --) threshold Y quantization + DTWT packing + threshold d quantization - b i Fig. 3. Block diagram of wavelet-based lossy compression. _) compressed Huffman run-length x(n) ’ decoding --) decoding - unpacking DTWT-’ data _) Fig. 4. Block diagram of wavelet-based lossy decompression. 114 I. Provaznik, J. Kozumplik j International Journal of Medical Injhnatics 45 (1997) I I I ~ 128 magnitude frequency responses IFOI=IHO( and IF1 =IHlI, wavelets Meyer I I I 1 I I I I I Fig. 5. Magnitude frequency responses of QMF lHhl = /Fhl and lHdl = 1~~1derived from Meyer’s wavelets. 2. Wavelet transform and its computer To decompose a signal x(t) on the base of realization functions, the projection of the signal on each wavelet g(e) can be computed. Thus, the To perform a time-frequency decomposi- wavelet transform (WT) of the signal x(t) is tion of a signal, we need to choose a set of defined as basic time functions depending on two parameters: time shift and frequency. A clas- sical approach is a short-time Fourier trans- form consisting of a Fourier transform computed on a windowed part of the signal. Eq. (1) represents cross-correlation of the The window length is fixed and it leads to signal to each wavelet, where 7 is time shift, R poor time resolution of either high or low is time dilation, and star denotes the complex frequency components. To overcome this conjugate. J2 is a normalization coefficient drawback, basic functions called wavelets are which allows each wavelet to keep the same used. The functions are defined again by two energy. The squared values of y(i,z) repre- but different parameters: time shift and time sent the energy of a part of the signal around dilation of a unique mother wavelet [7-lo]. time z in a frequency band related to A. I. Proaaznik, J. Kozumplik /International Journal of Medical Injbrmatics 45 (1997) 11 l-128 115 magnitude frequency responses IF01 and (Fll, wavelets blorthogonal 3.5 I I I I I 1 Fig. 6. Magnitude frequency responses of filters (a) IHhl, 117~1,and (b) (FAl, lFdl related to biorthogonal wavelets 3/S. For practical applications, we are faced with the problem of properly sampling in the f*(q)] =JFG(2mco)e-jor2mk, time-frequency plane, i.e. selecting the dis- (3) crete set of wavelet coefficients y(/i,z) to be used. One possibility leads to discrete wavelet where F[g(t)] = G(o) denotes the Fourier transform (DWT) of a continuous-time signal transform of the mother wavelet. DWT leads x(t) when /z= ZJ’ and r = $kr, , where A0> to an infinite discrete set of coefficients y(m,k) regularly spaced in cycle-octave do- 1, r, > 0, m, k are integers. DWT is often proposed with a discretization of the form main. In other words, the number of time 3,, = 2, r,, = 1, i.e. /z= 2”, r = 2”k [7]. Then, positions of y(m,k) doubles from m to m + 1. the DWT is called dyadic wavelet transform: From Eq. (3), it follows that time expansion of the mother wavelet corresponds to its 1 a spectrum compression in a dyadic way. The y(m,k) = - x(t)g(2 -“t - k) dt. (2) technique represented by DWT in Eq. (2) r 2” s - zc actually consists of applying cross-correlation With these notations, the Fourier transform to wavelets organized in octave bands having of m-th wavelet is given by a constant relative bandwidth. 116 I. Prooaznik, J. Kozumplik /International Journal of Medical Informatics 45 (1997) I1 l-128 Table 1 Orthogonal Meyer’s wavelet Biorthogonal wavelet 3/5 Wordlength [bit] -I 8 7 8 Average PRD [“/;I] 5.43 4.35 5.25 4.85 (4.43)” Average CR [-] 18.5 13.0 16.8 14.0 (14.3)” ACL [bit/sample] Average BR [bit/s] 0.65 324 462 0.92 0.71 351 0.86 (0.84>a 429 (420) Compression efficiency and reconstruction error. Thresholding 30 ,uV in three highest frequency bands. a Denotes results for data set after powerline noise removal (see the text). Table 2 Orthogonal Meyer’s wavelet Biorthogonal wavelet 315 Wordlength [bit] I 8 7 8 Average PRD [‘XI] 3.75 3.00 4.50 3.94 Average CR [-] 13.1 9.2 14.0 10.3 Average BR [bit/s] 458 652 429 583 Average ACL [bit/sample] 0.92 1.30 0.86 1.17 Compression efficiency and reconstruction error. Thresholding 30 ,uV in two highest frequency bands. To keep possibility of perfect recover of From a computational point of view, the original signal x(t) from wavelet coeffi- dyadic discrete-time wavelet transform cients, DWT should possesses orthogonality (DTWT) must be defined as: property. Then: x(t) y,(k) = f x(n)h,(2”k - n) n= --52 = F- ’ c fi m c y(m,k)e k pj’02’“kG(2Mw) 1 . (4) = f n= --?i h,(n)x(2”k-n). (7) We can also write Eq. (2) in term of convo- Eq. (7) represents an analysis of a discrete- lution the signal x(t) with time reverse func- time signal x(n) by a bank of digital octave tions: band filters. Coefficients yn?(k) represent the m undersampled output of m-th filter when Yhk) = x(t)h,(2”k - t) dt, (5) each 2”-th sample is taken, i.e. sampling fre- s -cc quency is 2”-times lower than sampling fre- where quency of x(n). 1 The principal block scheme of dyadic h,(2”k - t) = -g(2 - “t - k). DTWT is shown in Fig. 1. In Fig. 1, one J- 2” observes a repeating tree structure of two From Eq. (5) and Eq. (6), DWT consists of basic blocks of DTWT: two-channel analysis applying the signal to a bank of octave band and two-channel synthesis filter banks. The filters with impulse responses h,(O). blocks are depicted in Fig. 2. I. Prooaznik, J. Kozumplik /International Journal of Me&d Informatics 45 (1997) 11 I - 128 117 CSE Ilbmy. threshold 30. bns 8. wavelets: 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 orig- where PJz) is a low-pass filter and Pd( - z) inal signal, i.e. x(n) = x’(n - r), two condi- is a high-pass filter with symmetric fre- tions have to be fulfilled: quency response. PJz) is a halfband filter and is characterized by /Ham= ,+ = F&)&(z) + Fh(Z)H/JZ) = 22 - ?) (8) 0.5pqo)(, = 0. q=Kd - z) + f-/z(Z)ffh( -z) = 0, (9) To design the two-channel analysis and synthesis filter banks from Fig. 2, a suitable where r is time delay of a cascade Pd(z) must be selected such as Pd(z) = Hd(z)Fd(z) or IY~(z)F,~(z). Eq. (9) is met if: Hd(z)Fd(-7). Among the various systems Pd(z) proposed in the wavelet literature Fd(z) = Hh( - Z) and Fh(z) = - Hd( - z). [8,9], the following halfband low-pass filter (10) may be used for its relative simplicity: We can also write Eq. (8) as: Pd(Z) Fd(Z)&(Z) - Fd( - z)&( - Z)Pd(Z) =(- 1 +9zP2+ 16~-~+9z-~-~--)/32. -P(,(-Z)=2zrr, (11) (12) 118 I. Prooaznik, J. Kozumplik /international Journal of Medical Informatics 4.5 (I 997) Ill- 128 threshol& 30. bits: 8. wavelets. b135 I I I I 1 0 0 08 0 0 0 0 0 8 0 0 0 I I I I 1 5 10 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 Hh(z) = Hd( - z). (14) into a cascade of filters: Inserting Eq. (14) into Eq. (10) we have fW) 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 - Fh(z)=(1+2z-1-6z-2+2z-3+z-4)/8. l> of the filter Hd(z), we have (13) FJZ) = z -(A’- l)&(z- I), As a special case, quadrature mirror filters Fh(z) = zrcN-‘)Hh(zrl), (16) (QMF) can be used. Then Hd(z) and H,,(z) have the same magnitude frequency re- and Eq. (8) reduces after some straightfor- sponses but folded at o = n/2: ward calculations to I. Prouaznik, J. Kozumplik /international Journal of Medical informatics 45 (1997) Ill- 128 119 filtered CSE Ikbrary, threshold. 30. bits 8. wevelets: bi35 I , I I I I I I , I I I 5 10 15 20 25 30 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 nals, distortions occur at the beginning and (17) at the end of the signals, mainly due to the low-frequency wavelets. To overcome this The simplest QMFs are defined as: problem, periodic orthogonal wavelets are f&(z) = (1 + 2 - 9/z, employed [lo]. The basic concept relies on the use of DTWT obtained e.g. by cyclic convo- Fd(Z) = (1 + Z- ‘w, lution in frequency domain: l&(z) = (1 -z- ‘)/2, ~(4 = DFT,&mMWL,(~)1, (19) Fh(Z) = - (1 - z ~ ‘)/2. (18) where X(u) = DFT, [x(n)], DFT is discrete Fourier transform, m = 1,2, . . .. A4, k= The bank of filters Eq. (18) represents dyadic 0, 1, . ..) N/2”- I,u=O 713 *.., N- 1, M is a DTWT with Haar wavelets. number of wavelets, N = 2M is duration of Practical problems emerge when using the the signal and also a number of generated wavelet approach to analyze a discrete-time wavelet coefficients; and signal of finite duration. When computing the convolution formula Eq. (7) over such sig- H,(u) = J2mG*(2”‘v) (20) 120 I. Provaznik, J. Kozumplik /International Journal of‘ Medical Infbrmatics 4.5 (1997) I I I 128 origmalsrgnalwO06.vZ I 1 I I I I I I 1 I I I 0 02 04 06 08 1 12 14 16 1.8 2 t [set] decompressedsl gnal,threshold:30,~te:8, PRD=2 643,CR=12 89,waveiets h~eyer I I I I I I I c ‘n ' A' I 500 - / :; I I I I I - 08 1 12 1.4 16 18 t [set] 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 3. Wavelet-based compression of ECG to m-th wavelet. Subscript N,N/2” in Eq. signals (19) represents N-point inverse discrete Fourier transform (DFT- ‘) when each 2”-th 3.1. The compression concept sample is taken. Bandwidth of m-th filter is 2 - “fs 12, where Ji is sampling frequency, and The wavelet transform is a valuable tech- its output signal is fV, = 2 - “N samples long. nique for lossy compression of signals having Inverse DTWT is defined by: high-frequency components of relatively short duration. Such signals are, for instance, x(n) ECG records. Their QRS complexes contain = DFT- ‘[DFTN,,n,,~~,~(k>]Hm(v)l + Y(O)/N. the highest frequency components and cover about 10% of each heart cycle. This strictly (21) determines the lowest sampling frequency- A set of orthogonal wavelets can also be usually 500 Hz. Parts of the ECG signal constructed in frequency domain by sampling between QRS complexes are five or even their spectra. Detailed description of periodi- more times oversampled [ 111. Using dyadic cal orthogonal Meyer wavelets is in [lo]. wavelet approach to analyze the signal allows I. Prot’aznik, J. Kozumplik /International Journal of Medical Informatics 45 (1997) 1 I I - 128 121 ongmalsignalwO06.v2 0 02 04 06 08 1 12 1.4 16 1.8 2 t [set] decompressedslgnal,threshold:3O,bits:8,PRD=2694,CR=1457,wa~lets bi35 I I 1 ' ';.:I:+ -2ooot I 1 I I 0 02 04 06 08 1 12 14 16 1.8 2 t [sac] Fig. 11. Original signal No. 6, lead V2 (a) and related reconstructed signal (b). The compression algorithm used biorthogonal wavelets 315, threshold 30 PV in three highest frequency bands, quantization of wavelet coefficients to 8 bits. Resulted PRD = 2.69X, CR = 14.6. decomposition into frequency bands with erally real-valued numbers. The next step is perfect time localization of high-frequency then quantization. At this point, the data are components. Ideally, the highest bands (rep- prepared for the following lossless coding resented by wavelet coefficients y,(k) for the that increase compression efficiency. All lowest m) consist of long series of zeroes wavelet coefficients are packed into one se- alternated with short sequences of samples quence. Due to described properties of am, related to QRS complexes. In practice, the the zero-valued samples can be efficiently long series consist of ‘almost’ zero-valued stored by a run-length coder followed by a samples. standard entropic coder. From the above, the following compres- To recover the original data, the discussed sion scheme is proposed. Firstly, the original steps of compression scheme are applied in signal x(n) is decomposed by dyadic DTWT. reverse order: Huffman decoding, run-length The wavelet coefficients yrn(k) are threshol- decoding, data unpacking, and inverse ded, i.e. y,(k) of value lower than certain wavelet transform. The block diagram of the threshold are set to zero. This later leads to compression and decompression schemes are efficient coding. Wavelets coefficients are gen- in Figs. 3 and 4. 122 I. Provaznik, J. Kozumplik /International Journal of Medical Informatics 45 (1997) 1 I1 ~ 128 WOO6 v2 slgnel decomposition. lhreshold:30. wavelets- Meyer -4OC I I I I I I I I 1 1 50 100 150 200 250 300 350 400 450 500 -50 - " -100 - I I I 1 I 50 100 150 200 250 200 100 0 -100 -200 Y I I I I I 20 40 60 80 100 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 For our first experiments, Meyer’s periodic wavelet-based algorithm has been tested on wavelets have been used. The generating common standard in electrocardiology (CSE) Meyer wavelet g(t) is defined [lo] by its library of rest ECG signals. The library con- Fourier transform G(f) as follows: sists of 125 10-s records. Eight leads (I, II, Vl, V2, V3, V4, V5, V6) from each record IG(f)l = 0 for (fl~[O; 4 - e] have been chosen according to usual stan- dards for ECG data archiving. The signals (G(f)1 = SW--4) for ]fl~G - E;$ + E] have been sampled at frequency of 500 Hz with quantization step of 5 PV and repre- lG(f)l = 1 for lfl~[$-- E; 1 - 2~1 sented by integer numbers. 12-bit word- length has been assumed as suitable for data IG(f)l=S( -g&i) for If]~[l-22~;1+2~] storing. I. Provaznik, J. Kozumplik /International Journal of Medical Informatics 45 (1997) Ill- 128 123 wUO6v2 s~gnaldecomposi~on,tkeshold:30,wavelets: bt35 I -100 I I I I I 50 100 150 200 250 -500 , 1 I I I I 20 40 60 80 100 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 generat- ing wavelet g(n) is defined in time domain [9] where according to Eq. (13). Magnitude frequency responses related to mentioned biorthogonal W> = df) for fW 4 3/5 wavelets are in Fig. 6. Kf) = $733 forfd 5 01 3.4. Thresholding S(f) = l/G for f= 0 Of the various methods to modify the and wavelet coefficients (such as hard threshold- ing, soft thresholding, quantile thresholding, r(f) = 1 - expb*/(f- ~1~1, universal thresholding etc.), soft thresholding a = log[l - l/,:2], 0 < EI l/6. has been chosen. Using this method, values, of the wavelet coefficients are shifted to zero Magnitude frequency responses of QMFs level in selected high frequency bands. We derived from Meyer’s wavelets are in Fig. 5. thresholded the coefficients in two and three Also, the tests have been done with biorthog- frequency highest bands (125-250, 62.5-125, 124 I. Provaznik, J. Koxmplik /International Journal of Medical Informatics 45 (1997) 11 I ~ 128 origmal slgnal WOO6 i 250 decompressed slgnel,~reshold:30,b~ts 8, PRD=14 17,CR=ll 45,wevelets Meyer I I I I I I 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 QRS complexes are then represented by a 500 Hz). All bands have been thresholded at two-number code and thus efficiently stored. the same level of 30 ,uV. This value has been chosen as a compromise according to com- 3.6. En tropic coding pression efficiency and error of reconstruc- tion. A Huffman entropic coder has been used for further compression. This lossless coder 3.5. Quantization and run-length coding increased compression ratio 1.7 times in aver- age for the test data. The best results have The modified wavelet coefficients have been achieved by the Huffman coder using been stored with wordlength of 7 and 8 bits. no predictor. For lower values, the error of compression dramatically increased having a negative ef- fect on compression efficiency. Zero-valued 4. Results and discussion coefficients are stored by a run-length coder. Relatively long sequences of coefficients re- The proposed wavelet-based compression lated to parts of the original signal between algorithm has been tested for the following I. Prouaznik, J. Kozumplik /International Journal of Medical Informatics 45 (1997) 111~ I28 125 ongmalslgnalw006.1 300 250 200 7 E 150 f > 100 50 0 02 04 06 08 1 12 14 16 18 2 t [set] decompresseds~gnal,threshold30,b~ts8,PRD=1512,CR=1283,wavelets b135 I I I I I I I I I I 250 I I I / t I I I I I 0 02 04 06 08 1 12 14 16 1.8 2 t [set] 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, quantiza- tion of wavelet coefficients to 8 bits. Resulted PRD = 15.12%, CR = 12.8. parameters: Meyer’s and biorthogonal 3/5 and Table 2. generating wavelets, thresholding in two and One should note PRD and CR in both three highest frequency bands with threshold tables are average results. Due to distortion of 30 pV, quantization with 7 and 8 bit of some records by powerline noise and/or resolution. Compression efficiency has been myopotentials, PRD can significantly arise as evaluated by compression ratio (CR), bit rate the noise is suppressed by thresholding. The (BR), and average code length (ACL), (see remaining noise shorten sequences of zero- Appendix A). Error of reconstruction has valued samples in particular bands and CR is been evaluated by percent root-mean-square then typically low. This is well visible in Figs. difference (PRD), (see Appendix B). For ex- 7 and 8 in which distribution of resulting CR ample, one can choose biorthogonal 3/5 and PRD for all tested signals is shown. The wavelets, with 8 bit resolution, and threshold- results for the compression algorithm using ing in three bands. Then, resulted compres- Meyer’s wavelets are in Fig. 7, for the al- sion ratio CR for tests on all 1000 signals gorithm using biorthogonal wavelets 3/5 in CR = 14.0 (BR = 429 bit/s, ACL = 0.86 bit/ Fig. 8. The distribution of CR and PRD for sample) when reconstruction error PRD = a set of test signals with removed 60 Hz 4.85%. Results of other tests are in Table 1 powerline noise is in Fig. 9. 126 I. Pmvaznik, J. Kozumplik / Imternational Journal of Medical Informatics 45 (I 997) I1 1- 128 ongmal signal wUO6f I decompressed signal, threshoM30. b&8. PRD=B 849, CR=16 55. wavelets b135 0 02 04 0.6 08 1 1.2 1.4 16 18 2 t [set] 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 algorithm using biorthogonal wavelets 3/5. signal No. 6, lead V2, from CSE library and Original signal No. 6, lead I, filtered by a the related reconstructed signal for the com- notch filter at 60 and 180 Hz to suppress pression algorithm using Meyer’s wavelets powerline noise and related reconstructed sig- (Fig. 10) and for the algorithm using nal are shown in Fig. 16. The compression biorthogonal wavelets 3/5 (Fig. 11) are algorithm used biorthogonal wavelets 3/5, shown. Fig. 12 and Fig. 13 corresponding to threshold 30 ,uV in three highest frequency Figs. 10 and 11, respectively, show three bands, quantization of wavelet coefficients to highest frequency bands of the decomposed 8 bits. original signal where horizontal lines repre- Thresholding in three highest frequency sent the used thresholds. bands can sporadically disturb shape of P- An example of the original signal distorted waves when their spectral contents exceeds 30 by powerline noise (CSE signal No. 6, lead I) Hz. In order to avoid this shortcoming, we and the related reconstructed signal is in Fig. propose thresholding in two highest fre- 14 and Fig. 15. Fig. 14 represents coding by quency bands resulted in lower CR (see Table the compression algorithm using Meyer’s 2). Another solution consists of subdividing wavelets. Fig. 15 represents coding by the the third frequency band (31.25-62.5 Hz) I. Provaznik, J. Kozumplik /International Journal of Medical Informaties 45 (1997) I II- 128 127 into two bands of same bandwidth and Acknowledgements thresholding in a higher band (46.875-62.5 Hz). Then, CR would vary between values This project was partly supported by the from Tables 1 and 2. grant No. 102/94/0374 of the Czech State For modification of wavelet coefficients, Grant Agency. we have proposed soft thresholding using an empirically derived threshold. It should be noted that the use of optimum thresholding Appendix A. Compression efficiency where threshold value is derived from noise characteristics can increase compression effi- Compression ratio (CR) is defined by the ciency. Thus, the algorithm using such formula: thresholding can minimize distortion of QRS CR = amount of original data [bit] amount of compressed data [bit] L-1 complexes. More sophisticated algorithms could employ a QRS complex detector and (23) modify the wavelet coefficients between QRS complexes. For the best results, we have pro- Bit rate (BR) is defined by the formula: posed to remove powerline noise before com- pression. Further PRD reduction could be BR reached by using different quantization in = amount of compressed data [bit] @its - ‘1 particular frequency bands. Also, the results original signal duration [s] are influenced by the selection of the used (24) wavelets. Average code length (ACL) is defined by the formula: 5. Conclusions ACL = amount of compressed data [bit] In this paper, we propose to use a lossy amount of original data [sample] compression algorithm for ECG data pro- [bit.sample - ‘1 (25) cessing. The algorithm is based on orthogo- nal or biorthogonal wavelet transforms. Unimportant wavelet coefficients are set to Appendix B. Error of compression zero by soft thresholding (30 pV) in two or three highest frequency bands. Thresholding The error of reconstruction is characterized and quantization to 7 or 8 bits combined by percent root-mean-square difference with Huffman coding leads to average code (PRD) by the formula: length from 0.65 to 1.3 bit per sample with corresponding relative root-mean-square er- ccx<4 - x’W2 ror (PRD) between 5.43 and 3.00%. This PRD= in 100 [%] (26) compression is about 3-5 times more effi- ,j ox’ cient than direct lossless entropic coding. The v n proposed method was tested on the CSE where x(n) is the original signal, x’(n) is the ECG signal library. reconstructed signal. 128 I. Prouaznik. J. Kommplik / International Journal of’ Medical InJbrmatics 45 (1997) 111~ 128 References [9] G. Strang, T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge Press, 1996. [l] U.E. Ruttimann, H.V. Pipberger, Compression of [lo] 0. Bertrand, J. Bohorquez, J. Pernier, Time- the ECG by prediction or interpolation and en- frequency digital filtering based on an inver- tropy encoding, IEEE Trans. BME 26 (11) (1979) tible wavelet transform: An application to evoked 613-623. potentials, IEEE Trans. BME 41 (1) (1994) 77- [2] S.C. Tai, An extensive Markov system for ECG 88. exact coding, IEEE Trans. BME 42 (2) (1995) [ll] J. Kozumplik, I. Provaznik, ECG Data Compres- 230-232. sion by Nonuniform Sampling. Proc. of 4th Int. [3] C.S. Tai, Six-band sub-band coder on ECG Conf. SYMBIOSIS’95, Gliwice, Poland, 1995, pp. waveforms, Med. Biol. Eng. Comp. 30 (1992) 98- 100. 187-192. [12] J.A. Crowe, N.M. Gibson, M.S. Woolfson, M.G. [4] M.C. Aydin, A.E. Cetin, H. Koymen, ECG data Somekh, Wavelet transform as a potential tool compression by sub-band coding, Electron. Lett. for ECG analysis and compression, J. Biomed. 27 (1991) 359-360. Eng. 14 (5) (1992) 2699272. [5] A.E. Cetin, H. Koymen, M.C. Aydin, Multichan- [13] C. Li, C. Zheng, C. Tai, Detection of ECG char- nel ECG data compression by multirate signal acteristic points using wavelet transform, IEEE processing and transform domain coding tech- Trans. BME 42 (1) (1995) 21-28. niques, IEEE Trans. BME 40 (5) (1993) 4955499. [14] L. Sendhaji, G. Carrault, J.J. Belanger, G. Pas- [6] G. Nave, A. Cohen, ECG compression using sarielo, Comparing wavelet transforms for recog- long-term prediction, IEEE Trans. BME 40 (9) nizing cardiac patterns, IEEE Eng. Med. Biol. 13 (1993) 877-885. (2) (1995) 1677173. [7] I. Daubechies, Ten Lectures on Wavelet. Society [15] 0. Meste, H. Rix, P. Caminal, N. Thakor, Ven- for Industrial and Applied Mathematics, tricular late potentials characterization in time- Philadelphia, 1992. frequency domain by means of a wavelet [8] P.P. Vaidyanathan, Multirate Systems and Filter transform, IEEE Trans. BME 41 (7) (1994) 6255 Banksm, Prentice-Hall, 1993. 633.