If one poses the FIR filter design problem by requiring the maximum error
over certain bands of frequencies be minimized, we call the resulting
filter a Chebyshev filter or an equal ripple filter. The fact that the
minimization of the Chebyshev or L∞L∞ error results in an equal
ripple error comes from the alternation theorem. This very powerful
theorem allows one to minimize the Chebyshev error by directly
constructing an equal ripple approximation with the proper number of
ripples. That is the basis of several very effective algorithms, including the
Remez exchange algorithm.
There are several ways one could pose the Chebyshev FIR filter design
problem. For a simple length-N linear phase, lowpass filter with a transition band, if
one considers the length N, the passband ripple δpδp, the stopband
ripple δsδs, and the transition bandwidth Δ=ωs-ωpΔ=ωs-ωp, then one can fix or constrain any three of them and minimize
the fourth. Or, as Parks and McClellan do, fix the band edges, ωpωp and
ωsωs, and the ratio of δpδp
and δsδs and minimize one of them.
The Chebyshev error measure is often used for approximation in digital
filter design. This is particularly true when the signals and/or noise
are narrow band or single frequency or when one wants to minimize worst
case possibilities. Theoretical justification for its use has been given
by Weisburn, Parks, and Shenoy Entry 94. For FIR filter design,
the Parks-McClellan formulation of the filter design problem and
application of the Remez exchange algorithm is most commonly used
Entry 48, Entry 50. It is a particularly interesting and powerful method that
should be studied and understood to be fully utilized.
Linear programming was used earlier Entry 87, Entry 26, Entry 62 but
dropped out of favor when the Parks-McClellan algorithm was introduced.
It is now becoming more popular again because of more powerful computers,
better algorithms Entry 83, Entry 5, and linear programming's ability to
allow a variety of constraints Entry 80.
Still another approach to
achieving a Chebyshev approximation is to minimize the pthpth power of
the error using a large value of pp or to use an iterative scheme that
solves a weighted least squared error with the weights at each stage
determined by the error of the previous stage Entry 11. Still another
design method that produces an equal ripple error approximation uses a
constrained least squared error criterion Entry 78, Entry 76 which results in a
Chebyshev solution if tight constraints are imposed.
The early work by Herrmann and Schüssler Entry 27, Entry 31 and the algorithm
by Hofstetter, Oppenheim, and Siegel Entry 28, Entry 29 posed and solved a
similar problem but they had only approximate control of ωoωo (or
ωpωp or ωsωs) and always achieved the “extra ripple" design.
Given the proper specifications, the Parks-McClellan algorithm could
design any filter that the Hofstetter-Oppenheim-Siegel algorithm could,
but the opposite is not true. This seems to be one of the reasons the
Hofstetter-Oppenheim-Siegel algorithm is not commonly used.
The Chebyshev error is defined as the maximum difference between the
actual and desired response over a band or several bands of frequencies. This is
ε
=
max
ω
∈
Ω
|
A
(
ω
)
-
A
d
(
ω
)
|
ε
=
max
ω
∈
Ω
|
A
(
ω
)
-
A
d
(
ω
)
|
(1)
where ΩΩ is the union of the bands of frequencies that the
approximation is over Entry 16, Entry 20. The approximation problem in
filter design is to choose the filter coefficients to minimize
εε.
One way to minimize εε is to set up the frequency response in
four equations for the four types of linear phase FIR filters as done in
((Reference)), ((Reference)), and the corresponding sine expressions. An
alternative approach Entry 48 uses the fact that all four can be obtained
from the odd-length, even-symmetry type 1 and uses only ((Reference)). From
one of these frequency response representations together with powerful
Alternation Theorem several optimization schemes can be developed.
If the amplitude response for odd L is expressed as a sum of RR cosine terms
A
(
ω
)
=
∑
n
=
0
R
-
1
a
(
n
)
cos
(
ω
n
)
A
(
ω
)
=
∑
n
=
0
R
-
1
a
(
n
)
cos
(
ω
n
)
(2)
or for even L
A
(
ω
)
=
∑
n
=
1
R
a
(
n
)
cos
(
ω
(
n
-
1
/
2
)
)
A
(
ω
)
=
∑
n
=
1
R
a
(
n
)
cos
(
ω
(
n
-
1
/
2
)
)
(3)
with R=M+1=L+12R=M+1=L+12 for odd length-LL and R=L/2R=L/2 for even length-LL,
as derived in ((Reference)) and ((Reference)), then
Theorem 1
If A(ω)A(ω) is the linear combination of RR cosine functions given
in (Equation 2) or (Equation 3), the necessary and sufficient conditions for
A(ω)A(ω) to be the least Chebyshev error approximation to Ad(ω)Ad(ω)
over ω∈Ωω∈Ω are: The error function, ϵ(ω)=A(ω)-Ad(ω)ϵ(ω)=A(ω)-Ad(ω) have at leastR+1R+1 extremal frequencies in ΩΩ.
The extremal frequencies are ordered points ω1<ω2<⋯<ωR+1ω1<ω2<⋯<ωR+1 such that ϵ(ωk)=-ϵ(ωk+1)ϵ(ωk)=-ϵ(ωk+1)
and maxω∈Ω|ϵ(ω)|=|ϵ(ωk)|maxω∈Ω|ϵ(ω)|=|ϵ(ωk)| for k=1,2,⋯,R+1k=1,2,⋯,R+1.
The alternation theorem Entry 48, Entry 52 states that the minimum Chebyshev
error has at least R+1R+1 extremal frequencies. This is stated
mathematically by
A
(
ω
k
)
=
A
d
(
ω
k
)
+
(
-
1
)
k
δ
A
(
ω
k
)
=
A
d
(
ω
k
)
+
(
-
1
)
k
δ
(4)
for k=0,1,2,⋯,Rk=0,1,2,⋯,R, where the ωkωk are the ordered
extremal frequencies where the equal ripple error has maximum value. In
other words, the optimal solution to the linear phase FIR filter design
problem will have an equal ripple error function with a required number of
ripples. How all of these characteristics relate can be rather
complicated and good designs require experience Entry 30. When
applied to other approximation problems, care must be taken to ensure the
approximating functions satisfy the “Haar conditions" or other
restrictions Entry 16, Entry 50, Entry 20, Entry 52.
It is possible to pose the Chebyshev approximation problem in filter
design as a linear programming optimization problem
Entry 62, Entry 89, Entry 81, Entry 41. The error definition in (Equation 1) can be
written as an inequality by
A
d
(
ω
)
-
δ
≤
A
(
ω
)
≤
A
d
(
ω
)
+
δ
A
d
(
ω
)
-
δ
≤
A
(
ω
)
≤
A
d
(
ω
)
+
δ
(5)
where the scalar δδ is minimized.
The inequalities in (Equation 5) can be written as
A
≤
A
d
+
δ
A
≤
A
d
+
δ
(6)
-
A
≤
-
A
d
+
δ
-
A
≤
-
A
d
+
δ
(7)
or
A
-
δ
≤
A
d
A
-
δ
≤
A
d
(8)
-
A
-
δ
≤
-
A
d
-
A
-
δ
≤
-
A
d
(9)
which can be combined into one matrix inequality using ((Reference)) by
C
-
1
-
C
-
1
a
δ
≤
A
d
-
A
d
.
C
-
1
-
C
-
1
a
δ
≤
A
d
-
A
d
.
(10)
If δδ is minimized, the optimal Chebyshev approximation is achieved.
This is done by minimizing
ε
=
0
0
⋯
1
a
δ
ε
=
0
0
⋯
1
a
δ
(11)
which, together with the inequality of (Equation 10), is in the form of the
dual problem in linear programming Entry 19, Entry 43, Entry 82.
This can be solved using the lp()
command from the Optimization
Toolbox with Matlab Entry 23, the Meteor software system Entry 80,
CPlex Entry 7, or Karmarkar's algorithm Entry 5, Entry 35. The
Matlab lp
command is implemented in an m-file using a form of
quadratic programming algorithm that is not well suited to our filter
design problem. Meteor is a linear programming system using the simplex
algorithm written in Pascal by Ken Steiglitz especially for filter design.
It has been compiled on a variety of computers and efficiently designs
filters over 100 in length. The Karmarkar program written by Lang is a
relatively short m-file that is not particularly fast but is robust and
can design filters on the order of length-100. CPlex is a proprietary
program that can be used alone or called from Fortran programs and is
particularly robust and fast.
A Matlab program that applies its linear programming function
lp.m
to (Equation 10,Equation 11) for linear phase FIR filter design is
given by:
% lpdesign.m Design an FIR filter from L, f1, f2, and LF using LP.
% L is filter length, f1 and f2 are pass and stopband edges, LF is
% the number of freq samples. L is odd. Uses lp.m
% csb 5/22/91
L1 = fix(LF*f1/(.5-f2+f1)); L2 = LF - L1; %No. freq samples in PB, SB
Ad = [ones(L1,1); zeros(L2,1)]; %Samples of ideal response
f = [[0:L1-1]*f1/(L1-1), ([0:L2-1]*(.5-f2)/(L2-1) + f2)]'; %Freq samples
M = (L-1)/2;
C = cos(2*pi*(f*[0:M])); %Freq response matrix
CC = [C, -ones(LF,1); -C, -ones(LF,1)]; %LP matrix
AD = [Ad; -Ad];
c = [zeros(M+1,1);1]; %Cost function
x0 = [zeros(M+1,1);max(AD)+1]; %Starting values
x = lp(c,CC,AD,[],[],x0); %Call the LP
d = x(M+2); %delta or deviation
a = x(1:M+1); %Half impulse resp.
h = [a(M+1:-1:2);2*a(1);a(2:M+1)]./2; %Impulse response
This program has numerical problems for filters longer than 10 or 20 and
is fairly slow. The lp()
function uses an algorithm that seems not
well suited to the equations required by filter design. It would be nice
to have Meteor written in Matlab, both to show how the Simplex algorithm
works, and to have an efficient LP filter design system in Matlab. The
above program has been tested using Karmarkar's algorithm
Entry 5, Entry 66, Entry 83 as implemented in Matlab by Lang Entry 35.
It proved to be robust and reliable for lengths up to 100 or more. It was
faster than the Matlab function but slower than Meteor or CPlex. Its use
should be further investigated.
Direct use of quadratic programming and other optimization algorithms seem
promising
Entry 22, Entry 36, Entry 51, Entry 55, Entry 53, Entry 54, Entry 57, Entry 91, Entry 90, Entry 93
A very efficient algorithm which uses the results of the alternation
theorem is called the Remez exchange algorithm. Remez
Entry 63, Entry 16, Entry 52 showed that, under rather general conditions,
an algorithm that takes a starting estimate of the location of the
extremal frequencies and exchanges them with a new set calculated at each
iteration will converge to the optimal Chebyshev approximation. The
efficiency of this algorithm comes from finding the optimal solution by
directly constructing a function that satisfies the alternation theorem
rather than minimizing the Chebyshev error as done by the linear
programming technique. The Remez exchange algorithm has proven to be well
suited to the design of linear phase FIR filters Entry 44, Entry 47, Entry 32.
A particularly useful FIR filter design implementation of the Remez exchange
is called the Parks-McClellan algorithm and is described in
Entry 50, Entry 65, Entry 64, Entry 48. It has been implemented in Fortran in
Entry 49, Entry 64, Entry 18, Entry 48 and in Matlab in a program at the end of this
material. The Matlab program is particularly helpful in understanding how
the algorithm works, however, because it does not use any special tricks,
it is limited to lengths of 60 or so. Extensions and details can be found
in Entry 45, Entry 13, Entry 21, Entry 67, Entry 33, Entry 24, Entry 25, Entry 68, Entry 70, Entry 69, Entry 4.
This is a robust, efficient algorithm that significantly changed DSP when
Parks and McClellan first described it in 1972 and has undergone important
improvements. Examples are illustrated in Entry 64, Entry 46.
Parks and McClellan formulated the basic Chebyshev FIR filter design
problem by specifying the desired amplitude response A(ω)A(ω) and the
transition band edges, then minimizing the weighted Chebyshev error over
the pass and stop bands. For the basic lowpass filter illustrated in
Figure 1, the pass band edge ωpωp and the stop band edge ωsωs
are specified, the maximum passband error is related to the maximum stop
band error by δp=Kδsδp=Kδs and they are minimized.
Notice that if there is no transition band, i.e. ωp=ωsωp=ωs,
that δp+δs=1δp+δs=1 and no minimization is possible. While not
the case for a least squares approximation, a transition band is necessary
for the Chebyshev approximation problem to be well-posed. The effects of
a small transition band are large pass and stopband ripple as illustrated
in Figure 3b.
The alternation theorem states that the optimal approximation for this
problem will have an error function with R+1R+1 extremal points with
alternating signs. The theorem also states that there exists R+1R+1
frequencies such that, if the Chebyshev error at those frequencies are
equal and alternate in sign, it will be minimized over the pass band and
stop band. Note that there are nine extremal points in the length-15 example
shown in Figure 1, counting those at the band edges in addition
to those that are interior to the pass and stopbands. For this case,
R=(L+1)/2R=(L+1)/2 which agree with the example.
Parks and McClellan applied the Remez exchange algorithm Entry 50 to this
filter design problem by writing R+1R+1 equations using ((Reference)) and
(Equation 4) evaluated at the R+1R+1 extremal frequencies with RR unknown
cosine parameters a(n)a(n) and the unknown ripple value, δδ. In
matrix form this becomes
A
d
(
ω
0
)
A
d
(
ω
1
)
A
d
(
ω
2
)
A
d
(
ω
3
)
⋮
A
d
(
ω
R
)
=
cos
(
ω
0
0
)
cos
(
ω
0
1
)
⋯
cos
(
ω
0
(
R
-
1
)
)
1
cos
(
ω
1
0
)
cos
(
ω
1
1
)
⋯
cos
(
ω
1
(
R
-
1
)
)
-
1
cos
(
ω
2
0
)
cos
(
ω
2
1
)
⋯
cos
(
ω
2
(
R
-
1
)
)
1
cos
(
ω
3
0
)
cos
(
ω
3
1
)
⋯
cos
(
ω
3
(
R
-
1
)
)
-
1
⋮
⋮
cos
(
ω
R
0
)
cos
(
ω
R
1
)
⋯
cos
(
ω
R
M
)
±
1
a
(
0
)
a
(
1
)
a
(
2
)
⋮
a
(
R
-
1
)
δ
.
A
d
(
ω
0
)
A
d
(
ω
1
)
A
d
(
ω
2
)
A
d
(
ω
3
)
⋮
A
d
(
ω
R
)
=
cos
(
ω
0
0
)
cos
(
ω
0
1
)
⋯
cos
(
ω
0
(
R
-
1
)
)
1
cos
(
ω
1
0
)
cos
(
ω
1
1
)
⋯
cos
(
ω
1
(
R
-
1
)
)
-
1
cos
(
ω
2
0
)
cos
(
ω
2
1
)
⋯
cos
(
ω
2
(
R
-
1
)
)
1
cos
(
ω
3
0
)
cos
(
ω
3
1
)
⋯
cos
(
ω
3
(
R
-
1
)
)
-
1
⋮
⋮
cos
(
ω
R
0
)
cos
(
ω
R
1
)
⋯
cos
(
ω
R
M
)
±
1
a
(
0
)
a
(
1
)
a
(
2
)
⋮
a
(
R
-
1
)
δ
.
(12)
These equations are solved for a(n)a(n) and δδ using an initial guess
as to the location of the extremal frequencies ωiωi. This design is
optimal but only over the guessed frequencies, and we want optimality over
all of the pass and stopbands. Therefore, the amplitude response of the
filter is calculated over a dense set of frequency samples using
((Reference)) and a new set of estimates of the extremal frequencies is found
from the local minima and maxima and these are used to replace the initial
guesses (they are exchanged). This process is iteratively performed until
the guaranteed convergence is achieved and the optimal filter is designed.
The detailed steps of the Parks-McClellan algorithm are:
- Specify the ideal amplitude, Ad(ω)Ad(ω), the stop and pass band edges,
ωpωp and ωsωs, the error weight KK where δp=Kδsδp=Kδs, and the length LL.
- Choose R+1R+1 initial guesses for the extremal frequencies, ωiωi, in
the bands of approximation, ΩΩ.
This is often done uniformly over the pass and stop bands, including
ω=0,ωp,ωs,ω=0,ωp,ωs, and ππ.
- Calculate the cosine matrix at the current ωiωi and solve (Equation 12) for
a(n)a(n) and δδ which are optimal over these current extremal frequencies,
ωiωi.
- Using the a(n)a(n) or the equivalent h(n)h(n) from step 3, evaluate A(ω)A(ω)
over a dense set of frequencies. This amplitude response will interpolate
A(ωi)=Ad(ωi)±δA(ωi)=Ad(ωi)±δ at the extremal frequencies.
- Find R+1R+1 new extremal frequencies where the error has a local maximum
or minimum and has alternating sign. This includes the band edges.
- If the largest error is the same as δδ found in step 3, then
convergence has occured and the optimal filter has been designed, otherwise,
exchange the old extremal frequencies ωiωi used in step 2 and
return to step 3 for the next iteration.
- This iterative algorithm is guaranteed to converge to the unique optimal
solution using almost any starting points in step 2.
This iterative procedure is called a multiple exchange algorithm because
all of the extremal frequencies are up-dated each iteration. If only the
frequency of the largest error is up-dated each iteration, it is called a
single exchange algorithm which also converges but much more slowly. Some
modification of the Parks-McClellan method or the Remez exchange algorithm
will not converge as a multiple exchange, but will as a single exchange.
The Alternation theorem states that there will be a minimum of R+1R+1
extremal frequencies, even for multiband designs with arbitrary
Ad(ω)Ad(ω). If Ad(ω)Ad(ω) is piece-wise constant with TT
transition bands, one can derive the maximum possible number of extremal
frequencies and it is R+2TR+2T. This comes from the maximum number
of maxima and minima that a function of the form (Equation 2) or
(Equation 3) can have plus two at the edges of each transition band. For a
simple lowpass filter with one passband, one transition band, and one
stopband, there will be a minimum of R+1R+1 extremal frequencies and a
maximum of R+2R+2. For a bandpass filter, the maximum is R+4R+4. If a design has more than the minimum number of extremal frequencies, it is
called an extra ripple design. If it has the maximum number, it is
called a maximum ripple design.
It is interesting to note that at each iteration, the approximation is
optimal over that set of extremal frequencies and δδ increased
over the previous iteration. At convergence, δδ has increased to
the maximum error over ΩΩ and that is the minimum Chebyshev error.
At each iteration, the exchange of a proper set of extremal frequencies
with alternating signs of the errors is always possible. One can show
there will never be too few and if there are too many, one uses those
corresponding to the largest errors.
In step 4 it is suggested that the amplitude response A(ω)A(ω) be
calculated over a dense grid in the pass and stopbands and in step 5
the local extremes are found by searching over this dense grid. There
are more accurate methods that use bisection methods and/or Newton's
method to find the extremal points.
In step 2 it is suggested that the simultaneous equation of (Equation 12) be
solved. Parks and McClellan Entry 49 use a more efficient and
numerically robust method of evaluating δδ using a form of Cramer's
rule. With that δδ, an interpolation method can be used to find
a(n)a(n). This is faster and allows longer filters to be designed than with
the linear algebra based approach described here.
For the low pass filter, this formulation always has an extremal frequency
at both pass and stop band edges, ωpωp and ωsωs, and at
ω=0ω=0 and/or at ω=πω=π. The extra ripple filter has R+2R+2
extremal frequencies including both zero and pi. If this algorithm is
started with an incorrect number of extremal frequencies in the stop or
pass band, the iterations will correct this. It is interesting and
informative to plot the frequency response of the filters designed at each
iteration of this algorithm and observe how the correction takes place.
The Parks-McClellan algorithm starts with fixed pass and stop band edges
then minimizes a weighted form of the pass and stop band error ripple. In
some cases it may be more appropriate to fix one of the ripples and
minimize the other or to fix both ripples and minimize the transition band
width. Indeed Schüssler, Hofstetter, Tufts, and others Entry 31, Entry 27, Entry 28, Entry 29
formulated some of these ideas before Parks and McClellan developed their
algorithm. The DSP group at Rice has developed some modifications to
these methods and they are presented below.
Here we look at several examples of filters designed by the Parks-McClellan algorithm.
The examples here are length-15 with that shown in Figure 2a having a passband
0<f<0.30<f<0.3, a transition band 0.3<f<0.50.3<f<0.5, and a stopband 0.5<f<10.5<f<1. The
number of cosine terms in the frequency response formula is R=8R=8, therefore, the
alternation theorem says we must have at least R+1R+1 extremal points. There are four
in the passband, counting the one at zero frequency, the minimum, the maximum, and
the minimum at the bandedge. There are five in the stopband, counting the ones at the
bandedge and at f=1f=1. So, the number is nine which is at least R+1R+1. However, in
Figure 2c, there are ten extremal points but that is also at least 9, so it
also is optimal. For a low pass filter, the maximum number of extremal points is R+2R+2
and that is what this filter has. This special case is called the “maximum ripple" case.
It is possible to have ripples that do not touch the maximum value and, therefore,
are not considered extremal points. That is illustrated in Figure 3a. The
effects of a narrow transitionband are illustrated in Figure 3c. Note
the zero locations for these filters and how they relate to the amplitude response.
To illustrate some of the unexpected behavior that optimal filter designs can have,
consider the bandpass filter amplitude response shown in Figure 4. Here
we have a length-31 Chebyshev bandpass filter with a stopband
0<f<0.20<f<0.2, a transition band 0.2<f<0.250.2<f<0.25, a passband 0.25<f<0.50.25<f<0.5,
another transitionband 0.5<f<0.680.5<f<0.68, and a stopband 0.68<f<10.68<f<1. The
asymmetric transition bands cause large response in the transition band around
f=0.6f=0.6. However, this filter is optimal since the deviation occurs in part
of the frequency band that is not included in the optimization criterion. If
you think you don't care what happens in the transition bands, you may change
your mind with this kind of behavior.
If one wants to fix the pass band ripple and minimize the stop band ripple
Entry 70, equation (Equation 12) is changed so that the pass band ripple is
added to the appropriate top part of the vector AdAd