// **********
// 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);