function [x, e] = pulsed_signal(To,sig,par,ton,toff,tr,tf,Fs)
    // Pulsed signal generator
    //
    // Usage:
    //       [x, e] = pulsed_signal(To,sig,par,ton,toff,tr,tf,Fs)
    // where
    //       To:   Duration in s
    //       sig:  Type of carrier, one of the following strings:
    //             "sin","white","pink", "band"
    //       par:  Parameters according to type: 
    //             sin:   [frequency, amplitude]
    //             white: [RMS_value]
    //             pink:  [RMS_value]
    //             band:  [lower_freq, upper_freq, RMS_value]
    //       ton:  Duration of active part in s
    //       toff: Duration of passive part in s
    //       tr:   Fade-in (rise) time in s
    //       tf:   Fade-out (fall) time in s
    //       Fs:   Sample rate in Hz
    //       x:    Pulsed signal
    //       e:    Envelope of pulsed signal (low = 0, high = 1)
    //
    // NOTE: 1) Fade-in and fade-out envelopes are raised cosines,
    //       and their duration is computed from 0 to maximum.
    //       ton and toff is computed from 50 % to 50 % 
    //       2) All cycles are identical, even if the carrier
    //       is noise, so the signal is suitable to excite a moisy
    //       sistem in order to average the responses.
    //       3) Input arguments starting from ton are optional,
    //       being their default values: ton = 0.2, toff = ton, 
    //       tr = 0.02, tf = 0.02, Fs = 44100.
    //
    // -----------------------------
    // Author: Federico Miyara
    // Date:   2020-08-29
    
    // Missing input argument handling
    rhs = argn(2);
    if rhs<8
        Fs = 44100;
    end
    if rhs<7
        tf = 0.02;
    end
    if rhs<6
        tr = 0.02;
    end
    if rhs==4  
        toff = ton;
    end
    if rhs==3  
        ton = 0.2;
        toff = 0.2;
    end    
        
    // Signal (1 cycle)
    Kcycle = round((ton + toff)*Fs);
    tcycle = (0:Kcycle-1)/Fs;
    select sig
    case "sin"
        xcycle = par(2)*sin(2*%pi*par(1)*tcycle); 
    case "white"
        exec whitenoise.sci;
        xcycle = par(1) * whitenoise(ton + toff, Fs);
        xcycle = xcycle(1:Kcycle);
    case "pink"
        exec pinknoise.sci;
        xcycle = par(1) * pinknoise(ton + toff, Fs);
        xcycle = xcycle(1:Kcycle);
    case "band"
        exec bandnoise.sci;
        f1 = par(1);
        f2 = par(2);
        T = ton + toff;
        Xef = par(3);
        N = 32768;
        xcycle = bandnoise(f1, f2, T+N/Fs, Fs, Xef, N);
        xcycle = xcycle($-Kcycle+1:$);
    end
    
    // Envelope (1 cycle)
    // Initialize Envelope
    Ecycle = zeros(xcycle); 
    // Convert characteristic times into indices
    Kon = round(ton*Fs);
    Koff = round(toff*Fs);
    Kr = round(tr*Fs);
    Kf = round(tf*Fs);
    // Obtain the instants limiting analytic regions of
    // envelope
    k1 = 1;
    k2 = Kr;
    k3 = Kr + Kon - round((Kr+Kf)/2);
    k4 = Kon + round((Kr+Kf)/2);
    // Compute analytic regions of envelope 
    Ecycle(k1:k2) = 0.5*(1 - cos(%pi*(0:k2-k1)/Kr));
    Ecycle(k2+1:k3) = 1;
    Ecycle(k3+1:k4) = 0.5*(1 + cos(%pi*(0:k4-k3-1)/Kf));
    
    // Envelope * signal (1 cycle)
    Excycle = Ecycle .* xcycle;
    
    // Compute full length signal
    // Number of cycles
    M = ceil(To/(ton+toff))
    x = zeros(1,M*Kcycle);
    for k=1:M
        x((k-1)*Kcycle+1:k*Kcycle) = Excycle;
    end
    
    lhs = argn(1);
    if lhs==2
        for k=1:M
            e((k-1)*Kcycle+1:k*Kcycle) = Ecycle;
        end
    end          

               
endfunction
