ELEN E4810: Digital Signal Processing

Topic 8:
Filter Design: IIR
1. Filter Design Specifications
2. Analog Filter Design
3. Digital Filters from Analog Prototypes

Dan Ellis 2013-11-11 1

 1. Filter Design Specifications
 The filter design process:
Problem

Solution
performance G(z)
Analysis constraints Design Implement
transfer
function
• magnitude response • FIR/IIR • platform
• phase response • subtype • structure
• cost/complexity • order • ...

Dan Ellis 2013-11-11 2

 Performance Constraints
 .. in terms of magnitude response:
\G(ejW)\
0
passband
ripple
-10 Passband optimal minimum

transition band
filter stopband
will touch attenuation
-20 here
gain / dB

Stopband
-30

-40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 P
passband stopband frequencyW
edge edge
frequency frequency

Dan Ellis 2013-11-11 3

 Performance Constraints
 “Best” filter:
smallest
Passband Ripple
narrowest
Transition Band
greatest
Minimum SB Attenuation
 improving one usually worsens others
 But: increasing filter order (i.e. cost)
can improve all three measures

Dan Ellis 2013-11-11 4

 Passband Ripple
0
\G(ejW)\

-10

-20

-30

Passband \G(ejW)\
-40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 P

detail
1


√E

PB ripple
parameter
 Assume peak passband gain = 1 1
then minimum passband gain =
1+  2
 Or, ripple  max = 20 log10 1+  2
dB
Dan Ellis 2013-11-11 5
 Stopband Ripple
0
\G(ejW)\

-10

-20

-30

Stopband
-40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 P

detail

1
A

SB ripple 0
parameter freqW
 Peak passband gain is A× larger than
peak stopband gain
 Hence, minimum stopband attenuation

 s = 20 log10 A1 = 20 log10 A dB

Dan Ellis 2013-11-11 6
 Filter Type Choice: FIR vs. IIR
FIR IIR
 No feedback  Feedback
(just zeros) (poles & zeros)
 Always stable  May be unstable

 Can be  Difficult to control

linear phase phase
BUT  High order  Typ. < 1/10th
(20-2000) order of FIR (4-20)
 Unrelated to  Derive from
continuous- analog prototype
time filtering

Dan Ellis 2013-11-11 7

 FIR vs. IIR
 If you care about computational cost
→ use low-complexity IIR
(computation no object → linear phase FIR)

 If you care about phase response

→ use linear-phase FIR
(phase unimportant → go with simple IIR)

Dan Ellis 2013-11-11 8

 IIR Filter Design
 IIR filters are directly related to
analog filters (continuous time)
 via a mapping of H(s) (CT) to H(z) (DT) that
preserves many properties
 Analog filter design is sophisticated
signal processing research since 1940s


→ Design IIR filters via analog prototype

 need to learn some CT filter design

Dan Ellis 2013-11-11 9

 2. Analog Filter Design
 Decades of analysis of transistor-based
filters – sophisticated, well understood
 Basic choices:

 ripples vs. flatness in stop and/or passband

 more ripples → narrower transition band

Family PassBand StopBand

Butterworth flat flat
Chebyshev I ripples flat
Chebyshev II flat ripples
Elliptical ripples ripples
Dan Ellis 2013-11-11 10
 CT Transfer Functions
 Analog systems: s-transform (Laplace)
 frequency response still from a polynomial
Continuous-time Discrete-time

( )  ()
Transform H a s = ha t e dt
st
H d (z ) =  hd [n] z n

Frequency
response
H a ( j) ( )
Hd e j

Im{s} Im{z}
Pole/zero j≠ ej!
diagram Re{s} Re{z}
stable 1
stable s-plane z-plane
poles poles
s scales freq
Dan Ellis 2013-11-11 11
 Butterworth Filters
Maximally flat in pass and stop bands
 Magnitude
1 filter
response (LP): H a ( j) =
2
order
( )
2N
1+ c N

≠ ≪ ≠c,

\Ha(j7)\
1

0.8

|Ha(j≠)|2 →1 0.6

0.4
N
 ≠ = ≠c,
0.2
N
|Ha(j≠)|2 = 1/2 0
0 1 2 3
7/7c
3dB point analog freq

Dan Ellis 2013-11-11 12

 Butterworth Filters 6N dB/oct
rolloff
 ≠ ≫ ≠c, |Ha(j≠)|2 →(≠c/≠)2N

\Ha(j7)\/ dB
0
N
Log-log -10
-20
slope →
magnitude -30
-6N dB/oct
response -40
-50 N
-60
-70
-1 0 1

7/7c
10 10 10

n analog freq
d
H a ( j) = 0
2
 flat → n
d
@ ≠ = 0 for n = 1 .. 2N-1
Dan Ellis 2013-11-11 13
 Butterworth Filters
 How to meet design specifications?
1 1
\Ha(j7)\ 
( ) 1+2
p 2N
 1+
√E
c Design
Equation

( )
1 1
 1
2
A
1 log10  2
1+( )
 s 2N A2 N

A
7p 7s
7
c
( )
2 log10 s
p

 p
 k1 =  k=
A 1
2
s
=“discrimination”, ≪ 1 =“selectivity”, < 1
Dan Ellis 2013-11-11 14
 Butterworth Filters
1
H a ( j) =
2
 ... but what is Ha(s)?
1+ (c )2 N
 Traditionally, look it up in a table
 calculate N → normalized filter with ΩΩc = 1
 scale all coefficients for desired ΩΩ
c
1 Im{s}
 In fact, H a ( s ) = × ΩΩc
i(s  pi ) ×
Re{s}
j N +2 i1 ×  s 2N
where pi =  c e 2N
i = 1..N ×
  = 1
 c 
odd-indexed uniform divisions of c-radius circle
s-plane
Dan Ellis 2013-11-11 15
 Butterworth Example
0 dB
Design a Butterworth 

= -1 dB
√E
filter with 1 dB cutoff
at 1kHz and a
minimum attenuation
 = -40 dB
of 40 dB at 5 kHz A
1kHz 5kHz
7
1 7p 7s
1dB = 20 log10   = 0.259
2
1+  2
9999
40dB = 20 log10 A  A = 100
1 log
N  21 10 0.259
s log10 5
=5
p  N = 4  3.28
Dan Ellis 2013-11-11 16
 Butterworth Example
 Order N = 4 will satisfy constraints;
What are ΩΩc and filter coefficients?
 from a table, ΩΩ-1dB = 0.845 when ΩΩc = 1
ΩΩc = 1000/0.845 = 1.184 kHz
 from a table, get normalized coefficients for
N = 4, scale by 1184·2º

gain / dB
0

-10

Or, use Matlab:

-20
 -30

[B,A] = -40

-50

butter(N,Wc,’s’); -60
0 2000 4000 6000

freq / Hz
Dan Ellis 2013-11-11 17
M
 Chebyshev I Filter (Чебышёв)

 Equiripple in passband (flat in stopband)

→ minimize maximum error
ripple

gain / dB
0

1 depth
|H(j )| = 2 -10
N 
1 + 2 TN2 ( p
) -20

-30 N 
-40
0 0.5 1 1.5
7 2

Chebyshev cos(N cos 1 ) | | 1

polynomial TN ( ) =
of order N cosh(N cosh 1 ) | |>1

Dan Ellis 2013-11-11 18

 1
Chebyshev I Filter |H(j )|2 =
1 + 2 TN2 ( p
)

cos(N cos 1 ) | | 1
TN ( ) =
cosh(N cosh 1 ) | |>1
T (7)
5
2

T5 ( )
1
0

-1

-2

T25(7)
3

T52 ( )
2

-1

1/(1+0.1 T25(7))

1 1

0.5

1 + 0.1T52 ( ) 0

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

7

Dan Ellis 2013-11-11 19

 Chebyshev I Filter
 Design procedure:
 desired passband ripple → "
 min. stopband atten., ΩΩ , ΩΩ → N :
p s

1 1 1
= =
[
A 1+  TN ( p ) 1+  2 cosh N cosh 1  s
( )]
2 2 2 s 2
p

N
cosh 1
( A 2 1
 ) 1/k1, discrimination

( )
cosh-1 grows
1  s 1/k, slower than
cosh  p
selectivity log10

Dan Ellis 2013-11-11 20

 Chebyshev I Filter
 What is Ha(s)?
 complicated, get from a table
 .. or from Matlab cheby1(N,r,Wp,’s’)

 all-pole; can inspect them:

Im{s}
1

0.5

-0.5

-1
-1 -0.5 0 0.5 1
Re{s}
..like squashed-in Butterworth (why?)
Dan Ellis 2013-11-11 21
 Chebyshev II Filter
 Flat in passband, equiripple in stopband

gain / dB
0

1 N 
H a ( j) =
2 -10

-20
N 
constant
T ( s ) 
 2 -30
peak of
stopband
2  N p 
1+ 
-40
ripples
 TN (s ) 
-50

  
-60
~1/TN(1/ΩΩ) 0 0.5 1 1.5
7 2

zeros on imaginary axis

 Filter has poles and zeros (some )

 Complicated pole/zero pattern

Dan Ellis 2013-11-11 22

 Elliptical (Cauer) Filters
 Ripples in both passband and stopband
PB

gain / dB
0
ripple
-10 depth
1
H a ( j)
2
= -20
N  N 
1+  2 RN2 (p ) peak
-30

-40 SB
-50 ripples
function; satisfies -60
0 0.5 1 1.5
7 2

RN(ΩΩ-1) = RN(ΩΩ)-1 very narrow

zeros for ΩΩ<1 → poles for ΩΩ>1 transition band

 Complicated; not even closed form for Ν

Dan Ellis 2013-11-11 23
 Analog Filter Types Summary
gain / dB Butterworth Chebyshev I

gain / dB
0 0

-10 -10 r
-20 -20

-30 -30

-40 -40

-50
c
-50
p
N=6
-60 -60
0 0.5 1 1.5 2 0 0.5 1 1.5 2
r = 3 dB
Chebyshev II Elliptical
A = 40 dB
gain / dB

gain / dB
0 0

-10 -10 r
-20 -20

-30 -30

-40
A -40
A
-50 -50
s p
-60 -60
0 0.5 1 1.5 2 0 0.5 1 1.5 2

Dan Ellis 2013-11-11 24

 Analog Filter Transformations
 All filters types shown as lowpass;
other types (highpass, bandpass..)
derived via transformations
 i.e. s
ˆ = F 1
(s ) Desired alternate
response; still a
lowpass
prototype H LP (s )  H D (sˆ ) rational polynomial

 General mapping of s-plane

^
BUT keep LHHP & jΩΩ → jΩΩ;
poles OK, frequency response ‘shuffled’
Dan Ellis 2013-11-11 25
 Lowpass-to-Highpass
 Example transformation:
H HP (sˆ) = H LP ( s ) s=  p ˆ p
sˆ
 take prototype HLP(s) polynomial
 replace s with  p ˆp
sˆ
 simplify and rearrange
^
→ new polynomial HHP(s)

Dan Ellis 2013-11-11 26

 Lowpass-to-Highpass
 What happens to frequency response?
s = j  sˆ =
 p 
ˆp
 p
ˆp
j
=j (  p 

ˆp
) imaginary axis
stays on self...

ˆ =
 ...freq.→freq. ^
7
 = p  
ˆ =  HPF
 ˆp
 < p  
ˆ < 
ˆp LPF
7
LP passband HP passband 7p

 > p  
ˆ > 
ˆp ^
7p

LP stopband HP stopband
 Frequency axes inverted
Dan Ellis 2013-11-11 27
 Transformation Example
Design a Butterworth highpass filter
^
with PB edge -0.1dB @ 4 kHz (ΩΩp)
^
and SB edge -40 dB @ 1 kHz (ΩΩs)
 Lowpass prototype: make ΩΩ = 1
p
 p p
ˆ
  s = () ˆ = ()4
s
 Butterworth -0.1dB @ ΩΩ =1, -40dB @ ΩΩ =4
p s
1 log10 (  ) 1
A 12
0.1
N 
2
@ 0.1dB  = 10 10

2 log10 (  )  p 10
p

s
p
1+ ( c )
➝ N=5, 2N=10  c =  p /0.6866 = 1.4564
Dan Ellis 2013-11-11 28
 Transformation Example × ΩΩc
Im{s}

j N +2 1
×
 LPF proto has p = c e 2N
Re{s}
×
c N
 H LP ( s ) = N ×
=1 ( s  p )
 Map to HPF: H HP (sˆ) = H LP ( s ) s=  p ˆ p
sˆ
N zeros
cN cN sˆ N @ ^s = 0
 H HP (sˆ) = =
 =1(
N  p
sˆ
ˆp
 p ) =1 ( pˆ p  p sˆ)
N

^ /p
new poles @ s^ = ΩΩpΩΩ p l

Dan Ellis 2013-11-11 29

 Transformation Example
In Matlab:
ΩΩp ΩΩs Rp Rs
[N,Wc]=buttord(1,4,0.1,40,'s');
[B,A] = butter(N, Wc, 's');
[U,D] = lp2hp(B,A,2*pi*4000);
4
x 10

gain / dB
1.5 0
1 -10
Im{s}

0.5
-20
0
-30
-0.5
-40
-1
-50
-1.5
-60
-2 -1 0 1 0 1000 2000 3000 4000 5000
7 / Hz
4
Re{s} x 10

Dan Ellis 2013-11-11 30

 3. Analog Protos → IIR Filters
 Can we map high-performance CT
filters to DT domain?
 Approach: transformation H (s)→G(z)
a
i.e. G (z ) = H a (s ) s=F ( z )
where s = F(z) maps s-plane z-plane:
Im{s} s = F(z) Im{z}
Ha(s0) z = F-1(s) G(z0)
Every value of G(z) Re{s} Re{z}
1
is a value of Ha(s)
somewhere on the s-plane z-plane
s-plane & vice-versa
Dan Ellis 2013-11-11 31
 CT to DT Transformation
 Desired properties for s = F(z):
 s-plane jΩΩ axis z-plane unit circle
→ preserves frequency response values
 s-plane LHHP z-plane unit circle interior
→ preserves stability of poles
Im{s} Im{z}
Im u.c. ej!
jΩΩ
Re{s} Re{z}
1

s-plane LHHP UCI z-plane

Dan Ellis 2013-11-11 32

 Bilinear Transformation
 Solution: 1 z 1
z 1 Bilinear
s= =
1+z 1 z + 1 Transform
1+ s unique,
 Hence inverse: z =
1 s 1:1 mapping

|z| = 1 i.e.
 Freq. axis? s = j  z = 1+ j
1 j on unit circle
(1+ )+ j
 Poles? s =  + j  z = (1 ) j
2 1+ 2 +  2
+ 2
z = σ<0
1  2 +  + 
2 2
|z| < 1
Dan Ellis 2013-11-11 33
 Bilinear Transformation
 How can entire half-plane fit inside u.c.?
4 4

3 3

2 2

1 1

0 0

-1 -1

-2 -2

-3 -3

s-plane
-4
-4 -2 0 2 4
-4
-4 -2 0
z-plane
2 4

 Highly nonuniform warping!

 “Moebius Transformations Revealed”
[Link]

Dan Ellis 2013-11-11 34

 Bilinear Transformation
 What is CT DT freq. relation ΩΩ ! ?
j  j 2 j sin  / 2
z = e  s = 1+e j = 2 cos / 2 = j tan 2
1e
[Link] [Link]
W (DT)
P

 = tan ( 
)
0.5P

 i.e. 2 0 7
(CT)
 = 2 tan  1 -0.5P

-P
-5 0 5

 infinite range of CT frequency  <  < 

maps to finite DT freq. range  <  < 
 nonlinear; d    as ! → º pack it all in!
d
Dan Ellis 2013-11-11 35
 Frequency Warping
 Bilinear transform makes
Ge ( )
j
= H a ( j)  =2 tan1  for all !, ΩΩ
|G(ej )|
0.8

0.6
= 2tan 1( )  Same gain &
0.4 phase (", A...),
0.2
in same ‘order’,
0
-60 -40 -20 0

0
but with
-20
|Ha(j )| warped
-40
frequency axis
-60
0 1 2 3 4 5 6

Dan Ellis 2013-11-11 36

 Design Procedure
 Obtain DT filter specs:
 general form (LP, HP...),  p ,  s , 1 ,1
1+ 2 A
 ‘Warp’ frequencies to CT:
p s

 p = tan 2  s = tan 2
Old-  Design analog filter for  , , 1 ,1
p s 1+ 2 A
style  → Ha(s), CT filter polynomial
 Convert to DT domain: G (z ) = H a (s ) s=1z1
1+z 1
 → G(z), rational polynomial in z
 Implement digital filter!
Dan Ellis 2013-11-11 37
 Bilinear Transform Example
 DT domain requirements:
Lowpass, 1 dB ripple in PB, !p = 0.4º,
SB attenuation ≥ 40 dB @ !s = 0.5º,
attenuation increases with frequency
|G(ejW)| / dB

0
-1
PB ripples,
SB monotonic
→ Chebyshev I
-40

WP
0.4 0.5 1

Dan Ellis 2013-11-11 38

 Bilinear Transform Example
 Warp to CT domain:
p
 p = tan 2 = tan 0.2 = 0.7265 rad/sec
s
 s = tan 2 = tan 0.25 = 1.0 rad/sec
 Magnitude specs:
1 dB PB ripple
1/ 20
 1 = 10 = 0.8913   = 0.5087
2
1+

40 dB SB atten.
40 / 20
 A = 10
1 = 0.01  A = 100
Dan Ellis 2013-11-11 39
Bilinear Transform Example
 Chebyshev I design criteria:

N
cosh (
1 A 2 1

=
)
7.09 i.e. need N = 8
1  s
cosh  p ( )
 Design analog filter, map to DT, check:

\Ha(j7)\/ dB

\G(ejW)\/ dB
0 0

-10 -10

>> N=8; -20

-30
-20

-30

>> wp=0.7265; CT DT
-40

-50
-40

-50
freqs(B,A) freqz(b,a)
>> pbripple = 1.0; 7
-60
0
W/P0.5 1 1.5 2
-60
0 0.2 0.4 0.6 0.8 1

>> [B,A] = cheby1(N,pbripple,wp,'s');

>> [b,a] = bilinear(B,A,.5);
Dan Ellis 2013-11-11 40
M
 Other Filter Shapes
 Example was IIR LPF from LP prototype
 For other shapes (HPF, bandpass,...):

CT Bilinear
Bilinear Analog trans HD(s) transform
warp design
DT CT HLP(s)
specs specs GD(z)
Bilinear
transform GLP(z) DT
trans
 Transform LP→X in CT or DT domain...

Dan Ellis 2013-11-11 41

 DT Spectral Transformations
 Same idea as CT LPF→HPF mapping,
but in z-domain:
GD (zˆ ) = GL (z ) z=F ( zˆ ) = GL (F (zˆ))
 To behave well, z = F (zˆ ) should:
 map u.c. → u.c. (preserve G(ej!) values)
 map u.c. interior → u.c. interior (stability)
 i.e. F (zˆ) = 1 zˆ = 1 F (zˆ) < 1 zˆ < 1
 in fact, F ( z
ˆ ) matches the definition of an
1
allpass filter ... replace delays with F ( zˆ )
Dan Ellis 2013-11-11 42
 DT Frequency Warping
 Simplest mapping z = F (zˆ ) = 1
zˆ 
zˆ
has effect of warping frequency axis:
jˆ
e
zˆ = e jˆ
z=e =j
1 ejˆ Æ>0:
expand HF
P
W
0.8
A = 0.6

 tan (2 ) = 1 (2 )
0.6
1+ tan ˆ
 0.4
Æ<0:
A = -0.6 expand LF
0.2

0
0 0.2 0.4 0.6 0.8 P
^
W
Dan Ellis 2013-11-11 43
 Another Design Example
 Spec:
 Bandpass, from 800-1600 Hz (SR = 8kHz)
 Ripple = 1dB, min. stopband atten. = 60 dB

 8th order, best transition band

→ Use elliptical for best performance

 Full design path:
 prewarp spec frequencies
 design analog LPF prototype

 analog LPF → BPF

 CT BPF → DT BPF (Bilinear)

Dan Ellis 2013-11-11 44

 Another Design Example
 Frequency prewarping
 800 Hz @ SR = 8kHz
= 800/8000 x 2  rad/sample
→ L =0.2 , U =1600/8000 x 2  =0.4 
L = tan L/2 = 0.3249, U = 0.7265


 Note distinction between:

 application’s CT domain (800-1600 Hz)

and
 CT prototype domain (0.3249-0.7265 rad/s)

Dan Ellis 2013-11-11 45

 Another Design Example
Or, do it all in one step in Matlab:
[b,a] = ellip(8,1,60,


[800 1600]/(8000/2));
0

-20
Magnitude (dB)

-40

-60

-80
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized Frequency (¾P rad/sample)

1000

500
Phase (degrees)

-500

-1000
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized Frequency (¾P rad/sample)

Dan Ellis 2013-11-11 46