Programming in MATLAB; Periodic Patterns; Fourier Transforms

ccollins@math.utk.eduPreliminariesStart MATLAB (from Start Menu). Arrange your windows so you can see MATLAB and this webpage, or just print out this webpage.The Basic Plan1. 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 (assomething.m) clear your variables, and run it by typingsomething4. 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 ScriptHere'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 FunctionThe 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 asupdownf.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 PolynomialsYour assignment is to determine the Fourier Transform of some data.Basic ElementsWe will do all our work in the interval [0,2p] Start by creating a vector ofx-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 FunctionsAny function which can be built from these seven functions is defined by the seven coefficients you use. For example the function1 + 2cos(3x) - sin(2x)is defined by the coefficient list: (1, 0, 0, 2, 0, -1, 0). These are the coefficients of1, cos(x), cos(2x), cos(3x), sin(x), sin(2x), sin(3x)in order. Using the blockall7we 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)AssignmentCopy the following MATLAB code into a m-file and save it asbwave.min the folderwork. % 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 inYHints: 1. Once you think you know one part, subtract it fromYand work with the rest. For example, if you think the vector contains0.2 cos(x), setnY = Y - 0.2*c1;and plot and work withnY. 2. Compute the average value of all the vector. (Use the MATLAB functionmean). 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 valuesc0*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 sectionOptional Stuff on Orthogonality (or How to get the answer fast)Compute the dot products of the vectorsc0, 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, 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. 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 OrthogonalityBefore 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 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 VectorNow, use the same procedure we did above to determine the coefficients for your mystery vector 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)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 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 itz = fft(y(1001:1500)); N = length(z); plot(abs(z(1:N/2+1)))This graph shows thefrequenciespresent in this sound sample.