(Added 2013.02.06: Here are a few hints: pdf)
In this activity we study the some basic operations on signals and their interaction with the discrete Fourier transform.
First, we will write a function that will reverse a signal. That is, we want a function that takes a signal as input and outputs the signal defined by . Select File>New>Function to pull up an editor window in which we can write our function. Enter the following:
Here, rev (for reverse) is the name of our function. The x in the brackets represents the vector that we’re trying to reverse. The y is the output of the function; this should be the reverse of x. The two lines in the middle describe how to compute y from x. (At this point, you should figure out why this accomplishes the task.) The end indicates the end of our function. File>Save this file as rev.m. Now we can use rev in the command window just like any other MATLAB function. Try it out on a few vectors x to make sure that it works. Verify the identity we proved in class: .
Now it’s your turn. Select File>New>Function again. Please write a function — call it shift — that cyclically shifts the entries of a vector m places to the right. That is, given a signal , we want to calculate the signal defined by . Here, we always understand modulo , meaning that if , then means . The function shift that you are to write takes two inputs — x and m, so the first line of your function should be y=shift(x,m). Its should be end. As for how to actually get y from x, I want you to figure it out using the rev function as a guide.
The function you just wrote shifts components of a signal m places to the right. We’d also like the flexibility to shift components m units to the left. We could do one of two things. We could write another function called leftshift or something that does this, but that would be wasteful. Shifting each component of an N element signal m positions to the left (cycically!) is the same as shifting each component of our signal ??? units to the right. (What is ???). So you can use your existing shifting function to do left shifts, too.
We proved in class that if , then . Verify this formula for various x and m. For example:
Now let’s write a convolution function. Again, select File>New>Function to open a new editor window. We will write a function conv that takes as input two signals x and y and outputs the convolution z of x and y. I would suggest starting out by initializing your eventual answer z and filling in the first entry:
z=zeros(1,size(x)); (double check that)
z(1) = sum(x.*rev(y));
You could then proceed with a for loop:
You are not obligated to code it as I suggest. Just don’t do it in such a way that if defeats the purpose of the next part.
To check that our conv works, we can use the relationship between convolution, pointwise multiplication under the Fourier transform: conv(x,y) and ifft(fft(x).*fft(y)) should be equal. Check this for a few vectors x and y.
Now that we have a nice functional convolution function, we can use it to investigate some signals. First, let’s see what the effect of convolution in the frequency domain on a signal in the time domain. Let t=0:0.02:0.98; and let x=cos(2*pi*4*t)+3*cos(2*pi*10*t);. Have a look at x and fft(x) using a stem plot and explain why it looks how it looks. Let’s kill one of the peaks as follows. Set w=((1:50)<8)+((1:50)>42); and let xx=conv(x,ifft(w)). (This should give you the same vector as xx=ifft(fft(x).*w);.) Make a stem plot of xx. What function does the plot resemble? More precisely, what simple function, when sampled 50 times per second, gives xx? Check your guess (using MATLAB, of course). Explain what’s going on in a few sentences.
Please submit the following along as part of Assignment 2:
- your code for the shift and conv functions
- a few examples verifying that your code works
- your stem plots of x and xx and your few sentences of explanation.
Optional/time permitting: We just looked at an example of pointwise multiplication in the frequency domain and the corresponding convolution in the time domain. Let’s consider, on the other hand, pointwise mutliplication in the frequency domain. Set y=sin(2*pi*4*t); and v=(1:50)<26;. Have a look at stem(t,y) and stem(t,y.*v) as well as at stem(fft(y)) and stem(abs(fft(y.*v))). Since the Fourier transform maps pointwise multiplication to convolution, fft(y.*v)=conv(fft(y),fft(v)). The signal y.*v is a “zero padded” sine wave. Zero padding in the time domain corresponds to convolution with fft(v) in the frequency domain. Functions like fft(v) are important and we’ll be seeing them more later.