FIR Digital Filter Design
FIR filters are basic in spectral audio signal processing. In fact, the fastest way to implement long FIR filters in conventional CPUs^{5.1} is by means of FFT convolution. The convolution theorem for Fourier transforms (§2.3.5) states that the convolution of an input signal with a filter impulseresponse is given by the inverse DTFT of the product of the signal's spectrum times the filter's frequency response , i.e.,
(5.1) 
where
(5.2) 
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 (zeropadded) 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 socalled ``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 leastsquares, 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].
The Ideal Lowpass Filter
Consider the ideal lowpass filter, depicted in Fig.4.1. An ideal lowpass may be characterized by a gain of 1 for all frequencies below some cutoff frequency in Hz, and a gain of 0 for all higher frequencies.^{5.2} The impulse response of the ideal lowpass filter is easy to calculate:where denotes the normalized cutoff 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.3}In 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.
Lowpass Filter Design Specifications
Typical design parameters for a lowpass filter are shown in Fig.4.2. The design parameters are defined as follows: stopband ripple ( dB is common)
 passband ripple ( dB typical)
 stopband edge frequency
 passband edge frequency
 TW: transition width
 SBA: stopband attenuation
Ideal Lowpass Filter Revisited
The ideal lowpass filter of Fig.4.1 can now be described by the following specifications: The transition width TW is zero ( in Fig.4.2).
 The passband and stopband ripples are both zero
( in Fig.4.2, and ).
Optimal (but poor if unweighted)
LeastSquares
Impulse Response Design
Perhaps the most commonly employed error criterion in signal
processing is the leastsquares error criterion.
Let
denote some ideal filter impulse response, possibly
infinitely long, and let
denote the impulse response of a
length
causal FIR filter that we wish to design. The sum of
squared errors is given by
(5.4) 
where does not depend on . Note that . We can minimize the error by simply matching the first terms in the desired impulse response. That is, the optimal leastsquares FIR filter has the following ``tap'' coefficients:
The same solution works also for any norm (§4.10.1). That is, the error
(5.6) 
is also minimized by matching the leading terms of the desired impulse response. In the (leastsquares) case, we have, by the Fourier energy theorem (§2.3.8),
(5.7) 
Therefore, is an optimal leastsquares approximation to when is given by (4.5). In other words, the frequency response of the filter is optimal in the (leastsquares) sense.
Examples
Figure 4.3 shows the amplitude response of a length optimal leastsquares FIR lowpass filter, for the case in which the cutoff frequency is onefourth the sampling rate ( ).hh = firpm(L1,[0 0.5 0.6 1],[1 1 0 0]);where, in terms of the lowpass design specs defined in §4.2 above, we are asking for
 (passband edge frequency)^{5.5}
 (stopband edge frequency)
Frequency Sampling Method for
FIR Filter Design
The frequencysampling method for FIR filter design is perhaps
the simplest and most direct technique imaginable when a desired
frequency response has been specified. It consists simply of
uniformly sampling the desired frequency response, and
performing an inverse DFT to obtain the corresponding (finite)
impulse response [224, pp. 10523],
[198, pp. 25155]. The results are not optimal,
however, because the response generally deviates from what is desired
between the samples. When the desired frequencyresponse is
undersampled, which is typical, the resulting impulse response
will be time aliased to some extent. It is important to
evaluate the final impulse response via a simulated DTFT (FFT with
lots of zero padding), comparing to the originally desired frequency
response.
The frequencysampling method for FIR filter design is illustrated in
§4.6.2 below in the context of the window method for FIR
filter design, to which we now turn.
Window Method for FIR Filter Design
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. The window method consists of simply ``windowing'' a theoretically ideal filter impulse response by some suitably chosen window function , yielding(5.8) 
For example, as derived in Eq. (4.3), the impulse response of the ideal lowpass filter is the well known sinc function
sinc  (5.9) 
where is the total normalized bandwidth of the lowpass filter in Hz (counting both negative and positive frequencies), and denotes the cutoff 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 leastsquares 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 tradeoffs so as to maximize the filterdesign 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 timelimited, as needed for practical implementation. The window method always designs a finiteimpulseresponse (FIR) digital filter (as opposed to an infiniteimpulseresponse (IIR) digital filter). By the dual of the convolution theorem, pointwise multiplication in the time domain corresponds to convolution in the frequency domain. Thus, the designed filter has a frequency response given by
(5.10) 
where is the ideal frequency response and is the window transform. For the ideal lowpass filter, is a rectangular window in the frequency domain. The frequency response is thus obtained by convolving the rectangular window with the window transform . This implies several points which can be immediately seen in terms of this convolution operation:
 The passband gain is primarily the area under the main lobe of the window transform, provided the main lobe ``fits'' inside the passband (i.e., the total lowpass bandwidth is greater than or equal to the mainlobe width of ).
 The stopband gain is given by an integral over a portion of the side lobes of the window transform. Since sidelobes oscillate about zero, a finite integral over them is normally much smaller than the sidelobes themselves, due to adjacent sidelobe cancellation under the integral.
 The best stopband performance occurs when the cutoff frequency is set so that the stopband sidelobe 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 passband.
 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 stopband gain approaches the window sidelobe level, and the transition width approaches half the mainlobe width. Thus, for good results, the lowpass cutoff frequency should be set no lower than half the window's mainlobe width.
Matlab Support for the Window Method
Octave and the Matlab Signal Processing Toolbox have two functions implementing the window method for FIR digital filter design: fir1 designs lowpass, highpass, bandpass, and multibandpass filters.
 fir2 takes an arbitrary magnitude frequency response specification.
Bandpass Filter Design Example
The matlab code below designs a bandpass filter which passes frequencies between 4 kHz and 6 kHz, allowing transition bands from 34 kHz and 68 kHz (i.e., the stopbands are 03 kHz and 810 kHz, when the sampling rate is 20 kHz). The desired stopband attenuation is 80 dB, and the passband 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 bandpassfilter impulse response is computed in fir1, and the supplied Kaiser window is applied to shorten it to length .fs = 20000; % sampling rate F = [3000 4000 6000 8000]; % band limits A = [0 1 0]; % band type: 0='stop', 1='pass' dev = [0.0001 10^(0.1/20)1 0.0001]; % ripple/attenuation spec [M,Wn,beta,typ] = kaiserord(F,A,dev,fs); % window parameters b = fir1(M,Wn,typ,kaiser(M+1,beta),'noscale'); % filter designNote 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 passband edge has been moved to 6500 Hz instead of 6000 Hz, and the stopband 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 transitionwidths in filters designed by the window method must equal the windowtransform's mainlobe width. Therefore, the only way to achieve specs when there are multiple transition regions specified is to set the mainlobe width to the minimum transition width. For the others, it makes sense to center the transition within the requested transition region.
Under the Hood of kaiserord
Without kaiserord, we would need to implement Kaiser's formula [115,67] for estimating the Kaiserwindow required to achieve the given filter specs:where is the desired stopband attenuation in dB (typical values in audio work are to ). Note that this estimate for becomes too small when the filter passband width approaches zero. In the limit of a zerowidth passband, the frequency response becomes that of the Kaiser window transform itself. A nonzero passband width acts as a ``moving average'' lowpass filter on the sidelobes of the window transform, which brings them down in level. The kaiserord estimate assumes some of this sidelobe smoothing is present. A similar function from [198] for window design (as opposed to filter design^{5.7}) is
where now is the desired sidelobe attenuation in dB (as opposed to stopband attenuation). A plot showing Kaiser window sidelobe level for various values of is given in Fig.3.28. Similarly, the filter order is estimated from stopband attenuation and desired transition width using the empirical formula
(5.13) 
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.
Comparison to the Optimal Chebyshev FIR Bandpass Filter
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.M = 101; normF = [0 0.3 0.4 0.6 0.8 1.0]; % transition bands different %normF = [0 0.3 0.4 0.6 0.7 1.0]; % transition bands the same amp = [0 0 1 1 0 0]; % desired amplitude in each band [b2,err2] = firpm(M1,normF,amp); % optimal filter of length MFigure 4.7 shows the frequency response of the Chebyshev FIR filter designed by firpm, to be compared with the windowmethod 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. Remember the rule of thumb that the narrowest transitionband possible for a length FIR filter is on the order of , because that's the width of the mainlobe of a length rectangular window (measured between zerocrossings) (§3.1.2). Therefore, this value is quite exact for the transitionwidths of FIR bandpass filters designed by the window method using the rectangular window (when the mainlobe fits entirely within the adjacent passband and stopband). For a Hamming window, the windowmethod transition width would instead be . Thus, we might expect an optimal Chebyshev design to provide transition widths in the vicinity of , but probably not too close to or below In the example above, where the sampling rate was kHz, and the filter length was , we expect to be able to achieve transition bands circa Hz, but not so low as Hz. As we found above, Hz was underconstrained, while Hz was ok, being near the ``Hamming transition width.''
Hilbert Transform Design Example
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 symmetryproperties 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 negativefrequency 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.Primer on Hilbert Transform Theory
We need a Hilberttransform filter to compute the imaginary part of the analytic signal given its real part . That is,(5.14) 
where . In the frequency domain, we have
(5.15) 
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 phaseshift degrees). To pass the positivefrequency components unchanged, we would most naturally define for . However, conventionally, the positivefrequency Hilberttransform frequency response is defined more symmetrically as for , which gives and , i.e., the positivefrequency components of are multiplied by . In view of the foregoing, the frequency response of the ideal Hilberttransform filter may be defined as follows:
Note that the point at can be defined arbitrarily since the inverseFourier 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 discretetime Hilberttransform impulse response, as pursued in Problem 10. We will work with the usual continuoustime limit in the next section.
Hilbert Transform
The Hilbert transform of a real, continuoustime signal may be expressed as the convolution of with the Hilbert transform kernel:That is, the Hilbert transform of is given by
(5.18) 
Thus, the Hilbert transform is a noncausal linear timeinvariant filter. The complex analytic signal corresponding to the real signal is then given by
(5.19) 
To show this last equality (note the lower limit of 0 instead of the usual ), it is easiest to apply (4.16) in the frequency domain:
(5.20)  
(5.21) 
Thus, the negativefrequency components of are canceled, while the positivefrequency 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 finiteduration in some way.
Filtering and Windowing the Ideal
HilbertTransform Impulse Response
Let
denote the convolution kernel of the continuoustime
Hilbert transform from (4.17) above:
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 singlesideband filter needs a transition band between dc and , or higher, where denotes the mainlobe width (in rad/sample) of the window we choose, and a second transition band is needed between and . Note that we cannot allow a timedomain 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 noninteger advance or delay. We'll choose a halfsample delay. As a result, we'll need to delay the realpart filter by half a sample as well when we make a complete singlesideband filter. The matlab below illustrates the design of an FIR Hilberttransform filter by the window method using a Kaiser window. For a more practical illustration, the samplingrate assumed is set to Hz instead of being normalized to 1 as usual. The Kaiserwindow 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 stopband attenuation better than dB. The choice of , in setting the timebandwidth product of the Kaiser window, determines both the stopband rejection and the transition bandwidths required by our FIR frequency response.
M = 257; % window length = FIR filter length (Window Method) fs = 22050; % sampling rate assumed (Hz) f1 = 530; % lower passband limit = transition bandwidth (Hz) beta = 8; % beta for Kaiser window for decent sidelobe rejectionRecall 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 mainlobe 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.
Matlab, Continued
Given the above design parameters, we compute some derived parameters as follows:fn = fs/2; % Nyquist limit (Hz) f2 = fn  f1; % upper passband limit N = 2^(nextpow2(8*M)); % large FFT for interpolated display k1 = round(N*f1/fs); % lower band edge in bins if k1<2, k1=2; end; % cannot have dc or fn response kn = N/2 + 1; % bin index at Nyquist limit (1based) k2 = knk1+1; % highfrequency band edge f1 = k1*fs/N % quantized bandedge frequencies f2 = k2*fs/NSetting the upper transition band the same as the lowfrequency band ( ) provides an additional benefit: the symmetry of the desired response about cuts the computational expense of the filter in half, because it forces every other sample in the impulse response to be zero [224, p. 172].^{5.9}
Kaiser Window
With the filter length and Kaiser window as given above, we may compute the Kaiser window itself in matlab viaw = kaiser(M,beta)'; % Kaiser window in "linear phase form"The spectrum of this window (zeropadded by more than a factor of 8) is shown in Fig.4.9 (full magnitude spectrum) and Fig.4.10 (zoomin on the main lobe).
Windowing a Desired Impulse Response Computed by the
Frequency Sampling Method
The next step is to apply our Kaiser window to the ``desired'' impulse
response, where ``desired'' means a timeshifted (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
frequencysampling 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 nonobvious
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.
H = [ ([0:k12]/(k11)).^8,ones(1,k2k1+1),... ([k12:1:0]/(k11)).^8, zeros(1,N/21)];Figure 4.11 shows our desired amplitude response so constructed. Now we inverseFFT the desired frequency response to obtain the desired impulse response:
h = ifft(H); % desired impulse response hodd = imag(h(1:2:N)); % This should be zero ierr = norm(hodd)/norm(h); % Look at numerical roundoff error % Typical value: ierr = 4.1958e15 % Also look at time aliasing: aerr = norm(h(N/2N/32:N/2+N/32))/norm(h); % Typical value: 4.8300e04The real part of the desired impulse response is shown in Fig.4.12, and the imaginary part in Fig.4.13. Now use the Kaiser window to timelimit the desired impulse response:
% put window in zerophase form: wzp = [w((M+1)/2:M), zeros(1,NM), w(1:(M1)/2)]; hw = wzp .* h; % singlesideband FIR filter, zerocentered Hw = fft(hw); % for results display: plot(db(Hw)); hh = [hw(N(M1)/2+1:N),hw(1:(M+1)/2)]; % caual FIR % plot(db(fft([hh,zeros(1,NM)]))); % freq resp plotFigure 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.
More General FIR Filter Design
We have looked at more than just FIR filter design by the window method and frequencysampling technique. The general steps were Prepare a desired frequencyresponse that ``seems achievable''
 Inverse FFT
 Window the result (timelimit it)
 FFT that to see how it looks
Comparison to Optimal Chebyshev FIR Filter
Let's now compare the windowmethod 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 filterdesign software functions, such as
firpm have special modes for designing Hilberttransform
filters [224].
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.
The following Matlab command will try to design the FIR
Hilberttransform filter of the desired length with the desired
transition bands:
hri = firpm(M1, [f1,f2]/fn, [1,1], [1], 'Hilbert');Instead, however, we will use a more robust method [228] which uses the Remez exchange algorithm to design a lowpass filter, followed by modulation of the lowpass impulseresponse by a complex sinusoid at frequency in order to frequencyshift the lowpass to the singlesideband filter we seek:
tic; % remember the current time hrm = firpm(M1, [0,(f2fs/4)/fn,0.5,1], [1,1,0,0], [1,10]); dt = toc; % design time dt can be minutes hr = hrm .* j .^ [0:M1]; % modulate lowpass to singlesidebandThe weighting [1,10] in the call to firpm above says ``make the passband ripple times that of the stopband.'' For steadystate audio spectra, passband ripple can be as high as dB or more without audible consequences.^{5.11} The result is shown in Fig.4.16 (full amplitude response) and Fig.4.17 (zoomin on the dc transition band). By symmetry the highfrequency transition region is identical (but flipped):
Conclusions
We can note the following points regarding our singlesideband FIR filter design by means of direct Fourier intuition, frequencysampling, and the windowmethod: The passband ripple is much smaller than 0.1 dB, which is ``over designed'' and therefore wasting of taps.
 The stopband response ``droops'' which ``wastes'' filter taps when stopband attenuation is the only stopband specification. In other words, the first stopband ripple drives the spec ( dB), while all higherfrequency ripples are overdesigned. On the other hand, a highfrequency ``rolloff'' of this nature is quite natural in the frequency domain, and it corresponds to a ``smoother pulse'' in the time domain. Sometimes making the stopband attenuation uniform will cause small impulses at the beginning and end of the impulse response in the time domain. (The passband and stopband 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 passband is degraded by early rolloff. The passband 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 stopband is ideal, equiripple.
 The transition bandwidth is close to half that of the window method. (We already knew our chosen transition bandwidth was not ``tight'', but our ruleofthumb based on the Kaiserwindow mainlobe width predicted only about % excess width.)
 The passband is ideal, though overdesigned for static audio spectra.
 The computational design time is orders of magnitude larger than that for window method.
 The design fails to converge for filters much longer than 256 taps. (Need to increase working precision or use a different method to get longer optimal Chebyshev FIR filters.)
Generalized Window Method
Reiterating and expanding on points made in §4.6.3, often we need a filter with a frequency response that is not analytically known. An example is a graphic equalizer in which a user may manipulate sliders in a graphical user interface to control the gain in each of several frequency bands. From the foregoing, the following procedure, based in spirit on the window method (§4.5), can yield good results: Synthesize the desired frequency response as the smoothest possible interpolation of the desired frequencyresponse points. For example, in a graphic equalizer, cubic splines [286] could be used to connect the desired band gains.^{5.12}
 If the desired frequency response is real (as in simple band gains), either plan for a zerophase 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.
 Multiply by the window.
 Take an FFT (now with zero padding introduced by the window).
 Plot an overlay of the original desired response and the response retained after timedomain windowing, and verify that the specifications are within an acceptable range.
MinimumPhase Filter Design
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 [263] frequency response from an amplitude response. Let denote a desired complex, minimumphase frequency response in the digital domain ( plane):(5.23) 
and suppose we have only the amplitude response
(5.24) 
Then the phase response can be computed as the Hilbert transform of . This can be seen by inspecting the log frequency response:
(5.25) 
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 [263]. 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].
MinimumPhase and Causal Cepstra
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 minimumphase 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 [263], we may compute an approximate cepstrum as an inverse FFT of the log spectrum, and make it causal by ``flipping'' the negativetime cepstral coefficients around to positive time (adding them to the positivetime coefficients). That is , for and for . This effectively inverts all unstable poles and all nonminimumphase zeros with respect to the unit circle. In other terms, (if unstable), and (if nonminimum phase). The Laurent expansion of a differentiable function of a complex variable can be thought of as a twosided Taylor expansion, i.e., it includes both positive and negative powers of , e.g.,(5.26) 
In digital signal processing, a Laurent series is typically expanded about points on the unit circle in the plane, because the unit circleour frequency axismust lie within the annulus of convergence of the series expansion in most applications. The powerof terms are the noncausal terms, while the powerof terms are considered causal. The term in the above general example is associated with time 0, and is included with the causal terms.
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
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 Analysis of Sinusoids
Previous Section:
Spectrum Analysis Windows