IIR Filters

Analog Lowpass Filter
Specifications

Specifications
Typical magnitude response H a ( j) of an

analog lowpass filter may be given as
indicated below
In the passband, defined by 0 p , we

require
1 p H a ( j ) 1 + p , p
i.e., H a ( j) approximates unity within an
error of p
In the stopband, defined by s < , we
require
H a ( j ) s , s <
2
i.e., H a ( j) approximates zero within an

error of s
Copyright 2005, S. K. Mitra

Specifications

Specifications
p - passband edge frequency

s - stopband edge frequency
p - peak ripple value in the passband
s - peak ripple value in the stopband
Peak passband ripple
p = 20 log10 (1 p ) dB
Magnitude specifications may alternately be

given in a normalized form as indicated
below
Minimum stopband attenuation

s = 20 log10 ( s ) dB
3
4

Specifications
Analog Lowpass Filter Design

Two additional parameters are defined -
Here, the maximum value of the magnitude

in the passband assumed to be unity
(1) Transition ratio k =
1 / 1 + 2 - Maximum passband deviation,

given by the minimum value of the
magnitude in the passband
For a lowpass filter k < 1

(2) Discrimination parameter k1 =
Usually k1 << 1
1 - Maximum stopband magnitude

A
5
p
s
A2 1
6
Butterworth Approximation
The magnitude-square response of an N-th

order analog lowpass Butterworth filter
is given by
1
2
H a ( j ) =
1 + ( / c ) 2 N
2
First 2 N 1 derivatives of H a ( j) at = 0
are equal to zero
The Butterworth lowpass filter thus is said
to have a maximally-flat magnitude at = 0
Gain in dB is G () = 10 log10 H a ( j) 2
As G (0) = 0 and
G (c ) = 10 log10 (0.5) = 3.0103 3 dB
c is called the 3-dB cutoff frequency
8
Typical magnitude responses with c = 1
Two parameters completely characterizing a

Butterworth lowpass filter are c and N
These are determined from the specified
bandedges p and s , and minimum
passband magnitude 1 / 1 + 2, and
maximum stopband ripple 1 / A
Butterworth Filter
N=2
N=4
N = 10
Magnitude
1
0.8
0.6
0.4
0.2
0
10
c and N are thus determined from

2
1
= 1
H a ( j p ) =
1 + ( p / c ) 2 N 1 + 2
2
1
= 1
H a ( j s ) =
1 + ( s / c ) 2 N A2
Solving the above we get
log [( A2 1) / 2 ] log10 (1 / k1 )
N = 1 10
=
2 log10 ( s / p )
log10 (1 / k )
11
12
Since order N must be an integer, value

obtained is rounded up to the next highest
integer
This value of N is used next to determine c
by satisfying either the stopband edge or the
passband edge specification exactly
If the stopband edge specification is
satisfied, then the passband edge
specification is exceeded providing a safety
margin
Transfer function of an analog Butterworth
lowpass filter is given by
cN
cN
H a ( s) = C = N
=
DN ( s ) s + lN=01d s l lN=1 ( s pl )
l
Example - Determine the lowest order of a
Butterworth lowpass filter with a 1-dB cutoff
frequency at 1 kHz and a minimum attenuation of 40
dB at 5 kHz
Now
10 log10 1 2 = 1
1 +
which yields 2 = 0.25895
and
10 log10 12 = 40
A
which yields A2 = 10,000
where
pl = c e j[ ( N + 2l1) / 2 N ] , 1 l N
Denominator DN (s ) is known as the
Butterworth polynomial of order N
13
14
Chebyshev Approximation
A2 1 = 196.51334
1=
k1
Therefore
1 = s = 5
k p
and

order analog lowpass Type 1 Chebyshev filter
is given by
1
2
H a ( s) =
2 2
1 + TN ( / p )
Hence
where TN () is the Chebyshev polynomial

of order N:
log (1 / k1 )
N = 10
= 3.2811
log10 (1 / k )
cos( N cos 1 ),
TN () =
1
cosh( N cosh ),
We choose N = 4
15
1
>1
16
If at = s the magnitude is equal to 1/A,

then
2
1
= 12
H a ( j s ) =
2 2
1 + TN ( s / p ) A
Typical magnitude response plots of the

analog lowpass Type 1 Chebyshev filter are
shown below
Type 1 Chebyshev Filter
Magnitude
Solving the above we get

cosh 1 ( A2 1 / ) cosh 1 (1 / k1 )
=
N=
cosh 1 ( s / p )
cosh 1 (1 / k )
Order N is chosen as the nearest integer
greater than or equal to the above value
N=2
N=3
N=8
1
0.8
0.6
0.4
0.2
0
17
18

order analog lowpass Type 2 Chebyshev
(also called inverse Chebyshev) filter is
given by
1
2
H a ( j ) =
2
T ( / p )
1+ 2 N s
TN ( s / )
Typical magnitude response plots of the

analog lowpass Type 2 Chebyshev filter are
shown below
Type 2 Chebyshev Filter
Magnitude
0.6
0.4
0.2
where TN () is the Chebyshev polynomial

of order N
19
N=3
N=5
N=7
1
0.8
20
Elliptic Approximation
The square-magnitude response of an

elliptic lowpass filter is given by
1
2
H a ( j ) =
2 2
1 + RN ( / p )
where RN () is a rational function of order
N satisfying RN (1 / ) = 1 / RN () , with the
roots of its numerator lying in the interval
0 < < 1 and the roots of its denominator
lying in the interval 1 < <
The order N of the Type 2 Chebyshev filter

is determined from given , s , and A
using
cosh 1 ( A2 1 / ) cosh 1 (1 / k1 )
=
N=
cosh 1 ( s / p )
cosh 1 (1 / k )
Example - Determine the lowest order of a
Chebyshev lowpass filter with a 1-dB cutoff
frequency at 1 kHz and a minimum attenuation of
40 dB at 5 kHz -
N=
21
cosh 1 (1 / k1 )
= 2.6059
cosh 1 (1 / k )
22
For given p , s , , and A, the filter order

can be estimated using
2 log10 (4 / k1 )
N
log10 (1 / )
Example - Determine the lowest order of a elliptic

lowpass filter with a 1-dB cutoff frequency at 1
kHz and a minimum attenuation of 40 dB at 5 kHz
Note: k = 0.2 and 1 / k1 = 196.5134

Substituting these values we get
0 = 0.00255135,
k '= 0.979796,
where k ' = 1 k 2
0 = 1 k '
2(1 + k ')
= 0 + 2( 0 )5 + 15( 0 )9 + 150( 0 )13
23
= 0.0025513525
and hence N = 2.23308
Choose N = 3
24
Example - Design an elliptic lowpass filter

of lowest order with a 1-dB cutoff
frequency at 1 kHz and a minimum
attenuation of 40 dB at 5 kHz
Code fragments used
Typical magnitude response plots with p = 1

are shown below
Elliptic Filter
N=3
N=4
Magnitude
[N, Wn] = ellipord(Wp, Ws, Rp, Rs, s);
0.8
0.6
0.2
0
25
26
Design of Analog Highpass,

Highpass,
Bandpass and Bandstop Filters

Gain plot
Lowpass Elliptic Filter
Gain, dB
0
-20
-40
-60
[b, a] = ellip(N, Rp, Rs, Wn, s);

Wp = 2*pi*1000;
Ws = 2*pi*5000;
Rp = 1;
Rs = 40;
with
0.4
2000
4000
Frequency, Hz
6000
27
28
Design of Analog Highpass,

Highpass,
Bandpass and Bandstop Filters
Analog Highpass Filter Design

Spectral Transformation:
p
p
s=
s
where p is the passband edge frequency of
p is the passband edge
H LP ( s) and
frequency of H HP ( s)
On the imaginary axis the transformation is
p
p
=
Let s denote the Laplace transform variable

of prototype analog lowpass filter H LP (s)
and s denote the Laplace transform
variable of desired analog filter HD (s)
The mapping from s-domain to s -domain is
given by the invertible transformation
s = F (s)
H
(
s
)
=
H LP ( s ) s = F ( s)
Then
D
H LP ( s ) = HD ( s) s = F 1 ( s )
29
Steps involved in the design process:

Step 1 - Develop of specifications of a
prototype analog lowpass filter H LP (s)
from specifications of desired analog filter
HD (s ) using a frequency transformation
Step 2 - Design the prototype analog
lowpass filter
Step 3 - Determine the transfer function HD (s )
of desired analog filter by applying the
inverse frequency transformation to H LP (s)
30

=
p
p
p s
Stopband s p
Stopband
0
Passband
Stopband
Passband
p
s
Passband
Lowpass
Highpass
31
Example - Design an analog Butterworth

highpass filter with the specifications:
F p = 4 kHz, Fs = 1 kHz, p = 0.1 dB,
s = 40 dB
Choose p = 1
Then
2F p F p 4000
=
=
=4
s =
2Fs Fs 1000
Analog lowpass filter specifications: p = 1,
s = 4 , p = 0.1 dB, s = 40 dB
32
Analog Bandpass Filter

Design

Code fragments used
[N, Wn] = buttord(1, 4, 0.1, 40, s);
[B, A] = butter(N, Wn, s);
[num, den] = lp2hp(B, A, 2*pi*4000);
Spectral Transformation
o2
s 2 +
s = p
p2
p1 )
s(
Gain plots
-20
-20
-40
-60
33
-80
0
where p is the passband edge frequency

p1 and
p 2 are the lower
of H LP (s ), and
and upper passband edge frequencies of
desired bandpass filter H BP (s)
Highpass Filter
0
Gain, dB
Gain, dB
Prototype Lowpass Filter

0
-40
-60
10
-80
0
4
6
8
10
Frequency, kHz
34

Design
2
2
= p o
Bw
p2
p1 is the width of
where Bw =
passband and o is the passband center

frequency of the bandpass filter
Passband edge frequency p is mapped
p1 and
p 2, lower and upper
into m
passband edge frequencies
35

Design
= p
o2
2
Bw
p s
Stopband s p
Stopband
0
Passband
Stopband
Passband
Stopband
Passband
Stopband

s1 0
o
s1
s2
s2
p2 o
p1
p1
p2
Lowpass
Bandpass
36

Design

Design
Stopband edge frequency s is mapped

s 2 , lower and upper
s1 and
into m
stopband edge frequencies
Also,
o2 =
p1
p2 =
s1
s2
p1
p2 >
s1
s2
Case 1:
s 2 we can either
To make p1 p 2 = s1
increase any one of the stopband edges or
decrease any one of the passband edges as
shown below
If bandedge frequencies do not satisfy the

above condition, then one of the frequencies
needs to be changed to a new value so that
the condition is satisfied
37
Passband
Stopband
s1
Stopband
p1
p2
s2
38

Design

Design
s1
s2 /
p2
p1 to
(1) Decrease
larger passband and shorter
leftmost transition band
p1
p2 /
s2
s1 to
(2) Increase
o2 =
p1
p2 =
s1
s2
Note: The condition
p2
can also be satisfied by decreasing
which is not acceptable as the passband is
reduced from the desired value
Alternately, the condition can be satisfied
s 2 which is not acceptable
by increasing
as the upper stop band is reduced from the
desired value
No change in passband and shorter

39
40

Design

Design
p1
p2 <
s1
s2
Case 2:
p1
p2 =
s1
s 2 we can either
To make
decrease any one of the stopband edges or
increase any one of the passband edges as
shown below
Passband
Stopband
s1
p1
p 2 to
s1
s2 /
p1
(1) Increase
rightmost transition band
p1
p2 /
s1
s 2 to
(2) Decrease
No change in passband and shorter
rightmost transition band
p2
s2
Stopband
41
42

Design

Design
Example - Design an analog elliptic

bandpass filter with the specifications:
F p1 = 4 kHz, F p 2 = 7 kHz, Fs1 = 3 kHz
Fs 2 = 8 kHz, p = 1 dB, s = 22 dB
o2 =
p1
p2 =
s1
s2
Note: The condition
p1
can also be satisfied by increasing
which is not acceptable as the passband is
reduced from the desired value
s1 which is not acceptable
by decreasing
as the lower stopband is reduced from the
desired value
43
Now F p1F p 2 = 28 106 and Fs1Fs 2 = 24 106

Since F p1F p 2 > Fs1Fs 2 we choose
F p1 = Fs1Fs 2 / F p 2 = 3.571428 kHz
44

Design

Design
Code fragments used
[N, Wn] = ellipord(1, 1.4, 1, 22, s);

[B, A] = ellip(N, 1, 22, Wn, s);
[num, den]
= lp2bp(B, A, 2*pi*4.8989795, 2*pi*25/7);
We choose p = 1
Hence
24 9
s =
= 1.4
(25 / 7) 3
Gain plot
Prototype Lowpass Filter

s = 1.4 , p = 1 dB, s = 22 dB
45
0
Gain, dB
Gain, dB
46
-20
-40
-60
0
-20
-40
-60
0
5
10
15
Frequency,
kHz 2005, S. K. Mitra
Copyright
Analog Bandstop Filter Design

B
= s 2 w 2
o
where Bw = s 2 s1 is the width of

o is the stopband center
stopband and
frequency of the bandstop filter
Stopband edge frequency s is mapped
s1and
s 2 , lower and upper
into m
stopband edge frequencies
Spectral Transformation
s2
s1 )
s(
s = s
2
2o
s +
where s is the stopband edge frequency
s1 and
s 2 are the lower
of H LP ( s ) , and
and upper stopband edge frequencies of the
desired bandstop filter H BS ( s)
47
Bandpass Filter
48
Passband edge frequency p is mapped

p1 and
p 2 , lower and upper
into m
passband edge frequencies
Also,
p s
Stopband s p
Stopband
0
Passband
Passband
Stopband
Passband
Stopband
Passband
o
p1
o
p10
p2
p2
s2
s1
s1
s2
o2 =
p1
p2 =
s1
s2
If bandedge frequencies do not satisfy the

above condition, then one of the frequencies
needs to be changed to a new value so that
the condition is satisfied
Lowpass
Bandpass
49
50
p1
p2 >
s1
s2
Case 1:
p1
p2 =
s1
s 2 we can either
To make
increase any one of the stopband edges or
decrease any one of the passband edges as
shown below
Passband
Passband
Stopband
51
p 2 to
s1
s2 /
p2
(1) Decrease
larger high-frequency passband
and shorter rightmost transition band
p1
p2 /
s2
s 2 to
(2) Increase
No change in passbands and
shorter rightmost transition band
p1
s1
s2
p2
52
o2 =
p1
p2 =
s1
s2
Note: The condition
can also be satisfied by decreasing
p1
which is not acceptable as the low-
frequency passband is reduced from the
desired value
s1 which is not acceptable
by increasing
as the stopband is reduced from the desired
value
p1
p2 <
s1
s2
Case 1:
s 2 we can either
To make p1 p 2 = s1
decrease any one of the stopband edges or
increase any one of the passband edges as
shown below
53
Passband
54
Passband
Stopband
p1
s1
s2
p2

o2 =
p1
p2 =
s1
s2
Note: The condition
p2
can also be satisfied by increasing
which is not acceptable as the highfrequency passband is decreased from the
desired value
s 2 which is not acceptable
by decreasing
as the stopband is decreased
p1 to
s1
s2 /
p1
(1) Increase
s1 to
p1
p2 /
s1
(2) Decrease
No change in passbands and
shorter leftmost transition band
55
56
10
Digital Filter Design
Digital Filter Specifications

Usually, either the magnitude and/or the
phase (delay) response is specified for the
design of digital filter for most applications
In some situations, the unit sample response
or the step response may be specified
In most practical applications, the problem
of interest is the development of a realizable
approximation to a given magnitude
response specification
Objective - Determination of a realizable

transfer function G(z) approximating a
given frequency response specification is an
important step in the development of a
digital filter
If an IIR filter is desired, G(z) should be a
stable real rational function
Digital filter design is the process of
deriving the transfer function G(z)
1
We discuss in this course only the

magnitude approximation problem
There are four basic types of ideal filters
with magnitude responses as shown below
HLP(e j )
HHP (e j )
c 0
c2 c1
HBS (e j )
HBP (e j )
c1
c2
c2 c1
c1
c2
As the impulse response corresponding to

each of these ideal filters is noncausal and
of infinite length, these filters are not
realizable
In practice, the magnitude response
specifications of a digital filter in the
passband and in the stopband are given with
some acceptable tolerances
In addition, a transition band is specified
between the passband and stopband

j
As indicated in the figure, in the passband,

defined by 0 p , we require that
G (e j ) 1 with an error p, i.e.,
For example, the magnitude response G (e )

of a digital lowpass filter may be given as
indicated below
1 p G ( e j ) 1 + p ,
In the stopband, defined by s , we

require that G (e j ) 0 with an error s ,
i.e.,
G ( e j ) s , s
5
p - passband edge frequency

s - stopband edge frequency
p - peak ripple value in the passband
s - peak ripple value in the stopband
Since G (e j ) is a periodic function of ,
and G (e j ) of a real-coefficient digital
filter is an even function of
As a result, filter specifications are given
only for the frequency range 0

Specifications are often given in terms of
j
loss function A () = 20 log10 G (e ) in
dB
Peak passband ripple
p = 20 log10 (1 p ) dB
Minimum stopband attenuation
s = 20 log10 ( s ) dB
8
Magnitude specifications may alternately be

given in a normalized form as indicated
below
Here, the maximum value of the magnitude

in the passband is assumed to be unity
1 / 1 + 2 - Maximum passband deviation,
given by the minimum value of the
magnitude in the passband
1 - Maximum stopband magnitude
A

For the normalized specification, maximum
value of the gain function or the minimum
value of the loss function is 0 dB
Maximum passband attenuation max = 20 log10 1 + 2 dB
For p << 1, it can be shown that
max 20 log10 (1 2 p ) dB
11
10

In practice, passband edge frequency Fp
and stopband edge frequency Fs are
specified in Hz
For digital filter design, normalized
bandedge frequencies need to be computed
from specifications in Hz using
p 2 Fp
p =
=
= 2 FpT
FT
FT
2 Fs
s = s =
= 2 Fs T
FT
FT
12

Example - Let Fp = 7 kHz, Fs = 3 kHz, and
FT = 25 kHz
Then
2(7 103 )
p =
= 0.56
25 103
s =
2(3 103 )
= 0.24
25 103
13
Selection of Filter Type

The transfer function H(z) meeting the
frequency response specifications should be
a causal transfer function
For IIR digital filter design, the IIR transfer
function is a real rational function of z 1:
p0 + p1z 1 + p2 z 2 + L + pM z M
, MN
d 0 + d1z 1 + d 2 z 2 + L + d N z N
H(z) must be a stable transfer function and
must be of lowest order N for reduced
14 computational complexity
H ( z) =

For FIR digital filter design, the FIR
transfer function is a polynomial in z 1
with real coefficients:
Advantages in using an FIR filter (1) Can be designed with exact linear phase,
(2) Filter structure always stable with
quantized coefficients
Disadvantages in using an FIR filter - Order
of an FIR filter, in most cases, is
considerably higher than the order of an
equivalent IIR filter meeting the same
specifications, and FIR filter has thus higher
computational complexity
H ( z ) = h[n] z n
n =0
For reduced computational complexity,

degree N of H(z) must be as small as
possible
If a linear phase is desired, the filter
coefficients must satisfy the constraint:
h[n] = h[ N n]
15
16
Digital Filter Design:

Basic Approaches

Basic Approaches
Most common approach to IIR filter design (1) Convert the digital filter specifications
into an analog prototype lowpass filter
specifications
(2) Determine the analog lowpass filter
transfer function H a (s )
(3) Transform H a (s ) into the desired digital
transfer function G ( z )
17
18
This approach has been widely used for the

following reasons:
(1) Analog approximation techniques are
highly advanced
(2) They usually yield closed-form
solutions
(3) Extensive tables are available for
analog filter design
(4) Many applications require digital
simulation of analog systems

Basic Approaches

Basic Approaches
An analog transfer function to be denoted as

P ( s)
H a ( s) = a
Da ( s )
Basic idea behind the conversion of H a (s )

into G ( z ) is to apply a mapping from the
s-domain to the z-domain so that essential
properties of the analog frequency response
are preserved
Thus mapping function should be such that
Imaginary ( j ) axis in the s-plane be
where the subscript a specifically

indicates the analog domain
A digital transfer function derived from H a ( s )
shall be denoted as
P( z )
G( z) =
D( z )
19
20

Basic Approaches
Three commonly used approaches to FIR

filter design (1) Windowed Fourier series approach
(2) Frequency sampling approach
(3) Computer-based optimization methods
22
IIR Digital Filter Design: Bilinear

Transformation Method
Digital filter design consists of 3 steps:

(1) Develop the specifications of Ha ( s ) by
applying the inverse bilinear transformation
to specifications of G(z)
(2) Design Ha (s )
(3) Determine G(z) by applying bilinear
transformation to Ha ( s )
As a result, the parameter T has no effect on
G(z) and T = 2 is chosen for convenience
Above transformation maps a single point

in the s-plane to a unique point in the
z-plane and vice-versa
Relation between G(z) and Ha ( s ) is then
given by
G ( z ) = Ha ( s ) 2 1 z 1
s=
T 1+ z 1
Bilinear Transformation
Bilinear transformation 1
s = 2 1 z 1
T 1 + z
23

Basic Approaches
FIR filter design is based on a direct

approximation of the specified magnitude
response, with the often added requirement
that the phase be linear
The design of an FIR filter of order N may
be accomplished by finding either the
length-(N+1) impulse response samples {h[n]}
or the (N+1) samples of its frequency
response H (e j )
21
mapped onto the unit circle of the z-plane

A stable analog transfer function be mapped
into a stable digital transfer function
24
Inverse bilinear transformation for T = 2 is

z = 1+ s
1 s
Mapping of s-plane into the z-plane
For s = o + jo
2
2
(1 + o ) + jo
2 (1 + o ) + o
z=
z =
(1 o ) jo
(1 o ) 2 + o2
o = 0 z = 1
Thus,
o < 0 z < 1
o > 0 z > 1
25
26
For z = e with T = 2 we have

j e j / 2 (e j / 2 e j / 2 )
j = 1 e j = j / 2 j / 2 j / 2
1+ e
e
(e
+e
)
j 2 sin( / 2)
=
= j tan( / 2)
2 cos( / 2)
or = tan( / 2)
27
Mapping is highly nonlinear

Complete negative imaginary axis in the splane from = to = 0 is mapped into
the lower half of the unit circle in the z-plane
from z = 1 to z = 1
Complete positive imaginary axis in the splane from = 0 to = is mapped into the
upper half of the unit circle in the z-plane
from z = 1 to z = 1
28
29
= tan(/2)
Nonlinear mapping introduces a distortion
in the frequency axis called frequency
warping
Effect of warping shown below
Steps in the design of a digital filter (1) Prewarp ( p , s ) to find their analog
equivalents ( p , s )
(2) Design the analog filter H a (s )
(3) Design the digital filter G(z) by applying
bilinear transformation to H a (s )
Transformation can be used only to design
digital filters with prescribed magnitude
response with piecewise constant values
Transformation does not preserve phase
response of analog filter
30
IIR Digital Filter Design Using


Example - Consider
c
s + c
Applying bilinear transformation to the above
we get the transfer function of a first-order
digital lowpass Butterworth filter
Rearranging terms we get
H a ( s) =
1+ z 1
where
c (1 + z 1 )
=
(1 z 1 ) + c (1 + z 1 )
31
32


Example - Consider the second-order analog

notch transfer function
s 2 + o2
H a ( s) = 2
s + B s + o2
for which H a ( jo ) = 0
33
H a ( j 0) = H a ( j) = 1
o is called the notch frequency
If H a ( j 2 ) = H a ( j1 ) = 1 / 2 then
B = 2 1 is the 3-dB notch bandwidth
Then
=
35
1+
1 2 z + z
2 1 2 (1 + ) z 1 + z 2
1 + o2 B 1 tan( Bw / 2)
=
1 + o2 + B 1 + tan( Bw / 2)
1 o2
=
= cos o
1 + o2
34

Thus
0.954965 1.1226287 z 1 + 0.954965 z 2

1 1.1226287 z 1 + 0.90993 z 2
The gain and phase responses are shown below
G( z) =
2
0
-10
-20
-30
= 0.587785
where =
Gain, dB
= 0.90993
G ( z ) = H a ( s ) s =1 z
1+ z 1
2
2
(1 + o ) 2(1 o ) z 1 + (1 + o2 ) z 2
(1 + o2 + B) 2(1 o2 ) z 1 + (1 + o2 B ) z 2
1
2

Example - Design a 2nd-order digital notch
filter operating at a sampling rate of 400 Hz
with a notch frequency at 60 Hz, 3-dB notch
bandwidth of 6 Hz
Thus o = 2(60 / 400) = 0.3
Bw = 2(6 / 400) = 0.03
From the above values we get
1 c 1 tan(c / 2)
=
1 + c 1 + tan(c / 2)
Phase, radians
G ( z ) = H a ( s ) s =1 z 1
1
G ( z ) = 1 1 + z 1
2 1 z
36
-40
0
50
100
150
Frequency, Hz
200
1
0
-1
-2
0
50
100
150
Frequency, Hz
200
IIR Lowpass Digital Filter Design

Using Bilinear Transformation

Example - Design a lowpass Butterworth

digital filter with p = 0.25, s = 0.55 ,
p 0.5 dB, and s 15 dB
Thus
A2 = 31.622777
2 = 0.1220185
Prewarping we get
p = tan( p / 2) = tan(0.25 / 2) = 0.4142136
s = tan(s / 2) = tan(0.55 / 2) = 1.1708496
The inverse transition ratio is
1 s
=
= 2.8266809
k p
The inverse discrimination ratio is
If G (e j 0 ) = 1 this implies
20 log10 G (e j 0.25 ) 0.5
20 log10 G (e j 0.55 ) 15
37
38
N=
3rd-order lowpass Butterworth transfer

function for c = 1 is
1
Han ( s ) =
( s + 1)( s 2 + s + 1)
We choose N = 3
To determine c we use
2
1
1
=
2N
1 + ( p / c )
1 + 2
39
Denormalizing to get c = 0.588148 we

arrive at
s
Ha ( s ) = Han
0.588148
40

Design of IIR Highpass,

Highpass, Bandpass,
Bandpass,
and Bandstop Digital Filters
Applying bilinear transformation to H a (s )

we get the desired digital transfer function
G ( z ) = H a ( s ) s =1 z
First Approach (1) Prewarp digital frequency specifications

of desired digital filter GD (z ) to arrive at
frequency specifications of analog filter H D (s )
of same type
(2) Convert frequency specifications of H D (s )
into that of prototype analog lowpass filter
H LP (s )
1+ z 1
0
-10
Gain, dB
Magnitude
Magnitude and gain responses of G(z) shown

below:
0.8
0.6
0.4
41
-20
-30
0.2
0
0.2
0.4
0.6
/
0.8
-40
0
0.2
0.4
0.6
We then get
c = 1.419915 ( p ) = 0.588148
log10 (1 / k1 )
= 2.6586997
log10 (1 / k )
H a ( j p ) =
A2 1
= 15.841979


Thus
1
=
k1
0.8
Copyright
2005, S. K. Mitra
/
42
(3) Design analog lowpass filter H LP (s )


Highpass, Bandpass,
Bandpass,

Highpass, Bandpass,
Bandpass,
(4) Convert H LP (s ) into H D (s ) using

inverse frequency transformation used in
Step 2
(5) Design desired digital filter GD (z ) by
applying bilinear transformation to H D (s )
43
Second Approach (1) Prewarp digital frequency specifications

of desired digital filter GD (z ) to arrive at
frequency specifications of analog filter H D (s )
of same type
(2) Convert frequency specifications of H D (s )
into that of prototype analog lowpass filter
H LP (s )
44
IIR Highpass Digital Filter Design

Highpass, Bandpass,
Bandpass,
Design of a Type 1 Chebyshev IIR digital

highpass filter
Specifications: Fp = 700 Hz, Fs = 500 Hz,
p = 1 dB, s = 32 dB, FT = 2 kHz
Normalized angular bandedge frequencies
(3) Design analog lowpass filter H LP (s )

(4) Convert H LP (s ) into an IIR digital
transfer function GLP (z ) using bilinear
transformation
(5) Transform GLP (z ) into the desired
digital transfer function GD (z )
We illustrate the first approach
45
46
Prewarping these frequencies we get

p = tan( p / 2) = 1.9626105
s = tan(s / 2) = 1.0
MATLAB code fragments used for the design

[N, Wn] = cheb1ord(1, 1.9626105, 1, 32, s)
[B, A] = cheby1(N, 1, Wn, s);
[BT, AT] = lp2hp(B, A, 1.9626105);
[num, den] = bilinear(BT, AT, 0.5);
For the prototype analog lowpass filter choose

p =1

Using = p p we get s = 1.962105
0
Gain, dB
-10
s = 1.926105 , p = 1dB, s = 32 dB
2 Fp
2 700
=
= 0.7
FT
2000
2 Fs 2 500
s =
=
= 0.5
FT
2000
p =
47
-20
-30
-40
48
-50
0
0.2
0.4
0.6
/
0.8
1
IIR Bandpass Digital Filter Design
Design of a Butterworth IIR digital bandpass

filter
Specifications: p1 = 0.45 , p 2 = 0.65 ,
s1 = 0.3, s 2 = 0.75, p = 1 dB, s = 40 dB
Prewarping we get
p1 = tan( p1 / 2) = 0.8540807
p 2 = tan( p 2 / 2) = 1.6318517
p2
p1 = 0.777771
Width of passband Bw =
o2 =
p1
p 2 = 1.393733
s1
s 2 = 1.23010325
o2
s1and
We therefore modify s1 so that
exhibit
geometric
symmetry
with
s2
o2
respect to
s1 = 0.5773031
We set
For the prototype analog lowpass filter we
choose p = 1
s1 = tan(s1 / 2) = 0.5095254
s 2 = tan(s 2 / 2) = 2.41421356
49
50
2 2
Using = p o we get
Bw
s =
[N, Wn] = buttord(1, 2.3617627, 1, 40, s)

[B, A] = butter(N, Wn, s);
[BT, AT] = lp2bp(B, A, 1.1805647, 0.777771);
1.393733 0.3332788
= 2.3617627
0.5773031 0.777771
Specifications of prototype analog

Butterworth lowpass filter:
p = 1 , s = 2.3617627 , p = 1 dB,
s = 40 dB
51
0
Gain, dB
-10
-20
-30
-40
-50
52
0.2
0.4
0.6
0.8
/
IIR Bandstop Digital Filter Design
Design of an elliptic IIR digital bandstop filter

Specifications: s1 = 0.45 , s 2 = 0.65 ,
p1 = 0.3, p 2 = 0.75 , p = 1 dB, s = 40 dB
Prewarping we get
s 2 = 1.6318517,
s1 = 0.8540806,
p 2 = 2.4142136
p1 = 0.5095254,
s2
s1 = 0.777771
Width of stopband Bw =
p1and
p2
p1so that
We therefore modify
exhibit geometric symmetry with respect to
o2
53
o2 =
s 2
s1 = 1.393733
o2
p 2 p1 = 1.230103
p1 = 0.577303
We set
For the prototype analog lowpass filter we
choose s = 1
B
Using = s 2 w 2 we get
o
0.5095254 0.777771
p =
= 0.4234126
1.393733 0.3332787
54

[N, Wn] = ellipord(0.4234126, 1, 1, 40, s);
[B, A] = ellip(N, 1, 40, Wn, s);
[BT, AT] = lp2bs(B, A, 1.1805647, 0.777771);
0
Gain, dB
-10
-20
-30
-40
-50
0
55
0.2
0.4
0.6
0.8
/
10
Least Integral-Squared Error

Design of FIR Filters

Let H d (e j ) denote the desired frequency

response
Since H d (e j ) is a periodic function of
with a period 2, it can be expressed as a
Fourier series
In general, H d (e j ) is piecewise constant

with sharp transitions between bands
In which case, {hd [n]} is of infinite length
and noncausal
Objective - Find a finite-duration {ht [n]}
of length 2M+1 whose DTFT H t (e j )
approximates the desired DTFT H d (e j ) in
some sense
H d (e j ) = hd [n]e jn
n =
where
hd [n] =
1
1
j jn
H d (e )e d, n
2


Using Parsevals relation we can write
Commonly used approximation criterion Minimize the integral-squared error

1
j
j 2
=
H t (e ) H d (e ) d
2
where
n =
M
M 1
n =
n = M +1
= ht [n] hd [ n] + hd2 [n] + hd2 [n]

n= M
H t (e j ) = ht [n]e jn
n= M
= ht [n] hd [n]
It follows from the above that is

minimum when ht [n] = hd [n] for M n M
Best finite-length approximation to ideal
infinite-length impulse response in the
mean-square sense is obtained by truncation
Impulse Responses of Ideal

Filters

Ideal lowpass filter HLP(e j )
A causal FIR filter with an impulse response

h[n] can be derived from ht [n] by delaying:
h[n] = ht [n M ]
The causal FIR filter h[n] has the same
magnitude response as ht [n] and its phase
response has a linear phase shift of M
radians with respect to that of ht [n]
hLP [n] =
c 0
HHP (e j )
1
6
Ideal highpass filter -
sin c n
n , n
1 c ,
n=0
hHP [n] =
sin( n)
nc , n 0

Filters

Filters
Ideal bandstop filter -
Ideal bandpass filter -
HBS (e j )
HBP (e j )
1
1
c2 c1
c2 c1
c1
c2
c2
1 (c 2 c1 ) ,
n=0
hBS [n] =
sin( n) sin(c 2n)
nc1
n , n 0
sin(c 2 n) sin(c1n) , n 0
n
n
hBP [n] =
c 2 c1
n=0
,
c1
8

Filters

Filters
Ideal multiband filter HML (e
A5
A1
A4
A2
A3
Ideal discrete-time Hilbert transformer -
H ML (e j ) = Ak ,
j, < < 0
H HT (e j ) =
j , 0 < <
k 1 k ,
k = 1, 2,K, L
for n even
0,
hHT [n] =
2/ n, for n odd
sin( n)
hML [n] = ( Al Al +1) nL
L
l =1
10

Filters
Gibbs Phenomenon
Gibbs phenomenon - Oscillatory behavior in
the magnitude responses of causal FIR filters
obtained by truncating the impulse response
coefficients of ideal filters
Ideal discrete-time differentiator H DIF (e j ) = j,
1.5
11
Magnitude
n=0
0,
hDIF [ n] = cos n
n , n 0
12
N = 20
N = 60
0.5
0
0
0.2
0.4
0.6
0.8
/
Gibbs Phenomenon
13
As can be seen, as the length of the lowpass

filter is increased, the number of ripples in
both passband and stopband increases, with
a corresponding decrease in the ripple
widths
Height of the largest ripples remain the
same independent of length
Similar oscillatory behavior observed in the
magnitude responses of the truncated
versions of other types of ideal filters
Gibbs Phenomenon
Gibbs phenomenon can be explained by
treating the truncation operation as an
windowing operation:
ht [n] = hd [n] w[ n]
In the frequency domain
H t ( e j ) =
14
1
2
j
j ( )
) d
H d (e ) ( e
where H t (e j ) and (e j ) are the DTFTs

of ht [n] and w[n] , respectively
Gibbs Phenomenon
Gibbs Phenomenon
Thus H t (e j ) is obtained by a periodic

continuous convolution of H d (e j ) with
(e j )
15
If (e j ) is a very narrow pulse centered at

= 0 (ideally a delta function) compared to
variations in H d (e j ), then H t (e j ) will
approximate H d (e j ) very closely
Length 2M+1 of w[n] should be very large
On the other hand, length 2M+1 of ht [n]
should be as small as possible to reduce
computational complexity
16
Gibbs Phenomenon
Gibbs Phenomenon
A rectangular window is used to achieve

simple truncation:
1, 0 n M
wR [n] =
0, otherwise
Presence of oscillatory behavior in H t (e j )
is basically due to:
1) hd [n] is infinitely long and not absolutely
Oscillatory behavior can be explained by

j
examining the DTFT R (e ) of wR [n] :
Rectangular window
30
Amplitude
20
18
side lobe
M=4
0
-10
-1
summable, and hence filter is unstable

2) Rectangular window has an abrupt transition
to zero
17
main lobe
M = 10
10
-0.5
0
/
0.5
R (e ) has a main lobe centered at = 0

Other ripples are called sidelobes
Gibbs Phenomenon
19
Main lobe of R (e j ) characterized by its

width 4 /( 2 M + 1) defined by first zero
crossings on both sides of = 0
As M increases, width of main lobe
decreases as desired
Area under each lobe remains constant
while width of each lobe decreases with an
increase in M
Ripples in H t (e j ) around the point of
discontinuity occur more closely but with
no decrease in amplitude as M increases
Gibbs Phenomenon
Rectangular window has an abrupt transition
to zero outside the range M n M , which
results in Gibbs phenomenon in H t (e j )
Gibbs phenomenon can be reduced either:
(1) Using a window that tapers smoothly to
zero at each end, or
(2) Providing a smooth transition from
passband to stopband in the magnitude
specifications
20
Fixed Window Functions
Plots of magnitudes of the DTFTs of these

windows for M = 25 are shown below:
-20
-40
-40
-60
-100
0.2
0.4
0.6
0.8
-100
0.2
0.4
0.6
Hamming window
Blackman window
0
-20
-40
-40
Gain, dB
0
-20
-60
-100
-60
-80
-80
0.8
-60
-80
0.2
0.4
0.6
0.8
-100
0.2
0.4
0.6
0.8
1
/
Magnitude spectrum of each window

characterized by a main lobe centered at
= 0 followed by a series of sidelobes with
decreasing amplitudes
Parameters predicting the performance of a
window in filter design are:
Main lobe width
Relative sidelobe level
0
-20
-80
Hann window
Gain, dB
Gain, dB
Rectangular window
Gain, dB
Using a tapered window causes the height

of the sidelobes to diminish, with a
corresponding increase in the main lobe
width resulting in a wider transition at the
discontinuity
Hann:
w[n] = 0.5 + 0.5 cos( 2 n ), M n M
2M + 1
Hamming:
w[n] = 0.54 + 0.46 cos( 2 n ), M n M
2M + 1
Blackman:
w[n] = 0.42 + 0.5 cos( 2 n ) + 0.08 cos( 4 n )
2M + 1
2M + 1
Main lobe width ML - given by the

distance between zero crossings on both
sides of main lobe
Relative sidelobe level Asl - given by the
difference in dB between amplitudes of
largest sidelobe and main lobe

Distance between the locations of the
maximum passband deviation and minimum
stopband value ML
Observe H t (e j (c + ) ) + H t (e j (c ) ) 1
Thus,
H t (e jc ) 0.5
Passband and stopband ripples are the same
5
Width of transition band

= s p < ML
6
To ensure a fast transition from passband to

stopband, window should have a very small
main lobe width
To reduce the passband and stopband ripple
, the area under the sidelobes should be
very small
Unfortunately, these two requirements are
contradictory
In the case of rectangular, Hann, Hamming,

and Blackman windows, the value of ripple
does not depend on filter length or cutoff
frequency c , and is essentially constant
In addition,
c
M
where c is a constant for most practical
purposes
Rectangular window - ML = 4 /(2M + 1)

Asl = 13.3 dB, s = 20.9 dB, = 0.92 / M
Hann window - ML = 8 /(2M + 1)
Asl = 31.5 dB, s = 43.9 dB, = 3.11 / M
Hamming window - ML = 8 /(2M + 1)
Asl = 42.7 dB, s = 54.5 dB, = 3.32 / M
Blackman window - ML = 12 /( 2M + 1)
Asl = 58.1 dB, s = 75.3 dB, = 5.56 / M
9
Filter Design Steps (1) Set

c = ( p + s ) / 2
(2) Choose window based on specified s
(3) Estimate M using
c
M
10
FIR Filter Design Example

Lowpass filter of length 51 and c = / 2

Lowpass Filter Designed Using Hamming window
0
Gain, dB
Gain, dB
Lowpass Filter Designed Using Hann window

0
-50
An increase in the main lobe width is

associated with an increase in the width of
the transition band
A decrease in the sidelobe amplitude results
in an increase in the stopband attenuation
-50
-100
-100
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
Gain, dB
Lowpass Filter Designed Using Blackman window

0
-50
-100
11
0.2
0.4
0.6
/
0.8
12
Adjustable Window Functions

Dolph-Chebyshev window can be designed
with any specified relative sidelobe level
while the main lobe width adjusted by
choosing length appropriately
Filter order is estimated using
2.056 s 16.4
N=
2.85( )
Dolph-Chebyshev Window M
w[n] = 1 [1 + 2 Tk ( cos k ) cos 2nk ],
2M + 1
2M + 1
2M + 1
k =1
M nM
amplitude of sidelobe
where
=
main lobe amplitude
= cosh( 1 cosh 1 1 )
2M
and
cos(l cos 1 x ),
x 1
Tl ( x) =
1
cosh(l cosh x), x > 1
13
14
Properties of Dolph-Chebyshev window:

All sidelobes are of equal height
Stopband approximation error of filters
designed have essentially equiripple
behavior
For a given window length, it has the
smallest main lobe width compared to other
windows resulting in filters with the
smallest transition band
Gain response of a Dolph-Chebyshev

window of length 51 and relative sidelobe
level of 50 dB is shown below
Dolph-Chebyshev Window
Gain, dB
0
-20
-40
-60
0.2
0.4
0.6
0.8
-80
0
where is the normalized transition

bandwidth, e.g, for a lowpass filter
= s p
15
16
controls the minimum stopband

attenuation of the windowed filter response
is estimated using
Kaiser Window -
I 0{ 1 (n / M ) }
, M nM
I 0 ()
where is an adjustable parameter and I 0 (u )
is the modified zeroth-order Bessel function
of the first kind:
(u / 2) r
I 0 (u ) = 1 + [
]2
r!
r =1
Note I 0 (u ) > 0 for u > 0
20 (u / 2) r
In practice I 0 (u ) 1 + [
]2
r!
r =1
17
w[n] =
0.1102( s 8.7 ),
= 0.5842( s 21)0.4 + 0.07886( s 21),

0,
Filter order is estimated using
N=
for s > 50
for 21 s 50
for s < 21
s 8
2.285()
where is the normalized transition

bandwidth
18
sin(0.4 n)
Hence ht [ n] =
w[n], 12 n 12
n
where w[n] is the n-th coefficient of a
length-25 Kaiser window with = 3.3953
-20
Gain, dB
-20
-40
-80
0
20
Impulse Responses of FIR Filters

with a Smooth Transition
0.2
0.4
0.6
0.8
-80
0.2
0.4
0.6
0.8
Impulse Responses of FIR Filters

with a Smooth Transition
First-order spline passband-to-stopband

transition
Pth-order spline passband-to-stopband

transition
c = ( p + s ) / 2
c / ,
hLP [n] = 2 sin( n / 2 P) P sin(c n)

n / 2 P n
= s p
c / ,
hLP [n] = 2 sin( n / 2) sin(c n)

n
-40
-60
21
-60
Choose N = 24 implying M =12

19
Lowpass filter designed with Kaiser window
Kaiser Window
Gain, dB
Specifications: p = 0.3 , s = 0.5 ,

s = 40 dB
Thus c = ( p + s ) / 2 = 0.4
s = 10 s / 20 = 0.01
= 0.5842(19)0.4 + 0.07886 19 = 3.3953
32
N=
= 22.2886
2.285(0.2)
n=0
n=0
n >0
n >0
22
Lowpass FIR Filter Design

Example
Example
Magnitude
P = 1, N = 40
P = 2, N = 60
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
23
[dsp TIPS&TRICKS]
Greg Berchin
Precise Filter Design
ou have just been assigned

to a new project at work,
in which the objective is
to replace an existing analog system with a functionally equivalent digital system. Your
job is to design a digital filter that
matches the magnitude and phase
response of the existing systems analog
filter over a broad frequency range. You
are running out of ideas. The bilinear
transform and impulse invariance
methods provide poor matches to the
analog filter response, particularly at
high frequencies. Fast convolution
requires more computational resources
than you have and creates more
input/output latency than you can tolerate. What will you do?
In this article, we describe an obscure
but simple and powerful method for
designing a digital filter that approximates an arbitrary magnitude and phase
response. If applied to the problem
above, it can create a filter roughly comparable in computational burden and
latency to that created by the bilinear
transform method, with fidelity
approaching that of fast convolution. In
addition, the method described here can
also be applied to a wide variety of other
system identification tasks.
The filter design method we present
is called frequency-domain least-squares
(FDLS) [1][3]. The FDLS algorithm
produces a transfer function that approximates an arbitrary frequency response.
The input to the algorithm is a set of
magnitude and phase values at a large
number (typically thousands) of arbitrary
frequencies between 0 Hz and half the
sampling rate. The algorithms output is
a set of transfer function coefficients.
The FDLS algorithm is quite flexible in
that it can create transfer functions con-
taining poles and zeros (infinite response

filters), only zeros (finite response filters), or only poles (autoregressive networks). The algorithm uses nothing
more esoteric than basic linear algebra.
Before we can see how the technique
works, we need to review some basic linear algebra and matrix concepts.
BACKGROUND
First let us recall that, in order to
uniquely solve a system of equations, we
need as many equations as unknowns.
For example, the single equation with
one unknown 5x = 7 has the unique
solution x = 7/5. But the single equation with two unknowns 5x + 2y = 7
has multiple solutions x = (7 2y )/5
that depend on the unspecified y-value. If
another equation, 6x + 4y = 9, is
added to the equation above, there are
unique solutions for both x and y that
can be found algebraically or by matrix
inversion (denoted in the following by a
1 superscript):
5 2
6 4

x
7
=
y
9

x
5
=
y
6
1
= 38

16
7
8
=
1
7
9
1
16
7
5
9
32

9
16
.
21
45
+
16
32
2
4
Let us consider what happens if we

add another equation, x + y = 5, to the
pair that we already have above (we will
see later why we might want to do this).
There are no values of x and y that satisfy all three equations simultaneously. To
address this case, matrix algebra provides
the pseudoinverse, which determines
the values of x and y that come, in the
IEEE SIGNAL PROCESSING MAGAZINE [137] JANUARY 2007
least-squares sense, as close as possible

to satisfying all three equations. The
solution is then given by

5 2
7
x
6 4
= 9 or
y
1 1
5

1

5 2 T
5 2
x
6 4 6 4
y
1 1
1 1
T
5 2
7
for THIS case:
6 4 9
p_inv(X) =
1 1
5

[inv(X'X)][X']
0.3716
,
2.8491
where T denotes the matrix transpose. Of
course, the mathematical derivation of
the matrix inverse and pseudoinverse,
and the definition of least-squares, can
be found in any basic linear algebra text
[1]. And while our more mathematically
inclined readers will point out that there
are better ways than this to compute the
pseudoinverse, this method is adequate
for our example.
You may also remember that filter
specifications are commonly expressed in
terms of passband width and flatness,
transition band width, and stopband
attenuation. There may also be some general specifications about phase response
DSP Tips and Tricks introduces practical design and implementation signal processing algorithms that you
may wish to incorporate into your
designs. We encourage readers to
submit their contributions to
Associate Editors Rick Lyons
(r.lyons@ieee.org) or Britt Rorabaugh
(dspboss@aol.com).
1053-5888/07/$25.002007IEEE
[dsp TIPS&TRICKS]
continued
or time-domain performance, but the

exact magnitude and phase responses are
usually left to the designers discretion.
However, an important exception occurs
when a digital filter is to be used to emulate an analog filter. This is traditionally a
very difficult problem, because analog
systems are described by Laplace transforms using integration and differentiation, whereas digital systems are
described by z-transforms using delay.
Since the conversion between them is
nonlinear, the response of an analog system can only be approximated by a digital
system and vice-versa.
Let us assume that the transfer function of our digital filter (i.e., the mathematical description of the relationship
between the filters input and output) is
in a standard textbook form [2] given by
Y(z )
b0 + b1 z 1 + + bN z N
=
,
U(z )
1 + a1 z 1 + + aD z D
where U(z ) is the z-transform of the input
signal, Y(z ) is the z-transform of the output signal, and the a and b factors are realvalued coefficients. Furthermore, we
assume that the filter is causal, meaning
that its response to an input does not begin
until after the input is applied. Under these
assumptions, the time-domain difference
equation that implements our filter is
y(k) = a1 y(k 1) aD y(k D)
ues of which are not yet known. We also

know that the relationship between
input u and output y at any sample time
can be inferred from the frequency
response value A at frequency .
Combining these two ideas, we obtain
one equation in D + N + 1 unknowns
Let us conclude our background section with a comment on what a frequency response value means. In a simple
example, if the frequency response of a
system at a frequency 1 is given in
magnitude/phase form as A1 1 , the
output amplitude will be A1 times the
input amplitude and the output phase
will be shifted an angle 1 relative to the
input phase when a steady-state sine
wave of frequency 1 is applied to the
system. For instance, if the input to the
system described above at time k is
u1 (k) = cos(k1 ts ), where ts is the sampling period (equal to one over the sampling frequency), then the output will be
y1 (k) = A1 cos(k1 ts + 1 ) . The input
and output values at any sample time
can be determined in a similar manner.
For example, the input sample value N
samples
in
the
past
was
u1 (k N) = cos ((k N )1 ts ) and the
output sample value D samples in the
past was y1 (k D) = A1 cos ((k D)
1 ts + 1 ). For our purposes, since k
represents the current sample time, its
value can conveniently be set to zero.
y1 (0) = [y1 (1) . . . y1 (D)

a
1
..
.

aD
u1 (0) . . . u1 (N )]
b0 .

.
..
bN
(Note that the current-sample index k
has been set to zero.) If we repeat using
A2 2 at a different frequency 2 , we
obtain a second equation in D + N + 1
unknowns as shown in (a) at the bottom
of the page. And if we repeat at many
more different frequencies M than we
have unknowns D + N + 1, we know
from our review of linear algebra that the
pseudoinverse will compute values for
the set of coefficients a1 . . . aD and
b0 . . . bN that come as close as possible
to solving all of the equations, which is
exactly what we need to design our filter.
So now we can write (b), shown at the
bottom of the page. We can denote the
y1 (0) . . . yM (0) column vector above as
Y, the matrix as X, and the a1 . . . bN col-
FDLS FILTER APPROXIMATION

Based on our review of the pseudoinverse, transfer function, and frequency
response, we know that the output is a
combination of present and past input
and output values, each scaled by a set of
b or a coefficients (respectively), the val-
+ b0 u(k) + + bN u(k N ),
where the a and b coefficients are exactly
the same as in the transfer function
above, k is the time index, u(k) and y(k)
are the current values of the input and
output (respectively), u(k N ) was the
input value N samples in the past, and
y(k D) was the output value D samples
in the past. We can write the equation
above in matrix form as
y(k) = [y(k 1) . . . y(k D)
a
1
..
.

aD
u(k) . . . u(k N )]
b0 .

.
..

y1 (1) . . .
y1 (0)
=
y2 (0)
y2 (1) . . .
y1 (D)
y2 (D)
u1 (0) . . .
u2 (0) . . .
a
1
..
.

u1 (N)
aD .
u2 (N)
b
0

.
..
(a)
bN
y1 (1) . . .
y1 (0)
y2 (0) y2 (1) . . .

.. =
..
.
.
yM (0)
yM (1) . . .
y1 (D)
y2 (D)
..
.
u1 (0)
u2 (0)
..
.
yM (D)
uM (0) . . .
bN
...
...
a
1
..
u1 (N) .

u2 (N)
aD
..
b0 .

.
.
uM (N) ..
bN
(b)
We can now summarize our filter

design trick as follows:
1) Select the numerator order N and
the denominator order D, where N
and D do not have to be equal and
either one (but not both) may be zero.
(We have found no rule of thumb
for defining N and D; they are best
determined experimentally.)
2) Define the M separate input um
cosine sequences, each of length
(N + 1).
3) Compute the M separate output ym
cosine sequences, each of length D
(based on Am m ).
4) Fill the X matrix with the input um
and output ym cosine sequences.
5) Fill the Y vector with the M output
cosine values, ym (0) = Am cos( m ).
6) Compute the pseudoinverse; the
resulting vector contains the filter
coefficients.
A numerical example is shown in
Figures 1 and 2, which illustrate the
magnitude and phase, respectively, of a
real-world example analog system (black)
and of the associated bilinear transform
(blue), impulse invariance (green), and
FDLS (red) approximations. The sampling rate is equal to 240 Hz and
D = N = 12. The red FDLS graphs are
almost completely obscured by the black
analog system graphs. In this example,
the FDLS errors are often three to four
orders of magnitude smaller than those
of the other methods. (In Figure 2, the
bilinear transform curve is obscured by
the FDLS and analog curves at low frequencies and by the impulse invariance
curve at high frequencies.)
In terms of the computational complexity of an FDIS-designed filter, the
number of feedback and feed-forward
coefficients is determined by the variables D and N, respectively. As such, an
FDIS-designed
filter
requires
(N + D + 1) multiplies and (N + D)
additions per filter output sample.
Magnitude (dB)
(X T X )1 X T Y .
0
10
20
30
40
Black: Analog (Desired)

Blue: Bilinear Transform
Green: Impulse Invariance
Red: FDLS
50
102
101
100
Frequency (Hz)
101
102
[FIG1] Magnitude responses of the filter designed using the bilinear transform, impulse
invariance, and FDLS methods.
150
Phase (degrees)
umn vector as . With these notations,

Y = X, and the pseudoinverse solves
for the vector that contains the
desired filter coefficients
100
50
0
50
100
150
102
Black: Analog (Desired)

Blue: Bilinear Transform
Green: Impulse Invariance
Red: FDLS
101
100
Frequency (Hz)
101
102
[FIG2] Phase responses of the filter designed using the bilinear transform, impulse
invariance, and FDLS methods.
CONCLUSION
FDLS is a powerful method for designing digital filters. As is the case with all
approximation techniques, there are
circumstances in which the FDLS
method works well and others in which
it does not. The FDLS method does not
replace other filter design methods; it
provides one more method from which
to choose. FDLS is most useful in cases
where a specified frequency response
must be duplicated to within tight tolerances over a wide frequency range or
when the frequency response of an
existing system is known but the coefficients of the systems transfer function are unknown. It is up to the
designer to determine whether to use
it in any given situation. Detailed
examples and a MATLAB code implementation of the FDLS algorithm are
available at http://apollo.ee.columbia.
edu/spm/?i=external/tipsandtricks.
ACKNOWLEDGMENTS
My thanks to Jaime Andrs Aranguren
Cardona for providing the example
shown in Figures 1 and 2, originally
posted on the comp.dsp Newsgroup on
9 January 2005, which ultimately led to
this article.
AUTHOR
Greg Berchin (berchin@ieee.org) is a signal processing algorithm engineer who
provides contract engineering services
from Naperville, Illinois.
REFERENCES
[1] G. Strang, Linear Algebra and Its Applications,

2nd ed., Orlando, FL: Academic, pp. 103152, 1980.
[2] R. Lyons, Understanding Digital Signal
Processing, 2nd ed., Upper Saddle River, NJ:
Prentice-Hall, pp. 232240, 2004.
[3] G. Berchin, A new algorithm for system identification from frequency response information, masters thesis, University of California-Davis, 1988 .
[4] G. Berchin and M.A. Soderstrand, A transformdomain least-squares beamforming technique, in
Proc. IEEE Oceans 90 Conf., Arlington, VA, Sept.
1990.
[SP]
FDLS Examples
by Greg Berchin and Richard Lyons [January 2007]
The following material provides additional examples of the FDLS algorithm described in the
January 2007 IEEE Signal Processing magazine DSP Tips & Tricks column article "Precise Filter
Design" by Greg Berchin.
Algebraic Example
Recall the FDLS matrix expression
a1

y1 (0) y1 ( 1) K y1 ( D ) u1 (0) K u1 ( N ) M
y (0) y ( 1) K y ( D ) u (0) K u ( N ) a
2
2
2
2
D ,
2 =
b0
M
M
M
M
M

y M (0) y M ( 1) K y M ( D ) u M (0) K u M ( N ) M

bN
which we wrote as Y = X.
Each individual element in the Y column vector is of the form Amcos(m), and each
element in the X matrix is of the form A1cos(k1ts + 1) or cos(k1ts). Because all of these
elements are of the form of a product (Amplitude)[cos(angle)], each element in Y and X is equal
to a constant.
Now if, say, D = 10 and N = 9, then:
y1(0) = A1cos(1)
is the first element of the Y column vector and

[y1(1) ... y1(10) u1(0) ... u1(9)]
is the top row of the X matrix expression where
y1(1) = A1cos[(1)1ts + 1] = A1cos(1ts + 1)
y1(2) = A1cos[(2)1ts + 1] = A1cos(21ts + 1)
...
y1(10) = A1cos[(10)1ts + 1] = A1cos(101ts + 1)
and
u1(0) = cos[(0)1ts] = 1
u1(1) = cos[(1)1ts] = cos(1ts)
u1(2) = cos[(2)1ts] = cos(21ts)
...
u1(9) = cos[(9)1ts] = cos(91ts).
IEEE Copyright All Rights Reserved
Page 1 of 5
So the top row of the X matrix looks like:

[A1cos(1ts + 1) A1cos(21ts + 1) ... A1cos(101ts + 1) 1 cos(1ts)
cos(21ts) ... cos(91ts)].
The second row of the X matrix looks like:
[A2cos(2ts + 2) A2cos(22ts + 2) ... A2cos(102ts + 2) 1 cos(2ts)
cos(22ts) ... cos(92ts)].
And so on.
Numerical Example
Heres an example of the above expressions using actual numbers. Suppose we

need to approximate the transfer function coefficients for the system whose frequency
magnitude and phase response is that shown in Figure E1. Assume that our discretesystem sample rate is 1000 Hz, thus ts = 103 seconds, and N = D = 2. Also assume
M = 8 and we have the eight A1 -to- A8 magnitude sample values and the eight 1 -to- 8
phase samples, shown as dots in Figure E1, available to us as input values to the FDLS
algorithm.
M=8
m = 2 . fm :
1.5
A5
Linear
1
Desired magnitude response
0.5
A1
(a)
A8
0
0
100
f1 = 0 Hz
200
300
Frequency
f7 = 106.3 Hz
400
500
1 = 2(0.0)
2 = 2(19.69)
3 = 2(35.43)
4 = 2(51.18)
5 = 2(59.05)
6 = 2(66.93)
7 = 2(106.30)
8 = 2(389.76)
f8 = 389.76 Hz
f4 = 51.18 Hz
(b)
Phase (radians)
4
5
2
1
Desired phase response
0
2
8
0
100
200
300
Frequency
400
500
Fig E1
Page 2 of 5
In matrix form, the target analog system parameters are
0.0
19.6850
35.4331
fm = 51.1811
59.0551
66.9291
106.299
389.764
0.0
0.1237
0.2226
mts = 0.3216
0.3711
0.4205
0.6679
2.449
0.2172
0.2065
0.1696
Am = 0.0164
1.3959
0.6734
0.3490
0.3095
0.0
0.0156
0.0383
m = 3.0125
2.3087
0.955
0.0343
0.0031
where the fm vector is in Hz, the mts vector is in radians, and 1 m 8. The first two
elements of the Y vector are:
y1(0) = A1cos(1) = 0.2172cos(0) = 0.2172.
y2(0) = A2cos(2) = 0.2065cos(0.0156) = 0.2065.
The complete Y vector is:
Y =
A1cos(1)
A2cos(2)
A3cos(3)
A4cos(4)
A5cos(5)
A6cos(6)
A7cos(7)
A8cos(8)
0.2172
0.2065
0.1695
= 0.0162
0.9390
0.6605
0.3488
0.3095
The two elements of the "y1" part of the first row of the X vector are:
y1(1) = A1cos(1ts + 1) = 0.2172cos(0 + 0) = 0.2172.
y1(2) = A1cos(21ts + 1) = 0.2172cos(0 + 0) = 0.2172.
The two elements of the "y8" part of the eighth row of the X vector are:
y8(1) = A8cos(8ts + 8) = 0.3095cos(2.449 + 0.0031) = 0.2376.
y8(2) = A8cos(28ts + 8) = 0.3095cos(4.898 + 0.0031) = 0.562.
Page 3 of 5
The three elements of the "u1" part of the first row of the X matrix are:
u1(0) = cos(0) = 1
u1(1) = cos(1ts) = cos(0) = 1
u1(2) = cos(21ts) = cos(0) = 1.
The three elements of the "u8" part of the eighth row of the X matrix are:
u8(0) = cos(0) = 1
u8(1) = cos(8ts) = cos(2.449) = 0.7696
u8(2) = cos(28ts) = cos(4.898) = 0.1845.
The complete X matrix is:
A1cos(11ts+1)
A2cos(12ts+2)
A3cos(13ts+3)
X = A4cos(14ts+4)
A5cos(15ts+5)
A6cos(16ts+6)
A7cos(17ts+7)
A8cos(18ts+8)
0.2172
0.2045
0.1639
= 0.0147
0.5007
0.6564
0.2812
0.2376
0.2172
0.994
0.1502
0.0117
0.0059
0.5378
0.0928
0.0562
A1cos(21ts+1)
A2cos(22ts+2)
A3cos(23ts+3)
A4cos(24ts+4)
A5cos(25ts+5)
A6cos(26ts+6)
A7cos(27ts+7)
A8cos(28ts+8)
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.0000
0.9924
0.9753
0.9487
0.939
0.9129
0.7851
0.7696
cos(0)
cos(0)
cos(0)
cos(0)
cos(0)
cos(0)
cos(0)
cos(0)
cos(1ts)
cos(2ts)
cos(3ts)
cos(4ts)
cos(5ts)
cos(6ts)
cos(7ts)
cos(8ts)
cos(21ts)
cos(22ts)
cos(23ts)
cos(24ts)
cos(25ts)
cos(26ts)
cos(27ts)
cos(28ts)
1.0000
0.9696
0.9025
0.8002 .
0.7370
0.6667
0.2328
0.1845
Given the above Y vector and the X matrix, the FDLS algorithm computes the 2nd-order
(N = D = 2) transfer function coefficients vector M=8 as
M=8 =
1.8439
0.9842
0.3033
0.5762
0.3034
Page 4 of 5
Treated as filter coefficients, we can write vector M=8 as:

a0 = 1
a1 = 1.8439
a2 = 0.9842
b0 = 0.3033
b1 = 0.5762
b2 = 0.3034
implemented as the recursive filter network shown in Figure E2.

u(k)
y(k)
z 1
1.8439
z 1
0.3033
0.5762
0.3034
0.9842
Fig E2
(a)
Magnitude (dB)
The frequency-domain performance of the filter are the solid red curves shown in Figure E3.
There we see that the M=8 coefficients provide an accurate approximation to the desired
frequency response in Figure E1.
20
Desired magnitude response

Actual magnitude response
40
0
100
200
300
Frequency
400
500
400
500
(b)
Phase (radians)
4
2
Desired phase response

Actual phase response
0
2
0
100
200
300
Frequency
Fig E3
Page 5 of 5

IIR Filters

Uploaded by

Document Informationclick to expand document informationIIR_Filters

Document Informationclick to expand document information

Copyright:

Available Formats

IIR Filters

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

IIR Filters

Uploaded by

Copyright:

Available Formats

Analog Lowpass Filter

Analog Lowpass Filter

Typical magnitude response H a ( j) of an

In the passband, defined by 0 p , we

i.e., H a ( j) approximates zero within an

Copyright 2005, S. K. Mitra

Copyright 2005, S. K. Mitra

Analog Lowpass Filter

Analog Lowpass Filter

p - passband edge frequency

Magnitude specifications may alternately be

Minimum stopband attenuation

Copyright 2005, S. K. Mitra

Analog Lowpass Filter

Analog Lowpass Filter Design

Here, the maximum value of the magnitude

(1) Transition ratio k =

1 / 1 + 2 - Maximum passband deviation,

For a lowpass filter k < 1

1 - Maximum stopband magnitude

Copyright 2005, S. K. Mitra

The magnitude-square response of an N-th

Copyright 2005, S. K. Mitra

Copyright 2005, S. K. Mitra

Typical magnitude responses with c = 1

Two parameters completely characterizing a

Copyright 2005, S. K. Mitra

c and N are thus determined from

Since order N must be an integer, value

Copyright 2005, S. K. Mitra

The magnitude-square response of an N-th

where TN () is the Chebyshev polynomial

Copyright 2005, S. K. Mitra

If at = s the magnitude is equal to 1/A,

Typical magnitude response plots of the

Solving the above we get

Copyright 2005, S. K. Mitra

The magnitude-square response of an N-th

Typical magnitude response plots of the

where TN () is the Chebyshev polynomial

Copyright 2005, S. K. Mitra

The square-magnitude response of an

The order N of the Type 2 Chebyshev filter

Copyright 2005, S. K. Mitra

For given p , s , , and A, the filter order

Example - Determine the lowest order of a elliptic

Note: k = 0.2 and 1 / k1 = 196.5134

Copyright 2005, S. K. Mitra

Copyright 2005, S. K. Mitra

Analog Lowpass Filter Design

Example - Design an elliptic lowpass filter

Typical magnitude response plots with p = 1

[N, Wn] = ellipord(Wp, Ws, Rp, Rs, s);

Copyright 2005, S. K. Mitra

Design of Analog Highpass,

Analog Lowpass Filter Design

[b, a] = ellip(N, Rp, Rs, Wn, s);

Copyright 2005, S. K. Mitra