\chapter[IIR Filters]{Design of Infinite Impulse Response Filters}
\label{iir}
\index{IIR filters}

\section{Analog filters}
\index{analog filters}
In this section we review some of the most classical analog (or ``continuous time'') filters.
These are defined in frequency domain as rational transfer functions of the form:
$$H(s)=\displaystyle{ \frac{\sum_{i=0}^{m}{b_i s ^i}}{1+\sum_{i=1}^{n}{a_i s ^i}}}$$
The problem is then to determine the $a_i$ and $b_i$ coefficients or equivalently the zeros $z_i$ and poles $p_i$ of $H$ in order to meet given specifications of the (squared) magnitude response defined as: 
\begin{equation}
\label{hhh}
h^2(\omega)={\vert H(i\omega) \vert}^2={H(s)H(-s) \vert _{s=i\omega}}
\end{equation}

Thus, $h(\omega)$ is the spectrum of the output of a linear filter which admits a white noise as input.
We shall consider in this section only prototype lowpass filters, i.e., ideally we want to obtain $h(\omega) = 0$ for $\omega > \omega_c$ and $h(\omega) = 1$ for $\omega < \omega_c$.
Highpass, bandpass and stopband filters are then easily constructed by a simple change of variables.

The construction consists of finding a function $H_2$ of the complex variable $s$ such that $H_2(s) = H(s)H(-s)$ (i.e., 
is symmetric with respect to the imaginary axis and 
such that $H_2(i \omega) = h^2(\omega)$ along the imaginary axis).
Furthermore, the function $H_2(s)$ will be chosen to be rational and, 
consequently,  defined by its poles and zeros.

The transfer function of the filter, $H(s)$,  will then be 
defined by selecting all the poles and zeros which lie in the left hand
side of the complex $s$-plane.
In this way we obtain a stable and minimum phase filter.

\subsection{Butterworth Filters}
\index{analog filters!Butterworth}
\index{Butterworth filter}
The squared-magnitude response that we want to realize is given by:
\begin{equation}
\label{real}
h_n^2(\omega \vert \omega_c)=\displaystyle{\frac{1}{1+\displaystyle{{(\frac{\omega}{ \omega_c})}^{2n}}}}
\end{equation}

Here, $\omega_c$ is the cutoff frequency and $n$ the order.
A typical response can be plotted with the function {\tt buttmag} (see Figure~\ref{butt}):

The following code gives an example of the squared magnitude of a Butterworth filter of order 13 (see Figure~\ref{butt}).

\verbatok{\Diary Sanalog1.dia}
\end{verbatim}

\input{\Figdir analog1.tex}
\dessin{{\tt exec('analog1.code')} Magnitude in dB. $n=13,\omega_c=300$}
{butt}

From Figure~\ref{butt} we see that the magnitude is a decreasing function of the frequency. Also we see that $$h_n^2(\omega_c \vert \omega_c)=\frac{1}{2}$$ independently of the order $n$.

Let us  define now the stable transfer function $H(s)$ (Butterworth filter) 
by:
 $$H(s)=\frac{k_0}{\prod_{k=1}^n{(s-p_k)}}$$ where 
$$ p_{k}=e^{i\pi[1/2+(2k-1)/2n]} \qquad k=1,\ldots,n$$

This is done by the small function {\tt zpbutt} which computes the poles $p_k$ and the gain $k_0$:

For instance, with n=7 and $\omega_c = 3$, we obtain the following transfer function:

\verbatok{\Diary butt.dia}
\end{verbatim}

The poles $p_k$ of $H$ are located on the unit circle
and are symmetric with respect to 
the real axis as is illustrated in Figure~\ref{polesp}.
The figure was obtained as follows:
\verbatok{\Diary Sanalog2.dia}
\end{verbatim}

\input{\Figdir analog2.tex}
\dessin{{\tt exec('analog2.code')} Butterworth filter: pole positions. $n=13$}
{polesp}

We note the symmetry of the coefficients in $H(s)$ i.e., that $H(s)=\displaystyle{\tilde{H}(s)=s^n\,H(\frac{1}{s})}$ , which follows from the fact that for each
pole $p_k$ there corresponds the pole $\displaystyle{\frac{1}{p_k}=\overline{p_k}}$.
Also, we see that $H(-s)$ is obtained in the same way as $H(s)$ by selecting the (unstable) poles $-p_k$ instead of the $p_k$.  Since the set 
$\{ (p_k,-p_k) \qquad k=1,\ldots,n \}$ is made with the $2n$ roots of the polynomial $p(s)=1+{(-s^2)}^n$.  Thus,  we have:
$$H(s)H(-s)=\frac{1}{1+{(-s^2)}^n}$$
It follows immediately from (\ref{hhh}) that $H(s)$ realizes the response (\ref{real}).We see that $H(s)$ is obtained by a very simple spectral factorization of $p(s)$, which here can be done analytically.

\paragraph{Order determination:}
The filter order, $n$, completely specifies a Butterworth filter.  In general $n$ is determined by giving a desired attenuation $\frac{1}{A}$ at a specified ``normalized'' frequency $f=\frac{\omega_r}{\omega_c}$.
The filter order, $n$, is given by the solution of $\frac{1}{A^2}=h_n^2(f)$. We obtain immediately: 
\begin{equation}
n=\frac{\log_{10}(A^2-1)}{2\,\log_{10}(f)}
\end{equation}

\subsection{Chebyshev filters}
\index{analog filters!Chebyshev}
\index{Chebyshev filter}

The $n^{th}$ order Chebyshev polynomial $T_{n}(x)$ is defined recursively by:
$$\left\{
\begin{array}{ll}
T_{n+1}(x)=2xT_{n}(x)-T_{n-1}(x) & \\
T_{0}(x)=1 & T_1(x)=x
\end{array}
\right .$$
It may be shown that $T_{n}(x)$ is given more explicitly by:
\begin{equation}
T_{n}(x) = \left\{
\begin{array}{ll}
 \cos(n \cos^{-1}(x)) & \mbox{if $|x| < 1$} \\
 \cosh(n \cosh^{-1}(x)) & \mbox{otherwise}
\end{array}
\right . 
\label{expli}
\end{equation}
The recursive function {\tt chepol} implements this computation:
We note that the roots of $T_n$ are real and symmetric with respect
to the imaginary axis.
These polynomials are used to analytically define the squared-magnitude of the frequency response for a class of  analog filters.

\subsubsection{Type 1: Ripple in the passband}
\index{analog filters!Chebyshev!first type}
\index{Chebyshev filter!second type}

For type I Chebyshev filters the squared-magnitude function that we want to realize is:
\begin{equation}
 h_{1,n}^2(\omega \mid \omega_c , \epsilon)=\frac{1}{1+\epsilon^2 T_{n}^2(
\frac{\omega}{\omega_c})}
\label{eqch1}
\end{equation}
This function is completely specified once its three parameters $(\omega_c,\epsilon,n)$ are given.

Using {\tt chepol}, it is easy to compute the squared-magnitude response (\ref{eqch1}). The function {\tt cheb1mag} evaluates $h_{1,n}^2(\omega \mid \omega_c,\epsilon) $ for a given sample vector of $\omega 's $.
Note that for any value of $n$ one has $$h_{1,n}^2(\omega_c \mid \omega_c,\epsilon) = \frac{1}{1 + \epsilon^2}$$
The function $h_1$ is decreasing for $\omega > \omega_c$ with ``fast'' convergence to zero for ``large'' values of the order $n$.  The number of oscillations in the passband is also increasing with $n$.
If at $\omega = \omega_r > \omega_c$ , $h_{1,n}^2$ reaches the value $\frac{1}{A^2}$ then the parameters $(\omega_c,\epsilon,n)$ and $(A,\omega_r)$ are linked by the equation
$h_{1,n}^2(\omega_r | \omega_c,\epsilon) = \frac{1}{A^2}$ which may be written also as
\begin{equation}
A^2 = 1+\epsilon^2 T_{n}^{2}(\frac{\omega_r}{\omega_c})
\label{equ15}
\end{equation}
Using (\ref{expli}) this latter equation may be solved more explicitly:
$n \cosh^{-1}(f) = \cosh^{-1}(g)$ with $f = \displaystyle{\frac{\omega_r}{\omega_c}}$ and $g=\frac{A^2-1}{\epsilon}$.


Below is an example of the magnitude plotted in Figure~\ref{ch1mg}.

\verbatok{\Diary Sanalog3.dia}
\end{verbatim}

\input{\Figdir analog3.tex}
\dessin{{\tt exec('analog3.code')} Magnitude of a Type 1 Chebyshev filter}
{ch1mg}
Let us now calculate the transfer function of a Type 1 Chebyshev filter.
The transfer function is all-pole and the poles lie on an ellipse which can be determined totally by specifying the parameter $\epsilon$, the order $n$ and the cutoff frequency $\omega_c$.
The horizontal and vertical rays $a$ and $b$ of this ellipse are given by:
$a = \omega_c \displaystyle{\frac{\gamma - \gamma^{-1}}{2}}$ and 
$b = \omega_c \displaystyle{\frac{\gamma + \gamma^{-1}}{2}}$ where
$$\gamma = {(\frac{1+\sqrt{1+\epsilon^{2}}}{\epsilon})}^{1/n}$$

The poles $(p_k=\sigma_k + i\,\Omega_k , k=1,2,\ldots,n)$ are simple 
and, in order to have a stable filter, are regularly spaced over the 
left hand side of the ellipse.
The function {\tt zpch1} computes the poles (and the gain) of a Type 1 Chebyshev filter.

With the function {\tt zpch1} we can now calculate the transfer function which realizes the preceding example and recover the desired magnitude.
Compare Figures~\ref{ch1mg} and \ref{ch1mg2}, the latter figure being obtained as follows:

\verbatok{\Figdir analog4.dia}
\end{verbatim}

\input{\Figdir analog4.tex}
\dessin{{\tt exec('analog4.code')} Chebyshev filter: frequency response in magnitude}{ch1mg2}

\paragraph{Order determination:}
We have seen that the parameter ${\epsilon}$ specifies the size of the passband ripple and $\omega_c$ is the passband edge, thus,
one has $$ \frac{1}{1+\epsilon^2} \leq h_{1,n}^2(\omega \mid \omega_c) \leq 1 \qquad \mbox{for} \quad 0 \leq \omega \leq \omega_c$$
The order $n$ is generally determined by specifying a desired attenuation $\frac{1}{A}$ at a given ``normalized'' frequency $f = \frac{\omega_r}{\omega_c}$.
As in (\ref{equ15}), $n$ is given by the solution of $\frac{1}{A^2} = h_{1,n}^2(f | \omega_c,\epsilon)$:

$$n=\frac{\cosh^{-1}(\frac{\sqrt{A^2-1}}{\epsilon})}{\cosh^{-1}(f)} = \frac{\log(g+\sqrt(g^2-1))}{\log(f+\sqrt(f^2-1)}$$
where $g = \sqrt(\frac{A^2-1}{\epsilon^2})$


\subsubsection{Type 2: Ripple in the stopband}

\index{analog filters!Chebyshev!second type}
\index{Chebyshev filter!first type}

The squared-magnitude response of a Type 2 Chebyshev filter is defined by:
$$ h_{2,n}^{2}(\omega \mid \omega_r , A ) = \frac{1}{1+\frac{A^2-1}{T_{n}^{2}(\frac{\omega_r}{\omega})}}$$
Here $\omega_r$ is the passband edge and $A$ the attenuation at $\omega_r$.
The function {\tt cheb2mag} computes the squared-magnitude response.
Note that the sample vector must not include the value zero.  
Also, for any value of the order $n$ one has:
$$h_{2,n}^{2}(\omega_r \mid \omega_r,A) = \frac{1}{A^2} $$
The function is decreasing for $0 < \omega < \omega_r $ and displays ripple for $\omega > \omega_r$.

Also note that when the equation (\ref{equ15}) is satisfied both Type 1 and type 2 functions $h_{1,n}$ and $h_{2,n}$ take the same values at ${\omega_c}$ and ${\omega_r}$:
$$h_{1,n}^{2}(\omega_r | \omega_c, \epsilon) = h_{2,n}^{2}(\omega_r | \omega_r, A) = \frac{1}{A^2}$$
$$h_{2,n}^{2}(\omega_c | \omega_c, \epsilon) = h_{2,n}^{2}(\omega_c | \omega_r, A) = \frac{1}{1+\epsilon^2}$$


We can now plot for example a squared-magnitude response for the following specifications:
$1/A = 0.2 ,\omega_r = 6 ,$ and $n = 10$. The sample vector of $\omega$'s is chosen as {\tt 0:0.05:10}.
The magnitude plotted in dB is given by Figure~\ref{ch2mg}, and was 
generated by  the following code:

\verbatok{\Diary Sanalog5.dia}
\end{verbatim}

\input{\Figdir analog5.tex}
\dessin{{\tt exec('analog5.code')} Magnitude of a Type 2 Chebyshev filter}
{ch2mg}

The transfer function of a type 2 Chebyshev filter has both poles and zeros.
The zeros are imaginary and are located at
$$z_k = i \frac{\omega_r}{\cos(\frac{(2k-1)\pi}{2n})} \qquad k=1,2,\ldots,n$$
The poles $p_k = \sigma_k +i \Omega_k \qquad k=1,2,\ldots,n$ are found by solving for the singularities of the denominator of $h$.
They are defined by:
$$\sigma_k = \frac{\alpha_k}{{\alpha_k}^2+{\beta_k}^2}$$
$$\Omega_k = \frac{-\beta_k}{{\alpha_k}^2+{\beta_k}^2}$$
where
$$\alpha_k = -a\, \sin(\frac{(2k-1)\pi}{2n})$$
$$ \beta_k = b\, \cos(\frac{(2k-1)\pi}{2n})$$
and 
$$a  = \frac{\gamma - \gamma^{-1}}{2}$$
$$b  = \frac{\gamma + \gamma^{-1}}{2}$$
$$\gamma = {(A +\sqrt{A^2-1})}^{1/n}$$
The function {\tt zpch2} computes the poles and zeros of a type 2 Chebyshev filter, given its parameters $(\omega_r, A, n)$, according to the preceding formulas.

Let us consider the preceding example: we had $n=10$, $\omega_r=6$, $A=5$.

%exec('cheb2.code')

\verbatok{\Diary cheb2.dia}
\end{verbatim}

\paragraph{Order determination:}
We have seen that the parameter $A$ specifies the size of the passband ripple and $\omega_r$ is the stopband edge.  Thus, one has $$ 0 \leq h_{2,n}^{2}(\omega \mid \omega_r, A) \leq \frac{1}{A^2} \qquad \mbox{for} \quad \omega_r \leq \omega $$
The order $n$ is generally determined by specifying a desired attenuation $\frac{1}{1+\epsilon^2}$ at a given ``normalized'' frequency $f = \frac{\omega_r}{\omega_c}$ (remember that in order to define a type 2 Chebyshev filter $\omega_r$ must be given).
In a similar way as in (\ref{equ15}), $n$ is given by the solution of 
$\frac{1}{1+{\epsilon}^2} = h_{2,n}^{2}(f | \omega_r,A)$.  
Because of the symmetry in $\epsilon$ and $A$ we obtain the same solution as for type 1 filters:

$$n=\frac{\cosh^{-1}(\frac{\sqrt{A^2-1}}{\epsilon})}{\cosh^{-1}(f)} = \frac{\log(g+\sqrt(g^2-1))}{\log(f+\sqrt(f^2-1)}$$
where $g = \sqrt{\displaystyle{\frac{A^2-1}{\epsilon^2}}}$

\subsection{Elliptic filters}
\index{analog filters!elliptic}
\index{elliptic filter}
The elliptic filter presents ripple in both the passband and stopband.
It is based on the properties of the Jacobian elliptic function (see \cite{dieu},\cite{Abram}) that we briefly introduce now.
\subsubsection{Elliptic integral}
\index{elliptic integral}
Let us define for $z$ in the complex plane the function 
\begin{equation}
u(z)=\int_{0}^{z}{\frac{dt}{{(1-t^2)}^{1/2}{(1-m t^2)}^{1/2}}}
\label{ellint}
\end{equation}
where $m$ is a real constant \( 0 < m < 1 \).
(We will also use the notation $u(z,m)$ when necessary.)


We may assume that the functions $u_{1}={(1-t^2)}^{\frac{1}{2}}$ and $u_{2}={(1-m t^2)}^{\frac{1}{2}}$ are defined e.g. in the domain ${\cal D}$ made of the complex plane cut along the lines \( \{z;Re(z)=\mp 1 \) and  \(Im(z) < 0 \} \) and \( \{z;Re(z)=\mp \frac{1}{\sqrt{m}} \) and \(Im(z) < 0 \} \).
In other words we may choose for these complex functions the determination of the phase between $-\pi/2$ and $3 \pi/2$.
These functions are then completely specified in ${\cal D}$ by analytic continuation once we have fixed their values at $0$ as being $+1$.

Let us define now in the above open connected domain ${\cal D}$ the function
$$\Phi(t)=\frac{1}{u_{1}(t) u_{2}(t)}=\frac{1}{{(1-t^2)}^{\frac{1}{2}}{(1-m t^2)}^{\frac{1}{2}}}$$
and consider the path ${\cal S}$ which encircles the positive right quarter of the complex plane composed of the positive real axis, the point at infinity and the imaginary axis traversed from ${\infty}$ to zero.

As $t$ increases along the positive real axis the function $\Phi(t)$ first assumes positive real values for $0 \leq t <1$, then for $1 < t < 1/{\sqrt{m}}$, $\Phi(t)$ assumes purely imaginary values and for $t > 1/{\sqrt{m}}$, $\Phi(t)$ assumes real negative values.
At $\infty$ we have $\Phi(t)=0$ and as $t$ decreases from $\infty$ to $0$ along the imaginary axis $\Phi(t)$ assumes purely imaginary values.

Let us examine now the consequences of these properties on the function $u$.
When $z$ traces the  same path ${\cal S}$ the function $u(z)$ traces the border of a {\em rectangle} ${\cal R}_0$ with corners at $(O,K,K+iK',iK')$, each side being the image of the four parts of the strip ${\cal S}$ defined above.
This follows from the identities:

$$K(m)=\int_{0}^{1}{\Phi(t)dt}=-\int_{\frac{1}{\sqrt{m}}}^{\infty}{\Phi(t)dt}$$ and
$$iK'(m)=\int_{1}^{\frac{1}{\sqrt{m}}}{\Phi(t)dt}=\int_{0}^{\infty}{\Phi(it)dt}$$
Thus the points $(0,1,\frac{1}{\sqrt(m)},\infty)$ of the real axis are respectively mapped by $u$ into the points $(0,K,K+iK',iK')$ of the complex plane and the intervals $(0,1)$,$(1,\frac{1}{\sqrt{m}})$,$(\frac{1}{\sqrt{m}},\infty)$ are respectively mapped into the intervals $(0,K)$,$(K,K+iK')$,$(K+iK',iK')$.

It may be shown that $u$ realizes a conformal mapping of the first quadrant of the complex plane into the rectangle ${\cal R}_0$.
In particular any point $z$ with $Re(z) \geq 0$ and $Im(z) \geq 0$ admits an unique image in the rectangle.
In a similar way we can map, under  $u$, each of the four quadrants of the complex plane into a rectangle ${\cal R}$ (called the ``fundamental rectangle of periods'') with corners at
$(-K-iK',K-iK',K+iK',-K+iK')$, made of four smaller rectangles defined as ${\cal R}_0$ above and having the origin as common point.
The function $u$ defines a one-to-one correspondence of the complex plane into ${\cal R}$.

The function $u$ has been implemented as the {\tt \%asn} function.  This function may receive a real vector {\tt x} as argument.  The calculation is made componentwise.  A specific algorithm \cite{Carlson} is implemented in fortran.
We note that it sufficient to be able to compute $u(x)$ for $x \in (0,1)$ to have $u(z)$ for all nonnegative real and purely imaginary values $z$ thanks to the following changes of variable:
$$
\int_{1}^{x}{\frac{dt}{{(t^2-1)}^{1/2}{(1-m t^2)}^{1/2}}} = \int_{0}^{y}{\frac{dt}{{(1-t^2)}^{1/2}{(1-m_{1} t^2)}^{1/2}}}
$$
with $m_1=1-m$ , $y^2 = \frac{1}{m_1} \frac{x^2-1}{x^2}$ and $x \in (1,1/\sqrt{m})$
$$
\int_{\frac{1}{\sqrt{m}}}^{x}{\frac{dt}{{(t^2-1)}^{1/2}{(m t^2-1)}^{1/2}}} = \int_{0}^{y}{\frac{dt}{{(1-t^2)}^{1/2}{(1-m t^2)}^{1/2}}}
$$
with $y^2 = \frac{1}{m x^2}$ and $ x \in (1/\sqrt{m},\infty)$
$$
\int_{0}^{x}{\frac{dt}{{(1+t^2)}^{1/2}{(1+m t^2)}^{1/2}}} = \int_{0}^{y}{\frac{dt}{{(1-t^2)}^{1/2}{(1-m_{1} t^2)}^{1/2}}}
$$
with $m_1 = 1-m$ and $y^2 = \frac{x^2}{1 + x^2}$, which gives $u$ for purely imaginary values of its argument.

We can now for example use {\tt \%asn} to plot three sides of the rectangle ${\cal R}_0$ as the image by {\tt \%asn} of the three associated real intervals (see Figure~\ref{asn}). Here $m=0.8$ and we see that how a linear scale is transformed by {\tt \%asn} (parametric plot). In particular at $x=10$ the point $iK'$ is not reached (as expected).

\verbatok{\Diary Sanalog6.dia}
\end{verbatim}

\input{\Figdir analog6.tex}
\dessin{{\tt exec('analog6.code')} The rectangle ${\cal R}_0$ , image by $u$ of the positive real axis.}{asn}

The integrals $K$ and $K'$ are known as the ``complete'' elliptic integral: they may be calculated by {\tt \%asn} since we have $K(m)+iK'(m) = u(\frac{1}{\sqrt(m)},m)$.
 These numbers may also be calculated by the Arithmetic-Geometric-Mean algorithm which is implemented as the {\tt \%K} function.
Note that here $m$ can be vector-valued and {\tt \%K} computes $K(m)$ componentwise (this fact will be useful later).

\subsubsection{Elliptic function}
\index{elliptic function}

For $y$ inside the fundamental rectangle ${\cal R}$ the Jacobian elliptic function $sn$ is defined as the inverse function of $u$ i.e. we have for a fixed value of the parameter $m$:
$$u(z)=y \Leftrightarrow z=sn(y)$$

In particular we have for the corners of the rectangle ${\cal R}_0$ 
defined above:
$sn(0)=0$, $sn(K)=1$, $sn(K+iK') = \frac{1}{\sqrt{m}}$, $sn(iK') = \infty$.
In fact, the function $sn$ may be extended to the full complex plane as a meromorphic function.

Indeed,the ``symmetry principle'' states that if $f$ is an analytic function defined in an open set ${\cal D}$ whose boundary contains an {\em interval} $L$ of the complex plane and if in addition this interval is itself mapped by $f$ into another {\em interval} $L'=f(L)$ of the complex plane,
 then $f$ can be extended to the set $\sigma({\cal D})$, symmetric set of ${\cal D}$ with respect to $L$, by the formula $f(\sigma(z))=\sigma'(f(z))$ where $\sigma$ and $\sigma'$ are the symmetries with respect to $L$ and $L'$ respectively.

Since the sides of the fundamental rectangle are mapped into intervals of the real and imaginary axis of the complex plane by $sn$, it may be easily shown that the function $sn$ may be extended to the full complex plane by successive symmetries leading to a doubly periodic function with periods $4K$ and $2iK'$.

For real values of its argument, $sn(y)$ ``behaves'' like the sine function and for purely imaginary values of its argument, $sn(y)$ is purely imaginary and 
has poles at $\ldots,-3iK',-iK',iK',3iK',\ldots$.
For $y$ in the interval $(-iK',+iK')$, $sn(y)$ ``behaves'' like $i$ times the $\tan$ function and this pattern is repeated periodically.
For $m=0$, one has $K=\pi/2$, $K'=\infty$ and $sn$ coincides with the $\sin$ function.
For $m=1$, one has $K=\infty$, $K'=\pi/2$ and $sn$ coincides with the $\tanh$ function.

The function $sn$ has been implemented by the following function {\tt \%sn} which calls a fortran routine for real values of the argument, and use the addition formulas (\cite{Abram}) for imaginary values of the argument.
Note that x is allowed to be complex and vector-valued.

Let us plot for example the real and imaginary behavior of $sn$; this is done by the following commands which produce the Figures~\ref{realsn} and \ref{imagsn} and give respectively $sn(x)$ for $0 \leq x \leq 4K$ and $sn(iy)$ for $0 \leq y \leq 3K'/2$.
\verbatok{\Diary Sanalog7.dia}
\end{verbatim}

\input{\Figdir analog7.tex}
\dessin{{\tt exec('analog7.code')} Behavior of the sn function for real values}{realsn}

For imaginary values of the argument we must take care of the pole of $sn(z)$ at $z=iK'$:
\verbatok{\Diary Sanalog8.dia}
\end{verbatim}

\input{\Figdir analog8.tex}
\dessin{{\tt exec('analog8.code')} Behavior of the sn function for imaginary values}{imagsn}

\subsubsection{Squared Magnitude Response of Elliptic Filter}

The positions of poles and zeros of the $sn$ function will allow to define the squared magnitude response of a prototype lowpass elliptic filter.
The zeros of the $sn$ function are located at $2pK+2qiK'$ ,where $p$ and $q$ are arbitrary integers and its poles are located at $2pK+(2q+1)K'$.


For a fixed value of the parameter $m=m_1$, let us consider a path $\Sigma_n$ joining the points $(0,nK_1,nK_1+i{K'}_1,i{K'}_1)$ with n an {\em odd} integer and denoting $K_1 = K(m_1)$ and ${K'}_1 = K(1-m_1)$.
 From the discussion above we see that for $z \in (0,nK_1)$ the function $sn(z)$ oscillates between 0 and 1 periodically as shown in Figure~\ref{realsn}.
For $z \in (nK_1,nK+i{K'}_1)$, $sn(z)$ assumes purely imaginary values and increases in magnitude, with (real) limit values $sn(nK_1)=1$ and $sn(nK_1+i{K'}_1) = \frac{1}{\sqrt{m_1}}$.
Finally for $z \in (nK+i{K_1}',i{K_1}')$, $sn(z)$ oscillates periodically between $sn(nK_1+i{K_1}') = \frac{1}{\sqrt{m_1}}$ and $\infty$.

 For $z$ in the contour $\Sigma_n$, let us consider now the function:
\begin{equation} 
v(z)=\frac{1}{1+\epsilon^2 sn^{2}(z,m_1)}
\label{v(z)}
\end{equation}
Clearly, $v(z)$ oscillates between 1 and $\frac{1}{1+\epsilon^2}$ for $z \in (0,nK_1)$ and between 0 and $\frac{1}{1+\frac{\epsilon^2}{m_1}}$ for $z \in (nK_1+i{K_1}',i{K_1}')$.
Also, clearly $v(z)$ is a continuous function of $z$, which is real-valued for $z$ in the path $\Sigma_n$ and if we chose the parameter $m_{1} = \frac{\epsilon^2}{A^{2}-1}$ we can obtain an interesting behavior.
The function {\tt ell1mag} computes $v(z)$ for a given sample vector $z$ in the complex plane and for given parameters $\epsilon$ and $m_1$.


Now, we define the vector $z = [z_1,z_2,z_3]$ as a discretization of the path $\Sigma_n$, with $z_1$ a dicretization of $(0,nK_1)$, $z_2$ a discretization of $(nK_1,nK_1+i{k_1}')$ and $z_3$ a discretization of $(nK_1 + i{K_1}' , i{K_1}')$.
Then we can produce Figure~\ref{riri} which clearly shows the behavior of $v(z)$ for the above three parts of $z$.
\verbatok{\Figdir analog9.dia}
\end{verbatim}

\input{\Figdir analog9.tex}
\dessin{{\tt exec('analog9.code')} v(z) for z in $\Sigma_n$, with $n=9$}{riri}

Of course, many other paths than $\Sigma_n$ are possible, giving a large variety of different behaviors.
Note that the number of oscillations in both the passband and stopband is $n$, an arbitrary odd integer and the amplitudes of the oscillations are entirely specified by the choice of the parameter $m_1$.

  Thus, we have obtained a function $v$ whose behavior seems to correspond to a ``good'' squared magnitude for the prototype filter that we want to realize.
However, this function is defined along the path $\Sigma_n$ of the complex plane.
Since frequencies, $\omega$, are given as positive real values we must find a mapping of the positive real axis onto the path $\Sigma_n$.
But this is precisely done by the function $u(z)=sn^{-1}(z)$ given in (\ref{ellint}) or more generally, after scaling by the function $\alpha sn^{-1}(\frac{z}{\omega_c},m) + \beta$.
Indeed, as we have seen, such a function maps the positive real axis into the border of a rectangle ${\cal R}(\alpha,\beta,m)$.  The size of this rectangle depends on the parameter $m$ and we have to chose $m,\alpha,\beta$ such that ${\cal R}(\alpha,\beta,m) = \Sigma_n$.

To be more specific, we have to choose now the value of the parameter $m$ such that:
\begin{eqnarray}
\alpha \, sn(nK_1,m) +\beta & = & \omega_c \\
\alpha \, sn(nK_1+i{K'}_1 ,m) +\beta & = & \omega_r
\end{eqnarray}

Recalling that $sn^{-1}(1,m) = K(m)$ and $sn^{-1}(\frac{1}{\sqrt{m}},m) =K(m) + iK'(m)$ we see that to satisfy these equations we must chose 
\[ \alpha = \frac{nK_1}{K} = \frac{{K_1}'}{K'}\]
and 
\[ \beta = \left\{ \begin{array}{ll}
          0  & \mbox{if n is odd} \\
         K_1 & \mbox{if n is even}
                   \end{array}
           \right. \]

 In particular, we obtain that the parameters $n$,
$m_1 = \frac{\epsilon^2}{A^2-1}$ and $m = \frac{{\omega_c}^2}{{\omega_r}^2}$ 
of the filter cannot be chosen independently but that they must satisfy the equation:

\begin{equation}
n = \frac{{K'(m_1)}}{K(m_1)} \frac{K(m)}{K'(m)} = \frac{\chi_1}{\chi}
\label{nequal}
\end{equation}
  (We note $\chi_1 = \frac{K'(m_1)}{K(m_1)} = \frac{K(1-m_1)}{K(m_1)} $ and $\chi = \frac{K'(m)}{K(m)} = \frac{K(1-m)}{K(m)})$.
Usually $m_1$ is ``small'' (close to 0) which yields $\chi_1$ ``large'', and $m$ is ``large'' (close to 1) which yields $\chi$ ``large''.

In practice very good specifications are obtained with rather low orders.
In order to plot the frequency response magnitude of an elliptic prototype lowpass filter we can proceed as follows: first select the ripple parameters $\epsilon$ and $A$ and compute 
$ m_1 = \frac{\epsilon^2}{A^2-1}$
 and $\chi_1 = \frac{K'_1}{K1}$, then for various integer values of $n$ compute $m$ such that equation (\ref{nequal}) is satisfied or until the ratio $\frac{\omega_r}{\omega_c}$ is acceptable.

%The function {\tt findm} computes $m$ such that $\displaystyle{\frac{K(1-m)}{K(m)} = \chi}$ by 
%using a convergent iterative procedure.

%It is also possible to use the function {\tt findm}, for example, to evaluate the tradeoffs between the ripple specifications and the frequencies specifications.
 See Figure~\ref{isom.code}
\verbatok{\Figdir analog10.dia}
\end{verbatim}

\input{\Figdir analog10.tex}
\dessin{{\tt exec('analog10.code')} $\log(m)$ versus $\log(m_1)$ for order $n$ fixed}{isom.code}

Much flexibility is allowed in the choice of the parameters, provided that
equation (\ref{nequal}) is satisfied. 
The functions {\tt find\_freq} and {\tt find\_ripple} may be 
used to find the stopband edge $\omega_r$ when $\omega_c$, 
$\epsilon$, $A$, and $n$ are given and to find $\epsilon$ 
when the parameters $n$, $\omega_c$, $\omega_r$, and $A$ are given.

The following code shows how to find compatible parameters and produce 
Figure~\ref{ellmag.code}.

\verbatok{\Figdir analog11.dia}
\end{verbatim}

\input{\Figdir analog11.tex}
\dessin{{\tt exec('analog11.code')} Response of Prototype Elliptic Filter}
{ellmag.code}

\subsubsection{Construction of the elliptic filter}
The design of an elliptic filter is more complex than for the filters that we have defined until now.
First the parameters $n$, $m_1$, and $m$ which characterize the squared magnitude response cannot be chosen independently (see (\ref{nequal})). We have seen how to solve this difficulty.
Second, the squared magnitude response is not a rational function and moreover it has an infinite number of poles and zeros.

The construction is accomplished in two steps: first a transformation is made in the complex plane which maps the real axis to the imaginary axis and transforms the rectangular path $\Sigma_n$ to a rectangular path $\Sigma'_n$ in the LHS of the complex plane.
Then the elliptic function $sn$ is used to perform a transformation which maps the imaginary axis into $\Sigma'_n$.
Finally, only poles and zeros which are inside $\Sigma'_n$ are kept for the transfer function.
  
 Let us consider the pole-zero pattern of the function $v(z)$.
Clearly, the poles of $sn(z)$ become double zeros of $v(z)$ and the poles of $v(z)$ are found by solving the equation:
$$
1 + \epsilon^2 \, sn^2(z) = 0
$$

Thus the zeros of $v(z)$ in $\Sigma_n$ are located at $i\,K',i\,K'+2K,i\,K'+4K,\ldots,i\,K'+2p\,K$  and the poles of $v(z)$ in $\Sigma_n$ are located at $i\,u_0,i\,u_0+2K,i\,u_0+4K,\ldots,i\,u_0+2p\,K$ with $2p+1 = n$ and where we have noted $u_0 =sn^{-1}(\frac{i}{\epsilon},m_1)$.

Consider now the transformation $\lambda = i\,\frac{K'(m)}{K'(m_1)}u = i\,\frac{K(m)}{nK(m_1)}u$ ($n$ being given in (\ref{nequal})).
The above pole-zero pole pattern is mapped inside the LHS of the complex plane and the contour $\Sigma_n$ is mapped into $\Sigma'_n = (0,i\,K,-i\,K'+K,-K')$, and these points are respectively the image of the points $(0,i\,\omega_c,i\,\omega_r,i\,\infty)$ of the imaginary axis by the function $z \rightarrow i\,\omega_c\,sn(-iz,m)$.

The function {\tt zpell} performs these transformations and computes the poles and zeros of a prototype elliptic filter.

We illustrate the use of this function by the following example which uses the preceding numerical values of the parameters $\epsilon$, $A$, $\omega_c$, $\omega_r$.
The code produces Figure~\ref{zpellps}.
\verbatok{\Diary Sanalog12.dia}
\end{verbatim}

\input{\Figdir analog12.tex}
\dessin{{\tt exec('analog12.code')} Example of response of a filter obtained by {\tt zpell}}{zpellps}

\section{Design of IIR Filters From Analog Filters}
\label{desiir}
\index{IIR filter design}

	One way of designing IIR filters is by making
discrete approximations to analog filters.  If an approximation
method is to be used then it is desirable to verify that
important design characteristics of the analog filter are
preserved by the approximation.  A
stable, causal analog filter design should yield
a stable, causal digital filter under any approximation
technique.  Furthermore, essential characteristics of 
the analog filter should be preserved.  For example,
an analog low pass filter satisfies
certain design characteristics such as the location of
the cut-off frequency, the width of the transition band, and
the amount of error in the pass and stop bands.
The approximation of an analog filter should preserve 
these design specifications.

	One approach to approximating an analog
filter design is to sample the impulse response
of the analog filter.  This is known as the impulse invariance
approximation.  The relationship between
the analog and discrete transfer functions under this approximation is
%
\begin{equation}
H(z)|_{z=e^{sT}}=\frac{1}{T}\sum_{k=-\infty}^{\infty}H(s+j\frac{2\pi k}{T}).
\label{e.iir1}
\end{equation}
%
The approximation in (\ref{e.iir1}) takes $z=e^{sT}$.  Consequently,
the left half s-plane maps into the unit circle in the
z-plane, the right half s-plane maps outside the
unit circle, and the $j\omega$-axis in the s-plane maps 
to the unit circle in the z-plane.  Thus,
this approximation technique preserves causal, stable
filters.  However, since (\ref{e.iir1}) consists of a superposition
of shifted versions of $H(s)$ along the $j\omega$-axis,
aliasing can occur if the analog filter is not bandlimited.

	Because most analog filter design techniques
do not yield bandlimited filters aliasing is a problem.
For example, a high pass filter cannot be bandlimited.  Furthermore,
because of aliasing distortion, filter specifications pertaining
to band widths and errors are not necessarily preserved under the
impulse invariance approximation.

	In the following section two alternative
approximation techniques are discussed.  Each of these techniques avoids 
problems associated with aliasing.
\section{Approximation of Analog Filters}
\subsection{Approximation of the Derivative}

	Consider an analog filter which can be 
represented by a rational transfer function, $H(s)$,
where
%
\begin{equation}
H(s)=B(s)/A(s)
\label{e.iir2}
\end{equation}
%
and $A(s)$ and $B(s)$ are polynomial functions of $s$.
The relationship between the input, $X(s)$, and the output,
$Y(s)$, of the filter in (\ref{e.iir2}) can be expressed as
%
\begin{equation}
Y(s)=H(s)X(s)
\label{e.iir3}
\end{equation}
%
or because of the rational nature of $H(s)$
%
\begin{equation}
[\sum_{n=0}^{N}a_ns^n]Y(s)=[\sum_{m=0}^{M}b_ms^m]X(s)
\label{e.iir4}
\end{equation}
%
where the $\{a_n\}$ and the $\{b_m\}$ are the coefficients of the
polynomial functions $A(s)$ and $B(s)$, respectively.

	The relationship between the input and the output in the time domain 
is directly inferred from (\ref{e.iir4}),
%
\begin{equation}
\sum_{n=0}^{N}a_n\frac{d^n}{dt^n}y(t)=\sum_{m=0}^{M}b_m\frac{d^m}{dt^m}x(t).
\label{e.iir5}
\end{equation}
%
The differential equation in (\ref{e.iir5}) can be approximated
by using the Backward Difference Formula approximation to the 
derivative.  That is, for $T$ small we take
%
\begin{equation}
y'(t)|_{nT}\approx \frac{y(nT)-y(nT-T)}{T}.
\label{e.iir6}
\end{equation}
%
Because the operation in (\ref{e.iir6}) is linear and time-invariant
the approximation can be represented by
the z-transform,
%
\begin{equation}
Z\{y'(n)\}=(\frac{1-z^{-1}}{T})Z\{y(n)\}
\label{e.iir7}
\end{equation}
%
where $Z\{\cdot\}$ represents the z-transform operation and
$y'(n)$ and $y(n)$ are sampled sequences of the time functions
$y'(t)$ and $y(t)$, respectively.

	Higher order derivatives can be approximated by
repeated application of (\ref{e.iir6}) which in turn can be
represented with the z-transform by repeated multiplication
by the factor $(1-z^{-1})/T$.  Consequently, the result
in (\ref{e.iir5}) can be approximately represented by the z-transform
as
%
\begin{equation}
[\sum_{n=0}^{N}a_n(\frac{1-z^{-1}}{T})^n]Y(z)
=[\sum_{m=0}^{M}b_m(\frac{1-z^{-1}}{T})^m]X(z).
\label{e.iir8}
\end{equation}
%
Comparing (\ref{e.iir8}) to (\ref{e.iir4}) allows an identification of a
transform from the s-plane to the z-plane,
%
\begin{equation}
s=\frac{1-z^{-1}}{T}.
\label{e.iir9}
\end{equation}
%
Solving (\ref{e.iir9}) for $z$ yields
%
\begin{equation}
z=\frac{1}{1-sT}.
\label{e.iir10}
\end{equation}
%
which can be rewritten and evaluated for
$s=j\Omega$ as
%
\begin{equation}
z=\frac{1}{2}[1+\frac{1+j\Omega T}{1-j\Omega T}].
\label{e.iir11}
\end{equation}
%

	From (\ref{e.iir11}) it can be seen that the $j\Omega$-axis
in the s-plane maps to a circle of radius $1/2$ centered at
$1/2$ on the real axis in the z-plane.  The left half s-plane
maps to the interior of this circle and the right half s-plane
maps to the exterior of the circle.  Figure~\ref{f.iir1} illustrates
this transformation.
%
\input{\Figdir iir1.tex}
\dessin{{\tt exec('iir1.code')} Transform $s=(1-z^{-1})/T$}{f.iir1}
%

	The transform in (\ref{e.iir9}) yields stable causal discrete
filters when the analog filter is stable and causal.  However, since
the $j\Omega$-axis in the s-plane does not map to the unit circle
in the z-plane it is clear that the frequency response of the
digital filter will be a distorted version of the frequency
response of the analog filter.  This distortion may be acceptable
in the case where the frequency response of the analog
filter is bandlimited and the sampling period, $T$,
is much higher than the Nyquist rate.  Under these conditions
the transformed frequency response is concentrated on the small
circle in Figure~\ref{f.iir1} near the point $z=1$ and the 
frequency response on the unit circle is less distorted.

\subsection{Approximation of the Integral}

	An alternative approach to approximating the 
derivative of $y(t)$ is to approximate the integral of $y'(t)$.
For example if $y(t)|_{nT-T}$ is known, then $y(t)|_{nT}$ can
be approximated by the trapezoidal approximation rule
%
\begin{equation}
y(t)|_{nT}=\frac{T}{2}[y'(t)|_{nT}-y'(t)|_{nT-T}]+y(t)|_{nT-T}.
\label{e.iir12}
\end{equation}
%
Taking the z-transform of the sequences $y'(n)$ and $y(n)$
yields the relationship
%
\begin{equation}
Z\{y'(n)\}=\frac{2}{T}[\frac{1-z^{-1}}{1+z^{-1}}]Z\{y(n)\}
\label{e.iir13}
\end{equation}
%
and as before, we can substitute (\ref{e.iir13}) into (\ref{e.iir5}) and
then make a correspondence with (\ref{e.iir4}) yielding the
transform
%
\begin{equation}
s=\frac{2}{T}[\frac{1-z^{-1}}{1+z^{-1}}]
\label{e.iir14}
\end{equation}
%
between the s-plane and the z-plane.  The expression in (\ref{e.iir14}) is
known as the bilinear transform\index{bilinear transform}.

	Solving (\ref{e.iir14}) for $z$ yields
%
\begin{equation}
z=\frac{1+(sT/2)}{1-(sT/2)}
\label{e.iir15}
\end{equation}
%
and evaluating (\ref{e.iir15}) for $s=j\Omega$ gives
%
\begin{equation}
z=\frac{1+(j\Omega T/2)}{1-(j\Omega T/2)}.
\label{e.iir16}
\end{equation}
%
The expression in (\ref{e.iir16}) is an all-pass transformation
which has unit magnitude and phase which takes values from $-\pi$
to $\pi$ on the unit circle as $\Omega$ goes from $-\infty$ to
$\infty$.  The transformation in (\ref{e.iir15}) maps the left half
s-plane into the unit circle in the z-plane and
the right half s-plane outside of the unit circle
in the z-plane.  Consequently, stable causal analog filters
yield stable causal digital filters under this
transformation.  Furthermore, the $j\Omega$-axis in the
s-plane maps to the unit circle in the z-plane.
The mapping of the $j\Omega$-axis onto the unit circle is not linear
and, thus, there is a frequency warping distortion 
of the analog filter design when this 
transform is used.

	Because many filters of interest, such
as low pass, band pass, band pass, and stop band
filters have magnitudes which are piece-wise constant, frequency
warping distortion is of no consequence.  That is, the bilinear 
transformation maintains the characteristics of the
analog filter design.  However, if, in addition, the
phase of the analog filter is linear, the bilinear
transformation will destroy this property when used to obtain a
digital filter.
\section{Design of Low Pass Filters}

	For piece-wise constant specifications,	
the bilinear transform is the best of the 
three possible transforms discussed for converting
analog filter designs into digital filter designs.
Here we discuss how to use the bilinear transform
to design a standard digital low pass filter.  The next
section presents a series of other transformations
which can be used to convert a digital low pass filter into
a high pass, band pass, stop band, or another low pass filter.

	To effectively use the bilinear transform
to design digital filters, it is necessary to transform
the digital filter constraints into analog filter
constraints.  Evaluating (\ref{e.iir14}) at $z=e^{j\omega}$ yields
%
\begin{equation}
s=\frac{2j}{T}\tan(\omega/2)=\sigma+j\Omega.
\label{e.iir17}
\end{equation}
%
Thus, a digital filter constraint at $\omega$ corresponds
to an analog filter constraint at 
%
\begin{equation}
\Omega=\frac{2}{T}\tan(\omega/2).
\label{e.iir18}
\end{equation}
%
Consequently, to design a digital low pass filter with
cut-off frequency $\omega_c$ first requires
an analog low pass filter design with cut-off 
frequency
%
\begin{equation}
\Omega_c=2\tan(\omega_c/2).
\label{e.iir19}
\end{equation}
%
(where we have used (\ref{e.iir18}) with $T=1$).

	Any of the analog low pass filter design techniques
already discussed (such as the Butterworth, Chebyshev,
and elliptic filter designs) can be used to obtain
the digital low pass filter design.  The choice of model order
can be made by specifying additional constraints
on the filter design.  For example, specification of
a certain amount of attenuation at a specified frequency
in the stop band can be used to obtain the model order
of a Butterworth Filter.  Such a specification for
a  digital filter would be converted to an analog filter
specification using (\ref{e.iir18}) before designing the analog
filter.  More on filter order specification can be found in 
the section on analog filter design.
	
	An example of a typical digital low-pass filter 
design from a Chebyshev analog filter design of the first type
is as follows.
The digital low-pass filter is to have a cut-off frequency of
$\pi/2$.  This constraint is transformed to an analog constraint
using (\ref{e.iir19}).  The resulting analog constraint
takes the cut-off frequency to be $2\tan(\pi/4)=2$.
Now the function {\tt zpch1} is used to design the Chebyshev filter
of the first type of order 3 and passband ripple of .05.
Since the ripple of a Chebyshev filter is $1/(1+\epsilon^2)$
it follows that for a ripple of .05 in the passband that
$\epsilon=\sqrt{(1/.95)-1}=.22942$.  Thus, the call to the function
looks like
\verbatok{\Diary iir6.dia}
\end{verbatim}
where the transfer function {\tt hs} is calculated from the 
gain and the poles returned from the function.  The magnitude of the 
the transfer function can be plotted as follows
\verbatok{\Diary iir7.dia}
\end{verbatim}
which is displayed in Figure~\ref{f.iir2}.
%
\input{\Figdir iir2.tex}
\dessin{{\tt exec('iir2\_3.code')} Magnitude of Analog Filter}{f.iir2}
%

Now the analog low-pass filter design can be transformed to
a digital filter design using the bilinear transform as follows
\verbatok{\Diary iir8.dia}
\end{verbatim}
The result of the transform yields a filter which has a magnitude
as shown in Figure~\ref{f.iir3}.
%
\input{\Figdir iir3.tex}
\dessin{{\tt exec('iir2\_3.code')} Magnitude of Digital Filter}{f.iir3}
%
\section{Transforming Low Pass Filters}
\index{transforming low pass filters}
\label{s4.5}

	The previous section discussed the design of IIR low-pass
filters based on the bilinear transform approximation to an analog low-pass
filter.  It is possible to transform a digital low-pass filter
to a high-pass filter, band-pass filter, stop-band filter, or to
another low-pass filter by using transforms similar to the bilinear
transformation.  This section presents the necessary transformations
for each of the above filter types.  The development follows \cite{rabiner}.

	Assuming that the cut-off frequency of the original digital
low-pass filter is $\omega_c$ a new low-pass filter of cut-off
frequency $\omega_u$ can be created using the following transformation
%
\begin{equation}
z\rightarrow\frac{z-\alpha}{1-z\alpha}
\label{e.iir.a}
\end{equation}
%
where
%
\begin{equation}
\alpha=\frac{\sin[(\omega_c-\omega_u)/2]}{\sin[(\omega_c+\omega_u)/2]}.
\label{e.iir.b}
\end{equation}
%

	For a high-pass filter of cut-off frequency $\omega_u$ one can use
the transformation
%
\begin{equation}
z\rightarrow-\frac{z+\alpha}{1+z\alpha}
\label{e.iir.c}
\end{equation}
%
where
%
\begin{equation}
\alpha=-\frac{\cos[(\omega_c-\omega_u)/2]}{\cos[(\omega_c+\omega_u)/2]}.
\label{e.iir.d}
\end{equation}
%

	For a band-pass filter with $\omega_u$ and $\omega_l$
the upper and lower band edges, respectively, one would use the transformation
%
\begin{equation}
z\rightarrow-\frac{z^2-(2\alpha k/(k+1))z+((k-1)/(k+1))}{1-(2\alpha k/(k+1))z+((k-1)/(k+1))z^2}
\label{e.iir.e}
\end{equation}
%
where
%
\begin{equation}
\alpha=\frac{\cos[(\omega_u+\omega_l)/2]}{\cos[(\omega_u-\omega_l)/2]}
\label{e.iir.f}
\end{equation}
%
and
%
\begin{equation}
k=\cot[(\omega_u-\omega_l)/2]\tan(\omega_c/2).
\label{e.iir.g}
\end{equation}
%

	Finally, for a stop-band filter with $\omega_u$ and $\omega_l$
the upper and lower band edges, respectivly, the following transformation
is used
%
\begin{equation}
z\rightarrow\frac{z^2-(2\alpha/(k+1))z-((k-1)/(k+1))}{1-(2\alpha/(k+1))z-((k-1)/(k+1))z^2}
\label{e.iir.h}
\end{equation}
%
where
%
\begin{equation}
\alpha=\frac{\cos[(\omega_u+\omega_l)/2]}{\cos[(\omega_u-\omega_l)/2]}
\label{e.iir.i}
\end{equation}
%
and
%
\begin{equation}
k=\tan[(\omega_u-\omega_l)/2]\tan(\omega_c/2).
\label{e.iir.j}
\end{equation}
%
\section{How to Use the Function {\tt iir}}
\index{function syntax!iir@{\tt iir}}
\index{IIR filter design!function syntax}

	The call to the function {\tt iir} has the following syntax
\begin{verbatim}
--> [hz]=iir(n,ftype,fdesign,frq,delta)
\end{verbatim}
The argument {\tt n} is the filter order which must be a positive integer.
The argument {\tt ftype} is the filter type and can take values
{\tt 'lp'} for a low-pass filter, {\tt 'hp'} for a high-pass filter,
{\tt 'bp'} for a band-pass filter, or {\tt 'sb'} for a stop-band filter.

	The argument {\tt fdesign} indicates the type of analog filter design
technique is to be used to design the filter.  The value of {\tt fdesign}
can be {\tt 'butt'} for a Butterworth filter design, {\tt 'cheb1'}
for a Chebyshev filter design of the first type, {\tt 'cheb2'}
for a Chebyshev filter design of the second type, or {\tt 'ellip'}
for an elliptic filter design.  

	The argument {\tt frq} is a two-vector
which contains cut-off frequencies of the desired filter.  For low-pass
and high-pass filters only the first element of this vector is used.
The first element indicates the cut-off frequency of the desired filter.
The second element of this vector is used
for band-pass and stop-band filters.  This second element is the upper
band edge of the band-pass or stop-band filter, whereas the first
element of the vector is the lower band edge.

	The argument {\tt delta} is a two-vector of ripple values.
In the case of the Butterworth filter, {\tt delta} is not used.
For Chebyshev filters of the first 
type, only the first element of this vector is used and it serves
as the value of the ripple in the pass band.  Consequently,
the magnitude of a Chebyshev filter of the first type
ripples between 1 and 1-{\tt delta(1)} in the pass band.
For a Chebyshev filter of the second type only the second
element of {\tt delta} is used.  This value of {\tt delta}
is the ripple in the stop band of the filter.  Consequently,
the magnitude of a Chebyshev filter of the second type
ripples between 0 and {\tt delta(2)} in the stop band.
Finally, for the elliptic filter, both
the values of the first and second elements of the vector {\tt delta}
are used and they are
the ripple errors in the pass and stop bands, respectively.

	The output of the function, {\tt hz}, is a rational polynomial 
giving the coefficients of the desired filter.
\section{Examples}
\index{IIR filter design!examples}

	In this section we present two examples using the {\tt iir}
filter design function.  We remind the user that an 
important part of the filter design
process is that there is always a trade-off between the performance
and the expense of a filter design.  For
a filter with a small error in the pass and stop bands and with a
narrow transition band it will be necessary to implement a filter of
higher order (which requires more multiplies).  Consequently,
the filter design procedure is iterative.  The user specifies a
model order and then examines the magnitude of the
resulting filter to see if the design specifications are met.  If
specifications are not satisfied, then the user must start again
with a filter of higher model order.  Another important point to
keep in mind when using the function is that band pass and stop band
filters will generate transfer functions of twice the model order
specified.  This is due to that transformation of the prototype 
low pass filter using an all pass filter of order two 
(see Section~\ref{s4.5}).

	The first example is of a low-pass filter design using
a Chebyshev filter of the first type for the analog design.
The cut-off frequency of the digital filter is $\omega_c=.2$,
the filter order is $n=5$, and the ripple in the passband
is $\delta=.05$.  The call to the function is as follows
\verbatok{\Diary iir9.dia}
\end{verbatim}
The result of the filter design
is displayed in Figure~\ref{f.iir4}
%
\input{\Figdir iir4.tex}
\dessin{{\tt exec('iir4.code')} Digital Low-Pass Filter}{f.iir4}
%

	The second example is of a band-pass filter designed from
a third order analog elliptic filter with cut-frequencies $\omega_l=.15$
and $\omega_h=.25$ and ripples in the pass and stop bands, respectively,
as $\delta_p=.08$ and $\delta_s=.03$.  The call to $\Psi$Lab looks like
\verbatok{\Diary iir10.dia}
\end{verbatim}
and the resulting magnitude of the transfer function is illustrated in
Figure~\ref{f.iir5}.
%
\input{\Figdir iir5.tex}
\dessin{{\tt exec('iir5.code')} Digital Band-Pass Filter}{f.iir5}
%
Notice that the transfer function here is of order six whereas the
specified order was three.  For band pass and stop band filters
the user must specify a filter order of half the desired order to
obtain the desired result.

\section{Another Implementation of Digital IIR Filters}
\index{IIR filter design!alternate implementation}
\subsection{The {\tt eqiir} function}
\index{function syntax!eqiir@{\tt eqiir}}
\index{IIR filter design!function code}

The {\tt eqiir} function is an interface between $\Psi$Lab and the Fortran routine {\tt syredi} which is a modification of the well known {\tt eqiir} code \cite{ieee}. 
The {\tt eqiir} function allows one to  design four different types of filters, namely lowpass, highpass, symmetric stopband, and symmetric passband filters.  The algorithm is based on the bilinear transform of analog filters as described in the previous sections.  The filter obtained is a product of second order cells.
The order of the filter is computed automatically to meet the filter specifications.

The filter is given either by the set of its poles and zeros (output variables {\tt zpoles} and {\tt zzeros} of the {\tt eqiir} function) or equivalently by a the representation:
$$H(z)=fact \prod_{1}^{N}{n_i(z)/d_i(z)}$$
where the rational fraction $n_i(z)/d_i(z)$ is the i-th element of {\tt cells}.

\subsection{Examples}
\index{IIR filter design!examples}

\paragraph{Example 1 (Lowpass elliptic filter):}
Design of a lowpass elliptic filter with maximum ripples $\delta_p=0.02$, $\delta_s=0.001$ and cutoff frequencies $\omega_1 = 2 \frac{\pi}{10}$ and $\omega_2 = 4 \frac{\pi}{10}$.

\verbatok{\Diary eqiir1.dia}
\end{verbatim}

\paragraph{Example 2 (Lowpass Butterworth filter):}
Design of a lowpass Butterworth filter with the following specifications:\\
- 5dB passband attenuation at the normalized cutoff frequencies $.25 \pi/10$ and - 120dB attenuation at the stopband edge $0.5 \pi/10$.

\verbatok{\Diary eqiir2.dia}
\end{verbatim}

\paragraph{Example 3 (Lowpass type 1 Chebyshev filter):}
Design of a lowpass type 1 Chebyshev filter with 2dB ripple in the passband and -30 dB attenuation at the stopband edge.  The sampling frequency is assumed to be 3000Hz and the cutoff frequencies at 37.5Hz and 75Hz respectively.

\verbatok{\Diary eqiir3.dia}
\end{verbatim}

\paragraph{Example 4 (Elliptic symmetric bandpass filter):}
Design of a symmetric bandpass filter with edges at $\omega_1 = 0.251463, \omega_2 = \pi/10, \omega_3 = 2 \pi/10, \omega_4=0.773302$ and ripples in the passband and stopband given respectively by $\delta_p = 0.022763$, $\delta_s = 0.01$.

\verbatok{\Diary Seqiir4.dia}
\end{verbatim}

The preceding example shows how to compute the magnitude response by using the {\tt freq} primitive. The example is plotted in Figure~\ref{ell4}.

\input{\Figdir eqiir4.tex}
\dessin{{\tt exec('eqiir4.code')} Example of response obtained with {\tt eqiir}}{ell4}

%\section{Optimized IIR filters}