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

clear;
clf;
chdir("monchemin/")

// paramètres du lissage :
largeur = 9; // largeur de la fenêtre glissante (nb de pts)
degre = 3; // degré du polynôme

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

// Coefficients de convolution

function [a]=coefconvol(m, d)
    // m = largeur de la fenêtre
    // d : degré du polynôme
    l = (m-1)/2; // demie-largeur de la fenêtre
    z = (-l:l)'; // intervalle centré réduit
    J = ones(m,d+1);
    for i = 2:d+1
        J(:,i) = z.^(i-1); // matrice jacobienne
    end
    A = (J'*J)^(-1)*J';
    a = A(1:3,:);
endfunction

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

function [y, yprime, yseconde] = savitzkygolay(X, Y, larg, deg)
    // X, Y : nuage de points
    // larg : largeur de la fenêtre
    // deg : degré du polynôme
    n = size(X, "*");
    l = floor(larg/2);
    pas = (X($) - X(1))/(n-1);
//    y=Y;
//    Ymat = ones(l, n - larg);
//    yprime = zeros(Y);
//    yseconde = yprime;
    a = coefconvol(larg, deg);
    a(2, :) = 1/pas*a(2, :);
    a(3, :) = 2/pas^2*a(3, :);
    for i = 1:larg
        Ymat(i, :) = Y(i: n-larg+i)';
    end
    solution = a*Ymat;
    y = solution(1, :)';
    yprime = solution(2, :)';
    yseconde = solution(3, :)';
endfunction

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

// lecture des données

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

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

// Traitement des données

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

// affichage

decal = floor(largeur/2);
nbpts = size(Xinit, "*");
X1 = Xinit((decal+1):(nbpts-decal)); // on enlève les points non lissés

clf;

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

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

subplot(3,1,3)
plot(X1, Yseconde, "b")

// **********
// 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);
//        Yfoo=poldegtrois(Aopt,intervX);
//        subplot(3,1,1);plot(intervX, Yfoo, "r")
        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("pic_gaussien_bruite.txt", -1, 2)
Xinit = donnees(:,1);
Yinit = donnees(:,2);

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

// Traitement des données

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

// affichage

decal = floor(largeur/2);
nbpts = size(Xinit, "*");
vecteur1 = (decal+1):(nbpts-decal); // suppression des points non-lissés
X1 = Xinit(vecteur1);
Y0 = Yliss(vecteur1);
Y1 = Yprime(vecteur1);
Y2 = Yseconde(vecteur1);

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

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

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