% Fourier series computation
% by B. Boulet, 2004

% signal amplitude
A=1;
% signal period
T=1;
% number of harmonics to sum
N=50;
% time vector
finaltime=T;
starttime=-T
tstep=T/100;
t=[starttime:tstep:finaltime]';

% vector calc of each complex exponential
%a0=(exp(1)-1);
a0=0;
%First compute the ak's
spectrum=zeros(1,2*N+1);
for k=-N:N
 

 %Sawtooth going back to zero:
   if k>0
      eval(['a' num2str(k) '=(1/(j*k*2*pi)+(1-(-1)^k)/(2*(k*pi)^2))*((2*pi)/(1+j*4*pi))'])
   elseif k<0
      eval(['a_' num2str(abs(k)) '=(1/(j*k*2*pi)+(1-(-1)^k)/(2*(k*pi)^2))*((2*pi)/(1+j*4*pi))'])      


end
eval(['spectrum(N+1)=a0']);

end
%then compute the harmonic components of the signal
for k=-N:N
   if k>=0
      eval(['harm' num2str(k) '=a' num2str(k) '*exp(j*2*pi*k*t/T);'])
   elseif k<0
      eval(['harm_' num2str(abs(k)) '=a_' num2str(abs(k)) '*exp(j*2*pi*k*t/T);'])
end
end
%Compute the truncated Fourier series by summing the components
truncFS=zeros(length(t),1);
for k=-N:N
   if k>=0
      eval(['truncFS=truncFS+harm' num2str(k) ';'])
   elseif k<0
      eval(['truncFS=truncFS+harm_' num2str(abs(k)) ';'])
   end
   end
%plot truncated Fourier series
figure(1)
%plot(t,truncFS);
%hold on
plot(t,truncFS,t,harm1+harm_1,'--');
%hold off


%% Line spectrum
K=[-N:N];
figure(2)
subplot(211)
stem(K,abs(spectrum),'.-')
subplot(212)
stem(K,phase(spectrum),'.-')
