ECEn Matlab

 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).
 

 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?

A document to help gain an intuition behind convolution. convolutionguide.docx

The file below work together to create a MATLAB gui that simulates convolution

Please do not modify the following two files.
convolutiongui.m
convolutiongui.fig

You can modify the following three function files to change your signals x(t) and h(t)
x_t_gui.m This file modifies the signal x(t)
h_t_gui.m This file modifies the signal h(t)
time_parameters_gui.m This file returns the time parameters used for x(t) and h(t)

   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)?

 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.

 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)

Differential Equations-Euler's Method differential_equations_eulersmethod.docx

   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 Differential Equations-Euler's Method

 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

laplacetransformgui.fig laplacetransformgui.m

 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?

 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

 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?

 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. 

 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?

   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?

fft_plot.m

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

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

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.

   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))

   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?