# Connexions

You are here: Home » Content » Lab 5b - Digital Filter Design (part 2)

### 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?

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

#### 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

This module is included inLens: NSF Partnership in Signal Processing
By: Sidney BurrusAs a part of collection: "Purdue Digital Signal Processing Labs (ECE 438)"

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

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

• Featured Content

This module is included inLens: Connexions Featured Content
By: ConnexionsAs a part of collection: "Purdue Digital Signal Processing Labs (ECE 438)"

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

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

#### Also in these lenses

• UniqU content

This module is included inLens: UniqU's lens
By: UniqU, LLCAs a part of collection: "Purdue Digital Signal Processing Labs (ECE 438)"

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

• Lens for Engineering

This module is 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.

# Lab 5b - Digital Filter Design (part 2)

Module by: Charles A. Bouman. E-mail the author

Questions or comments concerning this laboratory should be directed to Prof. Charles A. Bouman, School of Electrical and Computer Engineering, Purdue University, West Lafayette IN 47907; (765) 494-0340; bouman@ecn.purdue.edu

## Introduction

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.

## Filter Design Using Standard Windows

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

h ( n ) = w ( n ) h i d e a l ( n ) , h ( n ) = w ( n ) h i d e a l ( n ) ,
(1)

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

h i d e a l ( n ) = ω c π sinc ω c π n - N - 1 2 , h i d e a l ( n ) = ω c π sinc ω c π n - N - 1 2 ,
(2)

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

The rectangular window is defined as

w ( n ) = 1 n = 0 , 1 , ... , N - 1 0 otherwise w ( n ) = 1 n = 0 , 1 , ... , N - 1 0 otherwise
(3)

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.

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.

1. Hanning window (as defined in Matlab, command hann(N)):
w(n)=0.5-0.5cos2πnN-1n=0,1,...,N-10otherwisew(n)=0.5-0.5cos2πnN-1n=0,1,...,N-10otherwise
(4)
2. Hamming window
w(n)=0.54-0.46cos2πnN-1n=0,1,...,N-10otherwisew(n)=0.54-0.46cos2πnN-1n=0,1,...,N-10otherwise
(5)
3. Blackman window
w(n)=0.42-0.5cos2πnN-1+0.08cos4πnN-1n=0,1,...,N-10otherwisew(n)=0.42-0.5cos2πnN-1+0.08cos4πnN-1n=0,1,...,N-10otherwise
(6)

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.

 Window (length N) Mainlobe width Peak-to-sidelobe amplitude (dB) R e c t a n g u l a r R e c t a n g u l a r 4 π / N 4 π / N - 13 d B - 13 d B H a n n i n g H a n n i n g 8 π / N 8 π / N - 32 d B - 32 d B H a m m i n g H a m m i n g 8 π / N 8 π / N - 43 d B - 43 d B B l a c k m a n B l a c k m a n 12 π / N 12 π / N - 58 d B - 58 d B

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.

### Note:

Use at least 512 sample points in computing the DTFT by typing the command DTFT(window,512). Type help DTFT for further information on this function.

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

1. Submit the figure containing the time domain plots of the four windows.
2. Submit the figure containing the DTFT (in decibels) of the four windows.
3. Submit the table of the measured and theoretical window spectrum parameters. Comment on how close the experimental results matched the ideal values. Also comment on the relation between the width of the mainlobe and the peak-to-sidelobe amplitude.
4. Submit the plots of your designed filter's impulse response and the magnitude of the filter's DTFT.

## Filter Design Using the Kaiser Window

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:

δ = min { δ p , δ s } δ = min { δ p , δ s }
(7)

The Kaiser window function of length NN is given by

w ( n ) = I 0 β n ( N - 1 - n ) N - 1 I 0 ( β ) n = 0 , 1 , ... , N - 1 0 otherwise w ( n ) = I 0 β n ( N - 1 - n ) N - 1 I 0 ( β ) n = 0 , 1 , ... , N - 1 0 otherwise
(8)

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:

β = 0 . 1102 ( A - 8 . 7 ) A > 50 0 . 5842 ( A - 21 ) 0 . 4 + 0 . 07886 ( A - 21 ) 21 A 50 0 . 0 A < 21 β = 0 . 1102 ( A - 8 . 7 ) A > 50 0 . 5842 ( A - 21 ) 0 . 4 + 0 . 07886 ( A - 21 ) 21 A 50 0 . 0 A < 21
(9)
N = 1 + A - 8 2 . 285 ( ω s - ω p ) N = 1 + A - 8 2 . 285 ( ω s - ω p )
(10)

where ·· is the ceiling function, i.e. xx 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 ββ:

• β = 0 β = 0
(11)
• β = 1 β = 1
(12)
• β = 5 β = 5
(13)

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

### Hint:

To create the Kaiser windows, use the Matlab command kaiser(N,beta) command where N is the length of the filter and beta is the shape parameter ββ. To insure at least 512 points in the plot use the command DTFT(window,512) when computing the DTFT.

### INLAB REPORT:

Submit the plots of the 3 Kaiser windows and the magnitude of their DTFT's in decibels. Comment on how the value ββ affects the shape of the window and the sidelobes of the DTFT.

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:

ω p = 1 . 8 ω c = 2 . 0 ω s = 2 . 2 δ p = 0 . 05 δ s = 0 . 005 ω p = 1 . 8 ω c = 2 . 0 ω s = 2 . 2 δ p = 0 . 05 δ s = 0 . 005
(14)

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

ω c = ω p + ω s 2 ω c = ω p + ω s 2
(15)

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.

### Hint:

To compute the passband ripple, find the value of the DTFT at frequencies corresponding to the passband using the command H(abs(w)<=1.8) where H is the DTFT of h(n)h(n) and w is the corresponding vector of frequencies. Then use this vector to compute the passband ripple. Use a similar procedure for the stopband ripple.

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:

1. Submit the values of ββ and NN that you computed.
2. Submit the three plots of the filter's magnitude response. Make sure the plots are labeled.
3. Submit the values of the passband and stopband ripple. Does this filter meet the design specifications?
4. Submit the magnitude plot of the DTFT in dB for the filtered signal. Compare this plot to the plot of Figure 4.
5. Comment on how the frequency content and the audio quality of the filtered signal have changed after filtering.

## FIR Filter Design Using Parks-McClellan Algorithm

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,

• Kaiser filters are not guaranteed to be the minimum length filter which meets the design constraints.
• Kaiser filters do not allow passband and stopband ripple to be varied independently.

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

y ( n ) = k = 0 N - 1 x ( n - k ) h ( k ) . y ( n ) = k = 0 N - 1 x ( n - k ) h ( k ) .
(16)

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

y ( n ) = k = 0 N / 2 - 1 x ( n - k ) + x ( n - N + 1 + k ) h ( k ) . y ( n ) = k = 0 N / 2 - 1 x ( n - k ) + x ( n - N + 1 + k ) h ( k ) .
(17)

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:

• f - vector containing an even number of band edge frequencies. For a simple low pass filter, f=[wp,ws] , where wp and ws are the passband and stopband frequencies, respectively.
• m - vector containing the ideal filter magnitudes of the filter in each band. For a simple low pass filter m=[1,0] .
• ripple - vector containing the allowed ripple in each band. For a simple low pass filter ripple=[delta_p,delta_s] , where delta_p and delta_s are the passband and stopband ripples, respectively.
• 2*pi - value, in radians, that corresponds to the sampling frequency.

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

H ( z ) = b ( 1 ) + b ( 2 ) z - 1 + + b ( n + 1 ) z - n H ( z ) = b ( 1 ) + b ( 2 ) z - 1 + + b ( n + 1 ) z - n
(18)

(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:

1. Submit the final measured values of filter length, passband ripple, and stopband ripple. How accurate was the filter order computation using Matlab's firpmord? How does the length of this filter compare to the filter designed using a Kaiser window?
2. Submit the plot of the filter's DTFT. How does the frequency response of the Parks-McClellan filter compare to the filter designed using the Kaiser window? Comment on the shape of both the passband and stopband.

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.

### INLAB REPORT:

Submit the plot of the DTFT magnitude for the filtered signal. Comment on how the audio quality of the signal changes after filtering. Also comment on any differences in audio quality between the Parks-McClellan filtered speech and the Kaiser filtered speech.

## Design of Discrete-Time IIR Filters Using Numerical Optimization

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

H θ ( z ) = θ 1 + θ 2 z - 1 + θ 3 z - 2 1 + θ 4 z - 1 + θ 5 z - 2 . H θ ( z ) = θ 1 + θ 2 z - 1 + θ 3 z - 2 1 + θ 4 z - 1 + θ 5 z - 2 .
(19)

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

H s h ( e j ω ) = sinc ω 2 π for | ω | < π H s h ( e j ω ) = sinc ω 2 π for | ω | < π
(20)

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ω):

Cost ( θ ) = - π π 1 H s h ( e j ω ) - H θ ( e j ω ) 2 d ω Cost ( θ ) = - π π 1 H s h ( e j ω ) - H θ ( e j ω ) 2 d ω
(21)

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:

Err d B ( ω ) = 20 log 10 H s h ( e j ω ) | H θ * ( e j ω ) | Err d B ( ω ) = 20 log 10 H s h ( e j ω ) | H θ * ( e j ω ) |
(22)

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:

• Write a Matlab function prefilter(w,theta) which computes the frequency response Hθ(ejω)Hθ(ejω) from equation Equation 19 for the vector of input frequencies w and the parameter vector theta.
• Write a Matlab function Cost(theta) which computes the total squared error of equation Equation 21. Use a sampling interval Δω=0.01Δω=0.01 for the functions Hθ(ejω)Hθ(ejω) and 1/Hsh(ejω)1/Hsh(ejω).
• Use the command fminsearch from Matlab's Optimization Toolbox to compute the value of the parameter θθ which minimizes Cost(theta). The function fminsearch has the syntax X = fminsearch('function_name',initial_value) where function_name is the name of the function being minimized (Cost), initial_value is the starting value for the unknown parameter, and X is the minimizing parameter vector. Choose an initial value of (θ1,θ2,θ3,θ4,θ5)=(1,0,0,0,0)(θ1,θ2,θ3,θ4,θ5)=(1,0,0,0,0) so that Hθ(ejω)=1Hθ(ejω)=1.
• Use the subplot command to plot the following three functions on the interval [-π,π][-π,π].
• The desired filter magnitude response 1/Hsh(ejω)1/Hsh(ejω).
• The designed IIR filter magnitude response |Hθ*(ejω)||Hθ*(ejω)|.
• The error in decibels, from equation Equation 22.

### INLAB REPORT

Do the following:

1. Submit the printouts of the code for the two Matlab functions prefilter.m and Cost.m.
2. Give an analytical expression for the optimized transfer function Hθ*(z)Hθ*(z) with the coefficients that were computed.
3. Submit the three plots. On the error plot, mark the frequency ranges where the approximation error is high.

## Footnotes

1. The advantages of using such symmetries varies considerably with the implementation and application. On many modern computing architectures the computational cost of adds and multiplies are similar, and the overhead of control loops may eliminate the advantages of reduced operations.

## References

1. Proakis, J. G., Manolakis, D. G. (1996). Digital Signal Processing. (3rd). New York: Prentice-Hall.

## Content actions

PDF | EPUB (?)

### What is an EPUB file?

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

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?

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