File:Savitzky-golay pic gaussien bruite.svg

// **********
// Constants and initialisation
// **********

clear;
clf;
chdir("mypath/")

// smoothing parameters :
width = 9; // width of the sliding window (number of pts)
poldeg = 3; // degree of the polynomial

// **********
// Functions
// **********

// Convolution coefficients

function [a]=convolcoefs(m, d)
    // m : width of the window (number of pts)
    // d : degree of the polynomial
    l = (m-1)/2; // half-width of the window
    z = (-l:l)'; // standardised abscissa
    J = ones(m,d+1);
    for i = 2:d+1
        J(:,i) = z.^(i-1); // jacobian matrix
    end
    A = (J'*J)^(-1)*J';
    a = A(1:3,:);
endfunction

// smoothing, determination of the first and second derivatives

function [y, yprime, ysecond] = savitzkygolay(X, Y, width, deg)
    // X, Y: collection of data
    // width: width of the window (number of pts)
    // deg: degree of the polynomial
    n = size(X, "*");
    l = floor(width/2);
    step = (X($)-X(1))/(n - 1);
    y = Y;
    yprime=zeros(Y);
    ysecond=yprime;
    a = convolcoefs(width, deg);
    a(2,:) = 1/step*a(2,:);
    a(3,:) = 2/step^2*a(2,:);
    for i = 1:width
        Ymat(i, :) = Y(i: n-width+i)';
    end
    solution = a*Ymat;
    y = solution(1, :)';
    yprime = solution(2, :)';
    ysecond = solution(3, :)';
endfunction

// **********
// Main program
// **********

// data reading

data = read("noisy_gaussian_peak.txt", -1, 2)
Xinit = data(:,1);
Yinit = data(:,2);

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

// Data processing

[Ysmooth, Yprime, Ysecond] = savitzkygolay(Xinit, Yinit, width, poldeg);

// Display

offset = floor(width/2);
nbpts = size(Xinit, "*");
X1 = Xinit((offset+1):(nbpts-offset)); // removal of non-smotthed points

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

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

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

// **********
// Constants and initialisation
// **********

clear;
clf;
chdir("mypath\")

// smoothing parameters
width = 9; // width of the sliding window (number of pts)

// **********
// Functions
// **********

// 3rd degree polynomial

function [y]=poldegthree(A, x)
    // Horner method
    y = ((A(1).*x + A(2)).*x + A(3)).*x + A(4);
endfunction

// regression with the 3rd degree polynomial

function [A]=regression(X, Y)
    // X, Y: column vectors of "width" values ;
    // determines the polynomial of degree 3
    // a*x^2 + b*x^2 + c*x + d
    // by regression on (X, Y)
    XX = [X.^3; X.^2; X];
    [a, b, sigma] = reglin(XX, Y);
    A = [a, b];
endfunction

// smoothing, determination of the first and second derivatives

function [y, yprime, ysecond] = savitzkygolay(X, Y, larg)
    // X, Y: collection of data
    n = size(X, "*");
    l = floor(larg/2); // halfwidth
    y=Y;
    yprime=zeros(Y);
    ysecond=yprime;
    for i=(l+1):(n-l)
        intervX = X((i-l):(i+l),1);
        intervY = Y((i-l):(i+l),1);
        Aopt = regression(intervX', intervY');
        x = X(i);
        y(i) = poldegthree(Aopt,x);
//        Yfoo=poldegthree(Aopt,intervX);
//        subplot(3,1,1) ; plot(intervX, Yfoo, "r")
        yprime(i) = (3*Aopt(1)*x + 2*Aopt(2))*x + Aopt(3); // Horner
        ysecond(i) = 6*Aopt(1)*x + 2*Aopt(2);
    end
endfunction

// **********
// Main program
// **********

// data reading

data = read("noisy_gaussian_peak.txt", -1, 2)
Xinit = data(:,1);
Yinit = data(:,2);

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

// Data processing

[Ysmooth, Yprime, Ysecond] = savitzkygolay(Xinit, Yinit, width);

// Display

offset = floor(width/2);
nbpts = size(Xinit, "*");
vector1 = (offset+1):(nbpts-offset); // suppresion des points non-lissés
X1 = Xinit(vector1);
Y0 = Yliss(vector1)
Y1 = Yprime(vector1);
Y2 = Ysecond(vector1);

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