IIR digital filter design for a low pass filter based on impulse invariance and bilinear transformation methods using butterworth analog filter

Nasser M. Abbasi
May 10, 2010

1 Introduction

This report derives a symbolic procedure to design a low pass IIR digital filter from an analog Butterworth filter using 2 methods: impulse invariance and bilinear transformation. Two numerical examples are then used to illustrate using the symbolic procedure. There are a total of 13 steps.

A Mathematica program with GUI is written to enable one to use this design for different parameters. A Matlab script is being written as well.

1.1 Filter specifications

Filter specifications are 5 parameters. The frequency specifications are analog frequencies, while the attenuations for the passband and the stopband are given in db

$\displaystyle % \begin{tabular}[c]{\vert l\vert l\vert}\hline $F_{s}$\ & The ... ...$\delta_{s}$\ & The attenuation in db at $f_{s}$\\ \hline \end{tabular} \ \$

This diagram below illustrates these specifications

The frequency specifications are in

and they must be first converted to digital frequencies $\omega$ where $0\leq\left\vert \omega\right\vert \leq\pi$ before using the attenuation specifications, The sampling frequency $F_{s}$ is used to do this conversion since $F_{s}$ corresponds to $2\pi$ on the digital frequency scale.

2 Analytical derivation of the design steps

These are the design steps

Convert specifications from analog to digital frequencies
Based on design method (impulse invariance of bilinear) apply the attenuation criteria to determine $\Omega_{c}$ and (the filter order)
Using $\Omega_{c}$ and find the locations of the poles of $H\left( s\right)$ , the butterworth analog filter.
Find $H\left( z\right)$ from $H\left( s\right)$ . The method of doing this depends if we are using impulse invariance or bilinear. This step is much simpler for the bilinear method as it does not require performing partial fractions decomposition on

Now we begin the analytical design procedure.

2.1 Impulse invariance method

We first convert analog specifications to digital specifications: $\frac {F_{s}}{2\pi}=\frac{f_{p}}{\omega_{p}}$ , hence $\omega_{p}=2\pi\frac{f_{p}% }{F_{s}}$ and $\omega_{s}=2\pi\frac{f_{s}}{F_{s}}$

Convert the criteria relative to the digital normalized scale:

$\displaystyle 20\log\left\vert H\left( e^{j\omega_{p}}\right) \right\vert$	$\displaystyle \geq \delta_{p}$
$\displaystyle 20\log\left\vert H\left( e^{j\omega_{s}}\right) \right\vert$	$\displaystyle \leq \delta_{s}%$

Therefore

$\displaystyle \left\vert H\left( e^{j\omega_{p}}\right) \right\vert$	$\displaystyle \geq 10^{\frac{\delta_{p}}{20}}T$	(A)
$\displaystyle \left\vert H\left( e^{j\omega_{s}}\right) \right\vert$	$\displaystyle \leq 10^{\frac{\delta_{s}}{20}}T%$	(B)

Notice that we added the

scaling factor. Now, Butterworth analog filter squared magnitude Fourier transform is given by

$\displaystyle \left\vert H_{a}\left( j\Omega\right) \right\vert ^{2}=\frac{T^{2}% }{1+\left( \frac{j\Omega}{j\Omega_{c}}\right) ^{2N}}%$

hence equations (A) and (B) above are now written in terms of the analog Butterworth amplitude frequency response and become

$\displaystyle \frac{T^{2}}{1+\left( \frac{\Omega_{p}}{\Omega_{c}}\right) ^{2N}}$	$\displaystyle \geq\left( 10^{\frac{\delta_{p}}{20}}\right) ^{2}T^{2}$
$\displaystyle \frac{T^{2}}{1+\left( \frac{\Omega_{s}}{\Omega_{c}}\right) ^{2N}}$	$\displaystyle \leq\left( 10^{\frac{\delta_{s}}{20}}\right) ^{2}T^{2}%$

Therefore the above becomes

$\displaystyle \frac{1}{1+\left( \frac{\Omega_{p}}{\Omega_{c}}\right) ^{2N}}$	$\displaystyle \geq\ 10^{\frac{\delta_{p}}{10}}$
$\displaystyle \frac{1}{1+\left( \frac{\Omega_{s}}{\Omega_{c}}\right) ^{2N}}$	$\displaystyle \leq10^{\frac{\delta_{s}}{10}}%$

Now, for impulse invariance, $\Omega_{p}=\frac{\omega_{p}}{T}\,\$ and $\Omega_{s}=\frac{\omega_{s}}{T}$ .

$\displaystyle 1+\left( \frac{\omega_{p}/T}{\Omega_{c}}\right) ^{2N}$	$\displaystyle \leq10^{-\frac {\delta_{p}}{10}}$	(1)
$\displaystyle 1+\left( \frac{\omega_{s}/T}{\Omega_{c}}\right) ^{2N}$	$\displaystyle \geq10^{-\frac {\delta_{s}}{10}}%$	(2)

Change inequalities to equalities and simplify

$\displaystyle \left( \frac{\omega_{p}/T}{\Omega_{c}}\right) ^{2N}$	$\displaystyle =10^{-\frac {\delta_{p}}{10}}-1$
$\displaystyle \left( \frac{\omega_{s}/T}{\Omega_{c}}\right) ^{2N}$	$\displaystyle =10^{-\frac {\delta_{s}}{10}}-1$

Divide the above 2 equations

$\displaystyle \left( \frac{\omega_{p}}{\omega_{s}}\right) ^{2N}$	$\displaystyle =\frac{10^{-\frac {\delta_{p}}{10}}-1}{10^{-\frac{\delta_{s}}{10}}-1}$
$\displaystyle 2N\left[ \log\left( \omega_{p}\right) -\log\left( \omega_{s}\right) \right]$	$\displaystyle =\log\left( 10^{-\frac{\delta_{p}}{10}}-1\right) -\log\left( 10^{-\frac{\delta_{s}}{10}}-1\right)$
$\displaystyle N$	$\displaystyle =\frac{1}{2}\frac{\log\left[ 10^{-\frac{\delta_{p}}{10}}-1\right]... ...}}-1\right] }{\log\left( \omega _{p}\right) -\log\left( \omega_{s}\right) }%$

We need to round the above to then nearest integer using the Ceiling function $\;$ i.e.

$\left\lceil N\right\rceil$

Now for impulse invariance method, use equation (1) above to solve for $\Omega_{c}$

$\displaystyle \left( \frac{\omega_{p}/T}{\Omega_{c}}\right) ^{2N}$	$\displaystyle =10^{-\frac {\delta_{p}}{10}}-1$
$\displaystyle 2N\left( \log_{10}\frac{\omega_{p}/T}{\Omega_{c}}\right)$	$\displaystyle =\log _{10}\left( 10^{-\frac{\delta_{p}}{10}}-1\right)$
$\displaystyle \log_{10}\frac{\omega_{p}/T}{\Omega_{c}}$	$\displaystyle =\frac{1}{2N}\log_{10}\left( 10^{-\frac{\delta_{p}}{10}}-1\right)$
$\displaystyle \frac{\omega_{p}/T}{\Omega_{c}}$	$\displaystyle =10^{\left( \frac{1}{2N}\log_{10}\left( 10^{-\frac{\delta_{p}}{10}}-1\right) \right) }$
$\displaystyle \Omega_{c}$	$\displaystyle =\frac{\omega_{p}/T}{10^{\left( \frac{1}{2N}\log_{10}\left( 10^{-\frac{\delta_{p}}{10}}-1\right) \right) }}%$

Now that we found and $\Omega_{c}$ , next find the poles of $H\left( s\right) .$ Since the butterworth magnitude square of the transfer function is

$\displaystyle \left\vert H_{a}\left( s\right) \right\vert ^{2}=\frac{T^{2}}{1+\left( \frac{s}{j\Omega_{c}}\right) ^{2N}}%$

Hence $H\left( s\right)$ poles are found by setting the denominator of the above to zero

$\displaystyle 1+\left( \frac{s}{j\Omega_{c}}\right) ^{2N}$	$\displaystyle =0$
$\displaystyle \left( \frac{s}{j\Omega_{c}}\right) ^{2N}$	$\displaystyle =-1$
	$\displaystyle =e^{j\left( \pi+2\pi k\right) }\ \ \ \ \ k=0,1,2,\cdots2N-1$
$\displaystyle \frac{s}{j\Omega_{c}}$	$\displaystyle =e^{j\left( \frac{\pi+2\pi k}{2N}\right) }$
$\displaystyle s$	$\displaystyle =j\Omega_{c}\ e^{j\left( \frac{\pi+2\pi k}{2N}\right) }$
	$\displaystyle =\Omega_{c}\ e^{j\frac{\pi}{2}}e^{j\left( \frac{\pi+2\pi k}{2N}\right) }$
	$\displaystyle =\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+2k+N\right) }{2N}\right) }%$

We only need to find the LHS poles, which are located at $i=0\cdots N-1$ , because these are the stable poles. Hence the $i^{th}$ pole is

$\displaystyle s_{i}=\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+2i+N\right) }{2N}\right) }%$

For example for

we get

$\displaystyle s_{0}=\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+N\right) }{2N}\right) }=\Omega_{c}\ e^{j\left( \frac{\pi7}{12}\right) }%$

Now we can write the analog filter generated based on the above selected poles, which is, for impulse invariance

$\displaystyle H_{a}\left( s\right) =\frac{T\ K}{% {\displaystyle\prod\limits_{i=0}^{N-1}} \left( s-s_{i}\right) }%$

(3)

is found by solving $H_{a}\left( 0\right) =T$ hence

$\displaystyle k=% {\displaystyle\prod\limits_{i=0}^{N-1}} \left( -s_{i}\right)$

Now we need to write poles in non-polar form and plug them into (3)

$\displaystyle s_{i}=\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+2i+N\right) }{2N}\... ...+N\right) }{2N}\right) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ i=0\cdots N-1$

Hence,

$\displaystyle H_{a}\left( s\right) =\frac{T\ K}{% {\displaystyle\prod\limits_{... ...+N\right) }{2N}% +j\sin\frac{\pi\left( 1+2i+N\right) }{2N}\right) \right) }%$

(4)

Where

$\displaystyle \Omega_{c}=\frac{\omega_{p}/T}{10^{\left( \frac{1}{2N}\log_{10}\left( 10^{-\frac{\delta_{p}}{10}}-1\right) \right) }}%$

and

$\displaystyle N=\left\lceil \frac{1}{2}\frac{\log\left[ 10^{-\frac{\delta_{p}}{... ...] }{\log\left( \omega_{p}\right) -\log\left( \omega_{s}\right) }\right\rceil$

Now that we have found $H\left( s\right)$ we need to convert it to $H\left( z\right)$

We need to make sure that we multiply poles of complex conjugates with each others to make the result simple to see.

Now that we have $H_{a}\left( s\right)$ , we do the $A\rightarrow D$ conversion. I.e. obtain $H\left( z\right)$ from the above $H\left( s\right)$ . When using impulse invariance, we need to perform partial fraction decomposition on (4) above in order to write $H\left( s\right)$ in this form

$\displaystyle H\left( s\right) =% {\displaystyle\sum\limits_{i=0}^{N-1}} \frac{A_{i}}{s-s_{i}}%$

For example, to obtain $A_{j}$ , we write

$\displaystyle A_{j}=\lim_{s\rightarrow s_{j}}H_{a}\left( s\right) =\frac{T\ k}{... ...style\prod\limits_{\substack{i=0\\ i\neq j}}^{N-1}} \left( s-s_{i}\right) }%$

Once we find all the $A^{\prime}s$ , we now write $H\left( z\right)$ as follows

$\displaystyle H\left( z\right) =% {\displaystyle\sum\limits_{i=0}^{N-1}} \frac{A_{i}}{1-\exp\left( s_{i}\ T\right) z^{-1}}%$

This completes the design. We can try to convert the above form of $H\left( z\right)$ to a rational expression as $H\left( z\right) =\frac{N\left( z\right) }{D\left( z\right) }$

2.2 bilinear transformation method

Convert the criteria relative to the digital normalized scale:

$\displaystyle 20\log\left\vert H\left( e^{j\omega_{p}}\right) \right\vert$	$\displaystyle \geq \delta_{p}$
$\displaystyle 20\log\left\vert H\left( e^{j\omega_{s}}\right) \right\vert$	$\displaystyle \leq \delta_{s}%$

Hence

$\displaystyle \left\vert H\left( e^{j\omega_{p}}\right) \right\vert$	$\displaystyle \geq 10^{\frac{\delta_{p}}{20}}$	(A)
$\displaystyle \left\vert H\left( e^{j\omega_{s}}\right) \right\vert$	$\displaystyle \leq 10^{\frac{\delta_{s}}{20}}%$	(B)

Butterworth analog filter squared magnitude Fourier transform is given by

$\displaystyle \left\vert H_{a}\left( j\Omega\right) \right\vert ^{2}=\frac{1}{1+\left( \frac{j\Omega}{j\Omega_{c}}\right) ^{2N}}%$

hence equations (A) and (B) above are now written in terms of the analog Butterworth amplitude frequency response and become

$\displaystyle \frac{1}{1+\left( \frac{\Omega_{p}}{\Omega_{c}}\right) ^{2N}}$	$\displaystyle \geq\left( 10^{\frac{\delta_{p}}{20}}\right) ^{2}=\ 10^{\frac{\delta_{p}}{10}}$
$\displaystyle \frac{1}{1+\left( \frac{\Omega_{s}}{\Omega_{c}}\right) ^{2N}}$	$\displaystyle \leq\left( 10^{\frac{\delta_{s}}{20}}\right) ^{2}=10^{\frac{\delta_{s}}{10}}%$

Now we assign values for $\Omega_{p}$ and $\Omega_{s}$ as follows, $\Omega _{p}=\frac{2}{T}\tan\left( \frac{\omega_{p}}{2}\right)$ , $\Omega_{s}% =\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right)$

$\displaystyle 1+\left( \frac{\frac{2}{T}\tan\left( \frac{\omega_{p}}{2}\right) }% {\Omega_{c}}\right) ^{2N}$	$\displaystyle \leq10^{\frac{\delta_{p}}{10}}$	(1)
$\displaystyle 1+\left( \frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right) }% {\Omega_{c}}\right) ^{2N}$	$\displaystyle \geq10^{\frac{\delta_{s}}{10}}%$	(2)

Change inequalities to equalities and simplify

$\displaystyle \left( \frac{\frac{2}{T}\tan\left( \frac{\omega_{p}}{2}\right) }{\Omega _{c}}\right) ^{2N}$	$\displaystyle =10^{\frac{\delta_{p}}{10}}-1$
$\displaystyle \left( \frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right) }{\Omega _{c}}\right) ^{2N}$	$\displaystyle =10^{\frac{\delta_{s}}{10}}-1$

Divide the above 2 equations

$\displaystyle \left( \frac{\tan\left( \frac{\omega_{p}}{2}\right) }{\tan\left( \frac{\omega_{s}}{2}\right) }\right) ^{2N}$	$\displaystyle =\frac{10^{\frac{\delta_{p}% }{10}}-1}{10^{\frac{\delta_{s}}{10}}-1}$
$\displaystyle 2N\left[ \log\left( \tan\left( \frac{\omega_{p}}{2}\right) \right) -\log\left( \tan\left( \frac{\omega_{s}}{2}\right) \right) \right]$	$\displaystyle =\log\left( 10^{\frac{\delta_{p}}{10}}-1\right) -\log\left( 10^{\frac {\delta_{s}}{10}}-1\right)$
$\displaystyle N$	$\displaystyle =\frac{1}{2}\frac{\log\left( 10^{\frac{\delta_{p}}{10}}-1\right)... ...ight) \right) -\log\left( \tan\left( \frac {\omega_{s}}{2}\right) \right) }%$

We need to round the above to then nearest integer using the Ceiling function $\;$ i.e.

$\left\lceil N\right\rceil$

Now for bilinear transformation we used equation (2) above to find $\Omega_{c},$ Hence we now solve for $\Omega_{c}$

$\displaystyle 1+\left( \frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right) }% {\Omega_{c}}\right) ^{2N}$	$\displaystyle =10^{\frac{\delta_{s}}{10}}$
$\displaystyle 2N\left( \log_{10}\frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right) }{\Omega_{c}}\right)$	$\displaystyle =\log_{10}\left( 10^{\frac{\delta_{s}}{10}% }-1\right)$
$\displaystyle \log_{10}\frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right) }% {\Omega_{c}}$	$\displaystyle =\frac{1}{2N}\log_{10}\left( 10^{\frac{\delta_{s}}{10}% }-1\right)$
$\displaystyle \frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right) }{\Omega_{c}}$	$\displaystyle =10^{\left( \frac{1}{2N}\log_{10}\left( 10^{\frac{\delta_{s}}{10}}-1\right) \right) }$
$\displaystyle \Omega_{c}$	$\displaystyle =\frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right) }{10^{\... ...frac{1}{2N}\log_{10}\left( 10^{\frac{\delta_{s}}{10}% }-1\right) \right) }}%$

Now that we found and $\Omega_{c}$ we find the poles of $H\left( s\right)$ . Since for bilinear the magnitude square of the transfer function is

$\displaystyle \left\vert H_{a}\left( s\right) \right\vert ^{2}=\frac{1}{1+\left( \frac {s}{j\Omega_{c}}\right) ^{2N}}%$

Hence $H\left( s\right)$ poles are found by setting the denominator of the above to zero

$\displaystyle 1+\left( \frac{s}{j\Omega_{c}}\right) ^{2N}$	$\displaystyle =0$
$\displaystyle \left( \frac{s}{j\Omega_{c}}\right) ^{2N}$	$\displaystyle =-1$
	$\displaystyle =e^{j\left( \pi+2\pi k\right) }\ \ \ \ \ k=0,1,2,\cdots2N-1$
$\displaystyle \frac{s}{j\Omega_{c}}$	$\displaystyle =e^{j\left( \frac{\pi+2\pi k}{2N}\right) }$
$\displaystyle s$	$\displaystyle =j\Omega_{c}\ e^{j\left( \frac{\pi+2\pi k}{2N}\right) }$
	$\displaystyle =\Omega_{c}\ e^{j\frac{\pi}{2}}e^{j\left( \frac{\pi+2\pi k}{2N}\right) }$
	$\displaystyle =\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+2k+N\right) }{2N}\right) }%$

We only need to find the LHS poles, which are located at $i=0\cdots N-1$ , because these are the stable poles. Hence the $i^{th}$ pole is

$\displaystyle s_{i}=\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+2i+N\right) }{2N}\right) }%$

For example for

we get

$\displaystyle s_{0}=\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+N\right) }{2N}\right) }=\Omega_{c}\ e^{j\left( \frac{\pi7}{12}\right) }%$

For bilinear,

is given by

$\displaystyle H_{a}\left( s\right) =\frac{K}{% {\displaystyle\prod\limits_{i=0}^{N-1}} \left( s-s_{i}\right) }%$

(3)

is found by solving $H_{a}\left( 0\right) =1$ , hence we obtain

$\displaystyle k=% {\displaystyle\prod\limits_{i=0}^{N-1}} \left( -s_{i}\right)$

We see that the same expression results for

for both cases.

Now we need to write poles in non-polar form and plug them into (3)

$\displaystyle s_{i}=\Omega_{c}\ e^{j\left( \frac{\pi\left( 1+2i+N\right) }{2N}\... ...+N\right) }{2N}\right) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ i=0\cdots N-1$

then

$\displaystyle H_{a}\left( s\right) =\frac{\ K}{% {\displaystyle\prod\limits_{i... ...+N\right) }{2N}% +j\sin\frac{\pi\left( 1+2i+N\right) }{2N}\right) \right) }%$

(4)

Where

$\displaystyle \Omega_{c}=\frac{\frac{2}{T}\tan\left( \frac{\omega_{s}}{2}\right... ...frac{1}{2N}\log_{10}\left( 10^{\frac{\delta_{s}}{10}% }-1\right) \right) }}%$

and

$\displaystyle N=\left\lceil \frac{1}{2}\frac{\log_{10}\left( 10^{\frac{\delta_{... ...log _{10}\left( \tan\left( \frac{\omega_{s}}{2}\right) \right) }\right\rceil$

Now that we have found $H\left( s\right)$ we need to convert it to $H\left( z\right)$ . After finding

as shown above, we simply replace

by $\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$ . This is much simpler than the impulse invariance method. Before doing this substitution, make sure to multiply poles which are complex conjugate of each others in the denominator of

. After this, then do the above substitution

2.3 Summary of analytical derivation method

We will now make a table with the derivation equations to follow to design in either bilinear or impulse invariance. Note that the same steps are used in both designs except for step . This table make it easy to develop a program.

$\displaystyle % \begin{tabular}[c]{\vert l\vert l\vert l\vert l\vert}\hline s... ... _{s=\frac{2}{T}\frac{1-z^{-1}% }{1+z^{-1}}}$\\ \hline \end{tabular} \ \ \$

3 Numerical design examples

3.1 Example 1

Sampling frequency $F_{s}=20khz$ , passband frequency $f_{p}=2khz$ , stopband frequency $f_{s}=3khz$ , with $\delta_{p}\geq-1db$ and $\delta_{stop}% \leq-15db$

3.1.1 using impulse invariance method (using T=1)

$\displaystyle % \begin{tabular}[c]{\vert l\vert l\vert}\hline step & Impulse ... ...8-j0.523\,68\allowbreak\right) z^{-1}}$\\ \hline \end{tabular} \ \ \ \ \ \$

3.1.2 bilinear method

Using the above design table, these are the numerical values: $T=\frac {1}{F_{s}}=\frac{1}{20000}$

$\displaystyle % \begin{tabular}[c]{\vert l\vert l\vert}\hline step & Bilinear... ...{-1}}}\rightarrow TODO$\\ \hline & $\cdots$\\ \hline \end{tabular} \ \ \ \$

$\allowbreak$

4 References

Alan Oppenheim, Ronald Schafer, Digital Signal Processing, Prentice-Hall, inc. 1975, Chapter 5.
Mostafa Shiva, Electrical engineering department, California state university, Fullerton, Lecture notes, handout H.
John Proakis, Dimitris Manolakis, digital signal processing, 3rd edition, section 8.3.

me 2012-05-09