Fourier Series

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.ccollins@math.utk.eduBasic Building BlocksWe 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];OrthogonalityOrthogonality is the mathematical generalization of perpindicularity. We say that two vectors areorthogonalif their dot product is zero. We say that two functions areorthogonalif 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 vectorsc0, c1, c2, c3, ...c0*c1' c0*c2' c0*c3' c0*s1'etc. To compute all the possibilities at once, useall7*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 listc0, 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 OrthogonalityBefore 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 withy, then divide the result by the appropriate value from the diagonal we computed above. For exampley*c0'/400produces an answer of 0.1000y*c1'/200produces an answer of -0.2000 Continue with the rest of the vectors.Mystery FunctionIf 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 incand 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 coefficientsBEFORE YOU LEAVE TODAY.Advanced ProcessingWe can simplify this process as follows: First, save the diagonal elements of the matrixall7*all7'D = diag(all7*all7')'Next, with any mystery vectorythe following computes the coefficients:(Y*all7')./DTry 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 itz = fft(y(1001:1500)); N = length(z); plot(abs(z(1:N/2+1)))This graph shows thefrequenciespresent 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