ECGR4124/5124 Digital Signal Processing Projects will NOT be accepted late. It is due at the beginning of class. Project 2: The Discrete Fourier Transform & Filtering Part I The Discrete Fourier Transform Point Allocation: (100pts 4124/ 150pts 5124) ECGR 4124 / ECGR5124 Task 1 In this task you must implement the matrix-vector form of the discrete fourier transform as described in class. To do this, you must write two MATLAB functions x_jw=myfft(x_n) and x_n=myifft(x_jw) which must reside in their own MATLAB files: myfft.m and myifft.m respectively. The function x_jw=myfft(x_n) implements the discrete fourier transform by computing the matrix WN and multiplying this matrix times the signal, x_n, which is assumed to be a column vector. The result is a column vector which is the discrete Fourier transform of the input, x_jw. The function x_n=myifft(x_jw) implements the inverse discrete Fourier transform by computing the matrix W−1 N and multiplying this matrix times the signal, x_jw, which is assumed to be a column vector. The result is a column vector which is the inverse discrete Fourier transform of the input, x_n. Test your DFT functions using a MATLAB script project2.m: 1. Generate a sinusoidal sequence x[n] = cos( π 3 n) 0 ≤ n ≤ N − 1 where N = 6. 2. Compute the DFT of x[n], X[k]. 3. Compute the IDFT of X[k], xˆ[n]. 4. Compute the sum of the point-by-point errors in the reconstruction: e = PN n=0 |x[n] − xˆ[n]|. If e > 10−15 you have a bug in your code, find it and fix it! To Turn in: You are responsible for turning in the following objects IN THE ORDER LISTED BELOW. Deductions will be made for reports which are messy and difficult to comprehend. Using your MATLAB program for steps 1-4 provide the following: 1. (10pts) Turn in your MATLAB code consisting of printouts of myfft.m, myifft.m, and project2.m. 2. (10pts) Using subplot() to generate a 3x1 matrix of plots, show the input signal x[n] from (1), the magnitude response of DFT |X[k]| from (2) and the phase response of the DFT from (2) θ[k]. Note: You may consider rounding the DFT values to their 6th digit of precision to prevent numerical instability, i.e., 0.12345678... rounds to 0.123457. 3. (15pts) Explain the following: For each index k, explain the (possibly complex) observed value of the DFT, X[k], in terms of the magnitude and phase of the complex number (there are 5 such numbers). 4. (5pts) Change the value of N such that N = 7 and use subplot() to generate a 2x1 matrix of plots showing the input signal x[n] from (1) the magnitude response of the DFT |X[k]| from (2) and the phase response of the DFT from (2) θ[k]. This DFT should look different! 5. (15pts) Explain why the DFT from steps (2) and (4) are different. 6. (10pts) Change the phase of x[n] by generating a new signal x[n] = cos( π 3 (n − 1)) 0 ≤ n ≤ N − 1 where N = 6 and recompute the DFT. 7. (10pts) Using subplot() to generate a 3x1 matrix of plots, show the input signal x[n] from (1), the magnitude response of the DFT |X[k]| from (2) and the phase response of the DFT from (2) θ[k]. 1 8. (15pts) Explain why the DFT from steps (2) and (7) are different. 9. (10pts) For a sinusoidal sequence x[n] = cos(n) 0 ≤ n ≤ N . Determine if there is a value for N which makes the magnitude response |X[k]| zero everywhere except for two indices of k as in Task 1, step (2). Write a paragraph explaining your solution and show a plot of your results. ECGR 5124 Task 2: Exploring the DFT This task asks students to explore modifications to the DFT in the frequency domain, i.e., frequency domain filtering, and to understand the implications of zero-padding for the DFT. Use your programs from task 1 to complete the list of questions below: 1. Generate a sinusoidal sequence x[n] = cos( π 3 n) 0 ≤ n ≤ N − 1 where N = 6. 2. Compute the DFT of x[n], X[k]. 3. Modify the values of X[k] to generate a new DFT, Xˆ[k], such that the IDFT of Xˆ[k] is complex exponential e −j π 3 n. Write a paragraph explaining your modification and why it generates a complex exponential in the sample (time) domain. 4. Using subplot() to generate a 5x1 matrix of plots, show the input signal x[n] from (1), the magnitude response,|X[k]|, from (2), the modified magnitude response, Xˆ[k] , from (3) and the real and complex parts of the time domain signal xˆ[n]. 5. Modify the values of X[k] to generate a new DFT, Xˆ[k], such that the IDFT of Xˆ[k] is a phase shifted version of the input with a phase shift of φ = π 4 . Write a paragraph explaining your modification and why it generates a phase shift of φ = π 4 in the sample (time) domain. 6. Using subplot() to generate a 4x1 matrix of plots, show the input signal x[n] from (1), the phase response, θ[k], from (2), the modified ˆθ[k] from (3) and the resulting time domain signal xˆ[n]. 7. Modify the values of X[k] to generate a new DFT, Xˆ[k], such that the IDFT of Xˆ[k] is frequency shifted version of the input, xˆ[n] = cos( 2π 3 n) 0 ≤ n ≤ N. Write a paragraph explaining your modification and why it generates a frequency shift to ω = 2π 3 in the sample (time) domain. 8. Using subplot() to generate a 4x1 matrix of plots, show the input signal x[n] from (1), the magnitude response,|X[k]|, from (2), the modified magnitude response, Xˆ[k] , from (3) and the time domain signal xˆ[n]. 9. Again, generate a sinusoidal sequence x[n] = cos( π 3 n) 0 ≤ n ≤ N − 1 where N = 6 as in (2) but now add Nz zeros onto the end of your signal. 10. Compute the DFT of x[n] for the following values of Nz = 0, 1, 10, 100, 1000. Using subplot() to generate a 5x1 matrix of plots, show |X[k]| for each of the values of Nz . Write a paragraph explaining what you see in these magnitude response plots. To turn in: 1. (15pts) Your MATLAB program to implement items 1-10 from task 2 above. 2. (35pts; 5pts per solution) Your solutions to questions 3,4,5,6,7,8 and 10. Note that 3,5,7, and 10 require written responses and 4,6, and 8 require plots which you should label and refer to from your written responses. 2 Part II Frequency Domain Filter Design with MATLAB Point Allocation: (100pts 4124/ 100pts 5124) The objective of this project is to allow students to gain some experience in designing filters using MATLAB’s filter design toolkit. 1 Prelab 1.1 Review of Filters – Frequency Domain The effect of a system on a signal input to it can be described in both time and frequency domains. Here we will study the frequency domain representation of the system as a filter. Different frequency components of the input signal experience different attenuation/amplification and delays as they pass through the filter thus a filter alters the signal input to it in frequency domain (Fig. 1). One of the main applications of a filter is to filter out the unwanted parts of an input signal. The unwanted part might be noise present out of the frequency band occupied by the desired signal (high frequency noise) or other interfering signals that are unwanted, e.g. when you tune to a radio station, the radio receiver will filter out other stations so you would hear a clear sound (Fig. 2). Figure 1: The alteration of frequency components of a signal input to a filter Simple filters can be categorized into four major categories: lowpass, highpass, bandpass and bandstop (or band-reject). Figure 3 shows the ideal forms of these four filters in the frequency domain. As the name suggests a lowpass filter for instance, passes the low frequency components of a signal while blocking the high frequency ones. As you have learned the convolution operation in time domain is equivalent to multiplication in the frequency domain. Thus, multiplying the frequency domain representation of the input signal with that of the filter, high frequency components of the input signal will be multiplied by zero and would thus vanish, while the lower frequency components will be multiplied by one and would pass without any change. 3 (a) (b) Figure 2: Some basic applications of filters: (a) noise reduction by removing out of band (high frequency noise), (b) removing interfering signals: passing of a tuned channel in a radio. In this project, you will learn how to use some tools in MATLAB signal processing toolbox to design filters with your desired characteristics. Figure 3: The ideal frequency response of the four basic filter types: 1) lowpass, 2) highpass, 3) bandpass and 4) bandstop (or band-reject) 1.2 Filter Design Before designing a filter, you should know what parameters need to be known for the design. The ideal filters as shown in figure 3, cannot be physically realized, since they are non-causal, i.e. to generate the output at instance n in time, they require the knowledge of future samples of the input signal, i.e. at n + 1, n + 2,... , which is impossible. The design techniques are as a result, approximations of the ideal case. The digital filter design methods fall into two main categories, FIR filter design and IIR filter design. You will learn more about these types in your future courses. In this project you will be introduced to some tools for basic IIR filter design. First of all you should know the parameters of the filter that you are going to design. Some of these parameters are described in this section. 4 1.3 dB Notation Something you will very often encounter when dealing with filters is the deci-Bell (dB) notation. The reason is that when you’re dealing with a wide range of variation in the magnitude of a quantity, say 5 orders of magnitude difference between maximum and minimum, it is hard to work with the numbers. If you use the logarithm of the magnitude, the variations will be too small! So Alexander Graham Bell, invented this new notation named after him by using 20 times the logarithm of the magnitude. Thus xdB = 20 log10 x Some figures that you should know without using your calculator are 3 dB and its multiples. 3 dB is almost equal to 2 (actually √ 2 = 1.414), thus if you have a quantity that is 3 dB greater than another, it is roughly two times greater than it. (Note that you are using logarithm, thus addition of the logarithms is equivalent to multiplication of the numbers themselves.) By analogy, -3 dB, 6 dB and -6 dB are roughly equivalent to 0.5 (actually √ 1 2 = .707), 4 and 0.25 respectively. 1.4 Passband Parameters The passband of the ideal filter is the band of frequencies where the filter frequency response is not zero. If a signal’s frequency content lies in this band, the signal will pass without any alteration. In case of practical filters, the frequency response does not have an abrupt change from zero to one or the value in passband (Fig 4). Instead, there is a transition from a higher level to a lower level. Hence, for a practical filter the band of frequencies over which the magnitude of frequency response is greater than a threshold is called the passband. This threshold is usually set at -3 dB of the maximum magnitude in passband. The frequency at which the filter’s frequency response magnitude reaches this threshold is called the cut-off frequency of the filter (fc). If the frequency response does not decrease monotonically in the passband the filter is said to have ripples in the passband. The difference between the maximum and minimum magnitude values in the passband defines the passband ripple factor Rp usually given in dB. The maximum value of the magnitude in passband is called the gain of the filter. 1.5 Stopband Parameters The stopband of the ideal filter is the band of frequencies where the filter frequency response is zero. Thus the frequency components of the signal lying in this band will not pass and will be fully removed. Again, since the transition from high to low is not abrupt in practical filters, a threshold level defines this region. So, the band of frequencies, where the magnitude of the frequency response of the filter drops below a threshold level is called stopband. The threshold value is dependent on the application of the filter usually about 40 dB and is usually denoted by As or Rs. If there are ripples in the stopband, the maximum value of the magnitude of frequency response in the stopband shall lie below the threshold. The frequency at which the magnitude of the filter response first falls below the threshold value is called the stopband frequency denoted by fs. 5 Figure 4: The filter parameters: 1) passband, 2) stopbands, 3) cut-off frequencies, 4) stopband frequencies, 5) passband ripple, 6) stopband level, 7) filter gain. 2 Four Classical IIR Filter Types In this section I will introduce and compare the four classical IIR filter types, namely Butterworth, Chebychev Type I and II, and Elliptic filters. Since straightforward algorithms for design and implementation of these filters exist, many software programs help you in designing these filters by simply specifying their parameters. Hence, for the sake of this project you need not know all the mathematical detail involved in the design of these filters. All filters here refer to the lowpass case. You can simply generalize the results to the other types by analogy. 2.1 Butterworth Filter The Butterworth filter of order N, also called the maximally flat filter, is an approximation of the ideal filter, which the first 2N − 1 derivatives of its magnitude squared are zero. As a result the frequency response of this filter decreases monotonically with frequency and |H(f = fc)| = √ 1 2 . The decrease is very slow in the passband and quick in the stopband. In a design problem where no ripple is acceptable in passband and stopband, Butterworth filter is a good choice. Figure 5: The frequency response of a Butterworth filter 6 2.2 Chebychev Type I Filter A Chebychev type I approximates the ideal filter by minimizing the absolute difference between the resulting filter and the ideal filter, over the passband. This filter results in appearance of ripples in the passband. The ripples have a constant magnitude of Rp which is one of the design parameters of this filter. The stopband response is maximally flat. The transition from passband to stopband occurs faster than the Butterworth filter: |H(f = fc)| = Rp. If the ripples in the passband are acceptable, a Chebychev filter usually require a lower-order transfer function than a Butterworth filter for the same specifications. Figure 6: The frequency response of a Chebychev type I filter 2.3 Chebychev Type II Filter A Chebychev type I approximates the ideal filter by minimizing the absolute difference between the resulting filter and the ideal filter, over the stopband. This filter results in appearance of ripples in the stopband. The ripples have a constant magnitude of Rs which is one of the design parameters of this filter. The passband response is maximally flat and |H(f = fc)| = Rs. The transition from passband to stopband occurs slower than the type I filter and for even filter order, it never reaches zero. Its advantage over type I is the absence of ripples in the passband. Figure 7: The frequency response of a Chebychev type II filter 2.4 Elliptic Filter The elliptic filter results in the steepest transition from passband to stopband by allowing ripples in both passband and stopband. This filter usually satisfies the desired filter specifications with the lowest order among the four types studied. The passband and stopband ripples Rp and Rs are also required as design parameters. 7 Figure 8: The frequency response of an elliptic filter 3 Four Classical FIR Filter Types In this section I will introduce and compare three classical FIR filter types: Rectangular Window, Raised Cosine Window, and Least Squares. Since straightforward algorithms for design and implementation of these filters exist, many software programs help you in designing these filters by simply specifying their parameters. Hence, for the sake of this project, you need not know all the mathematical detail involved in the design of these filters. All filters here refer to the lowpass case. You can simply generalize the results to the other types by analogy. 3.1 Windowed FIR Filter Design The generic concept of windowed filter design is to take an example of an ideal continuous filter having the desired filtering properties. This continuous filter is then sampled at the designed sampling rate to generate a (possibly infinite) sequence of filter coefficients. These coefficients are then truncated by selecting a window of samples from the complete sequence such that small coefficients of the filter are discarded. Unfortunately, sampling and subsequent truncation of the filter in this way causes the magnitude for some frequencies to be spread across multiple indices of the DFT in the Fourier domain; a phenomenon known as spectral leakage. Windowing of a simple waveform, like cos(ωt) causes its Fourier transform to develop non-zero values (commonly called spectral leakage) at frequencies other than ω. The leakage tends to be worst (highest) near ω and least at frequencies farthest from ω. If the signal under analysis is composed of two sinusoids of different frequencies, leakage can interfere with the ability to distinguish them spectrally. If their frequencies are dissimilar and one component is weaker, then leakage from the larger component can obscure the weaker’s presence. But if the frequencies are similar, leakage can render them unresolvable even when the sinusoids are of equal strength. –Wikipedia entry for Window Function Windowing functions seek to minimize spectral leakage for discrete FIR filters by attempting to attenuate the discontinuities introduced by truncating the filter. This is accomplished by re-weighting the filter coefficients within the window such that the coefficients at the end-points of the window lie very close to zero. A large variety of windows are available for FIR filter design and each windowing function has specific circumstances where they may be appropriate. Each windowing function strikes a compromise between resolution in the time domain, i.e., accurately preserving sharp variations in the input signal, and resolution in the frequency domain, i.e., preventing spectral leakage. We will examine the simplest window, the rectangular window, which provides no re-weighting at all and the raised cosine or Hann window. At the other extreme of dynamic range are the windows with the poorest resolution. These high-dynamicrange low-resolution windows are also poorest in terms of sensitivity; this is, if the input waveform contains random noise close to the signal frequency, the response to noise, compared to the sinusoid, will be higher than with a higher-resolution window. In other words, the ability to find weak sinusoids amidst the noise is 8 diminished by a high-dynamic-range window. High-dynamic-range windows are probably most often justified in wideband applications, where the spectrum being analyzed is expected to contain many different signals of various amplitudes. In between the extremes are moderate windows, such as Hamming and Hann. They are commonly used in narrowband applications, such as the spectrum of a telephone channel. In summary, spectral analysis involves a tradeoff between resolving comparable strength signals with similar frequencies and resolving disparate strength signals with dissimilar frequencies. That tradeoff occurs when the window function is chosen. 3.2 Rectangular Window Filter The rectangular window performs no re-weighting on the truncated filter coefficients at all. This preserves the performance of the filter in the time or sample domain. Yet, in the frequency domain there is significant leakage leading too poor frequency resolution. Figure 9: The frequency response of a rectangular windowed filter 3.3 Raised Cosine (Hann) Window Filter The raised cosine windowing function provides a compromise between frequency and time (sample) domain resolution. It re-weights the filter coefficients using a cosine as specified in equation (1) w(n) = 0.5 1 + cos 2πn N − 1 (1) where n denotes sample index and N denotes the length of the window. This re-weighting of the filter coefficients serves to improve the frequency resolution of the filter, however, it does cause sharp changes in the input signal to be smoothed due to the low-pass effect of the windowing function (even when the filter is intended to be a band-pass filter). 9 Figure 10: The frequency response of a raised cosine windowed filter 3.4 Least Squares Filter A least squares FIR filter design allows designers to specify a magnitude response having an arbitrary shape as a collection of points connected by straight lines. An optimization algorithm then computes the filter coefficients that have linear phase and best satisfy the specified design shape for the magnitude response by minimizing the sum of squared differences between the computed filter response and the specified filter response. Figure 11: The frequency response of a least squares filter 3.5 MATLAB filterDesigner: Filter Design and Signal Processing Tool In this section you will learn how to use a powerful signal processing tool in MATLAB, sptool which is part of the MATLAB’s signal processing toolbox and helps you with filter design and signal processing with an easy to use GUI. In what follows, you will design a bandpass filter that resembles the bandpass filter in a telephone system. You will then apply this filter to the sound file you used before so you will be able to hear the impact of the filtering operation on the signal by listening to the sound file before and after being filtered. 1. At MATLAB command prompt type filterDesigner 2. In the Frequency Specifications edit box enter the sampling frequency associated with your data file. The sampling frequency of a “.wav” sound file is provided as the second argument returned by the audioread() function, i.e., [ x, fs ]=audioread(’data.wav’) will place the sampling frequency (in Hz.) of the audio data in the ’data.wav’ file in the variable fs and the signal values in the vector x. 3. Select the Design Method: Butterworth IIR filter. 10 4. Select the Response Type: Bandpass. 5. Enter 200 Hz and 1200 Hz for the low and high stopband frequencies respectively. 6. Enter 250 Hz and 1000 Hz for the low and high passband frequencies respectively. 7. From the Magnitude Specifications box, enter 3 dB for the passband attenuation and 20dB for the stopband attenuation on either side of the passband. 8. Click the Design Filter button at the bottom. 9. Your filter is now designed and ready for use! 10. Export your filter to the MATLAB workspace. The representation of your filter will depend upon the design method (IIR, FIR) and type of filter you have selected. For IIR filters, exported filters may be a sequence of second order filters in cascade where each filter output is multiplied by a scalar gain factor, represented by the matrix SOS and the gain vector G. MATLAB will create these two variables in your MATLAB workspace when you export the filter. For FIR filters, exported filters may be a sequence of coefficients which directly represent the values of h[n]. 11. For IIR filters in (SOS,G) form, you can convert the (SOS,G) form of your filter into the more familiar transfer function using the [b, a] = sos2tf(SOS,G) MATLAB function. The resulting transfer function is a polynomial in z (or e jw for the Fourier Transform) where the equation for the numerator polynomial n(z) is given by N(z) = Pn i=1 biz 1−i and for the denominator polynomial, D(z) = Pn i=1 aiz 1−i . The resulting transfer function is then H(z) = b1 + b2z −1 + . . . + bn+1z −n a1 + a2z−1 + . . . + an+1z−n = N(z) D(z) To Turn In: 1. (15pts) For the Butterworth filter design of Steps 1-11, look at the coefficients of SOS and the coefficients of [a,b]. Indicate if one representation of the filter may be more desirable than the other. If so, why? 2. (5pts) Show a Bode plot of the magnitude spectrum for the unfiltered input signal. (50pts) Do the following steps for the Butterworth, Chebychev type II, elliptic, rectangular, raised cosine (Hann) and least squares filter types. Your report should have a 1-page analysis for each filter type generating a total of 6 pages. Each page should include the following information: 1. Plots of the values of the coefficients of the numerator and denominator polynomials ([b, a]) for your designed filter. For those filters having only a single coefficient for the denominator polynomial (FIR filters), there is no need to plot the vector of denominator coefficients a as this is a vector containing only the number 1. For this reason, FIR designs will provide you only the coefficients of the numerator polynomial, b. 2. Apply a representation for your filter to the provided input signal. To apply the transfer function filter, use the MATLAB command y = filter(b,a,x). To apply the (SOS,G) form of your filter, use the MATLAB command y = filtfilt(SOS,G,x). 3. Show a Bode plot of the designed filter magnitude response. Bode plots can be generated most easily using the MATLAB semilog() or loglog() function; note that log10 a 20 = 20 log10 a. Students should use this function to generate all of their Bode plots. Note that we do not visualize the phase response of these filters since all of the filter designs will have approximately linear phase in the passband if you examine the Bode plots of their phase response making their phase response behavior indistinguishable. 4. Show a Bode plot of the magnitude spectrum of the output/filtered signal. Describe how this plot can be predicted from the magnitude spectrum of the unfiltered signal and the magnitude response for your designed filter. 11 5. Listen to the output. Provide a comparative analysis of the filters you have examined. 1. (5pts) Indicate if there are any audible differences between the different filters. What did you notice and can you provide any explanation for the differences you may have noticed based on the analysis information for each filter. 2. (15pts) Compare the filters in terms of their performance using a table. There should be 6 columns of the table which are: (1) stopband ripple, (2) passband ripple, (3) transition band width, (4) passband gain and (5) stopband gain (6) total number of coefficients. Note the values of these parameters for each design and indicate which designs outperform competing filters for each of the respective (6) characteristics. There should be 6 rows for this table, one for each of the designed filters. 3. (10pts) In general, can you provide insight as to contexts where it may be appropriate to apply FIR filters in preference to IIR filters? Notes and supplemental material: Project 2 provides directions on specifying a bandpass filter that passes a specific frequency band. Yet, not all filter design situations will require the same set of parameters. Depending on the type of the filter (IIR, FIR) and the design method, e.g., Butterworth, Least-Squares, Chebyshev, etc., additional parameters may need to be specified. In such cases you are responsible for understanding what these parameters control and determining appropriate values of these parameters to achieve the goals of your filter design problem. Supplemental information for the Second Order Sections (SOS) or Biquadratic representation of an IIR filter are available online (see Wikipedia: Biquad Filter) and MATLAB’s specific representation for Biquadratic filters as provided in their documentation (Mathworks: Biquad Filter). The filter designs you are provides are returned as various forms of transfer functions, H(z). To plot the magnitude response of these transfer functions, one must consider the collection of values where z = e jωk where ωk = 2π N k and for some integer k = 0, 1, . . . , N − 1 and choice of N. In this case you will have the N-point DFT of your transfer function or the frequency response H[k] of your transfer function which is a sampling of the z-transform on the unit circle which also corresponds to a sampling of the DTFT H(z)|z=e jω = H(e jω)|ω=ωk = H[k]. Hence, you can substitute z = e jωk into your transfer function equation and compute values of the frequency response. You can sample this function as densely as you like by changing the value of N to be larger. This is also what the function freqz() in MATLAB does for you. 12