Because the integral of the square of a signal is a measure of its energy, there is some physical reason for minimizing the integral of the squared error , . Also, because of Parseval's theorem, a least squares approximation in the frequency domain is a least squares approximation in the time domain. However, minimizing the worst case squared error induces a minimum Chebyshev error problem in some formulations .

Discrete Frequency Samples of Error If we approximate the integral squared error by the sum of the squared error as given by q = 1 L ∑ k = 0 L - 1 ( A ( ω k ) - A d ( ω k ) ) 2 = 1 L ∑ k = 0 L - 1 e ( ω k ) 2 ≈ ∫ 0 π ( A ( ω ) - A d ( ω ) ) 2 d ω = ∫ 0 π e ( ω ) 2 d ω where the approximation error as a function of frequency is defined by e(ω)=A(ω)-Ad(ω) with A(ω) being the amplitude response of the filter and Ad(ω) being the desired amplitude response or the ideal response. The matrix statement for the error vector becomes ϵ = A - A d = C a - A d where C is the matrix of cosines from (), a is the vector of half of the filter coefficients from (), and Ad is the vector of samples of the ideal desired amplitude response. The number of samples of the amplitude response is L which should be five to twenty times the length of the filter to give a good approximation of the integral in most cases. The error to be minimized is q = ϵ T ϵ except for a scale factor of 1L. This could also be posed for the general phase problem by using H(ωk) rather than A(ωk) and h(n), the actual impulse response, rather than a(n), a nomralized half of the impulse response.

Truncated frequency sampling design using the inverse FFT or IDCT The design problem is posed by defining an error measure q as a sum of the squared differences between the actual and the desired frequency response over a set of L frequency samples. This error function is defined as q = 1 L ∑ k = 0 L - 1 | H ( ω k ) - H d ( ω k ) | 2 where Hd(ωk) are the L samples of the desired response. This problem is easier to formulate and solve if the frequency samples are equally spaced as in () which gives ω k = 2 π k / L and the problem is restricted to linear-phase filters where the real-valued amplitude A(ω) can be approximated rather than the complex frequency response H(ω). For approximations to a complex response, see "Complex L 2 and Minimum Phase Approximation". Linear phase and equally spaced samples cause () to become q = 1 L ∑ k = 0 L - 1 | A ( 2 π k / L ) - A d ( 2 π k / L ) | 2 or with a simpler notation q = 1 L ∑ k = 0 L - 1 | A k - A d k | 2 A very powerful property of the Fourier transform allows a straightforward design of least-squared-error FIR filters. Parseval's Theorem, which is based on the orthogonality of the DFT, states that the error defined by () in the frequency domain can also be calculated in the time domain by q = ∑ n = 0 L - 1 | h ( n ) - h d ( n ) | 2 where hd(n) is the length-L symmetric FIR filter that has the L frequency response amplitude samples Adk. This may be calculated by the frequency sampling method of using the special formulas such as () for length L or the inverse DFT. The filter to be designed has a length-N symmetric impulse response h(n) with L frequency response samples Ak. Because the filter h(n) is of length-N and symmetric, the error equation () can be split into two sums q = ∑ n = - M M | h ^ ( n ) - h ^ d ( n ) | 2 + 2 ∑ n = M + 1 ( L - 1 ) / 2 | h ^ d ( n ) | 2 where h^(n) and h^d(n) are the inverse DTFTs of Ak and Adk respectively, which means they are the h(n) and hd(n) shifted to be symmetic about n=0. This requires the number of frequency samples L must be odd. Equation clearly shows that to minimize q, the N values of h(n) are chosen to be equal to the equivalent N values of hd(n) making the first sum equal zero. In other words, h(n) is obtained by symmetrically truncating hd(n). The residual error is then given by the second summation above. An examination of the residual error as a function of N may aid in the choice of the filter length N. For the Type 1 linear-phase FIR filter (described in ) which has an odd length N and an even-symmetric impulse response, the L equally spaced samples of the frequency response from () gives A k = ∑ n = 0 M - 1 2 h ( n ) cos ( 2 π ( M - n ) k / L ) + h ( M ) for k=0,1,2,....,L-1, where M=(N-1)/2. This formula was derived as a special case of the DFT applied to the Type 1 real, even-symmetric FIR filter coefficients to calculate the sampled amplitude of the frequency response (perhaps better posed using a(n)). It was noted in that it is also a cosine transform and it can be shown that this transformation is orthogonal over the independent values of Ak, just as the DFT is. The desired ideal amplitude gives the ideal impulse response hd(n) from (0 by h d ( n ) = 1 N [ A d 0 + ∑ k = 1 M - 1 2 A d k cos ( 2 π ( n - M ) k / N ) ] . for n=0,1,⋯,L-1. This is used in (), and is the ideal impulse response that is truncated and shifted to give a causal, symmetric h(n). Use of the alternative equally-spaced sampling in (), which has no sample at zero frequency, requires hd(n) be calculated from () and (). The Type 2 filters with even N are developed in a similar way and use the design formulas () and (). These methods are summarized by: The filter design procedure for an odd-length Type 1 filter is to first design an odd-length-L FIR filter by the frequency sampling method from () or () or the IDFT, then to symmetrically truncate it to the desired odd-length N and shift it to make h(n) causal. To design an even-length Type 2 filter , start with an even-length-L frequency-sampling design from () or () or the IDFT and symmetrically truncate and shift. The resulting length-N FIR filters are optimal LS-error approximations to the desired frequency response over the L frequency samples. This approach can also be applied to the general arbitrary phase FIR filter design problem.

Weighted, Unevenly Sampled Discrete Least Squared Error Filter Design by Solving Simultaneous Equations It is sometimes desirable to formulate the least squared error design problem using unequally-spaced frequency samples and/or a weighting function on the error. This is not possible using the IDFT or derived formulas above and requires a different approach to the solution. Samples of the amplitude response derived for N odd in () are given by A ( ω k ) = ∑ n = 1 M 2 h ( M - n ) cos ( ω k n ) + h ( M ) for k=0,1,⋯,L-1. This relates the L frequency samples A(ωk) to the M+1 independent values of the symmetric length-N impulse response h(n). In the design problem where the Ak are given and the values for h(n) are to be found, this represents L equations with M+1 unknowns. Because of the symmetries of A(ω) shown in , only half of the L values of Ak are independent; however, in some cases, to have proper weights on all L samples, all must be calculated. sampled at L arbitrary frequencies can be written as a matrix equation C a = A where a is an M+1 length vector with elements which are the first half of h(n). C is an L by (M+1) matrix of the cosine terms from (), and A is a length-L vector of the frequency samples A(ωk). If the formula for the calculation of L values of the frequency response of a length-N FIR filter in () is used to define an error vector of differences as defined in () and the result is written in the matrix formulation of (), the error becomes C a = A = A d + e or C a - A d = e where e is a vector of differences between the actual and desired samples of the frequency response. The error measure defined in () becomes the quadratic form q = e T e For L>N, equation () is over determined and cannot, in general, be solved for a. The filter design error measure is the norm of e, as given in (). This error measure is minimized by making e orthogonal to the columns of C in (). Multiplying both sides of ( by the transpose of C gives C T C a = C T A d + C T e In order for q to be minimum, e must be orthogonal to the columns of C and, therefore, CTe must be zero. Hence, the optimal a must satisfy the “normal equations" , , which are C T C a = C T A d and which can be rewritten in terms of the pseudo-inverse , as a = [ C T C ] - 1 C T A d If L=N, this becomes the regular frequency-sampling problem and can be solved with zero error. For the case of interest in this section, where L>N, there are still only M+1 equations to be solved. For L>>N, equation () may be ill-conditioned, and () should not be used to solve them. Special methods will be necessary to avoid serious numerical problems . If a weighted error function is desired, () is modified to give q = 1 L ∑ k = 0 L - 1 W k | A ( ω k ) - A d ( ω k ) | 2 The normal equations of () become C T W C a = C T W A d where W is a positive-definite matrix of the weights. If zero weights are desired, the effect is be achieved by removing those frequencies from the set of L frequencies, not by using a zero value weight which would violate the vector-space conditions of a well-posed minimization problem. Although developed here for the linear-phase filter, () is a very general design approach for the FIR filter that allows arbitrary phase, as well as uneven frequency sampling and a weighting function in the error definition. For the arbitrary phase case, a complex F is obtained from sampling () and the full h(n) is used. For the special case of the equally-spaced frequency samples and linear- phase filter with unity weighting, the solution of () or () is the same as given by the frequency sampling design formulas. One of the important uses of the unequally spaced frequency samples is to create a transition band between the pass and stopbands where there are no samples. This “don't care" band does not contribute to the error measure q and allows better approximation to occur over the pass and stopbands. Of the many ways to solve () or (), one of the easiest and most reliable is the use of Matlab , which has a special command to solve this least-mean-squared error problem. Equation () should not be solved directly. For large L, it is ill-conditioned and a direct solution will probably have large errors. Matlab uses special algorithms to minimize these numerical errors. This approach was applied to the same problems that were solved by frequency sampling in the previous section. For N=L, the same results are obtained, thus verifying the theoretical prediction. As L becomes larger compared to N, more control is exerted over the behavior between the original sample points. As L becomes large compared to N, the solution approaches the same results as obtained where the error is defined as a continuous function of frequency and the integral of the squared error is minimized. Although the solution of the normal equations is a powerful and flexible technique, it can be slow, have numerical problems, and require large amounts of computer memory.

Examples of Discrete Least Squared Error Filter Design Here we will give examples of several least squared error designs of FIR filters.

Frequency Response of Length-15 FIR Filter Designed by Least Squared Error As for the frequency sampling design, we see a good lowpass filter frequency response with the actual amplitude interpolating the desired values at different points from the frequency sampling example in even though the length and band edge are the same. Notice there is less over shoot but more ripple near f=0. The Gibbs phenomenon is the same as for the Fourier series. If a transition band is introduces in the ideal amplitude response between f=0.4 and f=0.6 with a straight line, the overshoot is reduced significantly but with a slightly slower transition from the pass to stop band. This is illustrated in .

Frequency Response of Length-15 FIR Filter with a Transition Band Designed by Least Squared Error

Continuous Frequency Definition of Error Because the energy of a signal is the integral of the sum of the squares of the Fourier transform magnitude and because specifications are usually given in the frequency domain, a very reasonable error measure to minimize is the integral squared error given by q = 1 π ∫ 0 π | A d ( ω ) - A ( ω ) | 2 d ω where Ad(ω) is the desired ideal amplitude response, A(ω)=∑na(n)cos(ω(M-N)n) is the achieved amplitude response with the length h(n) related to h(n) by (). This integral squared error is approximated by the discrete squared error defined in () for L>>N which in some cases is much easier to minimize. However for some very useful cases, formulas can be found for h(n) that minimize () and that is what we will be considering in this section.

The Unweighted Least Integral Squared Error Approximation If the error measure is the unweighted integral squared error defined in (), Parseval's theorem gives the equivalent time-domain formulation for the error to be q = ∑ n = - ∞ ∞ | h d ( n ) - h ( n ) | 2 = 1 π ∫ 0 π | A d ( ω ) - A ( ω ) | 2 d ω . In general, this ideal response is infinite in duration and, therefore, cannot be realized exactly by an actual FIR filter. As was done in the case of the discrete error measure, we break the infinite sum in () into two parts, one of which depends on h(n) and the other does not. q = ∑ n = - M M | h d ( n ) - h ( n ) | 2 + 2 ∑ n = M + 1 ∞ | h d ( n ) | 2 Again, we see that the minimum q is achieved by using h(n)=hd(n) for -M≤n≤M. In other words, the infinitely long hd(n) is symmetrically truncated to give the optimal least integral squared error approximation. The problem then becomes one of finding the hd(n) to truncate. Here the integral definition of approximation error is used. This is usually what we really want, but in some cases the integrals can not be carried out and the sampled method above must be used.

Ideal Constant Gain Passband Lowpass Filter Here we assume the simplest ideal lowpass single band FIR filter to have unity passband gain for 0<ω<ω0 and zero stopband gain for ω0<ω<π similar to those in a and . This gives A d ( ω ) = { 0 ω 0 ω π 1 0 ω ω 0 as the ideal desired amplitude response. The ideal shifted filter coefficients are the inverse DTFT from () of this amplitude which for N odd are given by h ^ d ( n ) = 1 π ∫ 0 π A d ( ω ) cos ( ω n ) d ω = 1 π ∫ 0 ω 0 cos ( ω n ) d ω = ω 0 π sin ( ω 0 n ) ω 0 n which is sometimes called a “sinc" function. Note h^d(n) is generally infinite in length. This is now symmetrically truncated and shifted by M=(N-1)/2 to give the optimal, causal length-N FIR filter coefficients as h ( n ) = ω 0 π sin ( ω 0 ( n - M ) ) ω 0 ( n - M ) for 0 ≤ n ≤ N - 1 and h(n)=0 otherwise. The corresponding derivation for an even length starts with the inverse DTFT in () for a shifted even length filter is h ^ d = 1 π ∫ 0 π A d ( ω ) cos ( ω ( n + 1 / 2 ) ) d ω = ω 0 π sin ( ω 0 ( n + 1 / 2 ) ) ω 0 ( n + 1 / 2 ) which when truncated and shifted by N/2 gives the same formula as for the odd length design in () but one should note that M=(N-1)/2 is not an integer for an even N.

Ideal Linearly Increasing Gain Passband Lowpass Filter We now derive the design formula for a filter with an ideal amplitude response that is a linearly increasing function in the passband rather than a constant as was assumed above. This ideal amplitude response is given by and illustrated in For N odd, the ideal infinitely long shifted filter coefficients are the inverse DTFT of this amplitude given by A d ( ω ) = { 0 ω 0 ω π 1 π ω 0 ω ω 0 h ^ d ( n ) = 1 π ∫ 0 ω 0 ( ω π ) cos ( ω n ) d ω = cos ( ω 0 n ) - 1 π 2 n 2 + ω 0 sin ( ω 0 n ) π 2 n with the indeterminate h^d(0)=ω022π2. This is now truncated and shifted by M=(N-1)/2 to give the optimal, causal length-N FIR filter coefficients as h ( n ) = cos ( ω 0 ( n - M ) ) - 1 π 2 ( n - M ) 2 + ω 0 sin ( ω 0 ( n - M ) ) π 2 ( n - M ) for 0 ≤ n ≤ N - 1 and h(n)=0 otherwise. The corresponding derivation for an even length starts with the inverse DTFT for a shifted even length filter in () and after shifting by N/2 gives the same result as ().

Ideal Frequency Response of an FIR Filter with Increasing Gain in the Passband and Lowpass Cutoff

Ideal Differentiator plus Lowpass Filter Fortunately the inverse DTFT for an ideal differentiator combined with a lowpass filter can also be analytically evaluated. The ideal amplitude response is the same as () and but, since this case has an odd symmetric impulse response, the inverse DTFT uses sine functions which for odd N gives h ^ d ( n ) = 1 π ∫ 0 ω 0 ( 1 π ω ) sin ( ω n ) d ω = sin ( ω 0 n ) π 2 n 2 - ω 0 cos ( ω 0 n ) π 2 n with the indeterminate h^d(0)=0. This is now truncated and shifted by M=(N-1)/2 to give the optimal, causal length-N FIR filter coefficients as h ( n ) = sin ( ω 0 ( n - M ) ) π 2 ( n - M ) 2 - ω 0 cos ( ω 0 ( n - M ) ) π 2 ( n - M ) for 0 ≤ n ≤ N - 1 and h(n)=0 otherwise. Again the corresponding derivation for an even length gives the same result as in (). Note this very general single formula includes as special cases the odd and even length full band (ω0=π) differentiator given in . Also note that for a full band differentiator, an even length is much preferred because of the zero at ω=π for an odd length. However, for the differentiator with a lowpass filter, the zero aids in the lowpass filtering and, therefore, might be an advantage.

Hilbert Transformer The inverse DTFT for an ideal Hilbert transform combined with a lowpass filter can also be analytically evaluated. The ideal amplitude response is the same as () but with a constant phase shift of ϕ=π/2. Since this case has an odd symmetric impulse response, the inverse DTFT uses sine functions which for odd N which gives h ^ d ( n ) = 1 π ∫ 0 ω 0 sin ( ω n ) d ω = 1 - cos ( ω 0 n ) π n with the indeterminate h^d(0)=0. This is now truncated and shifted by M=(N-1)/2 to give the optimal, causal length-N FIR filter coefficients as h ( n ) = 1 - cos ( ω 0 ( n - M ) ) π ( n - M ) 0 ≤ n ≤ N - 1 and h(n)=0 otherwise. Again the corresponding derivation for an even length gives the same result as in (). The ideal amplitude response is shown in .

Ideal Frequency Response of an FIR Hilbert Transorm in the Passband and Lowpass Cutoff

Spline Transition Band Design All of the four lowpass filters described above exhibit the Gibbs phenomenon when truncated to a finite length. To remove this effect and to give a more explicit specification of the pass and stopband edges, a transition band is inserted between the pass and stopband. A transition function can be placed in this band to make the total desired amplitude response a continuous function. If we use a pth order spline as the transition function, the effect of adding this transition band to the basic lowpass filter ideal amplitude given in () is to multiply the ideal impulse response in () by a the Pth power of a sinc function to give h ^ d ( n ) = sin ( ω 0 n ) π n sin ( Δ n / p ) Δ n / p p where ω0=(ωs+ωp)/2 is the average band edge and Δ=(ωs-ωp)/2 is half the transition band width in radians per second normalized for one sample per second sampling rate , , . The spline produces a transition function which consists of p segments of pth order polynomials connected together so that p-1 derivatives are continuous at the junctions.

Ideal Lowpass Filter Amplitude with Order-p Spline Transition Function The optimal value of the exponent p is chosen as p=0.624(fs-fp)N (for a unity sampling rate) which minimizes the approximation error . Each of the four ideal lowpass filters derived in the previous can have a transition band added simply by multiplying their impulse response by the sinc weighting function as illustrated in (). shows an ideal unity gain filter amplitude response with examples of first, second, and tenth order spline transition functions. shows the ideal responses of the linear gain filter with fourth order spline transition function.

Ideal Increasing Amplitude Filter with Spline Transition Function

The Optimal Multiband Least Squared Error Design Method The optimal multiband design method consists of two somewhat independent parts. The first is the design of an optimal least squares lowpass filter with a transition band as described above or as calculated by an inverse FFT. The second part builds an optimal multiband filter from a combination of these optimal lowpass filters and is the main point of this . The unweighted least squared error linear phase FIR filter design problem is to find the filter coefficients that minimize the error defined by q = ∫ 0 π | A ( ω ) - A d ( ω ) | 2 d ω where A(ω) is the amplitude frequency response of the actual filter and Ad(ω) is the desired ideal amplitude response. This is done by truncating the inverse discrete time Fourier transform of Ad(ω). The difficulty is the analytical evaluation of the integral in the inverse transform . If a spline transition function is used, an analytical formula can be derived for the filter that minimizes (). The details of this result can be found in , . The infinitely long filters designed from the inverse discrete time Fourier transform of the ideal response have a frequency response which is the same as the ideal and, therefore, has no error. An ideal desired amplitude response can be formulated as the sum of simpler ideal lowpass filters, differentiators or Hilbert transformers together with their spline transition functions by A d ( ω ) = ∑ k K k A d k ( ω ) . where Adk(ω) is the desired lowpass response with a transition band in the kth band such as given in () or () and the Kk are chosen to build the desired Ad(ω). These Adk are the forms considered in the previous section along with any others that have analytical inverse DTFTs such as polynomials. Because of the linearity of the Fourier transform, a multiband ideal response can be constructed by simply adding and subtracting the impulse response of appropriate ideal lowpass filters. h ^ d ( n ) = ∑ k K k IDTFT { A d k ( ω ) } h ^ d ( n ) = ∑ k K k h ^ d k ( n ) Because of the orthogonality of the basis functions of the Fourier transform, the truncated sequence of the infinitely long impulse response h^d(n) will give an optimal approximation to Ad(ω) in a least squares sense. This argument allows no error weighting or “don't care" transition bands or traditional windowing methods. It does, however, allow the optimized spline transition functions , . Using these facts, an optimal multiband filter can be built up by successively adding and subtracting the impulse responses of optimal lowpass filters as done in (). For example a bandpass filter that approximates zero for 0<ω<ω1, has a spline transition band for ω1<ω<ω2, approximates one (or some other constant) for ω2<ω<ω3, has an independent second transition band for ω3<ω<ω4, and finally approximates zero for ω4<ω<π can be designed by first designing a simple lowpass filter with transition band ω3<ω<ω4 and then subtracting from its impulse response the impulse of a second lowpass filter designed with a transition band ω1<ω<ω2. A filter with two or more passbands can be designed by adding the impulse responses of two or more single passband filters. Indeed, a completely general design method can be formulated by alternately adding and subtracting lowpass filters starting with the highest frequency transition band and moving sequentially down to the lowest. If the ideal frequency response is not zero at ω=π, then one starts with a constant frequency response (an impulse in the time domain) and subtracts a lowpass filter (remember the length must be odd for this case). By scaling each lowpass filter, different gains are obtained in each band. Care must be taken that the constructed spline transition function properly fit the bands on both sides. This will not automatically happen if there are two adjacent bands with different slopes connected by one transitions function which are simply added together. It will automatically happen if each passband is separated by a stopband or if adjacent bands have the same slopes.

A Matlab Filter Design Program A Matlab program named fir3.m is given in the appendix of this book that will design optimal filters using the method described in the previous section. This particular program requires constant but arbitrary passband gains and uses a format for specifications similar to remez() in Matlab. It constructs the multiband filter from () by adding and subtracting optimal lowpass filters designed from the formula in () and calculated in the second program named fir3lp.m . The main program is given an even length vector f containing the normalized pass and stopband edges, including f=0 and f=1. It is also given an even length vector m containing the ideal response at each frequency in f . Because the lowpass filter has a constant passband, the ideal response of the multiband filter will have constant passbands. This means m will consist of adjacent terms that are equal. An example Matlab function call is given in the next section. The simple program listed in the appendix will design filters with constant gains in multiple passbands. From its construction it is easy to see how adding the use of the linear gain lowpass filter to the unity gain passband lowpass filter would allow designing optimal filters with linear gains in the passbands. By adding all four basic lowpass designs a calling program could be written that would automatically design one filter with a combination of all four characteristics. If the real and imaginary parts of a desired complex frequency response can be given in terms of the basic filters, nonlinear phase filters can be designed also. The programs are written to be consistent with Matlab's convention of normalizing for two samples per second sampling rate. The equations most of this book, however, are normalized for one sample per second.

Design Examples To show the results of using this new design approach, two examples of multiband filter design are presented here. The first is a filter with a stopband from ω=0 to ω=0.2, a transition band from ω=0.2 to ω=0.25, a passband with gain equal to 0.7 from ω=0.25 to ω=0.5, a transition band from ω=0.5 to ω=0.55, a passband with gain equal to 0.5 from ω=0.55 to ω=0.7, a transition band from ω=0.7 to ω=0.73, a stopband from ω=0.73 to ω=0.85, a transition band from ω=0.85 to ω=0.9, and a passband with gain equal one from ω=0.9 to ω=1. This is called with the Matlab program by

 
h  = fir3(51,[0 .2 .25 .5 .55 .7 .73 .85 .9 1],[0 0 .7 .7.5 .5 0 0 1 1])

and the amplitude response plot shown in a. The response for length of N=a01 is shown in Figures b and in c the zero locations are given. As an example of how versatile this approach can be, a length-101 linear phase multiband FIR filter was designed with different types of filtering being done in different bands. The signal frequencies in the band from 0<f<0.2 is differentiated,in the band from 0.25<f<0.4 is rejected, in 0.4<f<0.2 is hiltbert transformed, in 0<f<0.2 is rejected, and 0<f<0.2 is highpass filtered. In the transition bands between each of these processing bands, is an optimal spline transition function. This is a truly versatile multiband design technique with the only major limitation being that weighting is not possible. However, that limitation is removed in the next secession.

Frequency Response of Length-15 FIR Filter Designed by Least Squared Error

Weighted Least Integral Squares FIR Filter Design If the FIR filter design problem is posed as a weighted integral squared error approximation problem, a simple analytical design formula as in () or () is not possible (Recall that it is possible to easily introduce weights in the discrete approximation problem ()). In this section we consider a multiband generalization of an approach which is a mixture of analytical formulas and numerical solution of Toeplitz plus Hankel matrices which have been presented in , , . The most general definition of the linear phase weighted least squares FIR filter design problem , , defines the error measure as in () by q = 1 π ∫ Ω W ( ω ) |