Open access

Discrete Wavelet Transform Algorithms for Multi-Scale Analysis of Biomedical Signals

Written By

Juuso T. Olkkonen and Hannu Olkkonen

Published: 12 September 2011

DOI: 10.5772/38297

From the Edited Volume

Discrete Wavelet Transforms - Biomedical Applications

Edited by Hannu Olkkonen

Chapter metrics overview

1,975 Chapter Downloads

View Full Metrics

1. Introduction

The discrete wavelet transform (DWT) algorithms have a firm position in multi-scale processing of biomedical signals, such as EMG and EEG. The DWT algorithms were initially based on the compactly supported conjugate quadrature filters (CQFs) (Smith & Barnwell, 1986; Daubechies, 1988). However, a drawback in CQFs is due to the nonlinear phase effects such as spatial dislocations in multi-scale analysis. This is avoided in biorthogonal discrete wavelet transform (BDWT) algorithms, where the scaling and wavelet filters are symmetric and linear phase. The biorthogonal filters are usually constructed by a ladder-type network called lifting scheme (Sweldens, 1988; ITU-T, 2000). Efficient lifting BDWT structures have been developed for microprocessor and VLSI environment (Olkkonen et al. 2005; Olkkonen & Olkkonen, 2008). Only integer register shifts and summations are needed for implementation of the analysis and synthesis filters.

A severe obstacle in multi-scale DWT analysis is the dependence of the total energy of the wavelet coefficients in different scales on the fractional shifts of the analysed signal. If we have a discrete-time signal x [ n ] and the corresponding time shifted signal x [ n τ ] , where τ [ 0,1 ] , there occurs a notable difference in the energy of the wavelet coefficients as a function of the time shift. Kingsbury (2001) proposed a nearly shift invariant method, where the real and imaginary parts of the complex wavelet coefficients are approximately a Hilbert transform pair. The energy (absolute value) of the wavelet coefficients equals the envelope, which provides smoothness and approximate shift-invariance. Selesnick (2002) observed that using two parallel CQF banks, which are constructed so that the impulse responses of the scaling filters have half-sample delayed versions of each other: h [ n ] and h [ n 0.5 ] , the corresponding wavelet bases are a Hilbert transform pair. Selesnick (2002) proposed a spectral factorization method based on the half delay all-pass Thiran filters for design of the scaling filters. However, the scaling filters do not owe coefficient symmetry and the nonlinearity interferes with the spatial timing in different scales and prevents accurate statistical correlations.

In this book chapter we review the shift invariant DWT algorithms for multi-scale analysis of biomedical signals. We describe a dual-tree DWT, where two parallel CQF wavelet sequences form a Hilbert pair, which warrants the shift invariance. Next we review the construction of the shift invariant BDWT, which is based on the novel design of the Hilbert transform filter. Finally, we describe the FFT based computation of the analytic signal and the implementation of the shift invariant quadrature mirror filter (QMF) bank.

Advertisement

2. Shift invariant CQF bank

In the following we describe a shift invariant DWT algorithm based on two parallel real- valued CQF banks. The conventional CQF DWT bank consists of the H 0 ( z ) and H 1 ( z ) analysis filters and G 0 ( z ) and G 1 ( z ) synthesis filters for N odd (Fig. 1)

H 0 ( z ) = ( 1 + z 1 ) K P ( z ) H 1 ( z ) = z N H 0 ( z 1 ) G 0 ( z ) = H 1 ( z ) G 1 ( z ) = H 0 ( z ) E1

where P ( z ) is a polynomial in z 1 . The scaling filter H 0 ( z ) has the Kth order zero at ω = π . The wavelet filter H 1 ( z ) has the Kth order zero at ω = 0 , correspondingly. The filters are related via the perfect reconstruction (PR) condition

[ H 0 ( z ) H 1 ( x ) H 0 ( z ) H 1 ( z ) ] [ G 0 ( z ) G 1 ( z ) ] = [ 2 z N 0 ] E2

Figure 1.

The analysis and synthesis parts of the real-valued CQF DWT bank.

Let us denote the frequency response of the z-transform filter as

H ( z ) = n h n z n H ( ω ) = n h n e j ω n E3

Then we obtain the relations

H ( z ) H ( ω π ) H ( z 1 ) H ( ω π ) E4

where * denotes complex conjugation. The tree structured implementation of the two parallel real-valued CQF filter banks is described in Fig. 2. In M-stage CQF tree the frequency response of the wavelet sequence is

W M ( ω ) = H 1 ( ω 2 ) k = 2 M H 0 ( ω 2 k ) E5

Figure 2.

The implementation of two parallel real-valued CQF banks, which yields the wavelet sequences w 1 [ n ] , w 2 [ n ] and w ¯ 1 [ n ] , w ¯ 2 [ n ]

Next we construct a phase shifted parallel CQF filter bank consisting of the scaling filter H ¯ 0 ( z ) and the wavelet filter H ¯ 1 ( z ) (Fig. 2). Let us suppose that the scaling filters in parallel CQF trees are related as

H ¯ 0 ( ω ) = e j ϕ ( ω ) H 0 ( ω ) E6

where ϕ ( ω ) is a 2 π periodic phase function. Correspondingly, the CQF filters are related as

H 1 ( ω ) = e j ω N H 0 * ( ω π ) E7

We have

H ¯ 1 ( ω ) = e j ω N H ¯ 0 * ( ω π ) = e j ω N e j ϕ ( ω π ) H 0 * ( ω π ) = e j ϕ ( ω π ) H 1 ( ω ) E8

We may note that the phase shifted CQF bank (6,8) obeys the PR condition (2). The frequency response of the M-stage CQF wavelet sequence is

W ¯ M ( ω ) = H ¯ 1 ( ω 2 ) k = 2 M H ¯ 0 ( ω 2 k ) = e j ϕ ( ω 2 π ) H 1 ( ω 2 ) k = 2 M e j ϕ ( ω 2 k ) H 0 ( ω 2 k ) = e j ϕ ( ω 2 π ) e j k = 2 M ϕ ( ω 2 k ) H 1 ( ω 2 ) k = 2 M H 0 ( ω 2 k ) = e j θ W M ( ω ) E9

where the phase function

θ = ϕ ( ω 2 π ) k = 2 M ϕ ( ω 2 k ) E10

By selecting the phase function ϕ ( ω ) in (6) as

ϕ ( ω ) = ω 2 E11

the scaling filters (6) are half-sample delayed versions of each other. By inserting (11) in (10) we have

θ = ω 2 π 2 ω k = 2 M 1 2 k + 1 = π 2 + ω 2 M + 1 E12

The wavelet sequences (5,9) yielded by the CQF bank (1) and the phase shifted CQF bank (6,8) can be interpreted as real and imaginary parts of the complex wavelet sequence

W M C ( ω ) = W M ( ω ) + j W ¯ M ( ω ) E13

The requirement for the shift-invariance comes from

W ¯ M ( ω ) = { ψ M ( ω ) } E14

where denotes the Hilbert transform. The frequency response of the Hilbert transform operator is defined as ( ω ) = j sgn ( ω ) , where sgn ( ω ) = 1 for ω 0 and sgn ( ω ) = 0 for ω 0 . In this work we apply the Hilbert transform operator in the form

( ω ) = e j π / 2 sgn ( ω ) E15

The result (12) indicates that if the scaling filters are the half-sample delayed versions of each other, the resulting wavelet sequences are not precisely Hilbert transform pairs. There occurs a phase error term ω / 2 M + 1 , which depends both in frequency and the stage M of the wavelet sequence. However, the error term can transferred in front of the CQF tree by using the equivalence described in Fig. 3. Then the error term reduces to ω / 2 and the phase error term can be simply eliminated by prefiltering the analyzed signal by the half-sample delay operator, which has the frequency response D ( ω ) = e j ω / 2 . The total phase function is then

θ ( ω ) = D ( ω ) π / 2 + ω / 2 = π / 2 , which implies that the M-stage CQF wavelet sequence and the phase error corrected sequence are a Hilbert transform pair.

Figure 3.

The two equivalents for moving the phase function in front of the phase shifted CQF tree.

The two parallel BDWT trees can be considered to form a complex wavelet sequence by defining the Hilbert transform operator

a ( z ) = 1 + j ( z ) E16

By filtering the real-valued signal x [ n ] by the Hilbert transform operator results in an analytic signal

x a [ n ] = x [ n ] + j { x [ n ] } E17

whose magnitude response is zero at negative side of the frequency spectrum

X a ( ω ) = { 2 X ( ω ) 0 ω π 0 π ω 0 E18

Let us consider the complex wavelet sequence at the first stage (Fig. 6).The wavelet sequence is obtained by decimation of the high-pass filtered analytic signal

W 1 ( ω ) = [ X a ( ω ) H 1 ( ω ) ] 2 = W a ( ω ) 2 = 1 2 X a ( ω 2 ) H 1 ( ω 2 ) E19

The frequency spectrum of the undecimated wavelet sequence W a ( ω ) contains frequency components only in the range 0 ω π , but the frequency spectrum of the decimated analytic signal has the frequency band 0 ω 2 π . Hence, the decimation does not produce overlapping and leakage (aliasing) to the negative frequency range.

A key feature of the dual-tree wavelet transform is the shift invariance of the decimated analytic wavelet coefficients. The frequency spectrum of the decimated wavelet sequence of the fractionally delayed signal x [ n τ ] is 0.5 e j ω τ / 2 W a ( ω / 2 ) . The energy of the decimated wavelet coefficients is 0.5 | W ( ω / 2 ) | , which does not depend on the fractional delay.

Advertisement

2. Shift invariant BDWT filter bank

The two-channel BDWT filter bank is of the general form

H 0 ( z ) = ( 1 + z 1 ) L Q ( z ) H 1 ( z ) = ( 1 z 1 ) M R ( z ) G 0 ( z ) = H 1 ( z ) G 1 ( z ) = H 0 ( z ) E20

where the scaling filter H 0 ( z ) has the Lth order zero at ω = π . The wavelet filter H 1 ( z ) has the Mth order zero at ω = 0 , correspondingly. Q ( z ) and R ( z ) are polynomials in z 1 . The low-pass and high-pass reconstruction filters G 0 ( z ) and G 1 ( z ) are defined as in the CQF bank. For two-channel BDWT filter bank the PR relation is

[ H 0 ( z ) H 1 ( x ) H 0 ( z ) H 1 ( z ) ] [ G 0 ( z ) G 1 ( z ) ] = [ 2 z N 0 ] E21

An essential result is related to the modification of the BDWT bank (Olkkonen & Olkkonen, 2007a).

Lemma 1: If the scaling filter H 0 ( z ) , the wavelet filter H 1 ( z ) and the reconstruction filters G 0 ( z ) and G 1 ( z ) in BDWT filter bank (20) have a perfect reconstruction property (21), the following modified BDWT filter bank obeys the PR relation

H ¯ 0 ( z ) = F ( z ) H 0 ( z ) H ¯ 1 ( z ) = F 1 ( z ) H 1 ( z ) G ¯ 0 ( z ) = F 1 ( z ) G 0 ( z ) G ¯ 1 ( z ) = F ( z ) G 1 ( z ) E22

where F ( z ) is any polynomial in z 1 . Proof is yielded by direct insertion (22) to PR condition (21).

In the following we apply Lemma 1 for constructing the shift invariant BDWT filter bank. We describe a novel method for constructing the Hilbert transform filter ( z ) based on the half-sample delay filter D ( z ) = z 0.5 . The classical approach for design of the half-sample delay filter D ( z ) is based on the Thiran all-pass interpolator

D ( z ) = z 0.5 = k = 1 p c k + z 1 1 + c k z 1 = z N C ( z 1 ) C ( z ) = c N + c N 1 + + z N 1 + c 1 z 1 + + c N z N E23

where the c k coefficients are optimized so that the frequency response follows approximately D ( ω ) = e j ω / 2 . Correspondingly, the quadrature mirror filter D ( z ) has the frequency response

D ( ω π ) = e j ( ω π ) / 2 E24

The Hilbert transform filter is then obtained as

( ω ) = D ( ω ) D ( ω π ) = e j ω / 2 e j ( ω π ) / 2 = e j π / 2 ( z ) = D ( z ) D ( z ) E25

The Hilbert transform filter is inserted in the BDWT bank using the result of Lemma 1 (22). The modified prototype BDWT filter bank is

H ¯ 0 ( z ) = ( z ) H 0 ( z ) H ¯ 1 ( z ) = 1 ( z ) H 1 ( z ) G ¯ 0 ( z ) = 1 ( z ) G 0 ( z ) G ¯ 1 ( z ) = ( z ) G 1 ( z ) E26

A highly simplified BDWT filter bank can be obtained by noting that in (25) 1 ( z ) = ( z ) . We have

H ¯ 0 ( z ) = ( z ) H 0 ( z ) H ¯ 1 ( z ) = ( z ) H 1 ( z ) G ¯ 0 ( z ) = ( z ) G 0 ( z ) G ¯ 1 ( z ) = ( z ) G 1 ( z ) E27

The modified BDWT filter bank (27) can be realized by the Hilbert transform filter ( z ) , which works as a prefilter for the analysed signal (Fig. 4). The Hilbert transform filter ( z ) works as a postfilter in the reconstruction stage, respectively.

Figure 4.

The realization of the Hilbert transform filter.

An integer-valued Hilbert transform filter can be constructed by the B-spline transform (see details Olkkonen & Olkkonen, 2011b). The frequency response of the Hilbert transform filter shows a maximally flat magnitude spectrum. The phase spectrum corresponds to an ideal Hilbert transformer (15).

The Hilbert transform filter in Fig. 4 can be replaced by the Hilbert transform operator (16), which yields an analytic signal. This avoids the need for two parallel filter banks. In the following we describe a FFT based method for computation of the analytic signal and the implementation of the shift invariant quadrature mirror filter (QMF) bank.

Advertisement

3. FFT based computation of analytic signal

The fast Fourier transform of the signal x [ n ] , n = 0, 1, 2, …,N-1 is of the form

F F T N { x [ n ] } = Y k = n = 0 N 1 x [ n ] W N n k k = 0,..., N 1 E28

where W N = e 2 π j / N . The FFT coefficients Y k (k=N/2,…,N-1) represent the values in the negative frequency band ( π ω 0 ) . By zeroing those coefficients, the inverse fast Fourier transform (IFFT) yields an analytic signal. A more accurate result is obtained by weighting the FFT coefficients by a window

w k = { 2 k = 1,2,..., N / 2 1 0 k = N / 2 + 1,..., N 1 1 k = 0 a n d N / 2 E29

The analytic signal is then computed using the inverse FFT transform

x a [ n ] = 1 N k = 0 N 1 w k Y k W N n k = I F F T N { w k Y k } E30

The weighting sequence in (29) can be eliminated by writing

x a [ n ] = 1 N / 2 k = 0 N / 2 1 Y k W N n k Y 0 N + ( 1 ) n Y N / 2 N E31

Now, for even n we have

x a [ 2 n ] = 1 N / 2 k = 0 N / 2 1 Y k W N / 2 n k m e a n { x [ 2 n + 1 ] } = I F F T N / 2 { Y k } m e a n { x [ 2 n + 1 ] } E32

and for odd n

x a [ 2 n + 1 ] = 1 N / 2 k = 0 N / 2 1 W N k Y k W N / 2 n k m e a n [ x [ 2 n ] } = I F F T N / 2 { W N k Y k } m e a n { x [ 2 n ] } E33

For zero mean signal m e a n { x [ n ] } = 0 , which yields m e a n { x [ 2 n + 1 ] } = m e a n { x [ 2 n ] } . If the even points of the analytic signal are known, the FFT coefficients are solved from (32)

Y k = F F T N / 2 { x a [ 2 n ] m e a n { x [ 2 n ] } } k = 0,..., N / 2 1 E34

The odd points of the analytic signal are then computed from (33). We call this as the reconstruction property of the zero mean analytic signal. In the following we present a novel shift invariant QMF bank, which utilizes the reconstruction property of the analytic signal.

Advertisement

4. Shift invariant QMF bank

In QMF bank the scaling and wavelet filters obey the relation H 1 ( z ) = H 0 ( z ) , i.e. their frequency response is symmetric with respect to ω = π / 2 . In this work we define the scaling and wavelet filters as half band QMFs

H 0 ( z ) = 1 2 + z 1 A ( z 2 ) H 1 ( z ) = 1 2 z 1 A ( z 2 ) E35

The shift invariant tree structured QMF DWT is described in Fig. 5. The FFT based Hilbert transform operator a ( z ) produces an analytic signal, which is fed to the scaling H 0 ( z ) and wavelet H 1 ( z ) filters and decimated. If the original zero mean signal is x [ n ] , the decimated scaling and wavelet coefficients s [ n ] and w [ n ] are obtained from

s [ n ] = { h 0 [ n ] x a [ n ] } 2 w [ n ] = { h 1 [ n ] x a [ n ] } 2 E36

where denotes convolution. From (35) we have

H 0 ( z ) + H 1 ( z ) = 1 h 0 [ n ] + h 1 [ n ] = δ n E37

The reconstruction consists of the summation of the decimated signals. We obtain

s [ n ] + w [ n ] = { ( h 0 [ n ] + h 1 [ n ] ) x a [ n ] } 2 = { δ [ n ] x a [ n ] } 2 = x a [ 2 n ] E38

i.e. the summation of the decimated signals produces the even points x a [ 2 n ] of the analytic signal. The odd points x a [ 2 n + 1 ] of the analytic signal are then reconstructed from the even points x a [ 2 n ] via our results (32)-(34). The original signal is obtained from x [ n ] = r e a l ( x a [ n ] ) .

Figure 5.

The shift invariant tree structured QMF DWT.

Advertisement

5. Conclusion

The dual-tree DWT algorithms have appeared to outperform the real-valued DWTs in several applications such as denoising, texture analysis, speech recognition, processing of seismic signals and multiscale-analysis of neuroelectric signal analysis (Olkkonen et al. 2006; Olkkonen et al. 2007b, Olkkonen & Olkkonen, 2010, Olkkonen & Olkkonen 2011a).

Selesnick (2002) noted that a half-sample time-shift between the scaling filters in parallel CQF banks yields a nearly shift invariant DWT, where the wavelet bases form a Hilbert transform pair. However, the multi-scale analyses of neuroelectric signals have revealed that the first stages of wavelet sequences are quite poorly shift invariant. We reanalysed the condition and observed a phase-error term ω / 2 M + 1 (12) compared with the ideal phase response θ ( ω ) = π / 2 . The phase error attains s highest value at high frequency range and small stage M of the wavelet sequence. Fortunately the phase error term can be cancelled by adding a half-delay prefilter in front of the CQF chain. For this purpose the half-delay filter constructed by the B-spline transform (Olkkonen & Olkkonen, 2011b) is well suited. In addition, there exists many other design methods for half-delay filters (see e.g. Laakso et al. 1996; Johansson & Lowenborg, 2002; Pei & Tseng, 2003; Pei & Wang, 2004; Tseng, 2006).

In this book chapter we described a novel shift invariant dual-tree BDWT (27) based on Lemma 1 (22) and the Hilbert transform filter (25). In many respects the shift invariant BDWT bank outperforms the previous nearly shift invariant DWT approaches. The Hilbert transform filter assisted BDWT yields precisely shift invariant wavelet sequences, which permits the statistical analyses between scales in multi-scale analyses of biomedical signals such as EMG and EEG.

The Hilbert transform filter in Fig. 4 can be replaced by the Hilbert transform operator (16), which yields an analytic signal. This avoids the need for two parallel filter banks. In this work we described a FFT based method for computation of the analytic signal and the implementation of the shift invariant QMF bank. As a clear advantage of the half-band QMF structure is that the frequency responses of the scaling and wavelet filters are mirror symmetric with respect to ω = π / 2 . Hence, they split the energy of the signal precisely to the low-pass and high-pass fractions. The energy preservation property is of utmost importance in automated statistical signal processing of the multi-scale signals. In tree structured multi-scale analysis the linear phase of the QMFs is advantageous since the timing information of the wavelet coefficients in different scales is preserved. Without an exact timing of the subscale signals the statistical comparison of the wavelet coefficients in different scales is not relevant and may lead to misleading results. For example in EEG signal the neuroelectric discharge contains fast repetitive transients with related timing and overlapping waveforms. In multi-scale analysis different components can be separated due to their different timing and scale related intensification.

Advertisement

Acknowledgments

This work was supported by the National Technology Agency of Finland (TEKES).

References

  1. 1. Daubechies I. 1988Orthonormal bases of compactly supported wavelets. Commmun. Pure Appl. Math., 41
  2. 2. ITU-T(2000), , 800 -ISO DCD15444-1: JPEG2000 Image Coding System. Recommend. T.International Organization for Standardization, ISO/IEC JTC! SC29/WG1.
  3. 3. Johansson H. Lowenborg P. 2002Reconstruction of nonuniformy sampled bandlimited signals by means of digital fractional delay filters, IEEE Trans. Signal Process., 50 11 2757 2767
  4. 4. Kingsbury N. G. 2001Complex wavelets for shift invariant analysis and filtering of signals. J. Appl. Comput. Harmonic Analysis. 10
  5. 5. Laakso T. Valimaki V. Karjalainen M. Laine U. K. 1996Splitting the unit delay. Toolsfor fractional delay filter design, IEEE Signal Processing Magazine, 30 80
  6. 6. Olkkonen H. Pesola P. Olkkonen J. T. 2005Efficient lifting wavelet transform for micro-processor and VLSI applications. IEEE Signal Process. Lett. 12 2
  7. 7. Olkkonen H. Pesola P. Olkkonen J. T. Zhou H. 2006Hilbert transform assisted complex wavelet transform for neuroelectric signal analysis. J. Neuroscience Meth. 151
  8. 8. Olkkonen H. Olkkonen J. T. 2007aHalf-delay B-spline filter for construction of shift-invariant wavelet transform. IEEE Trans. Circuits and Systems II. 54 7
  9. 9. Olkkonen H. Olkkonen J. T. Pesola P. 2007bFFT-based computation of shift invariant analytic wavelet transform. IEEE Signal Process. Lett. 14 3
  10. 10. Olkkonen H. Olkkonen J. T. 2008Simplified biorthogonal discrete wavelet transform for VLSI architecture design. Signal, Image and Video Process. 2
  11. 11. Olkkonen H. Olkkonen J. T. 2010Shift-invariant B-spline wavelet transform for multi-scale analysis of neuroelectric signals. IET Signal Process. 4 6
  12. 12. Olkkonen J. T. Olkkonen H. 2011aShift invariant biorthogonal discrete wavelet transform for EEG signal analysis. Book chapter in: Discrete Wavelet Transforms- Theory and Applications, Edited by Juuso T. Olkkonen, Intech, 169 178
  13. 13. Olkkonen J. T. Olkkonen H. 2011bComplex Hilbert transform filter. J. Signal and Information Process. 2
  14. 14. Pei T. S. C. Tseng C. C. 2003An efficient design of a variable fractional delay filter using a first-order differentiator, IEEE Signal Processing Letters, 10 10 307 310
  15. 15. Pei S. C. &and Wang. P. H. 2004Closed-form design of all-pass fractional delay, IEEE Signal Processing Letters, 11 10 788 791
  16. 16. Selesnick I. W. 2002The design of approximate Hilbert transform pairs of wavelet bases. IEEE Trans. Signal Process. 50 5
  17. 17. Smith M. J. T. Barnwell T. P. 1986Exaxt reconstruction for tree-structured subband coders. IEEE Trans. Acoust. Speech Signal Process. 34
  18. 18. Sweldens W. 1988The lifting scheme: A construction of second generation wavelets. SIAM J. Math. Anal. 29
  19. 19. Tseng C. C. . 2006Digital integrator design using Simpson rule and fractional delay filter, IEEE Proc. Vision, Image and Signal Process., 153 1 79 85

Written By

Juuso T. Olkkonen and Hannu Olkkonen

Published: 12 September 2011