//x5.sci
clear;

T = 1;            // periode
N = 50;           // nombre d'echantillons dans T
KK = 10;          // nombre de periodes tracees
J = 50;
sigma = 2;        // ecart-type du bruit

K1 = KK + 1;
tChoix = [ "+ 1";  "+ " + string(KK); ..
           "+ " + string(10*KK); "+ " + string(100*KK) ];
nChoix = size(tChoix, "r");
tempos = [80000,0,20000,0];
tempo  = tempos(1);


function init_affichage()
    set("figure_style", "old"); 
    driver("X11");
    xset("default");
    xset("font", 0, 3);
    xset("pixmap", 1);
    xpause(50000);
    xbasc();
    xset("wdim", 940, 600);
    xset("wpos", 100,  0);
    xselect();
endfunction


function fin_affichage()
    xset("default");
    xset("pixmap", 0);
    driver("Rec");
endfunction


function z = normal(y)
    z = sigma*rand(y, "normal");
endfunction
 
 
function z = uniform(y)
    z = sigma*sqrt(12)*(rand(y, "uniform") - 0.5);
endfunction
 
 
function z = cauchy(y)
    z = tan( %pi*(rand(y, "uniform") - 0.5) )/4;
endfunction


function courbe(signal, cadre)
    plot2d(x, signal(:), 1, "011", " ", cadre, ticks);
    if pur 
        xset( "thickness", 2);
        plot2d( x, sp(:), 5, "000" );
        xset( "thickness", 1);
    end 
endfunction
 

function dessin(fin)
    xset("wwpc");
    debut = fin - KK*T;
    cadre = [debut, -ymax, fin, ymax];
    xsetech([0.0, 0.0, 1.0, 0.5]);
    courbe( sb, cadre);
    xsetech([0.0, 0.5, 1.0, 0.5]);
    courbe( sf, cadre);
    xtitle(sprintf("Signal bruité [en haut] et débruité [en bas] de " + ..
                   "%.0f a %.0f/%.0f", ..
                    debut, fin, K*T) );
    xset("wshow");
    if ( tempo > 0 )
        xpause(tempo);
    end
endfunction


init_affichage;

// function filtre()



x = linspace(0, T, N + 1)';
x($) = [];    
y  = sin(2*%pi*x/T);
ymax  = max(abs( y)) + 2*sigma;

cadre = [0, -ymax, KK*T, ymax];
ticks = [4, KK, 4, 10] ;
aRect = [0.05, 0.05, 0.10, 0.10];    
alea  = normal;

while %T;
    pur = %F;
    sp = y * ones(1,K1);
    sb = sp + alea(sp);
    sf = cumsum(sb, "c") ./ (  ones(N,1) * [1:K1] );
    
    x = linspace(0, K1*T, K1*N + 1)';
    x($) = [];
          
    xset("wwpc");
    xsetech(wrect = [0.0, 0.0, 1.0, 0.5], arect = aRect);
    courbe(sb, cadre);
    xsetech(wrect = [0.0, 0.5, 1.0, 0.5], arect = aRect);
    xtitle("Signal de période 1 + bruit (sur 10 périodes)");
    xset("color", 5);
    titlepage("  Cliquer pour débruiter ! ");
    xset("color", 1);
    xset("wshow");  
    xselect();
    xclick();
    
    K = KK; 
    k = K;  
    plusK = 1;       
    while %T
                    
        if  ( k == K )
            dessin(k*T);
            choix = x_choose( [tChoix; "Signal"], ..
                              "Periodes", "Stop");
            if ( choix == 0)
                break;
            elseif ( choix <=  nChoix )
                plusK = evstr( tChoix(choix) );
                K = K + plusK;
                tempo = tempos(choix);
            else
                pur = ~pur; 
                continue;
            end    
        end
        
        if ( plusK < 100 )
            for j = 1:J
                dessin((k + j/J)*T);
            end
        elseif ( plusK < 1000  )
            dessin(k*T);
        end
        
        k = k + 1;
        x = x + T;
        sb = [sb(:,2:$), y + alea(y)];
        sf = [sf(:,2:$), sf(:,$) + (sb(:,$) - sf(:,$))/k];   
    end

    select x_choose( ["de Gauss"; "uniforme"; "de Cauchy"], ..
                     ["Type de bruit"], "Quitter");     
    case 0 
        break;
    case 1
        alea = normal;
    case 2
        alea = uniform;
    case 3
        alea = cauchy;
    end 
        
end    
// endfunction

fin_affichage;