Skip to content Skip to navigation Skip to collection information

OpenStax-CNX

You are here: Home » Content » Digital Signal Processing and Digital Filter Design (Draft) » Least Squared Error Design of FIR Filters

Navigation

Lenses

What is a lens?

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

This content is ...

Affiliated with (What does "Affiliated with" mean?)

This content is either by members of the organizations listed or about topics related to the organizations listed. Click each link to see a list of all content affiliated with the organization.
  • NSF Partnership display tagshide tags

    This collection is included inLens: NSF Partnership in Signal Processing
    By: Sidney Burrus

    Click the "NSF Partnership" link to see all content affiliated with them.

    Click the tag icon tag icon to display tags associated with this content.

  • Featured Content

    This collection is included inLens: Connexions Featured Content
    By: Connexions

    Click the "Featured Content" link to see all content affiliated with them.

Also in these lenses

  • UniqU content

    This collection is included inLens: UniqU's lens
    By: UniqU, LLC

    Click the "UniqU content" link to see all content selected in this lens.

  • Lens for Engineering

    This module and collection are included inLens: Lens for Engineering
    By: Sidney Burrus

    Click the "Lens for Engineering" link to see all content selected in this lens.

Recently Viewed

This feature requires Javascript to be enabled.

Tags

(What is a tag?)

These tags come from the endorsement, affiliation, and other lenses that include this content.
 

Least Squared Error Design of FIR Filters

Module by: C. Sidney Burrus. E-mail the author

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 [8], [1]. 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 [22].

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 ω 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 ω
(1)

where the approximation error as a function of frequency is defined by e(ω)=A(ω)-Ad(ω)e(ω)=A(ω)-Ad(ω) with A(ω)A(ω) being the amplitude response of the filter and Ad(ω)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 ϵ = A - A d = C a - A d
(2)

where CC is the matrix of cosines from Equation 48 from FIR Digital Filters, aa is the vector of half of the filter coefficients from Equation 48 from FIR Digital Filters, and AdAd is the vector of samples of the ideal desired amplitude response. The number of samples of the amplitude response is LL 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 ϵ q = ϵ T ϵ
(3)

except for a scale factor of 1L1L.

This could also be posed for the general phase problem by using H(ωk)H(ωk) rather than A(ωk)A(ωk) and h(n)h(n), the actual impulse response, rather than a(n)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 qq as a sum of the squared differences between the actual and the desired frequency response over a set of LL frequency samples. This error function is defined as

q = 1 L k = 0 L - 1 | H ( ω k ) - H d ( ω k ) | 2 q = 1 L k = 0 L - 1 | H ( ω k ) - H d ( ω k ) | 2
(4)

where Hd(ωk)Hd(ωk) are the LL samples of the desired response. This problem is easier to formulate and solve if the frequency samples are equally spaced as in Equation 8 from FIR Filter Design by Frequency Sampling or Interpolation which gives

ω k = 2 π k / L ω k = 2 π k / L
(5)

and the problem is restricted to linear-phase filters where the real-valued amplitude A(ω)A(ω) can be approximated rather than the complex frequency response H(ω)H(ω). For approximations to a complex response, see "Complex and Minimum Phase Approximation".

Linear phase and equally spaced samples cause Equation 4 to become

q = 1 L k = 0 L - 1 | A ( 2 π k / L ) - A d ( 2 π k / L ) | 2 q = 1 L k = 0 L - 1 | A ( 2 π k / L ) - A d ( 2 π k / L ) | 2
(6)

or with a simpler notation

q = 1 L k = 0 L - 1 | A k - A d k | 2 q = 1 L k = 0 L - 1 | A k - A d k | 2
(7)

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 Equation 7 in the frequency domain can also be calculated in the time domain by

q = n = 0 L - 1 | h ( n ) - h d ( n ) | 2 q = n = 0 L - 1 | h ( n ) - h d ( n ) | 2
(8)

where hd(n)hd(n) is the length-L symmetric FIR filter that has the LL frequency response amplitude samples AdkAdk. This may be calculated by the frequency sampling method in the section Four Types of Linear-Phase FIR Filters using the special formulas such as Equation 8 from FIR Filter Design by Frequency Sampling or Interpolation for length LL or the inverse DFT. The filter to be designed has a length-N symmetric impulse response h(n)h(n) with LL frequency response samples AkAk.

Because the filter h(n)h(n) is of length-N and symmetric, the error equation Equation 8 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 q = n = - M M | h ^ ( n ) - h ^ d ( n ) | 2 + 2 n = M + 1 ( L - 1 ) / 2 | h ^ d ( n ) | 2
(9)

where h^(n)h^(n) and h^d(n)h^d(n) are the inverse DTFTs of AkAk and AdkAdk respectively, which means they are the h(n)h(n) and hd(n)hd(n) shifted to be symmetic about n=0n=0. This requires the number of frequency samples LL must be odd.

Equation 9 clearly shows that to minimize qq, the NN values of h(n)h(n) are chosen to be equal to the equivalent NN values of hd(n)hd(n) making the first sum equal zero. In other words, h(n)h(n) is obtained by symmetrically truncating hd(n)hd(n). The residual error is then given by the second summation above. An examination of the residual error as a function of NN may aid in the choice of the filter length NN.

For the Type 1 linear-phase FIR filter (described in the section Four Types of Linear-Phase FIR Filters) which has an odd length NN and an even-symmetric impulse response, the LL equally spaced samples of the frequency response from Equation ? from Fir Digital Filters gives

A k = n = 0 M - 1 2 h ( n ) cos ( 2 π ( M - n ) k / L ) + h ( M ) A k = n = 0 M - 1 2 h ( n ) cos ( 2 π ( M - n ) k / L ) + h ( M )
(10)

for k=0,1,2,....,L-1k=0,1,2,....,L-1, where M=(N-1)/2M=(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)a(n)). It was noted in the section Frequency Sampling Filter Design by Formulas that it is also a cosine transform and it can be shown that this transformation is orthogonal over the independent values of AkAk, just as the DFT is.

The desired ideal amplitude gives the ideal impulse response hd(n)hd(n) from Equation 29 from FIR Digital Filters by

h d ( n ) = 1 N [ A d 0 + k = 1 M - 1 2 A d k cos ( 2 π ( n - M ) k / N ) ] . h d ( n ) = 1 N [ A d 0 + k = 1 M - 1 2 A d k cos ( 2 π ( n - M ) k / N ) ] .
(11)

for n=0,1,,L-1n=0,1,,L-1. This is used in Equation 9, and is the ideal impulse response that is truncated and shifted to give a causal, symmetric h(n)h(n).

Use of the alternative equally-spaced sampling in Equation 9 from FIR Filter Design by Frequency Sampling or Interpolation, which has no sample at zero frequency, requires hd(n)hd(n) be calculated from Equation 11 from FIR Filter Design by Frequency Sampling or Interpolation and Equation 13 from FIR Filter Design by Frequency Sampling or Interpolation. The Type 2 filters with even NN are developed in a similar way and use the design formulas Equation 36 from FIR Digital Filters and Equation 37 from FIR Digital Filters. 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 Equation 5 from FIR Filter Design by Frequency Sampling or Interpolation or Equation 11 from FIR Filter Design by Frequency Sampling or Interpolation or the IDFT, then to symmetrically truncate it to the desired odd-length N and shift it to make h(n)h(n) causal. To design an even-length Type 2 filter , start with an even-length-L frequency-sampling design from Equation 7 from FIR Filter Design by Frequency Sampling or Interpolation or Equation 13 from FIR Filter Design by Frequency Sampling or Interpolation 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 LL 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 NN odd in Equation 2 from FIR Filter Design by Frequency Sampling or Interpolation are given by

A ( ω k ) = n = 1 M 2 h ( M - n ) cos ( ω k n ) + h ( M ) A ( ω k ) = n = 1 M 2 h ( M - n ) cos ( ω k n ) + h ( M )
(12)

for k=0,1,,L-1k=0,1,,L-1. This relates the LL frequency samples A(ωk)A(ωk) to the M+1 independent values of the symmetric length-N impulse response h(n). In the design problem where the AkAk are given and the values for h(n) are to be found, this represents LL equations with M+1 unknowns. Because of the symmetries of A(ω)A(ω) shown in Figure 5 from FIR Digital Filters, only half of the LL values of AkAk are independent; however, in some cases, to have proper weights on all LL samples, all must be calculated.

Equation 12 sampled at LL arbitrary frequencies can be written as a matrix equation

C a = A C a = A
(13)

where aa is an M+1M+1 length vector with elements which are the first half of h(n)h(n). CC is an LL by (M+1)(M+1) matrix of the cosine terms from Equation 12, and AA is a length-L vector of the frequency samples A(ωk)A(ωk).

If the formula for the calculation of LL values of the frequency response of a length-N FIR filter in Equation 10 is used to define an error vector of differences as defined in Equation 3 and the result is written in the matrix formulation of Equation 48 from FIR Filter Design by Frequency Sampling or Interpolation, the error becomes

C a = A = A d + e C a = A = A d + e
(14)

or

C a - A d = e C a - A d = e
(15)

where ee is a vector of differences between the actual and desired samples of the frequency response. The error measure defined in Equation 4 becomes the quadratic form

q = e T e q = e T e
(16)

For L>NL>N, equation Equation 13 is over determined and cannot, in general, be solved for aa. The filter design error measure is the norm of ee, as given in Equation 16. This error measure is minimized by making ee orthogonal to the columns of CC in Equation 15. Multiplying both sides of Equation 14 by the transpose of CC gives

C T C a = C T A d + C T e C T C a = C T A d + C T e
(17)

In order for qq to be minimum, ee must be orthogonal to the columns of CC and, therefore, CTeCTe must be zero. Hence, the optimal aa must satisfy the “normal equations" [10], [20], [9] which are

C T C a = C T A d C T C a = C T A d
(18)

and which can be rewritten in terms of the pseudo-inverse [10], [9] as

a = [ C T C ] - 1 C T A d a = [ C T C ] - 1 C T A d
(19)

If L=NL=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>NL>N, there are still only M+1M+1 equations to be solved. For L>>NL>>N, equation Equation 10 may be ill-conditioned, and Equation 19 should not be used to solve them. Special methods will be necessary to avoid serious numerical problems [9].

If a weighted error function is desired, Equation 4 is modified to give

q = 1 L k = 0 L - 1 W k | A ( ω k ) - A d ( ω k ) | 2 q = 1 L k = 0 L - 1 W k | A ( ω k ) - A d ( ω k ) | 2
(20)

The normal equations of Equation 18 become

C T W C a = C T W A d C T W C a = C T W A d
(21)

where WW 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 LL 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, Equation 21 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 FF is obtained from sampling Equation 28 from FIR Digital Filters and the full h(n)h(n) is used. For the special case of the equally-spaced frequency samples and linear- phase filter with unity weighting, the solution of Equation 18 or Equation 21 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 Equation 18 or Equation 21, 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 Equation 19 should not be solved directly. For large LL, 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=LN=L, the same results are obtained, thus verifying the theoretical prediction. As LL becomes larger compared to NN, more control is exerted over the behavior between the original sample points. As LL becomes large compared to NN, 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.

Figure 1: Frequency Response of Length-15 FIR Filter Designed by Least Squared Error
 This graph is labeled Length-15 FIR lowpass filter designed by least squared error. The x axis is labeled Normalized Frequency, and the y axis is labeled Amplitude Response, A. This figure consist of a box formed by the x and y axes and a line that extends perpendicularly from the y at y=14 and a line extending perpedicularly to the x axis at x=.5. In addition to this box there is a waveform that begins on the y axis just below the line perpendicular to the y axis. The waveform travels above the line and then back below and then back across the line further above the line. Then the line takes a very negative slope crossing the line perpendicular to the y axis and also the line perpendicular to the x axis. The line continues below the x axis crossing the axis just before x=.6. The line then undulates above and below the axis until the graph ends. There is an arrow pointing to the middle of the wave between existing inside the box. This arrow labels the area Achieved Amplitude and then below that is another arrow that labels the line perpendicular to the x axis Desired Amplitude.

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 Figure 1 even though the length and band edge are the same. Notice there is less over shoot but more ripple near f=0f=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.4f=0.4 and f=0.6f=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 Figure 2.

Figure 2: Frequency Response of Length-15 FIR Filter with a Transition Band Designed by Least Squared Error
Figure 2 (figDLSEtb.png)

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 ω q = 1 π 0 π | A d ( ω ) - A ( ω ) | 2 d ω
(22)

where Ad(ω)Ad(ω) is the desired ideal amplitude response, A(ω)=na(n)cos(ω(M-N)n)A(ω)=na(n)cos(ω(M-N)n) is the achieved amplitude response with the length h(n)h(n) related to h(n)h(n) by Equation 29 from FIR Digital Filters. This integral squared error is approximated by the discrete squared error defined in Equation 22 for L>>NL>>N which in some cases is much easier to minimize. However for some very useful cases, formulas can be found for h(n)h(n) that minimize Equation 22 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 Equation 22, 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 ω . q = n = - | h d ( n ) - h ( n ) | 2 = 1 π 0 π | A d ( ω ) - A ( ω ) | 2 d ω .
(23)

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 Equation 23 into two parts, one of which depends on h(n)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 q = n = - M M | h d ( n ) - h ( n ) | 2 + 2 n = M + 1 | h d ( n ) | 2
(24)

Again, we see that the minimum qq is achieved by using h(n)=hd(n)h(n)=hd(n) for -MnM-MnM. In other words, the infinitely long hd(n)hd(n) is symmetrically truncated to give the optimal least integral squared error approximation. The problem then becomes one of finding the hd(n)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<ω<ω00<ω<ω0 and zero stopband gain for ω0<ω<πω0<ω<π similar to those in Figure 8a from FIR Digital Filters and Figure 1. This gives

A d ( ω ) = { 1 0 < ω < ω 0 0 ω 0 < ω < π A d ( ω ) ={ 1 0 ω ω 0 0 ω 0 ω π
(25)

as the ideal desired amplitude response. The ideal shifted filter coefficients are the inverse DTFT from Equation 15 from Chebyshev or Equal Ripple Error Approximation Filters of this amplitude which for NN odd are given by

h ^ d ( n ) = 1 π 0 π A d ( ω ) cos ( ω n ) d ω h ^ d ( n ) = 1 π 0 π A d ( ω ) cos ( ω n ) d ω
(26)
= 1 π 0 ω 0 cos ( ω n ) d ω = ω 0 π sin ( ω 0 n ) ω 0 n = 1 π 0 ω 0 cos ( ω n ) d ω = ω 0 π sin ( ω 0 n ) ω 0 n
(27)

which is sometimes called a “sinc" function. Note h^d(n)h^d(n) is generally infinite in length. This is now symmetrically truncated and shifted by M=(N-1)/2M=(N-1)/2 to give the optimal, causal length-NN FIR filter coefficients as

h ( n ) = ω 0 π sin ( ω 0 ( n - M ) ) ω 0 ( n - M ) for 0 n N - 1 h ( n ) = ω 0 π sin ( ω 0 ( n - M ) ) ω 0 ( n - M ) for 0 n N - 1
(28)

and h(n)=0h(n)=0 otherwise. The corresponding derivation for an even length starts with the inverse DTFT in Equation 5 from Constrained Approximation and Mixed Criteria 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 ) h ^ d = 1 π 0 π A d ( ω ) cos ( ω ( n + 1 / 2 ) ) d ω = ω 0 π sin ( ω 0 ( n + 1 / 2 ) ) ω 0 ( n + 1 / 2 )
(29)

which when truncated and shifted by N/2N/2 gives the same formula as for the odd length design in Equation 28 but one should note that M=(N-1)/2M=(N-1)/2 is not an integer for an even NN.

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 Figure 3 For NN odd, the ideal infinitely long shifted filter coefficients are the inverse DTFT of this amplitude given by

A d ( ω ) = { 1 π ω 0 < ω < ω 0 0 ω 0 < ω < π A d ( ω ) ={ 1 π ω 0 ω ω 0 0 ω 0 ω π
(30)

and illustrated in Figure 3 For N odd, the ideal infinitely and shifted filter coefficients are the inverse DTFT of this amplitude given by

h ^ d ( n ) = 1 π 0 ω 0 ( ω π ) cos ( ω n ) d ω = cos ( ω 0 n ) - 1 π 2 n 2 + ω 0 sin ( ω 0 n ) π 2 n h ^ d ( n ) = 1 π 0 ω 0 ( ω π ) cos ( ω n ) d ω = cos ( ω 0 n ) - 1 π 2 n 2 + ω 0 sin ( ω 0 n ) π 2 n
(31)

with the indeterminate h^d(0)=ω022π2h^d(0)=ω022π2. This is now truncated and shifted by M=(N-1)/2M=(N-1)/2 to give the optimal, causal length-NN 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 h ( n ) = cos ( ω 0 ( n - M ) ) - 1 π 2 ( n - M ) 2 + ω 0 sin ( ω 0 ( n - M ) ) π 2 ( n - M ) for 0 n N - 1
(32)

and h(n)=0h(n)=0 otherwise. The corresponding derivation for an even length starts with the inverse DTFT for a shifted even length filter in Equation 15 from Chebyshev or Equal Ripple Error Approximation Filters and after shifting by N/2N/2 gives the same result as Equation 32.

Figure 3: Ideal Frequency Response of an FIR Filter with Increasing Gain in the Passband and Lowpass Cutoff
The graph is labeled Linearly Increasing Amplitude wiath a Lowpass Cutoff. The x axis is labeled Normalized Frequency and the y axis is Amplitude Response. The graph consist of a line extending from the origin with a positive slope and then the line proceeds straight down at x=.8 until it hits the x axis at a right angle.

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 Equation 30 and Figure 3 but, since this case has an odd symmetric impulse response, the inverse DTFT uses sine functions which for odd NN gives

h ^ d ( n ) = 1 π 0 ω 0 ( 1 π ω ) sin ( ω n ) d ω = sin ( ω 0 n ) π 2 n 2 - ω 0 cos ( ω 0 n ) π 2 n h ^ d ( n ) = 1 π 0 ω 0 ( 1 π ω ) sin ( ω n ) d ω = sin ( ω 0 n ) π 2 n 2 - ω 0 cos ( ω 0 n ) π 2 n
(33)

with the indeterminate h^d(0)=0h^d(0)=0. This is now truncated and shifted by M=(N-1)/2M=(N-1)/2 to give the optimal, causal length-NN 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 h ( n ) = sin ( ω 0 ( n - M ) ) π 2 ( n - M ) 2 - ω 0 cos ( ω 0 ( n - M ) ) π 2 ( n - M ) for 0 n N - 1
(34)

and h(n)=0h(n)=0 otherwise. Again the corresponding derivation for an even length gives the same result as in Equation 34. Note this very general single formula includes as special cases the odd and even length full band (ω0=πω0=π) differentiator given in [16]. 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 [15] combined with a lowpass filter can also be analytically evaluated. The ideal amplitude response is the same as Equation 25 but with a constant phase shift of ϕ=π/2ϕ=π/2. Since this case has an odd symmetric impulse response, the inverse DTFT uses sine functions which for odd NN which gives

h ^ d ( n ) = 1 π 0 ω 0 sin ( ω n ) d ω = 1 - cos ( ω 0 n ) π n h ^ d ( n ) = 1 π 0 ω 0 sin ( ω n ) d ω = 1 - cos ( ω 0 n ) π n
(35)

with the indeterminate h^d(0)=0h^d(0)=0. This is now truncated and shifted by M=(N-1)/2M=(N-1)/2 to give the optimal, causal length-NN FIR filter coefficients as

h ( n ) = 1 - cos ( ω 0 ( n - M ) ) π ( n - M ) 0 n N - 1 h ( n ) = 1 - cos ( ω 0 ( n - M ) ) π ( n - M ) 0 n N - 1
(36)

and h(n)=0h(n)=0 otherwise. Again the corresponding derivation for an even length gives the same result as in Equation 36. The ideal amplitude response is shown in Figure 4.

Figure 4: Ideal Frequency Response of an FIR Hilbert Transorm in the Passband and Lowpass Cutoff
The Cartesian graph contains in essence two squares formed by the presence of lines in both the bottom left and top left areas of the graph. In the bottom left square is formed by a line that runs along y=0, a line perpendicular to this at x=-.75. Another line runs parallel to the y=0 line from (-.75,-1) to (0,-1)Finally there is another line running vertically along x=0 from (0,-1) to (0,1). To the right of the top of this line there is a line running parallel to the y axis from (0,1) to (.75,1), and then another line extends down from (.75,1) to (.75,0). The First line at y=0 forms the bottom of this square.

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 pthpth order spline as the transition function, the effect of adding this transition band to the basic lowpass filter ideal amplitude given in Equation 25 is to multiply the ideal impulse response in Equation 27 by a the PthPth power of a sinc function to give

h ^ d ( n ) = sin ( ω 0 n ) π n sin ( Δ n / p ) Δ n / p p h ^ d ( n ) = sin ( ω 0 n ) π n sin ( Δ n / p ) Δ n / p p
(37)

where ω0=(ωs+ωp)/2ω0=(ωs+ωp)/2 is the average band edge and Δ=(ωs-ωp)/2Δ=(ωs-ωp)/2 is half the transition band width in radians per second normalized for one sample per second sampling rate [19], [16], [4]. The spline produces a transition function which consists of pp segments of pthpth order polynomials connected together so that p-1p-1 derivatives are continuous at the junctions.

Figure 5: Ideal Lowpass Filter Amplitude with Order-p Spline Transition Function
This graph is labeled Ideal Lowpass Filter Amplitude with Order-p Spline Transition Function. The x axis is labeled Normalized frequency and the y axis is labeled Amplitude Response, A_d. At y=0 there is another line running the length of the graph. The segment of the line between (.3,0) to (.7,0) is labeled Transition Band. This graph consist of three different lines all originating at (0,1). All three of the lines start out as horizontal lines. The first line to deviate from this path is labeled p=1 and  at the point where the Transition Band starts to where it ends the line takes a negative diagonal slope. The next line to deviate from the initial horizontal path is labeled p=2, and when it reaches the area of the transition band, the line takes a gentle negative slope and then ends by curving to the right before reaching the x axis. The final line labeled P=10. This line follows the initial horizontal path the longest and then takes the most drastic negative slope of the three down (almost a vertical line) and then curves to the right just a little before ending on the x axis.

The optimal value of the exponent pp is chosen as p=0.624(fs-fp)Np=0.624(fs-fp)N (for a unity sampling rate) which minimizes the approximation error [4]. 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 Equation 37. Figure 5 shows an ideal unity gain filter amplitude response with examples of first, second, and tenth order spline transition functions. Figure 6 shows the ideal responses of the linear gain filter with fourth order spline transition function.

Figure 6: Ideal Increasing Amplitude Filter with Spline Transition Function
The graph is labeled Ideal Increasing Amplitude with Spline Transition Function. The x axis is labeled Normalized Frequency and the y axis is Amplitude Response. The graph consist of a line extending from the origin with a positive slope and then the line proceeds straight down at about x=.75. Right before reaching the x axis the vertical line curves a little to the right just before (.8,0).

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 [5].

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 ω q = 0 π | A ( ω ) - A d ( ω ) | 2 d ω
(38)

where A(ω)A(ω) is the amplitude frequency response of the actual filter and Ad(ω)Ad(ω) is the desired ideal amplitude response. This is done by truncating the inverse discrete time Fourier transform of Ad(ω)Ad(ω). The difficulty is the analytical evaluation of the integral in the inverse transform [16]. If a spline transition function is used, an analytical formula can be derived for the filter that minimizes Equation 38. The details of this result can be found in [4], [5].

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 ( ω ) . A d ( ω ) = k K k A d k ( ω ) .
(39)

where Adk(ω)Adk(ω) is the desired lowpass response with a transition band in the kthkth band such as given in Equation 25 or Equation 30 and the KkKk are chosen to build the desired Ad(ω)Ad(ω). These AdkAdk 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 IDTFT { A d k ( ω ) }
(40)
h ^ d ( n ) = k K k h ^ d k ( n ) h ^ d ( n ) = k K k h ^ d k ( n )
(41)

Because of the orthogonality of the basis functions of the Fourier transform, the truncated sequence of the infinitely long impulse response h^d(n)h^d(n) will give an optimal approximation to Ad(ω)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 [4], [5].

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 Equation 40. For example a bandpass filter that approximates zero for 0<ω<ω10<ω<ω1, has a spline transition band for ω1<ω<ω2ω1<ω<ω2, approximates one (or some other constant) for ω2<ω<ω3ω2<ω<ω3, has an independent second transition band for ω3<ω<ω4ω3<ω<ω4, and finally approximates zero for ω4<ω<πω4<ω<π can be designed by first designing a simple lowpass filter with transition band ω3<ω<ω4ω3<ω<ω4 and then subtracting from its impulse response the impulse of a second lowpass filter designed with a transition band ω1<ω<ω2ω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 [11] 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 Equation 41 by adding and subtracting optimal lowpass filters designed from the formula in Equation 37 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=0f=0 and f=1f=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ω=0 to ω=0.2ω=0.2, a transition band from ω=0.2ω=0.2 to ω=0.25ω=0.25, a passband with gain equal to 0.7 from ω=0.25ω=0.25 to ω=0.5ω=0.5, a transition band from ω=0.5ω=0.5 to ω=0.55ω=0.55, a passband with gain equal to 0.5 from ω=0.55ω=0.55 to ω=0.7ω=0.7, a transition band from ω=0.7ω=0.7 to ω=0.73ω=0.73, a stopband from ω=0.73ω=0.73 to ω=0.85ω=0.85, a transition band from ω=0.85ω=0.85 to ω=0.9ω=0.9, and a passband with gain equal one from ω=0.9ω=0.9 to ω=1ω=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 Figure 7a. The response for length of N=101N=101 is shown in Figures Figure 7b and in Figure 7c 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 with frequencies in the band from 0<f<0.20<f<0.2 is differentiated, in the band from 0.23<f<0.40.23<f<0.4 is rejected, in 0.43<f<0.60.43<f<0.6 is Hilbert transformed, in 0.63<f<0.80.63<f<0.8 is rejected, and 0.83<f<1.00.83<f<1.0 is highpass filtered. In the transition bands between each of these processing bands, there is an optimal spline transition function. The amplitude response is shown in Figure 7d. 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.

Figure 7: Frequency Response and Zero Locations of FIR Filters Designed by Least Squared Error
This figure consist of four different graphs. The top left graph is labeled a. Length-51 Multiband FIR Filter. The y axis is labeled Magnitude Response. This graph consist of a wave form that starts with a squiggly section that is pretty much on top of the y=0 line and then the wave rises quickly at an almost vertical slope to y=.7. Then the wave squiggles til (.5 .7) where the slope takes an extremely negative slope and then it squiggles a little and the plummets back to the x axis and then squiggles a little and then takes a drastically positive slope almost to the top of the graph where is squiggles a tiny bit before running of the graph.  The second graph is labeled b. Length-201 Multiband FIR Filter. The graph is similar to the first graph but every place where there was squiggling the area is now a flat horizontal line. The bottom left graph is the most complicated of the four graphs. It is labeled c. Multiband Filter Zero Locations. the x axis is labeled Real Part of Z, and the y axis is labeled the Imaginary Part of z. In general this image consist of a large circle centered on the origin with little hollow circles in, on and around it. There three distinct areas around the perimeter of the circle where there are little hollow circles present on the perimeter of the circle. On the left hand side there are two groups of circles: one at the center of the bottom left quarter of the circle and then another at the center of the top left quarter of the circle. The other group is on the center of the right half of the circle. All of these circle groups consist of overlapping little circles. The two on the left half are groups of four circles and the one on the right half consist of 12 little circles. There are also little circles present inside the larger circle. There are three circles grouped together to the right of the left most area of the larger circle. There are six grouped together between the top left grouping of perimeter circles and the right half group and the same grouping mirror between the bottom left group of perimeter circles and the right group. There are also little circle on the outside of the larger circle. There are three to the left of the far left of the circle. There are also two groups of six circles above and below the top and bottom extremes of the circle. The final graph is more similar to first two graphs. It is labeled d. Multiband, Multitype FIR Filter. It consist of a waveform that has an extremely positive slope to begin with which peaks at (.25,1) and then immediately takes an equally negative slope down to the x axis. The wave then wavers a little along the x axis and then shoots straight up to the same peak as before and wavers a bit and then falls back the x axis wavers and repeats.

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 Equation 28 or Equation 40 is not possible (Recall that it is possible to easily introduce weights in the discrete approximation problem Equation 21). In this section we consider a multiband generalization [5] of an approach which is a mixture of analytical formulas and numerical solution of Toeplitz plus Hankel matrices which have been presented in [6], [21], [4].

The most general definition of the linear phase weighted least squares FIR filter design problem [6], [7], [2] defines the error measure as in Equation 10 from Constrained Approximation and Mixed Criteria by

q = 1 π Ω W ( ω ) | A d ( ω ) - A ( ω ) | 2 d ω . q = 1 π Ω W ( ω ) | A d ( ω ) - A ( ω ) | 2 d ω .
(42)

where ΩΩ is the set of frequencies that contribute to the error.

We set up the conditions for minimizing the error in Equation 42 for odd NN by using the same approach used in [4] which substitutes

A ( ω ) = n = 0 M a ( n ) cos ( ω n ) A ( ω ) = n = 0 M a ( n ) cos ( ω n )
(43)

from Equation 34 from FIR Digital Filters into Equation 42, differentiates qq in respect to each a(m)a(m), and then sets it equal to zero to give

1 π Ω W ( ω ) A d ( ω ) cos ( ω m ) d ω = n = 0 M a ( n ) 2 π Ω W ( ω ) cos ( ω n ) cos ( ω m ) d ω 1 π Ω W ( ω ) A d ( ω ) cos ( ω m ) d ω = n = 0 M a ( n ) 2 π Ω W ( ω ) cos ( ω n ) cos ( ω m ) d ω
(44)

where we can obtain the h(n)h(n) from the a(n)a(n) by the scaling and shifting in Equation 35 from FIR Digital Filters. We denote this in matrix form by

A w = C w a A w = C w a
(45)

with the elements of the M+1M+1 by M+1M+1 matrix CwCw as

c w ( m , n ) = 2 π Ω W ( ω ) cos ( ω n ) cos ( ω m ) d ω c w ( m , n ) = 2 π Ω W ( ω ) cos ( ω n ) cos ( ω m ) d ω
(46)

and the M+1M+1 by 1 vector of intermediate values awaw is given by an inverse DTFT of the weighted ideal amplitude response in

A w ( n ) = 1 π Ω W ( ω ) A d ( ω ) cos ( ω n ) d ω for 1 n N - 1 A w ( n ) = 1 π Ω W ( ω ) A d ( ω ) cos ( ω n ) d ω for 1 n N - 1
(47)

and

A w ( 0 ) = 1 2 π Ω W ( ω ) A d ( ω ) d ω . A w ( 0 ) = 1 2 π Ω W ( ω ) A d ( ω ) d ω .
(48)

Solving Equation 45 for the optimal aa which minimizes the integral weighted squared error Equation 42 is formally done by

a = C w - 1 A w a = C w - 1 A w
(49)

and more accurately done by special numerical algorithms. The case for even NN is easily derived by using Equation 40 from FIR Digital Filters to derive Equation 44. The actual length-NN filter coefficients h(n)h(n) are then found from a(n)a(n) using Equation 41 from FIR Digital Filters. Note that if the weighting is unity across the pass, transition, and stop bands, CwCw is the identity matrix. CwCw gives the effects of the weighting.

A similar formula was derived by Fleischer [6], Tufts, Rorabacher and Mosier [21], by Schüssler [17], by Oetken, Parks, and Schüssler [14], and by Burrus, Soewito, and Gopinath [4], [5] in addressing similar problems.

If the integrals in Equation 46 and Equation 47 can be analytically evaluated, the solution of the weighted squared error approximation is obtained by solving M+1M+1 equations. Fortunately these equations can be analytically evaluated for several interesting cases as was done for the unweighted case. The even length-NN case is derived in a similar way using Equation 40 from FIR Digital Filters and Equation 41 from FIR Digital Filters. An alternate formulation could modify the first column.

In most practical situations where specifications are set in the frequency domain, these filters are described in terms of frequency bands. We have already seen the idea of single pass, stop, and transition bands. We now allow multiple pass and stopbands separated by multiple transition bands. In order to obtain analytical solutions of Equation 46 and Equation 47 and to be consistent with usual practice, we restrict ourselves to constant weights over each separately defined frequency band. The error in Equation 42 now becomes

q = 1 π k W k ω k ω k + 1 | A d k ( ω ) - A ( ω ) | 2 d ω . q = 1 π k W k ω k ω k + 1 | A d k ( ω ) - A ( ω ) | 2 d ω .
(50)

where the weights are constant over each band and are given by W(ω)=WkW(ω)=Wk in the kthkth band defined by ωk<ω<ωk+1ωk<ω<ωk+1. The desired amplitude Adk(ω)Adk(ω) is likewise defined in the kthkth band and is hopefully simple enough to allow analytical evaluation of the formula Equation 26 for the ideal impulse response.

The form of Equation 50 causes Equation 46 to become

c w ( m , n ) = 2 π k W k ω k ω k + 1 cos ( ω n ) cos ( ω m ) d ω . c w ( m , n ) = 2 π k W k ω k ω k + 1 cos ( ω n ) cos ( ω m ) d ω .
(51)

which has an analytical solution as given by

c w ( m , n ) = 1 π k W k sin ( n - m ) ω k + 1 - sin ( n - m ) ω k ( n - m ) + sin ( n + m ) ω k + 1 - sin ( n + m ) ω k ( n + m ) c w ( m , n ) = 1 π k W k sin ( n - m ) ω k + 1 - sin ( n - m ) ω k ( n - m ) + sin ( n + m ) ω k + 1 - sin ( n + m ) ω k ( n + m )
(52)

which for FF band edges has terms that are indeterminate for n=m0n=m0 with values

c w ( n , n ) = 1 π k = 1 F - 2 W k ( ω k + 1 - ω k ) + sin ( 2 n ω k + 1 ) - sin ( 2 n ω k ) 2 n + W F - 1 π c w ( n , n ) = 1 π k = 1 F - 2 W k ( ω k + 1 - ω k ) + sin ( 2 n ω k + 1 ) - sin ( 2 n ω k ) 2 n + W F - 1 π
(53)

and for n=m=0n=m=0 as

c w ( 0 , 0 ) = 1 π k = 1 F - 2 W k 2 ( ω k + 1 - ω k ) + W F - 1 2 π c w ( 0 , 0 ) = 1 π k = 1 F - 2 W k 2 ( ω k + 1 - ω k ) + W F - 1 2 π
(54)

Since the matrix elements are functions of (n-m)(n-m) and (n+m)(n+m), CC is the sum of a Toeplitz and a Hankel matrix. This matrix can always be calculated and it simply depends on the set of band edges ωkωk and the band weights WkWk but not on the ideal amplitude response Ad(ω)Ad(ω). The case for even NN is similar but uses Equation 8 from Constrained Approximation and Mixed Criteria rather than Equation 7 from Constrained Approximation and Mixed Criteria with Equation 42 to derive an appropriate form of Equation 52.

If there are FF distinct band edges ωkωk, the first and last are ω1=0ω1=0 and ωF=πωF=π. This means part of the first term in the sum of Equation 52 is always zero and part of the last is zero except when n=m=0n=m=0 where it is ππ. Using these facts allows Equation 52 to be written

c w ( m , n ) = 1 π k = 1 F - 2 ( W k - W k + 1 ) sin ( n - m ) ω k + 1 n - m + sin ( n + m ) ω k + 1 n + m c w ( m , n ) = 1 π k = 1 F - 2 ( W k - W k + 1 ) sin ( n - m ) ω k + 1 n - m + sin ( n + m ) ω k + 1 n + m
(55)

which, together with appropriately modified Equation 53 and Equation 54, are good forms for programming. The Matlab program in the appendix contains the details.

Applying the form of Equation 50 to Equation 47 gives

a w ( n ) = 1 π k W k ω k ω k + 1 A d k ( ω ) cos ( ω n ) d ω . a w ( n ) = 1 π k W k ω k ω k + 1 A d k ( ω ) cos ( ω n ) d ω .
(56)

These integrals have been evaluated for the four basic filter types - constant gain passband lowpass filter, linear gain passband lowpass filter, differentiator plus lowpass filter, and Hilbert transformer plus lowpass filter - giving simple design formulas in Equation 28, Equation 32, Equation 34, and Equation 36.

Each basic filter type plus the effects of a transition band can be calculated and combined according to Equation 37. An example low pass filter with a weight of W1W1 in the passband and W2W2 in the transition band is given for odd NN gives for the intermediate coefficients h^w(n)h^w(n) from Equation 47 are

a w ( n ) = W 1 sin ( ω 2 n ) - sin ( ω 1 n ) π n + W 2 sin ( ω 0 n ) π n sin ( Δ n / p ) Δ n / p p - sin ( ω 1 n ) π n a w ( n ) = W 1 sin ( ω 2 n ) - sin ( ω 1 n ) π n + W 2 sin ( ω 0 n ) π n sin ( Δ n / p ) Δ n / p p - sin ( ω 1 n ) π n
(57)

A similar expression can be derived for even NN using Equation 8 from Constrained Approximation and Mixed Criteria.

This means the left hand vector in Equation 45 can be calculated as a weighted sum of inverse DTFTs such as in Equation 32 if the ideal desired amplitude can be constructed from the four basic types in FIR Digital Filters, each with optimal transition bands.

If one or more of the integrals in Equation 56 has no analytical solution, aw(n)aw(n) can be calculated numerically using a truncated weighted sum of inverse DFTs of a dense sampling of Adk(ω)Adk(ω) or made up of the passbands calculated from inverse DFTs and the transition bands added by multiplying appropriately by sinc functions since constructing an optimal spline transition function to be sampled would not be easy.

This gives a very powerful design method that allows multi band weighted least squares design of FIR filters. The calculation of the matrix CwCw in Equation 45 is always possible using Equation 52. Because using a true “don't care" transition band with a weight of zero might causes ill conditioning of Equation 45 for (fs-fp)N>12(fs-fp)N>12 as discussed in [4], one can add a spline transition function in Ad(ω)Ad(ω) to the definition in Equation 25 as done in Equation 57 and [4]. A very small weight used in the transition bands together with a spline transition function will improve the conditioning of Equation 45 with minor degradation of the optimality. This point needs further evaluation.

By using an inverse FFT perhaps plus a sinc induced transition function to calculate the components of Equation 56, this method can be used to design arbitrary shaped passbands. It can also be used for complex approximation by applying it to the real and imaginary parts of the desired Hd(ω)Hd(ω) separately and using the full, nonsymmetric h(n)h(n).

The form of the simultaneous equations Equation 45 that must be solved to design a filter by this method is interesting. If the weights in all pass, stop, and transition bands are unity, the CwCw matrix is the identity matrix and h^wh^w contains the filter coefficients. As the weights become less and less uniform or equal, the CwCw matrix becomes poorer conditioned. If the weights for the transition bands are zero, it is the smallest eigenvalues of CwCw that control the actual amplitude response A(ω)A(ω) in the transition bands. This explains why numerical errors in solving Equation 45 show up primarily in the transition bands. It also suggests this effect can be reduced by allowing a small weight in the transition bands. Indeed, one can design long filters by using spline transition functions with a small weight which then allows different pass and stopband weights.

Matlab Programs

A Matlab program for designing multiband FIR filters using the weighted least squares approximation described in Equation 55 and Equation 56 above is included in the Appendix. The program assumes passbands and stopbands that alternate with transition bands. The passbands are assumed to have constant gain for simplicity but that could be generalized with the results from FIR Digital Filters. The first for loop constructs the h^w(n)h^w(n) in Equation 47 by sequentially designing weighted bandpass filters and a separately weighted transition band similar to the example in Equation 57. These are added together in this loop to give the vector hwhw in Equation 45. The second for loop constructs the CC matrix in Equation 45 using the formula Equation 55. Care must be taken to correctly calculate the indeterminate values of sinc(0) and to properly include the effects of ωF=πωF=π.

When one uses zero weights in the transition bands, the CC matrix becomes ill conditioned when the product of the filter length NN and the sum of the transition band widths in Hertz is much above 12. This is an approximate rule which is somewhat affected by different passband and stopband weights, but it gives an indication of when numerical problems will occur. To reduce this problem, the program includes optimal spline transition functions so that a small weight can be used to improve the conditioning of CC with minimal effect on the optimality of the approximation.

Examples

To illustrate the effects of using a weighted least squared error design criterion, a simple length-21 linear phase lowpass FIR filter was designed, with unity weighting in the pass and stop bands and zero weighting in the transition band. The frequency response is shown in (Reference)a. The same filter is designed with a weight of 100 in the passband and the response is shown in (Reference)b and the case for a weight of 100 in the stopband is shown in (Reference)c. It is instructive to design many example filters and observe the effects of different weights, use of spline vs zero weight transition bands, and the effects of the transition band width on the pass and stopband performance.

Figure 8: Frequency Response of Length-21 FIR Filter Designed by Weighted Least Squared Error
This figure consist three graphs. The first graph is labeled Length-21 Lowpass FIR Filter Response, Least Squared error. It consist of a line extending to the right of the y axis and then intersecting another line that rises from the x axis. The lines along with the x and y axis form a rectangle. A wave form begins at the top of the rectangle. The wave wavers a little along the top of the rectangle and then falls to the x axis where is wavers as it makes it way to the left of the graph. The second graph is labeled Length 21 Weighted Passband Lowpass FIR Filter Response. The graph is essentially the same as the previous graph except for the absence of the rectangle and the waver of the wave along the x axis is much more pronounce. The third graph is labeled Length 21 Weighted Stopband Lowpass FIR Filter Response. This graph is similar to the previous to, but with this graph the waver of the initial segment is more pronounced and the waver of the section along the x axis is nearly non existent.

The same specifications that were used in the design using optimal spline transition functions of Figure 7a is used in the design here with unity weights in the stop and passbands and zero weights in the transition bands. The result is fairly similar to the spline function design and is shown in Figure 9a for N=51N=51. For lengths above around 131, numerical errors resulted in erratic performance in the transition bands. Indeed for this case, the use of the spline method would probably be superior. The advantage of the weighted least squares method is illustrated in Figure 9b where the same specifications are used but with a weight of 100 in the first passband, in Figure 9c where a weight of 100 is used in the second passband, and in Figure 9d where a weight of 100 is used in the third passband. This use of weights is impossible using the spline method or any windowing method.

Figure 9: Frequency Response of Length-15 FIR Filter Designed by Least Squared Error
This figure consist of four graphs. The top left graph is labeled Length-51 Multiband FIR Filter. It consist of a waveform that starts along the x axis with a slight waver and the shoots up almost vertically to about y=.7 where the wave wavers until it drops slightly to about y=.5 wavers a little and falls back to the x axis wavers a little and shoots all the back up to about y=1. The second graph labeled Length-51 Weighted Multiband FIR Filter progresses exactly the same way as the previous graph wave, but the first upper wavering area is actually smooth and flat in this graph. There is an arrow pointing to this section labeling it as weighted. The bottom left graph is labeled Length-51 Weighted Multiband FIR Filter and it also progresses the same way as the previous two graph except the second area of wavering is smoothed and flat and labeled with and arrow as weighted. The final graph the bottom right is the same except the last and highest portion is now flat and smooth and marked as weighted with an arrow.

Section Conclusions

This section has derived the four basic ideal lowpass filters: the constant gain passband lowpass filter, the linearly increasing gain passband lowpass filter, the differentiator with a lowpass filter, and the Hilbert transformer with a lowpass filter. It is shown that each of these can be modified to allow a spline transition function by a simple weighting function.

Because of using an L2L2 approximation error criterion and because of the orthogonality of the basis functions of the Fourier transform, it is shown that an optimal multiband filter can be built from the linear combination of these optimal building blocks. This new filter design method has the flexibility of the Parks-McClellan algorithm but the simplicity of the windowing methods. It is extremely fast and has no numerical problems. Unlike the windowing methods, the new method allows explicit independent control of multiple transition band edges and gives an optimal design. Its only limitation is not allowing error weighting.

We then derived a second method that likewise allowed multiple pass, stop, and transition bands with arbitrary band edges, but also allowed independent weighting of each frequency band. There are two limitations on this method. For long filters with wide transition bands with zero weights and where N(fp-fs)>12N(fp-fs)>12, the equations that must be solved are ill conditioned. This can be partially addressed using optimal spline functions with small weights in the transition bands. The second problem is that solving a large number of simultaneous equations can be slow and require considerable memory. These problems might be addressed by using special Toeplitz or Toeplitz plus Hankel algorithms [12] or some iterative method.

When should these methods be used? The second method which minimizes the weighted integral squared error should be used anytime the original problem dictates a squared error criterion and the product of the length and transition band width is less than twelve, N(fp-fs)<12N(fp-fs)<12. These conditions are often met because the squared error is a measure of the signal or noise energy and one seldom wants a long filter and a wide transition band. Even though this method requires solution of a set of simultaneous equations and is, therefore, slower than the spline transition function method, it executes in a few seconds on a PC or workstation and allows independent weighting of different frequency bands.

The first method which uses spline functions in the ideal response transition bands will design essentially arbitrarily long filters very quickly but it will not allow any error weighting. Although artificial transition functions are used in the ideal response, the optimized spline functions are very close to the response actually obtained by the second method with zero weighting in the transition band. This means the optimal approximation to the ideal response with spline functions transition bands is close to that obtained using the second numerical method. Comparisons of these effects for a single band can be found in [4]. If a Chebyshev approximation is desired, the Parks-McClellan method should be used although it too has numerical problems for long filters with wide transition bands. If different error measures are wanted in different bands, the iterative reweighted least squares (IRLS) algorithm [3] should be used. Recent research suggest that for many practical signal specifications, a mixture of Chebyshev and least squares is appropriate with no explicit transition bands [18].

If the equations that must be solved to obtain the optimal filter coefficients are ill-conditioned, an orthogonalization procedure can be used to improve the conditioning [13].

Characteristics of Optimal Filters

Gibbs phenomenon, transition band, pole-zero plots, etc.

Complex and Minimum Phase Approximation

Here we talk about which methods also solve the complex approximation problem. We talk about the minimum phase filter.

References

  1. Aaron, M. Robert. (1956, December). The Use of Least Squares in System Design. IRE Transaction on Circuit Theory, 3(4), 224–231.
  2. Algazi, V. Ralph and Suk, Minsoo. (1975, December). On the Frequency Weighted Least-Square Design of Finite Duration Filters. IEEE Transactions on Circuits and Systems, 22(12), 943–953.
  3. Burrus, C. S. and Barreto, J. A. and Selesnick, I. W. (1994, November). Iterative Reweighted Least Squares Design of FIR Filters. IEEE Transactions on Signal Processing, 42(11), 2926–2936.
  4. Burrus, C. S. and Soewito, A. W. and Gopinath, R. A. (1992, June). Least Squared Error FIR Filter Design with Transition Bands. IEEE Transactions on Signal Processing, 40(6), 1327–1340.
  5. Burrus, C. S. (1995, February). Multiband Least Squares FIR Filter Design. IEEE Transactions on Signal Processing, 43(2), 412–421.
  6. Fleischer, P. E. (1966, March). Digital Realization of Complex Transfer Functions. Simulation, 6, 171–180.
  7. Farden, David C. and Scharf, Louis L. (1974, June). Statistical Design of Nonrecursive Digital Filters. IEEE Transactions on ASSP, 22(3), 188–196.
  8. Guillemin, Ernst A. (1954, March). What is Nature's Error Criterion? [Follow-up articles in September 1954 issue.]. IRE Transaction on Circuit Theory, 1, 76.
  9. Lawson, C. L. and Hanson, R. J. (1974). Solving Least Squares Problems. Inglewood Cliffs, NJ: Prentice-Hall.
  10. Luenberger, D. G. (1969). Optimization by Vector Space Methods. New York: John Wiley & Sons.
  11. Moler, Cleve and Little, John and Bangert, Steve. (1989). Matlab User's Guide. South Natick, MA: The MathWorks, Inc.
  12. Merchant, G. A. and Parks, T. W. (1982, February). Efficient Solution of a Toeplitz-Plus-Hankel Coefficient Matrix System of Equations. IEEE Transactions on ASSP, 30(1), 40–44.
  13. Okuda, Masahiro and Ikehara, Masaaki and ichi Takahashi, Shin. (to appear 1998). Fast and Stable Least–Squares Approach for Linear Phase FIR Filters. IEEE Transactions on Signal Processing, 46,
  14. Oetken, G. and Parks, T. W. and Schüssler, H. W. (1975, June). New Results in the Design of Digital Interpolators. [also in Selected Papers in DSP II, IEEE Press, 1975, page 167.]. IEEE Trans. on ASSP, ASSP-23, 301–309.
  15. Oppenheim, A. V. and Schafer, R. W. (1989). Discrete-Time Signal Processing. Englewood Cliffs, NJ: Prentice-Hall.
  16. Parks, T. W. and Burrus, C. S. (1987). Digital Filter Design. New York: John Wiley & Sons.
  17. Schüssler, H. W. (1978, 1988). Digitale Signalverarbeitung II. [in German]. Universität Erlangen–Nürnberg, Erlangen, Germany: Lehrstuhl für Nachrichtentechnik.
  18. Selesnick, Ivan W. and Lang, Markus and Burrus, C. Sidney. (1996, August). Constrained Least Square Design of FIR Filters without Explicitly Specified Transition Bands. IEEE Transactions on Signal Processing, 44(8), 1879–1892.
  19. Schüssler, H. W. and Steffen, P. (1984). A Hybrid System for the Reconstruction of a Smooth Function from its Samples. Circuits, Systems, and Signal Processing, 3(3), 295–314.
  20. Strang, Gilbert. (1976). Linear Algebra and Its Applications. New York: Academic Press.
  21. Tufts, Donald W. and Rorabacher, Darold W. and Mosier, Wilbur E. (1970, June). Designing Simple, Effective Digital Filters. IEEE Transactions on Audio and Electroacoustics, 18(2), 142–158.
  22. Weisburn, B. A. and Parks, T. W. and Shenoy, B. G. (1994, April 19–22). Error Criteria for Filter Design. In Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing. (Vol. 3, p. III:565–568). IEEE ICASSP-94, Adelaide, Australia

Collection Navigation

Content actions

Download:

Collection as:

PDF | EPUB (?)

What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

Downloading to a reading device

For detailed instructions on how to download this content's EPUB to your specific device, click the "(?)" link.

| More downloads ...

Module as:

PDF | More downloads ...

Add:

Collection to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

Module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks