// **********
// Constantes et initialisation
// **********

clear;
clf;
chdir("monchemin\")

// paramètres du lissage :
largeur = 9; // largeur de la fenêtre glissante (nb de pts)

// **********
// Fonctions
// **********

// polynôme de degré 3

function [y]=poldegtrois(A, x)
    // méthode de Horner
    y = ((A(1).*x + A(2)).*x + A(3)).*x + A(4);
endfunction

// régression avec le polynôme de degré 3

function [A]=regression(X, Y)
    // X et Y : vecteurs colonne de 9 valeurs ;
    // détermine le polynôme de degré 3
    // a*x^2 + b*x^2 + c*x + d
    // par régression sur (X, Y)
    XX = [X.^3; X.^2; X];
    [a, b, sigma] = reglin(XX, Y);
    A = [a, b];
endfunction

// lissage, détermination de la dérivée et de la dérivée seconde

function [y, yprime, yseconde] = savitzkygolay(X, Y, larg)
    // X, Y : nuage de points
    // larg : largeur de fenêtre
    n = size(X, "*");
    decalage = floor(larg/2);
    y = Y;
    yprime = zeros(Y);
    yseconde = yprime;
    for i=(decalage+1):(n-decalage)
        intervX = X((i-decalage):(i+decalage),1);
        intervY = Y((i-decalage):(i+decalage),1);
        Aopt = regression(intervX', intervY');
        x = X(i);
        y(i) = poldegtrois(Aopt,x);
// running plot and picture saving to be introduced
        yprime(i) = (3*Aopt(1)*x + 2*Aopt(2))*x + Aopt(3); // Horner
        yseconde(i) = 6*Aopt(1)*x + 2*Aopt(2);
    end
endfunction

// **********
// Programme principal
// **********

// lecture des données

donnees = read("mes_donnees.txt", -1, 2)
Xinit = donnees(:,1);
Yinit = donnees(:,2);

// Traitement des données

[Yliss, Yprime, Yseconde] = savitzkygolay(Xinit, Yinit, largeur);