BioMedSigProcAna

An Iir-Filter Example: A Butterworth Filter Josef Goette Bern University of Applied Sciences, Biel Institute of Human Centered Engineering - microLab [email protected]

February 15, 2016

Contents 1 Introduction

1

2 Analog Butterworth Lowpass-Filters

4

3 Continuous-to-Discrete Transformations 3.1 Impulse Invariance Transformation . . . . . . . . 3.2 Bilinear Transformation . . . . . . . . . . . . . .

10 12 23

4 Discrete-Time Butterworth Filter Example 4.1 Design Using Impulse-Invariance Transformation 4.2 Design Using Bilinear Transformation . . . . . .

29 31 39

References

33 Iir Butterworth

46

i

2016

BioMedSigProcAna

c

Josef Goette, 2007–2016 All rights reserved. This work may not be translated or copied in whole or in part without the written permission by the author, except for brief excerpts in connection with reviews or scholarly analysis. Use in connection with any form of information storage and retrieval, electronic adaptation, computer software is forbidden.

33 Iir Butterworth

ii

2016

BioMedSigProcAna

1

Introduction

In the present document, we use, as usual, the following notation: Continuous-time signal frequencies come in “standard” characters such as f for the frequency in Hertz ([Hz]), and ω for the radian frequency in [rad/sec], ω = 2πf . The corresponding discrete-time signal frequencies use the corresponding symbols with hats: fˆ = ˆ f /fs = f Ts and ω ˆ = ωTs = ω/fs with fs being the sampling frequency in [Hz] and Ts the corresponding sampling interval Ts = 1/fs in [sec]. We abbreviate “continuous-time” by Ct and “discrete-time” by Dt. Further abbreviations are C2d for “continuous-to-discrete,” Tf for “transfer function,” and Lhp for “left-half plane.” We very closely follow [OS75] for the developments in the present document; see our additional remarks on details, which appear after the bibliography at the very end of the document on page 46.

33 Iir Butterworth

1

2016

BioMedSigProcAna

An Iir-Filter Design Approach Iir (recursive) filters are often design by what might be called “analog prototyping” • design a normalized Ct lowpass prototype filter • then: apply frequency-band transformations such as – – – –

lowpass lowpass lowpass lowpass

to to to to

lowpass highpass bandpass bandstop

• then: apply a continuous-to-discrete transformation • here: we work-out a Butterworth-filter example The acronym Iir stands for infinite impulse response; theoretically, the impulse response of this kind of filters never completely dies out. Of course, stable Iir filters have an impulse response that approaches more and more the zero line as time passes by; if the filter runs on a digital computer with floatingpoint arithmetic, then the impulse response will eventually become smaller than the smallest positive representable number; if it runs on a digital computer with fixed-point arithmetic, then the impulse response might begin to oscillate without end—see the discussion to the keyword “limit cycles” in our document [Goe16]. Recall that the acronym Ct means continuous time. Matlab supplies the following commands for frequencyband transformations: lp2lp() for the lowpass-to-lowpass transformation; lp2hp() for the lowpass-to-highpass transformation; 33 Iir Butterworth

2

2016

BioMedSigProcAna lp2bp() for the lowpass-to-bandpass transformation; and the command lp2bs() for the lowpass-to-bandstop transformation, respectively. Matlab has also commands for the continuous-to-discrete transformations, some of which stem from the Control System Toolbox (use help c2d), and some others from the Signal Processing Toolbox (use help impinvar and use help bilinear).1 As the table of contents reveals, we discuss in the present document the design based on the impulse-invariance transformation and the design based on the bilinear transformation. But also note that Matlab supplies “higher level” commands that integrate the mentioned “lower level” commands and thus considerably simplify the design of discrete-time and even digital filters2 from the perspective of the user; you might also want to try out the graphical filter design tool fdatool (use help fdatool to obtain a short description, and type fdatool in the Matlab command window to run the Gui).

1 The mentioned toolboxes just have different names for the commands; the algorithms behind these names are, however, very similar or even identical. 2 Recall that we name filters to be “discrete time” if they filter signals that have the time being discrete but the samples being still real numbers; “digital filters” filter then signals where not only the time is discrete, but where also the samples are represented be a finite number of bits. In the realm of Matlab simulations, we often also call filters and signals to be “discrete,” if the samples of the signals are represented by floating-point numbers—still finitely many bits, but . . . ; and we retain the notion “digital” for signals with samples that are represented by fixed-point numbers.

33 Iir Butterworth

3

2016

BioMedSigProcAna

2

Analog Butterworth Lowpass-Filters

Butterworth Filter Properties • magnitude response is maximally flat in passband ; for a N -th order lowpass, the first (2N − 1) derivatives of the squared magnitude function are zero at ω = 0 • the approximation to the ideal rectangular lowpass characteristic (brick-wall) is monotonic in passband as well as in stopband • squared magnitude function |HCt (jω)|2 =

1 1+



jω jωc

2N

• specified by just 2 parameters – the filter order N – the 3 dB cutoff-frequency ωc

We also often call these analog filters continuous-time (Ct) filters.

33 Iir Butterworth

4

2016

BioMedSigProcAna

Magnitude-Squared and Magnitude Responses .

|HCt (jω)|2 1

1/2

............................................................................. ........ ........ ..... ... ....... ..... .... .. ....... ..... .. ... N = 1, 2, 4, 8 ........ ......... .................. ........... ....................... ................................... ........ .......... ........................... ......... ... ....... .................. ...... ............ ........................................................................................... ......................................................................................................................................... ω

0

.

ωc

. |HCt (jω)| 1

√ 1/ 2

.

................................................................................... ............. .................. ......................... N = 1, 2, 4, 8 .................... ......................... ........................................ ........ .......... ........................... ............................... ... ...... ............. .......................... ...... ... ........ .... ......... ...................................... ....................................... ...... ................ ............... ......................................... ................................................................................................ ω

0

ωc

√ Recall that the √ linear gain 1/ 2 corresponds to −3 dB; we have 20 log10 (1/ 2) = −3.0103. Butterworth filters have a frequency magnitude response with no ripple, neither in the passband nor in the stopband. As the parameter N in the squared magnitude function on page 4 increases, the filter characteristics become sharper, meaning that if N increases, the characteristic in the passband stays closer to unity whereas the characteristic in the stopband comes close to zero more rapidly. Other well known analog filters are the Chebyshev filters, the inverse Chebyshev filters, and the Cauer filters. The Chebyshev 33 Iir Butterworth

5

2016

BioMedSigProcAna filters have a frequency magnitude response with a ripple in the passband but no ripple in the stopband; inverse Chebyshev filters have a dual characteristic to that of the Chebyshev filters: they have no ripple in the passband but a ripple in the stopband; Cauer filters—also called elliptic filters—have a ripple in both, in the passband as well as in the stopband. Chebyshev filters are often also called Chebyshev type I filters, whereas inverse Chebyshev filters are then called Chebyshev type II filters. For a comparable frequency filtering performance, Butterworth filters need the highest order, Chebyshev type I and Chebyshev type II filters need medium orders, and Cauer filters need the lowest order. A higher order means a higher number of operations needed to compute an output sample, and the need of a higher number of storage cells.

33 Iir Butterworth

6

2016

BioMedSigProcAna

Transfer Function Poles • denote by HCt (s) the transfer function of the Butterworth filter • from squared magnitude function we see that (jω → s) 1

HCt (s)HCt (−s) = 1+



s jωc

2N

• poles of squared magnitude function are spk = (−1)1/(2N ) · (jωc )    1 1 + 2k = ωc exp jπ , k = 0, 1, . . . , (2N −1) + 2 2N

33 Iir Butterworth

7

2016

BioMedSigProcAna

Transfer Function Poles: N = 3 Example for N = 3 ; 2N = 6:

1/6

spk = (−1) · j · ωc  1/6 ejπ/2 ωc = ej(π+k2π) jπ/2 = ejπ/6 |ejk2π/6 {z } |e {z } ωc 60◦ 90◦ spacing rotation

j ×

................................... ............ ........ ........ ....... ...... ...... ..... . ..... . . . ... .... . . .... .. . ... . ... ... . ... . . . ... ... . ... .. .... .. .... .. ... .. .. ... .. . ... ... ... .. . ... .. . ... ... ... ... ... .. ... ... ... . . ... .. .... ... ..... .... ...... ...... ....... ...... . . . ........ . . . ............. ..... ...............................

×

.............................. ............ ....... ....... ...... ...... ..... ..... .... . . . .... ... . . ... ... ... . ... .. . ... .... ... ... ... ... .. .. .. ... .. ... ... ... .. . ... . . . ... ... ... ... ... ... ... .... ... . . ..... . .... ...... ...... ....... ....... ......... ......................................

×

×

×

1

×

×

×

×

×ω c

×

2N = 6-th roots of (-1)

×

corresponding Butterworth poles

We see then that the poles of the Butterworth squared magnitude function have an angular spacing of 2π/(2N ) = π/N for an N -th order Butterworth filter. We further have for these poles • that they are distributed around a circle with radius ωc ; • that they are distributed symmetrically on either side of the imaginary axis; • that there are no poles on the imaginary axis itself; • and that there are poles on the real axis if N is odd, but not if N is even. 33 Iir Butterworth

8

2016

BioMedSigProcAna

Transfer Function • to determine the Tf HCt (s) from the Butterworth squaredmagnitude function: perform factorization HCt (s)HCt (−s) • observe that poles in squared magnitude function appear in pairs if s = sp is a pole then s = −sp is also a pole • to construct HCt (s) from squared magnitude function – choose one pole of each pair – choose stable pole (in Lhp) of each pair ....................................... .......... ........ ....... ...... ...... ...... ...... ..... . . . 1 ... ... . ... . .. ... . . ... ... ... . .. ... . ... .. . ... ..... .. ... .. ... ... ... .. ... ... .. ... 0 . c ... .. . ... ... ... ... ... .. ... .. . ... . .. ... ... .... .... ... ∗ ..... .... ...... 1 ..... . . . . . ...... ...... ........ ............ ........ ..................................

× p

poles of N = 3-rd order Butterworth transfer-function: HCt (s) =

×p = −ω

−p0 p1 p∗1 (s − p0 )(s − p1 )(s − p∗1 )

ωc

p ×

33 Iir Butterworth

9

2016

BioMedSigProcAna

3

Continuous-to-Discrete Transformations

C2d-Transformation Procedures there are various procedures to transform an analog prototypefilter design into a discrete-time filter: • impulse-invariance transformation • procedures based on numerical solution of differential equations – first forward-difference (forward Euler) – first backward-difference (backward Euler) – bilinear transformation

We discuss in more detail the impulse-invariance transformation in Subsection 3.1 below. For the procedures based on numerical solution of differential equations we note that the forward Euler transformation is the most simple transformation that has, however, the drawback that an unstable discrete-time system might result from a stable continuous-time (analog) prototype system. The transfer function H(z) of the discrete-time filter is obtained here by replacing the variable s in the continuous-time prototype transfer

33 Iir Butterworth

10

2016

BioMedSigProcAna function HCt (s) by (z − 1)/Ts , where Ts denotes the sampling interval, Ts = 1/fs with fs being the sampling frequency: forward Euler:

s=

 1 1 1 − z −1 z−1 = . Ts Ts z −1

The backward Euler transformation has the advantage over the forward Euler transformation that stable continuous-time designs are transformed into stable discrete-time designs. The transformation is given by backward Euler:

s=

 1 z−1 1 1 − z −1 . = Ts z Ts

The bilinear transformation is the most often used transformation in the design of discrete-time filters;3 therefore, we discuss it in more detail below in Subsection 3.2. For the interested reader we mention that the bilinear transformation is a member of the larger family of Moebius transformations, which are conformal mappings from complex plane to complex plane; a reference is [Hen74, Chapter 5].

3 In the design of simple discrete-time control algorithms, the Euler approximations are likewise very often used.

33 Iir Butterworth

11

2016

BioMedSigProcAna

3.1

Impulse Invariance Transformation

The Impulse Invariance Transformation given: continuous-time (Ct) prototype (analog prototype) HCt (s) •—• hCt (t) = ˆ impulse response of Ct filter find: for the discrete-time filter H(z) •—• h[n] = ˆ impulse response of Dt filter solution: sample the continuous-time impulse response h[n] = ˆ hCt (t = nTs ) where

Ts = ˆ sampling time interval

The above formulae state, reformulated in words, that the impulse response of the discrete-time filter is obtained from the impulse response of the continuous-time (analog) prototype filter through sampling of the latter.

33 Iir Butterworth

12

2016

BioMedSigProcAna

Impulse Invariance: Transfer Functions • it can be shown—see the sampling process—that – the z-transformation of h[n] •—• H(z) – is related to

– the Laplace transformation of hCt (t) •—• HCt (s) by

H(z)

z=esTs

  ∞ 2π 1 X HCt s + j k = Ts Ts k=−∞

• from z = esTs : – strips of width 2π/Ts is s-plane – map onto entire z-plane

33 Iir Butterworth

13

2016

BioMedSigProcAna

Impulse Invariance: Mapping of Planes s-plane ..ℑ(s)

.. .. .. .. . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . .. ... ... .. .. .. .. .. .

z-plane

3π Ts π Ts

− Tπs

ℜ(s)

=⇒

j ................................................................ . . . .................................................. ...................................................................... ......................................... ................................................. 1 ................................... ...................

− 3π Ts

• strips of width 2π/Ts is s-plane map onto entire z-plane • Lhp-part of each strip maps into unit circle

; stable poles of Ct filter go to stable poles of Dt filter

• imaginary axis in s-plane maps onto unit circle such that each segment of length 2π/Ts is mapped once around the circle

We should mention here that the mapping from the s-plane to the z-plane induced by the impulse-invariance procedure is no simple algebraic mapping, as are the mappings that result from the procedures based on numerical solutions of differential equations, an example of which—the bilinear transformation— we discuss in Subsection 3.2.

33 Iir Butterworth

14

2016

BioMedSigProcAna

Impulse Invariance: Frequency Responses • the Dt frequency response is expressed in terms of the Ct prototype frequency response as H z=e

jω ˆ




  ∞ 2π ω ˆ 1 X +j k HCt j = Ts Ts Ts k=−∞

• note: if and only if HCt (jω) ≡ 0 ,

for |ω| ≥

π Ts

then  
ω ˆ 1 HCt j H z = ej ωˆ = Ts Ts else aliasing

33 Iir Butterworth

15

2016

BioMedSigProcAna

Impulse Invariance: Aliasing HCt (jω) 1............ .. ..... . . . . . ..... .. ....... ..... . . . . ......... . . . . . . . . ................ .... . . . . . . . . . . . . ............................................. . . . . . . . . . ...................................... 2πfs ωc ω

1 ........ . .. Ts........ .......

jω ˆ

H(e )

.......... ................. .... ....... .... ... ... .... ... ... . ... ... ... ... ... ... ... ... ... ... ... ... ... . . . . ... . . . ... ... . . .... .. . . . . . . . . . . . .... .... .... .. .. ... . . . . . . . . . ..... . . ..... ..... .. ... ... . . . . . . ...... . . . . . . . . . . . ...... ...... ... ... ... ...... . . . . . . . . . . . . . . . . . . . . . ....... ...... ....... ....... ....... ..... . . . . . . ........ . . . . . . . . . . . . . ..... ..... .............. ..... .............. . . . . . . . . ......... . . . . . . . . . . ............. ............. ....... ....... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...................... .................................. .............. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ................................ ..................................................... ....................................................................................................................................................... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .......................... ..............

−2π

−π

0 ω ˆc

π

2π

ω ˆ

The above figure graphically shows the effects of aliasing: In the upper panel we show the frequency response of a continuoustime filter prototype which is not bandlimited; in the lower panel we show the first three aliasing parts giving, in sum together with all other aliasing parts, the discrete-time frequency response.

33 Iir Butterworth

16

2016

BioMedSigProcAna

Impulse Invariance: Interpretation of Mapping • Ct transfer function: partial fraction decomposition HCt (s) =

N X

Ak , s − sk

N X

Ak esk t u(t)

k=1

• | • hCt (t) =

sk = ˆ poles

k=1

• Dt impulse response and corresponding transfer function h[n] = hCt (nTs ) =

N X

Ak esk nTs u[n] =

k=1

k=1

• | • H(z) =

N X

k=1

N X

Ak eskTs


n

u[n]

Ak 1 − (eskTs ) z −1

We denote by u(t) the continuous-time unit-step function and, correspondingly, by u[n] the discrete-time unit-step function. Note that the above development is true if all poles have multiplicity 1; for poles with higher multiplicities we must use the corresponding partial fraction decompositions.

33 Iir Butterworth

17

2016

BioMedSigProcAna

Impulse Invariance: Interpretation of Mapping (2) s-plane poles sk • transformation of poles:

transform to z-plane poles zk = esk Ts

• coefficients Ak are equal for Ct and Dt • if Ct filter is stable, then Dt filter is also stable • we note however that – although poles are mapped by zk = esk Ts – the complete planes are not mapped by that relation – for example: zeros of Dt transfer function are functions of the coefficients Ak and of the poles sk

Concerning the stability of the filters we know that the continuous-time filter is stable if and only if all of its poles are in the left-half complex s-plane, real{sk } < 0. By the given mapping of the poles we then have for the discrete-time poles zk that |zk | = esk Ts < 1, meaning that the discrete-time poles are inside of the unit circle, and in turn, that the discrete-time filter is also stable. It is important to note for the impulse-invariant transformation that, although the poles are mapped from the s-plane to the z-plane by the relation zk = esk Ts , the planes themselves are not mapped by that relation. For example, zeros are not mapped in that way; the following simple example gives an illustration.

33 Iir Butterworth

18

2016

BioMedSigProcAna Example. Consider the second-order continuous-time transfer function HCt (s) = =

s+a (s + a)2 + b2 1/2 1/2 + , s + a + jb s + a − jb

having one zero at s = −a, the second zero at infinity, and a complex-conjugate pole pair at s = −a ± jb. The discrete-time transfer function obtained from the impulseinvariance transformation becomes H(z) =

1/2 1−

e−aTs e−jbTs z −1

+

1/2 1−

e−aTs e+jbTs z −1


+ e−jbTs /2 z −1 e 1−e = (1 − e−aTs e−jbTs z −1 ) (1 − e−aTs e+jbTs z −1 ) −aTs

+jbTs

1 − e−aTs cos(bTs )z −1 1 − e−aTs (e+jbTs + e−jbTs ) z −1 + e−2aTs z −2
z z − e−aTs cos(bTs ) = 2 . z − 2e−aTs cos(bTs )z + e−2aTs =

This discrete-time transfer function has one zero at the origin and the other zero at z = e−aTs cos(bTs ). We thus see that, although the poles are mapped by zk = esk Ts , the zeros are not mapped according to this formula.

33 Iir Butterworth

19

2016

BioMedSigProcAna

Impulse Invariance: Practical Advice • if Ct prototype is “sufficiently bandlimited,” then  
ω ˆ 1 H ej ωˆ ≈ HCt j Ts Ts • thus: “high” sampling rates (small Ts ) will minimize the aliasing, but the Dt filter may obtain a very high gain • advice: instead using H(z) =

N X

k=1

Ak 1 − (eskTs ) z −1

N X

Ak Ts 1 − (eskTs ) z −1

use H(z) =

k=1

• the Dt impulse response then is h[n] = Ts hCt (t = nTs )

33 Iir Butterworth

20

2016

BioMedSigProcAna

Impulse Invariance: Generalizations • one motivation to use the impulse-invariance procedure if Ct prototype filter is bandlimited then Dt filter frequency response closely approximates the Ct frequency response • another motivation to use the procedure: control some aspects of the time response of Dt filter – step invariance procedure – waveform invariance: extend the concept to preserve the output wave-shape for a variety of inputs • final remark: besides aliasing, the impulse-invariance approach transforms the frequency responses linearly

The mentioned step-invariance procedure just obtains the step response of the discrete-time filter by sampling the step response of the continuous-time prototype filter. The resulting discrete-time filter then might have desired step-response characteristics such as small rise time and low peak overshoot. An important feature of the impulse-invariance procedure is that, besides aliasing, the frequency responses of the continuoustime prototype filter and its discrete-time counterpart are linearly related, meaning that the shape of continuous-time frequency response is preserved in the discrete-time filter. This result is in contrast to the procedures which use algebraic transformations, an example of which—the bilinear transformation— we discuss in Subsection 3.2 below. We recall however, that the 33 Iir Butterworth

21

2016

BioMedSigProcAna impulse-invariance design technique is only appropriate for essentially bandlimited filters. Thus, the design of highpass or bandstop filters requires additional bandlimiting to avoid severe aliasing distortions.

33 Iir Butterworth

22

2016

BioMedSigProcAna

3.2

Bilinear Transformation

Bilinear Transformation Formulae given: Tf of analog prototype HCt (s) find: Tf H(z) of corresponding discrete-time filter solution: replace s in HCt (s) by s=

2 1 − z −1 Ts 1 + z −1

the inverse formula then is z=

1 + (Ts /2)s 1 − (Ts /2)s

You might want to verify the given formulae as follows: Start with a first-order continuous-time prototype system and first specify it by its differential equation; second, express the differential equation by its corresponding integral equation formulated for the n-th sampling time instant and its past neighbor (n − 1); third, approximate the involved integral by the trapezoidal rule; forth, replace the appearing derivatives by the expression given by the differential equation started with; fifth, rearrange terms to obtain the difference equation of the discretetime system, and from it obtain the transfer function H(z). Generalize to systems of order higher than one.

33 Iir Butterworth

23

2016

BioMedSigProcAna

Bilinear Transformation: Frequency Mapping • Ct frequencies ω ; on imaginary axis in s-plane: s = jω • Dt frequencies ω ˆ ; on unit circle in z-plane: z = ej ωˆ • unit circle maps to imaginary axis (and vice versa)
2 1 − e−j ωˆ s z = ej ωˆ = Ts 1 + e−j ωˆ = ··· 2 = j tan (ˆ ω/2) Ts = jω • therefore, frequency mapping is ω=

2 tan (ˆ ω /2) Ts

⇐⇒

ω ˆ = 2 arctan (ωTs /2)

To fill-in the dots in the above derivation, first use  ωˆ  ω ˆ ω ˆ 1−e−j ωˆ = e−j 2 ej 2 − e−j 2 ,

 ωˆ  ω ˆ ω ˆ 1+e−j ωˆ = e−j 2 ej 2 + e−j 2 ;

next use Euler’s formulae ω ˆ

ω ˆ

ej 2 − e−j 2 sin(ˆ ω /2) = , 2j

ω ˆ

ω ˆ

ej 2 + e−j 2 cos(ˆ ω /2) = , 2

to obtain j sin(ˆ ω /2)/ cos(ˆ ω /2) which is equal to j tan(ˆ ω /2). 33 Iir Butterworth

24

2016

BioMedSigProcAna

Bilinear Transformation: Frequency Mapping • ω ˆ= ˆ Dt frequency

←→

ω= ˆ Ct frequency

• the frequency mapping is ω ˆ = 2 arctan (ωTs /2) ω ˆ

....................................... ........................... . . . . . . . . . ...... ...... .. ... ω .. . .. ... . . . . .... ............. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ............... −π π

33 Iir Butterworth

25

2016

BioMedSigProcAna

Bilinear Transformation: Mapping of Planes • the mapping s 7→ z is z=

1 + (Ts /2)s 1 − (Ts /2)s

s-plane .

.. ..ℑ(s) .. .. .. .. .. .. ℜ(s) =⇒ .. .. .. .. .. .. .. ..

. .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . .. . . .. . . . . . . . . . . . . . . . .

33 Iir Butterworth

26

z-plane image of s .= .. jω

. ....................... ... .............................................................. ........................................ ................................................................................... ......................................1 ........................................ .............................. ... ................ ...j........

image of lefthalf plane

2016

BioMedSigProcAna

Bilinear Transformation: Mapping of Circles • in s-plane: Butterworth poles on a circle with radius ωc , equally spaced in angle • bilinear transformation is conformal mapping: Butterworth circle in s-plane maps to a circle in z-plane s-plane .

... .. .. .. .. .. .. .. .. .. .. .. .. .. .. .

z-plane

. .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......... . . . . . . . . . . . . . . . . . . ................... . . .............. ...... . . . . . . . . . . . . . . . ............ . . . . . . . .... . . . . . . . . . . . . . . ...... . . . . . . . . . . ... . . . . . . . . . . . . ....... . . . . . . . . . . . . ... . . . . . . . . . . . ..... . . . . . . . . . . . . . ... . . . . . . . . . . . .... . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . .. . . . . . . . . . . ..... . . . . . . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .... . . . . . . . . . . . . . . . . . . . . . . . . . . . .... . . . . . . . . . . . . . ... . . . . . . . . . . . . ...... . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . ........ . . . . . . . . . . . . . . . . . . . . . . . . . . . .......... . . . . . . . . ..... . . . . . . . . . . . . . . . . . . ...................... . .............. . . . . . . . . . . . . . . . . . . . . . . . . . ....... . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . .. . . . .. . . .. . . . .. . . .. . . . . . . . . . . . . . . . . . .

×

×

×

×

×ωc =⇒

×

Butterworth N = 3, 2N = 6

j .. ...... .. ......................................................................× . ........................ . . . . ............................................................ .. ............. ... .... . . ... . . . . . . . . ............× .. . .. . . . ... ...........................................................................× × ..... . .. . .. . ... ... ........................................................1 ... . . .............................×............................. ... ..... ................................... .........× ................... ..................................... 1−ωc Ts/2 1+ωc Ts/2

1+ωc Ts/2 1−ωc Ts/2

Note, however, that the Butterworth circle in the z-plane is neither centered at the origin, nor are the poles equally spaced in angle. But the left-half s-plane poles map into poles inside of the unit circle; therefore—as is always the case with the bilinear transformation—we obtain a stable discrete-time Butterworth filter from a continuous-time Butterworth prototype.

33 Iir Butterworth

27

2016

BioMedSigProcAna

0

π

.......... −3 dB H(ejωˆ )[dB] ... ... .... .... ..... ..... ...... ...... ...... ..... ..... .... ... ω ˆ 3dB ω ˆ

Bilinear Transformation: Frequency Responses

33 Iir Butterworth

ω ˆ π

0

................. ........................ . . . . . . . . . . . . . . . . . ............. ............. . . . . . . . ....... ..... ...... ... . . . . . ... . . . ω ˆ = 2 arctan (ωTs /2) ... . . . ... .... .... ω ......... −3 dB H Ct (jω)[dB] ... ... .... .... ..... ....... ........ ......... ........... ............... .................. .................... ............... ω3dB ω

28

2016

BioMedSigProcAna

4

Discrete-Time Butterworth Lowpass Filter Example

Frequency-Response Specifications • here: specifications in the Dt frequency domain • requirements: – passband magnitude constant to within 1 dB – passband is 0 ≤ ω ˆ≤ω ˆp = ˆ 0.2π – stopband attenuation > 15 dB

– stopband is 0.3π = ˆ ω ˆs ≤ ω ˆ ≤π • thus, if passband magnitude is normalized to 1 20 log10 H(ej0.2π ) ≥ −1 20 log10 H(ej0.3π ) ≤ −15

• for convenience we may assume that the sampling-time interval is Ts = ˆ 1

Note that we must distinguish between two different application scenarios: The first scenario designs a discrete-time filter which is specified in the discrete-time frequency domain. This is the situation we consider here and which allows to set the sampling time equal to one, because there exists no real continuoustime filter. 33 Iir Butterworth

29

2016

BioMedSigProcAna The second application scenario emulates a continuous-time filter with a discrete-time filter. Here the specifications are given in the continuous-time frequency domain, and our design has to select an appropriate sampling time interval Ts with which we plan to implement the emulating discrete-time filter. Obviously, in this scenario the sampling time interval cannot be set to unity.

33 Iir Butterworth

30

2016

BioMedSigProcAna

4.1

Design Using Impulse-Invariance Transformation

Analog Prototype Specifications • transform the Dt frequency specifications to corresponding specifications of a Ct prototype • recall: – impulse invariance design introduces aliasing – beside aliasing it linearly maps Ct to Dt frequencies • convenient procedure: – assume that aliasing is negligible – carry out design – verify performance of resulting filter • thus, mapping of critical frequencies is ω ˆ p = 0.2π

ω ˆ s = 0.3π

;

ωp = ω ˆ p /Ts = ω ˆp

;

20 log10 |HCt (jωp )| ≥ −1

;

ωs = ω ˆ s /Ts = ω ˆs

;

20 log10 |HCt (jωs )| ≤ −15

Recall that we have assumed that the sampling time interval Ts is unity. 33 Iir Butterworth

31

2016

BioMedSigProcAna

Analog Butterworth Parameters • the Ct Butterworth squared magnitude function is |HCt (jω)|2 =

1 1+



ω ωc

2N

• to do: determine the needed two parameters N and ωc • solving given specification inequalities with equality4  0.1  10 −1 1 log 101.5 −1
= 5.8858 N= 0.2 2 log 0.3 • to meet specifications: select N = 6 • inserting N = 6 into passband equation gives ωc = 0.7032

Note that because we must round up N to the nearest integer, not both specifications, passband and stopband, can be met exactly. If we insert N = 6 into the passband equation, the passband specifications are met exactly and the stopband specifications are exceeded for the Ct filter. Such a choice allows some margin for the aliasing that enters in the Dt filter. We next supply the detailed steps leading to the above results. We start with the two equations in the two unknowns 4 See

the derivations in the text below.

33 Iir Butterworth

32

2016

BioMedSigProcAna N and ωc that we obtain by setting the design specification inequalities to equality: |HCt (jωp )|2 =

1 1+



ωp ωc

−0.1 , 2N = 10

; |HCt (jωs )|2 =

1 1+



ωs ωc

1+



ωp ωc

2N

(1a)

= 100.1 . (1b)

−1.5 , 2N = 10

;

1+



ωs ωc

2N

(1c)

= 101.5 . (1d)

Taking in (1b) and (1d) logarithms (to any base) we next obtain
log 100.1 − 1 = = log (ωp ) − log (ωc ) , log 2N
  log 101.5 − 1 ωs log = = log (ωs ) − log (ωc ) . ωc 2N 

ωp ωc



(2a) (2b)

Subtracting (2b) from (2a) we obtain

 1  log (ωp ) − log (ωs ) = log 100.1 − 1 − log 101.5 − 1 2N | | {z } {z }  0.1  = log (ωp /ωs ) 10 − 1 = log 101.5 − 1

33 Iir Butterworth

33

2016

BioMedSigProcAna and in turn



100.1 − 1 log 1 101.5 − 1   N= ωp 2 log ωs



.

(3)

Finally, we insert ωp = ω ˆ p = 0.2π and ωs = ω ˆ s = 0.3π into (3). To numerically obtain the second filter parameter, ωc = 0.7032, via the passband specifications, we use the parameter N = 6 and ωp = 0.2π in (2a).

Analog Butterworth Poles • for N = 6 there are 2N = 12 poles of Butterworth squaredmagnitude function • these 12 poles are uniformly distributed in angle on a circle with radius ωc = 0.7032 • the Butterworth transfer function uses the N = 6 poles in the Lhp s-plane × ×

.............................. ........ ...... ...... ..... ..... ... ... ... . . ... ... . ... . . . .. .... .. 0.7032 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ... . . ........... .......... ... . . . .. . . . . . . ................ ..... .. . . . ....... .... ... . . . . . . . . . . . . . . . .... .. ... . ..... ... ..... ... . ... π ...... .................... .... ..... ..... .. ....... .......... .................................. 6 ....

×

×

×

×

×

×

×

33 Iir Butterworth

× ×

34

×

2016

BioMedSigProcAna

Analog Butterworth Transfer Function • the Lhp poles are p1 = −0.6792 + j · 0.1820

p∗1 = −0.6792 − j · 0.1820

p2 = −0.4972 + j · 0.4972

p∗2 = −0.4972 − j · 0.4972

p3 = −0.1820 + j · 0.6792

p∗3 = −0.1820 − j · 0.6792

• the corresponding transfer function becomes HCt (s) = =

p1 p∗1 p2 p∗2 p3 p∗3 (s − p1 ) (s − p∗1 ) (s − p2 ) (s − p∗2 ) (s − p3 ) (s − p∗3 ) 0.1209 (s2 + 1.3585s + 0.4945) (s2 + 0.9945s + 0.4945) ·

(s2

1 + 0.3640s + 0.4945)

We have normalized the above analog transfer function HCt (s) such that it has a Dc gain of one.

33 Iir Butterworth

35

2016

BioMedSigProcAna

Discrete-Time Butterworth Transfer Function • express HCt (s) as a partial fraction decomposition ; Ak , sk • apply H(z) =

N X

k=1

Ak 1 − (eskTs ) z −1

• to obtain H(z) =


1.8557 − 0.6304z −1 (1.0000 − 0.9973z −1 + 0.2571z −2)


−2.1428 + 1.1454z −1 + (1.0000 − 1.0691z −1 + 0.3699z −2)
0.2871 − 0.4466z −1 + (1.0000 − 1.2972z −1 + 0.6949z −2)

We have obtained the above parallel form of second-order sections by combining the terms of complex-conjugate pole pairs. Obviously, we might directly use this parallel form that naturally results from the impulse-invariant design procedure. If we desire a cascade form or a direct form, we must combine the separate second-order terms in an appropriate way.

33 Iir Butterworth

36

2016

BioMedSigProcAna

Discrete-Time Butterworth Magnitude-Response

1.0

0.5

0

H(ej ωˆ )

...................................... ... 0.8913 = ˆ −1 dB ... ... ... ... ... ... .... 0.1778 = ˆ −15 dB ........ ............................................................................................................................ 0.2π 0.3π

0 −15

−80

ω

π

H(ej ωˆ )

...................................dB ............ 0.8913 = ˆ −1 dB ....... ....... 0.1778 = ˆ −15 dB ........ ........ ......... ........... ............ .............. ................ ................. ........................ 0.2π 0.3π

33 Iir Butterworth

ω

π

37

2016

BioMedSigProcAna

Discrete-Time Butterworth Phase-Response π

0

−π

arg(H(ej ωˆ ))

.... ... ..... . . ... ..... .. ... ... ...... .... . .... ..... ω ... .... .0.2π 0.3π.............. π .... . ............ .... .. ................. .... .. ........................ ... .. ........................... ... . .................. ... .. ........... .... ........ .

Final Remarks to Our Impulse-Invariance Design. We recall that we have designed the filter to exactly meet—with the assumption that we have no aliasing—the passband specifications, that is, in the Butterworth case, to have a −1 dB attenuation at the passband edge frequency ω ˆ p = 0.2π; the design then exceeds the specification at the stopband edge frequency ω ˆ s = 0.3π. As we observe from the magnitude-response plots on page 37, it is true that at the passband edge the attenuation is slightly below −15 dB, indicating that the aliasing is not too strong, or, in other words, that the continuous-time filter is sufficiently bandlimited. In an other design this might not be true, such that the resulting discrete-time filter does not meet the specifications. If we have that situation, we may try to differently adjust the filter parameters (holding the order fixed), or we may try again with a higher-order filter.

33 Iir Butterworth

38

2016

BioMedSigProcAna

4.2

Design Using Bilinear Transformation

Analog Prototype Specifications • Dt frequency specifications must be pre-warped to corresponding analog frequencies, such that critical analog frequencies map to correct Dt frequencies • the frequency mapping function of the bilinear transformation is ω = (2/Ts ) tan (ˆ ω/2) = 2 tan (ˆ ω/2) • thus, mapping of critical frequencies is ω ˆ p = 0.2π

ω ˆ s = 0.3π

;

ωp = 2 tan (ˆ ωp /2) = 2 tan (0.1π)

;

20 log10 |HCt (jωp )| ≥ −1

;

ωs = 2 tan (ˆ ωs /2) = 2 tan (0.15π)

;

20 log10 |HCt (jωs )| ≤ −15

Recall that we conveniently have assumed that the sampling time interval is Ts = 1.

33 Iir Butterworth

39

2016

BioMedSigProcAna

Analog Butterworth Parameters • the Ct Butterworth squared magnitude function is |HCt (jω)|2 =

1 1+



ω ωc

2N

• to do: determine the needed two parameters N and ωc • solving given specification inequalities with equality5  0.1  10 −1 1 log 101.5 −1  = 5.3044  N= 2 log tan(0.1π) tan(0.15π)

• to meet specifications: select N = 6 • inserting N = 6 into stopband equation gives ωc = 0.7662

To supply detailed steps leading to the above results we may start with (3) on page 34. Using for the bilinear transformation ωp = 2 tan (ˆ ωp /2) = 2 tan (0.1π) and ωs = 2 tan (ˆ ωs /2) = 2 tan (0.15π), we obtain

N=

5 See

1 log 2



100.1 −1 101.5 −1

log



ωp ωs





 0.1  10 −1 log 1 101.5 −1  = 5.3044 .  = tan(0.1π) 2 log tan(0.15π)

the derivations in the accompanying text.

33 Iir Butterworth

40

2016

BioMedSigProcAna Because the order N of the Butterworth filter must be an integer, we must select N = 6 in order to meet the specifications. If we insert N = 6 into (2b), we can determine ωc to ωc = 0.7662. Note that by using the equation (2b)—which comes from the stopband constraint—to determine the ωc parameter, we meet the stopband specifications exactly and exceed the passband specifications. Such a choice is reasonable, because the bilinear transformation has no aliasing effects.

Analog Butterworth Poles • for N = 6 there are 2N = 12 poles of Butterworth squaredmagnitude function • these 12 poles are uniformly distributed in angle on a circle with radius ωc = 0.7662 • the Butterworth transfer function uses the N = 6 poles in the Lhp s-plane ...................................... ....... ...... ...... .... .... .... . . . ... ... ... .. ... . ... .... .. .... 0.7662 . . . . . . . . . . ... . ............ ................................................ . . . . . . ... . . . . ...... ..... . . . . . . . . . . . . .. . ....... . . . . . . . . . . . . . . ... .... ... .. . ..... ... ... .. .. . .. ... π ...... ................. ..... ..... .. ...... . . . . . . .... . ........ ................................... 6 .........

×

×

×

×

×

×

×

×

×

×

×

33 Iir Butterworth

×

41

2016

BioMedSigProcAna

Analog Butterworth Transfer Function • the Lhp poles are p1 = −0.7401 + j · 0.1983

p∗1 = −0.7401 − j · 0.1983

p2 = −0.5418 + j · 0.5418

p∗2 = −0.5418 − j · 0.5418

p3 = −0.1983 + j · 0.7401

p∗3 = −0.1983 − j · 0.7401

• the corresponding transfer function becomes HCt (s) = =

p1 p∗1 p2 p∗2 p3 p∗3 (s − p1 ) (s − p∗1 ) (s − p2 ) (s − p∗2 ) (s − p3 ) (s − p∗3 ) 0.2023 (s2 + 1.4802s + 0.5871) (s2 + 1.0836s + 0.5871) ·

(s2

1 + 0.3966s + 0.5871)

Again, we have normalized the above analog transfer function HCt (s) such that it has a Dc gain of one.

33 Iir Butterworth

42

2016

BioMedSigProcAna

Discrete-Time Butterworth Transfer Function • apply the bilinear transformation s=

2 1 − z −1 Ts 1 + z −1

• to the analog Butterworth transfer function HCt (s) to obtain
6 7.3769 · 10−4 1 + z −1 H(z) = (1.0000 − 0.9044z −1 + 0.2155z −2) · ·

33 Iir Butterworth

1 (1.0000 − 1.0106z −1 + 0.3583z −2) 1 (1.0000 − 1.2687z −1 + 0.7051z −2)

43

2016

BioMedSigProcAna

Discrete-Time Butterworth Magnitude-Response

1.0

0.5

0

H(ej ωˆ )

........................................ ... 0.9372 = ˆ −0.5634 dB ... ... ... ... ... ... .... 0.1778 = ˆ −15 dB ........ ........................................................................................................................... 0.2π 0.3π

0 −15

−80

ω

π

H(ej ωˆ )

...................................dB .............. 0.9372 = ˆ −0.5634 dB ....... ....... 0.1778 = ˆ −15 dB ....... ....... ....... ....... ....... ....... ....... ....... ....... 0.2π 0.3π

33 Iir Butterworth

ω π

44

2016

BioMedSigProcAna

Discrete-Time Butterworth Phase-Response π

0

−π

arg(H(ej ωˆ ))

.... ... ..... . . ... ..... .. ... ... ...... .... . .... ω ..... ... .... . 0.2π 0.3π......... π .... . ........ .... . ........... .... .. ............. . .... . ................ ... .. .................... ... .. ....................... .... .......................... ......

Final Remarks to Our Bilinear-Transformation Design. We recall that we have designed the filter to exactly meet the stopband specifications, that is, in the Butterworth case, to have a −15 dB attenuation at the stopband-edge frequency ω ˆ s = 0.3π; the design exceeds the specification at the passband-edge frequency ω ˆ p = 0.2π. As we observe from the magnitude-response plots on page 43, it is true that at the stopband edge the attenuation is at −15 dB, whereas at the passband edge the attenuation is only about −0.5634 dB, leaving a certain margin to the −1 dB required by the specifications. Comparison of the two Designs. If we compare the bilineartransformation based design on page 43 to the design based on the impulse-invariance transformation on page 37, we see

33 Iir Butterworth

45

2016

BioMedSigProcAna that the magnitude function of the bilinear transformation design falls off more rapidly than the magnitude function of the impulse-invariance transformation design. This is because the bilinear transformation maps the entire jω axis of the s-plane onto the unit circle, and the continuous-time Butterworth filter of order 6 has a zero of multiplicity 6 at s → ∞; the resulting discrete-time filter then has a zero of multiplicity 6 at z = −1.

References [Goe16] Josef Goette. Biomedical Signal Processing and Analysis—On Fixed-Point Filter Realizations. Bern University of Applied Sciences, Script at the Bfh-ti Biel/Bienne, HuCE-microLab, February 2016. [Hen74] Peter Henrici. Applied and Computational Complex Analysis, volume 1 of Pure and Applied Mathematics. John Wiley & Sons, New York, 1974. [OS75] Alan W. Oppenheim and Ronald W. Schafer. Digital Signal Processing. Prentice-Hall Inc., Englewood Cliffs, N.J., 1975. [Ran02] Rangaraj M. Rangayyan. Biomedical Signal Analysis: A Case-Study Approach. IEEE Press, New York, 2002. Bfh-ti Biel/Bienne Library 57.08 RANGA. In our development, we have mainly followed [OS75, Chapter 5, Sections 5.1 and 5.2]; there you also find design examples for Chebyshev and Cauer lowpass filters, as well as examples for using frequency transformations to design filters with highpass, bandpass, and bandstop characteristics. On Butterworth filters you might also want to consult [Ran02, Chapter 3, pp. 118 ff.]. 33 Iir Butterworth

46

2016