Download Butterworth Filter Design: A Comprehensive Guide with Examples and more Lecture notes Digital Signal Processing in PDF only on Docsity!
Digital Signal Processing – 5
V. IIR Digital Filters
(Deleted in 2007 Syllabus). (Added in 2007 Syllabus).
2007 Syllabus: Analog filter approximations – Butterworth and Chebyshev, Design of IIR digital filters from analog filters, Bilinear transformation method, Step and Impulse invariance techniques, Spectral transformations, Design examples: Analog-Digital transformations
Contents: 5.1 Introduction 5.2 The normalized analog, low pass, Butterworth filter 5.3 Time domain invariance 5.4 Bilinear transformation 5.5 Nonlinear relationship of frequencies in bilinear transformation 5.6 Digital filter design – The Butterworth filter 5.7 Analog design using digital filters 5.8 Frequency transformation 5.9 The Chebyshev filter 5.10 The Elliptic filter
5.1 Introduction
Nomenclature With a 0 = 1 in the linear constant coefficient difference equation, a 0 y(n) + a 1 y(n–1) + … + aN y(n–N) = b 0 x(n) + b 1 x(n–1) + … + bM x(n–M) , a 0 0 we have,
H(z) =
N
i
i i
M
i
i i
a z
bz
1
0
1
This represents an IIR filter if at least one of a 1 through aN is nonzero, and all the roots of the denominator are not canceled exactly by the roots of the numerator. In general, there are M finite zeros and N finite poles. There is no restriction that M should be less than or greater than or equal to N. In most cases, especially digital filters derived from analog designs, M ≤ N. Systems of this type are called Nth^ order systems. This is the case with IIR filter design in this Unit. When M > N, the order of the system is no longer unambiguous. In this case, H(z) may be taken to be an Nth^ order system in cascade with an FIR filter of order ( M – N ). When N = 0, as in the case of an FIR filter, according to our convention the order is 0. However, it is more meaningful in such a case to focus on M and call the filter an FIR filter of M stages or ( M +1) coefficients.
Example The system H(z) = ( 1 z ^8 ) ( 1 z ^1 )is not an IIR filter. Why (verify)?
IIR filter design An analog filter specified by the Laplace transfer function, Ha(s) , may be designed to either frequency domain or time domain specifications. Similarly, a digital filter, H(z) , may be required to have either (1) a given frequency response, or (2) a specific time domain response to an impulse, step, or ramp, etc.
Analog design using digital filters, ωi = ΩiT Another possibility is that a digital filter may be required to simulate a continuous-time (analog) system. To simulate an analog filter the discrete- time filter is used in the A/D – H(z) – D/A structure shown below. The A/D converter can be thought of roughly as a sampler and coder, while the D/A converter, in many cases, represents a decoder and holder followed by a low pass filter ( smoothing filter ). The A/D converter may be preceded by a low pass filter, also called an anti-aliasing filter or pre-filter.
We will usually be given a set of analog requirements with critical frequencies Ω 1 , Ω 2 ,…, ΩN in radians/sec., and the corresponding frequency response magnitudes K 1 , K 2 , …, KN in dB. The sampling rate 1/ T of the A/D converter will be specified or can be determined from the input signals under consideration.
Equivalent Analog Filter
x(n) y(n)
xa(nT) ya(nT)
(1/ T ) samples/sec (1/ T ) samples/sec
A/D H(z) D/A
xa(t) ya(t)
5.2 The normalized analog, low pass, Butterworth filter
As a lead-in to digital filter design we look at a simple analog low pass filter – an RC filter, and its frequency response.
Example 5.2.1 Find the transfer function, Ha(s) , impulse response, ha(t) , and frequency response, Ha(jΩ) , of the following system.
Solution This is a voltage divider. The transfer function is given by
Ha(s) = ()
X s
Ys
R sC
sC 1
s RC
RC
s
Taking the inverse Laplace transform gives the impulse response,
ha(t) = 5 e ^5 tu ( t )
The frequency response is
Ha(jΩ) = H (^) a ( s ) s j = 5
j
= 2 tan( /^5 )
1 1 ( / 5 )
e^ j
The cut-off frequency is Ωc = 5 rad/sec. The gain at Ω = 0 is 1.
The MATLAB plots of frequency response of Ha(jΩ) = 5 ( j 5 )are shown below. We
use the function fplot. The analog frequency, Omega ( Ω) , extends from 0 to ∞; however, the plots cover the range 0 to 6 π rad/sec.
subplot(2, 1, 1), fplot('abs(5/(5+jOmega))', [-6pi, 6pi], 'k'); xlabel ('Omega, rad/sec'), ylabel('|H(Omega)|'); grid; title ('Magnitude') % subplot(2, 1, 2), fplot('angle(5/(5+jOmega))', [-6pi, 6pi], 'k'); xlabel ('Omega, rad/sec'), ylabel('Phase of H(Omega)'); grid; title ('Phase')
Ω , rad./sec.
|Ha(jΩ)|
Ωc = 5
Input = x(t) C^ = 20μF Output = y(t)
R = 10kΩ
-15 -10 -5 0 5 10 15
1
Omega, rad/sec
|H(Omega)|
Magnitude
-15 -10 -5 0 5 10 15
0
1
2
Omega, rad/sec
Phase of H(Omega)
Phase
If we adjust the values of the components R and C so that 1 RC = 1, we would have Ha(s)
= 1
s
which is a normalized filter with cut-off frequency Ωc = 1 rad/sec and gain of 1 at Ω = 0.
Such a normalized LP filter could be transformed to another LP filter with a different cut-off frequency of, say, Ωc = 10 rad/sec by the low pass to low pass transformation s → s / 10 . The
transfer function then becomes
Ha(s) = (^1) (/ 10 )
s s s
s
which still has a dc gain of 1. The gain of this filter could be scaled by a multiplier, say, K , so that
Ha(s) = 10
s
K
which has a dc gain of K and a cut-off frequency of Ωc = 10 rad/sec. The frequency response of the normalized filter Ha(s) = 1 /( s 1 )is Ha(jΩ) = 1 ( j 1 ). The
corresponding MATLAB plots are shown below using the function plot. Omega is a vector , consequently we use “./” instead of “/” etc.
Omega = -6pi: pi/256: 6pi; H = 1./(1.+ j .*Omega); subplot(2, 1, 1), plot(Omega, abs(H), 'k'); xlabel ('Omega, rad/sec'), ylabel('|H(Omega)|'); grid; title ('Magnitude') subplot(2, 1, 2), plot(Omega, angle(H), 'k'); xlabel ('Omega, rad/sec'), ylabel('Phase of H(Omega)'); grid; title ('Phase')
0 5 10 15
1
Omega, rad/sec
|H(Omega)|
Magnitude
1st Order, Cut-off = 5 rad/sec
Frequency response analysis The frequency response analysis in the analog frequency ( Ω ) domain is given by the following equations which may be used to illustrate, qualitatively, the effect of LP, HP or BP analog filters on a signal.
Y(s) = H(s) X(s) Y(jΩ) = H(jΩ) X(jΩ) = H ( j ) ej H ( j ) X ( j ) ej X ( j )= H ( j ) X ( j ) ej ^ H ( j ) X ( j Y () = H () X () and Y ()= H ()+ X ()
Similarly, if we have a digital filter H(z) the frequency response analysis in the digital frequency ( ω ) domain is given by the following equations which may be used to illustrate, qualitatively, the effect of LP, HP or BP digital filters on a signal.
Y(z) = H(z) X(z)
Y(jω) = H(jω) X(jω) = H^ ( j^ ^ )^ ej H ( j^ ) X^ ( j^ ^ )^^ ej X ( j^ )= H^ (^ j ^ ) X ( j^ ^ )^ ej ^ H ( j ) X ( j^
Y ( ) = H ( ) X ( ) and Y ( )= H ( )+ X ( )
X(s) Y(s) H(s)
X(z) Y(z) H(z)
Writing down the Nth^ order filter transfer function H(s) from the pole locations Let us look at some analog Butterworth filter theory.
(1) |H(jΩ)| = 1 / c ^2 N
decreases monotonically with Ω. No ripples.
(2) Poles lie on the unit circle in the s -plane (for the Chebyshev filter, in contrast, they lie on an ellipse). (3) The transition band is wider (than in the case of the Chebyshev filter). (4) For the same specifications the Butterworth filter has more poles (or, is of higher order) than the Chebyshev filter. This means that the Butterworth filter needs more components to build.
The normalized analog Butterworth filter has a gain of |H(jΩ)| = 1 at Ω = 0, and a cut-off frequency of Ωc = 1 rad/sec. Given the order, N , of the filter we want to be able to write down its transfer function from the pole locations on the Butterworth circle.
Example 5.2.2 Given the order N of the filter, divide the unit circle into 2 N equal parts and place poles on the unit circle at (360^0 /2 N ) apart. The H(s) will be made up of the N poles in the left half plane only. Remember complex valued poles must occur as complex conjugate pairs. There will
be no poles on the imaginary axis. Since the N poles must lie on the left half semicircle, when N is odd the odd pole must be at s = – 1. Thus, for N = 1, there is one pole, s 1 = – 1, and H(s) is given by
Radius = Ωc = 1
ζ
jΩ
s -plane
Ω , rad./sec.
|H(jΩ)|
Ωc = 1
Direction of increasing order of filter
H(s) = 1
s s
s
s
Example 5.2.3 Filter order N = 2 so that 2 N = 4 and 360^0 /4 = 90^0. The pole plot is shown above. The poles are at
s 1 =
j and s 2 =
j ,
so that
H(s) = ( )( )
s s 1 s s 2
s j s j
Denominator is
Dr. =
2
2
s –
2
2
j = 2
s s j = s^2 2 s 1
So
H(s) = 2 1
s^2 s
Example 5.2.4 Filter order N = 3, so that 2 N = 6 and 360^0 /6 = 60^0. Poles are at
s1, 2, 3 = – 1,
j , and
j s + 2 s+ 1
2
Radius = Ωc = 1
1
ζ
jΩ
s -plane
s 3 = s * 2
s 2
s 1
- s 3
Radius = Ωc = 1
ζ
jΩ s -plane
Determining the order and transfer function from the specifications A typical magnitude response specification is sketched below. The magnitudes at the critical frequencies Ω 1 and Ω 2 are A and B , respectively. Typically Ω 1 is in the pass band or is the edge of the pass band and Ω 2 is in the stop band or is the edge of the stop band. For illustrative purposes we have arbitrarily
taken A = 0.707 (thus Ω 1 is the cut-off frequency, but this need not be the case) and B = 0.25. The log-magnitude specification is diagrammed below. Note that (20 log A ) = K 1 dB and (20 log B ) = K 2 dB. Thus the analog filter specifications are
0 20 log 10 H ( j ) K 1 for all Ω Ω 1
20 log 10 H ( j ) K 2 for all Ω Ω 2
With the magnitude |H(jΩ)| given by the Butterworth function,
|H(jΩ)| = 1 / c ^2 N
and using the equality condition at the critical frequencies in the above specifications the order, N , of the filter is given by
Ω , rad./sec.
|H(jΩ)|
A = 0.
Pass band
Transition band
Stop band
B = 0.
dB (= 20 log 10 |H(jΩ)| )
K 1
K 2
N =
2
1 10
/ 10
/ 10 10
2 log
log (^2)
1 K
K
→ (3.16, Ludeman)
The result is rounded to the next larger integer. For example, if N = 3.2 by the above calculation then it is rounded up to 4, and the order of the required filter is N = 4. In such a case the resulting filter would exceed the specification at both Ω 1 and Ω 2. The cut-off frequency Ωc is determined from one of the two equations below.
Ωc = 2 N K / 10
1 10 1 1
or Ωc = 2 N K / 10
2 10 2 1
→ (3.17 Ludeman)
The equation on the left will result in the specification being met exactly at Ω 1 while the specification is exceeded at Ω 2. The equation on the right results in the specification being met exactly at Ω 2 and exceeded at Ω 1. Note that the design equation for N may be written in the alternative form
N =
1
2 10
/ 10
/ 10 10
2 log
log (^1)
2 K
K
Example 5.2.6 What is the order and transfer function of the analog Butterworth filter that satisfies the following specification?
1 = 200 rad/sec K 1 = – 1 dB
2 = 600 rad/sec K 2 = – 30 dB
Solution The order N is given by
N =
2
1 10
/ 10
/ 10 10
2 log
log (^2)
1 K
K
2 log
log
10
( 30 )/ 10
( 1 )/ 10 10 =
2 log
log
10
3
- 1 10
2 log
log
10
10
2 log
log
10
10
2 log
log
10
10
2 log 0. 3333333
log 0. 00025918 10
Now, as in an earlier example, locate 8 poles uniformly on the unit circle, making sure to satisfy all the requirements … and write down the transfer function, H ( s ), of the normalized
Butterworth filter (with a cut-off frequency of 1 rad/sec),
H ( s ) = 1. 848 1 0. 765 1
s^2 s s^2 s
Example 5.3.1 [Low pass filter] For the analog filter Ha(s) = s
A
find the H(z) corresponding
to the impulse invariant design using a sample rate of 1/ T samples/sec.
Solution The analog system’s impulse response is ha ( t )= ℒ–^1
s
A
= Ae ^ tu ( t ). The
corresponding h(n) is then given by
h(n) = ha ( t ) t nT = (^) Ae ^ nTu ( nT )= A e u ( n ) Tn = Aa nu ( n )
where, as previously, we have set e ^ T= a. The discrete-time filter, then, is given by the z - transform of h(n)
H(z) = ʓ{ h(n) } = ʓ A e u ( n ) Tn = T z e
Az (^) = z a
Az
= 1 1 e z
A
(^) T = 1 az 1
A
which has a pole at z = e ^ T= a. In effect, the pole at s = –α in the s -plane is mapped to a pole at
z = e ^ T= a in the z -plane. ( HW What is the difference equation?)
X z
Yz = 1 1 az
A
→ Y ( z )( 1 az ^1 )= AX ( z )
→ y ( n ) ay ( n 1 )= Ax ( n ) → y ( n )= Ax ( n ) ay ( n 1 )
Relationship between the s -plane and the z -plane ( z esT ) We can extend the above procedure
to the case where Ha(s) is given as a sum of N terms with distinct poles as
Ha(s) =
N
k (^) k
k s
A
For this case the impulse invariant design, H(z) , is given by
H(z) = ʓ
t nT
N
k (^) k
k s
A
L
1
1
=
N
k
T
k z e^ k
Az 1
(^) =
N
k
T
k e z
A
11 ^ k^1
where L–^1 means Laplace inverse. We observe that a pole at s = – k in the s -plane gives rise to a
pole at z = e ^ kT in the z -plane and the coefficients in the partial fraction expansion of Ha(s) and H(z) are equal. If the analog filter is stable, corresponding to – k being in the left half plane,
then the magnitude of e ^ kT will be less than unity, so that the corresponding pole of the digital filter is inside the unit circle, and as a result the digital filter also is stable. While the poles in the s -plane “map” to poles in the z -plane according to the relationship
z = esT , it is important to recognize that the impulse invariance design procedure does not correspond to a mapping ( transformation ) of the s -plane to the z -plane by that relationship or in fact by any relationship. (An example of a transformation is where we actually make a
substitution , say, s =
1
1
1
z
z T
which, of course, is the bilinear transformation). For example,
the zeros of Ha(s) do not map to zeros of H(z) according to this relation. See also matched z - transform later.
We can explore the relationship z esT keeping in mind that it only applies to poles and
that it is not a transformation. Set s = ζ + jΩ and z = r ej in z = esT to get r ej = e (^ j ) T =
e T^ ej T so that r = e^ T and ω = ΩT.
The above relations can be used to show that poles in the left half of the primary strip in the s -plane map into poles within the unit circle in the z -plane as shown in the figure for s = s 1.
Mapping of poles, z = esT s -plane pole s = ζ + jΩ
z -plane pole z = esT= r ej r ω
0 1 1 0 jΩs/ 2 – 1 1 π
- ∞ + jΩs/ 2 – 0 0 π
- ∞ – jΩs/ 2 – 0 0 π
For s 1 = ζ 1 + jΩ 1 we have r = e^ ^1 T and ω = Ω 1 T. However, poles at s 2 and s 3 (which are a
distance Ωs from s 1 ) also will be mapped to the same pole that s 1 is mapped to. In fact, an infinite number of s -plane poles will be mapped to the same z -plane pole in a many-to-one relationship. These frequencies differ by Ωs = 2 πFs = 2 π/T ( Fs is the sampling frequency in Hertz). This is called aliasing (of the poles) and is a drawback of the impulse-invariant design. The analog system poles will not be aliased in this manner if, in the first place, they are confined to the “primary strip” of width Ωs = 2πFs = 2π/T in the s -plane. In a similar fashion poles located in the right half of the primary strip in the s -plane will be mapped to the outside of the unit circle in the z -plane. Here again the mapping of the s -plane poles to the z -plane poles is many-to-one.
Owing to the aliasing, the impulse invariant design is suitable for the design of low pass and band pass filters but not for high pass filters.
(Omit) Matched z-transform In this method we apply the mapping z = esT not only to the poles but also to the zeros of Ha(s). As a result the observations made above are valid for the matched z -transform method of filter design.
σ Primary strip, width = Ωs = 2π/T
s -plane
jΩ
- jπ/T
jΩs/2 = jπ/T
j3π/T
- j3π/T
s 2
s 1 = ζ 1 +jΩ 1
s 3
Ωs
z -plane
Im
Re 1
Ω 1 T
Example 5.3.2 [Impulse invariant design of 2nd^ order Butterworth filter] Obtain the impulse
invariant digital filter corresponding to the 2nd^ order Butterworth filter Ha(s) = 2 2 4
s^2 s with sampling time T = 1 sec. Solution Note that we are using T = 1 sec. simply for the purpose of comparing with the bilinear design done later with T = 1 sec. It is important to remember that in bilinear design calculations the value of T is immaterial since it gets cancelled in the design process ; but in impulse invariant design there is no such cancellation, so the value of T is critical (the smaller, the better).
Ha(s) = 2 2 4
(^) s (^2) s = 2 2 2 2 2 2 2
s s
2 2 2 2
s
2 2 2 2
s
The expression in braces is in familiar form and can be converted to its impulse invariant digital filter equivalent. See 3(c) in HW.
Step invariant design Here the response of the digital filter to the unit step sequence, u(n) , is chosen to be samples of the analog step response. In this way, if the analog filter has good step response characteristics, such as small rise-time and low peak over-shoot, these characteristics would be preserved in the digital filter. Clearly this idea of waveform invariance can be extended to the preservation of the output wave shape for a variety of inputs.
(Omit) Problem Given the analog system Ha(s) , let ha(t) be its impulse response and let pa(t) be its step response. The system Ha(s) is given to be continuous-time linear time-invariant. Let also
For T = 0.1, valid for Ω < 10 π
For T = 1, valid for Ω < π /
|Ha(jΩ) |
|Heq(jΩ)|
π
π (^) 10π
10π
pa(t)
t
h(n) be the unit sample response, p(n) be the step response, and, H(z) be the system function,
of a discrete-time linear shift-invariant filter. Then,
(a) If h(n) = ha(nT) , does p(n) =
n
k
ha ( kT )?
(b) If p(n) = pa(nT) , does h(n) = ha(nT)?
Solution (a) If h(n) = ha(nT) , does p(n) =
n
k
ha ( kT )? We know that
u(n) =
n
k
( k )
This is seen to be true by writing it out in full as
u(n) = δ(– ) + δ(– +1) + …+ δ(–1) + δ(0) + δ(1) +…+ δ(n)
where n is implicitly some positive integer. Take, for instance, n = 3; then, from the above equation u(3) = δ(0) = 1, all the other terms being zero. In other words u(n) is a linear combination of unit sample functions. And, since the response to δ(k) is h(k) , therefore, the response to u(n) is a linear combination of the unit sample responses h(k). That is,
p(n) =
n
k
h ( k )=
n
k
ha ( kT )
Therefore, the answer to the above question is, Yes.
(b) If p(n) = pa(nT) , does h(n) = ha(nT)? Since δ(n) = u(n) – u(n–1) , the response of the digital system, H(z) , to the input δ(n) is
h(n) = p(n) – p(n–1) = pa(nT) – pa(nT–T) ha(nT)
Therefore, the answer to the above question is, No.
Example 5.3.3 [LP filter] [Step Invariance] [Problem 5.4b, O & S] Consider the continuous-
time system Ha(s) =
s
A
with unit step response pa(t). Determine the system function, H(z) , i.e.,
the z -transform of the unit sample response h(n) of a discrete-time system designed from this system on the basis of step-invariance, such that p(n) = pa(nT) , where
p(n) =
n
k
h ( k ) and pa(t) =
t
ha ( ) d
Solution Since pa(t) =
t ha ( ) d we have
ℒ[ pa(t) ] = Pa(s) = s
H (^) a ( s )
s ( s )
A
s
K 1
s
K 2
where K 1 and K 2 are the coefficients of the partial fraction expansion, given by
