issue with calculation of derivative using FFT

Here is my issue: I am trying an algorithm using FFT to fit and get a function's approximative in order and calculate it's derivative in an more easy way ( for any data), but I get some discontinuity. at first I am using the initial fonction:

but as you can see on my figure, some noise appears near $\sigma = 0$ for the first derivative and it became worse for the second one.

Does anyone have a solution to fix that issue? I am kind of lost now

thanks in advance

Here is my code

    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 

Answer 1

[David Goodmanson]

Hi Yilmaz,

The problem is that your sine wave is not really periodic. With linspace, the first point is sigma = -pi and the last point is sigma = pi. These two are the same as far as sine is concerned, so you are repeating the same point twice in the array. For correct periodicity, you need the last point in the sigma array to be one point short of equaling pi. With the array

sigma = 2*pi*((0:N-1)/N) - pi;

the problems go away.

[Yilmaz Ayaz]

Hi David, thank you I just tried that and the problem is solved as you said. but now there is a new issue which is that I get 2 period instead of 1.

the code is as follow now

    N=512
    sigma = 2*pi*((0:N-1)/N) - pi;

    cn = 2/N *fft(y);      % FFT to retriece the complex Fourier coeff  end

    bk = real(cn(2:end));            ak = -1*imag(cn(2:end)); % fourier coefficient for y
        a0=cn(1);

         %%%calculation of the approximated function y and phi with fourrier
         %%%coefficient
        for i=1:length(sigma)
            yf(i) = a0; 
            for j=2:length(ak)
                jj=j-1;
                       yf(i) = (yf(i) +( -ak(j)*sin(jj*sigma(i) )+ bk(j)*cos(jj*sigma(i)) ));          % approxiamated y
            end 
        end

% derivative  of approximated functions using Fourrier coeff
 for j = 1:length(sigma)

        Yo(j) = 0;          Yoo(j) = 0; 

        for l=2:length(ak)
            k=l-1;
            kk=k;
             Yo(j) = 1/(1)*(Yo(j) +kk*( -ak(l) * cos(k*sigma(j)) - bk(l) * sin(k*sigma(j)) ) );              % 1st derivative of y
             Yoo(j) = -1/(N*N)*(Yoo(j) - kk^2 *( ak(l) * sin(k*sigma(j)) + bk(l) * cos(k*sigma(j)) ) );     % 2nd derivative of y         
        end
   end