George Mason University Signals and Systems I Spring 2018 Laboratory Project #4 Assigned: March 5, 2018 Due Date: Laboratory Section on Week of March 26, 2018 Lab Report Your report for this lab will consist of all the analytical (i.e, pencil/paper) work, MATLAB plots and code, and relevant explanations. Each student must do his or her own work on this lab. However, you may ask other students or any of the teaching staff for advice. 1 Prelab (a) We have seen that an important class of linear time-invariant systems is one in which the input and output are related by a linear constant coefficient differential equation, such as d2 y(t) dy(t) +2 + 10y(t) = x(t) 2 dt dt or, more generally, N M X dk y(t) X dk x(t) ak = bk (1) dtk dtk k=0 k=0 We have also seen that if the input to a linear time-invariant system is an exponential, x(t) = est where s = σ + jω is a complex number, then the output will be an exponential of exactly the same form, but scaled by a complex number whose value depends on s and the impulse response of the system, h(t). More specifically, the output is y(t) = H(s)est where Z∞ H(s) = h(t)e−st dt −∞ is the system function. Thus, exponentials are eigenfunctions of linear time-invariant systems. In the special case when s = jω, x(t) is a complex exponential of frequency ω, x(t) = ejωt and y(t) = H(jω)ejωt In this case, Z∞ H(jω) = −∞ h(t)e−jωt dt and H(jω) is called the frequency response of the system (filter). Clearly, H(jω) may be found from the system function by setting s = jω, H(jω) = H(s) s=jω For an LTI system described by a LCCDE of the form given in Eq.(1), we have seen that H(s) is a ratio of polynomials M X H(s) = sM sM −1 B(s) bM + bM −1 + · · · + b1 s + b0 = = k=0 N N −1 N A(s) aN s + aN −1 s + · · · + a1 s + a0 X bk sk = ak sk P (s) Q(s) k=0 where the coefficients ak and bk are the coefficients of the differential equation. From this, the frequency response is found by setting s = jω. The magnitude of H(jω) at a given frequency ω tells us how much a complex exponential is attenuated or amplified in amplitude, and the phase (angle) of H(jω) indicates how much the complex exponential is shifted in phase (delayed or advanced in time). MATLAB has a number of useful m-files to find and plot the frequency response of a system and to filter signals. One of these that you will be using in this lab is freqs, which returns the frequency response H(jω) of a continuous-time system that is specified by the coefficients ak and bk in the numerator and denominator polynomials of H(jω). For example, if we have a linear time-invariant system defined by the LCCDE d3 y(t) d2 y(t) dy(t) + 2 + 10 + y(t) = x(t) 3 2 dt dt dt and if we enter the following commands in MATLAB , >> a = [1 2 10 1]; b = [1]; >> w=linspace(0,5,100); >> H=freqs(b,a,w); then H will contain the values of the frequency response in rad/s at the frequencies specified by the vector w. In this case, H contains 100 samples of H(jω) within the range 0 ≤ ω ≤ 5 rad/s. The magnitude of the frequency response may then be plotted as follows, >> plot(w,abs(H)) >> xlabel(’Frequency (rad/s)’) >> ylabel(’|H(jw)|’); and the plot will be as shown in Figure 1. There are other ways to use MATLAB and these may be found by typing help freqs. Other useful MATLAB functions that you will be using in this lab that you have used before include tf(b,a), impulse(sys,tfinal,dt), and lsim(sys,u,t,x0). Again, for documentation on any of these use the help function. Figure 1: Plot of the magnitude of the frequency response of a filter with coefficients b=[1] and a=[1 2 10 1] . (b) Knowing the systems’ frequency response H(ω) is particularly helpful when your input is a continuous periodic signal. As you have studied, a periodic signal can be represented as a Fourier Series ∞ X x(t) = Xn e nωo t (2) n=−∞ where Xn = 1 To Z x(t)e−nωo t dx. (3) To Notice that x(t) is now made up of a constant, Xn , times an ”everlasting” exponential, enωo t . It is an easy process then to calculate the output to a system knowing its frequency response: ∞ X y(t) = H( nωo ) Xn e nωo t (4) n=−∞ 2 Periodic Signals In this section we will focus on the output to an LTI system defined by its system differential equation when the input is periodic. This will involve determining the Xn coefficients of a Fourier Series expansion for use as an approximation to the input signal. The next step will be to calculate the frequency response H( ω) for use in determining the output y(t). 2.1 Pulse Wave Signal You are given a system whose linear constant coefficient differential equation is d3 y(t) d2 y(t) dy(t) + 59 + 1951 + 9179 y(t) = 9179 x(t) 3 2 dt dt dt (5) The input to this system is the periodic signal shown in Figure 2. Pulse Wave 2 1 t -1 1 2 3 4 -1 Figure 2: Plot of Periodic Input (a) Calculate via hand analysis the Fourier Series coefficients Xn . Now write a MATLAB m-script file to plot the Xn coefficients using the stem function >> nmax=20; >> n=[-nmax:nmax]; >> L_n=length(n); >> >> >> >> >> >> % % % % omega_0=?; % % X_n = zeros(1,L_n); % for i=1:1:L_n X_n(i) = provide end X_n((L_n+1)/2)=?; % % You determine this number Identifies the coefficients to calculate; Coefficients to calculate. Note: zero is at (L_n+1)/2 The radial frequency of the square wave shown in Figure 2. Initialize the coefficients you calculation; Provide a specific answer for X_o if the calculation is faulty (e.g., 1/0). >> figure(1) >> stem(n,abs(X_n)); %Plot the magnitude of X_n From the stem plot, determine nmax where Xnmax is greater than 1% of the largest Xn . How fast do the Xn0 s fall off compared to n? (Hint: you may want to use the ”data cursor” under the ”tools” menu) Instructor Verification (separate page) (b) Approximate the input x(t) using the Fourier System coefficients X−nmax to Xnmax that you calculated above. >> >> >> >> >> >> t=0:0.005:2*To; % Determine To from Figure 2. x=X_n*exp(j*omega_0*n’*t); figure(2) plot(t,real(x)),xlabel(’t’), ylabel(’Partial Sum’) % Plot real(x) to remove any small imaginary part due to % computational error axis([0 2*To -1.5 2.5]) %You provide the period To text(0.05,-0.25,[’max. har. = ’ num2str(nmax)]) How does the Fourier Series approximation to the square wave compare with the actual square wave? What discrepancies do you notice? Instructor Verification (separate page) (c) Determine H( ω) from Eq. (5). For the Fourier Series, ω will be replaced with n ωo so that H( ω) = H( n ωo ). Now calculate the output of your system y(t) given the input x(t) >> >> >> >> >> >> >> for i=1:1:L_n H_n(i)= ? ; Use your derived solution from Equation (5). end y=(H_n.*X_n)*exp(j*omega_0*n’*t); figure(3) plot(t,real(y)) % Be sure to label your plot! axis([0 6 -1.5 2.5]) How does the output compare to the input? What type of filter is h(t)? Instructor Verification (separate page) Triangular Wave 2 1 -1 1 2 3 4 5 t Figure 3: Plot of Second Periodic Input. 2.2 Triangular Wave Signal In this section consider the triangular wave signal shown in Figure 3. (a) Using hand analysis calculate the Fourier Series coefficients Xn . (b) Plot the Xn coefficients using the stem function. (c) What value of nmax did you find so that Dnmax is greater than 1% of the highest Xn ? How does this compare with Part 2.1(a) above? How fast do the Xn0 s of the triangle wave fall off compared to those of the pulse wave? Instructor Verification (separate page) 3 Power Spectrum A power spectrum shows where the power of a signal is distributed across the frequency spectrum. A single sinusoid would be expected to have a majority of its power at the frequency of the sinusoid while other signals may be more uniformly distributed. In this section we’ll look at several wave forms to see how their power is distributed. 3.1 Sinusoid Consider a vector x consisting of a unit-amplitude sinusoid with frequency fo = 330 Hz. Given the sampling frequency fs = 33, 100 Hz, the spacing between samples would then be dt = 1/fs = 3.02115x10−5 seconds. Set the number of samples of your waveform to L = 216 = 65, 536. (a) Write a MATLAB program to create the sinusoid x and plot the first 500 points. Next, listen to the signal x using the soundsc command. The format for this function is soundsc(x,fs). Failure to add the fs will mean you will be listening for a long time! >> >> >> >> >> >> >> >> >> >> >> >> >> fs = 33100; % Sampling frequency fo = 330; % Sinusoid frequency dt = 1/fs; % Sampling interval L = 2^(16); % Number of points in sinusoid t = [0:dt:(L-1)*dt]; x = sin(fo*2*pi*t); figure(5) plot(t(1:500),x(1:500)); % Plot the first 500 points xlabel(’Time (sec)’); title(’Sine’); ylim([-1.25 1.25]) % Expand y-axis grid soundsc(x,fs) (b) The MATLAB code below computes the spectrum of the sinusoid x(t) and graphs the spectrum from 0 to 3000 Hz. >> >> >> >> >> >> >> >> >> Fx = 20*log10(abs(fft(x,L))); % Spectrum of signal x ff = fs*linspace(0,1,L); % Frequency values figure(6) plot(ff,Fx) xlim([0 3000]) % Limits frequencies to 0 -> 3,000 Hz. xlabel(’Frequency (Hz)’); ylabel(’Spectrum Magnitude (dB)’) title(’Sine’) grid What can you say about the distribution? Why isn’t all the power located at fo ? 3.2 Square Wave (a) Next consider a square wave. MATLAB provides a function titled square(x) which creates a unit square wave with frequency 2π Hz. Using the same parameters as given in Section 3.1 above, use the code below to plot the first 500 points of the square wave. Also, listen to the square wave signal using the MATLAB command soundsc. >> >> >> >> >> >> >> >> >> >> fs = 33100; % Sampling frequency fo = 330; % Sinusoid frequency dt = 1/fs; % Sampling interval L = 2^(16); % Number of points in signal t = [0:dt:(L-1)*dt]; x = square(fo*2*pi*t); figure(7) plot(t(1:500),x(1:500)); % Plot the first 500 points xlabel(’Time (sec)’); title(’Sine’); >> ylim([-1.25 1.25]) >> grid >> soundsc(x,fs) % Expand y-axis How does the square wave signal sound when compared to the sinusoidal signal in Section 3.1(a)? (b) Use the MATLAB code below compute and graph the spectrum of the square wave. >> >> >> >> >> >> >> >> Fx = 20*log10(abs(fft(x,L))); % Spectrum of signal x ff = fs*linspace(0,1,L); % Frequency values figure(8), plot(ff,Fx) xlim([0 3000]) % Limits frequencies to 0 -> 3,000 Hz. xlabel(’Frequency (Hz)’); ylabel(’Spectrum Magnitude (dB)’) title(’Square Wave’) grid How does this spectrum compare to that found for the sinusoid in Section 3.1(b)? 3.3 Speech In this section we’ll look at the spectrum of a noisy signal. Noise may be defined as an unknown and undesirable signal, n(t), that is added to a desired signal x(t), i.e., y(t) = x(t) + n(t) An important application of signal processing is to reduce the noise in y(t) by filtering out as much of the noise as possible. Of course, in able to do so, something must be known about the properties of the noise and/or signal. In this lab, you will investigate ho one might remove additive noise from a speech. (a) Load lab04_signal.mat into MATLAB as shown below. The following two variables will be loaded, signal and fs . load lab04_signal.mat To listen to the signal, use soundsc(signal,fs) being sure to include fs . Describe what you hear. (b) Plot the magnitude spectrum of the signal using the same format as before. >> >> >> >> >> >> L=length(signal); Fx_filtered=20*log10(abs(fft(signal,L))); ff=fs*linspace(0,1,L); figure(9) plot(ff,Fx_filtered); % Be sure to label your plot! xlim([0 fs/2]); Compare this spectrum to both the sinusoid and square wave. What can you say about the frequency distribution of this signal? (c) There are a number of standard filters that have been created over the years. The Butterworth filter is one that has also been made into a function in MATLAB called, appropriately, butter. It consists of three parameters butter(n,wc,’ftype’) where n is the order of the differential equation, wc is the cutoff frequency, and 0 f type0 is the type of filter: ’low’,’high’, or ’stop’. We will be using the ’low’ pass filter. Use the code below to plot H(jω) for the filter being used. Next apply the filter as shown in the code and listen to the results. >> >> >> >> >> >> wc = 0.5; % Cutoff frequency is the range 0 < wc < 1. [b,a] = butter(5,wc,’low’) ; figure(10) freqs(b,a,100); % Plots to characteristics of H(jw) signal_filtered = filter(b,a,signal); soundsc(signal_filtered,fs); How does this compare with the original signal? Try several cutoff frequencies (wc) and note any differences. (d) Finally, calculate and plot the spectrum of the filtered signal using the code below. >> >> >> >> >> >> L=length(signal_filtered); Fx_filtered=20*log10(abs(fft(signal_filtered,L))); ff=fs*linspace(0,1,L); figure(11) plot(ff,Fx_filtered); xlim([0 fs/2]); How does this differ from the original signal’s spectrum? Lab #4 ECE 220: Spring 2018 Instructor Verification Sheet For each verification, be prepared to explain your answer and respond to other related questions that the lab TA’s or professors might ask. Turn this page in along with your lab report. Name: Date of Lab: Part 2.1(a): What value of nmax did you find so that Xnmax is greater than 1% of the highest Xn ? Part 2.1(b): How does the Fourier Series approximation to the pulse wave compare with the actual square wave? What discrepancies do you notice? Part 2.1(c): What was the impact of passing the pulse wave through the system defined by the differential equation in (5)? Part 2.2: What value of nmax did you find so that Xnmax is greater than 1% of the highest Xn ? How does this compare with Part 2.1(a) above? How fast do the Xn0 s of the triangle wave fall off compared to those of the square wave? Lab Report Format Signals and Systems 1 ECE 220 Spring 2018 Each lab report that you submit for this class is a formal report that is well-prepared, carefully written, and follows the format given below. The report must include descriptions and analysis of all assigned portions of the lab. Lab reports are to be typed and submitted through the class Blackboard website. Once you have completed your report, please submit it via Blackboard prior to the start of the next lab. Thus, Lab #1’s report is due prior to the start of Lab #2. General Guidelines • Your report should be typewritten, neat and well-organized. • The report must follow the format given below. All analytical work and calculations made should be clearly explained. • MATLAB code is to be placed in report appendices. • It is expected that the report will be grammatically correct with no spelling errors. Points will be deducted for poor grammar and/or spelling errors. You are encouraged to use a spell checker. • Explanations that describe your work must be included in the report. Plots that are submitted within the report must include properly labeled axes and a title. Each plot should be given a figure number and a caption. Points will be deducted if you include unnecessary graphs and plots. • Within the text of your report, when referring to a particular plot, refer to it by number. Detailed Lab Report Format Your report must contain the following sections. Subsections may be added for purposes of clarification. 1. Title Page: This page contains an identification of the Laboratory by number and title. It also must include the name and G-Number of the author and the date of submission. 2. Introduction: This section contains a description of the purpose and objectives of the laboratory project. Do not simply copy text from the assignment. The Introduction should summarize the topics covered in the laboratory and a brief summary of the key results obtained. 3. Main Body: This section is the most detailed part in your report and it will generally contain subsections. Note that it should not contain any MATLAB code. What it does contain is a description of all results obtained during the lab as well as any figures that were generated with MATLAB or other sources. If theoretical calculations are required, these should be detailed and, where appropriate, compared with the experimental portion of the lab. 4. Conclusions: In this Section, summarize the conclusions that you made once the results were obtained. The conclusions should be tied back to the objectives of the Lab Project. This should be a concise section that focuses on the important results obtained and lists the conclusions that resulted from these results. 5. References/Appendices: Include a description of any external references that you used (e.g. web pages, technical papers, etc.) during the lab. The appendices should also contain listings of your MATLAB code. You should document your code with comment lines, where appropriate. This holds for both the primary code as well as any scripts and subroutines that you have written. Remember, if a grader has trouble understanding your MATLAB code, you will not receive full credit for that portion of your report. Adding carefully chosen comments within the code will make it much more likely that your code is fully understood. 1
Purchase answer to see full attachment