Optimal FIR Digital Filter Design
We now look briefly at the topic of optimal FIR filter design. We saw examples above of optimal Chebyshev designs (§4.5.2). and an oversimplified optimal leastsquares design (§4.3). Here we elaborate a bit on optimality formulations under various error norms.
Lp norms
The norm of an dimensional vector (signal) is defined as
(5.27) 
Special Cases

norm
(5.28)
 Sum of the absolute values of the elements
 ``City block'' distance

norm
(5.29)
 ``Euclidean'' distance
 Minimized by ``Least Squares'' techniques

norm
In the limit as , the norm of is dominated by the maximum element of . Optimal Chebyshev filters minimize this norm of the frequencyresponse error.
Filter Design using Lp Norms
Formulated as an norm minimization, the FIR filter design problem can be stated as follows:
(5.31) 
where
 FIR filter coefficients
 suitable discrete set of frequencies
 desired (complex) frequency response
 obtained frequency response (typically fft(h))
 (optional) error weighting function
Optimal Chebyshev FIR Filters
As we've seen above, the defining characteristic of FIR filters optimal in the Chebyshev sense is that they minimize the maximum frequencyresponse errormagnitude over the frequency axis. In other terms, an optimal Chebyshev FIR filter is optimal in the minimax sense: The filter coefficients are chosen to minimize the worstcase error (maximum weighted errormagnitude ripple) over all frequencies. This also means it is optimal in the sense because, as noted above, the norm of a weighted frequencyresponse error is the maximum magnitude over all frequencies:
(5.32) 
Thus, we can say that an optimal Chebyshev filter minimizes the norm of the (possibly weighted) frequencyresponse error. The norm is also called the uniform norm. While the optimal Chebyshev FIR filter is unique, in principle, there is no guarantee that any particular numerical algorithm can find it.
The optimal Chebyshev FIR filter can often be found effectively using the Remez multiple exchange algorithm (typically called the ParksMcClellan algorithm when applied to FIR filter design) [176,224,66]. This was illustrated in §4.6.4 above. The ParksMcClellan/Remez algorithm also appears to be the most efficient known method for designing optimal Chebyshev FIR filters (as compared with, say linear programming methods using matlab's linprog as in §3.13). This algorithm is available in Matlab's Signal Processing Toolbox as firpm() (remez() in (Linux) Octave).^{5.13}There is also a version of the Remez exchange algorithm for complex FIR filters. See §4.10.7 below for a few details.
The Remez multiple exchange algorithm has its limitations, however. In particular, convergence of the FIR filter coefficients is unlikely for FIR filters longer than a few hundred taps or so.
Optimal Chebyshev FIR filters are normally designed to be linear phase [263] so that the desired frequency response can be taken to be real (i.e., first a zerophase FIR filter is designed). The design of linearphase FIR filters in the frequency domain can therefore be characterized as real polynomial approximation on the unit circle [229,258].
In optimal Chebyshev filter designs, the error exhibits an equiripple characteristicthat is, if the desired response is and the ripple magnitude is , then the frequency response of the optimal FIR filter (in the unweighted case, i.e., for all ) will oscillate between and as increases. The powerful alternation theorem characterizes optimal Chebyshev solutions in terms of the alternating error peaks. Essentially, if one finds sufficiently many for the given FIR filter order, then you have found the unique optimal Chebyshev solution [224]. Another remarkable result is that the Remez multiple exchange converges monotonically to the unique optimal Chebyshev solution (in the absence of numerical roundoff errors).
Fine online introductions to the theory and practice of Chebyshevoptimal FIR filter design are given in [32,283].
The window method (§4.5) and Remezexchange method together span many practical FIR filter design needs, from ``quick and dirty'' to essentially ideal FIR filters (in terms of conventional specifications).
LeastSquares LinearPhase FIR Filter Design
Another versatile, effective, and oftenused case is the weighted least squares method, which is implemented in the matlab function firls and others. A good general reference in this area is [204].
Let the FIR filter length be samples, with even, and suppose we'll initially design it to be centered about the time origin (``zero phase''). Then the frequency response is given on our frequency grid by
(5.33) 
Enforcing even symmetry in the impulse response, i.e., , gives a zerophase FIR filter that we can later rightshift samples to make a causal, linear phase filter. In this case, the frequency response reduces to a sum of cosines:
(5.34) 
or, in matrix form:
Recall from §3.13.8, that the Remez multiple exchange algorithm is based on this formulation internally. In that case, the lefthandside includes the alternating error, and the frequency grid iteratively seeks the frequencies of maximum errorthe socalled extremal frequencies.
In matrix notation, our filterdesign problem can be stated as (cf. §3.13.8)
(5.36) 
where these quantities are defined in (4.35). We can denote the optimal leastsquares solution by
(5.37) 
To find , we need to minimize
This is a quadratic form in . Therefore, it has a global minimum which we can find by setting the gradient to zero, and solving for .^{5.14}Assuming all quantities are real, equating the gradient to zero yields the socalled normal equations
(5.39) 
with solution
(5.40) 
The matrix
(5.41) 
is known as the (MoorePenrose) pseudoinverse of the matrix . It can be interpreted as an orthogonal projection matrix, projecting onto the columnspace of [264], as we illustrate further in the next section.
Geometric Interpretation of Least Squares
Typically, the number of frequency constraints is much greater than the number of design variables (filter coefficients). In these cases, we have an overdetermined system of equations (more equations than unknowns). Therefore, we cannot generally satisfy all the equations, and are left with minimizing some error criterion to find the ``optimal compromise'' solution.
In the case of leastsquares approximation, we are minimizing the Euclidean distance, which suggests the geometrical interpretation shown in Fig.4.19.
Thus, the desired vector
is the vector sum of its
best leastsquares approximation
plus an orthogonal error
:
(5.42) 
In practice, the leastsquares solution can be found by minimizing the sum of squared errors:
(5.43) 
Figure 4.19 suggests that the error vector is orthogonal to the column space of the matrix , hence it must be orthogonal to each column in :
(5.44) 
This is how the orthogonality principle can be used to derive the fact that the best least squares solution is given by
(5.45) 
In matlab, it is numerically superior to use ``h= A h'' as opposed to explicitly computing the pseudoinverse as in ``h = pinv(A) * d''. For a discussion of numerical issues in matrix leastsquares problems, see, e.g., [92].
We will return to leastsquares optimality in §5.7.1 for the purpose of estimating the parameters of sinusoidal peaks in spectra.
Matlab Support for LeastSquares FIR Filter Design
Some of the available functions are as follows:
 firls  leastsquares linearphase FIR filter design
for piecewise constant desired amplitude responses  also designs
Hilbert transformers and differentiators
 fircls  constrained leastsquares linearphase FIR
filter design for piecewise constant desired amplitude responses 
constraints provide lower and upper bounds on the frequency response
 fircls1  constrained leastsquares linearphase FIR
filter design for lowpass and highpass filters  supports relative
weighting of passband and stopband error
For more information, type help firls and/or doc firls, etc., and refer to the ``See Also'' section of the documentation for pointers to more relevant functions.
Chebyshev FIR Design via Linear Programming
We return now to the norm minimization problem of §4.10.2:
and discuss its formulation as a linear programming problem, very similar to the optimal window formulations in §3.13. We can rewrite (4.46) as
(5.47) 
where denotes the th row of the matrix . This can be expressed as
s.t.  (5.48) 
Introducing a new variable
(5.49) 
then we can write
(5.50) 
and our optimization problem can be written in more standard form:
s.t.  (5.51) 
Thus, we are minimizing a linear objective, subject to a set of linear inequality constraints. This is known as a linear programming problem, as discussed previously in §3.13.1, and it may be solved using the matlab linprog function. As in the case of optimal window design, linprog is not normally as efficient as the Remez multiple exchange algorithm (firpm), but it is more general, allowing for linear equality and inequality constraints to be imposed.
More General Real FIR Filters
So far we have looked at the design of linear phase filters. In this case, , and are all real. In some applications, we need to specify both the magnitude and phase of the frequency response. Examples include
 minimum phase filters [263],
 inverse filters (``deconvolution''),
 fractional delay filters [266],
 interpolation polyphase subfilters (Chapter 11)
NonlinearPhase FIR Filter Design
Above, we considered only linearphase (symmetric) FIR filters. The same methods also work for antisymmetric FIR filters having a purely imaginary frequency response, when zerocentered, such as differentiators and Hilbert transformers [224].
We now look at extension to nonlinearphase FIR filters, managed by treating the real and imaginary parts separately in the frequency domain [218]. In the nonlinearphase case, the frequency response is complex in general. Therefore, in the formulation Eq. (4.35) both and are complex, but we still desire the FIR filter coefficients to be real. If we try to use ' ' or pinv in matlab, we will generally get a complex result for .
Problem Formulation
(5.52) 
where , , and . Hence we have,
(5.53) 
which can be written as
(5.54) 
or
(5.55) 
which is written in terms of only real variables.
In summary, we can use the standard leastsquares solvers in matlab and end up with a real solution for the case of complex desired spectra and nonlinearphase FIR filters.
Matlab for General FIR Filter Design
The cfirpm function (formerly cremez) [116,117] in the Matlab Signal Processing Toolbox performs complex FIR filter design (``Complex FIR ParksMcClellan''). Convergence is theoretically guaranteed for arbitrary magnitude and phase specifications versus frequency. It reduces to ParksMcClellan algorithm (Remez second algorithm) as a special case.
The firgr function (formerly gremez) in the Matlab Filter Design Toolbox performs ``generalized'' FIR filter design, adding support for minimumphase FIR filter design, among other features [254].
Finally, the fircband function in the Matlab DSP System Toolbox designs a variety of real FIR filters with various filtertypes and constraints supported.
This is of course only a small sampling of what is available. See, e.g., the Matlab documentation on its various toolboxes relevant to filter design (especially the Signal Processing and Filter Design toolboxes) for much more.
SecondOrder Cone Problems
In SecondOrder Cone Problems (SOCP), a linear function is minimized over the intersection of an affine set and the product of secondorder (quadratic) cones [153,22]. Nonlinear, convex problem including linear and (convex) quadratic programs are special cases. SOCP problems are solved by efficient primaldual interiorpoint methods. The number of iterations required to solve a problem grows at most as the square root of the problem size. A typical number of iterations ranges between 5 and 50, almost independent of the problem size.
Resources
 LIPSOL: Matlab code for linear programming using interior point methods.
 Matlab's linprog (in the Optimization Toolbox)
 Octave's lp (SourceForge package)
Nonlinear Optimization in Matlab
There are various matlab functions available for nonlinear optimizations as well. These can be utilized in more exotic FIR filter designs, such as designs driven more by perceptual criteria:
 The fsolve function in Octave, or the Matlab Optimization Toolbox, attempts to solve unconstrained, overdetermined, nonlinear systems of equations.
 The Octave function sqp handles constrained nonlinear optimization.
 The Octave optim package includes many additional functions such as leasqr for performing LevenbergMarquardt nonlinear regression. (Say, e.g., which leasqr and explore its directory.)
 The Matlab Optimization Toolbox similarly contains many functions for optimization.
Next Section:
Spectrum of a Sinusoid
Previous Section:
MinimumPhase and Causal Cepstra