compute the numerator variables. Therefore most of the iterations solve a
smaller linear system than the methods involved up to this point.

####  III-B11 Steiglitz-McBride iterative algorithm

In 1965 K. Steiglitz and L. McBride presented an algorithm [34, 33] that has
become quite popular in statistics and engineering applications. The
Steiglitz-McBride method (SMB) considers the problem of deriving a transfer
function for either an analog or digital system from their input and output
data; in essence it is a time-domain method. Therefore it is mentioned in this
work for completeness as it closely relates to the methods by Levy, SK and
Sid-Ahmed, yet it is far better known and understood.

The derivation of the SMB method follows closely that of SK. In the Z-domain,
the transfer function of a digital system is defined by

$H(z)=\frac{B(z)}{A(z)}=\frac{b_{0}+b_{1}z^{-1}+\ldots+b_{N}z^{-N}}{1+a_{1}z_{-1}+\ldots+a_{N}z^{-N}}$

Furthermore

$Y(z)=H(z)X(z)=\frac{B(z)}{A(z)}X(z)$

Steiglitz and McBride define the following problem,

$\min\;\varepsilon_{s}=\sum_{i}\mathcal{E}_{i}(z)^{2}=\frac{1}{2\pi
j}\oint\left|X(z)\frac{B(z)}{A(z)}-D(z)\right|^{2}\frac{dz}{z}$ (95)

where $X(z)=\sum_{j}x_{j}z^{-j}$ and $D(z)=\sum_{j}d_{j}z^{-j}$ represent the
z-transforms of the input and desired signals respectively. Equation (95) is
the familiar nonlinear solution error function expressed in the Z-domain.
Steiglitz and McBride realized the complexity of such function and proposed
the iterative solution (95) using a simpler problem defined by

$\min\;\varepsilon_{e}=\sum_{i}\mathcal{E}_{i}(z)^{2}=\frac{1}{2\pi
j}\oint\left|X(z)B(z)-D(z)A(z)\right|^{2}\frac{dz}{z}$ (96)

This linearized error function is the familiar equation error in the Z-domain.
Steiglitz and McBride proposed a two-mode iterative approach. The SMB Mode 1
iteration is similar to the SK method, in that at the $k$-th iteration a
linearized error criterion based on (96) is used,

$\begin{split}\mathcal{E}_{k}(z)&=\frac{B_{k}(z)}{A_{k-1}(z)}X(z)-\frac{A_{k}(z)}{A_{k-1}(z)}D(z)\\\
&=W_{k}(z)\left[B_{k}(z)X(z)-A_{k}(z)D(z)\right]\end{split}$ (97)

where

$W_{k}(z)=\frac{1}{A_{k-1}(z)}$

Their derivation444For more details the reader should refer to [34, 33]. leads
to the familiar linear system

$\mathbf{C}\boldsymbol{x}=\boldsymbol{y}$

with the following vector definitions

$\boldsymbol{x}=\left\\{\begin{array}[]{c}b_{0}\\\ \vdots\\\ b_{N}\\\ a_{1}\\\
\vdots\\\ a_{N}\end{array}\right\\}\hskip
42.67912pt\boldsymbol{q}_{j}=\left\\{\begin{array}[]{c}x_{j}\\\ \vdots\\\
x_{j-N+1}\\\ d_{j-1}\\\ \vdots\\\ d_{j-N}\end{array}\right\\}$

The vector $\boldsymbol{q}_{j}$ is referred to as the input-output vector.
Then

$\displaystyle\mathbf{C}$ $\displaystyle=$
$\displaystyle\sum_{j}\boldsymbol{q}_{j}\boldsymbol{q}_{j}^{T}$
$\displaystyle\boldsymbol{y}$ $\displaystyle=$
$\displaystyle\sum_{j}d_{j}\boldsymbol{q}_{j}$

SMB Mode 2 is an attempt at reducing further the error once Mode 1 produces an
estimate close enough to the actual solution. The idea behind Mode 2 is to
consider the solution error defined by (95) and equate its partial derivatives
with respect to the coefficients to zero. Steiglitz and McBride showed [34,
33] that this could be attained by defining a new vector

$\boldsymbol{r}_{j}=\left\\{\begin{array}[]{c}x_{j}\\\ \vdots\\\ x_{j-N+1}\\\
y_{j-1}\\\ \vdots\\\ y_{j-N}\end{array}\right\\}$

Then

$\displaystyle\mathbf{C}$ $\displaystyle=$
$\displaystyle\sum_{j}\boldsymbol{r}_{j}\boldsymbol{q}_{j}^{T}$
$\displaystyle\boldsymbol{y}$ $\displaystyle=$
$\displaystyle\sum_{j}d_{j}\boldsymbol{r}_{j}$

The main diference between Mode 1 and Mode 2 is the fact that Mode 1 uses the
desired values to compute its vectors and matrices, whereas Mode 2 uses the
actual output values from the filter. The rationale behind this is that at the
beggining, the output function $y(t)$ is not accurate, so the desired function
provides better data for computations. On the other hand, Mode 1 does not
really solve the desired problem. Once Mode 1 is deemed to have reached the
vicinity of the solution, one can use true partial derivatives to compute the
gradient and find the actual solution; this is what Mode 2 does.

It has been claimed that under certain conditions the Steiglitz-McBride
algorithm converges. However no guarantee of global convergence exists. A more
thorough discussion of the Steiglitz-McBride algorithm and its relationships
to other parameter estimation algorithms (such as the Iterative Quadratic
Maximum Likelihood algorithm, or IQML) are found in [62, 63, 64].

####  III-B12 Jackson’s method

The following is a recent approach (from 2008) by Leland Jackson [36] based in
the frequency domain. Consider vectors $\boldsymbol{a}\in\mathbb{R}^{N}$ and
$\boldsymbol{b}\in\mathbb{R}^{M}$ such that

$H(\omega)=\frac{B(\omega)}{A(\omega)}$

where $H(\omega),B(\omega),A(\omega)$ are the Fourier transforms of
$\boldsymbol{h},\boldsymbol{b}$ and $\boldsymbol{a}$ respectively. For a
discrete frequency set one can describe Fourier transform vectors
$\boldsymbol{B}=\mathbf{W}_{b}\boldsymbol{b}$ and
$\boldsymbol{A}=\mathbf{W}_{a}\boldsymbol{a}$ (where
$\mathbf{W}_{b},\mathbf{W}_{a}$ correspond to the discrete Fourier kernels for
$\boldsymbol{b},\boldsymbol{a}$ respectively). Define

$H_{a}(\omega_{k})=\frac{1}{A(\omega_{k})}$

In vector notation, let
$\mathbf{D}_{a}=\mbox{diag}(\boldsymbol{H}_{a})=\mbox{diag}(1/\boldsymbol{A})$.
Then

$H(\omega)=\frac{B(\omega)}{A(\omega)}=H_{a}(\omega)B(\omega)\Rightarrow\boldsymbol{H}=\mathbf{D}_{a}\boldsymbol{B}$
(98)

Let $H_{d}(\omega)$ be the desired complex frequency response. Define
$\mathbf{D}_{d}=\mbox{diag}(\boldsymbol{H}_{d})$. Then one wants to solve

$\min\;\boldsymbol{E}^{*}\boldsymbol{E}=\lVert\boldsymbol{E}\rVert_{2}^{2}$

where $\boldsymbol{E}=\boldsymbol{H}-\boldsymbol{H}_{d}$. From (98) one can
write $\boldsymbol{H}=\boldsymbol{H}_{d}+\boldsymbol{E}$ as

$\boldsymbol{H}=\mathbf{D}_{a}\boldsymbol{B}=\mathbf{D}_{a}\mathbf{W}_{b}\boldsymbol{b}$
(99)

Therefore

$\boldsymbol{H}_{d}=\boldsymbol{H}-\boldsymbol{E}=\mathbf{D}_{a}\mathbf{W}_{b}\boldsymbol{b}-\boldsymbol{E}$
(100)

Solving (100) for $\boldsymbol{b}$ one gets

$\boldsymbol{b}=(\mathbf{D}_{a}\mathbf{W}_{b})\backslash\boldsymbol{H}_{d}$
(101)

Also,

$\boldsymbol{H}_{d}=\mathbf{D}_{d}\hat{\boldsymbol{I}}=\mathbf{D}_{d}\mathbf{D}_{a}\boldsymbol{A}=\mathbf{D}_{a}\mathbf{D}_{d}\boldsymbol{A}=\mathbf{D}_{a}\mathbf{D}_{d}\mathbf{W}_{a}\boldsymbol{a}$

where $\hat{\boldsymbol{I}}$ is a unit column vector. Therefore

$\boldsymbol{H}-\boldsymbol{E}=\boldsymbol{H}_{d}=\mathbf{D}_{a}\mathbf{D}_{d}\mathbf{W}_{a}\boldsymbol{a}$

From (100) we get

$\mathbf{D}_{a}\mathbf{W}_{b}\boldsymbol{b}-\boldsymbol{E}=\mathbf{D}_{a}\mathbf{D}_{d}\mathbf{W}_{a}\boldsymbol{a}$

or

$\mathbf{D}_{a}\mathbf{D}_{d}\mathbf{W}_{a}\boldsymbol{a}+\boldsymbol{E}=\mathbf{D}_{a}\mathbf{W}_{b}\boldsymbol{b}$

which in a least squares sense results in

$\boldsymbol{a}=(\mathbf{D}_{a}\mathbf{D}_{d}\mathbf{W}_{a})\backslash(\mathbf{D}_{a}\mathbf{W}_{b}\boldsymbol{b})$
(102)

From (101) one gets

$\boldsymbol{a}=(\mathbf{D}_{a}\mathbf{D}_{d}\mathbf{W}_{a})\backslash(\mathbf{D}_{a}\mathbf{W}_{b}[(\mathbf{D}_{a}\mathbf{W}_{b})\backslash\boldsymbol{H}_{d}])$

As a summary, at the $i$-th iteration one can write (100) and (102) as
follows,

$\displaystyle\boldsymbol{b}_{i}$
$\displaystyle=(\mbox{diag}(1/\boldsymbol{A}_{i-1})\mathbf{W}_{b})\backslash\boldsymbol{H}_{d}$
$\displaystyle\boldsymbol{a}_{i}$
$\displaystyle=(\mbox{diag}(1/\boldsymbol{A}_{i-1})\mbox{diag}(\boldsymbol{H}_{d})\mathbf{W}_{a})\backslash(\mbox{diag}(1/\boldsymbol{A}_{i-1})\mathbf{W}_{b}\boldsymbol{b}_{i})$

####  III-B13 Soewito’s quasilinearization method

Consider the equation error residual function

$\displaystyle\begin{split}e(\omega_{k})&=B(\omega_{k})-D(\omega_{k})\cdot
A(\omega_{k})\\\
&=\sum_{n=0}^{M}b_{n}e^{-j\omega_{k}n}-D(\omega_{k})\cdot\left(1+\sum_{n=1}^{N}a_{n}e^{-j\omega_{k}n}\right)\\\
&=b_{0}+b_{1}e^{-j\omega_{k}}+\cdots+b_{M}e^{-j\omega_{k}M}\cdots\\\ &\hskip
28.45274pt-D_{k}-D_{k}a_{1}e^{-j\omega_{k}}-\cdots-
D_{k}a_{N}e^{-j\omega_{k}N}\\\ &=\left(b_{0}+\cdots
b_{M}e^{-j\omega_{k}M}\right)-\cdots\\\ &\hskip
28.45274ptD_{k}\left(a_{1}e^{-j\omega_{k}}+\cdots
a_{N}e^{-j\omega_{k}N}\right)-D_{k}\end{split}$

with $D_{k}=D(\omega_{k})$. The last equation indicates that one can represent
the equation error in matrix form as follows,

$\boldsymbol{e}=\mathbf{F}\boldsymbol{h}-\boldsymbol{D}$

where $\mathbf{F}=[\mathbf{F}_{1}\quad\mathbf{F}_{2}]$ and

$\mathbf{F}_{1}=\left[\begin{array}[]{ccccccc}1&e^{-j\omega_{0}}&\cdots&e^{-j\omega_{0}M}\\\
\vdots&\vdots&&\vdots\\\
1&e^{-j\omega_{L-1}}&\cdots&e^{-j\omega_{L-1}M}\end{array}\right]$
$\mathbf{F}_{2}=\left[\begin{array}[]{ccccccc}-D_{0}e^{-j\omega_{0}}&\cdots&-D_{0}e^{-j\omega_{0}N}\\\
\vdots&&\vdots\\\
-D_{L-1}e^{-j\omega_{L-1}}&\cdots&-D_{L-1}e^{-j\omega_{L-1}N}\end{array}\right]$

and

$\boldsymbol{h}=\left[\begin{array}[]{c}b_{0}\\\ b_{1}\\\ \vdots\\\ b_{M}\\\
a_{1}\\\ \vdots\\\
a_{N}\end{array}\right]\quad\mbox{and}\quad\boldsymbol{D}=\left[\begin{array}[]{c}D_{0}\\\
\vdots\\\ D_{L-1}\end{array}\right]$

Consider now the solution error residual function

$\displaystyle s(\omega_{k})$
$\displaystyle=H(\omega_{k})-D(\omega_{k})=\frac{B(\omega_{k})}{A(\omega_{k})}-D(\omega_{k})$
$\displaystyle=\frac{1}{A(\omega_{k})}\left[B(\omega_{k})-D(\omega_{k})\cdot
A(\omega_{k})\right]$ $\displaystyle=W(\omega_{k})e(\omega_{k})$

Therefore one can write the solution error in matrix form as follows

$\boldsymbol{s}=\mathbf{W}(\mathbf{F}\boldsymbol{h}-\boldsymbol{D})$ (103)

where $\mathbf{W}$ is a diagonal matrix with $\frac{1}{A(\omega)}$ in its
diagonal. From (103) the least-squared solution error
$\varepsilon_{s}=\boldsymbol{s}^{*}\boldsymbol{s}$ can be minimized by

$\boldsymbol{h}=(\mathbf{F}^{*}\mathbf{W}^{2}\mathbf{F})^{-1}\mathbf{F}^{*}\mathbf{W}^{2}\boldsymbol{D}$
(104)

From (104) an iteration555Soewito refers to this expression as the Steiglitz-
McBride Mode-1 in frequency domain. could be defined as follows

$\boldsymbol{h}_{i+1}=(\mathbf{F}^{*}\mathbf{W}_{i}^{2}\mathbf{F})^{-1}\mathbf{F}^{*}\mathbf{W}_{i}^{2}\boldsymbol{D}$

by setting the weights $\mathbf{W}$ in (103) equal to $A_{k}(\omega)$, the
Fourier transform of the current solution for $\boldsymbol{a}$.

A more formal approach to minimizing $\varepsilon_{s}$ consists in using a
gradient method (these approaches are often referred to as Newton-like
methods). First one needs to compute the Jacobian matrix $\mathbf{J}$ of
$\boldsymbol{s}$, where the $pq$-th term of $\mathbf{J}$ is given by
$\mathbf{J}_{pq}=\frac{\partial\boldsymbol{s}_{p}}{\partial\boldsymbol{h}_{q}}$
with $\boldsymbol{s}$ as defined in (103). Note that the $p$-th element of
$\boldsymbol{s}$ is given by

$s_{p}=H_{p}-D_{p}=\frac{B_{p}}{A_{p}}-D_{p}$

For simplicity one can consider these reduced form expressions for the
independent components of $\boldsymbol{h}$,

$\displaystyle\frac{\partial s_{p}}{\partial b_{q}}$
$\displaystyle=\frac{1}{A_{p}}\frac{\partial}{\partial
b_{q}}\sum_{n=0}^{M}b_{n}e^{-j\omega_{p}n}=W_{p}e^{-j\omega_{p}q}$
$\displaystyle\frac{\partial s_{p}}{\partial a_{q}}$
$\displaystyle=B_{p}\frac{\partial}{\partial
a_{q}}\frac{1}{A_{p}}=\frac{-B_{p}}{A_{p}^{2}}\frac{\partial}{\partial
a_{q}}\left(1+\sum_{n=1}^{N}a_{n}e^{-j\omega_{p}n}\right)$
$\displaystyle=\frac{-1}{A_{p}}\cdot\frac{B_{p}}{A_{p}}\cdot
e^{-j\omega_{p}q}=-W_{p}H_{p}e^{-j\omega_{p}q}$

Therefore on can express the Jacobian $\mathbf{J}$ as follows,

$\mathbf{J}=\mathbf{W}\mathbf{G}$ (105)

where $\mathbf{G}=[\mathbf{G}_{1}\quad\mathbf{G}_{2}]$ and

$\mathbf{G}_{1}=\left[\begin{array}[]{cccc}1&e^{-j\omega_{0}}&\cdots&e^{-j\omega_{0}M}\\\
\vdots&\vdots&&\vdots\\\
1&e^{-j\omega_{L-1}}&\cdots&e^{-j\omega_{L-1}M}\end{array}\right]$
$\mathbf{G}_{2}=\left[\begin{array}[]{ccc}-H_{0}e^{-j\omega_{0}}&\cdots&-H_{0}e^{-j\omega_{0}N}\\\
\vdots&&\vdots\\\
-H_{L-1}e^{-j\omega_{L-1}}&\cdots&-H_{L-1}e^{-j\omega_{L-1}N}\end{array}\right]$

Consider the solution error least-squares problem given by

$\min_{\boldsymbol{h}}f(\boldsymbol{h})=\boldsymbol{s}^{T}\boldsymbol{s}$

where $\boldsymbol{s}$ is the solution error residual vector as defined in
(103) and depends on $\boldsymbol{h}$. It can be shown [21, pp. 219] that the
gradient of the squared error $f(\boldsymbol{h})$ (namely
$\nabla\boldsymbol{f}$) is given by

$\nabla\boldsymbol{f}=\mathbf{J}^{*}\boldsymbol{s}$ (106)

A necessary condition for a vector $\boldsymbol{h}$ to be a local minimizer of
$f(\boldsymbol{h})$ is that the gradient $\nabla\boldsymbol{f}$ be zero at
such vector. With this in mind and combining (103) and (105) in (106) one gets

$\nabla\boldsymbol{f}=\mathbf{G}^{*}\mathbf{W}^{2}(\mathbf{F}\boldsymbol{h}-\boldsymbol{D})=\mathbf{0}$
(107)

Solving the system (107) gives

$\boldsymbol{h}=(\mathbf{G}^{*}\mathbf{W}^{2}\mathbf{F})^{-1}\mathbf{G}^{*}\mathbf{W}^{2}\boldsymbol{D}$

An iteration can be defined as follows666Soewito refers to this expression as
the Steiglitz-McBride Mode-2 in frequency domain. Compare to the Mode-1
expression and the use of $G_{i}$ instead of $F$.

$\boldsymbol{h}_{i+1}=(\mathbf{G}_{i}^{*}\mathbf{W}_{i}^{2}\mathbf{F})^{-1}\mathbf{G}_{i}^{*}\mathbf{W}_{i}^{2}\boldsymbol{D}$
(108)

where matrices $\mathbf{W}$ and $\mathbf{G}$ reflect their dependency on
current values of $\boldsymbol{a}$ and $\boldsymbol{b}$.

Atmadji Soewito [35] expanded the method of quasilinearization of Bellman and
Kalaba [65] to the design of IIR filters. To understand his method consider
the first order of Taylor’s expansion near $H_{i}(z)$, given by

$\displaystyle H_{i+1}(z)$ $\displaystyle=H_{i}(z)+$
$\displaystyle\frac{[B_{i+1}(z)-B_{i}(z)]A_{i}(z)-[A_{i+1}(z)-A_{i}(z)]B_{i}(z)}{A_{i}^{2}(z)}$
$\displaystyle=H_{i}(z)+$
$\displaystyle\frac{B_{i+1}(z)-B_{i}(z)}{A_{i}(z)}-\frac{B_{i}(z)[A_{i+1}(z)-A_{i}(z)]}{A_{i}^{2}(z)}$

Using the last result in the solution error residual function $s(\omega)$ and
applying simplification leads to

$\begin{split}s(\omega)&=\frac{B_{i+1}(\omega)}{A_{i}(\omega)}-\frac{H_{i}(\omega)A_{i+1}(\omega)}{A_{i}(\omega)}+\frac{B_{i}(\omega)}{A_{i}(\omega)}-D(\omega)\\\
&=\frac{1}{A_{i}(\omega)}[B_{i+1}(\omega)-H_{i}(\omega)A_{i+1}(\omega)+B_{i}(\omega)\ldots\\\
&\hskip 56.9055pt-A_{i}(\omega)D(\omega)]\end{split}$ (109)

Equation (109) can be expressed (dropping the use of $\omega$ for simplicity)
as

$s=W\Big{(}\big{(}\left[B_{i+1}-H_{i}(A_{i+1}-1)\right]-H_{i}\big{)}+\big{(}\left[B_{i}-D(A_{i}-1)\right]-D\big{)}\Big{)}$
(110)

One can recognize the two terms in brackets as
$\mathbf{G}\boldsymbol{h}_{i+1}$ and $\mathbf{F}\boldsymbol{h}_{i}$
respectively. Therefore (110) can be represented in matrix notation as
follows,

$\boldsymbol{s}=\mathbf{W}[\mathbf{G}\boldsymbol{h}_{i+1}-(\boldsymbol{D}+\boldsymbol{H}_{i}-\mathbf{F}\boldsymbol{h}_{i})]$
(111)

with $\boldsymbol{H}=[H_{0},H_{1},\cdots,H_{L-1}]^{T}$. Therefore one can
minimize $\boldsymbol{s}^{T}\boldsymbol{s}$ from (111) with

$\boldsymbol{h}_{i+1}=(\mathbf{G}_{i}^{*}\mathbf{W}_{i}^{2}\mathbf{G}_{i})^{-1}\mathbf{G}_{i}^{*}\mathbf{W}_{i}^{2}(\boldsymbol{D}+\boldsymbol{H}_{i}-\mathbf{F}\boldsymbol{h}_{i})$
(112)

since all the terms inside the parenthesis in (112) are constant at the
$(i+1)$-th iteration. In a sense, (112) is similar to (108), where the desired
function is updated from iteration to iteration as in (112).

It is important to note that any of the three algorithms can be modified to
solve a weighted $l_{2}$ IIR approximation using a weighting function
$W(\omega)$ by defining

$V(\omega)=\frac{W(\omega)}{A(\omega)}$ (113)

Taking (113) into account, the following is a summary of the three different
updates discussed so far:

SMB Frequency Mode-1: $\displaystyle\hskip
8.5359pt\boldsymbol{h}_{i+1}=(\mathbf{F}^{*}\mathbf{V}_{i}^{2}\mathbf{F})^{-1}\mathbf{F}^{*}\mathbf{V}_{i}^{2}\boldsymbol{D}$
SMB Frequency Mode-2: $\displaystyle\hskip
8.5359pt\boldsymbol{h}_{i+1}=(\mathbf{G}_{i}^{*}\mathbf{V}_{i}^{2}\mathbf{F})^{-1}\mathbf{G}_{i}^{*}\mathbf{V}_{i}^{2}\boldsymbol{D}$
Soewito’s quasilinearization: $\displaystyle\hskip
8.5359pt\boldsymbol{h}_{i+1}=(\mathbf{G}_{i}^{*}\mathbf{V}_{i}^{2}\mathbf{G}_{i})^{-1}\ast\ldots$
$\displaystyle\hskip
51.21504pt\mathbf{G}_{i}^{*}\mathbf{V}_{i}^{2}(\boldsymbol{D}+\boldsymbol{H}_{i}-\mathbf{F}\boldsymbol{h}_{i})$

###  III-C $l_{p}$ approximation

Infinite Impulse Response (IIR) filters are important tools in signal
processing. The flexibility they offer with the use of poles and zeros allows
for relatively small filters meeting specifications that would require
somewhat larger FIR filters. Therefore designing IIR filters in an efficient
and robust manner is an inportant problem.

This section covers the design of a number of important $l_{p}$ IIR problems.
The methods proposed are consistent with the methods presented for FIR
filters, allowing one to build up on the lessons learned from FIR design
problems. The complex $l_{p}$ IIR problem is first presented in Section
III-C1, being an essential tool for other relevant problems. The $l_{p}$
frequency-dependent IIR problem is also introduced in Section III-C1. While
the frequency-dependent formulation might not be practical in itself as a
filter design formulation, it is fundamental for the more relevant magnitude
$l_{p}$ IIR filter design problem, presented in Section III-C2.

Some complications appear when designing IIR filters, among which the
intrinsic least squares solving step clearly arises from the rest. Being a
nonlinear problem, special handling of this step is required. It was detemined
after thorough experimentation that the quasilinearization method of Soewito
presented in Section III-B13 can be employed successfully to handle this
issue.

Figure 20: Block diagram for complex $l_{p}$ IIR algorithm.

####  III-C1 Complex and frequency-dependent $l_{p}$ approximation

Chapter II introduced the problem of designing $l_{p}$ complex FIR filters.
The complex $l_{p}$ IIR algorithm builds up on its FIR counterpart by
introducing a nested structure that internally solves for an $l_{2}$ complex
IIR problem. Figure 20 illustrates this procedure in more detail. This method
was first presented in [66].

Figure 21: Results for complex $l_{100}$ IIR design.

Figure 22: Maximum error for $l_{2}$ and $l_{100}$ complex IIR designs.

Figure 23: Error curve for $l_{100}$ complex IIR design.

Compared to its FIR counterpart, the IIR method only replaces the weighted
linear least squares problem for Soewito’s quasilinearization algorithm. While
this nesting approach might suggest an increase in computational expense, it
was found in practice that after the initial $l_{2}$ iteration, in general the
$l_{p}$ iterations only require from one to only a few internal weighted
$l_{2}$ quasilinearization iterations, thus maintaining the algorithm
efficiency. Figures 21 through 23 present results for a design example using a
length-5 IIR filter with $p=100$ and transition edge frequencies of 0.2 and
0.24 (in normalized frequency).

Figure 21 compares the $l_{2}$ and $l_{p}$ results and includes the desired
frequency samples. Note that no transition band was specified. Figure 22
illustrates the effect of increasing $p$. The largest error for the $l_{2}$
solution is located at the transition band edges. As $p$ increases the
algorithm weights the larger errors heavier; as a result the largest errors
tend to decrease. In this case the magnitude of the frequency response went
from 0.155 at the stopband edge (in the $l_{2}$ case) to 0.07 (for the $l_{p}$
design). Figure 23 shows the error function for the $l_{p}$ design,
illustrating the quasiequiripple behavior for large values of $p$.

Another fact worth noting from Figure 21 is the increase in the peak in the
right hand side of the passband edge (around $f=0.22$). The $l_{p}$ solution
increased the amplitude of this peak with respect to the corresponding $l_{2}$
solution. This is to be expected, since this peak occurs at frequencies not
included in the specifications, and since the $l_{p}$ algorithm will move
poles and zeros around in order to meet find the optimal $l_{p}$ solution
(based on the frequencies included for the filter derivation). The addition of
a specified transition band function (such as a spline) would allow for
control of this effect, depending on the user’s preferences.

The frequency-dependent FIR problem was first introduced in Section II-E.
Following the FIR approach, one can design IIR frequency-dependent filters by
merely replacing the linear weighted least squares step by a nonlinear
approach, such as the quasilinearization method presented in Section III-B13
(as in the complex $l_{p}$ IIR case). This problem illustrates the flexibility
in design for $l_{p}$ IRLS-based methods.

Figure 24: Block diagram for magnitude $l_{p}$ IIR method.

####  III-C2 Magnitude $l_{p}$ IIR design

The previous sections present algorithms that are based on complex
specifications; that is, the user must specify both desired magnitude and
phase responses. In some cases it might be better to specify a desired
magnitude response only, while allowing an algorihm to select the phase that
optimally minimizes the magnitude error. Note that if an algorithm is given a
phase in addition to a magnitude function, it must then make a compromise
between approximating both functions. The magnitude $l_{p}$ IIR approximation
problem overcomes this dilemma by posing the problem only in terms of a
desired magnitude function. The algorithm would then find the optimal phase
that provides the optimal magnitude approximation. A mathematical formulation
follows,

$\min_{\boldsymbol{a},\boldsymbol{b}}\;\;\left\|\;\left|D(\omega)\right|-\left|\frac{B(\omega;\boldsymbol{b})}{A(\omega;\boldsymbol{a})}\right|\;\right\|_{p}^{p}$
(114)

A critical idea behind the magnitude approach is to allow the algorithm to
find the optimum phase for a magnitude approximation. It is important to
recognize that the optimal magnitude filter indeed has a complex frequency
response. Atmadji Soewito [35] published in 1990 a theorem in the context of
$l_{2}$ IIR design that demonstrated that the phase corresponding to an
optimal magnitude approximation could be found iteratively by updating the
desired phase in a complex approximation scenario. In other words, given a
desired complex response $D_{0}$ one can solve a complex $l_{2}$ problem and
take the resulting phase to form a new desired response $D^{+}$ from the
original desired magnitude response with the new phase. That is,

$D_{i+1}=|D_{0}|e^{j\phi_{i}}$

where $D_{0}$ represents the original desired magnitude response and
$e^{j\phi_{i}}$ is the resulting phase from the previous iteration. This
approach was independently suggested [36] by Leland Jackson and Stephen Kay in
2008.

This work introduces an algorithm to solve the magnitude $l_{p}$ IIR problem
by combining the IRLS-based complex $l_{p}$ IIR algorithm from Section III-C1
with the phase updating ideas from Soewito, Jackson and Kay. The resulting
algorithm is robust, efficient and flexible, allowing for different orders in
the numerator and denominator as well as even or uneven sampling in frequency
space, plus the optional use of specified transition bands. A block diagram
for this method is presented in Figure 24.

The overall $l_{p}$ IIR magnitude procedure can be summarized as follows,

  1. 1.

Experimental analysis demonstrated that a reasonable initial solution for each
of the three main stages would allow for faster convergence. It was found that
the frequency domain Prony method by Burrus [28] (presented in Section III-B4)
offered a good initial guess. In Figure 24 this method is iterated to update
the specified phase. The outcome of this step would be an equation error
$l_{2}$ magnitude design.

  2. 2.

The equation error $l_{2}$ magnitude solution from the previous step
initializes a second stage where one uses quasilinearization to update the
desired phase. Quasilinearization solves the true solution error complex
approximation. Therefore by iterating on the phase one finds at convergence a
solution error $l_{2}$ magnitude design.

  3. 3.

The rest of the algorithm follows the same idea as in the previous step,
except that the least squared step becomes a weighted one (to account for the
necessary $l_{p}$ homotopy weighting). It is also crucial to include the
partial updating introduced in Section II-B2. By iterating on the weights one
would find a solution error $l_{p}$ magnitude design.

Figures 25 through 29 illustrate the effectiveness of this algorithm at each
of the three different stages for length-5 filters $\boldsymbol{a}$ and
$\boldsymbol{b}$, with transition edge frequencies of 0.2 and 0.24 (in
normalized frequency) and $p=30$. A linear transition band was specified.
Figures 25, 25 and 25 show the equation error $l_{2}$, solution error $l_{2}$
and solution error $l_{p}$. Figure 28 shows a comparison of the magnitude
error functions for the solution error $l_{2}$ and $l_{p}$ designs. Figure 29
shows the phase responses for the three designs.

Figure 25: Equation error $l_{2}$ magnitude design.

Figure 26: Solution error $l_{2}$ magnitude design.

Figure 27: Solution error $l_{p}$ magnitude design.

Figure 28: Comparison of $l_{2}$ and $l_{p}$ IIR magnitude designs

Figure 29: Phase responses for $l_{2}$ and $l_{p}$ IIR magnitude designs.

From Figures 28 and 29 one can see that the algorithm has changed the phase
response in a way that makes the maximum magnitude error (located in the
stopband edge frequency) to be reduced by approximately half its value.
Furthermore, Figure 28 demonstrates that one can reach quasiequiripple
behavior with relatively low values of $p$ (for the examples shown, $p$ was
set to 30).

##  IV Conclusions

Digital filters are essential building blocks for signal processing
applications. One of the main goals of this work is to illustrate the
versatility and relevance of $l_{p}$ norms in the design of digital filters.
While popular and well understood, $l_{2}$ and $l_{\infty}$ filters do tend to
accentuate specific benefits from their respective designs; filters designed
using $l_{p}$ norms as optimality criteria can offer a tradeoff between the
benefits of these two commonly used criteria. This work presented a number of
applications of $L_{p}$ norms in both FIR and IIR filter design, and their
corresponding design algorithms and software implementation.

The basic workhorse for the methods presented in this document is the
Iterative Reweighted Least Squares algorithm, a simple yet powerful method
that sets itself naturally adept for the design of $l_{p}$ filters. The notion
of converting a mathematically complex problem into a series of significantly
easier optimization problems is common in optimization. Nevertheless, the
existence results from Theorem 1 strongly motivate the use of IRLS methods to
design $l_{p}$ filters. Knowing that optimal weights exist that would turn the
solution of a weighted least squared problem into the solution of a least-$p$
problem must at the very least captivate the curiosity of the reader. The
challenge lies in finding a robust and efficient method to find such weights.
All the methods presented in this work work under this basic framework,
updating iteratively the weighting function of a least squares problem in
order to find the optimal $l_{p}$ filter for a given application. Therefore it
is possible to develop a suite of computer programs in a modular way, where
with few adjustments one can solve a variety of problems.

Throughout this document one can find examples of the versatility of the IRLS
approach. One can change the internal linear objective function from a complex
exponential kernel to a sinusoidal one to solve complex and linear phase FIR
filters respectively using the same algorithm. Further adaptations can be
incorporated with ease, such as the proposed adaptive solution to improve
robustness.

Another important design example permits to make $p$ into a function of
frequency to allow for different $p$-norms in different frequency bands. Such
design merely requires a few changes in the implementation of the algorithm,
yet allows for fancier, more elegant problems to be solved, such as the
Constrained Least Squares (CLS) problem. In the context of FIR filters, this
document presents the CLS problem from an $l_{p}$ prespective. While the work
by John Adams [16] set a milestone in digital filter design, this work
introduces a strong algorithm and a different perspective to the problem from
that by Adams and other authors. The IRLS $l_{p}$-based approach from this
work proves to be robust and flexible, allowing for even and uneven sampling.
Furthermore, while a user can use fixed transition bands, one would benefit
much from using a flexible transition band formulation, where the proposed
IRLS-based algorithm literally finds the optimal transition band definition
based on the constraint specifications. Such flexibility allows for tight
constrains that would otherwise cause other algorithms to fail meeting the
constraint specifications, or simply not converging at all. Section II-F
introduced two problem formulations as well as results that illustrate the
method’s effectiveness at solving the CLS problem.

While previous work exists in the area of FIR design (or in linear $l_{p}$
approximation for that matter), the problem of designing $l_{p}$ IIR filters
has been far less explored. A natural reason for this is the fact that $l_{2}$
IIR design is in itself an open research area (and a rather complicated
problem as well). Traditional linear optimization approaches cannot be
directly used for either of these problems, and nonlinear optimization tools
often prove either slow or do not converge.

This work presents the $l_{p}$ IIR design problem as a natural extension of
the FIR counterpart, where in a modular fashion the linear weigthed $l_{2}$
section of the algorithms is replaced by a nonlinear weighted $l_{2}$ version.
This problem formulation allows for the IIR implementation of virtually all
the IRLS FIR methods presented in Chapter II. Dealing with the weighted
nonlinear $l_{2}$ problem is a different story.

The problem of rational $l_{2}$ approximation has been studied for some time.
However the sources of ideas and results related to this problem are scattered
across several areas of study. One of the contributions of this work is an
organized summary of efforts in rational $l_{2}$ optimization, particularly
related to the design of IIR digital filters. The work in Section III-B also
lays down a framework for the IIR methods proposed in this work.

As mentioned in Section III-C, some complications arise when designing IIR
$l_{p}$ filters. Aside from the intrinsic $l_{2}$ problem, it is necessary to
properly combine a number of ideas that allowed for robust and efficient
$l_{p}$ FIR methods. A design algorithm for complex $l_{p}$ IIF filters were
presented in Section III-C1; this algorithm combined Soewito’s
quasilinearization with ideas such as $l_{p}$ homotopy, partial updating and
the adaptive modification. In practice, the combination of these ideas showed
to be practical and the resulting algorithm remained robust. It was also found
that after a few $p$-steps, the internal $l_{2}$ algorithm required from one
to merely a few iterations on average, thus maintaining the algorithm
efficient.

One of the main contributions of this work is the introduction of an IRLS-
based method to solve $l_{p}$ IIR design problems. By properly combining the
principle of magnitude approximation via phase updating (from Soewito, Jackson
and Kay) with the complex IIR algorithm one can find optimal magntiude $l_{p}$
designs. This work also introduced a sequence of steps that improve the
efficiency and robustness of this algorithm, by dividing the design process
into three stages and by using suitable initial guesses for each stage.

Some of the examples in this document were designed using Matlab programs. It
is worth to notice the common elements between these programs, alluding to the
modularity of the implementations. An added benefit to this setup is that
further advances in any of the topics covered in this work can easily be
ported to most if not all of the algorithms.

Digital filter design is and will remain an important topic in digital signal
processing. It is the hope of the author to have motivated in the reader some
curiosity for the use of $l_{p}$ norms as design criteria for applications in
FIR and IIR filter design. This work is by no means comprehensive, and is
meant to inspire the consideration of the flexibility of IRLS algorithms for
new $l_{p}$ related problems.

## References

  * [1] R. E. Ziemer, W. H. Tranter, and D. R. Fannin, _Signals and Systems: Continuous and Discrete_ , 4th ed. Prentice Hall, 1998.
  * [2] J. L. Walsh and T. S. Motzkin, “Polynomials of Best Approximation on an Interval,” _Proceeedings of the National Academy of Sciences, USA_ , vol. 45, pp. 1523–1528, October 1959.
  * [3] T. S. Motzkin and J. L. Walsh, “Polynomials of Best Approximation on a Real Finite Point Set I,” _Trans. American Mathematical Society_ , vol. 91, no. 2, pp. 231–245, May 1959.
  * [4] C. L. Lawson, “Contributions to the theory of linear least maximum approximations,” Ph.D. dissertation, UCLA, 1961.
  * [5] J. R. Rice and K. H. Usow, “The Lawson Algorithm and Extensions,” _Mathematics of Computation_ , vol. 22, pp. 118–127, 1968.
  * [6] J. R. Rice, _The Approximation of Functions_. Addison-Wesley, 1964, vol. 1.
  * [7] C. S. Burrus, J. A. Barreto, and I. W. Selesnick, “Iterative Reweighted Least-Squares Design of FIR Filters,” _IEEE Transactions on Signal Processing_ , vol. 42, no. 11, pp. 2926–2936, November 1994.
  * [8] L. A. Karlovitz, “Construction of Nearest Points in the $L^{p}$, $p$ even and $L^{\infty}$ norms, I.” _Journal of Approximation Theory_ , vol. 3, pp. 123–127, 1970.
  * [9] S. W. Kahng, “Best $L_{p}$ Approximations,” _Mathematics of Computation_ , vol. 26, no. 118, pp. 505–508, April 1972.
  * [10] R. Fletcher, J. A. Grant, and M. D. Hebden, “The Calculation of Linear Best $L_{p}$ Approximations,” _The Computer Journal_ , vol. 14, no. 118, pp. 276–279, Apr 1972.
  * [11] M. Aoki, _Introduction to Optimization Techniques_. The Macmillan Company, 1971.
  * [12] J. A. Barreto, “$L_{p}$ Approximation by the Iterative Reweighted Least Squares Method and the Design of Digital FIR Filters in One Dimension,” Master’s thesis, Rice University, 1992.
  * [13] C. S. Burrus and J. A. Barreto, “Least $p$-power Error Design of FIR Filters,” in _Proc. IEEE Int. Symp. Circuits, Syst. ISCAS-92_ , vol. 2, San Diego, CA, May 1992, pp. 545–548.
  * [14] J. Nocedal and S. J. Wright, _Numerical Optimization_ , ser. Springer series in operations research. New York, NY: Springer-Verlag, 1999.
  * [15] J. W. Adams and J. L. Sullivan, “Peak-Constrained Least-Squares Optimization,” _IEEE Trans. on Signal Processing_ , vol. 46, no. 2, pp. 306–321, Febr. 1998.
  * [16] J. W. Adams, “FIR Digital Filters with Least-Squares Stopbands Subject to Peak-Gain Constraints,” _IEEE Trans. on Circuits and Systems_ , vol. 39, no. 4, pp. 376–388, April 1991.
  * [17] I. W. Selesnick, M. Lang, and C. S. Burrus, “Constrained Least Square Design of FIR Filters without Specified Transition Bands,” _IEEE Transactions on Signal Processing_ , vol. 44, no. 8, pp. 1879–1892, August 1996.
  * [18] M. Lang, I. W. Selesnick, and C. S. Burrus, “Constrained Least Square Design of 2-D FIR Filters,” _IEEE Transactions on Signal Processing_ , vol. 44, no. 5, pp. 1234–1241, May 1996.
  * [19] I. W. Selesnick, M. Lang, and C. S. Burrus, “A Modified Algorithm for Constrained Least Square Design of Multiband FIR Filters Without Specified Transition Bands,” _IEEE Transactions on Signal Processing_ , vol. 46, no. 2, pp. 497–501, Feb. 1998.
  * [20] R. Fletcher and M. J. D. Powell, “A Rapidly Convergent Descent Method for Minimization,” _Computer Journal_ , vol. 6, no. 2, pp. 163–168, 1963.
  * [21] J. E. Dennis and R. B. Schnabel, _Numerical Methods for Unconstrained Optimization and Nonlinear Equations_. Philadelphia, PA: SIAM, 1996.
  * [22] L. S. e. a. T. P. Krauss, _Signal Processing Toolbox User’s Guide_. The MathWorks, 1994, ch. 2, pp. 143–145.
  * [23] J. T. Spanos and D. L. Mingori, “Newton Algorithm for Fitting Transfer Functions to Frequency Measurements,” _Journal of Guidance, Control and Dynamics_ , vol. 16, no. 1, pp. 34–39, January 1993.
  * [24] D. C. Sorensen, “Newton’s Method with a Dodel Trust Region Modification,” _SIAM Journal of Numerical Analysis_ , vol. 16, pp. 409–426, 1982.
  * [25] R. Kumaresan and C. S. Burrus, “Fitting a Pole-Zero Filter Model to Arbitrary Frequency Response Samples,” _Proc. ASILOMAR_ , pp. 1649–1652, 1991.
  * [26] F. B. Hildebrand, _Introduction to Numerical Analysis_. McGraw-Hill, 1974.
  * [27] E. C. Levy, “Complex-Curve Fitting,” _IRE Transactions on Automatic Control_ , vol. AC-4, no. 1, pp. 37–43, May 1959.
  * [28] T. W. Parks and C. S. Burrus, _Digital Filter Design_. John Wiley and Sons, 1987.
  * [29] B. G. C. F. M. R. de Prony, “Essai Expérimental et Analytique: Sur les lois de la Dilatabilité des fluides élastiques et sur celles de la Force expansive de la vapeur de l’eau et de la vapeur de l’alkool, à différentes températures.” _Journal de l’École Polytechnique (Paris)_ , vol. 1, no. 2, pp. 24–76, 1795.
  * [30] H. E. Padé, “Sur la Représentation Approchée d’une Fonction par des Fractions Rationnelles,” _Annales Scientifiques de LÉcole Normale Supérieure (Paris)_ , vol. 9, no. 3, pp. 1–98, 1892.
  * [31] C. K. Sanathanan and J. Koerner, “Transfer Function Synthesis as a Ratio of Two Complex Polynomials,” _IEEE Transactions on Automatic Control_ , vol. AC-8, pp. 56–58, January 1963.
  * [32] A. C. M. A. Sid-Ahmed and G. A. Jullien, “Computational Techniques for Least-Square Design of Recursive Digital Filters,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. ASSP-26, no. 5, pp. 477–480, October 1978.
  * [33] H. W. S. L. E.McBride and K. Steiglitz, “Time-Domain Approximation by Iterative Methods,” _IEEE Transactions on Circuit Theory_ , vol. CT-13, no. 4, pp. 381–87, December 1966.
  * [34] K. Steiglitz and L. E. McBride, “A Technique for the Identification of Linear Systems,” _IEEE Transactions on Automatic Control_ , vol. AC-10, pp. 461–64, October 1965.
  * [35] A. W. Soewito, “Least square digital filter design in the frequency domain,” Ph.D. dissertation, Rice University, December 1990.
  * [36] L. B. Jackson, “Frequency-Domain Steiglitz-McBride Method for Least-Squares IIR Filter Design, ARMA Modeling, and Periodogram Smoothing,” _IEEE Signal Processing Letters_ , vol. 15, pp. 49–52, 2008.
  * [37] A. Antoniou, _Digital Filters: Analysis, Design, and Applications_ , 2nd ed. McGraw-Hill, 1993.
  * [38] E. Cunningham, _Digital Filtering: An Introduction_. Houghton-Mifflin, 1992.
  * [39] E. W. Cheney, _Introduction to Approximation Theory_ , ser. Intl. Series in Pure and Applied Mathematics. McGraw-Hill, 1966.
  * [40] V. Chvatal, _Linear Programming_. Freeman and Co., 1980.
  * [41] G. Strang, _Introduction to Applied Mathematics_. Wellesley-Cambridge Press, 1986.
  * [42] S. A. Ruzinsky, “$L_{1}$ and $L_{\infty}$ Minimization via a Variant of Karmarkar’s Algorithm,” _IEEE Trans. on Acoustics, Speech and Signal Processing_ , vol. 37, no. 2, pp. 245–253, Febr. 1989.
  * [43] A. V. Oppenheim and R. W. Schafer, _Discrete-Time Signal Processing_. Englewood Cliffs, N.J.: Prentice-Hall, 1989.
  * [44] C. S. Burrus, J. H. McClellan _et al._ , _Computer-based Exercises for Signal Processing_. Prentice Hall, 1994\.
  * [45] R. A. Vargas and C. S. Burrus, “Adaptive Iterative Reweighted Least Squares Design of $L_{p}$ FIR Filters,” _Proc. ICASSP_ , vol. 3, Signal Processing, Theory and Methods, pp. 1129–32, 1999.
  * [46] K. S. Shanmugan and A. M. Breipohl, _Random Signals: Detection, Estimation and Data Analysis_. John Wiley & Sons, 1988.
  * [47] J. G. Proakis and D. G. Manolakis, _Digital Signal Processing_. Macmillan Publishing Co., 1988.
  * [48] S.-P. Wu, S. Boyd, and L. Vandenberghe, “Fir filter design via spectral factorization and convex optimization,” to Appear in Applied Computational Control, Signal and Communications, Biswa Datt editor, Birkhauser, 1997.
  * [49] K. Steiglitz, “Computer-Aided Design of Recursive Digital Filters,” _IEEE Transactions on Audio and Electroacoustics_ , vol. AU-18, no. 2, pp. 123–129, June 1970.
  * [50] A. Deczky, “Synthesis of Recursive Digital Filters Using the Minimum $p$-error Criterior,” _IEEE Transactions on Audio and Electroacoustics_ , vol. AU-20, no. 4, pp. 257–263, October 1972.
  * [51] R. Kumaresan, “Identification of Rational Transfer Functions from Frequency Response Samples,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 26, no. 6, pp. 925–934, Novermber 1990.
  * [52] R. E. Mickens, _Difference equations_. Van Nostrand Reinhold, 1987.
  * [53] S. N. Elaydi, _An Introduction to Difference Equations_ , ser. Undergraduate texts in mathematics. New York: Springer, 1996.
  * [54] D. F. T. Jr., “On Fluids, Networks, and Engineering Education,” in _Aspects of Network_ , R. E. Kalman and N. DeClaris, Eds. Hol, Rinehart and Winston, Inc, 1971, pp. 591–612.
  * [55] M. L. V. Blaricum, “A Review of Prony’s Method techniques for Parameter Estimation,” in _Air Force Statistical Estimation Workshop_ , May 1978, pp. 125–135.
  * [56] L. Weiss and R. N. McDonough, “Prony’s method, Z-transforms , and Padé approximation,” _SIAM Review_ , vol. 5, no. 2, pp. 145–149, April 1963\.
  * [57] I. Barrondale and D. D. Olesky, “Exponential Approximation using Prony’s Method,” in _The Numerical Solution of Nonlinear Problems_ , C. T. H. Baker and C. Phillips, Eds. New York: Oxford University Press, 1981, ch. 20, pp. 258–269.
  * [58] S. L. Marple, Jr, _Digital Spectral Analysis with Applications_ , ser. Signal Processing Series. Englewood Cliffs, NJ: Prentice-Hall, 1987.
  * [59] C. F. Gerald and P. O. Wheatley, _Applied Numerical Analysis_. Reading, MA: Addison-Wesley, 1984.
  * [60] H. Cabannes, Ed., _Pade Approximates Method and its Applications to Mechanics_ , ser. Lecture notes in physics. Springer-Verlag, 1976, no. 47.
  * [61] G. H. Golub and C. F. V. Loan, _Matrix Computations_. Johns Hopkins University Press, 1996.
  * [62] H. Fan and M. Doroslovački, “On ”Global Convergence” of Steiglitz-McBride Adaptive Algorithm,” _IEEE Transactions on Circuits and Systems II_ , vol. 40, no. 2, pp. 73–87, February 1993.
  * [63] P. Stoica and T. Söderström, “The Steiglitz-McBride Identification Algorithm Revisited-Convergence Analysis and Accuracy Aspects,” _IEEE Transactions on Automatic Control_ , vol. AC-26, no. 3, pp. 712–17, June 1981\.
  * [64] J. H. McClellan and D. Lee, “Exact Equivalence of the Steiglitz-McBride Iteration and IQML,” _IEEE Transactions on Signal Processing_ , vol. 39, no. 2, pp. 509–12, February 1991.
  * [65] R. E. Bellman and R. E. Kalaba, _Quasilinearization and Nonlinear Boundary-Value Problems_. New York: American Elsevier, 1965.
  * [66] R. A. Vargas and C. S. Burrus, “On the Design of $L_{p}$ IIR Filters with Arbitrary Frequency Bands,” in _Proc. ICASSP_ , vol. 6, 2001, pp. 3829–3832.