% Copyright INRIA

%\documentstyle{book}
%\begin{document}

\def\verbatok#1{\expandafter \begin{verbatim}
\input{#1}}

\chapter{Optimization in filter design}

In this chapter, some optimization techniques are investigated for the design 
of IIR as well as FIR filters. Those techniques are particuliarly usefull 
when non standard responses are desired.

%\include{lpiir}

\label{optimized IIR filters}
\section{Optimized IIR filters}
\index{optimized!IIR filters}

In a previous chapter on the design of IIR filters, several 
methods have been investigated which make use of closed-form 
expressions for this purpose. In this section the desired 
specifications are achieved with the help of optimization 
techniques: \index{optimization} the set of convenient filter parameters 
is computed by minimizing some error criterion \cite{rabiner}. This technique 
has to be used when nonconventional filter
specifications are to be achieved.
\subsection{Minimum Lp design}
\label{minlp}
\index{minimum Lp design}
The criterion that has been retained here is the minimum ${\Bbb L}_p$ 
error\index{minimum!Lp error} and the set of parameters to be optimized, the set of poles 
and zeros the cardinal of which being specified by the user. 
The reason for such a choice of parameters is that specifications 
on the group delay\index{group delay} are much more easily 
written for this set than 
for the usual filter coefficients - especially the computations 
of gradients - One may note also that the minimum ${\Bbb L}_p$ 
error criterion admits the well-known minimum mean-square 
error\index{mean-square error} criterion as a particular case 
by setting $p$ to two.

Now, let $H(z)$ be the transfer function of the filter given as a 
cascade of $K$ second-order sections:
\begin{eqnarray}
H(z)=A\prod_{k=1}^K\frac{z^2-2r_{0k}cos\theta_{0k}z+r_{0k}^2}
{z^2-2r_{pk}cos\theta_{pk}z+r_{pk}^2}
\label{e.minlp.1}
\end{eqnarray}
The set of parameters to be optimized is then given by the following 
vector: 
\begin{eqnarray}
\phi=(r_{0k}, \theta_{0k}, r_{pk}, \theta_{pk}, A) & k=1,K
\label{e.minlp.2}
\end{eqnarray}
where index $0$ stands for zeros and index $p$ for poles, no confusion 
being to be expected with index $p$ in the ${\Bbb L}_p$ error criterion\index{error criterion}. 
Usually the specifications for a filter are given separately for 
the magnitude $|H(e^{j\omega})|$ and/or the group delay\index{group delay}  
$\tau(\omega)$; 
the corresponding expressions are:
\begin{equation}
\begin{array}{l}
|H(e^{j\omega})|\triangleq a(\phi , \omega ) \nonumber 
\end{array}
\end{equation}
\begin{eqnarray}
 = 
A\prod_{k=1}^K\frac{(1-2r_{0k}cos(\omega-\theta_{0k})+r_{0k}^2)^{1/2} 
(1-2r_{0k}cos(\omega+\theta_{0k})+r_{0k}^2)^{1/2}}
{(1-2r_{pk}cos(\omega-\theta_{pk})+r_{pk}^2)^{1/2} 
(1-2r_{pk}cos(\omega+\theta_{pk})+r_{pk}^2)^{1/2}}
\label{e.minlp.3}
\end{eqnarray}
%\end{equation}
and
\begin{eqnarray}
& &\tau(\phi,\omega)=                \nonumber           \\& & \sum_{k=1}^K
\{
\frac{1-r_{pk}cos(\omega-\theta_{pk})}
{(1-2r_{pk}cos(\omega-\theta_{pk})+r_{pk}^2)^{1/2}} + 
\frac{1-r_{pk}cos(\omega+\theta_{pk})}
{(1-2r_{pk}cos(\omega+\theta_{pk})+r_{pk}^2)^{1/2}} \nonumber \\ & &-  
\frac{1-r_{0k}cos(\omega-\theta_{0k})}
{(1-2r_{0k}cos(\omega-\theta_{0k})+r_{0k}^2)^{1/2}} - 
\frac{1-r_{0k}cos(\omega+\theta_{0k})}
{(1-2r_{0k}cos(\omega+\theta_{0k})+r_{0k}^2)^{1/2}}
\}
\end{eqnarray}
Defining the desired magnitude response\index{magnitude response} 
$a_d(\omega_j)$ and 
group delay\index{group delay} $\tau_d(\omega)$, 
the minimum ${\Bbb L}_p$-design problem 
can be formulated by mean of the following error function: 
\begin{eqnarray}
E(\phi)&=&\lambda\sum_{j=1}^Jw_a(\omega_j)[a(\phi,\omega_j)-a_d(\omega_j)]^{2p} \nonumber
\\ &+&(1-\lambda)\sum_{j=1}^Jw_{\tau}(\omega_j)[\tau(\phi,\omega_j)-
\tau_d(\omega_j)]^{2p}
\label{e.minlp.5}
\end{eqnarray}
where $w_a(\omega_j)$ and $w_{\tau}(\omega_j)$ are weighting functions 
defined on a dense grid of frequencies $\{\omega_j /0\leq\omega_j 
\leq \pi\}$
and $\lambda$ is a real number belonging to $[0,1]$ that reflects the 
importance of the specifications on the magnitude relative to that 
on the group delay\index{group delay} in a straightforward fashion.
One seek after a vector $\phi^{*}$ such that $E(\phi^{*})$ is 
minimum: this problem is easily solved in Scilab with the help 
of the function {\tt optim} the purpose of which is the resolution 
of general nonlinear optimization problems . 
We refer the reader to
the relevant documentation \cite{delebecque} for a detailed explanation 
of its use.
The {\tt optim} function needs some 
input parameters, the main of which being what is called a {\em simulator}: 
it may be given as a Scilab function and its purpose is to give the cost function\index{cost function} 
of the relevant problem as well as its gradient relative to the 
specified parameters. For the minimum ${\Bbb L}_p$ design problem, the 
simulator is named {\tt iirlp} and makes use of two other macros: 
{\tt iirmag} and {\tt iirgroup}; it 
gives $E(\phi)$ together with its gradient relative to $\phi$.

The following example will give an idea of what can be achieved 
with this optimization\index{optimization} technique: we are given a low-pass elliptic filter 
with normalized cut-off frequency 0.25, transition\index{transition} 
bandwidth equal to
0.15, ripple in the passband of 0.1 
and ripple in the stopband of 0.01 (i.e. 40dB of attenuation); 
with the Scilab function 
{\tt eqiir} we have obtained a filter of the fourth order together 
with its zeros and poles. 

\verbatok{\Diary optiir2.dia}
\end{verbatim}
Now we want to inverse this filter, that is 
we seek after a filter the magnitude reponse of which times that of 
the original elliptic filter equals one in the passband, while in the 
stopband the total attenuation is less than 80dB. This situation 
appears for example when a digital filter is needed, after analog-to-digital 
conversion, to compensate the 
deformations of a signal by an anti-aliasing analog filter.  
 The corresponding 
specifications are obtained the following way:
\verbatok{\Diary optiir1.dia}
\end{verbatim}
Although the constraints on the parameters can be easily specified, 
an unconstrained optimization procedure has been applied because the 
initial values are unknown (hence may be far away from the 
optimal solution, if any) and the order is unknown too. Instead, as 
indicated in \cite{ieee}, the 
Scilab function {\tt optim} will be applied several times and 
when some pole or zero goes outside the unit circle, it will be 
replaced by the symmetric (with respect to the unit circle) complex 
number and a new optimization runned. To see the results obtained with 
a constrained optimization, it may be interesting to run the 
following command, recalling 
that the constraints on the poles and zeros are:
\begin{equation}
\left\{
\begin{array}{ll}
0\leq r_{0k} \leq 1 & 0\leq r_{pk} \leq 1 \\
0\leq \theta_{0k} \leq 2\pi & 0\leq \theta_{pk} \leq 2\pi
\end{array}
\right.
\label{e.minlp.6}
\end{equation}
\verbatok{\Diary optiir3.dia}
\end{verbatim}
Another method to solve this problem would be to run an optimization  
with penalties on the constraints, in order to keep the poles and zeros 
inside the unit circle: we did not try it. 
Now, back to the unconstrained optimization, after several runs of {\tt optim} without constraints, an optimal solution has 
been achieved for the chosen filter order. Nevertheless it is 
seen on Figure~\ref{f.minlp.1} that this solution is not satisfactory:
%%%%%%%%%%%
\input{\Figdir optiir1.tex}
\dessin{{\tt exec('optiir.1.code')} Minimum mean-square design. Fourth order IIR filter}{f.minlp.1}
%%%%%%%%%%
Figure~\ref{f.minlp.2} shows that the product of the two magnitude 
responses is far from being equal to one in the passband and that 
the total prescribed attenuation is not achieved (the horizontal line
is at -80 dB):
%%%%%%%%%%%
\input{\Figdir optiir2.tex}
\dessin{{\tt exec('optiir.2.code')} Resulting magnitude response. Log scale}{f.minlp.2}
%%%%%%%%%%
So a second-order block has been added (four 
more parameters) to the transfer function found at the preceding step, 
leading to a filter of order six:
\verbatok{\Diary optiir4.dia}
\end{verbatim}
The same optimization procedure has been applied with this initial value, 
resulting in the following solution vector:
\verbatok{\Diary optiir5.dia}
\end{verbatim}
the desired magnitude reponse and the one achieved
with that solution appear in Figure~\ref{f.minlp.3}, while the product 
of the log-magnitude responses is in Figure~\ref{f.minlp.4}.
%%%%%%%%%%%
\input{\Figdir optiir3.tex}
\dessin{{\tt exec('optiir.3.code')} Minimum mean-square design. Sixth order IIR filter}{f.minlp.3}
%%%%%%%%%%
%%%%%%%%%%%
\input{\Figdir optiir4.tex}
\dessin{{\tt exec('optiir.4.code')} Resulting magnitude response. Log scale}{f.minlp.4}
%%%%%%%%%%
As we are interested in what happens in the passband, focusing on it 
is needed: this is done in Figure~\ref{f.minlp.5} and we see that for  
$\omega \in [0,.45]$ the ripple is equal to 0.07 dB. 
The reader may convince himself that better approximation may be 
obtained with an increase of the filter order; we mention too that 
the specification of $a_d$ at the beginning of the transition is not 
likely to be that of a real filter (it has not a monotone decreasing 
behaviour in that region !) and that a more realistic desired response 
could have been best approximated with the same number of parameters.
%%%%%%%%%%%
\input{\Figdir optiir5.tex}
\dessin{{\tt exec('optiir.5.code')} Log-magnitude response. $\omega \in [0,0.45]$}{f.minlp.5}
%%%%%%%%%%d

%\include{optfir}

\label{optimized FIR filters}
\section{Optimized FIR filters}
\index{optimized!FIR filters}

As for IIR filters, optimization techniques may be used to achieve 
particuliar specifications for FIR filters \cite{rabiner} . In this framework, 
the design problem formulation leads to a simple linear programming
\index{linear programming} 
problem, which makes this approach attractive, as opposed to nonlinear 
methods used in the case of optimization based, IIR filter design. 

As the approach in this section is based on the frequency sampling 
technique \ref{s1}, we first refer to 
the frequency response of a linear phase FIR filter as given by formula \ref{e1}. In fact, because the linear phase term can 
be ignored for the design purpose, we rather consider the real function: 
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
H^{*}(e^{j\omega}) &=& \sum_{k=0}^{N-1}H(k)S(\omega,k)
\label{e.optfir.1}
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
where 
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
S(\omega,k) &=& e^{-jk\pi/N}\frac{\sin(N\omega/2)}{\sin(\omega/2-k\pi/N)}\nonumber\\
&=&\pm e^{-jk\pi/N}\frac{\sin(N(\omega/2)-k\pi/N)}{\sin(\omega/2-k\pi/N)}
\label{e.optfir.2}
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
are the interpolating functions.
Usually in filter design, specifications are given in passbands and stopbands 
while absent in transition bands for which the width rather is given. 
For that reason, $H^{*}(e^{j\omega})$ can be written: 
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
H^{*}(e^{j\omega}) &=& B(\omega)+\sum_{i=1}^pT_iA_i(\omega)
\label{e.optfir.3}
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
where $B(\omega)$ gives the contribution to $H^{*}(e^{j\omega})$ of all the 
fixed frequency samples (that is those in the passbands and the stopbands) 
and the $A_i(\omega)$ the contribution of all the unconstrained samples 
(that is the ones in the transitions) with respective magnitude $T_i$, these 
being to be optimized. In the sequel, the union of passbands will be 
called region 1,noted $R_1$ and that of passbands region 2, noted $R_2$. 
We now want, for a fixed approximation error in $R_1$, to find the 
linear phase FIR filter giving the maximum attenuation in $R_2$ - note 
that the converse is another possible issue - This can be formulated 
as follows: 

For some fixed $\epsilon$ and desired frequency response $H_{d}(e^{j\omega})$, 
find the set of ${T_i}$, solution of: 
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
%\stackrel{\textstyle min}{{T_i}}\; \stackrel{\textstyle max}{\omega \in R_2} 
%\; \stackrel{\textstyle |H^{*}(e^{j\omega}) - H_{d}(e^{j\omega})|}{}
\min_{T_i}\max_{\omega \in R_2}|H^{*}(e^{j\omega}) - H_{d}(e^{j\omega})|
\label{e.optfir.4}
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
and subject to the constraints: 
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
|H^{*}(e^{j\omega}) - H_{d}(e^{j\omega})| \leq \epsilon
\label{e.optfir.5}
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
Because of the linearity of $H_{d}(e^{j\omega})$ with respect to the $T_i$, 
we are led to a linear programming problem, the cost function to minimize 
being the maximum attenuation in $R_2$, which will be denoted by $T_{p+1}$ 
for notational convenience and the constraints on the $T_i , i=1 \ldots p$, 
being formulated with the help of (\ref{e.optfir.5}). The optimization 
problem may be formulated as to find the set of $T_i's$ such that 
$T_{p+1}$ is minimum subject to the constraints:
%%%%%%%%%%%%%%%%%%%%
\[
\begin{array}{c}
\left.
\begin{array}{cc}
\sum_{i=1}^pT_iA_i(\omega) & \leq  \epsilon-B(\omega)+H_{d}(e^{j\omega})\\ \\
-\sum_{i=1}^pT_iA_i(\omega) & \leq  \epsilon+B(\omega)-H_{d}(e^{j\omega})
\end{array}
\right \} \omega \in R_1\\ \\
\left .
\begin{array}{cc} 
\sum_{i=1}^pT_iA_i(\omega)-T_{p+1} & \leq  -B(\omega)+H_{d}(e^{j\omega})\\ \\
-\sum_{i=1}^pT_iA_i(\omega)-T_{p+1} & \leq  B(\omega)-H_{d}(e^{j\omega})
\end{array}
\right \} \omega \in R_2\\
\end{array}
\]
%%%%%%%%%%%%%%%%%%%%
Now the problem is in a suitable form for solution via the classical 
$Simplex$ $Method$. Let us mention too that, with minor changes, the 
problem -and the associated macro- may be stated as to find the filter 
with minimum ripple in 
the passbands for a given attenuation in the stopbands.
\index{simplex method}
In the following, only an example of the standard lowpass filter type
is treated although any other frequency 
response can be approximated with this method. 

\paragraph{Example 1}: figure \ref{optfir.fig1} 
 shows the frequency response 
of a lowpass type 1 filter with the following specifications: $n$=64;
cut-off frequency, $f_c$=.15; $\epsilon=0.01$; three samples in the transition.
%%%%%%%%%%%
\input{\Figdir optfir1.tex}
\dessin{{\tt exec('optfir1.code')} Linear programming design. 64-point lowpass FIR filter}{optfir.fig1}
%%%%%%%%%%
Figure \ref{optfir.fig2} shows the log magnitude squared of the initial filter 
defined by the rectangular window and the optimized filter. 
$0.28$ and $0.30$ $\epsilon=0.001$ and three samples in each transition.
%%%%%%%%%%%
\input{\Figdir optfir2.tex}
\dessin{{\tt exec('optfir2.code')} Linear programming design. }{optfir.fig2}
%%%%%%%%%%

 
%\end{document}