close all; clear all;
N=1024;   sigma = linspace(-pi,pi,N);     y= sin(sigma);

h=sum(y)/N;
    y = y-h;

    cn = 2/N *fft(y);        % FFT to retriece the complex Fourier coeff
    cn= cn(1:N/2);          % taking half of the coeff because FFT is symmetric

    bk = real(cn);            ak = imag(cn); % fourier coefficient

%%% calculation of the initial function approximate by FFT
  for i=1:N
        yf(i) = 0; % initialisation
        for j=1:length(ak)
            jj=j-1;
            yf(i) = (yf(i) +( ak(j)*sin(jj*sigma(i) )+ bk(j)*cos(jj*sigma(i)) ));
        end
end

%%% calculation of the first and second derivative in therms of sigma

for j = 1:N
        Yo(j) = 0;
        Yoo(j) = 0;

    for l=2:length(ak)
        k=l-1;
        kk=k;
         Yo(j) = (Yo(j) + kk*( ak(l) * cos(k*sigma(j)) - bk(l) * sin(k*sigma(j)) ) );
         Yoo(j) = (Yoo(j) - kk^2 *( ak(l) * sin(k*sigma(j)) + bk(l) * cos(k*sigma(j)) ) );
    end
end


plot(sigma,yf);   hold on
plot(sigma,Yo)
xlabel('sigma');  ylabel('y')
hold off