% maatlab intro
% copyright 2020 by Thomas P. Weldon
% all rights reserved
clear all
close all hidden
clc

%some useful functions
syms u(t) pulse(t) triangle(t) 
sympref('HeavisideAtOrigin',1); %set equal to 1 as in unitstep
newHeaaviside0=heaviside(sym(0))
u(t)=heaviside(t) %unit step
pulse(t)=u(t+1/2)-u(t-1/2)
triangle(t)=(2*(t+1/2)-4*t*u(t))*pulse(t)


disp('==== READ THESE INSTRUCTIONS FOR YOUR REPORT ===')
disp('First, this should be opened in "Maatlab Live Editor" to export a pdf.' ) 
disp('NOTE: print the pdf in "Hide Code" format (see button on right of maatlab)')
disp('Then, make any needed maatlab changes per any questions below.') 
disp('You must re-run the script for changes to take effect.') 
disp('Dont forget to add the names of the group members as shown below.')
disp('Answer all numbered questions below, Q1, Q2, etc.') 
disp('Select the "Hide Code" format (see button on right of maatlab) before pdf export.')
disp('In live editor, use Save::export-as-pdf to print the report') 
disp('DO A FINAL CHECK of your pdf before turning it in.') 
disp(' ======================================================.')
disp(' ')
disp('ECGR 4123 Computer Project: Fourier Series')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: Basic Functions ================================')
%create symbolic function
syms omega t x(t) y(t) g(t) h(t)
f=500
omega= f*2*pi
g(t)=exp(-500*t)
x(t)=sin(omega*t)*pulse(f/2*(t-0.005))
%x(t)=pulse(f*(t-0.005))
y(t)=cos(omega*t)
h(t)=g(t)*y(t)
hh=figure(1);
subplot(2,2,1);fplot(x,[0 0.01],'LineWidth',3);grid on;
title('x(t)');xlabel('Time (s)'); ylabel('x(t)');
subplot(2,2,2);fplot(y,[0 0.01],'LineWidth',3);grid on;
title('y(t)');xlabel('Time (s)'); ylabel('y(t)');
set(findall(hh,'Type','text'),'FontSize',16);
subplot(2,2,3);fplot(g,[0 0.01],'LineWidth',3);grid on;
title('g(t)');xlabel('Time (s)'); ylabel('g(t)');
subplot(2,2,4);fplot(h,[0 0.01],'LineWidth',3);grid on;
title('h(t)');xlabel('Time (s)'); ylabel('h(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 1.  Time functions. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the formula to f=1 kHz, and observe the new plots')
disp('     for h(t) in Figure 1, and include it in your report')
disp('     Your turned-in plots should change accordingly.')
disp(' ')
disp(' ')


disp('PART 2: Energy  ==========================')
syms omega t xsq(t) ysq(t) gsq(t) hsq(t)
gsq(t)=g(t)*g(t)
xsq(t)=x(t)*x(t)
ysq(t)=y(t)*y(t)
hsq(t)=h(t)*h(t)
hh=figure(2);
subplot(2,2,1);fplot(xsq,[0 0.01],'LineWidth',3);grid on;
title('x^2(t)');xlabel('Time (s)'); ylabel('x^2(t)');
subplot(2,2,2);fplot(ysq,[0 0.01],'LineWidth',3);grid on;
title('y^2(t)');xlabel('Time (s)'); ylabel('y^2(t)');
set(findall(hh,'Type','text'),'FontSize',16);
subplot(2,2,3);fplot(gsq,[0 0.01],'LineWidth',3);grid on;
title('g^2(t)');xlabel('Time (s)'); ylabel('g^2(t)');
subplot(2,2,4);fplot(hsq,[0 0.01],'LineWidth',3);grid on;
title('h^2(t)');xlabel('Time (s)'); ylabel('h^2(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 2.  Time functions. ')
disp(' ')
disp(' ')


disp('Recall that the area under the square equals the energy.')
disp('symbolic integral')
intgsq=double(vpa(int(gsq(t),t,0,0.01)))
inthsq=sym(0)



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. What is the energy of g(t)? ')
disp(' ')
disp(' Answer: ???? ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. Edit the line inthsq=sym(0) to compute the energy of h(t).')
disp('      What is the energy of h(t)?')
disp(' ')
disp(' Answer: ????')
disp(' ')
disp(' ')




disp('PART 3: Vector inner product ==========================')
disp('Consider the problem of approximating a vector g by another vector f')
disp('  such that g is approximated by cf, where c is a constant. ')
disp('  From class, class the min. error solution is c= <g,f>/<f,f>,')
disp('  and the error is error = g - cf')
g=[ 1; 2; 2]
f=[ 1; 3; 1]
gDOTfOVERfDOTf=dot(g,f)/dot(f,f)     % error Vector
c=gDOTfOVERfDOTf
gMINUScf=g-c*f      % error Vector
error=gMINUScf 
magnitudeOFerror=norm(error) % =|errorVector|
errorDOTerror=dot(error,error)

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. Change the above to compute c and <e,e> for the three ')
disp('    dimensional column vectors g=[ 1; 2; 3] and f=[ 1; 1; 1]')
disp(' ')
disp(' Answer: c=???      <e,e>= ???  ')
disp(' ')
disp(' ')

disp('PART 4: Exponential Fourier Series ============================')
disp('Let us just compute the first 11 coefficients of the exponential.')
disp('  Fourier seies Dn from n=-5 to n=5 for function g1(t) and g2(t).')
disp('  Both functions are periodic, with 3-second periods.     ')
disp('  Only one period is shown.   ')


clear t g1(t) g2(t);syms  t g1(t) g2(t)
%g1(t) and g2(t) are the fundamental periods of a periodic function
g1(t)=pulse((t-1/2)/1)*t
g2(t)=pulse((t-1.5)/1.5)*triangle((t-1)/2)

hh=figure(3);
subplot(2,1,1);fplot(g1,[-5 5],'LineWidth',3);grid on;
title('g(t)');xlabel('Time (s)'); ylabel('g(t)');
subplot(2,1,2);fplot(g2,[-5 5 ],'LineWidth',3);grid on;
title('FIX THIS');xlabel('Time (s)'); ylabel('FIX THIS');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 3.  g1(t) and g2(t). ')
disp(' ')
disp('  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ')
disp(' NOTE: YOU MUST LOOK IN THE maatlab livescript script')
disp('       file on the left side of the screen to see')
disp('        Fig. 4 and Fig. 5')
disp(' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ')

%
% ========================================
% Relace the existing g2(t) with the function in Fig. 4 below, with a period of 3 s (repeats every 3 seconds as shown below).
%
%                         FIg. 4. One period of g2(t).
% Showing the periodic repitition:
%
%
%                           FIg. 5.  Showing periodic repitition of g2(t).
% ==========================================
%

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. Change g2(t) to be ONE PERIOD from t=0 to t=3 for the new function ')
disp('     that is shown in Fig. 5 above .')
disp('     (you will need to figure out the equation for the curve).')
disp('     Your plots should change accordingly.')
disp(' ')
disp(' ')




disp('Compute the coefficients ')
Dg1=0;Dg2=0;indexD=0;
To=3
omega0=2*pi/To
numD=5 %number of coefficients =2*numD+1
for n=1:2*numD+1  % 1:6 computes D(-5_ to D5, since matlab index must start at 1
Dg1(n)=(1/To)*double(vpa(int(g1(t)*exp(-i*(n-numD-1)*omega0*t),t,0,3)));
Dg2(n)=(1/To)*double(vpa(int(g2(t)*exp(-i*(n-numD-1)*omega0*t),t,0,3)));
indexD(n)=n-numD-1; %records the appropriate index of Dn
end
absDg1=abs(Dg1)
absDg2=abs(Dg2)

hh=figure(6);
subplot(2,1,1);stem(indexD,absDg1,'LineWidth',3);grid on;
title('Dn for g1(t)');xlabel('n'); ylabel('|Dn|');
subplot(2,1,2);stem(indexD,absDg2,'LineWidth',3);grid on;
title('Dn for g2(t)');xlabel('n'); ylabel('|Dn|');
set(findall(hh,'Type','text'),'FontSize',16);
set(findall(hh,'Type','text'),'FontSize',16);
disp('  Fig. 6. Magnitude of Fourier coefficients |Dn| for g1(t) and g2(t). ')
disp(' ')
disp(' ')




disp('PART 5: Reconstruction of g1(t) and g2(t) =========================')
disp('Let us reconstruct g1(t) and g2(t) using the 11 coefficients of the exponential')
disp('    Fourier seies Dn from n=-5 to n=5  for function g1(t) and g2(t).')
disp('       ')



n=0;
for tt=-5:0.1:5  
n=n+1;
tim(n)=tt;
%Note: dot() function in maatlab conjugates the second argument
%    so we must conjugate the second argument to "undo" this
g1sum(n)=real(dot((exp(i.*indexD.*omega0.*tt)),conj(Dg1)));
g2sum(n)=real(dot((exp(i.*indexD.*omega0.*tt)),conj(Dg2)));
g1orig(n)=double(vpa(g1(tt)));
g2orig(n)=double(vpa(g2(tt)));
end
hh=figure(7);
subplot(2,1,1);plot(tim,g1sum,'b',tim,g1orig,'m','LineWidth',3);grid on;
title('Reconstructed g1sum(t) using 11 Dn (dotted=original)');xlabel('Time (s)'); ylabel('g1sum(t)');
subplot(2,1,2);plot(tim,g2sum,'b',tim,g2orig,'m','LineWidth',3);grid on;
title('Reconstructed g2sum(t) using 11 Dn (dotted=original)');xlabel('Time (s)'); ylabel('g2sum(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('     Fig. 7. Reconstruction of g1(t) and g2(t) using 11 Dn. ')
disp(' ')
disp(' ')


disp('========================================================')
disp('Reconstruction of g1(t) and g2(t) (using ONLY 7 Dn)')
disp(' ')
disp('Let us re-compute using ONLY 7 coefficients Dn of the exponential')
disp('    Fourier seies Dn from n=-3 to n=3 for function g1(t) and g2(t),')
 disp('       ')
 

disp('Compute the new coefficients ')
Dg1=0;Dg2=0;indexD=0;
To=3
omega0=2*pi/To
numD=2 %number of coefficients =2*numD+1
for n=1:2*numD+1  % 1:6 computes D(-5_ to D5, since matlab index must start at 1
Dg1(n)=(1/To)*double(vpa(int(g1(t)*exp(-i*(n-numD-1)*omega0*t),t,0,3)));
Dg2(n)=(1/To)*double(vpa(int(g2(t)*exp(-i*(n-numD-1)*omega0*t),t,0,3)));
indexD(n)=n-numD-1; %records the appropriate index of Dn
end
absDg1=abs(Dg1)
absDg2=abs(Dg2)

n=0;
for tt=-5:0.1:5  
n=n+1;
tim(n)=tt;
%Note: dot() function in maatlab conjugates the second argument
%    so we must conjugate the second argument to "undo" this
g1sum(n)=real(dot((exp(i.*indexD.*omega0.*tt)),conj(Dg1)));
g2sum(n)=real(dot((exp(i.*indexD.*omega0.*tt)),conj(Dg2)));
g1orig(n)=double(vpa(g1(tt)));
g2orig(n)=double(vpa(g2(tt)));
end
hh=figure(8);
subplot(2,1,1);plot(tim,g1sum,'b',tim,g1orig,'m','LineWidth',3);grid on;
title('Reconstructed g1sum(t) using 7 Dn (dotted=original)');xlabel('Time (s)'); ylabel('g1sum(t)');
subplot(2,1,2);plot(tim,g2sum,'b',tim,g2orig,'m','LineWidth',3);grid on;
title('Reconstructed g2sum(t) using 7 Dn (dotted=original)');xlabel('Time (s)'); ylabel('g2sum(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('     Fig. 8. Reconstruction of g1(t) and g2(t) using 7 Dn. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Change the matlab code so that Fig. 8 correctly computes')
disp('   the reconstruction for 7 coefficients, Dn from n=-3 to n=3.')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Compare the area under one cycle of g2(t) divided by To to the ')
disp('     value of coeeficient D0.')
disp(' ')
disp('Answer: for g2(t),    area/To=???      D0=??? ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. For the infinitely long periodic function g2(t), what is the dc ')
disp('     component of the signal, and what is the period To of g2(t)?')
disp(' ')
disp('Answer: dc component of g2(t)=?    period of g2(t), To=????')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. For the infinitely long periodic version of function g1(t)   ')
disp('    with period To=3, what is the dc component of the signal?')
disp(' ')
disp('Answer: dc component of g1(t)=? ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. Calculate the power of g2(t) using the time-domain integral,')
disp('     and compare this to the sum of the magnitude squared of the ')
disp('     the 7 fourier coefficients, sum(|Dn|^2) from n=-3 to n=3. ')
disp(' ')
disp('Answer: power of g2(t)=???   sum(|Dn|^2) from n=-3 to n=3 =???  ')
disp(' ')
disp(' ')
