clear;
chdir('mypath/')

// Parameters of the function
 
step = 0.02;
X = [0:step:10]';
speak = 0.1; // sigma peaks
varnoise = 0.01; // variance of the noise
 
// Gauss function
 
function [y]=gauss(x, sigma)
   y=1/(sqrt(2*%pi*sigma))*exp(-x^2/(2*sigma));
endfunction
 
// global function
 
function [y]=f(x, sigma, variancenoise)
    for i = 1:size(x,'*')
        y (i) = gauss(x(i) - 2, sigma) + 0.5*gauss(x(i) - 6.2, sigma)...
        + gauss(x(i)-7, sigma) + variancenoise*rand(1,'normal');
    end
endfunction
 
// Data generation

Y = f(X, speak, varnoise);

// plot(X,Y);

// Saving

write('three_gaussian_peaks.txt', [X,Y]);

// **********
// Constants and initialisation
// **********

clear;
clf;
chdir('mypath/')

// smoothing parameters :
width = 9; // width of the sliding window (number of pts)
poldeg = 3; // degree of the polynomial

// **********
// Functions
// **********

// Convolution coefficients

function [a]=convolcoefs(m, d)
    // m : width of the window
    // d : degree of the polynomial
    l = (m-1)/2; // half-width of the window
    a = ones(3,m);
    z = (-l:l)'; // standardised abscissa
    J = ones(m,d+1);
    for i = 2:d+1
        J(:,i) = z.^(i-1); // jacobian matrix
    end
    A = (J'*J)^{-1}*J';
    a = A(1:3,:);
endfunction

// smoothing, determination of the first and second derivatives

function [y, yprime, ysecond] = savitzkygolay(X, Y, wdt, deg)
    // X, Y: collection of data
    // wdt: width of the window
    // deg : degree of the polynomial
    n = size(X,'*');
    l = floor(wdt/2);
    step = (X($) - X(1))/(n - 1);
    y=Y;
    yprime = zeros(Y);
    ysecond = yprime;
    a = convolcoefs(wdt, deg);
    a(2,:) = 1/step*a(2,:);
    a(3,:) = 2/step^2*a(3,:);
    for i=(l+1):(n-l)
        y(i) = a(1,:)*Y(i-l:i+l);
        yprime(i) = a(2,:)*Y(i-l:i+l);
        ysecond(i) = a(3,:)*Y(i-l:i+l);
    end
endfunction

// **********
// Main program
// **********

// data reading

data = read('three_gaussian_peaks.txt', -1, 2)
Xinit = data(:,1);
Yinit = data(:,2);

//subplot(3,1,1)
//plot(Xdef, Ydef, "b")

// Data processing

[Ysmooth, Yprime, Ysecond] = savitzkygolay(Xinit, Yinit, width, poldeg);

// Display

offset = floor(width/2);
nbpts=size(Xinit,'*');
X1=Xinit((offset+1):(nbpts-offset));
X2=Xinit((2*offset+1):(nbpts-2*offset));
Y1=Yprime((offset+1):(nbpts-offset));
Y2=Ysecond((2*offset+1):(nbpts-2*offset));

subplot(3,1,1)
plot(Xinit, Yinit, "b")
plot(Xinit, Ysmooth, "r")

subplot(3,1,2)
plot(X1, Y1, "b")

subplot(3,1,3)
plot(X2, Y2, "b")