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

```Preliminaries

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)
and run it by typing something
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

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)

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:

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