\documentclass{slides}
\usepackage{epsfig}
\newcommand{\maintitle}[1]{\begin{center}\Large #1 \end{center}}
\newcommand{\slidetitle}[1]{\begin{par}\Large #1\end{par}}

\begin{document}

\maintitle{Solving the Selesnick-Burrus Filter Design Equations}


\vfil
\begin{center}

{\large
John B. Little
}

\vfil

Department of Mathematics and Computer Science

{\small
College of the Holy Cross

Worcester, MA 01610

{\tt email: little@mathcs.holycross.edu}}
\end{center}

\begin{center}
{\large AMS Sectional Meeting}

{\large Stevens Institute of Technology}
\vfil
{\small April 29, 2001}
\end{center}

%%%%%%%%%%%%%%%%
\begin{slide}

\slidetitle{{Outline of Talk}}

We will discuss an application of computational commutative
algebra to a question in {\it signal processing}.

\begin{enumerate}
\item Background on digital filters
\item A filter design problem
\item Results
\item Methods
\end{enumerate}
\end{slide}

\begin{slide}
\slidetitle{Digital Signals}

A {\em digital} signal is a quantized function of a 
{\em discrete} variable.  
\begin{itemize}
\item only have signal values at 
$n\Delta t$, $n$ an integer, $\Delta t$ called the 
{\em sampling interval} 
\item signal values are elements of a discrete set.
\end{itemize}

For simplicity, we ignore quantization effects
and study {\em discrete-time, real-valued} signals.

Can write a signal as a real-valued function of 
an integer variable, $x[n]$, or equivalently, a 
list of real numbers indexed by the integers:
\[
 \cdots, x[-2], x[-1], x[0], x[1], x[2], \cdots
\]
\end{slide}


\begin{slide}
\slidetitle{Filters}

Signal processing operations can be described 
mathematically by means of functions $\Gamma$:
\[
\Gamma : \{ {\rm signals} \} \to \{ {\rm signals} \}.
\]

We consider only operations with the following properties

\begin{enumerate}
\item (Linearity)  If $z[n] = c x[n] + y[n]$, then
\[
\Gamma(z[n]) = c \Gamma(x[n]) + \Gamma(y[n])
\]
\item (Shift-invariance)  If $y[n] = x[n + k]$ for all 
$n$, $k$ fixed, and the outputs are
$w[n] = \Gamma(x[n])$, $z[n] = \Gamma(y[n])$, then
\[
z[n] = \Gamma(y[n]) = \Gamma(x[n+k]) = w[n+k]
\]
\end{enumerate}

A $\Gamma$ satisfying both of these properties is said
to be a {\em filter}.
\end{slide}

\begin{slide}
\slidetitle{Example -- A ``Moving Average'' Filter}

For any input signal $x[n]$, define 
\begin{eqnarray*}
y[n] &=& \Gamma(x[n])\\
     &=& (x[n] + x[n-1] + x[n-2])/3
\end{eqnarray*}

If $x[n] = \cos(n\pi/2)$, 
\[
\cdots,1,0,-1,0,1,0,-1,0,1,\cdots
\]
for instance, then $\Gamma(x[n]) = $
\[
\cdots, 0,1/3,0,-1/3,0,1/3,0,-1/3,0,\cdots
\]

It is easy to see that this satisfies the linearity
and shift-invariance properties, so this defines
an filter.


\end{slide}
\begin{slide}
\slidetitle{Impulse Response of a Filter}

A nice result about filters:  Let $\delta[n]$ be 
\[
\cdots, 0,0,1,0,0, \cdots
\]
($1$ at time $n = 0$) -- called a {\em unit
impulse} input.  $\Gamma(\delta[n])$ is called the 
{\em impulse response} of the filter.

{\em If we know $h[n] = \Gamma(\delta[n])$, then we
also know the output of the filter on any other
signal $x[n]$.}  Reason: We can write
\[
x[n] = \sum_{k=-\infty}^\infty x[k]\delta[n-k]
\]
Hence by linearity 
\[
\Gamma(x[n]) = \sum_{k=-\infty}^\infty x[k] \Gamma(\delta[n-k]).
\]
Then by shift-invariance 
\[
\Gamma(x[n]) = \sum_{k=-\infty}^\infty x[k] h[n-k].
\]
\end{slide}

\begin{slide}
\slidetitle{The Transfer Function}

Package $x[n]$ by its ``$z$-transform'', the series
\[
X(z) = \sum_{n=-\infty}^\infty x[n] z^{-n}
\]
Do same for the impulse response of a filter $h[n]$
and the output $y[n]$ for input $x[n]$:
\[
H(z) = \sum_{n=-\infty}^\infty h[n] z^{-n},\ 
Y(z) = \sum_{n=-\infty}^\infty y[n] z^{-n} 
\]

Coefficient of $z^{-n}$ in the product 
$H(z)X(z)$ is:
\[
\sum_{k=-\infty}^\infty h[k]x[n-k],
\]
same as $y[n]$ (!)  

That is: $Y(z) = H(z)X(z)$.  $H(z)$ is called the 
{\em transfer function}.
\end{slide}

\begin{slide}
\slidetitle{Filter Design}

Note:  $H(e^{i\omega})$ is the (discrete) Fourier transform 
of $h$.

Goal is to design filters to perform various operations on 
signals, by designing $H(z)$.

Typical problems: design
a {\em low-pass} or {\em high pass} filters to meet 
certain specifications.

\begin{itemize}
\item Low-pass filters:  remove high-frequency components
of signals -- typically {\em smooth out} or {\em blur} signals
(can be used to remove high-frequency ``noise'')
\item High-pass filters:  remove low-frequency components --
typically focus on fine details, rapid changes
(can be used to detect features)
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Selesnick-Burrus Design Paradigm}

{\em reference:} Selesnick and Burrus, {\sl IEEE Trans. Circuits
and Systems II}, v. 45(1), January 1998.

Given the low-pass filter $H(z)$, consider

\begin{itemize}
\item {\em magnitude response:}
$F(\omega) =  |H(e^{i\omega})|^2$
\item {\em group delay:}
$G(\omega) = {d\over d\omega}\left({\rm arg}(H(e^{i\omega}))\right)$
\end{itemize}

We would like to, simultaneously, 
\begin{enumerate}
\item make the passband {\em flat},
\item make the stopband {\em flat}, and
\item make the group delay $G(\omega)$  {\em flat}
across the passband.
\end{enumerate}


\end{slide}

\begin{slide}
\slidetitle{Specifications}

Provided $N \ge K + L + M + 1$, can find
$H(z)$ satisfying:
\begin{itemize}
\item  
$h[0] = 1$ and $F^{(2i)}(0) = 0$, for $i = 1, \ldots M$.
({\em flatness in passband}).
\item
$z = -1$ a zero of $H(z)$ of order $K$ ({\it flatness in stopband}
$\omega = \pi$),  
\item
$G^{(2j)}(0) = 0$, for $j = 1, \ldots, L$ ({\it flat group delay})
\end{itemize}

($F$ and $G$ are even functions, so odd-order derivatives at 0 
automatically zero.)
\end{slide}

\begin{slide}
\slidetitle{The System of Equations}

\begin{itemize}
\item Computing derivatives of $F$ and $G$ explicitly give
{\em polynomial equations} in the filter coefficients $h[n]$, or 
in their {\em moments}
\begin{equation}
m[k] = \sum_{n=0}^{N-1} n^k h[n], k \ge 0 
\end{equation}
\item $h[n]$ and $m[k]$ are related
by the linear equations (1), so can use either set of variables.
\item Moments are somewhat more convenient.
\item Need to understand the solutions of a rather large,
quite messy system of polynomial equations in $m[k]$.
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Form of The Equations}

\begin{itemize}
\item Flatness equations (quadratic in the $m[i]$):
\[
{2i\choose i} m[i]^2 + 2\sum_{\ell=0}^{i-1} {2i\choose \ell}(-1)^{i+\ell} m[\ell] m[2i-\ell]
\]
$ = 0$, $i = 1,\ldots, M$.
\[
\sum_{\ell=0}^j (1-{2\ell\over 2j+1}) {2j+1\choose \ell} (-1)^\ell m[\ell] m[2j+1-\ell]
\]
$= 0$, $j = 1,\ldots, L$.
\item Normalization:  $m[0]=1$
\item Zeros at $z = -1$: remainder of $H(z)$ on division by $(1+z)^K$ is
zero -- $K$ linear equations on $m[i]$ 
\item Additional equations:  Can express $m[k]$ for $k \ge N$ in terms 
of $m[0],\ldots,m[N-1]$ by solving linear equations, and substitute.
\end{itemize}

\end{slide}

\begin{slide}
\slidetitle{To Summarize}

Get  $L + M + K = N-1$ equations in $N-1$ variables $m[1], \ldots, m[N-1]$,
and we expect only finitely many solutions.  (Mainly interested in 
real ones.)

Notation:  

\begin{itemize}
\item $I_{K,L,M}$ is the ideal generated by all these, 
\item $J_{L,M}$ is the ideal generated by the quadratics 
(doesn't depend on $K$)
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{``Easy'' and ``Hard'' Cases}

Selesnick and Burris: 

\begin{itemize}
\item  ("Region I") The {\em easier} cases are 
\[
\left\lfloor {M-1\over 2}\right\rfloor \le L \le M
\]
All variables except $m[1]$ can be eliminated
using  linear algebra techniques.
\item ("Region II")
The {\em harder} cases are
\[
0 \le L \le \left\lfloor {M-1\over 2}\right\rfloor - 1
\]
\end{itemize}

Example: if $L = 0$, there are $2^M$ solutions of the
design equations in these cases(!)  (One way to see this:
the generators for $I_{K,0,M}$ 
are already ``almost'' a Gr\"obner basis for $I_{K,0,M}$
with respect to a {\em grevlex} order with
$m[N-1] > m[N-2] > \cdots > m[1]$.)
\end{slide}

\begin{slide}
\slidetitle{Some Results}


\[\# {\rm sols} = \dim {\bf C}[m[1],\ldots,m[N-1]]/I_{K,L,M}\]
in ``Region II''-- doesn't depend on $K \ge 1$.  
$L = 0$ column omitted: $2^M$.

$$\vbox{\tabskip=0pt \offinterlineskip
\halign{\strut
\vrule#&\hfil#\quad
&\vrule#&\quad\hfil#\quad
&\vrule#&\quad\hfil#\quad
&\vrule#&\quad\hfil#\quad
&\vrule#&\quad\hfil#\quad
&\vrule#&\quad\hfil#\quad
&\vrule#&\quad\hfil#\quad
&\vrule#\cr
\multispan{12}\hrulefill\cr
&M/L&&1&&2&&3&&4&&5&&\cr\multispan{12}\hrulefill\cr
& 5&& 16&&  &&  &&  &&  &\cr\multispan{12}\hrulefill\cr
& 6&& 26&&  &&  &&  &&  &\cr\multispan{12}\hrulefill\cr
& 7&& 48&&24&&  &&  &&  &\cr\multispan{12}\hrulefill\cr
& 8&& 78&&38&&  &&  &&  &\cr\multispan{12}\hrulefill\cr
& 9&&152&&66&&32&&  &&  &\cr\multispan{12}\hrulefill\cr
&10&&   &&112&&50&&  &&  &\cr\multispan{12}\hrulefill\cr
&11&&   &&  &&86&&40&&   &\cr\multispan{12}\hrulefill\cr
&12&&   &&  &&142&&62&&  &\cr\multispan{12}\hrulefill\cr
&13&&   &&  &&   &&106&&48&\cr\multispan{12}\hrulefill\cr
}}
$$

Cases $M \le 7$ check with Selesnick-Burrus.
$M\ge 8$ new; blank space below Region II boundary 
means (?)  
\end{slide}

\begin{slide}
\slidetitle{Methods Used}

\begin{itemize}
\item For the cases $M$ odd, $L = (M-3)/2$ at the ``corners''
of the Region II boundary, {\em surprisingly(?)}, a {\em lex} GB
is not too hard to compute if the variables are ordered
$m[N-1] > m[N-2] > \cdots > m[1]$.  The dimension reported is then
the degree of the univariate polynomial in $m[1]$, and this can
be solved by numerical root-finding.  
\item For some cases with very small $L$, a {\em grevlex} GB
and eigenvalue methods are effective.
\item Best general method I have found is a ``hybrid'' approach
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{The ``Best'' Solution Method}

\begin{itemize}
\item The quadratic flatness conditions are homogeneous if 
${\rm weight}(m[i]]) = i$, and equations have a very special form. 
\item In fact, the variety defined by
the quadratics, $V(J_{L,M})$, is a special type of {\em rational
scroll} -- it has a parametrization via rational (monomial) curves
with ``directrix'' a linear variety.  
\item It is always easy to 
find this parametrization (a {\em lex} GB of $J_{L,M}$
with $m[N-1] > m[N-2] > \cdots > m[1]$ is easy to compute and shows
it clearly).  
\item Substituting from this parametrization reduces the 
number of variables considerably to $M - L$.
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{``Best'' Solution Method, cont.}

\begin{itemize}
\item The equations left have quite low degrees in all variables except 
$m[1]$.
\item By {\em either} sparse resultants, or pairwise ordinary
(Sylvester) resultants, can eliminate all the variables except $m[1]$
and get a univariate polynomial in the ideal $I_{K,L,M}$.  
\item In cases I have looked at this way, the actual computation
is fast, but I don't have an algorithm for forming the 
pairwise ordinary (Sylvester) resultants (did it ``by hand''
trial and and error), and there are always
extraneous factors in the result.
\item Sparse resultants are relatively slower, and also have
large extraneous factors.
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Things That Make These Systems Hard}

\begin{itemize}
\item High degrees of univariate polynomials
  
\item For large $M,L$, the univariate polynomial has {\em gigantic}
coefficients.

\item A {\em BIG} numerical analysis problem.  Need to use very 
many decimal digits in floating-point computation to get accurate roots.  

\item In fact this is the 
current ``bottle-neck'' in the cases I have worked out.  
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Conclusions}

CAS's like Maple, etc. where 
the number of digits in floating point computation can be set
are very convenient for this!  (Arbitrary precision
arithmetic systems would also be well-adapted, of course.)

These systems of equations are a real challenge to the ``state
of the art'' in polynomial system solving.  It's possible
to ``nibble around the edges'' of Region II a bit, but 
solving them requires a lot of computer time.  A ``new idea''
is probably needed for a general solution method.
\end{slide}
\end{document}