The final project for my Mechanical Engineering Analysis class was to apply learned numerical methods to engineering systems using Matlab. The final assignment submission was the video made by our Matlab script, and the code as well. Below is the final video, the full script and a little bit about the project.
I used Matlab to see how radio waves decay through space. I put together a small transmitter and connected it to a 12V bench power supply. I then played a 423 Hz tone (A note) through the radio, which would be easy to curve-fit using an r-squared method.
The data was recorded with a microphone at various distances, imported, and then the script turned it into a vector from the audio file using Matlab’s audioread function. I then used a loop to find the sine wave that would best fit the audio vector.
I added code for a graph at each distance, with both the sine wave of best fit, overlayed with the original audio data. I then graphed the least-squares value for all distances to see which had the lowest deviation, and hence, the clearest sound. Unsurprisingly, it was the closest distance to the transmitter.
I decided to make a simple stop motion video as well, to add some creativity to an otherwise dry presentation. I used Matlab’s imread function to import some drawings into the media library. I then made them display for only a fraction of a second, or for one frame, as Matlab’s frame rate is fairly slow. The result was a very fun animation of a radio wave transmitting from a tower. I then displayed the graphs for about three seconds each, and then made another small stop motion video for the closing. Below is the full code. As a note, no AI was used to write this code, and all drawings are original and were done by John (Jack) Kennedy.
Full Code:
clc close all clear all %% Setting up variables phiA=linspace(0,360,180); %adding a phase shift variable for the loop ampA=linspace(0.04,0.07,10);% Making an ampltude variable for the loop w=432*2*pi;% Adding the frequency used to the line of best fit equation and converting to radians. %% Distance 1 minerrora=1e5;% Making a large minimum error value to start with [a,Fa]=audioread('Distance 1.m4a',[2500,7000]);% Importing the file received at the first distance t=[0:length(a)-1]*1/Fa;% Making a vector for the time variable to have something to plot the sound values against AA=transpose(a);% Transposes into a more useable row vector for i=1:10% Going through the number of amplitudes & making a loop for using the brute force method for finding a line of best fit for j=1:180% Making a loop for the phase shift variable sinfitA=ampA(i)*sin(w*t+phiA(j));% Equation for a sine wave ra=sinfitA-AA;% Calculate R2 value; R2A(i,j)=sum(ra.^2);% Squaring the r value for positive number if R2A(i,j)◁minerrora% Seeing if the r-squared is less than the last one ampbestA=ampA(i);% Save if better phibestA=phiA(j);% Save if better minerrora=R2A(i,j);% Save if better end% End if statement end% End loop end% End loop ybestA=ampbestA*sin(w*t+phibestA);% this value makes the line of best fit with the best values selected from the previous loop figure(1) % Adding a figure plot to make the plots separate plot(t,a,t,ybestA)% Plotting the line of best fit and the original audio data against time xlabel('Time (seconds)')% Adding an xlabel ylabel('Ampiltude (Volts)')% Adding a y label title('Signal at 5 feet')% Adding a title legend('Received Signal','Line of Best Fit')% Adding a legend F1=getframe(gcf);% Gets the frame for the movie part R2StaticA=sum((ybestA-AA).^2);% Finds the final r-squared value for the final data plot%% Distance 2 minerrorb=1e5;% big boy [b,Fb]=audioread('Distance 2.m4a',[2500,7000]);% Importing the file received at the second distance BB=transpose(b);% Taking the transpose of b to make a row vector for i=1:10% Going through the number of amplitudes & making a loop for using the brute force method for finding a line of best fit for j=1:180% Making a loop for the phase shift variable sinfitB=ampA(i)*sin(w*t+phiA(j));% Get sine wave of best fit rb=sinfitB-BB;% Calculate R2 value; R2B(i,j)=sum(rb.^2);% Squaring the r value to get a positive number if R2B(i,j)◁minerrorb% Seeing if the r-squared is less than the last one ampbestB=ampA(i);% Save if better phibestB=phiA(j);% Save if better minerrorb=R2B(i,j);% Save if better end% End if statement end% End loop end% End loop ybestB=ampbestB*sin(w*t+phibestB);% this value makes the line of best fit with the best values selected from the previous loop figure(2)% Adding a figure plot to make the plots separate plot(t,b,t,ybestB)% Plotting the line of best fit and the original audio data against time xlabel('Time (s)')% Adding a label for x ylabel('Amplitude (Volts)')% Adding a label for y title('Signal at 10 feet')% Adding a title legend('Received Signal','Line of Best Fit') F2=getframe(gcf);% Gets the frame for the movie part R2StaticB=sum((ybestB-BB).^2);% Finds the final r-squared value for the final data plot %% Distance 3 minerrorc=1e5;% big [c,Fc]=audioread('Distance 3.m4a',[2500,7000]);%Importing the audio file received at the next distance CC=transpose(c);% The audio file is read as a column vector, so this make is a more useable row vector for i=1:10% Going through the number of amplitudes & making a loop for using the brute force method for finding a line of best fit for j=1:180% Making a loop for the phase shift variable sinfitC=ampA(i)*sin(w*t+phiA(j));% Get sine wave of best fit rc=sinfitC-CC;% Calculate R2 value; R2C(i,j)=sum(rc.^2);% Squaring the r value to get a positive number if R2C(i,j)◁minerrorc% Seeing if the r-squared is less than the last one ampbestC=ampA(i);% Save if better phibestC=phiA(j);% Save if better minerrorc=R2C(i,j);% Save if better end% End if statement end% End loop end% End loop ybestC=ampbestC*sin(w*t+phibestC);% this value makes the line of best fit with the best values selected from the previous loop figure(3)% Making a new figure plot(t,c,t,ybestC)% Plotting the line of best fit and the original audio data against time xlabel('Time (s)')% Adding a label for x ylabel('Amplitude (Volts)')% Adding a label for y title('Signal at 15 feet')% Adding a title legend('Received Signal','Line of Best Fit') F3=getframe(gcf);% Takes the plot as a frame for the movie part R2StaticC=sum((ybestC-CC).^2);% Finds the final r-squared value for the final data plot %% Distance 4 minerrord=1e5;% big [d,Fd]=audioread('Distance 4.m4a',[2500,7000]);% Importing the audio file received at the next distance DD=transpose(d);% The audio file is read as a column vector, so this make is a more useable row vector for i=1:10% Going through the number of amplitudes & making a loop for using the brute force method for finding a line of best fit for j=1:180% Making a loop for the phase shift variable sinfitD=ampA(i)*sin(w*t+phiA(j));% Get sine wave of best fit rd=sinfitD-DD;% Calculate R2 value; R2D(i,j)=sum(rd.^2);% Squaring the r value to get a positive number if R2D(i,j)◁minerrord% Seeing if the r-squared is less than the last one ampbestD=ampA(i);% Save if better phibestD=phiA(j);% Save if better minerrord=R2D(i,j);% Save if better end% End if statement end% End loop end% End loop ybestD=ampbestD*sin(w*t+phibestD);% this value makes the line of best fit with the best values selected from the previous loop figure(4)% Making a new figure plot(t,d,t,ybestD)% Plotting the line of best fit and the original audio data against time xlabel('Time (s)')% Adding a label for x ylabel('Amplitude (Volts)')% Adding a label for y title('Signal at 20 feet')% Adding a title legend('Received Signal','Line of Best Fit') F4=getframe(gcf);% Takes the plot as a frame for the movie part R2StaticD=sum((ybestD-DD).^2);% Finds the final r-squared value for the final data plot %% Distance 5 minerrore=1e5;% big [e,Fe]=audioread('Distance 5.m4a',[2500,7000]);% Importing the audio file received at the following distance EE=transpose(e);% The audio file is read as a column vector, so this make is a more useable row vector for i=1:10% Going through the number of amplitudes & making a loop for using the brute force method for finding a line of best fit for j=1:180% Making a loop for the phase shift variable sinfitE=ampA(i)*sin(w*t+phiA(j));% Get sine wave of best fit re=sinfitE-EE;% Calculate R2 value; R2E(i,j)=sum(re.^2);% Squaring the r value to get a positive number if R2E(i,j)◁minerrore% Seeing if the r-squared is less than the last one ampbestE=ampA(i);% Save if better phibestE=phiA(j);% Save if better minerrore=R2E(i,j);% Save if better end% End if statement end% End loop end% End loop ybestE=ampbestE*sin(w*t+phibestE);% this value makes the line of best fit with the best values selected from the previous loop figure(5)% Making a new figure plot(t,e,t,ybestE)% Plotting the line of best fit and the original audio data against time xlabel('Time (s)')% Adding a label for x ylabel('Amplitude (Volts)')% Adding a label for y title('Signal at 25 feet')% Adding a title legend('Received Signal','Line of Best Fit') F5=getframe(gcf);% Takes the plot as a frame for the movie part R2StaticE=sum((ybestE-EE).^2);% Finds the final r-squared value for the final data plot %% Distance 6 minerrorf=1e5; [f,Ff]=audioread('Distance 6.m4a',[2500,7000]);% Importing the audio file received at the final distance FF=transpose(f);% The audio file is read as a column vector, so this make is a more useable row vector for i=1:10% Going through the number of amplitudes & making a loop for using the brute force method for finding a line of best fit for j=1:180% Making a loop for the phase shift variable sinfitF=ampA(i)*sin(w*t+phiA(j));% Get sine wave of best fit rf=sinfitF-FF;% Calculate R2 value; R2F(i,j)=sum(rf.^2);% Squaring the r value to get a positive number if R2F(i,j)◁minerrorf% Seeing if the r-squared is less than the last one ampbestF=ampA(i);% Save if better phibestF=phiA(j);% Save if better minerrorf=R2F(i,j);% Save if better end% End if statement end% End loop end% End loop ybestF=ampbestF*sin(w*t+phibestF);% this value makes the line of best fit with the best values selected from the previous loop figure(6)% Making a new figure plot(t,f,t,ybestF)% Plotting the line of best fit and the original audio data against time xlabel('Time (seconds)')% Adding a label for x ylabel('Amplitude (Volts)')% Adding a label for y title('Signal at 30 feet')% Adding a title legend('Received Signal','Line of Best Fit') F6=getframe(gcf);% Takes the plot as a frame for the movie part R2StaticF=sum((ybestF-FF).^2);% Finds the final r-squared value for the final data plot %% Original transmitted signal [y,Fr]=audioread('432.mp3',[2500,7000]);% Importing the original transmitted audio file figure(7)%Making a new figure y1=(y(:,1));% Taking data from only one chanel(mono) plot(y1)%plotting the original mono signal xlabel('Sample Number')% Adding an x label ylabel('Amplitude (volts)')% Adding a y label title('Original Mono signal (432 Hertz)')% Adding a title F7=getframe(gcf);% Getting the final frame for the video %% Finding the r-squared values Res=[R2StaticF,R2StaticE,R2StaticD,R2StaticC,R2StaticB,R2StaticA];% Making a vector in order to plot the residuals Dis=[5,10,15,20,25,30];% Making a distance vector figure(8) plot(Dis,Res,'r*')% Making a plot for these datapoints hold on xlabel('Distance from Antenna (ft)') ylabel('Clarity') title('Distance from Antenna vs Amount of Static') [yax,index]=max(Res);%this line finds the minimum value of fx and assigs it to two variabels xmaz=Dis(index);%this line indexes the xvalue that corresponds to minimum fx value plot(xmaz,yax,'bo')%this line then plots the minimum values with a red circle. legend('Residual of each distance','Loud and Clear') F8=getframe(gcf); %% Video part Im1=imread('frame1.png'); % Read image into MATLAB Im1=imresize(Im1,[420,560]);% Resize image Frame1=im2frame(Im1);% Making the first frame Im2=imread('frame2.png'); % Read image into MATLAB Im2=imresize(Im2,[420,560]);% Resize image Frame2=im2frame(Im2);% Making the second frame Im3=imread('frame3.png'); % Read image into MATLAB Im3=imresize(Im3,[420,560]);% Resize image Frame3=im2frame(Im3);% Making the third frame Im4=imread('frame4.png'); % Read image into MATLAB Im4=imresize(Im4,[420,560]);% Resize image Frame4=im2frame(Im4);% Making the fourth frame cmap=gray(256);% Because the final frames were scanned in black and white, the colormap was empty Im5=imread('endframe1.png'); % Read image into MATLAB Im5=imresize(Im5,[420,560]);% Resize image endframe1=im2frame(Im5,cmap);%This adds a colormap to the image so it can be displayed in the video Im6=imread('endframe2.png'); % Read image into MATLAB Im6=imresize(Im6,[420,560]);% Resize image endframe2=im2frame(Im6,cmap);%This adds a colormap to the image so it can be dsplayed in the vid v=VideoWriter('Final Project','MPEG-4');%Initializes the videowriter to make the video open(v)% Begins the video for i=1:80 % does these four frames 80 times, making up roughly 10 seconds of a very simple animation writeVideo(v,Frame1)% Frame 1 of the intro page writeVideo(v,Frame2)% Frame 2 of the inro page writeVideo(v,Frame3)% Frame 3 of the introduction writeVideo(v,Frame4)% Frame 4 of the introduction end% Ends the introduction for i=1:110% Writing the video for about 3 seconds writeVideo(v,F1)% write the video with the first frame end% Ends this part of the video for i=1:110% Writing the video for about 3 seconds writeVideo(v,F2)% write the video with the second frame end for i=1:110% Writing the video for about 3 seconds writeVideo(v,F3)% Writing the video with the third plot end for i=1:110% Writing the video for about 3 seconds writeVideo(v,F4)% Writing the video with the fourth plot end for i=1:110% Writing the video for about 3 seconds writeVideo(v,F5)% Writing the video with the fifth plot end for i=1:110% Writing the video for about 3 seconds writeVideo(v,F6)% Writing the video with the sixth plot end for i=1:110% Writing the video for about 3 seconds writeVideo(v,F7)% Writing the video with the eventh plot end for i=1:110% Writing the video for about 3 seconds writeVideo(v,F8)% Writing the video with the eigth plot end for i=1:40% Writing the video for about 10 seconds because each of these images loop for 3 frames because I got lazy on the ending animation and didn't want to draw 4 frames :( for j=1:3% loops this image for 3 frames writeVideo(v,endframe1)% Adds the image end% ends this loop for j=1:3% loops this image for 3 frames writeVideo(v,endframe2)% Addsthe image end% Ends the loop end% Ends the lazy loop close(v)% Ends video