function [y, t] = velvetnoise(T, Nd, Fs)
    // Velvet noise generator
    //
    // Usage:
    //       [y, t] = velvetnoise(T, Nd, Fs)
    //
    // where T:  Duration in s. The actual number of samples
    //           is floor(F*Ts) + 1
    //       Nd: Pulse density (number of pulses per second).
    //           It may be a scaler if it is intended to be
    //           constant or a two component vector if the
    //           density is to increase quadratically from
    //           Nd(1) to Nd(2). Once Nd(2) is reached it 
    //           remains constant until the target duration
    //           is reached 
    //       Fs: Sampling rate in Hz
    //       y:  Velvet noise output row vector
    //       t:  Optional time vector 
    //
    // This function generates velvet random noise of given
    // duration of constant or variable pulse density. No 
    // attempt to get unity RMS is done since that would
    // saturate the audio interface.
    //
    // If Fs is not provided, its default value is 44100. If 
    // Fs = 1 then the number of samples to generate is 
    // floor(T)+1.
    //
    // ------------------------------------
    // Author: Federico Miyara
    // Date:   2023-11-20

    // Provide default values for missing arguments
    rhs = argn(2);
    if rhs < 3
        Fs = 44100;
    end
    if rhs < 2 
        Nd = 1000;
    end
    
    if 1==2
        // For testing purposes --REMOVE--
        T = 3
        Fs = 44100
        Nd(1) = 1/0.01
        Nd(2) = 5000
    end
    
    // Signal length
    L = floor(Fs*T) + 1;
    
    // Velvet noise pulse position
    // Initialization
    k = zeros(1, L);
    if length(Nd)==1
        // Constant pulse density
        // Mean time between pulses in samples
        Td = Fs/Nd;
        // Indices where the pulses are located
        for m=1:T*Nd
            k(m) = round(m*Td + rand()*(Td-1));
        end
    else
        // Variable pulse density
        // Initial and final mean time between pulses
        // in samples
        Td1 = Fs/Nd(1);
        Td2 = Fs/Nd(2);
        // Index of the first pulse
        m = 1;
        k(m) = round(Td1 + rand()*(Td1-1));
        while k(m)<=L
            // Update mean time between pulses 
            Ndt = (k(m)/Fs)^2*Nd(1)^3;
            if Ndt>Nd(2)
                break
            end
            Td = Fs/Ndt; 
            m = m + 1;
            k(m) = k(m-1) + round(Td + rand()*(Td-1)); 
        end
        Td = Fs/Nd(2);
        while k(m)<=L
            m = m + 1;
            k(m) = k(m-1) + round(Td + rand()*(Td-1)); 
        end
    end
    // Limit k to the defined values
    k = k(1:m);
    // Remove values above the target size of signal
    k(find(k>L)) = []; 
    
    // Randomized +1 and -1 sequence
    u = 2*double(rand(1:length(k)) > 0.5) - 1;
    // Initialize y 
    y = zeros(1,L);
    // Apply sequence to the pulse indices
    y(k) = u;
    
    // If required as an output argument, create
    // time vector
    if argn(1)==2
        t = (0:L-1)/Fs;
    end

endfunction
