hold off
echo on
axis([0 100 0 4])
B=A(1500:1600,2);
N=size(B,1)
tt=(1:N)';
axis([0 N 0 3])
plot(tt,B)
title('Measured data')
pause
c0=ones(size(tt));
a0=sum(c0.*B)/N
hold on
axis([0 100 0 3])
app=a0*c0;
plot(app,'r-')
title('Zeroth order Fourier approximation')
hold off
M=50
aa=zeros(1,M);
bb=zeros(1,M);
for k=1:M
    pause
    ck=cos(tt*2*pi*k/N);
    sk=sin(tt*2*pi*k/N);
    ak=sum(ck.*B)*2/N
    bk=sum(sk.*B)*2/N
    aa(k)=ak;
    bb(k)=bk;
    app=app+ak*ck+bk*sk;
    plot(tt,B,'b-',tt,app,'r-')
    title('Fourier approximation of order')
    text(100,1,'Fourier approximation of order','FontSize',18)
    axis([0 N 1 3])
    end
 pause
 hold off
 plot([1:50],aa,'b-',[1:50],bb,'r-')
title('Fourier coefficients = spectrum of the signal')
pause
plot(tt,B,'b-',tt,a0*c0+aa(38)*cos(tt*2*pi*38/N),'r-')
title('only mean plus one frequency, one cosine')    
pause
plot(tt,B-a0*c0-aa(38)*cos(tt*2*pi*38/N),'r-')
title('difference: signal minus mean minus one frequency, one cosine')    
pause
plot(tt,B-a0*c0-aa(38)*cos(tt*2*pi*38/N),'r-')