Matlab Assignments

# MATLAB Exercises

### Lab1 - Plotting magnitude and phase of complex functions of a real variable

Lab 1 should introduce students to MATLAB, m files, command window, workspace, arrays, multiplication, powers, exp, sum, component-wise operations, defining complex numbers, complex arrays, plot, abs, phase, for loop and repeated addition for computing sums.

Exercise 1 - Plot the magnitude and phase of the function $X(j \omega) = \frac{j \omega}{1+j \omega}$ as a function of $\omega$ for $-10 \pi \leq \omega \leq 10 \pi$

Exercise 2 - Plot the magnitude and phase of the function $x(t) = imag(3 - e^{(1-j2\pi)t})$ as a function of $t$ for $-10 \leq t \leq 10$

Exercise 3 - Plot the magnitude and phase of the function $X(j \omega) = e^{j 3 \omega} + e^{j 5 \omega}$ as a function of $\omega$ for $-10 \pi \leq \omega \leq 10 \pi$. For what values of $\omega$ is $X(j \omega) = 0$? Solve tihs problem analytically and verify your answer.

Exercise 4 - Consider the complex function $\displaystyle{x(t) = \sum_{k=0}^6 e^{j \frac{2 \pi k t}{8}}}$. Evaluate this sum in MATLAB either by using the sum function or using a for loop. Plot the real and imaginary parts of $x(t)$ versus $t$ for $-10 \leq t \leq 10$. Now, analytically determine a closed-form expression for $x(t)$. Plot this function in MATLAB and verify that both of them coincide.

What to turn in? - You should turn in your MATLAB code and plots for all the exercises. When plotting magnitude and phase, use subplot in MATLAB to plot them as two subplots. Use xlabel and ylabel to clearly label what you are plotting on the X and Y axes.

### Lab2 - Manipulating signals

This lab is meant to go with lectures that provide a basic introduction to signals such as continuous and discrete-time signals, operations on signals such as scaling of the time axis, shift in time. In this Lab, you will learn to use wavrecord to record an audio clip or look for an audio clip online in wav form and try to plot that signal as a function of time. You will simulate echo in an auditorium and listen to the sound.

Exercise 1 - Create your own audio file that is at least 4 seconds long (the exact time duration is not really important, but do not make the file too long). Use the wavrecord command in MATLAB and a sampling frequency of 10000 Hz. You can also try to find a wav file online. Here is one that I like samplewavfile . Use the sound command in MATLAB to play the sound. Make sure the recording is fine.

Exercise 2 - Using the wavread command, read the signal into a vector called x. Also read the sampling frequency in to a variable called Fs. Make sure you understand what this sampling frequency means. Plot the received signal as a function of time. Your time axis must have units in seconds.

Exercise 3 - Even though the audio signal is originally a continuous-time signal, we cannot represent continuous-time signals on a digitial computer directly. In this example, you have sampled the continuous-time signal and you are representing it as a discrete-time signal. Create the signal $x[2n]$ and play the signal by using the sound command and use the same sampling frequency as before. See that it sound shriller. Think about why this should be the case.

Exercise 4 - Sum of sinusiods of the same frequency. One of the useful results about sinusoids and cosinusoids is that sums of sinusoids and cosinusoids of the same frequency but with arbitrary amplitude and phase result in another sinsuoid or cosinusoid of the same frequency. In this exercise you will observe this using MATLAB. Let $x_1(t) = A_1 \sin(\omega t + \theta_1)$ and let $x_2(t) = A_2 \sin(\omega t + \theta_2)$. For $A_1 = 1, A_2 = \sqrt{3}, \omega = 4 \pi, \theta_1 = 0, \theta_2 = \pi/4$, plot $x_1(t)$ and $x_2(t)$ as a function of $t$ for $-10 \leq t \leq 10$ . Use 200 samples along the time axis, i.e., time should be incremented in steps of 0.01. Observe that each of the signals is periodic with the time period that you expect. Now plot the sum of the signals and observe that it is also a sinusoid with the exact same frequency (time period). Now, analytically find an expression for $x_3(t)$ in the form $B \sin(\omega t + \phi)$, i.e., find $B$ and $\phi$. Plot this signal for the same range of values of $t$ and observe that the signals coincide. Use the hold on command in MATLAB and plot the two signals in different colors. We can easily extend this result to see that a linear combination of any number of sinusoids of the same frequency results in another sinusoid of the same frequency.

Exercise 5 - Sum of sinusoids of different frequencies. Let us now consider the same example as above, except that we consider the case when the frequencies of sinusoids are different. Let $x_1(t) = A_1 \sin(\omega_1 t + \theta_1)$, $x_2(t) = A_2 \sin(\omega_2 t + \theta_2)$ and $x_3(t) = A_3 \sin(\omega_2 t + \theta_3)$. For $A_1 = 1, A_2 = 1, A_3 =1$, $\omega_1 = 4 \pi, \omega_2 = 2 \sqrt{5} \pi, \omega_3 = 2 \sqrt{3} \pi$, $\theta_1 = 0, \theta_2 = \pi/4, \theta_3 = \pi/6$, plot $x_1(t),x_2(t)$ and $x_3(t)$ as a function of $t$ for $-10 \leq t \leq 10$. Let $x_4(t)$ be the sum of the three signals. Is $x_4(t)$ periodic? Why or why not? Adding sinusoids of different frequencies arises in modeling wireless communication channels when the receiver or transmitter is moving. In wireless channels, typically the transmitting signals propagate via several reflections (similar to echo) and the net result is that if a sinusoid is transmitted, the received signal will be the sum of sinusoids with different amplitudes and phases. When there is relative mobility between the transmitter and the receiver, these reflected paths also undergo different Doppler shifts and, hence, their frequencies will be slightly different.

Exercise 6 - Consider the CT signal $x(t) = e^{-t} \ \sin(3t) U(t+2)$. Plot this signal as a function of $t$ for $-10 \leq t \leq 10$ where the t axis is incremented in steps of 0.01. Compute the energy of this signal by first using the symbolic toolbox in MATLAB analytically. Then use a Riemann sum approximation to computing the energy using MATLAB.

Exercise 7 - Consider the signal $x[n] = e^{-0.1n} \ U[n+20]$. Let $y[n] = \sum_{-\infty}^{n} x[m]$. Write a MATLAB program to compute $y[n]$ and plot it. Use the stem function to plot discrete-time signals. Notice that $y[n] = y[n-1]+x[n]$. Using this can you simplify your program? Read the help on the filter command and implement this using the filter command.

What to turn in? - Hard copy of programs for Exercises 1 - 4. Exercises 5-7 are compulsory for honors students

### Lab3 - Systems

Exercise 1 - Simulate a room with echo. Suppose this speech signal is played in a large room which results in significant echo. Let us model the echo as being composed of a direct path which has no delay and two reflected paths, one corresponding to a delay of $\tau$ seconds and another corresponding to a delay of $2 \tau$ seconds. Let the reflection coefficient for these two paths be 0.9 and 0.81 respectively. Simulate the received signal $y[n]$ and play it using the sound command. Play around with different values of $\tau$ and in your opinion, what is the smallest value of $\tau$ for which you can clearly start to hear some distortion in the signal?

Exercise 2 - Just for your interest - If you know the delays and the reflection coefficients exactly, can you think of a way to recover the signal $x[n]$ from $y[n]$?

Exercise 3 - (Optional) There is a simple way to recover $x[n]$ from $y[n]$. Unfortunately, in practice this simple technique will not work well since what you get to observe will always have some noise in addition to just the superimposed echoes. One way to model the received signal is to add a Gaussian noise with zero mean and variance $\sigma^2$ to the signal $y[n]$. For $\sigma^2 = 0.01$ simulate the received signal and run your recovery algorithm on this. How does the recovery algorithm perform?

### Lab 4 - Continuous-Time Fourier Series

Fourier's results says that a large class of periodic signals can be represented as linear combinations of sinusoids and cosinusoids or, equivalently, complex exponentials. More precisely, let $x(t)$ be a continuous-time periodic signal which satisfies some conditions known as the Dirichlet conditions (let us worry about these conditions now). Let the time period of the signal be $T$. We will refer to $\omega_0 = \frac{2\pi}{T}$ as the fundamental frequency of this signal. It turns out that the signal $x(t)$ can be expressed as

(1)
\begin{align} x(t) = \sum_{k=-\infty}^{\infty} X[k] e^{j k \omega_0 t} \end{align}

where

(2)
\begin{align} X[k] = \frac{1}{T} \int_0^T x(t) e^{j k \omega_0 t} \ dt \end{align}

In this lab, we will see this to be true for some periodic signals. Specifically, for a given set of $X[k]$, you will reconstruct $x(t)$ using Eq. (1)

Part 1 - First let us consider the case when

(3)
\begin{align} X[k] = \left\{ \begin{array}{ll} 1/2 & k = 0 \\ \frac{2 \sin(2\pi k T_1/T)}{2 \pi k}, & k \neq 0. \end{array} \right. \end{align}

with $T_1 = 1/4, T=1$. Write a progam in MATLAB to compute the signal $x(t)$ given in Eq. (1) for $-2 \leq t \leq 2$. It is not possible to implement infinite summations in MATLAB. Therefore, we will use an approximation given by

(4)
\begin{align} x_N(t) = \sum_{k=-N}^{N} X[k] e^{j k \omega_0 t} \end{align}

First set $N=1$ and plot $x_1(t)$. Using the hold on command in MATLAB, change $N$ to 2 and plot $x_2(t)$ on top $x_1(t)$. Repeat for $N = 3, 4, 5, 10, 20, 50$. What do you observe?

Part 2 - Now start a new figure and change $X[k]$ to $X[k] = \frac{1-e^{-4}}{4+j2 \pi k}$. Set $T=2$ and plot $x_N(t)$ for $-4 \leq t \leq 4$ for $N = 1, 2, 3, 5, 10, 20,50$. What do you observe?

Think about what you have accomplished and how this means that different signals can be written as linear combinations of complex exponentials. What questions come to your mind? Note these questions down so that we can talk about them on Tuesday.

Part 3 - For the honors students - By now you have probably realized that the signal in Part 1 is $x(t) = rect(t/2T_1)$ ($T_1 = 1/4$). Let $e_N(t)$ denote the approximation error defined by $e_N(t) = x(t)-x_N(t)$. Plot the energy in approximation error $e_N(t)$ as a function of $N$. What do you observe? Read Section 3.4 from the book and write a few words about whether your result shows convergence of the Fourier series.

What to turn in ? - MATLAB code and a plot of $x_N(t)$ for the different values of $N$ for Part 1 and Part 2. For each part, there should be only one plot and the different figures for different values of $N$ should be superposed on the same plot. Honors section students should turn in the code for Part 3 and the plot of the energy in the approximation error as a function of $N$.

### Lab 5 - Discrete-Time Fourier Series

Write a program in MATLAB that will compute the discrete-time Fourier series (DTFS) coefficients given a discrete-time signal $x[n]$ and its time period $N$. Using this program compute the DTFS coefficients of the signals in Problem 3.28 a (parts a and c only). Compare with the analytical expressions you got when working out the homework for parts a and c. Plot the magnitude and phase of $X[k]$ as a function of $k$ for $0 \leq k \leq N-1$. Use the stem function in MATLAB for the plot.

The assignment will be due on Tuesday along with the rest of the homework.
What to turn in ? - MATLAB code and a total of 4 plots (magnitude and phase of the Fourier series coefficients for part a and part c).

### Project - Detection of Signals (Gravity Waves) in Noise

You may have heard the exciting news that scientists were able to record gravity waves from the collision of two black holes. What a historic moment for science! In my opinion, it is as much a victory for engineering as it is for science, since the engineering that has gone behind the development of detectors to detect extremely weak signals and the algorithms that are needed to detect the weak signal in the presence of noise are remarkable! One way for people to get a feel for the gravity wave is to convert the gravity wave into a sound signal and such as signal is shown below.

The type of signal recorded is called a chirp signal. A chirp signal is a sinusoidal signal whose frequency linearly increases with time. Assume that a chirp signal is of length $T$ seconds and the frequency of the sinusoid increases from $\omega_1$ at the $t=0$ to $\omega_2$ at $t=T$. Express such a signal mathematically. Let $x_L(t)$ and $x_W(t)$ be the two signals recorded at the LIGO observatories in Livingston, Louisiana and Hangford, Washington. Mathematically, what is the relationship between the two signals?

Let $x(t)$ denote a 10 second long signal recorded by the LIGO detectors. A reasonable model for $x(t)$ is to assume that $x(t)$ contains three segments, i.e., $x(t) = [0_{\tau}(t) \ x_c(t) \ 0_{8-\tau}(t)] + n(t)$, where $x_c(t)$ is a chirp signal of duration 2 seconds, whose frequency increases linearly from 32 Hz to 256 Hz and whose amplitude is 1, $0_{\tau}(t)$ is a signal which is zero for $\tau$ seconds, and $n(t)$ is a noise signal. Let $x[n]$ denote the discre-time signal obtained by sampling $x(t)$ at a sampling rate of 10,000 samples per second. Create the vector $x[n]$ by assuming the sampled version of $n(t)$ to be Gaussian random variables with zero mean and a variance of 2.

When "looking for" gravity waves, a priori we do not know whether a chirp signal is present in $x(t)$ and if it is, we do not know when it starts within the 10 second span. How can you implement a detector that detects if a chirp is present in $x(t)$ and if it is when it starts? Read about correlation functions. See Problem 1.37 and 2.65 from Oppenheim, Wilsky and Young.

I would like to emphasize that the signal recorded in the LIGO detectors is *not* the audio signal created by the collision of black holes. The signal actually corresponds to the expansion and shrinking of distances between objects, i.e., warping of the space. We can nevertheless pretend that these signals are audio signals and listen" to them through speakers. This should also emphasize the utility of the tools we are learning in this course to model a variety of phenomena from gravity waves to audio signals!

### Project - Interference Rejection

Communication systems are often subject to interference from a variety of sources. In this project, you will use MATLAB to read a wav file, simulate the effect of narrowband interference and process the distorted signal to recover the original signal. You should write a MATLAB program that performs the following functions

1. (25 points) Read a speech file which is available in wav format. The file anykey.wav is available for download. But I encourage you to find your own wav file on the internet or you can also create your own speech file if you prefer. Use wavread function in MATLAB to read the file. Play the sound file through the speakers and observe that you can hear the speech file. Plot the magnitude of the Fourier transform of the speech file. The x-axis should be in rad/s.
2. (25 points) Simulate the effect of interference. We will model the interference as being composed of 10 cosinusoid functions at frequencies 4500,4550,..,4950 Hz whose amplitudes are all Gaussian random variables with zero mean and variance 1/10. Play this through the speakers and observe that you can see the effect of interference. Plot the magnitude of the Fourier transform of the speech file with interference. You should see that the interference is dominant.
3. (50 points) Design an LTI system to remove the interference by designing a filter that cuts off all frequencies greater than 4000 Hz. Process the received signal by passing it through a LTI system to mitigate the effect of interference. Play this through the speakers and observe that the effect of interference has been largely mitigated. Notice here that you are processing the received signal in the discrete-time domain. The conv function in MATLAB can be used to perform convolution. Another alternative is to use the filter command. Think about how you may have to truncate any impulse responses that are of infinite duration. Plot the magnitude of the Fourier transform of the processed signal as a function of $\omega$. The units for $\omega$ should be rad/s.

You can work in a team of two students. You should submit your MATLAB code and three plots of the magnitude of Fourier transforms plotted as a function of $\omega$.

The following commands in MATLAB will be useful. Do a help on these commands or read about them by searching on google. The syntax for writing a for loop will also be useful.
In this project, we are interested in plotting the magnitude of the Fourier transforms as a function of $\omega$. Normally, what we would like to get is the function $X(j \omega)$ for all values of $\omega$. However, on a digital computer, we cannot do this. We have to evaluate $X(j \omega)$ for a some finitely many values of $\omega$. A built in function in MATLAB called fft can be used to obtain the DTFT of a signal $x[n]$ of length $N$, namely $X(e^{j\Omega})$ for $N$ uniformly spaced values in the interval $[-\pi,\pi]$. Be sure to read the help section for the fft function and the ifft function. This will tell you why you need the fftshift and ifftshift functions.