Source: https://www.dsprelated.com/freebooks/sasp/FIR_Digital_Filter_Design.html
Timestamp: 2019-04-20 03:02:35+00:00

Document:
As usual with the DTFT, the sampling rate is assumed to be . In practice, the DTFT is used in sampled form, , replacing the DTFT with a (zero-padded) FFT, as will be discussed in Chapter 8. To make the most of FFT processors for FIR filter implementation, we will need flexible ways to design many kinds of FIR filters.
This chapter provides a starting point in the area of FIR digital filter design. The so-called ``window method'' for FIR filter design, also based on the convolution theorem for Fourier transforms, is discussed in some detail, and compared with an optimal Chebyshev method. Other methods, such as least-squares, are discussed briefly to provide further perspective. Tools for FIR filter design in both Octave and the Matlab Signal Processing Toolbox are listed where applicable. For more information on digital filter design, see, e.g., the documentation for the Matlab Signal Processing Toolbox and/or [263,283,32,204,275,224,198,258].
Consider the ideal lowpass filter, depicted in Fig.4.1.
Figure 4.1: Desired amplitude response (gain versus frequency) for an ideal lowpass filter.
where denotes the normalized cut-off frequency in radians per sample. Thus, the impulse response of an ideal lowpass filter is a sinc function.
Unfortunately, we cannot implement the ideal lowpass filter in practice because its impulse response is infinitely long in time. It is also noncausal; it cannot be shifted to make it causal because the impulse response extends all the way to time . It is clear we will have to accept some sort of compromise in the design of any practical lowpass filter.
The subject of digital filter design is generally concerned with finding an optimal approximation to the desired frequency response by minimizing some norm of a prescribed error criterion with respect to a set of practical filter coefficients, perhaps subject also to some constraints (usually linear equality or inequality constraints) on the filter coefficients, as we saw for optimal window design in §3.13.5.3In audio applications, optimality is difficult to define precisely because perception is involved. It is therefore valuable to consider also suboptimal methods that are ``close enough'' to optimal, and which may have other advantages such as extreme simplicity and/or speed. We will examine some specific cases below.
Perhaps the most commonly employed error criterion in signal processing is the least-squares error criterion.
is also minimized by matching the leading terms of the desired impulse response.
Therefore, is an optimal least-squares approximation to when is given by (4.5). In other words, the frequency response of the filter is optimal in the (least-squares) sense.
Figure 4.3 shows the amplitude response of a length optimal least-squares FIR lowpass filter, for the case in which the cut-off frequency is one-fourth the sampling rate ( ).
Figure: Amplitude response of a length FIR lowpass-filter obtained by truncating the ideal impulse response.
We see that, although the impulse response is optimal in the least-squares sense (in fact optimal under any norm with any error-weighting), the filter is quite poor from an audio perspective. In particular, the stop-band gain, in which zero is desired, is only about 10 dB down. Furthermore, increasing the length of the filter does not help, as evidenced by the length 71 result in Fig.4.4.
Figure: Amplitude response of a length FIR lowpass-filter obtained by the Remez Exchange Algorithm (function firpm in the Matlab Signal Processing Toolbox).
We see that the Chebyshev design has a stop-band attenuation better than 60 dB, no corner-frequency resonance, and the error is equiripple in both stop-band (visible) and pass-band (not visible). Note also that there is a transition band between the pass-band and stop-band (specified in the call to firpm as being between normalized frequencies 0.5 and 0.6).
The main problem with the least-squares design examples above is the absence of a transition band specification. That is, the filter specification calls for an infinite roll-off rate from the pass-band to the stop-band, and this cannot be accomplished by any FIR filter. (Review Fig.4.2 for an illustration of more practical lowpass-filter design specifications.) With a transition band and a weighting function, least-squares FIR filter design can perform very well in practice. As a rule of thumb, the transition bandwidth should be at least , where is the FIR filter length in samples. (Recall that the main-lobe width of a length rectangular window is (§3.1.2).) Such a rule respects the basic Fourier duality of length in the time domain and ``minimum feature width'' in the frequency domain.
The window method for digital filter design is fast, convenient, and robust, but generally suboptimal. It is easily understood in terms of the convolution theorem for Fourier transforms, making it instructive to study after the Fourier theorems and windows for spectrum analysis. It can be effectively combined with the frequency sampling method, as we will see in §4.6 below.
where is the total normalized bandwidth of the lowpass filter in Hz (counting both negative and positive frequencies), and denotes the cut-off frequency in Hz. As noted earlier, we cannot implement this filter in practice because it is noncausal and infinitely long.
Since sinc decays away from time 0 as , we would expect to be able to truncate it to the interval , for some sufficiently large , and obtain a pretty good FIR filter which approximates the ideal filter. This would be an example of using the window method with the rectangular window. We saw in §4.3 that such a choice is optimal in the least-squares sense, but it designs relatively poor audio filters. Choosing other windows corresponds to tapering the ideal impulse response to zero instead of truncating it. Tapering better preserves the shape of the desired frequency response, as we will see. By choosing the window carefully, we can manage various trade-offs so as to maximize the filter-design quality in a given application.
Window functions are always time limited. This means there is always a finite integer such that for all . The final windowed impulse response is thus always time-limited, as needed for practical implementation. The window method always designs a finite-impulse-response (FIR) digital filter (as opposed to an infinite-impulse-response (IIR) digital filter).
The pass-band gain is primarily the area under the main lobe of the window transform, provided the main lobe ``fits'' inside the pass-band (i.e., the total lowpass bandwidth is greater than or equal to the main-lobe width of ).
The stop-band gain is given by an integral over a portion of the side lobes of the window transform. Since side-lobes oscillate about zero, a finite integral over them is normally much smaller than the side-lobes themselves, due to adjacent side-lobe cancellation under the integral.
The best stop-band performance occurs when the cut-off frequency is set so that the stop-band side-lobe integral traverses a whole number of side lobes.
The transition bandwidth is equal to the bandwidth of the main lobe of the window transform, again provided that the main lobe ``fits'' inside the pass-band.
For very small lowpass bandwidths , approaches an impulse in the frequency domain. Since the impulse is the identity operator under convolution, the resulting lowpass filter approaches the window transform for small . In particular, the stop-band gain approaches the window side-lobe level, and the transition width approaches half the main-lobe width. Thus, for good results, the lowpass cut-off frequency should be set no lower than half the window's main-lobe width.
fir1 designs lowpass, highpass, bandpass, and multi-bandpass filters.
fir2 takes an arbitrary magnitude frequency response specification.
The default window type is Hamming, but any window can be passed in as an argument. In addition, there is a function kaiserord for estimating the parameters of a Kaiser window which will achieve the desired filter specifications.
The matlab code below designs a bandpass filter which passes frequencies between 4 kHz and 6 kHz, allowing transition bands from 3-4 kHz and 6-8 kHz (i.e., the stop-bands are 0-3 kHz and 8-10 kHz, when the sampling rate is 20 kHz). The desired stop-band attenuation is 80 dB, and the pass-band ripple is required to be no greater than 0.1 dB. For these specifications, the function kaiserord returns a beta value of and a window length of . These values are passed to the function kaiser which computes the window function itself. The ideal bandpass-filter impulse response is computed in fir1, and the supplied Kaiser window is applied to shorten it to length .
Note the conciseness of the matlab code thanks to the use of kaiserord and fir1 from Octave or the Matlab Signal Processing Toolbox.
Figure 4.6 shows the magnitude frequency response of the resulting FIR filter . Note that the upper pass-band edge has been moved to 6500 Hz instead of 6000 Hz, and the stop-band begins at 7500 Hz instead of 8000 Hz as requested. While this may look like a bug at first, it's actually a perfectly fine solution. As discussed earlier (§4.5), all transition-widths in filters designed by the window method must equal the window-transform's main-lobe width. Therefore, the only way to achieve specs when there are multiple transition regions specified is to set the main-lobe width to the minimum transition width. For the others, it makes sense to center the transition within the requested transition region.
Figure 4.6: Amplitude response of the FIR bandpass filter designed by the window method.
where is the desired stop-band attenuation in dB (typical values in audio work are to ). Note that this estimate for becomes too small when the filter pass-band width approaches zero. In the limit of a zero-width pass-band, the frequency response becomes that of the Kaiser window transform itself. A non-zero pass-band width acts as a ``moving average'' lowpass filter on the side-lobes of the window transform, which brings them down in level. The kaiserord estimate assumes some of this side-lobe smoothing is present.
where now is the desired side-lobe attenuation in dB (as opposed to stop-band attenuation). A plot showing Kaiser window side-lobe level for various values of is given in Fig.3.28.
where is in radians between 0 and .
Without the function fir1, we would have to manually implement the window method of filter design by (1) constructing the impulse response of the ideal bandpass filter (a cosine modulated sinc function), (2) computing the Kaiser window using the estimated length and from above, then finally (3) windowing the ideal impulse response with the Kaiser window to obtain the FIR filter coefficients . A manual design of this nature will be illustrated in the Hilbert transform example of §4.6.
To provide some perspective on the results, let's compare the window method to the optimal Chebyshev FIR filter (§4.10) for the same length and design specifications above.
The following Matlab code illustrates two different bandpass filter designs. The first (different transition bands) illustrates a problem we'll look at. The second (equal transition bands, commented out), avoids the problem.
Figure 4.7 shows the frequency response of the Chebyshev FIR filter designed by firpm, to be compared with the window-method FIR filter in Fig.4.6. Note that the upper transition band ``blows up''. This is a well known failure mode in FIR filter design using the Remez exchange algorithm [176,224]. It can be eliminated by narrowing the transition band, as shown in Fig.4.8. There is no error penalty in the transition region, so it is necessary that each one be ``sufficiently narrow'' to avoid this phenomenon.
Figure 4.7: Amplitude response of the optimal Chebyshev FIR bandpass filter designed by the Remez exchange method.
Figure 4.8: Amplitude response of the optimal Chebyshev FIR bandpass filter as in Fig.4.7 with the upper transition band narrowed from 2 kHz down to 1 kHz in width.
We will now use the window method to design a complex bandpass filter which passes positive frequencies and rejects negative frequencies.
Since every real signal possesses a Hermitian spectrum , i.e., , it follows that, if we filter out the negative frequencies, we will destroy this spectral symmetry, and the output signal will be complex for every nonzero real input signal (excluding dc and half the sampling rate). In other terms, we want a filter which produces a ``single sideband'' (SSB) output signal in response to any real input signal. The Hermitian spectrum of a real signal can be viewed as two sidebands about dc (with one sideband being the ``conjugate flip'' of the other). See §2.3.3 for a review of Fourier symmetry-properties for real signals.
An ``analytic signal'' in signal processing is defined as any signal having only positive or only negative frequencies, but not both (typically only positive frequencies). In principle, the imaginary part of an analytic signal is computed from its real part by the Hilbert transform (defined and discussed below). In other words, one can ``filter out'' the negative-frequency components of a signal by taking its Hilbert transform and forming the analytic signal . Thus, an alternative problem specification is to ask for a (real) filter which approximates the Hilbert transform as closely as possible for a given filter order.
where denotes the frequency response of the Hilbert transform . Since by definition we have for , we must have for , so that for negative frequencies (an allpass response with phase-shift degrees). To pass the positive-frequency components unchanged, we would most naturally define for . However, conventionally, the positive-frequency Hilbert-transform frequency response is defined more symmetrically as for , which gives and , i.e., the positive-frequency components of are multiplied by .
Note that the point at can be defined arbitrarily since the inverse-Fourier transform integral is not affected by a single finite point (being a ``set of measure zero'').
The ideal filter impulse response is obtained by finding the inverse Fourier transform of (4.16). For discrete time, we may take the inverse DTFT of (4.16) to obtain the ideal discrete-time Hilbert-transform impulse response, as pursued in Problem 10. We will work with the usual continuous-time limit in the next section.
Thus, the Hilbert transform is a non-causal linear time-invariant filter.
Thus, the negative-frequency components of are canceled, while the positive-frequency components are doubled. This occurs because, as discussed above, the Hilbert transform is an allpass filter that provides a degree phase shift at all negative frequencies, and a degree phase shift at all positive frequencies, as indicated in (4.16). The use of the Hilbert transform to create an analytic signal from a real signal is one of its main applications. However, as the preceding sections make clear, a Hilbert transform in practice is far from ideal because it must be made finite-duration in some way.
Convolving a real signal with this kernel produces the imaginary part of the corresponding analytic signal. The way the ``window method'' for digital filter design is classically done is to simply sample the ideal impulse response to obtain and then window it to give . However, we know from above (e.g., §4.5.2) that we need to provide transition bands in order to obtain a reasonable design. A single-sideband filter needs a transition band between dc and , or higher, where denotes the main-lobe width (in rad/sample) of the window we choose, and a second transition band is needed between and .
Note that we cannot allow a time-domain sample at time 0 in (4.22) because it would be infinity. Instead, time 0 should be taken to lie between two samples, thereby introducing a small non-integer advance or delay. We'll choose a half-sample delay. As a result, we'll need to delay the real-part filter by half a sample as well when we make a complete single-sideband filter.
The matlab below illustrates the design of an FIR Hilbert-transform filter by the window method using a Kaiser window. For a more practical illustration, the sampling-rate assumed is set to Hz instead of being normalized to 1 as usual. The Kaiser-window parameter is set to , which normally gives ``pretty good'' audio performance (cf. Fig.3.28). From Fig.3.28, we see that we can expect a stop-band attenuation better than dB. The choice of , in setting the time-bandwidth product of the Kaiser window, determines both the stop-band rejection and the transition bandwidths required by our FIR frequency response.
Recall that, for a rectangular window, our minimum transition bandwidth would be Hz, and for a Hamming window, Hz. In this example, using a Kaiser window with ( ), the main-lobe width is on the order of Hz, so we expect transition bandwidths of this width. The choice above should therefore be sufficient, but not ``tight''.5.8 For each doubling of the filter length (or each halving of the sampling rate), we may cut in half.
w = kaiser(M,beta)'; % Kaiser window in "linear phase form"
The spectrum of this window (zero-padded by more than a factor of 8) is shown in Fig.4.9 (full magnitude spectrum) and Fig.4.10 (zoom-in on the main lobe).
Figure 4.9: The Kaiser window magnitude spectrum.
Figure 4.10: Kaiser window magnitude spectrum near dc.
The next step is to apply our Kaiser window to the ``desired'' impulse response, where ``desired'' means a time-shifted (by 1/2 sample) and bandlimited (to introduce transition bands) version of the ``ideal'' impulse response in (4.22). In principle, we are using the frequency-sampling method (§4.4) to prepare a desired FIR filter of length as the inverse FFT of a desired frequency response prepared by direct Fourier intuition. This long FIR filter is then ``windowed'' down to length to give us our final FIR filter designed by the window method.
If the smallest transition bandwidth is Hz, then the FFT size should satisfy . Otherwise, there may be too much time aliasing in the desired impulse response.5.10 The only non-obvious part in the matlab below is ``.^8'' which smooths the taper to zero and looks better on a log magnitude scale. It would also make sense to do a linear taper on a dB scale which corresponds to an exponential taper to zero.
Figure 4.11 shows our desired amplitude response so constructed.
Figure 4.11: Desired single-sideband-filter frequency response.
The real part of the desired impulse response is shown in Fig.4.12, and the imaginary part in Fig.4.13.
Figure 4.12: Real part of the desired single-sideband-filter impulse response.
Figure 4.13: Desired Hilbert-transform-filter impulse response.
Figure 4.14 and Fig.4.15 show the normalized dB magnitude frequency response of our final FIR filter consisting of the nonzero samples of hw.
Figure 4.14: Frequency response of the Kaiser-windowed impulse response.
Figure 4.15: Zoomed frequency response of the Kaiser-windowed impulse response.
In general, the end result will be a somewhat smoothed version of what we originally asked for in the frequency domain. However, this smoothing will be minimal if we asked for a truly ``doable'' desired frequency response. Because the above steps are fast and non-iterative, they can be used effectively in response to an interactive user interface (such as a set of audio graphic-equalizer sliders that are interpolated to form the desired frequency response), or even in a real-time adaptive system.
Let's now compare the window-method design using the Kaiser window to the optimal equiripple FIR filter design given by the Remez multiple exchange algorithm.
Note, by the way, that many filter-design software functions, such as firpm have special modes for designing Hilbert-transform filters .
It turns out that the Remez exchange algorithm has convergence problems for filters larger than a few hundred taps. Therefore, the FIR filter length was chosen above to be small enough to work out in this comparison. However, keep in mind that for very large filter orders, the Remez exchange method may not be an option. There are more recently developed methods for optimal Chebyshev FIR filter design, using ``convex optimization'' techniques, that may continue to work at very high orders [218,22,153]. The fast nonparametric methods discussed above (frequency sampling, window method) will work fine at extremely high orders.
Figure 4.16: Frequency response of the optimal Chebyshev FIR filter designed by the Remez exchange algorithm.
Figure 4.17: Transition region of the optimal Chebyshev frequency response.
In this case we did not normalize the peak amplitude response to 0 dB because it has a ripple peak of about 1 dB in the pass-band. Figure 4.18 shows a zoom-in on the pass-band ripple.
Figure 4.18: Pass-Band ripple for optimal Chebyshev frequency response.
The pass-band ripple is much smaller than 0.1 dB, which is ``over designed'' and therefore wasting of taps.
The stop-band response ``droops'' which ``wastes'' filter taps when stop-band attenuation is the only stop-band specification. In other words, the first stop-band ripple drives the spec ( dB), while all higher-frequency ripples are over-designed. On the other hand, a high-frequency ``roll-off'' of this nature is quite natural in the frequency domain, and it corresponds to a ``smoother pulse'' in the time domain. Sometimes making the stop-band attenuation uniform will cause small impulses at the beginning and end of the impulse response in the time domain. (The pass-band and stop-band ripple can ``add up'' under the inverse Fourier transform integral.) Recall this impulsive endpoint phenomenon for the Chebyshev window shown in Fig.3.33.
The pass-band is degraded by early roll-off. The pass-band edge is not exactly in the desired place.
The filter length can be thousands of taps long without running into numerical failure. Filters this long are actually needed for sampling rate conversion [270,218].
The stop-band is ideal, equiripple.
The pass-band is ideal, though over-designed for static audio spectra.
The computational design time is orders of magnitude larger than that for window method.
Practical application of Hilbert-transform filters to phase-vocoder analysis for additive synthesis is discussed in §G.10.
If the desired frequency response is real (as in simple band gains), either plan for a zero-phase filter in the end, or synthesize a desired phase, such as linear phase or minimum phase (see §4.8 below).
Perform the inverse Fourier transform of the (sampled) desired frequency response to obtain the desired impulse response.
Plot an overlay of the desired impulse response and the window to be applied, ensuring that the great majority of the signal energy in the desired impulse response lies under the window to be used.
Take an FFT (now with zero padding introduced by the window).
Plot an overlay of the original desired response and the response retained after time-domain windowing, and verify that the specifications are within an acceptable range.
In summary, FIR filters can be designed nonparametrically, directly in the frequency domain, followed by a final smoothing (windowing in the time domain) which guarantees that the FIR length will be precisely limited. As we'll discuss in Chapter 8, it is necessary to precisely limit the FIR filter length to avoid time-aliasing in an FFT-convolution implementation.
Above, we used the Hilbert transform to find the imaginary part of an analytic signal from its real part. A closely related application of the Hilbert transform is constructing a minimum phase  frequency response from an amplitude response.
If is computed from by the Hilbert transform, then is an ``analytic signal'' in the frequency domain. Therefore, it has no ``negative times,'' i.e., it is causal. The time domain signal corresponding to a log spectrum is called the cepstrum . It is reviewed in the next section that a frequency response is minimum phase if and only if the corresponding cepstrum is causal [198, Ch. 10], [263, Ch. 11].
To show that a frequency response is minimum phase if and only if the corresponding cepstrum is causal, we may take the log of the corresponding transfer function, obtaining a sum of terms of the form for the zeros and for the poles. Since all poles and zeros of a minimum-phase system must be inside the unit circle of the plane, the Laurent expansion of all such terms (the cepstrum) must be causal. In practice, as discussed in , we may compute an approximate cepstrum as an inverse FFT of the log spectrum, and make it causal by ``flipping'' the negative-time cepstral coefficients around to positive time (adding them to the positive-time coefficients). That is , for and for . This effectively inverts all unstable poles and all non-minimum-phase zeros with respect to the unit circle. In other terms, (if unstable), and (if non-minimum phase).
In digital signal processing, a Laurent series is typically expanded about points on the unit circle in the plane, because the unit circle--our frequency axis--must lie within the annulus of convergence of the series expansion in most applications. The power-of- terms are the noncausal terms, while the power-of- terms are considered causal. The term in the above general example is associated with time 0, and is included with the causal terms.
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.
In the limit as , the norm of is dominated by the maximum element of . Optimal Chebyshev filters minimize this norm of the frequency-response error.
An especially valuable property of FIR filter design under norms is that the error norm is typically a convex function of the filter coefficients, rendering it amenable to a wide variety of convex-optimization algorithms . The following sections look at some specific cases.
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).
Fine online introductions to the theory and practice of Chebyshev-optimal FIR filter design are given in [32,283].
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).
Another versatile, effective, and often-used 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 .
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.
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 the case of least-squares approximation, we are minimizing the Euclidean distance, which suggests the geometrical interpretation shown in Fig.4.19.
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., .
We will return to least-squares optimality in §5.7.1 for the purpose of estimating the parameters of sinusoidal peaks in spectra.
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.
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.
In general, these all involve complex desired frequency-response samples ( complex), and the matrix remains a more general Fourier transform matrix that cannot be reduced to a matrix of cosines.
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 .
which is written in terms of only real variables.
In summary, we can use the standard least-squares solvers in matlab and end up with a real solution for the case of complex desired spectra and nonlinear-phase FIR filters.
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.
The firgr function (formerly gremez) in the Matlab Filter Design Toolbox performs ``generalized'' FIR filter design, adding support for minimum-phase FIR filter design, among other features .
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.
See §3.13 for examples of optimal FFT window design using linprog.
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 Matlab Optimization Toolbox similarly contains many functions for optimization.
Again, this is only a sampling of what is available.

References: §3
 §4
 §4
 §4
 §2
 §4
 §4
 §4
 §3
 §4
 §3
 §5
 §3
 §3