In this section a frequencysampling design method is
developed such that the frequency response of the IIR filter will interpolate or
pass through the given samples of a desired response. This development is
parallel to that used for Prony's method in the time domain.
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 LSerror approximation can be done as
for the FIR filter [57]. 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.
Nevertheless, it is a useful noniterative design method. Finally,
a general discussion of iterative design methods for LSfrequency
response error is given.
The method for calculating samples of the frequency response
of an IIR filter can be reversed to design
a filter much the same way it was for the FIR filter using frequency sampling [57].
The ztransform 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
.
(18)The frequency response of the filter is given by setting z=ejωz=ejω. Using the notation
H
(
ω
)
=
H
(
z
)

z
=
e

j
ω
.
H
(
ω
)
=
H
(
z
)

z
=
e

j
ω
.
(19)Equallyspaced samples of the frequency response are chosen so
that the number of samples is equal to the number of unknown
coefficients in Equation 18. 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
)
(20)and can be calculated from the length(L+1)(L+1)) DFTs of the numerator
and denominator (padded with zeros to the proper length).
H
k
=
DFT
{
b
n
}
DFT
{
a
n
}
=
B
k
A
k
H
k
=
DFT
{
b
n
}
DFT
{
a
n
}
=
B
k
A
k
(21)where the indicated division is termbyterm division for each value
of kk. Multiplication of both sides of Equation 21 by AkAk gives
B
k
=
H
k
A
k
B
k
=
H
k
A
k
(22)If the length(L+1)(L+1) inverse DFT of HkHk is denoted by the length
(L+1)(L+1) sequence hnhn, equation Equation 22 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
(23)Note that the hnhn in Equation 23 are not the impulse response values
of the filter (they are an aliased version of it) as used in the FIR case or
in Equation 6. Using the same approach as used for Prony's method in the time
domain, a more compact matrix notation is
b
0
=
H
a
0
b
0
=
H
a
0
(24)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 LNLN terms of
the righthand vector of Equation 23 are zero, the HH matrix can be
reduced by deleting the rightmost LNLN columns to give H0H0
which causes Equation 24 to become
b
0
=
H
0
a
b
0
=
H
0
a
(25)Because the first element of aa is unity, it is partitioned
to remove the unity term and the remaining lengthNN vector is
denoted a*a*. The simultaneous equations represented by Equation 25
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
*
(26)where H1H1 is (M+1)(M+1) by (N+1)(N+1), h1h1 is
length(LM)(LM), and H2H2 is (LM)(LM) by NN. The lower (LM)(LM)
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
(10) are written
which allows the calculation of bb.
If L=N+ML=N+M, 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 [42], [11] and there
are multiple solutions, a lower order problem can be posed. If
there are no solutions, the approximation methods must be
used and/or the assumed order increased.
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
[57]. Also note that there is no control
over the stability of the filter designed by this method.
This approach uses of the DFT therefore does not allow the
possibility of unequally spaced frequency samples as was possible
for FIR filter design.
The frequencysampling 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 frequencyresponse samples are not
consistent with what an IIR filter can achieve. One solution to
this problem is the same as for the FIR case
[57], the use of more frequency samples than the number of
filter coefficients and the definition of an approximation error
function that can be minimized. Another solution is choose another
desired frequency that is closer to what can be achieved.
There is also no simple restriction that
will guarantee stable filters. If the frequencyresponse samples
are consistent with an unstable filter, that is what will be
designed.