This is a continuation of Lab 9. The main difference is
that instead of guesswork, we will develop a systematic
way to determine the composition of our mystery function.
Basic Building Blocks
We will do all our work in the interval [0,2p]
Copy and paste the following MATLAB code:
x = 2*pi*[0:399]/400;
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);
all7 = [c0;c1;c2;c3;s1;s2;s3];
Orthogonality is the mathematical generalization of perpindicularity.
We say that two vectors are orthogonal if their dot product
is zero. We say that two functions are orthogonal if their
inner product is zero.
x1 = [1 0 0];
x2 = [0 1 0];
x3 = [0 0 1];
These three vectors are mutually orthogonal. We knew this already
as they represent the directions of the x, y and z-axes, respectively.
Note that since in MATLAB these are 1x3 vectors, computing x1*x2',
multiplies a 1x3 with a 3x1 to get a 1x1, or scalar, result.
Now, look at the dot products of the vectors c0, c1, c2, c3, ...
To compute all the possibilities at once, use
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.
The end result is that these vectors are mutually orthogonal.
Note the values on the diagonal, we will need them later.
Before we work with our mystery function, 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.
If you don't have the script which produces the mystery function on
your computer, copy it from below and run it with your name (or any
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;
Mail me email@example.com the name or code you used and your guess for the
coefficients BEFORE YOU LEAVE TODAY.
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
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 Analysis. 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).
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;
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)
Let's grab a part of this sound and analyze it
z = fft(y(1001:1500));
N = length(z);
This graph shows the frequencies present in this sound
% script: BWAVE.M - builds a random function
name = input('Enter your name: ','s');
name = double(name);
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))
Y = cz*b7;
plot(xz,Y), axis([0 2*pi, min(Y)-0.1, max(Y)+0.1])
disp('Your mystery function is stored in Y')
clear cz b7 xz sh name
% end of script