This is the second part of a two week laboratory in digital filter design. The first week of the laboratory covered some basic examples of FIR and IIR filters, and then introduced the concepts of filter design. In this week we will cover more systematic methods of designing both FIR and IIR filters.

Download DTFT.m for the following section.

We can generalize the idea of truncation by using different windowing functions to truncate an ideal filter's impulse response. Note that by simply truncating the ideal filter's impulse response, we are actually multiplying (or “windowing”) the impulse response by a shifted rect()rect() function. This particular type of window is called a rectangular window. In general, the impulse reponse h(n)h(n) of the designed filter is related to the impulse response hideal(n)hideal(n) of the ideal filter by the relation

where w(n)w(n) is an NN-point window. We assume that

where ωcωc is the cutoff frequency and NN is the desired window length.

The rectangular window is defined as

The DTFT of w(n)w(n) for N=21N=21 is shown in Figure 1. The rectangular window is usually not preferred because it leads to the large stopband and passband ripple as shown in Figure 2.

Figure 1: DTFT of a rectangular window of length 21.

Figure 2: Frequency response of low-pass filter, designed using the truncation method.

More desirable frequency characteristics can be obtained by making a better selection for the window, w(n)w(n). In fact, a variety of raised cosine windows are widely used for this purpose. Some popular windows are listed below.

In filter design using different truncation windows, there are two key frequency domain effects that are important to the design: the transition band roll-off, and the passband and stopband ripple (see Figure 3 below). There are two corresponding parameters in the spectrum of each type of window that influence these filter parameters. The filter's roll-off is related to the width of center lobe of the window's magnitude spectrum. The ripple is influenced by the ratio of the mainlobe amplitude to the first sidelobe amplitude (or difference if using a dB scale). These two window spectrum parameters are not independent, and you should see a trend as you examine the spectra for different windows. The theoretical values for the mainlobe width and the peak-to-sidelobe amplitude are shown in Table 1.

Plot the rectangular, Hamming, Hanning, and Blackman window functions of length 21 on a single figure using the subplot command. You may use the Matlab commands hamming, hann, and blackman. Then compute and plot the DTFT magnitude of each of the four windows. Plot the magnitudes on a decibel scale (i.e., plot 20log10|W(ejω)|20log10|W(ejω)|). Download and use the function DTFT.m to compute the DTFT.

Measure the null-to-null mainlobe width (in rad/sample) and the peak-to-sidelobe amplitude (in dB) from the logarithmic magnitude response plot for each window type. The Matlab command zoom is helpful for this. Make a table with these values and the theoretical ones.

Now use a Hamming window to design a lowpass filter h(n)h(n) with a cutoff frequency of ωc=2.0ωc=2.0 and length 21. Note: You need to use Equation 1 and Equation 2 for this design. In the same figure, plot the filter's impulse response, and the magnitude of the filter's DTFT in decibels.

INLAB REPORT

Download nspeech2.mat for the following section.

Figure 3: Tolerance specifications for the frequency response of a filter.

The standard windows of the "Filter Design Using Standard Windows" section are an improvement over simple truncation, but these windows still do not allow for arbitrary choices of transition bandwidth and ripple. In 1964, James Kaiser derived a family of near-optimal windows that can be used to design filters which meet or exceed any filter specification. The Kaiser window depends on two parameters: the window length NN, and a parameter ββ which controls the shape of the window. Large values of ββ reduce the window sidelobes and therefore result in reduced passband and stopband ripple. The only restriction in the Kaiser filter design method is that the passband and stopband ripple must be equal in magnitude. Therefore, the Kaiser filter must be designed to meet the smaller of the two ripple constraints:

The Kaiser window function of length NN is given by

where I0(·)I0(·) is the zero'th order modified Bessel function of the first kind, NN is the length of the window, and ββ is the shape parameter.

Kaiser found that values of ββ and NN could be chosen to meet any set of design parameters, (δ,ωp,ωs)(δ,ωp,ωs), by defining A=-20log10δA=-20log10δ and using the following two equations:

where ⌈·⌉⌈·⌉ is the ceiling function, i.e. ⌈x⌉⌈x⌉ is the smallest integer which is greater than or equal to xx.

To further investigate the Kaiser window, plot the Kaiser windows and their DTFT magnitudes (in dB) for N=21N=21 and the following values of ββ:

For each case use at least 512 points in the plot of the DTFT.

Figure 4: DTFT of a section of noisy speech.

Next use a Kaiser window to design a low pass filter, h(n)h(n), to remove the noise from the signal in nspeech2.mat using equations Equation 1 and Equation 2. To do this, use equations Equation 9 and Equation 10 to compute the values of ββ and NN that will yield the following design specifications:

The low pass filter designed with the Kaiser method will automatically have a cut-off frequency centered between ωpωp and ωsωs.

Plot the magnitude of the DTFT of h(n)h(n) for |ω|<π|ω|<π . Create three plots in the same figure: one that shows the entire frequency response, and ones that zoom in on the passband and stopband ripple, respectively. Mark ωpωp, ωsωs, δpδp, and δsδs on these plots where appropriate. Note: Since the ripple is measured on a magnitude scale, DO NOT use a decibel scale on this set of plots.

From the Matlab prompt, compute the stopband and passband ripple (do not do this graphically). Record the stopband and passband ripple to three decimal places.

Filter the noisy speech signal in nspeech2.mat using the filter you have designed. Then compute the DTFT of 400 samples of the filtered signal starting at time n=20000n=20000 (i.e. 20001:20400 ). Plot the magnitude of the DTFT samples in decibels versus frequency in radians for |ω|<π|ω|<π. Compare this with the spectrum of the noisy speech signal shown in Figure 4. Play the noisy and filtered speech signals back using sound and listen to them carefully.

INLAB REPORT

Do the following:

Click here for help on the firpm function for Parks-McClellan filter design. Download the data file nspeech2.mat for the following section.

Kaiser windows are versatile since they allow the design of arbitrary filters which meet specific design constraints. However, filters designed with Kaiser windows still have a number of disadvantages. For example,

Minimizing filter length is important because in many applications the length of the filter determines the amount of computation. For example, an FIR filter of length NN may be directly implemented in the time domain by evaluating the expression

For each output value y(n)y(n) this expression requires NN multiplies and N-1N-1 additions.

Oftentimes h(n)h(n) is a symmetric filter so that h(n)=h(N-1-n)h(n)=h(N-1-n). If the filter h(n)h(n) is symmetric and NN is even, then Equation 16 may be more efficiently computed as

This strategy reduces the computation to N/2N/2 multiplies and N-1N-1 adds for any value of NN1. Note that the computational effort is linearly proportional to the length of the filter.

The Kaiser filters do not guarantee the minimum possible filter length. Since the filter has equal passband and stopband ripple, it will usually exceed design requirements in one of the two bands; this results in an unnecessarily long filter. A better design would allow the stopband and passband constraints to be specified separately.

In 1972, Parks and McClellan devised a methodology for designing symmetric filters that minimize filter length for a particular set of design constraints {ωpωp, ωsωs, δpδp, δsδs}. The resulting filters minimize the maximum error between the desired frequency response and the actual frequency response by spreading the approximation error uniformly over each band. The Parks and McClellan algorithm makes use of the Remez exchange algorithm and Chebyshev approximation theory. Such filters that exhibit equiripple behavior in both the passband and the stopband, and are sometimes called equiripple filters.

As with Kaiser filters, designing a filter with the Parks and McClellan algorithm is a two step process. First the length (i.e. order) of the filter must be computed based on the design constraints. Then the optimal filter for a specified length can be determined. As with Kaiser windows, the filter length computation is approximate so the resulting filter may exceed or violate the design constraints. This is generally not a problem since the filter can be redesigned for different lengths until the constraints are just met.

The Matlab command for computing the approximate filter length is

[n,fo,mo,w] = firpmord(f,m,ripple,2*pi)

where the inputs are:

The outputs of the command are n = filter length - 1 , and the vectors fo , mo , and w which are intermediate filter parameters.

Once the filter length, n , is obtained, the Matlab command for designing a Parks-McClellan filter is b = firpm(n,fo,mo,w) . The inputs n , fo , mo , and w are the corresponding outputs of firpmord, and the output b is a vector of FIR filter coefficients such that

(What is the impulse response of this filter?)

For further information, read the help document on using Matlab to implement the Parks-McClellan algorithm.

Now design a symmetric FIR filter using firpmord and firpm in Matlab to meet the design specifications given in the "Filter Design Using the Kaiser Window" section. Compute the DTFT of the filter's response for at least 512 points, and use this result to compute the passband and stopband ripple of the filter that was designed. Adjust the filter length until the minimum order which meets the design constraints is found. Plot the magnitude of the DTFT in dB of the final filter design.

INLAB REPORT

Do the following:

Use the filter you have designed to remove the noise from the signal nspeech2.mat . Play the noisy and filtered speech signals back using sound and listen to them carefully. Compute the DTFT of 400 samples of the filtered signal starting at time n=20,001n=20,001 (i.e. 20001:20400 ). Plot the magnitude of the DTFT in decibels versus frequency in radians for |ω|<π|ω|<π. Compare this with the spectrum of the noisy speech signal shown in Figure 4, and also with the magnitude of the DTFT of the Kaiser filtered signal.

In this section, we consider the design of discrete-time IIR filters through the direct search of filter parameters that will minimize a specific design criterion. Such “brute force” approaches to filter design have become increasingly more popular due to the wide availability of high speed computers and robust numerical optimization methods.

Typically, numerical approaches to filter design have two parts. First, they design a cost, or error criterion. This criterion is a measure of the difference between the ideal filter response and the response of the computed or “approximate” filter. The goal is to find the approximate filter with the lowest cost (error). Mean square error is a popular cost criterion. The second part is to minimize the cost with respect to the filter parameters. We will perform the required numerical optimization with the fminsearch function in Matlab's Optimization Toolbox.

In order to formulate a cost criterion, we must first select a model for the discrete-time filter of interest. There are many ways of doing this, but we will use the coefficients of a rational transfer function to model (or parameterize) the set of second order IIR filters. In this case, the elements of the vector θ=[θ1,θ2,θ3,θ4,θ5]θ=[θ1,θ2,θ3,θ4,θ5] are the coefficients of the transfer function

Using this parameterization, we may then define a function Cost(θ)Cost(θ) which is the “cost” of using the filter Hθ(z)Hθ(z).

To illustrate this numerical optimization approach, we will design a digital filter that compensates for the roll-off due to the sample-and-hold process in an audio CD player. In lab 4, we saw that the sample-and-hold operation in a conventional D/A converter causes the reconstructed signal to be filtered by the function

One method of reducing this distortion is to digitally “pre-filter” a signal with the inverse transfer function, 1/Hsh1/Hsh. This filter 1/Hsh1/Hsh pre-distorts the audio signal so the reconstructed signal has the desired frequency response. We would like to approximate the filter 1/Hsh1/Hsh using the second order filter of Equation 19.

For an audio CD player, the magnitude of the frequency response is generally considered to be more important than the phase. This is because we are not perceptually sensitive to phase distortion in sound. Therefore, we may choose a cost function which computes the total squared error between the magnitudes of the desired pre-filter response, 1/Hsh(ejω)1/Hsh(ejω), and the second order filter Hθ(ejω)Hθ(ejω):

The θθ parameters that minimize this cost function will be the parameters of our designed filter. A more complex approach might also account for the filter phase, but for simplicity we will only try to match the filter magnitudes.

After the filter is designed, we may compute the difference between the CD player's frequency response in dB and the ideal desired response in dB that the CD player should have:

where θ*θ* is the optimized value of θθ and Hθ*(ejω)Hθ*(ejω) is the optimum second-order filter.

Do the following to perform this filter design:

INLAB REPORT

Do the following: