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

clear;
clf;
chdir('monchemin/')

// paramètres du lissage :
largeur = 41; // 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_gauss_dissym_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,'*');
X1=Xinit((decal+1):(nbpts-decal));
X2=Xinit((2*decal+1):(nbpts-2*decal));
Y1=Yprime((decal+1):(nbpts-decal));
Y2=Yseconde((2*decal+1):(nbpts-2*decal));

[foo, pinfg] = max(Yprime); // position du point d'inflexion à gauche
[foo, pinfd] = min(Yprime); // " droite
[foo, pos] = min(Yseconde); // position du sommet
ymax = Yliss(pos); // valeur du maximum
print(%io(2), Xinit(pinfg))
print(%io(2), Xinit(pinfd))
print(%io(2), Xinit(pos))
print(%io(2), ymax)

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

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

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