File:Savitzky-golay pic gaussien bruite norme l1.svg
Summary
| Description |
Français : Lissage d'un pic gaussien bruité, avec une verison modifiée de l'algorithme de Savitzky-Golay : le polynôme de degré 3 est ajusté par régression en utilisant la norme ℓ1 (somme des valeurs absolues des résidus) au lieu des moindres carrés.
Cette méthode est réputée moins sensible aux points aberrants. Créé avec Scilab, modiufié avec Inkscape.English: Smoothing of a noisy gaussian peak using a modified version of the Savitzky-Golay algorithm: the 3rd degree polynomial is determined by regression using the ℓ1 norm (sum of the absolute values of the residuals) instead of the least squares.
This method is supposed to be less sensitive to outliers. Created with Scialb, modified with Inkscape. |
| Date | |
| Source | Own work. ame data as File:Savitzky-golay pic gaussien bruite.svg |
| Author | Cdang |
Scilab source
 |
This media was created with Scilab, a free open-source software. Here is a listing of the Scilab source used to create this file. |
For the generation of the data, see File:Savitzky-golay pic gaussien bruite.svg.
For the use of the ℓ1 norm, we use the standard non-linear least square regression dunction, leastsq(), but the residual is defined as the square root of the difference.
// **********
// 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")
Licensing
I, the copyright holder of this work, hereby publish it under the following license:
This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license.
- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
