function y = envelope(x, tau, Fs, n, graph, fig)
    // Envelope detection
    // 
    // Detects the energy envelope of a signal by low-pass
    // filtering the squared signal and then taking its
    // square root. It is equivalent to the short-time RMS 
    // value of the signal, also called exponential average.
    //
    // Usage:
    //       y = envelope(x, tau, Fs, n, graph, fig)
    //
    // where x:     input signal
    //       y:     output signal (envelope of x) 
    //       tau:   time constant of detector in s
    //       Fs:    sampling rate in Hz
    //       n:     filter order 
    //       graph: graphic output ("nof", "lin", "log", "dB")
    //       fig:   number of figure where output is plotted
    //
    // fig, graph, n, Fs and tau are optional arguments. Their
    // default values are fig = gcf().figure_id+1, "nof" (no 
    // graphic output), n = 1, Fs = 44100, tau = 0.01. The 
    // argument graph affects the magnitude axis of the 
    // envelope graph. If "dB", the reference value is 1.
    // The envelope is computed by filterig the squared signal 
    // with a Butterworth filter of cutoff angular frequency 
    // 1/tau and order n. Typically n = 1 but a higher order
    // might be used to smooth out residual high frequency 
    // ripple, at the expense of slightly increased lag caused
    // by group delay.
    // The envelope retains the original sampling rate. It 
    // may be downsampled to a lower sampling rate by means 
    // of the function intdec().
    // 
    // --------------------------------------------------
    // Author:   Federico Miyara
    // Date:     2006-07-23
    // Revision: 2006-08-14
    // Revision: 2007-11-24 
    // Revision: 2007-12-10 
    // Revision: 2008-02-26
    // Revision: 2008-08-04
    // Revision: 2009-09-14
    // Revision: 2010-06-07
    // Revision: 2020-03-28
    // Revision: 2022-04-01 
    
    // Default arguments if not given
    rhs = argn(2);
    if rhs < 6
       fig = gcf().figure_id + 1;
    end
    if rhs < 5
       graph = "nof";
    end
    if rhs < 4
       n = 1;
    end    
    if rhs < 3
        Fs = 44100;
    end
    if rhs < 2
        tau = 0.01;
    end

    // Find low-pass filter denominator and numerator
    // coefficients using bilinear transformation of 
    // a Butterworth lowpass prototype    
    hs = analpf(n,"butt",[0 0], 1/tau);
    hz = horner(hs, 2*Fs*(%z-1)/(%z+1));
    Af = coeff(hz.den)($:-1:1);
    Bf = coeff(hz.num)($:-1:1);

    // Filter squared signal and take square root
    y = sqrt(filter(Bf, Af, x.^2));

    // Graphic output 
    if rhs >= 5
        t = [0:length(y)-1]/Fs;
        if rhs == 5 
            figure;
            gcf().background = -2;
        else
            figure(fig);
            gcf().background = -2;
        end
        select graph
        case 'lin'
            plot2d(t, y)
            xlabel('t  [s]')
            ylabel('y')
        case 'log'
            plot2d("nl", t, y)
            xlabel('t  [s]')
            ylabel('log(y)')
        case 'dB'
            plot2d(t, 20*log10(y))
            xlabel('t  [s]')
            ylabel('y  [dB]')
        end
        title("Envelope (tau = " + string(tau) + "s)")
        gca().grid = [1 1];
    end

endfunction

