Revision as of 12:50, 18 December 2008 by Mcwalker (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)


DT Fourier Series in Matlab with ONE Command

DT Fourier Series with a single MATLAB command!


Calculating fourier series by hand can often become time consuming and error prone. Matlab has an easy and fast built-in fuction for computing discrete time fourier series coefficents. Unfortunely, this wont help you on exams, but it might save you considerable time on homework assignemnts.

The command is ifft. It takes in a vector representing your signal and produces a vector of the fourier series coefficients. Two examples are provided below:

Example 1: The signal is represented by the graph below and is periodic for all time:

This signal can be represented by a vector. Each element in the vector corresponds to the value that the signal takes at each time interval. At time 0, the value is 2. At time 1 the value is 1 and for time 2 and 3 the value is 0. This can be represented by the vector below:

[2,1,0,0]?

To find the fourier series coefficients, we would use the following matlab code:

signal = [2,1,0,0]; fouriercoefs = ifft(signal)

The output gives:

fouriercoefs = 0.7500 0.5000 + 0.2500i 0.2500 0.5000 - 0.2500i

This means that the fourier series coefficients are: a0 = .75 a1 = .5+.25j a2 = .25 a3 = .5-.25j.

Example 2: The signal is represented by the graph below and is periodic for all time:

This signal can be represented by a vector like before. You may recognize this signal as x(t) = t for 0<=t<=2 in continuous time. Since we're in discrete time, the vector below represents the signal:

[0,1,2]?

Note that this is NOT the same as [0,.5,1,1.5,2]?. In that case, the continuous time example would be x(t) = .5t for 0<=t<=4. To find the fourier series coefficients, we would use the following matlab code:

signal = [0,1,2]; fouriercoefs = ifft(signal)

The output gives:

fouriercoefs = 1.0000 -0.5000 - 0.2887i -0.5000 + 0.2887i

This means that the fourier series coefficients are: a0 = 1 a1 = -.5 - .2887j a2 = -.5 + .2877j

Hopefully this will help you take DT fourier series in matlab easier and faster. These are examples which you can easily verify by hand.


comments:

... --shaun.p.greene.1, Wed, 26 Sep 2007 22:46:19 Good find in the functions library.

Although I'm not sure that the ifft() does what we really want, I tried it on the homework, and it gives me pretty much the same answer, so I definately think its close to what we need.

I found a function called dftmtx() that will generate the k and n matrix that is needed for finding the ak values. I apologize in advance because I'm not good in latex.

we have

$ \displaystyle a_k = \frac{1}{N}\cdot \sum_{k=0}^{k=N-1} \left( signal \cdot e^{(\jmath\frac{2\pi}{N} k n)} \right) $

the dftmtx() command should give you the matrix that has your exp(......).

Then, all you have left to find the aks is to

ak = signal * dfftmtx(N);

which gives a matrix of ak values that is almost identical to what the ifft() command gives you.

Thanks again for the good start with the ifft function, it got me moving on this homework.

... --shaun.p.greene.1, Wed, 26 Sep 2007 22:43:05 sry, that paragraph is really hard to read.

help? --andrew.c.daniel.1, Thu, 27 Sep 2007 00:04:44 if you used ifft() in the homework and it didn't work you could try transposing the [1xn]? vector wavread gives you to a [nx1]? horizontal vector matlab syntax: A_transposed = A';

... --tom.l.moffitt.1, Thu, 27 Sep 2007 00:15:22 If you do it on the homework, make sure you take it over one period. Since with voice recordings each period wont be exactly the same, it's a good idea to just do it over one.

A note about fft vs ifft --ross.a.howard.1, Thu, 27 Sep 2007 09:31:01 If you look at the help fft page it gives the equations it represents with fft and ifft. There are some differences between them and our book. fft has the correct sign in the complex exponential, but is not multiplied by 1/N. ifft finds the conjugate of the aks (which when you plot abs(ak) it does not matter) and has the 1/N term. For finding the aks using fft: aks=aks./length(aks); One should also note the helpful function shiftfft. It is useful because the matrix returned from fft starts at index 1 but contains a0, index 2 contains a1, and so on until it reaches the highest k, then it starts counting down. This means that when you plot the aks, they will not be in the right order. The shiftfft function correctly puts the negative aks on the left of the a0. (aks=shiftfft(aks);) Note: if you used ifft to find the aks, then use the shiftifft instead. Hope this helps. :)

Tutorial: Vector/Matrix Manipulation in Matlab

The Best Way

The best way to reverse a vector is probably as follows:

Say we have a vector,

test = [1,2,3,4,5]

and we want to reverse the numbers to provide us with a new vector:

newVector = [5,4,3,2,1]

To accomplish this without using a while loop or for loop, we can simply use the variable 'end', that is predefined in matlab:

newVector = test(end:-1:1)

What this statement say is take the last element of test ('end'), place that in the first element of newVector, then take the next last element of test ('end - 1') and put it into the second element of newVector, and so on, until newVector is a reversed copy of test. This method works with both row and column vectors.

Other Ways

Use the function flipud(vector) to reverse a column vector quickly.

>> flipud([1;2;3])
ans =
3
2
1
Use the function fliplr(vector) to reverse a row vector quickly.
>> fliplr([1,2,3])
ans =
3 2 1

Fourier in Matlab

This is some matlab code to check your fourier ak's. All you have to change is the expression for ak, and any particular ak's that are not in the formula (usually a0 is one of these). The other thing you must change is in the expression for sum, the last number (in this case 6) is the period (T). Note: This is what I got for problem 3.22(e).

t=-6:.001:6;
k=-40:40;
ak=0;
f=0;
sum=0;
ak=1./(j.*k.*pi).*(cos(2.*pi./3.*k)-cos(pi./3.*k)); ak(41)=0;
for num=-40:40
sum=ak(num+41).*exp(j*num*t*2*pi/6);
f=f+sum;
end plot(t,f)

Sound in Matlab

Recording sound:

WAVRECORD(N,FS,CH) Records N samples at frequency FS. CH is number of channels. Use 1 for the homeworks. Returns a matrix of size N, contiaining the samples
WAVPLAY(Y,FS) Plays audio stored in vector Y at a sampling rate of FS Hz.
Plot(Y) Plots vector Y (useful to visualize an audio signal). Plotting at different frequencies:
x=y ( 1 : N : length(y) ); will create a vector x, which is the signal Y at 1/N of its original frequency. (takes every N element of y and puts in in x)

Plotting in Matlab

plot(x,y,z)

  • Will plot X versus Y. If Y is omitted, the function will plot X versus size(X).
  • Z parameter will change the format of the plot. It is a character string. For example, Z = 'b*' will plot using blue stars. For all possible options, check the help page for plot

semilogx(x,y,z)

  • Creates plot with logarithmic scale for x

semilogy(x,y,z)

  • Creates plot with logarithmic scale for y

loglog(x,y,z)

  • Creates plot with logarithmic scales for both x and y

plot(x1,y1,z1,x2,y2,z2...xn,yn,zn)

  • Will plot multiple curves on one graph where X, Y and Z each correspond to one individual curve.

subplot(y,x,z)

  • Will allow you to plot many graphs at once. X and Y are the number of graphs in each row and column, respectively. Z is the current graph that the plot() command will affect.
  • To change which plot to modify, you must reuse the subplot command with the same X and Y, and change the Z value to the plot you wish to modify

xlabel('label'),ylabel('label')

  • Changes the label for the x and y axis on the current graph.

title('label')

  • Changes the title of the current graph.

grid on/ grid off/ grid minor

  • Turns a grid on or off. grid minor adds more lines to the plot.

Functions in Matlab

Writing functions in matlab can save you the trouble of typing the same lines of code over and over again. Functions are written in m-files. Open MATLAB, then go to file->new->m-file. Let's begin with an empty function named test that does nothing.

function test()

That was easy. Now, save the file as test.m. Every function you write must have the same name as the name of the m-file that contains it. A function foo() would be in foo.m, bar() would be in bar.m, and so on. Go to the matlab input and type test() to execute your function. Nothing happens, since we only made an empty function.

Now, let's try and make our function do something. Let's modify our file like so:

function test()

fprintf('Hello World\n')

Now save the function and execute it again. As you can see, anything that happens within your function will execute every time you call it from the main window. Whats the point of this? Let's give it an input, and get an output instead. We'll add two numbers:

function x = test(a,b)

x = a+b;

As you can see by running this code, you can give your functions inputs, and they will give you outputs. You can give your functions as many inputs as you want. In the prompt, try typing y = test(5,6). You will see that y is now equal to 11. The function will return the value of x as whatever it is equal to when the function finishes.

What if you want multiple outputs, you ask? Easy. Lets change our code to return a number of outputs. You can return an arbitrary number of outputs with your function, just like inputs.

function [x, y, z]? = test(a,b)

x = a+b;

y = a-b;

z = a*b;

Go back to the main prompt, and type [a,b,c]? = test(6,4). You will see that now a = 10, b = 2, and c = 24.

Alumni Liaison

Correspondence Chess Grandmaster and Purdue Alumni

Prof. Dan Fleetwood