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 least-squares design (§4.3). Here we elaborate a bit on optimality formulations under various error norms.
- Sum of the absolute values of the elements
- ``City block'' distance
- ``Euclidean'' distance
- Minimized by ``Least Squares'' techniques
- 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 frequency-response error-magnitude 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 worst-case error (maximum weighted error-magnitude ripple) over all frequencies. This also means it is optimal in the sense because, as noted above, the norm of a weighted frequency-response error is the maximum magnitude over all frequencies:
Thus, we can say that an optimal Chebyshev filter minimizes the norm of the (possibly weighted) frequency-response 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 Parks-McClellan algorithm when applied to FIR filter design) [176,224,66]. This was illustrated in §4.6.4 above. The Parks-McClellan/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.13There 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  so that the desired frequency response can be taken to be real (i.e., first a zero-phase FIR filter is designed). The design of linear-phase 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 characteristic--that 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 . Another remarkable result is that the Remez multiple exchange converges monotonically to the unique optimal Chebyshev solution (in the absence of numerical round-off errors).
The window method (§4.5) and Remez-exchange method together span many practical FIR filter design needs, from ``quick and dirty'' to essentially ideal FIR filters (in terms of conventional specifications).
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
Enforcing even symmetry in the impulse response, i.e., , gives a zero-phase FIR filter that we can later right-shift samples to make a causal, linear phase filter. In this case, the frequency response reduces to a sum of cosines:
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 left-hand-side includes the alternating error, and the frequency grid iteratively seeks the frequencies of maximum error--the so-called extremal frequencies.
where these quantities are defined in (4.35). We can denote the optimal least-squares solution by
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.14Assuming all quantities are real, equating the gradient to zero yields the so-called normal equations
is known as the (Moore-Penrose) pseudo-inverse of the matrix . It can be interpreted as an orthogonal projection matrix, projecting onto the column-space of , as we illustrate further in the next section.
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 practice, the least-squares solution can be found by minimizing the sum of squared errors:
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 :
This is how the orthogonality principle can be used to derive the fact that the best least squares solution is given by
In matlab, it is numerically superior to use ``h= A h'' as opposed to explicitly computing the pseudo-inverse as in ``h = pinv(A) * d''. For a discussion of numerical issues in matrix least-squares problems, see, e.g., .
Some of the available functions are as follows:
- firls - least-squares linear-phase FIR filter design
for piecewise constant desired amplitude responses -- also designs
Hilbert transformers and differentiators
- fircls - constrained least-squares linear-phase FIR
filter design for piecewise constant desired amplitude responses --
constraints provide lower and upper bounds on the frequency response
- fircls1 - constrained least-squares linear-phase FIR
filter design for lowpass and highpass filters -- supports relative
weighting of pass-band and stop-band 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
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
where denotes the th row of the matrix . This can be expressed as
Introducing a new variable
then we can write
and our optimization problem can be written in more standard form:
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.
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 ,
- inverse filters (``deconvolution''),
- fractional delay filters ,
- interpolation polyphase subfilters (Chapter 11)
Nonlinear-Phase FIR Filter Design
Above, we considered only linear-phase (symmetric) FIR filters. The same methods also work for antisymmetric FIR filters having a purely imaginary frequency response, when zero-centered, such as differentiators and Hilbert transformers .
We now look at extension to nonlinear-phase FIR filters, managed by treating the real and imaginary parts separately in the frequency domain . In the nonlinear-phase 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 .
where , , and . Hence we have,
which can be written as
which is written in terms of only real variables.
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 Parks-McClellan''). Convergence is theoretically guaranteed for arbitrary magnitude and phase specifications versus frequency. It reduces to Parks-McClellan algorithm (Remez second algorithm) as a special case.
Finally, the fircband function in the Matlab DSP System Toolbox designs a variety of real FIR filters with various filter-types 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.
In Second-Order Cone Problems (SOCP), a linear function is minimized over the intersection of an affine set and the product of second-order (quadratic) cones [153,22]. Nonlinear, convex problem including linear and (convex) quadratic programs are special cases. SOCP problems are solved by efficient primal-dual interior-point 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.
- LIPSOL: Matlab code for linear programming using interior point methods.
- Matlab's linprog (in the Optimization Toolbox)
- Octave's lp (SourceForge package)
- 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 Levenberg-Marquardt nonlinear regression. (Say, e.g., which leasqr and explore its directory.)
- The Matlab Optimization Toolbox similarly contains many functions for optimization.
Spectrum of a Sinusoid
Minimum-Phase and Causal Cepstra