% This function illustrates reconstruction of a periodic function via 
% the fft.
function main_fourier_series

%DRIVER1
DRIVER2

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: This function uses the "incorrect" interpolation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER1

N  = 15;
tt = linspace(0,2*pi*(1 - 1/N),N)';
uu = ufun(tt);
cc = fft(uu);

ttfine = linspace(0,2*pi,1000);
uufine = ufun(ttfine);
vvfine = zeros(size(uufine));
for k = 0:(N-1)
  vvfine = vvfine + (1/N)*cc(k+1)*exp(1i*k*ttfine);
end

figure(1)
subplot(1,2,1)
plot(tt,uu,'r.-','MarkerSize',15)
title('The function samples')

subplot(1,2,2)
plot(0:(N-1),real(cc),'r.',...
     0:(N-1),imag(cc),'b.',...
     'MarkerSize',15)
legend('real(fft(uu))','imag(fft(uu))')
title('DFT coefficients')

figure(2)
plot(ttfine,uufine,'k',...
     ttfine,real(vvfine),'b',...
     tt,uu,'r.',...
     'LineWidth',2,'MarkerSize',20)
legend('exact function','interpolated function','sample points')
xlabel('t')
keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER2

N  = 51;
tt = linspace(0,2*pi*(1 - 1/N),N)';
uu = ufun(tt);
cc = fft(uu);

ttfine = linspace(0,2*pi,1000);
uufine = ufun(ttfine);
vvfine = (1/N)*cc(1)*ones(size(ttfine));
for k = 1:floor(N/2)
  vvfine = vvfine + (1/N)*cc(k+1)*exp(1i*k*ttfine) ...
                  + (1/N)*cc(N+1-k)*exp(-1i*k*ttfine);
end

figure(1)
subplot(1,2,1)
plot(tt,uu,'r.-','MarkerSize',15)
title('The function samples')

subplot(1,2,2)
ind = [(ceil(N/2)+1):N,1:ceil(N/2)];
kk  = (-floor(N/2)):floor(N/2);
plot(kk,real(cc(ind)),'r.',...
     kk,imag(cc(ind)),'b.',...
     'MarkerSize',15)
legend('real(fft(uu))','imag(fft(uu))')
title('DFT coefficients')

figure(2)
plot(ttfine,uufine,'k',...
     ttfine,real(vvfine),'b',...
     tt,uu,'r.',...
     'LineWidth',2,'MarkerSize',20)
legend('exact function','interpolated function','sample points')
xlabel('t')
    
keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uu = ufun(tt)

%uu = cos(3*tt).*(1 - sin(2*tt).^3);
%uu = 1-abs(1-tt/pi);
%uu = 1 + cos(4*tt).*sin(tt);
uu = exp(-10*(tt-pi).^2);
%uu = exp(-10*(tt-pi).^2).*sin(14*tt);
%uu = bump(tt-pi);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uu = bump(tt)

uu      = zeros(size(tt));
ind     = find(abs(tt) < 1);
uu(ind) = exp(-1./(1-tt(ind).*tt(ind)));

return

