GS 2003/Math - Collins
Lab 2
Programming in MATLAB; Periodic Patterns; Fourier Transforms

Preliminaries

Start MATLAB (from Start Menu).
Arrange your windows so you can see MATLAB and
this webpage, or just print out this webpage.

The Basic Plan

1. Go under the file menu and choose New M-File.
   This will bring up a new window which you can type in.
   Move/resize windows so you can see all two/three at once.
2. Create the variables and try different commands in the command window.
   When you find something that works, copy it into your M-file.
3. At times, save your M-file (as something.m)
   clear your variables,
   and run it by typing something
4. Fix any problems (error messages or bad answers),
   keep adding features, and repeat steps 2-4.
5. For final test, run through all different types
   of inputs to see if your program performs as you expected.

First Script

Here's the script from class today. Copy and paste it
into the M-File editor window.  Save it as updown.m

x(1) = input(‘Enter initial value:’);
k = 1;
while (x(k) ~= 1)
   if (mod(x(k),2) == 0)
      x(k+1) = x(k)/2;
   else
      x(k+1) = 3*x(k)+1;
   end
   k = k + 1;
end


Run it by typing updown in the
command window.  After it is done, type disp(x)
to see the values computed.

First Function

The main two differences between a script and a function are that
(1) all the calculations and variables in a function are only available
    in the function; while in the script they are available in the main
    part of MATLAB.  This is called the scope of the variables.
(2) because of the scope, functions start with a line telling which
    variables will be passed in and which will be passed out of the
    function

Modify updown.m by replacing the first line by these three lines:
    
    function x = updownf(a)
    % computes the up/down sequence
    x(1) = a;
    
Save the result as updownf.m (The name of the file
must match the name of the function)

Now try typing updownf(11) in the command window.


Assignment (Part 1) 1. Which starting value from 1 to 20 produces the longest sequence before it gets back to 1? 2. Challenge: which starting value from 1 to 500 produces the largest value in its sequence?
Trigonometric Polynomials Your assignment is to determine the Fourier Transform of some data. Basic Elements We will do all our work in the interval [0,2p] Start by creating a vector of x-values in this interval: x = 2*pi*[0:399]/400; Our functions are: 1, cos(x), cos(2x), cos(3x), sin(x), sin(2x), sin(3x) c0 = cos(0*x); c1 = cos(1*x); c2 = cos(2*x); c3 = cos(3*x); s1 = sin(1*x); s2 = sin(2*x); s3 = sin(3*x); Plot these functions: plot(x,c1,x,c2,'--',x,c3,':') axis([0, 2*pi, -1 1]) title('cos(x), cos(2x), cos(3x)') plot(x,s1,x,s2,'--',x,s3,':') axis([0, 2*pi, -1 1]) title('sin(x), sin(2x), sin(3x)') To make further computations easier, let's combine all these functions into one block: all7 = [c0;c1;c2;c3;s1;s2;s3]; Building Functions Any function which can be built from these seven functions is defined by the seven coefficients you use. For example the function 1 + 2cos(3x) - sin(2x) is defined by the coefficient list: (1, 0, 0, 2, 0, -1, 0). These are the coefficients of 1, cos(x), cos(2x), cos(3x), sin(x), sin(2x), sin(3x) in order. Using the block all7 we can easily evaluate and plot this function: c = [1, 0, 0, 2, 0, -1, 0]; y = c*all7; plot(x,y) Try your hand and creating and plotting the following functions: a. -1 + cos(x) - sin(2x) b. cos(2x) + 1/2*cos(3x) - sin(2x) + 1/4*sin(3x) c. 1 + 3/5*cos(x) - 1/5*cos(2x) + 1/7*cos(3x) d. -sin(x) + 2*sin(2x) - 1/10*sin(3x) Assignment Copy the following MATLAB code into a m-file and save it as bwave.m in the folder work. % script: BWAVE.M - builds a random function name = input('Enter your name: ','s'); name = double(name); rand('state',round(sum(name)*pi)) sh = 64 + floor(60*rand(1)); name = [name,64+floor(60*rand(1,7))]; cz = (name(1:7)-sh); cz = cz/norm(cz); cz = round(cz*10)/10; xz = 2*pi*[0:399]/400; b7 = [ones(size(xz)) sin(xz) cos(2*xz) sin(3*xz) cos(xz) sin(2*xz) cos(3*xz)]; Y = cz*b7; plot(x,Y), axis([0 2*pi, min(Y)-0.1, max(Y)+0.1]) disp('Your mystery vector is stored in Y') clear cz b7 xz sh name % end of script Run this program and type in your name when prompted. Your assignment is to determine what set of coefficients were used to generate the vector with values in Y Hints: 1. Once you think you know one part, subtract it from Y and work with the rest. For example, if you think the vector contains 0.2 cos(x), set nY = Y - 0.2*c1; and plot and work with nY. 2. Compute the average value of all the vector. (Use the MATLAB function mean). 3. Count the peaks to find the highest frequency. 4. Look at the values at the endpoints and at p, and other interesting values. 5. Once you have an idea, you can use trial and error to figure out the possible coefficients. Try plotting your guess and the actual vector on the same graph: c = [ ]; % your guess y = c*all7; plot(x,y,x,Y) 6. All the coefficients are rounded to 1 decimal place and should be between -2 and 2. 7. Compute the values c0*Y', c1*Y', ..., and see if they tell you anything.
Assignment (Part 2) 3. Write out your name and the best set of coefficients you found by this trial and error method.
Do the Assignment above (#3) before you read this section Optional Stuff on Orthogonality (or How to get the answer fast) Compute the dot products of the vectors c0, c1, c2, c3, ... c0*c1' c0*c2' c0*c3' c0*s1' etc. The dot product of two vectors is proportional to the cosine of the angle between the two vectors. In particular if the dot product is 0, it means the angle is 90 degrees, i.e. the two vectors are perpindicular (or orthogonal). To compute all the possibilities at once, use all7*all7' This produces a 7x7 matrix, with the value in row i, column j is the dot product of the i-th vector with the j-th vector in the list c0, c1, c2, c3, s1, s2, s3. You should see that basically each vector is orthogonal to every other vector. Note the values on the diagonal, we will need them later. Using Orthogonality Before we work with our mystery vector, let's work with one where we know the composition: y = 0.1*c0 - 0.2*c1 + 0.3*c2 -0.4*c3 + 0.5*s1 - 0.6*s2 + 0.7*s3; Now, take each vector and form the dot product of it with y, then divide the result by the appropriate value from the diagonal we computed above. For example y*c0'/400 produces an answer of 0.1000 y*c1'/200 produces an answer of -0.2000 Continue with the rest of the vectors. Mystery Vector Now, use the same procedure we did above to determine the coefficients for your mystery vector When you get an answer, store the result in c and plot the result compared with the mystery function: c = [1, 0, 0, 2, 0, -1, 0]; y = c*all7; plot(x,y,x,Y) Advanced Processing We can simplify this process as follows: First, save the diagonal elements of the matrix all7*all7' D = diag(all7*all7')' Next, with any mystery vector y the following computes the coefficients: (Y*all7')./D Try it! If you like it, try it with different mystery functions. This process of using a standard set of orthogonal functions to decompose another function is called Fourier Transform. What we have been doing by hand, i.e. computing the coefficients, is called the Fourier Transform. Since it is so useful there is a special version of it called the Fast Fourier Transform (or FFT). Try this: z = fft(Y); N = length(z); a = 2*real(z(1:N/2+1))/N; a(1) = a(1)/2; b = 2*imag(z(1:N/2+1))/N; a(1:5) b(1:5) Normally FFT produces a vector containing complex numbers so we have to do extra work to convert these complex numbers into real numbers that we can understand. (That's what all the extra stuff is). Now, for fun, try this: load handel % this loads a sound file from Handel's Messiah sound(y) % this plays it (if you have speakers) plot(y) Let's grab a part of this sound and analyze it z = fft(y(1001:1500)); N = length(z); plot(abs(z(1:N/2+1))) This graph shows the frequencies present in this sound sample.
ccollins@math.utk.edu
Last Modified: June 12, 2003