This is the first part of a two week laboratory in digital filter design. The first week of the laboratory covers some basic examples of FIR and IIR filters, and then introduces the concepts of FIR filter design. Then the second week covers systematic methods of designing both FIR and IIR filters.

In digital signal processing applications, it is often necessary to change the relative amplitudes of frequency components or remove undesired frequencies of a signal. This process is called filtering. Digital filters are used in a variety of applications. In Laboratory 4, we saw that digital filters may be used in systems that perform interpolation and decimation on discrete-time signals. Digital filters are also used in audio systems that allow the listener to adjust the bass (low-frequency energy) and the treble (high frequency energy) of audio signals.

Digital filter design requires the use of both frequency domain and time domain techniques. This is because filter design specifications are often given in the frequency domain, but filters are usually implemented in the time domain with a difference equation. Typically, frequency domain analysis is done using the Z-transform and the discrete-time Fourier Transform (DTFT).

In general, a linear and time-invariant causal digital filter with input x(n)x(n) and output y(n)y(n) may be specified by its difference equation

where bibi and akak are coefficients which parameterize the filter. This filter is said to have NN zeros and MM poles. Each new value of the output signal, y(n)y(n), is determined by past values of the output, and by present and past values of the input. The impulse response, h(n)h(n), is the response of the filter to an input of δ(n)δ(n), and is therefore the solution to the recursive difference equation

There are two general classes of digital filters: infinite impulse response (IIR) and finite impulse response (FIR). The FIR case occurs when ak=0ak=0, for all kk. Such a filter is said to have no poles, only zeros. In this case, the difference equation in Equation 2 becomes

Since Equation 3 is no longer recursive, the impulse response has finite duration NN.

In the case where ak≠0ak≠0, the difference equation usually represents an IIR filter. In this case, Equation 2 will usually generate an impulse response which has non-zero values as n→∞n→∞. However, later we will see that for certain values of ak≠0ak≠0 and bibi, it is possible to generate an FIR filter response.

The Z-transform is the major tool used for analyzing the frequency response of filters and their difference equations. The Z-transform of a discrete-time signal, x(n)x(n), is given by

where zz is a complex variable. The DTFT may be thought of as a special case of the Z-transform where zz is evaluated on the unit circle in the complex plane.

From the definition of the Z-transform, a change of variable m=n-Km=n-K shows that a delay of KK samples in the time domain is equivalent to multiplication by z-Kz-K in the Z-transform domain.

We may use this fact to re-write Equation 1 in the Z-transform domain, by taking Z-transforms of both sides of the equation:

From this formula, we see that any filter which can be represented by a linear difference equation with constant coefficients has a rational transfer function (i.e. a transfer function which is a ratio of polynomials). From this result, we may compute the frequency response of the filter by evaluating H(z)H(z) on the unit circle:

There are many different methods for implementing a general recursive difference equation of the form Equation 1. Depending on the application, some methods may be more robust to quantization error, require fewer multiplies or adds, or require less memory. Figure 1 shows a system diagram known as the direct form implementation; it works for any discrete-time filter described by the difference equation in Equation 1. Note that the boxes containing the symbol z-1z-1 represent unit delays, while a parameter written next to a signal path represents multiplication by that parameter.

Figure 1: Direct form implementation for a discrete-time filter described by a linear recursive difference equation.

Download the files, nspeech1.mat and DTFT.m for the following section.

Figure 2: Location of two zeros for simple a FIR filter.

To illustrate the use of zeros in filter design, you will design a simple second order FIR filter with the two zeros on the unit circle as shown in Figure 2. In order for the filter's impulse response to be real-valued, the two zeros must be complex conjugates of one another:

where θθ is the angle of z1z1 relative to the positive real axis. We will see later that θ∈[0,π]θ∈[0,π] may be interpreted as the location of the zeros in the frequency response.

The transfer function for this filter is given by

Use this transfer function to determine the difference equation for this filter. Then draw the corresponding system diagram and compute the filter's impulse response h(n)h(n).

This filter is an FIR filter because it has impulse response h(n)h(n) of finite duration. Any filter with only zeros and no poles other than those at 0 and ±∞±∞ is an FIR filter. Zeros in the transfer function represent frequencies that are not passed through the filter. This can be useful for removing unwanted frequencies in a signal. The fact that Hf(z)Hf(z) has zeros at e±jθe±jθ implies that Hf(e±jθ)=0Hf(e±jθ)=0. This means that the filter will not pass pure sine waves at a frequency of ω=θω=θ.

Use Matlab to compute and plot the magnitude of the filter's frequency response |Hf(ejω)||Hf(ejω)| as a function of ωω on the interval -π<ω<π-π<ω<π, for the following three values of θθ:

Put all three plots on the same figure using the subplot command.

In the next experiment, we will use the filter Hf(z)Hf(z) to remove an undesirable sinusoidal interference from a speech signal. To run the experiment, first download the audio signal nspeech1.mat, and the M-file DTFT.m Load nspeech1.mat into Matlab using the command load nspeech1. This will load the signal nspeech1 into the workspace. Play nspeech1 using the sound command, and then plot 101 samples of the signal for the time indices (100:200).

We will next use the DTFT command to compute samples of the DTFT of the audio signal. The DTFT command has the syntax

[X,w]=DTFT(x,M)

where x is a signal which is assumed to start at time n=0n=0, and MM specifies the number of output points of the DTFT. The command [X,w]=DTFT(x,0) will generate a DTFT that is the same duration as the input; if this is not sufficient, it may be increased by specifying M . The outputs of the function are a vector X containing the samples of the DTFT, and a vector w containing the corresponding frequencies of these samples.

Compute the magnitude of the DTFT of 1001 samples of the audio signal for the time indices (100:1100). Plot the magnitude of the DTFT samples versus frequency for |ω|<π|ω|<π. Notice that there are two large peaks corresponding to the sinusoidal interference signal. Use the Matlab max command to determine the exact frequency of the peaks. This will be the value of θθ that we will use for filtering with Hf(z)Hf(z).

Write a Matlab function FIRfilter(x) that implements the filter Hf(z)Hf(z) with the measured value of θθ and outputs the filtered signal (Hint: Use convolution). Apply the new function FIRfilter to the nspeech1 vector to attenuate the sinusoidal interference. Listen to the filtered signal to hear the effects of the filter. Plot 101 samples of the signal for the time indices (100:200), and plot the magnitude of the DTFT of 1001 samples of the filtered signal for the time indices (100:1100).

INLAB REPORT

For both the original audio signal and the filtered output, hand in the following:

Download the file pcm.mat for the following section.

Figure 3: Location of two poles for a simple IIR filter.

While zeros attenuate a filtered signal, poles amplify signals that are near their frequency. In this section, we will illustrate how poles can be used to separate a narrow-band signal from adjacent noise. Such filters are commonly used to separate modulated signals from background noise in applications such as radio-frequency demodulation.

Consider the following transfer function for a second order IIR filter with complex-conjugate poles:

Figure 3 shows the locations of the two poles of this filter. The poles have the form

where rr is the distance from the origin, and θθ is the angle of p1p1 relative to the positive real axis. From the theory of Z-transforms, we know that a causal filter is stable if and only if its poles are located within the unit circle. This implies that this filter is stable if and only if |r|<1|r|<1. However, we will see that by locating the poles close to the unit circle, the filter's bandwidth may be made extremely narrow around θθ.

This two-pole system is an example of an IIR filter because its impulse response has infinite duration. Any filter with nontrivial poles (not at z = 0 or ±∞±∞) is an IIR filter unless the poles are canceled by zeros.

Calculate the magnitude of the filter's frequency response |Hi(ejw)||Hi(ejw)| on |ω|<π|ω|<π for θ=π/3θ=π/3 and the following three values of rr.

Put all three plots on the same figure using the subplot command.

In the following experiment, we will use the filter Hi(z)Hi(z) to separate a modulated sinusoid from background noise. To run the experiment, first download the file pcm.mat and load it into the Matlab workspace using the command load pcm . Play pcm using the sound command. Plot 101 samples of the signal for indices (100:200), and then compute the magnitude of the DTFT of 1001 samples of pcm using the time indices (100:1100). Plot the magnitude of the DTFT samples versus radial frequency for |ω|<π|ω|<π. The two peaks in the spectrum correspond to the center frequency of the modulated signal. The low amplitude wideband content is the background noise. In this exercise, you will use the IIR filter described above to amplify the desired signal, relative to the background noise.

The pcm signal is modulated at 3146Hz and sampled at 8kHz. Use these values to calculate the value of θθ for the filter Hi(z)Hi(z). Remember from the sampling theorem that a radial frequency of 2π2π corresponds to the sampling frequency.

Write a Matlab function IIRfilter(x) that implements the filter Hi(z)Hi(z). In this case, you need to use a for loop to implement the recursive difference equation. Use your calculated value of θθ and r=0.995r=0.995. You can assume that y(n) is equal to 0 for negative values of nn. Apply the new command IIRfilter to the signal pcm to separate the desired signal from the background noise, and listen to the filtered signal to hear the effects. Plot the filtered signal for indices (100:200), and then compute the DTFT of 1001 samples of the filtered signal using the time indices (100:1100). Plot the magnitude of this DTFT. In order to see the DTFT around ω=θω=θ more clearly, plot also the portion of this DTFT for the values of ωω in the range [θ-0.02,θ+0.02][θ-0.02,θ+0.02]. Use your calculated value of θθ.

INLAB REPORT

For both the pcm signal and the filtered output, submit the following:

Download the file nspeech2.mat for the following section.

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

Oftentimes it is necessary to design a good approximation to an ideal lowpass, highpass or bandpass filter. Figure 4 illustrates the typical characteristics of a real low-pass filter. The frequencies |ω|<ωp|ω|<ωp are known as the passband, and the frequencies ωs<|ω|≤πωs<|ω|≤π are the stopband. For any real filter, ωp<ωsωp<ωs. The range of frequencies ωp≤ω≤ωsωp≤ω≤ωs is known as the transition band. The magnitude of the filter response, H(ejω)H(ejω), is constrained in the passband and stopband by the following two equations

where δpδp and δsδs are known as the passband and stopband ripple respectively. Most lowpass filter design techniques depend on the specification of these four parameters: ωpωp, ωsωs, δpδp, and δsδs.

Figure 5: DTFT of a section of noisy speech.

To illustrate the selection of these parameters consider the problem of filtering out background noise from a speech signal. Figure 5 shows the magnitude of the DTFT over a window of such a signal, called nspeech2. Notice that there are two main components in nspeech2: one at the low frequencies and one at the high. The high frequency signal is noise, and it is band limited to |ω|>2.2|ω|>2.2. The low frequency signal is speech and it is band limited to |ω|<1.8|ω|<1.8. Download the file nspeech2.mat. and load it into the Matlab workspace. It contains the signal nspeech2 from Figure 5. Play the nspeech2 using the sound command and note the quality of the speech and background noise.

In the following sections, we will compute low-pass filters for separating the speech and noise using a number of different methods.

Filter Design Using Truncation

Ideally, a low-pass filter with cutoff frequency ωcωc should have a frequency response of

H i d e a l ( e j w ) = 1 | ω | ≤ ω c 0 ω c < | ω | ≤ π H i d e a l ( e j w ) = 1 | ω | ≤ ω c 0 ω c < | ω | ≤ π

(21)

and a corresponding impulse response of

h i d e a l ( n ) = ω c π sinc ( ω c n π ) for - ∞ < n < ∞ h i d e a l ( n ) = ω c π sinc ( ω c n π ) for - ∞ < n < ∞

(22)

However, no real filter can have this frequency response because hideal(n)hideal(n) is infinite in duration.

One method for creating a realizable approximation to an ideal filter is to truncate this impulse response outside of n∈[-M,M]n∈[-M,M].

h t r u n c ( n ) = ω c π sinc ( ω c π n ) n = - M , ... , 0 , 1 , ... , M 0 otherwise h t r u n c ( n ) = ω c π sinc ( ω c π n ) n = - M , ... , 0 , 1 , ... , M 0 otherwise

(23)

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

Figure 6 shows the magnitude response of the lowpass filter with cutoff frequency ωc=2.0ωc=2.0, with the impulse response truncated to n∈[-10,10]n∈[-10,10]. Notice the oscillatory behavior of the magnitude response near the cutoff frequency and the large amount of ripple in the stopband.

Due to the modulation property of the DTFT, the frequency response of the truncated filter is the result of convolving the magnitude response of the ideal filter (a rect) with the DTFT of the truncating window. The DTFT of the truncating window, shown in Figure 7, is similar to a sinc function since it is the DTFT of a sampled rectangular window. Notice that this DTFT has very large sidelobes, which lead to large stopband ripple in the final filter.

Figure 7: DTFT of a rectangular window of size 21.

A truncated impulse response is of finite duration, yet the filter is still noncausal. In order to make the FIR filter causal, it must be shifted to the right by MM units. For a filter of size N=2M+1N=2M+1 this shifted and truncated filter is given by

h ( n ) = ω c π sinc ω c π ( n - N - 1 2 ) n = 0 , 1 , ... , N - 1 0 otherwise . h ( n ) = ω c π sinc ω c π ( n - N - 1 2 ) n = 0 , 1 , ... , N - 1 0 otherwise .

(24)

This time shift of (N-1)/2(N-1)/2 units to the right corresponds to multiplying the frequency response by e-jω(N-1)/2e-jω(N-1)/2. It does not affect the magnitude response of the filter, but adds a factor of -jω(N-1)/2-jω(N-1)/2 to the phase response. Such a filter is called linear phase because the phase is a linear function of ωω.

It is interesting to see that the filter formula of Equation 24 is valid for NN both even and odd. While both of these filters are linear phase, they have different characteristics in the time domain. When NN is odd, then the value at n=(N-1)/2n=(N-1)/2 is sampled at the peak of the sincsinc function, but when NN is even, then the two values at n=N/2n=N/2 and n=(N/2)-1n=(N/2)-1 straddle the peak.

To examine the effect of filter size on the frequency characteristics of the filter, write a Matlab function LPFtrunc(N) that computes the truncated and shifted impulse response of size NN for a low pass filter with a cutoff frequency of ωc=2.0ωc=2.0. For each of the following filter sizes, plot the magnitude of the filter's DTFT in decibels. Hints: The magnitude of the response in decibels is given by |HdB(ejω)|=20log10|H(ejω)||HdB(ejω)|=20log10|H(ejω)|. Note that the log command in Matlab computes the natural logarithm. Therefore, use the log10 command to compute decibels. To get an accurate representation of the DTFT make sure that you compute at least 512 sample points using the command [X,w]=DTFT(filter_response,512) .

• N = 21 N = 21

  (25)
• N = 101 N = 101

  (26)

Now download the noisy speech signal nspeech2.mat , and load it into the Matlab workspace. Apply the two filters with the above sizes to this signal. Since these are FIR filters, you can simply convolve them with the audio signal. Listen carefully to the unfiltered and filtered signals, and note the result. Can you hear a difference between the two filtered signals? In order to hear the filtered signals better, you may want to multiply each of them by 2 or 3 before using sound.

INLAB REPORT

• Submit the plots of the magnitude response for the two filters (not in decibels). On each of the plots, mark the passband, the transition band and the stopband.
• Submit the plots of the magnitude response in decibels for the two filters.
• Explain how the filter size effects the stopband ripple. Why does it have this effect?
• Comment on the quality of the filtered signals. Does the filter size have a noticeable effect on the audio quality?