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

clear;
clf;
chdir('C:\Documents and Settings\christophe.dang\My Documents\documentation\algo\regression\')

// 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 [R] = residual(A, X, Y)
    // calcul du résidus selon la norme l1
    // sqrt : pour pouvoir utiliser leastsq(), car
    // (sqrt())^2 = abs()
    R = sqrt(abs(Y - poldegtrois(A, X)))
endfunction

function [Aopt]=regression(Ainit, 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)
    [S, Aopt] = leastsq(list(residual, X, Y), Ainit);
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
    Ainit = [0, 0, 0, 0];
    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(Ainit, 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,'*');
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));

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")