The preceding design methods have been based on designing an
analog prototype filter and then converting it to a digital
filter. This approach is appropriate for the class of
approximations where analytical solutions are possible, but not
for many others. In the remaining part of this chapter, methods
will be developed that directly design the desired digital
filter. Most approaches are extensions of methods used for FIR
filters, but they are more complicated for the IIR case where
rational approximation is being performed rather than polynomial
approximation.
In this section a frequency-sampling design method is
developed such that the frequency response of the IIR filter will pass
through the given samples of a desired response. Since an IIR filter
cannot have linear phase, the sampled response must contain both magnitude
and phase. The extension of the frequency- sampling method to a LS-error
approximation is not as simple as for the FIR filter. The method presented
in this section uses a criterion based on the equation error rather than
the more common error between the actual and desired frequency
response[1]. Nevertheless, it is a useful noniterative design
method. Finally, a general discussion of iterative design methods for
LS-frequency response error is given.
The method for calculating samples of the frequency response
of an IIR filter presented in the section on Properties of IIR Filters
can be reversed to design
a filter much the same way it was for the FIR filter using frequency sampling.
The z-transform transfer function for an IIR filter is given by
H
(
z
)
=
B
(
z
)
A
(
z
)
=
b
0
+
b
1
z
-
1
+
⋯
+
b
M
z
-
M
1
+
a
1
z
-
1
+
⋯
+
a
N
z
-
N
.
H
(
z
)
=
B
(
z
)
A
(
z
)
=
b
0
+
b
1
z
-
1
+
⋯
+
b
M
z
-
M
1
+
a
1
z
-
1
+
⋯
+
a
N
z
-
N
.
(1)
The frequency response of the filter is given by setting z=e-jωz=e-jω. Using the notation
H
(
ω
)
=
H
(
z
)
|
z
=
e
-
j
ω
.
H
(
ω
)
=
H
(
z
)
|
z
=
e
-
j
ω
.
(2)
Equally-spaced samples of the frequency response are chosen so
that the number of samples is equal to the number of unknown
coefficients in (Equation 1). These L+1L+1 = M+N+1M+N+1 samples of this
frequency response are given by
H
k
=
H
(
ω
k
)
=
H
(
2
π
k
L
+
1
)
H
k
=
H
(
ω
k
)
=
H
(
2
π
k
L
+
1
)
(3)
and can be calculated from the length-(L+1)(L+1)) DFTs of the numerator
and denominator.
H
k
=
DFT
{
b
n
}
DFT
{
a
n
}
=
B
k
A
k
H
k
=
DFT
{
b
n
}
DFT
{
a
n
}
=
B
k
A
k
(4)
where the indicated division is term-by-term division for each value
of kk. Multiplication of both sides of (Equation 4) by AkAk gives
B
k
=
H
k
A
k
B
k
=
H
k
A
k
(5)
If the length-(L+1)(L+1) inverse DFT of HkHk is denoted by the length-
(L+1)(L+1) sequence hnhn, equation (Equation 5) becomes cyclic convolution
which can be expressed in matrix form by
b
0
b
1
⋮
b
M
0
⋮
0
=
h
0
h
L
h
L
-
1
⋯
h
1
h
1
h
0
h
L
h
2
h
1
h
0
⋮
⋮
h
L
⋯
h
0
1
a
1
⋮
a
N
0
⋮
0
b
0
b
1
⋮
b
M
0
⋮
0
=
h
0
h
L
h
L
-
1
⋯
h
1
h
1
h
0
h
L
h
2
h
1
h
0
⋮
⋮
h
L
⋯
h
0
1
a
1
⋮
a
N
0
⋮
0
(6)
Note that the hnhn in (Equation 6) are not the impulse response values
of the filter as used in the FIR case. A more compact matrix notation of is
b
0
=
H
a
0
b
0
=
H
a
0
(7)
where HH is (L+1)(L+1) by (L+1)(L+1), bb is length-(M+1)(M+1),
and aa is length-(N+1)(N+1). Because the lower L-NL-N terms of
the right-hand vector of (Equation 6) are zero, the HH matrix can be
reduced by deleting the right-most L-NL-N columns to give H0H0
which causes (Equation 7) to become
b
0
=
H
0
a
b
0
=
H
0
a
(8)
Because the first element of aa is unity, it is partitioned
to remove the unity term and the remaining length-NN vector is
denoted a*a*. The simultaneous equations represented by (Equation 8)
are uncoupled by further partitioning of the HH matrix as
shown in
b
0
=
H
1
h
1
H
2
1
a
*
b
0
=
H
1
h
1
H
2
1
a
*
(9)
where H1H1 is (M+1)(M+1) by (N+1)(N+1), h1h1 is
length-(L-M)(L-M), and H2H2 is (L-M)(L-M) by NN. The lower (L-M)(L-M)
equations are written
0
=
h
1
+
H
2
a
*
0
=
h
1
+
H
2
a
*
(10)
or
h
1
=
-
H
2
a
*
h
1
=
-
H
2
a
*
(11)
which must be solved for a*a*. The upper M+1M+1 equations of
(10) are written
b
=
H
1
a
b
=
H
1
a
(12)
which allows the calculation of bb.
If L=N+ML=N+M, H2H2 is square. If H2H2 is
nonsingular, (Equation 11) can be solved exactly for the denominator
coefficients in a*a*, which are augmented by the unity term
to give aa. From (Equation 12), the numerator coefficients in bb are found. If H2H2 is singular [16], (Reference) and there
are multiple solutions, a lower order problem can be posed. If
there are no solutions, the approximation methods must be
used.
Note that any order numerator and denominator can be prescribed. If
the filter is in fact an FIR filter, aa is unity and a*a* does not exist. Under these conditions, (Equation 12) states that bn=hnbn=hn, which is one of the cases of FIR frequency sampling covered
[19]. Also note that there is no control
over the stability of the filter designed by this method.
ummary
In this section, an interpolation design method was developed and
analyzed. Use of the DFT converted the frequency- domain
specifications to the time domain. A matrix partitioning allowed
uncoupling the solution for the numerator from the solution of the
denominator coefficients. The use of the DFT prevents the
possibility of unequally spaced frequency samples as was possible
for FIR filter design. The solution of simultaneous equations would
allow unequal spacing which is not as troublesome as with the FIR
filter because IIR filters are usually of lower order.
The frequency-sampling design of IIR filters is somewhat more
complicated than for FIR filters because of the requirement that
H2H2 be nonsingular. As for the FIR filter, the samples of
the desired frequency response must satisfy the conditions to insure
that hnhn are real. The power of this method is its ability to
interpolate arbitrary magnitude and phase specification. In contrast
to most direct IIR design methods, this method does not require any
iterative optimization with the accompanying convergence problems.
As with the FIR version, because this design approach is an
interpolation method rather than an approximation method, the
results may be poor between the interpolation points. This usually
happens when the desired frequency-response samples are not
consistent with what an IIR filter can achieve. One solution to
this problem is the same as for the FIR case
[19], the use of more frequency samples than the number of
filter coefficients and the definition of an approximation error
function that can be minimized. There is no simple restriction that
will guarantee stable filters. If the frequency-response samples
are consistent with an unstable filter, that is what will be
designed.
In order to obtain better practical filter designs, the
interpolation scheme of the previous section is extended to give an
approximation design method [19]. It should be noted at the
outset that the method developed in this section minimizes an
equation-error measure and not the usual frequency-response error
measure.
The number of frequency samples specified, L+1L+1, will be made
larger than the number of filter coefficients, M+N+1M+N+1. This means
that H2H2 is rectangular and, therefore, (Equation 11) cannot in
general be satisfied. To formulate an approximation problem, a
length-(L+1)(L+1) error vector εε is introduced in (Equation 8)
and (Equation 9) to give
b
0
=
H
0
a
+
ε
b
0
=
H
0
a
+
ε
(13)
Equation (Equation 11) becomes
h
1
-
ε
=
-
H
2
a
*
h
1
-
ε
=
-
H
2
a
*
(14)
where now H2H2 is rectangular with (L-M)>N(L-M)>N. Using
the same methods as used to derive (Equation 11), the error εε is minimized in a least-squared error sense by the
solution of the normal equations [16]
H
2
T
h
1
=
-
H
2
T
H
2
a
*
H
2
T
h
1
=
-
H
2
T
H
2
a
*
(15)
If the equations are not singular, the solution is
a
*
=
-
[
H
2
T
H
2
]
-
1
H
2
T
h
1
.
a
*
=
-
[
H
2
T
H
2
]
-
1
H
2
T
h
1
.
(16)
If the normal equations are singular, the pseudo-inverse
[16], (Reference) can be used to obtain a minimum norm or reduced
order solution.
The numerator coefficients are found by the same techniques as
before in (Equation 12)
b
=
H
1
a
b
=
H
1
a
(17)
which results in the upper M+1M+1 terms in εε
being zero and the total squared equation error being minimum.
As is true for LS-error design of FIR filters, (Equation 15) is often
numerically ill-conditioned and (Equation 16) should not be used to solve for
a*a*. Special algorithms such as those used by Matlab and
LINPACK [18], [7] should be employed.
The error εε defined in (Equation 13) can better be
understood by considering the frequency-domain formulation. Taking
the DFT of (Equation 13) gives
B
k
=
H
k
A
k
+
ε
B
k
=
H
k
A
k
+
ε
(18)
where εε is the error in trying to satisfy (Equation 8) when the
equations are over-specified. This can be reformulated in terms of
EE, the difference between the frequency response samples of
the designed filter and the desired response samples, by dividing
(Equation 8) by AkAk to give
E
k
=
B
k
A
k
-
H
k
=
ε
k
A
k
E
k
=
B
k
A
k
-
H
k
=
ε
k
A
k
(19)
EE is the error in the solution of the approximation
problem, and εε is the error in the equations defining
the problem. The usual statement of a frequency-domain approximation
problem is in terms of minimizing some measure of EE, but
that results in solving nonlinear equations. The design procedure
developed in this section minimizes the squared error εε,
thus only requiring the solution of linear equations. There is an
important relation between these problems. Equation (Equation 19) shows that
minimizing εε is the same as minimizing EE
weighted by AA. However, AA is unknown until after the problem is
solved.
Although this is posed as a frequency-domain design method, the
method of solution for both the interpolation problem and the LS
equation-error problem is the same as the time-domain Prony's
method, discussed in Section 7.5 of reference [19].
Numerous modifications and extensions can be made to this method. If the
desired frequency response is close to what can be achieved by an IIR
filter, this method will give a design approximately the same as that of a
true least-squared solution-error method. It can be shown that
ε=0↔E=0ε=0↔E=0. In some cases, improved
results can be obtained by estimating AkAk and using that as a weight on
εε to approximate minimizing EE. There are iterative
methods based on solving (Equation 16) and (Equation 17) to obtain values
for AkAk. These values are used as weights on εε to solve for
a new set of AkAk used as a new set of weights to solve again for AkAk[19][26]. We found this approach to converge slowly, but a recent paper
using the log-magnitude [14] was more successful. Other
approaches are given in [24], [23], [11]. The solution of
(Equation 16) and (Equation 17) is sometimes used to obtain starting
values for other iterative optimization algorithms that need good starting
values for convergence.
An interesting iterative design algorithm that can design to approximate
complex or magnitude frequency responses has be recently proposed by
Jackson [11]. A different approach to the same problem was posed
by Soewito [26], [28].
To illustrate this design method a sixth-order lowpass filter was
designed with 41 frequency samples to approximate. The magnitude of
those less than 0.2 Hz is one and of those greater than 0.2 is zero.
The phase was experimentally adjusted to result in a good magnitude
response. The design was performed with Program 9 in the appendix
of [19] and the frequency response is shown in Figure 7-33 of
[19]. Matlab programs have recently been written which are
smaller and easier to understand than those in FORTRAN.
ummary
In this section an LS-error approximation method was posed to
design IIR filters. By using an equation-error rather than a
solution-error criterion, a problem resulted that required only the
solution of simultaneous linear equations.
Like the FIR filter version, the IIR frequency sampling design
method and the LS equation-error extension can be used for complex
approximation and, therefore, can design with both magnitude and
phase specifications.
If the desired frequency-response samples are close to what an IIR
filter of the specified order can achieve, this method will produce
a filter very close to what a true least-squared error method would.
However, when the specifications are not consistent with what can be
achieved and the approximating error is large, the results can be
very poor and in some cases, unstable. It is particularly difficult
to set realistic phase response specifications. With this method, it
is even more important to have a design environment that will allow
easy trial-and-error procedure.
Newly published works which will be discussed here are
[14], [13], [22], [10], [20], [17], [12]. Other references can
be found in [19], [14], [6], [26], [28], [5]. The Matlab command invfreqz()
which is an inverse to the freqz()
command gives a
similar or, perhaps, the same result as the method described in this
note but uses a different formulation [15], [25].
Practical problems occur in the design of a filter to
separate signals according to their energy. Because the energy
content of a signal is the integral or sum of the square of the
signal, a mean-squared-error measure is natural. Unfortunately,
for the IIR filter design problem, the optimization procedure is
nonlinear. This was pointed out in the last section where the
equation error was used in order to have a linear problem.
Because of the nonlinear nature of the least-squared-error
minimization, the method of solution becomes dependent on the
desired frequency response, and therefore, there is no single
method for design. The mean-squared error for magnitude
approximation is defined as
q
(
x
)
=
∑
i
=
0
L
|
H
(
ω
i
)
|
-
|
H
d
(
ω
i
)
|
2
q
(
x
)
=
∑
i
=
0
L
|
H
(
ω
i
)
|
-
|
H
d
(
ω
i
)
|
2
(20)
where xx is a vector of filter parameters chosen to minimize q,
and the error is sampled at L+1L+1 frequencies ωiωi. Steiglitz
[19] chose the parameter vector x to be the coefficients of a
cascade structure in order to best fit an iterative optimization
scheme. He applied a standard optimization algorithm, the
Fletcher-Powell method, to the minimization of (Equation 19). Other
methods which are more directly related to a squared-error measure
can also be used.
Practical difficulties exist in solving this approximation
problem. In some cases, local minima are found rather than the
global minimum. In other cases, convergence of the minimization
algorithm is slow or does not occur at all. Numerical problems can
result from ill-conditioned equations, and there is no guarantee
that the designed filter will be stable.
An important factor is the choice of a desired frequency-
response function Hd(ω)Hd(ω) that does not result in the optimum
approximation having a large error. This often means not having an
abrupt discontinuity between the passband and stopband.
Another factor is the starting of the iterative optimization
algorithm with a set of coefficients in xx that is close to the
optimum. This can be accomplished by using the frequency sampling
method to give a design that can be
used to start a least-squares algorithm. Because the error defined
in (Equation 19) is in terms of magnitudes, an unstable design can be
converted to a stable one by moving the unstable pole at a radius of
rr in the zz-plane to a radius of 1/r1/r. This does not change the
magnitude frequency response and does stabilize the effect of that
pole [19].
A generalization of the idea of a squared-error measure is
defined by raising the error to the pp power where pp is a
positive integer. This error is defined by
q
(
x
)
=
∑
i
=
0
L
|
H
(
ω
)
-
H
d
(
ω
)
|
p
q
(
x
)
=
∑
i
=
0
L
|
H
(
ω
)
-
H
d
(
ω
)
|
p
(21)
Deczky [8] developed this approach and used the Fletcher-Powell
method to minimize (Equation 21). He also applied this method to the
approximation of a desired group-delay function. An important
characteristic of this formulation is that the solution approaches the
Chebyshev or mini-max solution as p becomes large. Initial work shows the
method of iteratively reweighted least squared error (IRLS) as was applied
to the FIR filter design in (Reference) can also be used for LpLp and
constrained least squared error optimal design of IIR filters
[29].
The error measure that often best meets filter design
specifications is the maximum error in the frequency response
that occurs over a band. The filter design problem becomes the
problem of minimizing the maximum error (the min-max problem).
Among several approaches to this error minimization, one is
by Deczky which minimizes the p-power error of (Equation 21) for large
p. Generally, p=10p=10 or greater approximates a Chebyshev result
[19]. Another is by Dolan and Kaiser which uses a
penalty-function approach.
Linear programming can be applied to this error measure
by linearizing the equations in much the same way as
in (Equation 15) [21]. In contrast to the FIR case, this can be a
practical design method because the order of practical IIR
filters is generally much lower than for FIR filters. A scheme
called differential correction has also proven to be effective.
Although the rational approximation problem is nonlinear, an
application of the Remes exchange algorithm can be implemented
[19]. Since the zeros of the numerator of the transfer
function mainly control the stopband characteristics of a filter,
and the zeros of the denominator mainly control the passband, the
effects of the two are somewhat uncoupled. An application of the
Remes exchange algorithm, alternating between the numerator and
denominator, gives an effective method for designing IIR filters
with a Chebyshev error criterion. If the order of the
numerator and denominator are the same and the desired filter is
an ideal lowpass filter, the Remes exchange should give the same
result as the elliptic function filter. However,
this approach allows any order numerator or denominator to be set
and any shape passband to be approximated. There are cases where
a lower-order denominator than numerator results in a filter with
fewer required muliplications than an elliptic-function filter
[19].
The problem of designing an IIR digital filter with a
prescribed time-domain response is addressed in this section.
Most formulations of time-domain design of IIR filters result in
nonlinear equations for the same reasons as for frequency-domain
design. Prony, in 1790, derived a special formulation for the
analysis of elastic properties of gases, which resulted in linear
equations. A more general form of Prony's method can be applied
to the IIR filter design by use of a matrix description [2], [19].
The transfer function of an IIR filter is given by
H
(
z
)
=
B
(
z
)
A
(
z
)
=
b
0
+
b
1
z
-
1
+
⋯
+
b
M
z
-
M
1
+
a
1
z
-
1
+
⋯
+
a
N
z
-
N
=
h
0
+
h
1
z
-
1
+
h
2
z
-
2
+
⋯
H
(
z
)
=
B
(
z
)
A
(
z
)
=
b
0
+
b
1
z
-
1
+
⋯
+
b
M
z
-
M
1
+
a
1
z
-
1
+
⋯
+
a
N
z
-
N
=
h
0
+
h
1
z
-
1
+
h
2
z
-
2
+
⋯
(22)
and the impulse response h(n)h(n) is related to H(z)H(z) by the zz
transform.
H
(
z
)
=
∑
n
=
0
∞
h
(
n
)
z
-
n
H
(
z
)
=
∑
n
=
0
∞
h
(
n
)
z
-
n
(23)
Equation (Equation 23) can be written
B
(
z
)
=
H
(
z
)
A
(
z
)
B
(
z
)
=
H
(
z
)
A
(
z
)
(24)
which is the z-transform version of convolution. This convolution
can be written as a matrix multiplication. Using the first K+1
terms of the impulse response, this is written
b
0
b
1
⋮
b
M
0
⋮
0
=
h
0
0
0
⋯
0
h
1
h
0
0
h
2
h
1
h
0
⋮
⋮
h
L
⋯
h
0
1
a
1
⋮
a
N
0
⋮
0
b
0
b
1
⋮
b
M
0
⋮
0
=
h
0
0
0
⋯
0
h
1
h
0
0
h
2
h
1
h
0
⋮
⋮
h
L
⋯
h
0
1
a
1
⋮
a
N
0
⋮
0
(25)
In order to uncouple the calculations of the anan and the bnbn,
the matrices are partitioned to give
b
0
=
H
1
h
1
H
2
1
a
*
b
0
=
H
1
h
1
H
2
1
a
*
(26)
where bb is the vector of the M+1M+1 numerator coefficients of (Equation 22),
a*a* is the vector of the NN denominator coefficients (ao=1ao=1), h1h1 is the
vector of the last (K-M)(K-M) terms of the impulse response, H1H1 is the
M+1M+1 by N+1N+1 partition of (Equation 25), and H2H2 is the (K-M)(K-M) by NN remaining
part. The lower K-MK-M equations are written
0
=
h
1
+
H
2
a
*
0
=
h
1
+
H
2
a
*
(27)
or
h
1
=
-
H
2
a
*
h
1
=
-
H
2
a
*
(28)
which must be solved for a*a*. The upper M+1M+1 equations of
(Equation 26) are written
b
=
H
1
a
b
=
H
1
a
(29)
which allows the calculation of bb.
If L=N+ML=N+M, then H2H2 is square. If H2H2 is nonsingular, (Equation 28)
can be solved exactly for the denominator coefficients in a*a*, which are augmented by the unity term to give aa.
From (Equation 29), the numerator coefficients in bb are found.
If H2H2 is singular [16], (Reference) and there are multiple
solutions, a lower order problem can be posed. If there are no
solutions, the methods of the next section must be used.
Note that any order numerator and denominator can be prescribed. If
the filter is in fact an FIR filter, aa is unity and a*a* does not exist. Under these conditions, (Equation 29) states
that bn=hnbn=hn, which is one of the cases of FIR frequency
sampling covered in Section 3.1 of [19]. Also note that there
is no control over the stability of the filter designed by this
method.
Although Prony's method, applied to the time-domain design
problem here, is similar to the solution of the frequency-sampling
IIR design problem, there are important differences. The
inverse DFT is used to obtain the matrix in the frequency domain problem, which is
cyclic convolution. Equation (Equation 25) is noncyclic convolution and
the K+1K+1 terms of h(n)h(n), used to form HH, result from a truncation
of the infinitely long sequence.
In order to obtain better practical filter designs, the interpolation
scheme of the previous section is extended to give an approximation design
method [19]. It should be noted at the outset that the method developed
in this section minimizes an equation-error measure and not the usual
frequency-response error measure.
The number of samples specified, L+1L+1, will be made larger than
the number of filter coefficients, M+N+1M+N+1. This means that H2H2 is
rectangular and, therefore, ((Reference)) cannot in general be satisfied. To
formulate an approximation problem, a length-(L+1)(L+1) error vector
εε is introduced in ((Reference)) and ((Reference)) to give
b
0
=
H
0
a
+
ε
b
0
=
H
0
a
+
ε
(30)
Equation (Equation 28) becomes
h
1
-
ε
=
-
H
2
a
*
h
1
-
ε
=
-
H
2
a
*
(31)
where now H2H2 is rectangular with (L-M)>N(L-M)>N. Using the same
methods as used to derive (Equation 28), the error εε is minimized in a
least-squared error sense by the solution of the normal equations [16]
H
2
T
h
1
=
-
H
2
T
H
2
a
*
H
2
T
h
1
=
-
H
2
T
H
2
a
*
(32)
If the equations are not singular, the solution is
a
*
=
-
[
H
2
T
H
2
]
-
1
H
2
T
h
1
.
a
*
=
-
[
H
2
T
H
2
]
-
1
H
2
T
h
1
.
(33)
If the normal equations are singular, the pseudo-inverse [16], (Reference)
can be used to obtain a minimum norm or reduced order solution.
The numerator coefficients are found by the same techniques as before in
(Equation 29)
b
=
H
1
a
b
=
H
1
a
(34)
which results in the upper M+1M+1 terms in εε being zero
and the total squared equation error being minimum.
As is true for LS-error design of FIR filters, (Equation 32) is often numerically
ill-conditioned and (Equation 33) should not be used to solve for a*a*.
Special algorithms such as those used by Matlab and LINPACK
[18], [7] should be employed.
Various modifications can be made to the form of Prony's
method presented. After the denominator is found by minimizing
the equation error, the numerator can be found by minimizing the
solution error. It is possible to mix the exact and approximate
methods. The details can be found in [19], [3], [4].
Several modifications to Prony's method have been made to
use it to minimize the solution error. Most of these iteratively
minimize a weighted-equation error with Prony's method and update
the weights from the previous determination of aa[9], [27].
If an LS-error, time-domain approximation is the desired
result, a minimization technique can be applied directly to the
solution error. The most successful method seems to be the Gauss-
Newton algorithm with a step-size control. This combined with
Prony's method to find starting parameters is an effective design
tool.
-
Berchin, Greg. (2007, January). Precise Filter Design. IEEE Signal Processing Magazine, 24(1), 137–139.
-
Burrus, C. S. and Parks, T. W. (1970, June). Time Domain Design of Recursive Digital Filters. [also in IEEE Press DSP Reprints, 1972]. IEEE Transactions on Audio and Electroacoustics, 18(2), 137–141.
-
Brophy, F. and Salazar, A. C. (1973). Considerations of the Pade Approximant Technique in the Synthesis of Recursive Digital Filters. IEEE Transactions on Audio and Electroacoustics, 21, 500–505.
-
Brophy, F. and Salazar, A. C. (1974). Recursive digital filter Systhesis in the Time Domain. IEEE Transactions on Acoustics, Speech, and Signal Processing, 22, 45–55.
-
Burrus, C. Sidney and Vargas, Ricardo A. (2008, March). The Direct Design of Recursive or IIR Digital filters. In Proceedings of International Symposium on Communications, Control, and Signal Processing. ISCCSP-08, Malta
-
Cuyt, Annie and Wuytack, Luc. (1987). Nonlinear Methods in Numerical Analysis. [North–Holland Mathematics Studies: 136]. Amsterdam: Elsevier.
-
Dongarra, J. J. and Bunch, J. R. and Moler, C. B. and Stewart, G. W. (1979). LINPACK User's Guide. Philadelphia, PA: SIAM.
-
Deczky, A. G. (1972). Equiripple and Minimum Chebyshev Approximations for Recurxive Digital Filters. IEEE Transactions on Acoustics, Speech, and Signal Processing, 20, 257–263.
-
Evans, A. G. and Fischl, R. (1973). Optimal Least Squares Time-Domain Systhesis of Recursive Digital Filters. IEEE Transactions on Audio and Electroacoustics, 21, 61–65.
-
Garcia, C. B. and Zangwill, W. I. (1981). Pathways to Solutions, Fixed Points, and Equilibria. Englewood Cliffs, NJ: Prentice-Hall.
-
Jackson, L. B. (2008). Frequency-Domain Steiglitz-McBride Method for Least-Squares IIR Filter Design, ARMA Modeling, and Periodogram Smoothing. IEEE Signal Processing Letters, 15, 49–52.
-
Kumaresan, R. and Burrus, C. S. (1991, May 14–17). Fitting a Pole-Zero Filter Model to Arbitrary Frequency Response Samples. In Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing. (p. 1649–1652). IEEE ICASSP-91, Toronto
-
Kobayashi, Takao and Imai, Satoshi. (1990, April). Complex Chebyshev Approximation for IIR Digital Filters in Iterative WLS Technique. In Proceedings of the ICASSP-90. (p. 13.D3.12). Albuquerque, NM
-
Kobayashi, Takao and Imai, Satoshi. (1990, February). Design of IIR Digital Filters with Arbitrary Log Magnitude Function by WLS Techniques. IEEE Trans. on ASSP, 38(2), 247–252.
-
Levy, E. C. (1959, May). Complex-Curve Fitting. IRE Transactions on Automatic control, AC-4(1), 37–43.
-
Lawson, C. L. and Hanson, R. J. (1974). Solving Least Squares Problems. Inglewood Cliffs, NJ: Prentice-Hall.
-
McClellan, J. H. and Lee, D. W. (1991, February). Exact Equivalence of the Steiglitz–McBride Iteration and IQML. IEEE Transactions on Signal Processing, 39(2), 509–512.
-
Moler, Cleve and Little, John and Bangert, Steve. (1989). Matlab User's Guide. South Natick, MA: The MathWorks, Inc.
-
Parks, T. W. and Burrus, C. S. (1987). Digital Filter Design. New York: John Wiley & Sons.
-
Pintelon, Rik. (1991, June). Comments on Design of IIR Filters in the Complex Domain. IEEE Transactions on Signal Processing, 39(6), 1454–1455.
-
Rabiner, L. R. and Gold, B. (1975). Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall.
-
Stonick, V. L. and Alexander, S. T. (1990, April). Global Optimal IIR Filter Design and ARMA Estimation Using Homotopy Continuation Methods. In Proceedings of the ICASSP-90. (p. 13.D3.13). Albuquerque, NM
-
Shaw, Arnab K. (1995, November). Optimal Design of Digital IIR Filters by Model-Fitting Frequency Response Data. IEEE Transactions on Circuits and Systems–II, 42, 702–710.
-
Sanathanan, C. K. and Koerner, J. (1963, January). Transfer Function Synthesis as a Ratio of Two Complex Polynomials. IEEE Transactions on Automatic Control, 8, 56–58.
-
Smith, Julius O. (1983, June). Techniques for Digital Filter Design and System Identification with Application to the Violin. [IIR filter design on page 50]. Ph. D. Thesis. EE Dept., Stanford University.
-
Soewito, Atmadji W. (1990, December). Least Square Digital Filter Design in the Frequency Domain. Ph. D. Thesis. Rice University, ECE Department, Houston, TX 77251.
-
Steiglitz, K. (1970). Computer-Aided Design of Recursive Digital Filters. IEEE Transactions on Audio and Electroacoustics, 18, 123–129.
-
Vargas, Ricardo Arturo. (2008, May). Iterative Design of Digital Filters. Ph. D. Thesis. Rice University, ECE Department, Houston, TX 77251.
-
Vargas, Ricardo A. and Burrus, C. Sidney. (2001, May 7-11). On the Design of IIR Filters with Arbitrary Frequency Response. In Proceedings of IEEE International Conference on Acoustics Speech and Signal Processing. ICASSP-01, Salt Lake City