function [x,t]=ctfs_synthesis(T0,ak,t)
%%% Function to implement CT Fourier series synthesis equation
%%%
%%%   [x,t]=ctfs_synthesis(T0,ak,t)
%%%    
%%%   Variables:
%%%        T0 = fundamental period
%%%        ak = Fourier series coefficients
%%%         t = vector of times [OPTIONAL; if no times are specified
%%%                              function computes values for 2 periods]
%%%         x = x(t), samples of the periodic time signal
%%%

Kmax=(length(ak)-1)/2;

if nargin<2
  error('CTFS_SYNTHESIS:  not enough input arguments')
elseif nargin<3 | isempty(t)
  if Kmax>0
    maxfreq=Kmax/T0;
  else
    maxfreq=1/T0;
  end
  t=-T0:1/(8*maxfreq):T0; % if t not specified, sample at four times Nyquist
end

w0=2*pi/T0;				% Fundamental frequency
k=-Kmax:Kmax;				% Vector of k's

Ek=exp(j*t'*k*w0);			% Matrix of complex exponentials

ak=ak(:);				% Make sure ak is column vector

x=Ek*ak;				% Compute x as matrix multiply

%%% Note:  100*eps is on order of 2*10^(-14).  This is needed to deal
%%% with impulse train calculations.  You'll still get imaginary values
%%% if you use lots of coefficients for the impulse train.

if max(abs(imag(x)))<100*eps		% Make sure x is real
  x=real(x);				% (eliminate small imaginary part)
end

return