fit curve with fft and calculate it's derivatives

274 Views Asked by At

Here is my issue: I am trying an algorithm using FFT to fit and get a function's approximative in order and calculate its derivative in an easier way (for any data), but I get some discontinuity. At first I am using the initial function: $$ y(\sigma, t_0) = \sin(\sigma) $$ with $-\pi \le \sigma \le \pi$ to try it.

After the FFT the function can be approximated by: $$ y\left(\sigma ,t\right)=\sum _{k=0}^{N/2}[a_k\left(t\right)\sin\left(k\sigma \right)+b_k\left(t\right)\cos\left(k\sigma \right)]. $$

The first and second derivatives should be :

$$ y_{\sigma }=\sum _{k=1}^{N/2}k\left[a_k\left(t\right)\cos\left(k\sigma\right)-b_k\left(t\right)\sin\left(k\sigma \right)\right] $$

$$ y_{\sigma \sigma}=-\sum _{k=1}^{N/2}k^2\left[a_k\left(t\right)\sin\left(k\sigma \right)+b_k\left(t\right)\cos\left(k\sigma\right)\right],$$

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

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

Thanks in advance.

figure1

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