====== ECEn 380 MATLAB Assignments ====== Welcome to the 380 section for the ECEn MATLAB wiki. For a brief refresher on MATLAB, [[ matlab_guide | look here.]] The syntax for the problem numbers is such:


[Language][Class].[Book authors][Chapter].[Section].[Problem]

Thus for help on problem "MT380.UY6.8.1", for example, you can refer to section 6-8 in the textbook (page 279). [[ EcEn 380 Full Table of Content ]] ^ Chapter ^ Assignments ^ Course Topic ^ MATLAB Topic ^ | Ch 1 | [[ecen_380_assignments#Chapter 1: Signals | Signals ]] ||| | ::: | [[ecen_380_assignments#Signal Transformations | MT380.UY1.2.1 ]] | Signal Transformations | [[matlab_guide#Array Operations| Array Operations]],[[matlab_guide#Graph Function| plot]] ,[[matlab_guide#Annotation| Graph Annotation]] ,[[matlab_guide#Figures| Figures]] | | ::: | [[ecen_380_assignments#Signal Power and Energy | MT380.UY1.5.1 ]] | Signal Power and Energy | sum (used to evaluate numerical integrals) | | Ch 2 | [[ecen_380_assignments#Chapter 2: LTI Systems | Linear Time-Invariant (LTI) Systems ]] ||| | ::: | [[ecen_380_assignments#MT380.UY2.3.0 | MT380.UY2.3.0 ]] | Convolution Aid | | | ::: | [[ecen_380_assignments#MT380.UY2.3.1 | MT380.UY2.3.1 ]] | Convolution with function files | [[matlab_guide#Function Files| Function Files]], [[matlab_guide#For-loop| For-loop]] | | ::: | [[ecen_380_assignments#MT380.UY2.3.1 | MT380.UY2.3.2 ]] | Convolution with arrays | [[matlab_guide#For-loop| For-loop]], [[matlab_guide#Ones and Zeros| Zeros]] | | ::: | [[ecen_380_assignments#MT380.UY2.5.1 | MT380.UY2.5.1 ]] | Convolution Properties | Review | | ::: | [[ecen_380_assignments#MT380.UY2.8.1 | MT380.UY2.8.1 ]] | Imp.Resp. of 2nd Ord LCCDEs | [[matlab_guide#Roots| Roots]] | | ::: | [[ecen_380_assignments#MT380.UY2.9.1 | MT380.UY2.9.1 ]] | Car Suspension | Review | | Ch 3 | [[ecen_380_assignments#Chapter 3: Laplace Transform | Laplace Transform ]] ||| | ::: | [[ecen_380_assignments#MT380.UY3.2.1 | MT380.UY3.2.1 ]] | Poles and Zeros | [[matlab_guide#Mesh| Mesh]] | | Ch 4 | [[ecen_380_assignments#Chapter 4: Applications of the Laplace Transform | Applications of the Laplace Transform ]] ||| | ::: | [[ecen_380_assignments#MT380.UY4.5.1 | MT380.UY4.5.1 ]] | Op-Amp Circuits | Review | | Ch 5 | [[ecen_380_assignments#Chapter 5: Fourier Analysis Techniques | Fourier Analysis Techniques ]] ||| | ::: | [[ecen_380_assignments#MT380.UY5.4.1 | MT380.UY5.4.1 ]] | Fourier Series Coefficients | [[matlab_guide#For-loop| For-loop]] | | ::: | [[ecen_380_assignments#MT380.UY5.12.1 | MT380.UY5.12.1 ]] | Circuit Analysis with Fourier | Review | | Ch 6 | [[ecen_380_assignments#Chapter 6: Applications of the Fourier Transform | Applications of the Fourier Transform ]] ||| | ::: | [[ecen_380_assignments#MT380.UY6.8.1 | MT380.UY6.8.1 ]] | Butterworth Filters | Review | | ::: | [[ecen_380_assignments#MT380.UY6.12.1 | MT380.UY6.12.1 ]] | Sampling Theorem | [[matlab_guide#fopen| fopen]],[[matlab_guide#fscanf| fscanf]],[[matlab_guide#fclose| fclose]], [[matlab_guide#interp| Interpolate]],[[matlab_guide#sound| sound]] | | Ch 7 | [[ecen_380_assignments#Chapter 7: Discrete-Time Signals and Systems | Discrete-Time Signals and Systems ]] ||| | ::: | [[ecen_380_assignments#MT380.UY7.2.1 | MT380.UY7.2.1 ]] | D-Time Signal Functions | rat | | ::: | [[ecen_380_assignments#MT380.UY7.13.1 | MT380.UY7.13.1 ]] | D-Time Fourier Series | Review | | ::: | [[ecen_380_assignments#MT380.UY7.15.1 | MT380.UY7.15.1 ]] | Discrete Fourier Transform | Review | | Extras | [[ecen_380_assignments#Music Harmonics | Music Harmonics ]] | FFT,Fourier Series Coeff | FFT Plot ====== Chapter 1: Signals ====== ===== Signal Transformations ===== ==== MT380.UY1.2.1 ==== == Description == Signal Transformations Objective: Help student gain a visual understanding of several basic signal transformations. Commands: array operations, plot, labels, figure. Background Information: In MATLAB and other programming environments continuous signals cannot be represented by vectors or arrays, since these are inherently discrete. However, it is possible to approximate a continuous signal by taking samples of the signal, where the time interval between samples is small enough that the signal exhibits only small changes in amplitude between adjacent samples. For example, consider the signal x(t) = exp(-0.5*t)*u(t). This signal "turns on" at t = 0, and then exponentially decays towards 0 with an initial slope of -0.5 (at t = 0). Given this initial slope, the time interval between samples should be significantly less than 0.5. A nice choice might be 100 times smaller, although this will clearly depend on the application. Exercise: Consider the signal x(t) = sin(2*pi/4*t). This is a continuous signal, "turning on" at time t = negative infinity and continuing towards t = positive infinity. Clearly we can't plot this signal for all time in Matlab. In this exercise, you should limit your time axis from -6 s to 6 s, and use a sampling period (interval between time points) of 0.01 s. First, set up a vector t holding your time values, and create the signal x(t) = sin(2*pi/4*t). This can be done with the following code: t_start = -6; t_end = 6; t_interval = 0.01; t = t_start:t_interval:t_end; x = sin(2*pi/4*t); Then do the following: (1) Plot the signal x(t) (2) Create the following transformations of the signal x(t), and plot each of them in the same figure but different plots using the subplot command: (a) y1(t) = x(2*t) (b) y2(t) = x(0.5*t) (c) y3(t) = x(t-1) (d) y4(t) = x(-t) (e) y5(t) = x(4-2*t) Give each plot a title to distinguish the original signal from its transformations. (3) Answer the following questions: (a) Explain the behavior you observe in 2(a) and 2(b). (b) Explain in words how the signal x(t) is transformed in 2(e). == Solution Image == {{:380matlab:ch1:ul_1_2_si_signaltransformationssolutionimage.jpg?400|}} == Solution == {{:380matlab:solutions:ch1:ul_1_2_signaltransformations.m|}} ===== Signal Power and Energy ===== ==== MT380.UY1.5.1 ==== ==Description == UY_1_5_1 Signal Power and Energy Objective: Help students gain a better understanding of the relation between frequency, total energy, and average power for sinusoidal signals of finite time duration.Also,to teach students how to approximate an integral in MATLAB. Commands: sum (used to evaluate numerical integrals) Exercise: Consider the signal x(t) = cos(2*pi*f*t)*(u(t) - u(t-6)), where t represents time and f represents frequency in Hz. (1) Create four signals of this form with frequencies of 1, 10, 50, and 1000 Hz over the internal 0 <= t <= 6s (the interval where they are non-zero), as shown below: f = [1 10 50 1000]; t_start = 0; t_end = 6; delta_t = 0.0001; t = t_start:deta_t:t_end-delta_t; y = cos(2*pi*f'*t); (2) Calculate the total energy and average power (over the time interval where the signal is non-zero) of each of the four signals at different frequencies. Do this by using MATLAB to approximate (or numerically evaluate) the appropriate integrals. For example, the total energy in the signal will be given by the integral of the absolute value of the signal squared over the interval 0 <= t <= 6s. You can approximate this integral by summing up the areas of rectangles that are delta_t wide and |y(t)|^2 across the entire interval. DO NOT USE THE INTEGRAL FUNCTION. (3) Report both the total energy and average power (over the interval 0 to 6s for all four signals. (4) Answer the following question: (a) What impact does frequency have on the total energy and average power of a sinusoidal signal analyzed over a finite time duration? ==Solution Image == {{:380matlab:ch1:uy_1_5_1_si1_poweravgsolutionimage1.jpg?400|}} {{:380matlab:ch1:uy_1_5_1_si2_poweravgsolutionimage2.jpg?400|}} == Solution == {{:380matlab:solutions:ch1:uy_1_5_1_poweravg.m|}} ====== Chapter 2: LTI Systems ====== ===== Convolution Aid ===== ==== MT380.UY2.3.0 ==== A document to help gain an intuition behind convolution. {{:380matlab:ch2:convolutionguide.docx|}} The file below work together to create a MATLAB gui that simulates convolution Please do not modify the following two files. \\ {{:380matlab:ch2:convolutiongui.m|}} \\ {{:380matlab:ch2:convolutiongui.fig|}} \\ You can modify the following three function files to change your signals x(t) and h(t) \\ {{:380matlab:ch2:x_t_gui.m|}} This file modifies the signal x(t) \\ {{:380matlab:ch2:h_t_gui.m|}} This file modifies the signal h(t) \\ {{:380matlab:ch2:time_parameters_gui.m|}} This file returns the time parameters used for x(t) and h(t) \\ ===== Convolution Function File ===== ==== MT380.UY2.3.1 ==== MT380_UY2_3_1_Convolution Function File Objective: Implement an approximation to continuous-time convolution in Matlab. Commands: creating function files Background Information: The convolution between a function x(t) and h(t) is given by the "convolution integral": y(t) = integral(x(tau)*h(t-tau)*dtau) from tau = -infinity to +infinity It is clearly impossible to numerically approximate an integral with infinite bounds in Matlab. However, if both x(t) and h(t) are of finite duration, then the integral over tau is only nonzero for a range of tau values equal to the sum of the widths of the two signals. As we have done with other continuous time signals, we can approximate y(t), x(t) and h(t) with discrete approximations in Matlab, and perform the convolution integral numerically by summing trapezoids. Exercise: In this exercise, you will write code to approximate continuous- time convolution for two signals of finite duration: x(t) = u(t) - u(t-4) and h(t) = 2*(u(t) - u(t-2)). To do this, you will create two functions. The first will be a function called x_t, which takes as a parameter the value t and returns x(t) (the signal x(t) evaluated at time point t). function [x] = x_t(t) The second function will be called h_t, and, just like x_t, it will take as a parameter the value t and return h(t) (the signal h(t) evaluated at time point t). function [h] = h_t(t) Convince yourself that since x(t) is only non-zero from 0 < t < 4 and h(t) is only non-zero from 0 < t < 2, y(t) will only be non-zero from 0 < t < 6. Once you have created the function x_t and h_t, you will write code that evaluates y(t) over the interval -2 < t < 8 (just so you can see where the convolution goes to zero as well) by numerically performing the convolution integral. (1) Create function files for x(t) and h(t). The code for the function x_t is provided below. This code should be put in a file called x_t.m. function [x] = x_t(t) if t < 0 || t > 4 x = 0; else x = 1; end end Again note that this function takes a value t and returns x(t) (the signal x(t) evaluated at the time t that is passed in). For example, x_t(-1) will return 0, and x_t(1.5) will return 1. (2) Create a time array that contains the time points from -2 < t < 8 over which you are going to evaluate the convolution integral, and create a vector y of the same length as t to hold your output: t_beg = -2; d_t = 0.01; t_end = 8; t = t_beg:d_t:t_end; y = zeros(1, length(t)); (3) Now cycle through all of the time values in t using a for loop and evaluate the convolution integral at each value of t (using a nested for loop over a new time variable tau). We need to make sure we evaluate the convolution integral over a suitable range of times tau. In this case, we are going to use -10 < tau < 10. This is a bigger interval than we need. In the questions, you will be required to determine what the smallest tau interval that can be used in this case is. Here is some code to get you started: tau_beg = -10; d_tau = 0.01; tau_end = 10; for m = 1:length(t) for tau = tau_beg:d_tau:tau_end y(m) = y() + INSERT YOUR CODE HERE end end (3) Plot y(t) over the interval -2 < t < 8. Question: a) What is the smallest range of values you could have used for tau in this case? b) This code can take a few seconds to run. Could you make the code run more quickly by only evaluating the convolution integral over tau values where x(t) * h(t - tau) is non-zero? Over what bounds of tau does the tau integral need to run in terms of t? Integrate this change in your code. Hint: use the array t that is given to you. for m = 1:length(t) for tau = CHANGE THE RANGE OF INTEGRATION IN TERMS OF t y(m) = ... end end c) What effect does using a larger or smaller value of dtau have when numerically convolving x(t) with h(t)? == Solution Image == {{:380matlab:ch2:mt380_uy2_3_1_si_convolutionffsolutionimage.jpg?400|}} == Solution == {{:380matlab:solutions:ch2:h_t.m|}} \\ {{:380matlab:solutions:ch2:x_t.m|}} \\ {{:380matlab:solutions:ch2:mt380_uy2_3_1_convolutionff.m|}} \\ ===== Convolution With Arrays ===== ==== MT380.UY2.3.2 ==== MT380.UY2.3.2 Convolution of two arrays Objective: Understand the mechanics of convolution in order to create a function that will approximate continuous time convolution. Commands: for-loop, zeros, functions Exercise: Create a function that can approximate continuous convolution. The functions parameters should consist of two arrays representing signals, x(t) and h(t), and a time_interval that represents the time increments in an array that represents time. ex: t = 0:time_interval:5. Your time interval will serve as dtau in y(t) = integral (x(tau)*h(t-tau)dtau, -inf, +inf). The function declaration should look like function [y] = convolution(x, h, time_interval). with x, and h representing two signals. a) To test your function let x(t) = h(t) = u(t)-u(t-4) and convolve the two signals together. b) Plot the computed signal. Some helpful information is below. Convolution is defined as y(t) = integral (x(tau)*h(t-tau)dtau, -inf, +inf) Working with causal signals, convolution can be rewritten as y(t) = integral (x(tau)*h(t-tau)dtau, 0, t) with 0 being the lower bound of integration and t being the upper bound of integration. In order to do Convolution in MATLAB the formula has to change since MATLAB can only approximate a continuous signal. Time is lost and is replaced by indexing. For example consider the code: t= 0:.1:1-.1; t(s) = [0 .1 .2 .3 .4 .5 .6 .7 .8 .9]; x = sin(pi/0.2*t); x = [0 1 0 1 0 1 0 1 0 1] sin(pi/o.2*t(2)) = x(2) To visually demonstrate this consider the two arrays below. let t = [0 1 2 3 ] then t(1) represents 0s let x(t) = [1 2 3 4] then x(1) = 1 at t = 0s let h(t) = [3 2 1 0] then h(4) = 0 at t = 3s let y(t) = [3 8 14 20 11 4 0] then y(1) = 3 at t = 0s Because indexing in MATLAB begins with 1 instead of 0, the convolution formula needs to be adapted for proper indexing. y(index) = integral (x(tau)*h(index-tau)dtau, 1, index) tau cannot start at 0 since x(0) does not exist in MATLAB. This means that tau must start at 1. However, that creates a problem since both index and tau starts at 1. Thus h( (index = 1) - (tau =1 ) ) = h(0) is out of range. To fix this, the equation needs to be adapted even further. y(index) = integral (x(tau)*h(index+1 - tau)dtau, 1, index) Below is an example of this equation. Assume that dtau = 1. at index < 1 y(index <1) = integral (x(tau)*h(index+1<2-tau)dtau, 1, index) | | x[tau] | 1 2 3 4 h[index-tau] 1 2 3 | ------------------------------------------------------------- tau -4 -3 -2 -1 0 1 2 3 4 5 6 During this period, the two signals do not come in contact with each other, thus y[index <1] must be zero. This is the property of a causal signal convolved with a causal system. at index = 1 y(1) = integral (x(tau)*h(2-tau)dtau, 1, index) | | x[tau] | 1 2 3 4 h[index-tau] 1 2 | 3 ------------------------------------------------------------- tau -3 -2 -1 0 1 2 3 4 5 6 At index = 1, and with tau ranging from 1 to index, the two signals only have non-zero values at tau = 1. This indicates that y(1) = 1*3. at index = 2 y(2) = integral (x(tau)*h(3-tau)dtau, 1, index) | | x[tau] | 1 2 3 4 h[t-tau] 1 | 2 3 ------------------------------------------------------------- tau -3 -2 -1 0 1 2 3 4 5 6 At index = 2, and with tau ranging from 1 to index, the two signals both have non-zero values when tau = 1, and 2. This indicates that y(2) = 1*2 + 3*2. This is because y(t) depends on the integration (summation in discrete time) of the two signals being multiplied together. at index = 3 y(3) = integral (x(tau)*h(4-tau)dtau, 1, index) | | x[tau] | 1 2 3 4 h[t-tau] | 1 2 3 ------------------------------------------------------------- tau -4 -3 -2 -1 0 1 2 3 4 5 6 At index = 3, and with tau ranging from 1 to index, the two signals both have non-zero values when tau = 1,2, and 3. This indicates that y(3) = 1*1 + 2*2 + 3*3. at index = 4 y(4) = integral (x(tau)*h(5-tau)dtau, 1, index) | | x[tau] | 1 2 3 4 h[t-tau] | 1 2 3 ------------------------------------------------------------- tau -4 -3 -2 -1 0 1 2 3 4 5 6 At index = 4, and with tau ranging from 1 to index, the two signals both have non-zero values when tau = 2,3 and 4. This indicates that y(4) = 2*1 + 3*2 + 4*3. at index = 5 y(5) = integral (x(tau)*h(6-tau)dtau, 1, index) | | x[tau] | 1 2 3 4 h[t-tau] | 1 2 3 ------------------------------------------------------------- tau -4 -3 -2 -1 0 1 2 3 4 5 6 At index = 5, and with tau ranging from 1 to index, the two signals both have non-zero values when tau = 3 and 4. This indicates that y(5) = 3*1 + 2*4. at index = 6 y(6) = integral (x(tau)*h(7-tau)dtau, 1, index) | | x[tau ] | 1 2 3 4 h[t-tau] | 1 2 3 ------------------------------------------------------------- tau -4 -3 -2 -1 0 1 2 3 4 5 6 At index = 6, and with tau ranging from 1 to index, the two signals both have non-zero values when tau = 4. This indicates that y(6) = 4*1. at index > 6 y(index>6) = integral (x(tau)*h(index+1>7-tau)dtau, 1, index) | | x[tau] | 1 2 3 4 h[t-tau] | 1 2 3 ------------------------------------------------------------------ tau -4 -3 -2 -1 0 1 2 3 4 5 6 7 At index > 6, and with tau ranging from 1 to index, the two signals only have non-zero values. This indicates that y(>6) = 0. Because of the time reversal, we need to ensure that MATLAB doesn't try to index a part of an array that is out of bounds. There are many ways of doing this. The easiest for most people will be an if-statement that checks to see if the index is out of bounds. Another way of doing this is by zero padding. Zero padding is when you concatenate an array with zeros to ensure that MATLAB doesn't index out of bounds. == Solution Image == {{:380matlab:ch2:mt380_uy2_3_1_si_convolutionffsolutionimage.jpg?400|}} == Solution == {{:380matlab:solutions:ch2:convolution.m|}} ===== Convolution Properties ===== ==== MT380.UY2.5.1 ==== == Description == UY2_5_1 Convolution Properties Objective: Gain a visual understanding of the following convolution properties: Commutative, Associative, Distributive, Causal*Causal, Time-shift, Width, and Area. Commands: subplot, figure. Exercise: Use the following functions in this exercise: x1(t) = 2u(t) - 2u(t-2) x2(t) = 2r(t) -4r(t-1) + 2r(t-2) x3(t) = 2u(t) - 3u(t-1) + 3u(t-t) -2u(t-4) x4(t) = r(t)u(t)u(-(t-1)) x5(t) = u(t) - u(t-4) Create two figures. In one figure plot each signal in its own subplot. In the other figure, plot the results to parts a to e. The '*' symbol denotes convolution. a) Compute y1(t) = x2(t)*x3(t) and y2(t) = x3(t)*x2(t) You should have four subplots in this first figure. What is the difference between y1(t) and y2(t)? If there is no difference, state that. b) Compute y3(t) = [x2(t)*x3(t)]*x4(t) and y4(t) = x2(t)*[x3(t)*x4(t)] What is the difference between y3(t) and y4(t)? If there is no difference, state that. c) Compute y5(t) = x3(t)*[x1(t) + x2(t)] and y6(t) = x3(t)*x1(t) + x3(t)*x2(t) What is the difference between y5(t) and y6(t)? If there is no difference, state that. d) Compute y7(t) = x2(t-1)*x3(t-2). Compare y7(t) with the answer in part a. What is the delay in y(t), and how does it relate to the delay in the two signals x2(t-1) and x3(t-2)? e) Compute y8(t) = x1(t)*x5(t). How does the area of y8(t) relate to the areas of x1(t) and x5(t). Use the previous plots created to answer the following questions. What is the relation between the width of the convolved signal (y1(t)-y8(t)) and the width of the signals that were convolved (x1(t)-x5(t)). Signals x1(t)-x5(t) are all causal signals. If you convolved any two or more signals together, the output signal must also be (causal / non-causal) == Solution Image == {{:380matlab:ch2:uy2_5_1_si1convolutionpropertiessolutionimage1.jpg?400|}} {{:380matlab:ch2:uy2_5_1_si2convolutionpropertiessolutionimage2.jpg?400|}} == Solution == {{:380matlab:solutions:ch2:uy2_5_1_convolutionproperties.m|}} ===== Imp.Resp. of 2nd Ord LCCDEs Aid ===== ==== MT380.UY2.8.0 ==== Differential Equations-Euler's Method {{:380matlab:ch2:differential_equations_eulersmethod.docx|}} ===== Imp.Resp. of 2nd Ord LCCDEs ===== ==== MT380.UY2.8.1 ==== == Description == UY2_8_1 Impulse Response of 2nd Order LCCDEs Objectives: Become familiar with Euler's method for solving 1st order differential equations; the effect of overdamped, critically damped, and underdamped systems to a sinusoidal frequency; and approximating derivations in matlab. Commands: review roots. Exercise: Use the figure below for this problem. The value of the capacitor is 2e-3F, the value of the inductor is 1H, and the resistor will take on four different values R = [1,10,20,40]. The input signal is Vs = | cos(2*pi*f*t) for 0 <= t < 1(s) | 0 for 1(s) <= t <= 3(s) and f represents the frequency in Hz, and is 22 Hz. This frequency is close to the resonate frequency of the system. Vout is Vo. You will create an impulse response of the system for every value of R. When you plot, create four different figures for each different resistor value. in each figure, create subplots: 1 for h1(t), 1 for h2(t), 1 for hc(t), 1 for h(t), and 1 for h(t)*Vs (Convolution) When you are creating the impulse responses you will have to create a time array for each of them. t = 0:time_interval:time_end. For the time_interval use 2/1000, and for time_end use 5*tau. The behavior of the impulse response can be captured in a time duration of 5*tau. The value of tau will change depending on the value of R. Remember, tau is the real part of the root of the characteristic equation. Begin by creating your 2nd order differential equation based on the image provided. a) Split the second order differential equation into to coupled first-order equations. Use Euler's method to approximate h1(t) and h2(t). Compare them to the analytical solutions provided in the book by plotting them. Remember to do this for each value of R - analytic approach. Used for comparison - h1A = exp(p(1).*t); - h2A = exp(p(2).*t); b) Convolve your h1(t) and h2(t) to get hc(t). Compare it to the analytical solution provided in the book by plotting them. - analytical solution to hc - hcA = (1/(p(1)-p(2)))*(exp(p(1).*tc)-exp(p(2).*tc)); - tc is the time array for hcA c) Approximate the derivative of hc to get h(t). Compare it to the analytical solution provided in the book by plotting them. - analytical solution - hA = (1/(p(1)-p(2)))*(p(1)*exp(p(1).*tc)-p(2)*exp(p(2).*tc)); d) Convolve your four h(t)s, one for each resistor value, with the input signal and plot them. e) Questions: 1) How close of an approximation is Euler's method, and for what value of R is the approximation the worse and why? 2) Using the graph that shows Vout, look at the time after the input signal stops (t = 1s). What is the relationship between the value of R and Vout after (t = 1s). Why does it behave like this? For further exploration, change the resonate frequency of the input signal and see what happens. For help with this assignment see the document [[ecen_380_assignments#MT380.UY2.8.0 | Differential Equations-Euler's Method]] == Figure == {{:380matlab:ch2:uy2_8_1_impresp2ndordercircuit.png?300|}} == Solution Image == {{:380matlab:ch2:uy2_8_1_si1_impresp2ndordersolutionimage1.jpg?400|}} {{:380matlab:ch2:uy2_8_1_si2_impresp2ndordersolutionimage2.jpg?400|}} {{:380matlab:ch2:uy2_8_1_si3_impresp2ndordersolutionimage3.jpg?400|}} {{:380matlab:ch2:uy2_8_1_si4_impresp2ndordersolutionimage4.jpg?400|}} == Solution == {{:380matlab:solutions:ch2:uy2_8_1_impresp2ndorder.m|}} ===== Car Suspension System ===== ==== MT380.UY2.9.1 ==== Euler's method 2nd order Objective: Use Euler's method for solving 2nd order differential equations to approximate the displacement of a car as it drives over a pot hole. Exercise: You are driving a car at 10m/s when you run over a pot hole that is 0.2 meters deep and 0.5 meters wide. The car weighs 1800 kg (450 kg per wheel), and the suspension system has a spring constant of 10^5 N/m and a Viscous damper of 5*10^3 Ns/m. a) Use Euler's method to plot the path of the car starting from when it first hits the pothole. Hint- you will need the cart to start at level ground. b) compare the results from part a with the analytical solution. Look at table 2-3 on page 71 and the example of driving over a pothole on page 72 to compute the analytical solution. Plot it on the same graph used in part a. c) Compare and contrast the two plots == Solution == {{:380matlab:solutions:ch2:uy2_9_1_eulermethod2ndorder.m|}} ====== Chapter 3: Laplace Transform ====== ===== Laplace Transform Aid ===== ==== MT380.UY3.0 ==== {{:380matlab:ch3:laplacetransformgui.fig|}} {{:380matlab:ch3:laplacetransformgui.m|}} ===== Poles and Zeros ===== ==== MT380.UY3.2.1 ==== == Description == S Domain Graph Objective: Learn how to plot the magnitude and phase of a Laplace transform in order to gain a visual understanding of the s-domain Command: review mesh-grid, mesh and abs. Exercise: Consider the system with impulse response h(t) = exp(-0.9*t)cos(4*pi*t)u(t). The transfer function of the system is (just look in the table!): H(s) = (s+0.9)/((s+0.9)^2 + ((4*pi)^2); We are going to evaluate the transfer function H(s) for different values of s in the complex plane. The variable s can be split into its real part and its imaginary part: s = sigma + j*omega; with sigma representing the real part and omega representing the imaginary part. We'll evaluate H(s) for values of sigma in the range -3 <= sigma <= 3 and values of omega in the range from -6*pi < omega < 6*pi: sigma = -3:0.02:3; omega = 2*pi*(-3:0.02:3); Note that we have made it so the sigma and omega arrays have equal numbers of elements (in this case 301 elements). We now want a matrix of s values that will be a 301x301 matrix with the complex value of s at each point. The following code using the "meshgrid" command accomplishes this: [sigmaGrid,omegaGrid] = meshgrid(sigma,omega); sgrid = sigmaGrid + j*omegaGrid; Now sgrid should be a 301x301 matrix, where every entry corresponds to a different value of s in the complex plain. (a) Evaluate the Laplace transform at every value of s contained in sgrid using the following code: Hs = (sgrid + 0.9) ./ ((sgrid + 0.9).^2 + (4*pi)^2); (b) Use the 'mesh' command to plot the magnitude and phase (two plots) of Hs. The x-axis will be sigma (the real part of s), the y-axis will be omega (the imaginary part of s), and the z-axis (the "height" of the 3D plot) will be the magnitude or phase of H(s) at the corresponding value of s. Note that the command abs(Hs) returns the magnitude of Hs, and the command angle(Hs) returns the phase. Use the following code to plot, but put each of these in its own figure window. mesh(sigma, omega, abs(Hs)); mesh(sigma, omega, angle(Hs)); (c) Questions: (1) Using the transfer function H(s), determine the poles and zeros of the system. (2) Verify that your magnitude plot shows the poles of the system. Don't worry about seeing the zeros, it will be difficult since the plot is on a very large scale. Is it obvious where the poles are? (3) Orient your plots so that the omega axis is facing you, and the magnitude and phase are going vertical. As a function of omega, is phase an even/odd function? is magnitude an even/odd function? == Solution Image == {{:380matlab:ch3:uy3_1_1_si_1_sdomaingraphsolutionimage1.jpg?400|}} {{:380matlab:ch3:uy3_1_1_si_2_sdomaingraphsolutionimage2.jpg?400|}} == Solution == {{:380matlab:solutions:ch3:uy3_1_1_sdomaingraph.m|}} ====== Chapter 4: Applications of the Laplace Transform ====== ===== Op-Amp Circuits ===== ==== MT380.UY4.5.1 ==== Op-amp circuit Objective: Graph the s-domain transfer function of an op-amp circuit to gain a visual understanding. Also, create a bode plot of the s-domain transfer function with sigma = 0 to introduce students to the Fourier domain. Functions to learn: none Exercise: Using the figure associated with the problem, solve for the op-amp's transfer function H(s). The values are shown below L = 10e-2; % inductor C = 10e-6; % capacitor R1 = 10; R2 = 2000; a)Evaluate H(s) at all possible combinations of s with s = sigma + j*omega, sigma = -20:1:20, and omega = 2*pi*(-1000:10:1000). b) graph the magnitude of H(s) with the frequency axis in units of Hz. c) where are the zeros ? d) Evaluate H(s) with sigma = 0, and omega = 2*pi*(-1000:10:1000). Graph the magnitude in terms of decibels, dB = 20log(Vo/Vin); and with the frequency axis in units of Hz. Also, graph the phase in terms of degrees on the same plot. e) In part d you plotted the bode plot of the Op-Amp. This is to introduce you into the Fourier domain == Figure == {{:380matlab:ch3:uy4_5_1_op_amp.png?400|}} == Solution == {{:380matlab:solutions:ch4:uy4_5_1_op_ampcircuit.m|}} ====== Chapter 5: Fourier Analysis Techniques ====== ===== Fourier Series Coefficients ===== ==== MT380.UY5.4.1 ==== UY5_4_1_ComputationOfFourierSeriesCoefficients Objective: Use MATLAB to calculate the Fourier series coefficients of a periodic, continuous signal. Once the Fourier series coefficients are calculated, compare the Fourier series expansion with a finite number of harmonics to the original signal. Commands: numerical integration using sum Exercise: The signal x(t) that we will be using in this problem is periodic, with a period of 5 seconds. However, we will only be representing a single period of this periodic signal x(t) in Matlab. Use the following code to generate an array containing a single period of the signal x(t) over the time interval from t = -1 to t = 4: T = 5; period of signal w = 2*pi/T; fundamental frequency time parameters t_beg = -1; t_end = t_beg + T; dt = .001; time step t = t_beg:dt:t_end-dt; time array (covers one period of x(t) creation of x(t) signal x = zeros(1,length(t)); for m = 1:length(t) if t(m) <0 x(m) = 0; elseif t(m) < 1 x(m) =1; elseif t(m) < 2 x(m) = t(m)-1; elseif t(m) < 2.5 x(m) = 3-t(m); elseif t(m) < 3 x(m) = -2+t(m); else x(m) = 4-t(m); end end plot x(t) to see what it looks like (over one period) close all; figure(1); plot(t, x); xlabel('time (s)'); ylabel('x(t)'); a) As discussed in class, any periodic continuous signal can be described by its Fourier series expansion. The Fourier series coefficients x_n describe how much of each harmonic of the fundamental frequency there is in the signal. The equation to calculate the Fourier series coefficients x_n is shown below: x_n = (1/To)*integral(x(t)*exp(-j*n*w*t)*dt ,0, To) (1) Use MATLAB to numerically compute the equation above to find the Fourier series coefficients for a fixed number of harmonics of the fundamental frequency. (Let's start with using the first 10 harmonics, and then you can experiment with adding more later. Note that if we use the first 10 harmonics, we are preserving 21 terms in the Fourier series expansion summation: the DC component x_0 and then x_-10 up to x_10.) The following code gets you most of the way there: set up our n index (the index of the Fourier Series coefficient x_n) highest_harmonic = 50; This is the highest harmonic we'll use n = -highest_harmonic:highest_harmonic; x_n = zeros(1,length(n)); Calculate the Fourier Series Coefficients x_n: x_n = (1/To)*integral(x(t)*exp(-j*n*w*t)*dt ,0, To) for k = 1:length(n) x_n(k) = INSERT CODE HERE TO DO THE NUMERICAL INTEGRATION USING 'sum' end b) Now let's find the Fourier series expansion approximation of x(t) using the x_n that we calculated. We'll call it x_FS(t). The infinite Fourier series expansion equation is: x(t) = sum(x_n*exp(j*n*w*t), -inf, +inf) (2) But we're only keeping 2*highest_harmonic+1 terms in the summation as follows: x_FS(t) = sum(x_n*exp(j*n*w*t), -highest_harmonic, +highest_harmonic) (3) Use MATLAB to create x_FS(t), the approximate of x(t), over one period. The code below will get you started: x_FS = zeros(1,length(x)); x_FS(t) will be our Fourier Series expansion approximation of x(t) Calculate x_FS(t), our Fourier Series expansion approximation of x(t) for k = 1:length(n) x_FS = x_FS + INSERT YOUR CODE HERE end c) Now let's plot x_FS(t) and compare it to x(t). Note that we are plotting only the real part of x_FS(t). This is because x_FS(t) will have a small imaginary part since it isn't a perfect representation of the actual (all real) signal x(t): figure(2); plot(t,real(x_FS),'r','LineWidth',2); hold on plot(t,x,'k:','LineWidth',2); title('Comparison of Fourier Series approximation of x(t) with x(t)'); legend('x_F_S(t)','x(t)'); xlabel('time (s)'); d) And finally, let's plot the Fourier series coefficients x_n using stem plots. Since the x_n are complex, we will create one stem plot showing the magnitude of x_n and one stem plot showing the phase of x_n: phase_x_n = angle(x_n)*180/pi; phase mag_x_n = abs(x_n); magnitude of Fourier series coefficients figure(3); subplot(1,2,1); stem(n, mag_x_n); title('Fourier Series Coefficients - Magnitude'); ylabel('|x_n|'); xlabel('n'); subplot(1,2,2); stem(n, phase_x_n); title('Fourier Series Coefficients - Phase'); ylabel('angle(x_n) (degrees)'); xlabel('n'); e) Questions: 1) What affect does the number of harmonics used have on the reconstruction of the original signal? For this question, try using different values of the variable highest_harmonic. Maybe try 5, 10 and 100. 2) The reconstructed signal experiences Gibbs Phenomena. What signal characteristic causes Gibbs Phenomena? == Figure == {{:380matlab:ch5:uy5_4_1_computationoffourierseriescoefficients_figure.jpg?400|}} == Solution Images == {{:380matlab:ch5:uy5_4_1_si_computationoffourierseriescoefficientssolutionimage.jpg?400|}} {{:380matlab:ch5:uy5_4_1_si2_computationoffourierseriescoefficientssolutionimage2.jpg?400|}} == Solution == {{:380matlab:solutions:ch5:uy5_4_1_fourierseriescoefficients.m|}} ===== Circuit Analysis with Fourier ===== ==== MT380.UY5.12.1 ==== UY5_12_1_Circuit Analysis Fourier Transform Objective: Use Fourier analysis to calculate the output to a system when stimulated. Exercise: This problem is based off of problem 5.61 that is in the book. You will still need to do part a by hand. The book will direct you to the needed figures. a) Derive the impulse function for the circuit. b) Create a function file that models the impulse function obtained in part a. c) Create another function file that models the input signal. Remember that A = 5mA and T = 3. d) Using the functions obtained in part b and c find and plot Vout. Hint convolve the functions. e) Take a moment to appreciate how much easier that problem became by using MATLAB compared to if you had to convolve the two functions by hand. == Solution == {{:380matlab:solutions:ch5:uy5_12_1_circuitanalysisfouriertransform.m|}} {{:380matlab:solutions:ch5:h_t.m|}} {{:380matlab:solutions:ch5:x_t.m|}} ====== Chapter 6: Applications of the Fourier Transform ====== ===== Butterworth Filters ===== ==== MT380.UY6.8.1 ==== UY6_3_1_Active Filters Objective: Model a second order butterworth filter and test it by passing signals at various frequencies through the filter. Since this is a design problem, and you have already covered butterworth filters in lab, less code will be given. functions: review roots Exercise: This exercise will be broken down into several parts part 1) Creating the impulse response of the low pass filter. You will model a second order lowpass filter using the butterworth design. The filter will be created in a separate function file because you will use this filter in another MATLAB exercise. The function will take two arguments: the corner frequency in Hz and the time-step (which we'll call "dt") used in creating the signals. The function will return an array that represents the impulse response of the 2nd order lowpass filter. It is crucial that the dt used in creating the filter is the same dt used to create the signal that it is filtering. Below is the function declaration that is placed in the function file that creates the filter. function [h] = h_t(fc,dt) Inside the function file you will need to convert the passed in frequency in Hertz into radians per second: wc = 2*pi*fc; Section 6-8 in your text book shows how to create the frequency response, H(jw), for a second order butterworth lowpass filter using the specified frequency. For aid, refer to Table 6-3 (page 285) to obtain the pole locations and coefficients a1 and a2 that you will use to design your butterworth filter. Look at example 6-10 in your book if you need more help. The frequency response should look like the equation below. (Remember that you can get the frequency response from H(s) by substituting s = jw.) H(jw) = wc^2/(s^2 + a1*wc*s + a2*wc^2) (1) The frequency response needs to be changed into its time equivalent impulse response h(t). This will be done in two steps. In the first step the transfer function needs to be changed into a second order differential equation. If you have forgotten how to do this, refer to section 2-7.3. The second order differential equation should be of the form below. d^2y/dt^2 + a1*dy/dt + a2*y(t) = b2*x(t) (2) - note, the coefficients a1 and a2 are not the same a1 and a2 coefficients in equation (1). The second step is to transform the second order differential equation into the impulse response h(t). Refer back to section 2-8. First calculate the poles (p1, and p2) as shown in equation 2.126 (pg 65) in your book. p = roots([1 a1 a2]); - note, a1 and a2 refer to the coefficients from equation (2). Second, form the time array that will be used. The time array will start from t =0s, increment by dt, and end at t = 10*tau. Tau is the time constant obtained from the roots. The roots have the basic form p = -sigma +/- j*wd where sigma represents the real part and wd represents the complex part. Sigma is related to your time constant, tau. tau = 1/sigma (3) So figure out what your sigma is, and insert in your Matlab code: tau = INSERT CODE HERE; this is the code Once tau is obtained, the time array can be created. time parameters time_beg = 0; time_end = 10*tau t = time_beg:dt:time_end - dt; time array The reason why we are using tau is because the impulse response in non-finite. However, its limit as t -> inf+ is 0. lim h(t) = 0 t -> inf+ After about 5*tau the value of h(t) (or h(5*tau)) is almost zero since the impulse response will decay exponentially. We are going to t = 10*tau just to make sure that the behavior response is completely captured. Also, since h(t) is causal, the time array can start at t =0. Now that we have the coefficients, roots, and the time array, h(t) can be constructed using equation 2.136 in your book (pg 66). To make things simple, below is a compilation of the above code that you will put in your function file that creates the low pass filter impulse response. function [h] = h_t(fc,dt) wc = 2*pi*fc; the corner frequency in rads/s a1_bwf and a2_bwf are coefficients obtained from table 6-3 in the book (pg 285). a1_bwf = INSERT VALUE; a2_bwf = INSERT VALUE; a1 and a2 are coefficients associated with the 2nd order differential equation. a1 = INSERT CODE; a2 = INSERT CODE; r = roots([1, a1, a2]); obtain the roots of the polynomial. b2 is a coefficient according to equation 2.121 in the book it represents b2 b2 = wc^2; Now let's set up the impulse response tau = INSERT CODE; time parameters time_beg = 0; time_end = 10*tau; t = time_beg:dt:time_end - dt; time array impulse response notice that the equation is simplifies because b1 = 0 h = INSERT CODE; part 2) Now that you have constructed your filter, let's test it. Copy the code below in a separate script and run the script. The script will call your function file h_t.m to create the lowpass filter and then pass three signals of amplitude 5 and frequencies: 100,400,4000 Hz. The script will then plot the filtered signals. fc = 400; corner frequency of filter(Hz) time parameters time_beg = 0; time_end_x = 0.1; dt = 1/100000; time step tx = time_beg:dt:time_end_x - dt; time array for signals impulse response of low pass filter h = h_t(fc,dt); th = time_beg:dt:length(h)*dt-dt; time array for filter frequencies for the two signals f1 = 100; f2 = 400; f3 = 4000; signals to be filtered x1 = 5*cos(2*pi*f1*tx); x2 = 5*cos(2*pi*f2*tx); x3 = 5*cos(2*pi*f3*tx); filtered signals y1 = conv(x1,h)*dt; y2 = conv(x2,h)*dt; y3 = conv(x3,h)*dt; time array for y ty = time_beg:dt:length(y1)*dt-dt; figure(1); plots the impulse response subplot(2,2,1); plot(th,h); xlabel('time (s)'); ylabel('amplitude'); title('impulse response h(t)'); plots the filtered signal with Hz = 100 subplot(2,2,2) plot(ty,y1); xlabel('time (s)'); ylabel('amplitude'); title('filtered signal with Hz = 100'); plots the filtered signal with Hz = 400 subplot(2,2,3) plot(ty,y2); xlabel('time (s)'); ylabel('amplitude'); title('filtered signal with Hz = 400'); plots the filtered signal with Hz = 4000 subplot(2,2,4) plot(ty,y3); xlabel('time (s)'); ylabel('amplitude'); title('filtered signal with Hz = 4000'); part 3) You will now analyse the plots to verify that your filter is working. a) You designed your filter with an fc of 400 Hz and a dc gain of 1. Verify that the filtered signal with f1 = 100 Hz still has a max amplitude of 5. b) Now look at the filtered signal with f2 = 400 Hz. Since the signal has the same frequency as the corner frequency what should its amplitude approximately be? d) Looking at the last signal with f3 = 4000 Hz, ignore the first part of the filtered signal (the transient response) which has an amplitude larger than the rest of the signal. f3 is ten times the frequency of f2, and you are using a second order filter. How much smaller should the amplitude of the signal with frequency 4000 Hz be compare to the signal with frequency 400 Hz? == Solution Image == {{:380matlab:ch6:uy6_3_1si_butterworthfilterssolutionimage.jpg?400|}} == Solution == {{:380matlab:solutions:ch6:uy6_8_1_butterworthfilter.m|}} \\ The function file below creates the impulse response of a second order lowpass filter. \\ {{:380matlab:solutions:ch6:h_t.m|}} \\ ===== Sampling Theorem ===== ==== MT380.UY6.12.1 ==== UY6_12_1 Sampling Theorem With Noise Objective: This exercise will expose students to the process of sampling data, and the effects of aliasing. Functions: fopen, fscan, fclose, interp, and sound. Caution: The command 'sound' will clip at +/- 1V. Before playing any array using the sound command, make sure that every value is within the range. You should not have to worry about it in this lab, but it is good to verify so that you develop the habit. You can easily verify the max and min value in your workspace by specifying it to show those values, or you can do it by inspection of the graphs that you will plot. Exercise: For this exercise you will be given a file to read in. You will take the file and read it in. This file contains an array of data that represents sound. You will filter, sample, and filter again to gain an idea on the effects of aliasing. This problem is broken down into 7 parts to help simplify and organize it. Some of the descriptions will be long, but read them because they will help save you time. Part 1) Get Data From the File You will use the command fopen to read in the file 'PianoDataWithNoise.txt'. Everything in the file represents data used in generating the sound; just open up the file in notepad to see its format. The file is long and its length is 5512500 which represents 5 seconds of data. i.e. 1102500 elements represents 1 second. When you read in the file, the format of the file is exponential notation. After you read the file, make sure you close it. Reading the file will take MATLAB about 15 seconds, but it only needs to read it in once. After you have the file's content stored in a variable, comment it out so it doesn't read it every time you run the program. The variable will remain in the workspace unless you clear it. Part 2) Creating the filter. Before you sample a signal you always want to filter it with a lowpass filter. In modulo MT380.UY6.8.1 you created a function file that generates a second order filter. Create a filter with a corner frequency of 500 Hz. Remember that the samples per second is 1102500, and the time interval/step used for creating the filter is 1/1102500. Part 3) Filter Data a) You will have two arrays of data. One array will not be pre-filtered and the other array will be filtered. Use the filter created in part 2 to filter the data obtained in part 1, and store the data in another array. b) Plot the non-pre-filtered data (NPFD) and the pre-filtered data (PFD) as a function of time. I have given you a function file called fft_plot. If you have not downloaded it from the wiki page, download it now. This function will plot the frequency magnitude spectrum. It takes in two parameters: the array containing the data, and the number of samples that represent a second. Ex fft_plot(NPFD, 1102500). Use this function to plot the frequency magnitude of NPFD and PFD. Magnify the plots so you can see the frequencies between 0 and 1500 Hz. Notice that each signal has a frequency at 1200 and 1500 Hz; However, the PFD frequency at 100 and 1500 Hz is much smaller because of the filter. These frequencies represents possible noise that can get into your system. Even though the noise shows up in the PFD signal, it is negligibly small, and can be ignored. Also, other than the noise, the other frequencies are about 500 Hz or below. Part 4) Sample Data a) The data that you have been given represents a continuous signal. In this part you will simulate sampling a continuous signal. Sample the NPFD signal at 1050 samples/second, and sample the PFD at 1050 sample/second and 500 samples/second. At this point you should have three different data arrays. 1) a pre-filtered signal sampled at 1050 samples/s or Hz 2) a pre-filtered signal sampled at 500 samples/s or Hz 3) a non pre-filtered signal sampled at 1050 sample/s or Hz b) These three signals will need to be modified in order to be played using the sound command because the sound card can only play data at specific frequencies, and they are 8000,11025, 22050, 44100, 48000, and 96000 Hz. To do this, you will need to interpolate the signals using the interp command. The signals sampled at 1050 Hz will need to be interpolated to have 22050 samples/s and the signal sampled at 500 Hz will need to be interpolated to have 8000 samples/s. Interpolating a signal does not change the frequencies. c) After you have interpolated it, plot all three signals as a function of time, and plot their frequency magnitude spectrum using the fft_plot function file provided. Compare and contrasts these plots to the plots created in part. d) Questions 1) After sampling the NPFD at 1050, what frequencies did the noise alias to? 2) Why are some of the aliased frequencies attenuated? 3) Why does the PFD signal sampled at 1050 Hz not suffer from aliasing? 4) Why did the PFD signal sampled at 500 Hz become really distorted? Part 5) Create Filters for sampled data For the signals sampled at 1050 Hz, create a filter will a cut off frequency of 500 Hz. For the signal sampled at 500 Hz, create a filter with a cut off frequency of 250 Hz using the filter function file you created in modulo MT380.UY6.8.1. Part 6) Filter the sampled data a) Filter the sampled data with their corresponding filters. b) Plot each signal as a function of time, and plot their frequency magnitude using the fft_plot function file provided. Remember to pass in the correct samples/s for each signal. c) Questions 1) Why is it important to filter after sampling data? 2) If you have a Nyquist sampling rate of 1000 Hz, what should the maximum corner frequency of the low pass filter be and why? Part 7) Play the data a) You now get to play your data using the sound command. Remember to make sure that all values in any signal array is within +/- 1. b) Play all signals with there corresponding frequencies. The signals sampled at 1050 Hz should be played at 22050 Hz, and the signal sampled at 500 Hz should be played at 8000 Hz c) Compare your results with the sound files provided on the wiki page. d) Questions 1) You should hear a low hum in the NPFD signal. It is most apparent at the beginning of the sound file. Why is it important to pre-filter before sampling. 2) You signal sampled at 500 Hz doesn't resemble the original sound. Why is it important to sample at a frequency greater than the Nyquist sampling rate? == Sound File == == fft_plot function file == {{:380matlab:solutions:ch6:fft_plot.m|}} == Original Music File == {{:380matlab:solutions:ch6:pianofinal_5s.wav|}} == Music Files After Sampling == This file was sampled at 1050 samples/sec, then filtered again with fc of 500 Hz. \\ You should be able to hear a low buz that is not in the original file \\ {{:380matlab:solutions:ch6:pianos1nff.wav|}} This file was first filtered with fc of 500 Hz, sampled at 1050 samples/sec, then filtered again with fc of 500 Hz. \\ This file should sound similar to the original just quieter due to filtering \\ {{:380matlab:solutions:ch6:pianos1ff.wav|}} This file was first filtered with fc of 500 Hz, sampled at 500 samples/sec, then filtered again with fc of 250 Hz.\\ This file shows the effects of aliasing. \\ {{:380matlab:solutions:ch6:pianos8ff.wav|}} ====== Chapter 7: Discrete-Time Signals and Systems ====== ===== D-Time Signal Functions ===== ==== MT380.UY7.2.1 ==== MT380.UY7.2.1 Discrete Time Signal Functions Objective: gain a visual understanding of the fundamental period. Functions: rat Exercise: Consider the causal signal x = cos(2*pi*f*t) with f = 5 Hz. a) Calculate the period of the signal. b) Write a function file that can calculate the fundamental period of a sampled signal for any period (To) and sample rate (Ts). This function must be able to indicate if a fundamental period does not exist. For example, in my program, if the fundamental period does not exist i return a 0. Hint: for this function file, consider using the rat command. c) You will sample signal x at the following three sample rates (Ts):0.035, 0.25, 0.01, and .01/pi. Using the function file created in part b, determine the fundamental period of signal x sampled at the three different sample rates. For example, if I were to sample signal x at a sample rate of 0.077 my fundamental period would be 200 samples. d) Sample the signal x at the indicated sample rates for the duration of 1 fundamental period. Continuing my example, my time duration would be (No*Ts) or 15.4 seconds. You should have three sampled signals of different time duration e) Approximate three continuous time signals (x) by sampling at a very high frequency and with a time duration equal to the fundamental period of all three sampling rates.Continuing my example, I would create a signal that approximates a continuous signal that begins at t = 0s and ends at t = 15.4s. f) Plot each sampled signal x with its corresponding continuous approximated signal. See solution graphs for a visual understanding. If the fundamental period does not exist, do not plot the data points, just indicate it in the plot's title. g) Questions: 1) For each sampled signal, how many periods of the continuous approximated signal fit within one fundamental period of the sampled signal? 2) How do the numbers obtained in part 1 relate to the book equation 7.21 (No = k*(To/Ts)) == Solution == {{:380matlab:solutions:ch7:mt380_uy7_2_1_discrete_time_signal_functions.m|}} {{:380matlab:solutions:ch7:fundamentalperiod.m|}} ===== D-Time Fourier Series ===== ==== MT380.UY7.13.1 ==== == Description == MT380_UY7_13_1_DiscreteTimeFourierSeries Objective: Create a Program that simulates the Discrete-Time Fourier Series of a sampled signal. Commands: none Exercise: The file, signal.txt, contains one period of a sampled periodic signal. a)Open the file, import the data, and plot the original signal. b) Find the expansion coefficients of the given signal. c) Plot the magnitude line spectrum and phase line spectrum of the imported signal using the expansion coefficients obtained in part b d) Reconstruct the original signal using the expansion coefficients and plot it. e) Question: What is the highest frequency in the signal? == Signal == == Solution Image == == Solution == {{:380matlab:solutions:ch7:mt380_uy7_13_1_discretetimefourierseries.m|}} ===== Discrete Fourier Transform ===== ==== MT380.UY7.15.1 ==== MT380 UY7 15 1 Discrete Fourier Transform Objective: The objective of the assignment is to create a function that can implement the Discrete Fourier Transform and plot it. We have provided all of the code except for the actual implementation of the DFT, but please go through the code carefully and try to understand it. Commands: Review Exercise: You are given a file that represents a non-periodic signal. In this exercise you will (1) window the signal, (2) apply the DFT to the windowed signal, (3) and then plot the signals. a) Begin by importing the signal file (file located below), and storing its contents to a variable. This variable will hold the data points of the non-periodic signal. % this command will open the file. You will need to replace % filename with the name of the file. fileID = fopen('filename.txt'); % the variable signal will contain the data points that represents the % signal. % the variable num_of_samples contains the number of data points in the % file. [original_signal,num_of_samples] = fscanf(fileID,'%e'); % transpose the signal original_signal = original_signal'; % Since you have what you need, you can close the file. fclose(fileID); % Below is the sample rate at which the original signal % was created, in samples/second. sample_rate = 1200; % plot the signal as a function of n(index) figure(1); subplot(2,2,1); plot(original_signal); title('original signal'); Notice how the bulk of the information is around the center of the plot, and that there is little non-zero data at other locations. We will window the signal to isolate the bulk of the data. This limits the number of sampled values, and thus reducing the number of computations. %the two parameters below will index 'original signal' around the bulk of %the information. n_window_lower = 320; n_window_upper = 2119; % window the original signal windowed_signal = original_signal(n_window_lower:n_window_upper); Plot the windowed_signal and notice the difference between the 'original_signal' and the 'windowed_signal'. % plot the windowed_signal as a function of n subplot(2,2,2); plot(windowed_signal); title('windowed signal'); b) Create a function file that will perform the DFT. A skeleton of the function file is provided below. I have created the function file to operate in two modes depending on how many parameters you pass into the function file. All you need to do is modify the code that is between the comment lines to actually implement the DFT summation equation. Comment line %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Xk, f] = myDFT(signal,sample_rate) % Xk is an array of the DFT complex coefficients % f is an array of the frequencies that correspond to Xk % signal is a variable that contains the data points of your signal % sample_rate is the rate at which the signal was samples/created % This function file supports two modes. The first mode is the generic DFT, % and is triggered when only yhe signal is passed to the function. Ex % myDFT(windowed_signal) % If both parameters are passed, then the function file will execute the % second mode % Ex. myDFT(windowned_signal, sample_rate); % Checks to see how many parameters were passed in. If only one parameter % is passed in, then the original DFT is calculated. if(nargin < 2) sample_rate = 1; end % The equation for the DFT is % Xk = sum(x[n]*exp(-1j*k*omega_o*n), from n = 0 to n = No-1) % No represents the number of sampled values No = length(signal); % fundamental frequency gamma_o = 2*pi/No; % n is an array that indexes the data points n = 1:length(signal); % Xn are the DFT complex coefficients Xk = zeros(1,length(signal)); % f is an array that will contain the frequencies (Hz) at which the DFT is % applied f = zeros(1,length(signal)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Student will only modify the code in between the comment lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % iterates through the harmonics of gamma_o to find the DFT complex % coefficients and the frequencies (Hz) for k = 0:No-1; Xk(k+1) = sum( %INSERT CODE HERE); f(k+1) = k*gamma_o/(2*pi)*sample_rate; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Do not modify the code below %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % In the DFT, the values ranged from 0 to pi represent the positive % frequencies while the values from pi to 2pi represent the negative % frequencies Since the magnitude of the DFT is an even function we can % get all of the necessary information from looking at the values of Xn % that represent the frequencies from 0 to pi if(nargin == 2) Xk = 2*Xk(1:round(length(Xk))/2); f = f(1:round(length(f))/2); end c) Call the function to perform the DFT on your signal passing in only your signal to myDFT function. Call the function again, but pass in both the signal and the sample rate. % call function to perform the DFT [Xn_mode1,f_mode1] = myDFT(windowed_signal); [Xn_mode2,f_mode2] = myDFT(windowed_signal,sample_rate); % plot the DFT transform subplot(2,2,3); plot(f_mode1,abs(Xn_mode1)); title('DFT mode1 complex coefficients'); xlabel('frequency (hz)'); ylabel('DFT complex coefficients magnitude'); subplot(2,2,4); plot(f_mode2,abs(Xn_mode2)); title('DFT mode 2 complex coefficients'); xlabel('frequency (hz)'); ylabel('DFT complex coefficients magnitude'); d) Question: 1) Looking at the DFT mode2 plot (where you pass in the sample rate, what frequencies are in the original signal. Hint: there are only four. Remember to submit your plots and code. == Signal File == {{:380matlab:ch7:mt380_uy7_15_1_s_discretefouriertransform_signal.txt|}} == Solution Image == {{:380matlab:ch7:mt380_uy7_15_1_si_discretefouriertransform_solutionimage.jpg?400|}} == Signal Creation == {{:380matlab:solutions:ch7:mt380_uy7_15_1_sc_discretefouriertransform_signalcreation.m|}} == Solution == DFT Function File. {{:380matlab:solutions:ch7:mt380_uy7_15_1_dft.m|}} Associated Script File {{:380matlab:solutions:ch7:mt380_uy7_15_1_discretefouriertransform.m|}} ====== Extras ====== ===== Music Harmonics ===== This MATLAB assignment explores the relationship between music and the frequency domain using the FFT. The file provided gives you a sample of a cello playing the note D3. An FFT of note D3 is taken in order to approximate Discrete Fourier Transform coefficients. These coefficients are used to approximate other notes pertaining to a Cello. All you need to do is download the file and explore what is provided to you. Note, all the code is provided to you; however, you are encouraged to modify it and play around. Be creative. As of now, you do not need to turn in anything in order to receive credit for this assignment. == File == {{:380matlab:funapplications:music_harmonics.zip|}}