Wavelet transform in electrocardiography-data compression by bdr85067

VIEWS: 77 PAGES: 18

									                                                                                                                                                       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.

								
To top