plot(tt,V,'b-',tt,V,'ro')
title('Measured data')
%
% Number of data points
%
N=size(tt,2)
%
% Period
%
t1=tt(N)
t0=tt(1)
p=t1-t0
tttt=[0:p/300:p];
ttt=tttt+t0*ones(size(tttt));
%
% time increments Delta-t's
%
dt=tt(2:N)-tt(1:N-1);
%
% f_0(t) = 1, a_0 zero-th order Fourier coeff.
% using trapezoidal rule
%
a0=sum(((V(2:N)+V(1:N-1))/2).*dt)/p
hold on
c0=ones(size(ttt));
app=a0*c0;
plot(tt,V,'k-',tt,V,'bo',ttt,app,'r-',tt+p,V,'k-',tt+p,V,'bo',ttt+p,app,'r-')
title('Zeroth order Fourier approximation')
pause
hold off
%
% maximal order of Fourier approximation
%
M=10
aa=zeros(1,M);
bb=zeros(1,M);
for k=1:M
    pause(1.2)
    %
    % Evaluate the shifted sines and cosines at sample times
    %
    ck=cos(2*pi*k*(tt-t0*ones(size(tt)))/p);
    sk=sin(2*pi*k*(tt-t0*ones(size(tt)))/p);
    %
    % Multiply the data against the weight functions cos(kt) and sin(kt)
    %
    Vck=ck.*V;
    Vsk=sk.*V;
    %
    % Trapezoidal rule
    %    
    ak=2*sum(((Vck(2:N)+Vck(1:N-1))/2).*dt)/p
    bk=2*sum(((Vsk(2:N)+Vsk(1:N-1))/2).*dt)/p
    %
    % Store newly computed coefficients for later access in arrays
    %
    aa(k)=ak;
    bb(k)=bk;
    app=app+ak*cos(tttt*2*pi/p*k)+bk*sin(tttt*2*pi/p*k);
    plot(tt,V,'k-',tt,V,'bo',ttt,app,'r-',tt+p,V,'k-',tt+p,V,'bo',ttt+p,app,'r-')
    mytitle=strcat('Fourier approximation of order ',num2str(k));
    title(mytitle)
    pause
    end
 hold off
 plot([1:M],aa,'b-',[1:M],bb,'g-',[1:M],sqrt(aa.*aa+bb.*bb),'r-',[1,M],[0,0],'k-')
 title('Fourier coefficients = spectrum of the signal')