mode(0)
clear all;

// This code has been freely adapted from the one of Eli Billauer
// http://www.billauer.co.il/peakdet.html 
//
// Axes of progress : 
//      - determination if the first extrema is a peak or a valley / currently the assumption has been made that is a peak (use of filtering to smooth the curve and then determine if versmax is initialy %t)
//      - does vectorization possible (at least a difficult to implement) ?
//      - use of moving average curves?


function [max_vector,min_vector,status] = peak_valley_searh(vect, tol, lower_limit, upper_limit)

    status = 0; 

    // verifications
    [nl,nc] = size(vect);
    N = size(vect,"*");

    if ( argn(2)  < 2) then
        printf("Not enough input argument => the calculation stops !\n");
        status = 1;

    elseif ( (nl <> N) & (nc <> N) ) then
        printf("the input vector has not the correct dimension - it must be a row or column vector !\n");
        status = 1;

    elseif (tol <= 0.) then
        printf("Tolerance must strictly be positive i.e. greater than 0 - the calculation stops !\n");
        status = 1;

    else

        // program beginning
        // variables initialization
        // min_vector and max_vector have initially the same size than vect (avoiding dynamic allocation) ; lzeros will be removed at the last stage
        max_vector = zeros(N,2); min_vector = zeros(N,2);
        tempo_min = %inf; tempo_max = -%inf;
        position_min = %nan; position_max = %nan;
        versmax = %t; // Mean assymption here (but not necessary true: the first extremum is a peak => need to be improved

        for i = 1 : N
            current_point = vect(i);

            if (current_point > tempo_max ) then
                tempo_max = current_point;
                position_max = i;
            end
            if (current_point < tempo_min) then
                tempo_min = current_point;
                position_min = i;
            end

            if (versmax == %t) then
                if (current_point < (tempo_max - tol)) then
                    max_vector(i,1) = position_max;
                    max_vector(i,2) = tempo_max;
                    tempo_min = current_point;
                    position_min = i;
                    versmax = %f; // now the next extremum is a valley
                end
            else
                if (current_point > (tempo_min + tol)) then
                    min_vector(i,1) = position_min;
                    min_vector(i,2) = tempo_min;
                    tempo_max = current_point;
                    position_max = i;
                    versmax = %t; // // now the next extremum is a peak                  
                end
            end
            status = 0;
        end

        // zeros values are removed
        ind = find(max_vector(:,1) == 0 );
        max_vector(ind,:) = []; clear ind;
        ind = find(min_vector(:,1) == 0 );
        min_vector(ind,:) = []; clear ind;

        // all the peaks above the upper limit are removed
        if (upper_limit <> %f) then
            ind = find(max_vector(:,2) > upper_limit );
            max_vector(ind,:) = []; clear ind; 
        end
        //all the valleys bellow the upper limit are removed
        if (lower_limit <> %f) then
            ind = find(min_vector(:,2) < lower_limit ); 
            min_vector(ind,:) = []; clear ind;
        end


    end
endfunction


ni = 0; ns=20; stp = 0.001;
t=[ni:stp:ns]';
n = (ns/stp)+1;
x=0.3*sin(t) + sin(1.3*t) + 0.9*sin(4.2*t) + 0.02*rand(n,1,"normal");

//lower_limit = %f;
//upper_limit =%f;
lower_limit = -0.4;
upper_limit = 0.8;
tol = 0.5;

// function call
[max_vector,min_vector,status] = peak_valley_searh(x, tol,lower_limit,upper_limit);


scf(1);
drawlater();
xgrid(3);
f = gcf(); 
//f                                          
f.figure_size = [1000, 1000];                
f.background = color(255,255,255);           
a = gca();                                   
//a                                          
a.font_size = 2;                             
a.x_label.text = "X axis" ;                  
a.x_location="bottom";                       
a.x_label.font_angle=0;                      
a.x_label.font_size = 4;                     
a.y_label.text = "Y axis";
a.y_location="left";
a.y_label.font_angle=-90;
a.Y_label.font_size = 4;
a.title.text = "Title";                     
a.title.font_size = 5;
a.line_style = 1;	

// début des courbes
plot(t,x);
e1 = gce(); 
p1 = e1.children;                           
p1.thickness = 1;
p1.line_style = 1;
p1.foreground = 2;

n = size(min_vector,"*");
if  (n <> 0) then
    plot(t(min_vector(:,1)), min_vector(:,2));
    e2 = gce();                            
    p2 = e2.children;                      
    p2.thickness = 2;
    p2.line_mode = "off";                  // line is desactivated -> only marks are plotted
    p2.mark_mode="on";
    p2.mark_style = 14;
    p2.mark_size_unit="tabulated";
    p2.mark_size = 5;
    p2.mark_foreground = 3;
    p2.mark_background = 0;                 // 0 = transparent
end

n = size(max_vector,"*");
if  (n <> 0) then
    plot(t(max_vector(:,1)), max_vector(:,2));
    e3 = gce();                              
    p3 = e3.children;                        
    p3.thickness = 2;
    p3.line_mode = "off";
    p3.mark_mode="on";
    p3.mark_style = 10;
    p3.mark_size_unit="tabulated";
    p3.mark_size=3;
    p3.mark_foreground = 5;
    p3.mark_background = 0;                  // 0 = transparent
end

if (lower_limit <> %f) then
    lineinf = [t(1) lower_limit ; t($) lower_limit];
    plot(lineinf(:,1), lineinf(:,2));
    e4 = gce();
    p4 = e4.children;                           
    p4.line_mode = "on";
    p4.thickness = 2;
    p4.line_style = 3;                          // dash line
    p4.foreground = 6;
end


if (upper_limit <> %f) then
    linesup = [t(1) upper_limit ; t($) upper_limit];
    plot(linesup(:,1), linesup(:,2));
    e5 = gce();
    p5 = e5.children;                           
    p5.line_mode = "on";
    p5.thickness = 2;
    p5.line_style = 3;                          // dash lines
    p5.foreground = 1;
end

noms_courbes = ["Courbe" "Valleys" "Peaks" "Lower limit" "Upper limit"];

// add of the legend
legend(noms_courbes,pos=1)
drawnow();
xs2png(1,'peak_search.png');