PROGRAMMING EXERCISES FOR CHAPTERS 2-12 OF: "THE SCIENTIST AND ENGINEER'S GUIDE TO DIGITAL SIGNAL PROCESSING" COMMENTS AND SUGGESTIONS ON THESE EXERCISES ARE WELCOME! SEND THEM TO THE AUTHOR AT: SMITH@DSPGUIDE.COM CHAPTER 2: STATISTICS, PROBABILITY AND NOISE Use the following algorithm to find the value of each sample in a test signal, where the function "rnd" returns a random number between 0 and 1: x = rnd if x < 0.6 then x = x-0.2 1. Calculation of the histogram. a. Sketch the pdf of the process that generated this signal. b. Generate and plot a signal containing 300 points from this algorithm. c. Generate and plot the histogram of 1000 points acquired from this algorithm. Use a bin width of 0.01. d. Repeat (b) using 10,000 samples. e. Repeat (b) using 100,000 samples. f. Estimate the peak-to-peak noise on the histograms in (c) to (e). Give your answer as the actual number, not a percentage. g. Which of these three histograms is the best estimate of the pdf? The worst estimate? h. Which histogram has the lowest peak-to-peak noise? The highest peak-to- peak noise? i. Explain the apparent contradiction between your answers in (g) and (h). 2. Calculation of the mean and standard deviation. a. Generate a signal containing 10,000 points from this algorithm. Calculate the mean and standard deviation, using a method of your choice. b. Generate three signals as in (a), and add them together. Calculate the mean and standard deviation of the resulting signal. c. Repeat (b) adding 10 signals. d. Repeat (b) adding 30 signals. e. Make a graph of your data in (a)-(d), proving that when random signals are added, their means add and their standard deviations add in quadrature. In your graphs, the x-axis should be "number of points," and the y-axis should be "mean" or "standard deviation." 3. The Central Limit Theorem. a. Calculate and plot the histogram of the signal in 2(c). b. On this same graph, show a Gaussian curve with the mean and standard deviation determined in 2(c). c. Comment on the match between the two in terms of accuracy and precision. d. Would increasing the number of samples improve the accuracy or precision of the match? e. Would increasing the number of signals added improve the accuracy of the match? CHAPTER 3: ADC AND DAC 1. Generate a 1000 point digital signal that simulates an analog signal, consisting of a sine wave at 30.76 hertz, being sampled at 1 kHz for 1 second. Plot this digital signal with the x-axis labeled "time", running from 0 to 1 second. 2. Simulate sampling the analog signal at 200 Hz, by discarding every 4 out of 5 samples from the digital signal you created in step 1. a. Plot this resampled signal, with the x-axis labeled as "time", running from 0 to 1 second. b. Has aliasing occurred in this signal? Explain. c. What is the digital frequency of this signal? d. Could the original analog signal be reconstructed from this digital signal? 3. Repeat #2 with the sampling at 50 Hz, 33.33 Hz, and 28.57 Hz. CHAPTER 4: DSP SOFTWARE This exercise looks at the problem of adding numbers that are very different is size. Write a computer program(s) to complete the following. 1. Generate a 256 samples long sine wave with an amplitude of one, and a frequency such that it completes 3 full periods in the 256 samples. Represent each of the samples using single precision floating point. We will call this signal, x[ ] 2. Add a constant, k = 300,000 to each of the samples in the signal. 3. Subtract the same constant, k, from each of the samples. Call this reconstructed signal, y[ ] 4. Find the difference (i.e., the reconstruction error) between x[ ] and y[]. Call this signal, d[ ]. 5. Plot the reconstructed signal, y[ ], and the difference signal, d[ ]. 6. Repeat steps 1-5 for k = 3,000,000. 7. Repeat steps 1-5 for k = 30,000,000. 8. Answer the following questions: a. "When floating point numbers of very different size are added, the quantization noise on the [fill in the blank] number destroys the information contained in the [fill in the blank] number. " b. The results of step 7 show that the information in the reconstructed signal is completely destroyed with k is equal to 30 million, or greater. If this exercise were repeated with the sine wave of amplitude 0.001, how large of value of k would be needed to destroy the information in the reconstructed signal? c. If this exercise were repeated using double precision, how large of value for k would be needed to destroy the information in the reconstructed signal? CHAPTER 5: LINEAR SYSTEMS These computer exercises will help you understand "The Fundamental Concept in DSP", as illustrated in Fig. 5-11. 1. Generate and plot the following four signals, each 500 samples long: a. x1[n] = sin(2*pi*n/100) b. x2[n] = 4*exp(-(n-150)^2/300) - exp(-(n-150)^2/2500) c. x3[n] = 1 for 240 < n < 300 -2 for 299 < n < 380 0 otherwise d. x4[n] = rnd+rnd+rnd+rnd+rnd+rnd-3 (where the function rnd returns a random number uniformly distributed between 0 and 1) 2. Generate and plot the sum of these four signals: x[n] = x1[n] + x2[n] + x3[n] + x4[n] 3. Pass each of the five signals through the linear system defined below, creating the signals: y[n], y1[n], y2[n], y3[n] and y4[n]. Using x[n] and y[n] as an example, the system is defined by the difference equation: y[n] = 0.05 x[n] + 0.95 y[n-1] This can be carried out with the following program: y[0] = 0 for n = 1 to 499 y[n] = 0.05*x[n] + 0.95*y[n-1] next n As will be discussed in Chapter 19, this is called a "single-pole recursive low-pass filter." It operates in the same way as an RC circuit in electronics, smoothing the waveform and removing high-frequency noise. 4. Compare y[n] with the sum of y1[n], y2[n], y3[n], and y4[n]. Are they identical? Test this by subtracting one from the other. Plot the result and explain any difference between the two signals. CHAPTER 6: CONVOLUTION 1. Write a program that calculates y[n] = x[n]*h[n], where y[n] is 598 samples, x[n] is 500 samples, and h[n] is 99 samples. 2. Generate an impulse response, h[n], according to the algorithm below. As discussed in Chapter 16, this is the filter kernel for a "low-pass windowed-sinc" filter. When convolved with an input signal, this filter passes sinusoids that have fewer than 25 cycles in 500 samples, and blocks sinusoids with a higher frequency. Make a plot of this signal. for i = 0 to 98 h[i] = 0.31752 * sin(0.314159 * (i-49.00001)) / (i-49.00001) h[i] = h[i] * (0.54 - 0.46 * cos(0.0641114 * i)) next i 3. Test your program by convolving h[n] with the signal described below. What should the output of your program be in response to this signal? Why? x[n] = 1 for n = 0 x[n] = 0 otherwise 4. Generate a more complicated test signal, x[n], that consists of two sinusoids added together. The first sinusoid will have an amplitude of 1, and make 6 complete cycles in the 500 samples. The second sinusoid will have an amplitude of 0.5 and make 44 complete cycles in the 500 samples. Make of plot of this signal. 5. Test your convolution program by filtering the signal you created in 4, with the filter kernel you created in 2. Plot this signal. Has the filter passed the lower frequency signal, while blocking the higher frequency signal? Comment on its effectiveness. 6. Turn the filter kernel into a high-pass filter by changing the sign of all the samples, and then adding one to sample 49 (you can look ahead to Fig. 14-5 if you want more information on this procedure). Test this high-pass filter in the same way as step 4. CHAPTER 7: PROPERTIES OF CONVOLUTION The exercise looks at two different ways of detect a known waveform in a noisy signal. The waveform to be detected is an exponentially decaying pulse. The first detection method is to threshold the first difference of the signal. The second method is to threshold the correlation of the signal with the known waveform. 1. Generate a 600 sample test signal, containing three target pulses: x[n] = exp(-(n-100)/15), if 99 < n < 160 = exp(-(n-300)/15), if 299 < n < 360 = exp(-(n-500)/15), if 499 < n < 560 = 0, otherwise 2. Generate a 600 sample signal of normally distributed random noise with mean = 0, and SD = 1. 3. Generate four test signals with signal-to-noise ratios (SNRs) of 50, 10, 5, and 2.5. Do this by adding the signal from step 1, with an appropriately scaled version of the noise from step 2. For this problem, we will define the SNR to be equal to the peak amplitude of the target waveform (which is one), divided by the standard deviation of the noise. Plot these four test signals. 4. Calculate and plot the first difference of the four test signals. 5. Calculate and plot the correlation of each of the four signals with the 60 point target signal. (Check your program by making sure that the peaks have the proper symmetry). 6. Based on your results, estimate the minimum SNR that the signal must have in order to reliably detect the target by: a. a visual inspection of the waveform b. thresholding the first difference c. thresholding the correlation signal CHAPTER 8: THE DISCRETE FOURIER TRANSFORM 1. Generate and plot a 512 point test signal: x[n] = (sin(2 pi n 0.08) + 2 sin(2 pi n 0.3)) exp(-(n-200)^2 / 60^2) This signal is composed of two sinusoids, with frequencies of 0.08 and 0.3, multiplied by a Gaussian, with a standard deviation of 60. Each sinusoid produces a peak in the frequency domain, while the Gaussian makes these peaks wider and more uniform (more about this in Chapter 9). 2. Take the DFT of this signal. Plot the real and imaginary parts. These plots contain the same information that is contained in the time domain. Is the information in a form that humans can easily understand? Explain, using these plots as an example. 3. Convert the frequency spectrum into polar form. Plot the magnitude on a linear amplitude scale. Is this information in a form that humans can easily understand? Explain. 4. Plot the magnitude on a log amplitude scale. Do the samples between the two peaks have a value of zero? Explain. 5. Plot the phase signal. Identify those sections of the phase that are meaningful, and those sections that are nothing more than meaningless noise. 6. Unwrap the phase in the meaningful sections, and plot. Does this unwrapping make the phase easier to understand? Explain. 7. Convert the polar frequency spectrum back into rectangular notation, and then take the Inverse DFT. Compare the resulting time domain signal with the original? Are they identical? Plot the difference between the two signals, and explain why it is not entirely zero. CHAPTER 9: APPLICATIONS OF THE DFT 1. Create a 600 sample test signal: x[n] = cos(2 pi n 8 / 600) exp(n/200). Add random noise to this signal, with mean = 0 and SD = 1. Plot this signal. 2. Create a 9 point impulse response, h[n]: 1/25, 2/25, 3/25, 4/25, 5/25, 4/25, 3/25, 2/25, 1/25. Plot this signal. What kind of filter is this? 3. Calculate and plot the convolution of x[n] and h[n]. How has this convolution improved the signal? 4. Pad x[n] with zeros to form a 1024 sample signal, calculate the spectrum, and plot the magnitude. 5. In the spectrum of x[n], identify the portion that is mostly signal, and the portion that is mostly noise. What characteristic of the noise generated in step 1 insures that the noise is white? Why does the noise appear irregular in this spectrum, instead of perfectly flat? How could you make the noise appear flatter? 6. Pad h[n] with zeros to form a 1024 sample signal, calculate the frequency response, and plot the magnitude. Identify the frequencies that are passed through the filter (> 90% amplitude), those that are partially passed, and those that are mostly blocked (< 10% amplitude). 7. Multiply the frequency spectrum by the frequency response (see Eq. 9-1), and take the inverse DFT. Is this signal identical to that obtained by direct convolution? Test this by subtracting the two signals, and plotting the result. Explain the result. CHAPTER 10: PROPERTIES OF THE DFT Time domain shifting, time domain aliasing. 1. Generate and plot a 512 sample signal containing a Gaussian curve: x[n] = exp(-(n-200)^2/900) 2. Calculate the DFT, convert to polar form, and plot the magnitude and phase. 3. Modify the frequency spectrum such that the time domain signal will be shifted by 270 samples to the right. 4. Take the Inverse DFT of the modified spectrum. Plot the resulting time domain signal. Has the signal shifted as expected? Has aliasing occurred? Explain. Modulation and frequency domain aliasing. For each of the following steps, plot the time domain signal, calculate the DFT, and plot the magnitude. 5. Generate and plot the following 512 sample signal, calculate the DFT, and plot the magnitude. x[n] = exp(-(n-200)^2/900) sin(2 pi n 0.027) + 0.08 6. Generate and plot the following 512 sample signal, calculate the DFT, and plot the magnitude. x[n] = sin(2 pi n 0.3125) 7. Multiply the time domain signals created in steps 5 and 6, plot, calculate the DFT, and plot the magnitude. On your plot, identify the upper sideband, the lower sideband, and the carrier wave. Has aliasing occurred? 8. Repeat (6), except using a frequency of 0.4922 for the sinusoid (instead of 0.3125). Has aliasing occurred? Why are their two peaks to the left of the carrier wave, instead of only one? Explain. CHAPTER 11: FOURIER TRANSFORM PAIRS 1. Generate the following waveforms, each 512 samples long, centered on sample 200. Normalize the amplitude of each waveform to give it unity area. a. Rectangular pulse, 45 samples wide (i.e., 45 nonzero values) b. Triangular pulse, 45 samples wide (i.e., 43 nonzero values) c. Gaussian, standard deviation = 13.5 e. Hamming window, M = 66 (see Eq. 16-1 and Fig. 16-2a) 2. Calculate the DFT of each waveform and convert to polar form. 3. To allow a fair comparison of the four magnitudes, each of these signals has a "cutoff frequency" equal to 0.01. That is, the width of each time domain waveform is selected to make the magnitude have a value of 0.707 at sample 5 (i.e., 5/512 = 0.01). However, one of these four signals is incorrect, it does not have a cutoff frequency of 0.01. Identify which of the four is incorrect, and modify it so that it does have a cutoff frequency of 0.01. Explain your modification. 4. Plot each of the four time domain waveforms (the correct versions), each of the four magnitudes, and each of the four magnitudes on a logarithmic amplitude scale. 5. Suppose that these four waveforms were being used as filter kernels. The goal is to pass a signal at frequency 0.005, while blocking noise in the frequency band of 0.03 to 0.4. Answer the following. a. Why must the four time domain waveforms be normalized to have the same area in this comparison? What would be the affect on the magnitudes if they were allowed to have different areas? b. Are any of these four filters better or worse than the others in how well the signal is passed? Explain. c. Are any of these four filters better or worse than the others in how well the noise is blocked? Explain. d. Which of these four filters is the best for this application? e. Which of these four filters is the worst for this application? CHAPTER 12: THE FAST FOURIER TRANSFORM 1. Write (or copy) an FFT and an IFFT subroutine. 2. Test that the subroutines are operating correctly: a. Generate two random signals, each containing 256 samples, and place them in the real and imaginary parts of the time domain. Take the FFT and then the IFFT. Calculate the difference between the original signals, and the reconstructed signals. Is the difference zero? Explain. Make a plot of the difference for the real part. b. Place a random signal into the real part of the time domain, and zeros into the imaginary part. Take the FFT and convert to polar form. Plot the real & imaginary parts. (To make the frequency spectra easier to graph, you will probably want to use random numbers that have a mean of zero). Does the frequency domain have the proper symmetry? Explain 3. Use the FFT in a spectral analysis problem: a. Generate a 1,024,000 sample signal containing normally distributed random noise with zero mean and a standard deviation of one, plus a sine wave with amplitude 0.1, and a frequency of 0.2. b. Plot 1024 samples from this signal. Is the sine wave visible? c. Take the FFT of these 1024 samples and plot the magnitude. Is the sine wave visible? d. Use the method of Chapter 9 to calculate the average frequency spectrum of the 1,024,000 point signal (break the signal into 1024 sample segments, calculate the FFTs, and average the magnitudes). Is the sine wave visible in the averaged spectrum? 4. Answer the following questions. a. If the DFT by correlation were used in step 2a instead of the FFT, approximately how much larger in amplitude would the difference signal be? b. In step 2b, why will the frequency spectra be easier to plot if the time domain signals have a mean of zero? c. In step 2b, if the random signal were placed in the imaginary part, and all zeros placed in the real part, how would the frequency domain be changed?