A Low-power Analogue SC-CMOS Filter suitable to implement the Wavelet Algorithm to analyse ECG signals in Pacemakers

Contents

Acknowledgments

Introduction

Chapter I Evaluation of the filter expression in s-domain and z-domain
1.1) Laplace transform and Padè approximation
1.2) Impulse invariance technique

Chapter II Filter optimization – dynamic range analysis
2.1) State-Space realization
2.2) State-Space optimization
2.3) Hessenberg and Schur decompositions
2.4) Optimal capacitance distribution

Chapter III Comparison of discrete and continuous time approaches
3.1) State-Space realization in continuous time

Chapter IV Circuital realization of the filter
4.1) Synthesis of the state-space filter
4.2) Filter settings
4.3) First approach to the circuit
4.4) Impulse transient response of the circuit with ideal amplifiers
4.5) Step transient response of the circuit with ideal amplifiers
4.6) Sine transient response of the circuit with ideal amplifiers
4.7) Clock feedthrough

Chapter V Two-stage CMOS operational amplifier
5.1) Two-stage CMOS operational amplifier topology
5.2) Analysis of the amplifier
5.3) Frequency compensation
5.4) Integrator bandwidth
5.5) Considerations on the amplifier
5.6) Calculus of the transconductance by noise-optimisation
5.7) Calculus of the biasing current and design of the amplifier

Chapter VI Final schematic of the wavelet filter
6.1) Transmission gate and dummy gate
6.2) Impulse transient response
6.3) Step transient response
6.4) Sine transient response
6.5) Power consumption
6.6) Noise analysis

Chapter VII Considerations on the wavelet filter

References

Acknowledgments

To my parents, Carmela and Sebastiano Salvo, for their constant support, help and patience during all these years.

To my supervisors W. A. Serdijn and S. A. P. Haddad, for receiving me very well in Delft, for their constant help and important advising.

Introduction

The QRS complex represents ventricular depolarization and consists of three waveforms. The normal complex begins with a downward deflection known as the Q wave, followed by an upward deflection called the R wave. The next downward deflection will be the S wave. All ventricular complexes are known as QRS complexes even if every wave is not present in all complexes. The normal human QRS complex lasts about 0 . 04 to 0 . 11 seconds and its waveform can be seen in Fig. 1, which shows a typical external electrocardiogram (ECG).

illustration not visible in this excerpt

Fig. 1. Example of a typical external electrocardiogram.

The detection of the QRS complex is very important because it is related to different heart dysfunctions such as:

- Nonspecific intraventricular conduction delay
- Some cases of left anterior or posterior fascicular block
- Incomplete right or left bundle branch block
- Ectopic rhythms originating in the ventricles (e.g., ventricular tachycardia, pacemaker rhythm)
- Presence of necrotic heart tissue
- Ventricular hypertrophy

The wavelet transformation [1] is a very promising tool to characterize non-stationary signals such as the QRS complex because it gives good estimation of time and frequency localization. In fact, the analysis of the signal is performed at various resolutions and accomplished by decomposition into elementary functions that are well localized both in time and frequency domain. The idea is to introduce a so-called mother wavelet to analyse a given signal which will be represented in terms of a wavelet expansion using a linear combination of the coefficients of the wavelet functions derived from the mother wavelet by using different scales.

Let Ψ(t) be a real or complex valued function in L2(R), the Hilbert space of measurable, square integrable, one-dimensional functions [2].

The function Ψ(t) is said to be a wavelet if and only if its Fourier transform Ψ'(ω) satisfies:

illustration not visible in this excerpt

This admissibility condition implies that:

illustration not visible in this excerpt

which means that Ψ(t) is oscillatory and its average is zero.

Let

illustration not visible in this excerpt

be the dilation of Ψ(t) by a scale factor of a > 0. The factor in the above expression is used for energy normalization.

The wavelet transform of a function f(t) L2(R) at a scale a and position ζ is given by:

illustration not visible in this excerpt

where * denotes the complex conjugation.

This transformation satisfies energy conservation and the original signal f(t) can be reconstructed from the wavelet transform.

For smaller values of scale a and the position ζ the wavelet is contracted in the time domain and the wavelet transform gives information about the finer details of the signal. For larger values of a, the wavelet expands and its transform gives a global view of the signal.

If the scale parameter is the set of integral powers of [Abbildung in dieser Leseprobe nicht enthalten], the wavelet is called a dyadic wavelet. To cover the whole frequency domain, the Fourier transform of must satisfy the relation:

illustration not visible in this excerpt

The detection of the QRS complex is based on the modulus maxima of the wavelet transform defined as any point Wf (2j, τ0) such that

illustration not visible in this excerpt

The modulus maxima and the zero-crossing of the wavelet transform correspond to the sharp edges in the signal [3]. In the QRS complex there are two modulus maxima with opposite signs of Wf (2j, τ) with a zero-crossing between them; therefore a good method to detect the QRS complex is to apply thresholds to the wavelet transform of the ECG signal.

However, the advantageous properties derived from the wavelet transform also have a limitation due to the time-frequency representation which cannot do better than what follows from the Heisenberg principle:

illustration not visible in this excerpt

This relation means that it is not possible to achieve at the same time good performances in time and frequency. For this reason, a good choice for the wavelet is the first derivative Gaussian function, which has the minimum time-frequency product of 1/2. Moreover, the waveform of the first derivative Gaussian function resembles the shape of the QRS complex.

The filter can be implemented either in continuous time or in discrete time. Implementing linear conductances as resistors is not possible because of the large values which should be used in low power applications. In continuous time it is possible, for example, using a log-domain approach that bypasses the need for large resistances. Nevertheless, in this research work the switched-capacitor (SC) approach in discrete time is preferred because it allows for a large signal voltage swing at the integrator output and consequently a better dynamic range.

In order to implement the filter in the discrete time (in the z-domain), the impulse invariance technique is adopted to preserve the impulse response of the continuous time analysis by defining

illustration not visible in this excerpt

where T is the sample period. The analysis in the z-domain is very useful to easily implement the filter with the SC technique.

Starting from the mathematical expression of the Gaussian function, after the Laplace Transform, the Padè Approximation will be used to arrive at a rational form of the function suitable for a filter implementation. It will be shown that a SC realization is preferable instead of a direct implementation of the filter in Continuous Time in order to avoid high-value resistor elements. The optimization of the sensitivity, noise and dynamic range of the filter will be carried out following a State-Space methodology, as described in [4], both for Continuous and Discrete Time. The resulting filter is a 7th order complementary metal-oxide semiconductor (CMOS) SC filter designed implementing a two-stage operational amplifier. The proposed technique will be proven to be suitable for pacemaker applications where low power consumption is a necessary requirement.

Chapter I Evaluation of the filter expression in s-domain and z-domain

1.1) Laplace transform and Padè approximation

A Gaussian function has generally the following parametric form:

illustration not visible in this excerpt

Let . The previous expression can be rewritten as:

illustration not visible in this excerpt

To guarantee the causality of the filter in the evaluation of the Laplace transform, the Gaussian function has been delayed by choosing. In this way, as it is possible to see from Fig. 2, the Gaussian function lays on the right half of the plane.

illustration not visible in this excerpt

Fig. 2. Gaussian function g centred in t=3.

Its first derivative is (Fig. 3):

illustration not visible in this excerpt

Fig. 3. First derivative of the Gaussian function.

After applying the Laplace transform, it is necessary to apply the Padè approximation to arrive at a rational function which can lead us to the filter implementation [5].

The Padè approximation can be thought as a generalization of the Taylor polynomial. To be more precise, a Padè approximation of order (m, n) of a function f(x) at a point x0 is the rational function where p(x) is a polynomial of degree m while q(x) is a polynomial of degree n and the formal power series of f(x)q(x)-p(x) in the point x0 begins with the term xm+n+1.

A first attempt to implement the filter has been carried out by using the Chebyshev polynomials because this allows approximating the Gaussian function in an interval and not only in a specific point. Unfortunately, the resulting filter, c(s), is instable as proven by the poles calculation (performed in Maple software).

illustration not visible in this excerpt

The poles are:

illustration not visible in this excerpt

As the real parts of s4 and s5 are positive, the filter is not realizable.

Thus, it has been chosen to apply directly the Padè approximation.

Different orders of approximation have been tried at the point x0 = 0; increasing the order yields a better approximation but it also requires a higher filter order. In Figs. 4-7 it is possible to see the results for approximations of orders (3,6), (3,7), (3,8) and (3,9).

illustration not visible in this excerpt

Fig. 4. Padè approximation of order (3,6) of g'.

illustration not visible in this excerpt

Fig. 5. Padè approximation of order (3,7) of g'.

illustration not visible in this excerpt

Fig. 6. Padè approximation of order (3,8) of g'.

illustration not visible in this excerpt

Fig. 7. Padè approximation of order (3,9) of g'.

A good compromise seems to be the order (3,7), which is chosen for the filter. The mathematical expression is

illustration not visible in this excerpt

The poles are:

illustration not visible in this excerpt

For all the poles lay on the left side of the s-plane, the filter is stable. Moreover, starting with a stable filter in the s-domain will give a stable filter in the z-domain.

1.2) Impulse invariance technique

(1.1) gives the transfer function of an IIR filter. The impulse invariance technique is applied to transform the function into the z-domain.

The impulse response, h(t), of an analogue filter of order m is defined as the inverse Laplace transform of its system function Ha(s) [6], given as

illustration not visible in this excerpt

where the (-si) are the poles of the analogue transfer function while Ai are constants. Thus,

illustration not visible in this excerpt

where L-1 denotes the inverse Laplace transform.

If we require that

illustration not visible in this excerpt

then

illustration not visible in this excerpt

Taking the z-transform of this last equation, it follows that

illustration not visible in this excerpt

or

illustration not visible in this excerpt

It must be noted that the complex-plane mapping produced by the z-transform, although it gives the same time response, provides a different frequency response.

The frequency response of a sampled data filter function is equal to the frequency response of a continuous function plus contributions of the response displaced by multiples of 2π/T. This addition is called "aliasing".

Aliasing can be reduced if we sample faster than twice the bandwidth of the function (Nyquist's theorem):

illustration not visible in this excerpt

where fs is the sampling frequency and B the bandwidth of the function.

To respect Nyquist's theorem and to have a good resolution, a normalized sampling frequency of 3 Hz is chosen. A higher value of the sampling frequency is not advisable because it renders the poles to be too close to the unit circle, which causes the instability of the filter.

For the calculus of the impulse invariance of (1.1), Matlab® has been used. We have that

illustration not visible in this excerpt

where

and

illustration not visible in this excerpt

Fig. 8 shows the position of the poles. It is possible to see that they are inside the unit circle to confirm the stability of the filter.

illustration not visible in this excerpt

Fig. 8. Zero-pole plot of g3,7. The filter is stable as the poles sufficiently lay inside the unit circle.

Using Simulink®, the discrete step response of the filter can be evaluated. The result is shown in Fig. 9.

illustration not visible in this excerpt

Fig. 9. Step response of g3,7(z). The step amplitude is 1.

If we compare Fig. 9 to Fig. 5, we can see the validity of the invariance impulse method.

In next chapter, the (1.2) will be optimized for power and noise applying a state-space realization.

Chapter II Filter optimization – dynamic range analysis

2.1) State-Space realization

Dynamic range is formally defined as the ratio of the maximum to the minimum signal magnitudes that a system is able to process individually, not necessarily at the same time [4]. It is different from the instantaneous signal-to-noise ratio, which refers to the ratio of the simultaneous measure of signal and noise at a given instant of time. For linear systems both definitions coincide. However, in nonlinear systems, cross-products of signal and noise appear at the output, making the instantaneous noise magnitude dependent on the correspondent signal magnitude.

A system's dynamic range is essentially determined by:

- system's maximum signal representation capacity
- internally generated noise

As they are independent, it is desirable to optimize them orthogonally. Following this idea, the controllability and observability gramians are important in the orthogonalization, which can be achieved by a state-space realization.

Given the discrete function G(z), with indexes M and N, , written as

illustration not visible in this excerpt

a state-space representation of G(z) is defined as

illustration not visible in this excerpt

where u is the input, x the state vector and y the output. The matrices A, B, C, D are defined as

illustration not visible in this excerpt

This form is called canonical form.

The transfer function can be expressed in terms of the state matrices as

illustration not visible in this excerpt

As it is well known, the limited word-length of filter coefficients yields deviations from the required transfer function. To reduce these deviations, which can be very relevant because of the high order of the filter, the long E number representation has been used in Matlab®.

The controllability gramian K is defined as

illustration not visible in this excerpt

and is referred to the state power of the signal.

The observability gramian W is defined as

illustration not visible in this excerpt

and is referred to the system's output noise.

Assuming that all the integrators of the filter have the same output amplitude AM, which should not be exceeded for more than a limited fraction of time and limited magnitude to avoid excessive distortion, the following scaling constraint holds:

illustration not visible in this excerpt

where

- δ(p) is a nonlinear, monotonically increasing function of the fraction of time p.
- σ is the variance of the signal.

Thus, the upper limit of the signal power is given by

illustration not visible in this excerpt

where Q is the weighting matrix and is equal to CTC, and kii are the main diagonal elements of K.

The dynamic range can be expressed in terms of K and W by the relation [4, 7]:

[...]