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
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
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 [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 [48], [50]. It is a particularly interesting and powerful method that should be studied and understood to be fully utilized.
Linear programming was used earlier [87], [26], [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 [83], [5], and linear programming's ability to allow a variety of constraints [80].
Still another approach to
achieving a Chebyshev approximation is to minimize the
The early work by Herrmann and Schüssler [27], [31] and the algorithm
by Hofstetter, Oppenheim, and Siegel [28], [29] posed and solved a
similar problem but they had only approximate control of
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
where
One way to minimize
If the amplitude response for odd L is expressed as a sum of
or for even L
with
If
The alternation theorem [48], [52] states that the minimum Chebyshev
error has at least
for
It is possible to pose the Chebyshev approximation problem in filter design as a linear programming optimization problem [62], [89], [81], [41]. The error definition in Equation 1 can be written as an inequality by
where the scalar
The inequalities in Equation 5 can be written as
or
which can be combined into one matrix inequality using Equation 48 from FIR Digital Filters by
If
which, together with the inequality of Equation 10, is in the form of the dual problem in linear programming [19], [43], [82].
This can be solved using the lp() command from the Optimization Toolbox with Matlab [23], the Meteor software system [80], CPlex [7], or Karmarkar's algorithm [5], [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
[5], [66], [83] as implemented in Matlab by Lang [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 [22], [36], [51], [55], [53], [54], [57], [91], [90], [93]
A very efficient algorithm which uses the results of the alternation theorem is called the Remez exchange algorithm. Remez [63], [16], [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 [44], [47], [32].
A particularly useful FIR filter design implementation of the Remez exchange is called the Parks-McClellan algorithm and is described in [50], [65], [64], [48]. It has been implemented in Fortran in [49], [64], [18], [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 [45], [13], [21], [67], [33], [24], [25], [68], [70], [69], [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 [64], [46].
Parks and McClellan formulated the basic Chebyshev FIR filter design
problem by specifying the desired amplitude response
 |
Notice that if there is no transition band, i.e.
The alternation theorem states that the optimal approximation for this
problem will have an error function with
Parks and McClellan applied the Remez exchange algorithm [50] to this
filter design problem by writing
These equations are solved for
The detailed steps of the Parks-McClellan algorithm are:
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
It is interesting to note that at each iteration, the approximation is
optimal over that set of extremal frequencies and
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
In step 2 it is suggested that the simultaneous equation of Equation 12 be
solved. Parks and McClellan [49] use a more efficient and
numerically robust method of evaluating
For the low pass filter, this formulation always has an extremal frequency
at both pass and stop band edges,
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 [31], [27], [28], [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
 |
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
 |
If one wants to fix the pass band ripple and minimize the stop band ripple
[70], equation Equation 12 is changed so that the pass band ripple is
added to the appropriate top part of the vector
Iteration of this equation will keep the pass band ripple
This algorithm [28], [29], [70] came into existence in order to design
the filters posed by Herrmann and Schüssler [31], [27] where both
the pass and stop band ripple sizes,
This algorithm is iterated as a multiple exchange, keeping the number of ripples in the pass and stop band constant, to give an optimal extra ripple filter. The location and width of the transition band is controlled only by the choice of how the number of initial ripples are divided between the pass and stop band. The final filter may not have the transition located where you want it. Indeed, no solution may exist with the desired location of the transition band.
The designs produced by the HOS algorithm are always maximum ripple but this comes with a loss of accurate control over the location of the transition band. The algorithm is not, strictly speaking, an optimization algorithm. It is an interpolation algorithm. The Chebyshev error is not minimized, the designed amplitude interpolates the specified error ripples. However, although not directly minimized, the transition band width of these designs seems to be minimized [65], [61], [64]. Extra or maximum ripple designs seem to be efficient in using all the zeros to produce small ripple size and narrow transition bands, however, the loss of accurate control over the location of the transition bands becomes even more problematic with multiple transition band designs. Perhaps some compromise methods can be devised that use some of the efficiency of the maximum ripple approximations with some of the control of other methods. The next two design methods are of that type.
Shpak and Antoniou [67] propose decoupling the size of the pass
and stopband ripple sizes in order to have control over the pass and stop
band edges and have an extra ripple design. The Parks-McClellan design
has the ripple sizes related with a fixed weight
This method has two ripple sizes,
Because the arguments in the Weisburn, Parks, and Shenoy paper
[94] require the assumption of no signal or
noise energy in the transition band, it is now
more obvious that a narrow transition band is very desirable. For this
reason it may be better to fix the pass and stop band peak error,
For problems where the signal and noise spectra are such that a specific
frequency
The optimal Chebyshev filter designed with this new algorithm
is generally not extra ripple and, therefore, will have an extremal
frequency at
The exchange algorithm is done as by Parks and McClellan finding new
extremal frequencies at each iteration, but with fixed ripple sizes in both
pass and stop bands. This new algorithm reduces the transition band
width as done by the Hofstetter, Oppenheim, and Siegel method but with the
transition band location controlled and without requiring the extra ripple
solution. Note that any transition band frequency could be fixed. It
could be
Extending this formulation and algorithm to the multiple transition band case complicates the problem as the solution may not be unique or may have anomalous behavior in one of the transition bands. Details of the solution to this problem are given in [70].
All of the design methods discussed so far have assumed that
Fortunately, emperical formulas have been derived that give a good estimate of the relationship of the length of an optimal Chebyshev FIR filter for given pass and stopband ripple and transition band edges [64], [65]. Kaiser's formula is
and it is fairly accurate for average filter specifications (neither wide nor narrow bands).
In order to better understand the nature of an optimal Chebyshev and to see the power of the Parks-McClellan algorithm, we present the design of a length-21 linear phase FIR bandpass filter. To see the effects of the design specifications, we will fix the two pass band edges and the upper stop band edge, then look at the effects of varying the lower stop band edge. The Matlab program that generated the designs is:
% ChebyPlot9.m generates Chebyshev figures.
% Change in opt frequency response as band edge is changed, csb 1/26/07
N = 20;
M = [0 0 1 1 0 0];
W = [7.5 10 7.5];
ff = [0:512]/512; k=0;
%for fk = .10:.02:.34
% k = k+1;
clf;
for k = 1:6
fk = .1 + .02*(k-1);
F = [0 fk .35 .8 .85 1];
b = firpm(N,F,M,W);
%clf;
axis([0 1 0 1.2]);
AA = abs(fft(b,1024)); AA = AA(1:513);
dd = max(AA(1:50));
ddd = dd*(W(1)/W(2));
subplot(3,2,k); plot(ff,AA,'r'); hold;
plot([0 F(2) F(2) F(5) F(5) 1],[dd dd 0 0 dd dd],'b');
plot([0 F(3) F(3) F(4) F(4) 1],[0 0 1-ddd 1-ddd 0 0],'b');
plot([0 F(3) F(3) F(4) F(4) 1],[0 0 1+ddd 1+ddd 0 0],'b');
title('L-21 Chebyshev Filter, f_s = 0.1');
ylabel('Magnitude |H(\omega)|');
pause;
end; hold off;
The results are shown in Figures Figure 5 and Figure 6.
 |
 |
Note the large transmission peaks in the transition band of Figures Figure 5a, b, and c that result from the two transition bands being very different in width. As the lower transition band narrows, this peak grows smaller and eventually disappears in Figure 5f. Note that there are two extremal points in the lower stop band of Figure 5b and seven in the pass band, while there are three in the lower stop band of Figures Figure 5c and d and six in the pass band. But, there are always twelve total (thirteen for a case between Figures Figure 5b and c). In Figure 6d, there are only five extremal points in the pass band but twelve total. The same filter is optimal for the conditions given in Figures Figure 6a, b, and c. Much can be learned about optimal filters by running experiments in Matlab. Remember, all of these are optimal for the specifications given.
It is possible to approximate the effects of Chebyshev approximation by
minimizing the
Examples of expected and unexpected results of optimality. Rabiner's work will be used here. The non-unique designs for certain multiband designs will be illustrated.
Algorithms that directly use the alternation theorem, such as the standard Remez multiple exchange algorithm, are difficult to apply to the complex approximation or 2-D approximation problem because the concept of “alternation" is difficult to define and the number of ripples in an optimal solution is more difficult to determine [92], [79], [84], [12], [40], [39], [85]. Work has been done on the complex approximation problem at Rice by Parks and Chen [17] and by Burrus, Barreto, and Selesnick [11], [9], at Erlangen by Schuessler, Preuss, Schulist, and Lang [59], [60], [71], [73], [74], [72], at MIT by Alkhairy et al [1], [2], at USC by Tseng and Griffiths [88], [86], at Georgia Tech by Karam and McClellan [34], at Cornell by Burnside and Parks [14], and by Potchinkov and Reemtsen at Cottbus [53], [54], [57], [58], [56]. The work done by Adams which uses an implementation of a constrained quadratic programming algorithm might be useful here [3], [6]. Lang has extended and further developed this constrained approach [37], [38], [42] and Selesnick is applying it to IIR filter design [75]. Tseng gives a good summary of complex approximation in [86].
By adding the Chebyshev filter design methods described above to the Parks-McClellan algorithm, one has a rather complete set of approaches to equal ripple filter designs that allows a wide variety of specifications. The new exchange algorithm which minimizes the transition band width while allowing the specification of the center or either edge of the transition band edge may fit many design environments better than the traditional Parks-McClellan. An alternative approach which specifies the pass and stop band peak error yet has no zero weighted transition band will be presented in Constrained Least Squares Design[77], [78]. Matlab programs are available for the Parks-McClellan algorithm, the modified Parks-McClellan algorithm, the Hofstetter-Oppenheim-Siegel algorithm, the new minimum transition band design algorithm, and the constrained least squares algorithm. They are written with a common format and notation to easily see how they are programmed and how they are related. This book generally presents the lowpass case. The bandpass and multi-band cases use the same ideas but are a bit more complicated and are discussed in more detail in the references.
More about this module: Metadata | Downloads | Version History
This work is licensed by Daniel Williamson under a Creative Commons Attribution License (CC-BY 3.0), and is an Open Educational Resource.
Last edited by C. Sidney Burrus on Nov 17, 2012 5:58 pm -0600.