% 
%  Complex Fourier Series
%
%  inputs:  Tlow and Thigh: The period goes from Tlow to Thigh
%           N = the number of terms to use in the Fourier series
%           N_Periods = the number of periods to plot
%
function Fourier_Series(Tlow,Thigh,N,N_Periods)
global T m
%
%  set the fundamental period
%
  T = Thigh-Tlow;
%
%##################################################################
%##################################################################
%
%  first numerically integrate to get the dc value. c0
%
  c0 = (1/T)*quadl(@fcn,Tlow,Thigh);
%
%  numerically integrate to get the remaining coefficients, cn
%
  for i = 1:N
      m = i;  
      c(i) = quadl(@integrand,Tlow,Thigh);
  end;
%
%  Now compute the true function
%
  t = linspace(Tlow,Thigh,10000);
  true_value = feval(@fcn,t);
%
%  replicate the function over the desired periods
%
  tp = t;
  tv = true_value;
  for kk = 2: N_Periods
      tp = tp + T;
      t = [t tp(2:end)];
      true_value = [true_value tv(2:end)];
  end;
%
% we need to look at negative and positive values for m
% so we use peoperties of Fourier Series, 
%
% The Fourier Series estimate will end up in the variable est
%
  m = 1;
  est = T*c(1)*conj(feval(@periodic,t))+T*conj(c(1))*feval(@periodic,t);
%  
  for i = 2:N
    m = i; 
    est = est + T*c(i)*conj(feval(@periodic,t));
    est = est + T*conj(c(i))*feval(@periodic,t);
  end;
  est = est+c0;
%
%##############################################################
%#############################################################
%
%  This is where you need to add your calls to various functions
%  you will be writing
%
%  at this point the following variables have been defined:
%
%  the Fourier series representation is in the variable 'est'
%  the time vector is in the variable 't'
%  the true (original) function is in variable 'true_value'
%  the Fourier series coefficients are in 'c0' and in the array 'c'
%   However, the c array goes from c(1) to c(N), c = [c1 c2 ... cN]
%  the period is in the variable 'T'
%
%  put calls to your functions below here
%





%
%  plot the true value and the Fourier Series representation
%
  plot(t,true_value,'-',t,real(est),'--',t,imag(est),'.-'); grid; 
  title(['Fourier Series, N = ' num2str(N) ]);
  legend('True','Real Estimate', 'Imaginary Estimate');
%
%  we are all done with the main program now
%
return;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  This is the function we want to determine the Fourier series representation of
%  You will modify this for each of your functions...
%  The peroid is in the variable T, which is you can use if you need to
%
function value = fcn(t);
global T m
%
%  return the function in the variable `value'
%
%  the default function
%
  value = 0.135*exp(t).*(t<2)+1*((t>=2)&(t<3))+(4-t).*((t>=3)&(t<4))+0*(t>=4);
%
  return;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  add your functions here
%




%
% Don't change anything below this line
%
%###############################################################
%###############################################################
%
%  Integrate the function times the appropriate complex exponential
%
function c = integrand(x)
  c = feval(@fcn,x).*feval(@periodic,x);
return;
%
%  Determine the proper complex exponential. In order to use this with the
%  integrand functions, the harmonic and fundamental period are global
%  parameters ...
%
function value = periodic(x)
global T m
  value = (1/T)*exp(-j*2*pi*m*x/T);
return;

