// **********
// Initialisation
// **********

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

// **********
// Constants
// **********

// Initial parameter of the model
// obtained by Savitzky-Golay
// x position estimated by the minima of the second derivative
// height H by the height of the curve at these positions
// full-width at half maximum FWHM by the length of the chord
// that links the inflection points (extrema of the derivative)
// sigma = FWHM/2,35
Ainit = [1.96 ; 1.26 ; 0.29
6.12 ; 0.64 ; 0.20
7.06 ; 1.27 ; 0.26];

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

// Gaussian function

function [y] = gauss(A, x)
    // A(1): peak position
    // A(2): peak height
    // A(3): peak "width"
    y = A(2)*exp(-(x - A(1)).^2/A(3));
endfunction

function [y] = profil(A, x)
    // A(1): 1st peak position
    // A(2): 1st peak height
    // A(3): 1st peak "width"
    // A(4)-A(6): same, 2nd pic
    // A(7)-A(9): same, 3rd pic
    G1 = A(1:3);
    G2 = A(4:6);
    G3 = A(7:9);
    y = gauss(G1, x) + gauss(G2, x) + gauss(G3, x);
endfunction

// Residual

function [e] = res(A, x, y)
    e = profil(A, x) - y;

endfunction

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

// Data reading

data = read('three_gaussian_peaks.txt', -1, 2)
X = data(:,1);
Y = data(:,2);

[S, Aopt] = leastsq(list(res, X, Y), Ainit)
Yopt = profil(Aopt, X);
YGopt1 = gauss(Aopt(1:3),X);
YGopt2 = gauss(Aopt(4:6),X);
YGopt3 = gauss(Aopt(7:9),X);

// Display

plot2d(X, Y, style=0)
//plot(X, Yopt, "-r")
plot(X, YGopt1, "-c")
plot(X, YGopt2, "-g")
plot(X, YGopt3, "-b")

height = max(Y);
xmin = min(X)

xstring(xmin, height,...
 'Pic 1 : X1 = '+string(Aopt(1))+' pour 2'...
+' ; Al = '+string(Aopt(2))+' pour 1,26'...
+' ; H1 = '+string(Aopt(3))+' pour 0,2.')
xstring(xmin, 0.95*height,...
 'Pic 2 : X2 = '+string(Aopt(4))+' pour 6,2'...
+' ; A2 = '+string(Aopt(5))+' pour 0,63'...
+' ; H2 = '+string(Aopt(6))+' pour 0,2.')
xstring(xmin, 0.9*height,...
 'Pic 3 : X3 = '+string(Aopt(7))+' pour 7'...
+' ; A3 = '+string(Aopt(8))+' pour 1,26'...
+' ; H3 = '+string(Aopt(9))+' pour 0,2.')
xstring(xmin, 0.85*height,...
 'variance : S = '+string(S))