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

  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.

  Try this:

      x1 = [1 0 0];
      x2 = [0 1 0];
      x3 = [0 0 1];

      x1*x2'
      x1*x3'
      x2*x3'
      x1*x1'

  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, ...

      c0*c1'
      c0*c2'
      c0*c3'
      c0*s1'
      etc.

  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.

  The end result is that these vectors are mutually orthogonal.

  Note the values on the diagonal, we will need them later.

Using Orthogonality

  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.

Mystery Function

  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
  code word).

  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)

  Mail me ccollins@math.utk.edu the name or code you used and your guess for the 
  coefficients BEFORE YOU LEAVE TODAY.

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 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).

   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.


   % 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(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