function [y] = bandfilter(x, fo, typ, Fs, nd, notch)
    // Several types of bandpass filtering
    // This function takes an audio signal and filters it according
    // to a pass band of specified center frequency and type. It 
    // can also be used as a notch filter, filtering out the 
    // specified band
    // 
    // Usage:
    //       y = bandfilter(x, fo, typ, Fs, nd, notch)
    //
    // where x:     vector containing input signal
    //       fo:    center frequency in Hz or two-component vector
    //              containing lower and upper band limits fmin,
    //              fmax; in this case typ is irrelevant
    //       typ:   critical (0), octave (1), one-third octave (3), etc.
    //       Fs:    sampling rate 
    //       nd:    no-downsampling flag; if equal to 1 no downsampling
    //              is applied to low frequency bands
    //       notch: notch flag; if equal to 1 instead of a bandpass
    //              filter a notch or bandreject filter is obtained 
    //       y:     filtered otput
    //      
    // Fs, typ, nd, notch are optional, being their default values 
    // 44100, 1 (octave), 0 (enable downsampling) and 0 (bandpass). 
    // fo should be smaller than Fs/2.
    //
    // See also: fftfilter, zkrab, zbark, fcritband, critenv,
    //           trackfilter, critband
    //
    // --------------------------
    // Author: Federico Miyara
    // Date: 2008-07-30
    //       2009-02-11
    //       2009-11-16
    //       2016-03-18
    //       2019-03-17 
    //       2020-03-26

    // Dealing with arguments

    rhs = argn(2);

    // If less than 2 or more than 6 input arguments exit with error 
    // message
    if rhs < 2 | rhs > 6
        error(msprintf(gettext(..
        "%s: Wrong number of input arguments: %d to %d expected.\n"), ..
        "bandfilter", 2, 4));
    end

    // If x is not a vector exit with error message
    if ~isvector(x) 
        error(msprintf(gettext(..
        "%s: Wrong size for argument #%d: Vector expected.\n"), ..
        "bandfilter", 1));
    end

    // If fo is not a scalar or two-element real vector exit with
    // error message
    if length(fo) < 1 | length(fo) > 2
        error(msprintf(gettext("%s: Argument #%d: %s or %s expected.\n"), ..
        "bandfilter", 2, "scalar", "vector of length 2"));
    end

    // If notch is not given, it is set by default to 0
    // (normal bandpass filter)
    if rhs<6
        notch = 0;
    end

    // If nd is not given, it is set by default to 0
    // (enable downsampling for low frequency bands)
    if rhs<5
        nd = 0;
    end

    // If Fs is not given, it is set by default to 44100 Hz
    if rhs<4
        Fs = 44100;
    end

    // If typ is not given, it is set by default to 1 (octave band)
    if rhs<3
        typ = 1;
    end

    // Signal size and length
    xsiz = size(x);
    xlen = length(x);

    // Determination of low and high cut-off frequencies,
    // finf and fsup
    if length(fo)==2
        finf = fo(1);
        fsup = fo(2);
    else
        if typ == 0
            // If it is a critical band 
            exec("zkrab.sci"); 
            exec("zbark.sci"); 
            finf = zkrab(zbark(fo) - 0.5);
            fsup = zkrab(zbark(fo) + 0.5);
        else
            // If it is a fractional octave band
            finf = fo / 2^(1/typ/2);
            fsup = fo * 2^(1/typ/2);
        end
    end

    if finf > Fs/2 then
        // If band is completely above Nyquist frequency
        mprintf("fmin = " + string(finf) + " Hz" + "\n")
        mprintf("Band exceeds Nyquist frequency")
        return    
    end

    if fsup > Fs/2
        // If band is partially above Nyquist frequency 
        mprintf("fmax = " + string(fsup) + " Hz" + "\n")
        mprintf("Upper frequency exceeds Nyquist frequency. ")
        mprintf("Band will be limited to Nyquist" + "\n" )
        fsup = Fs/2;
    end

    // Load fftfilter
    if ~exists("fftfilter")
        exec("fftfilter.sci");
    elseif type(fftfilter)<>13
        mprintf("fftfilter exists as a non-function variable. " + ... 
        "It has been saved as %%fftfilter")
        %fftfilter = fftfilter;
        exec("fftfilter.sci");
    end

    // Computation is performed by first downsampling the signal,
    // as this turns out to be computationally more efficient.
    // However, the signal is only downsampled if its upper
    // cut-off frequency is at least 6 times lower than the
    // original sampling rate. This ensures that with the
    // minimum downsampling (by a factor of 2) the the new
    // Nyquist frequency will be at least 1.5 times higher than 
    // the upper cut-off frequency, which in turn guarantees
    // that the error due to the antialiasing filter contained 
    // in the resampling function (intdec) will be small (it is
    // about 0.3 % in the worst case).
    // Even if the downsampling could be done by any integer 
    // factor, being the maximum possible
    //      q = floor(Fs/3/fsup),
    // it is preferable that it be by a power of 2 since the 
    // downsampled frequencies belong to the set of original
    // frequencies.
    // The factor can be obtained from the condition:
    //      q = 2^r <= Fs/(3*fsup)
    // Solving for the largest integer exponent r,
    //      r <= log2(Fs/(3*fsup))
    //      r = floor(log2(Fs/(3*fsup)))
    //      q = 2^r = 2^floor(log2(Fs/(3*fsup)))
    // EXAMPLE 1: Original sampling rate, 44100 Hz. For an upper
    // cut-off frequency of 7350 Hz it turns out that the signal
    // is downsampled by a factor of 2, i.e., the new sampling 
    // rate is 22050 Hz. The new Nyquist frequency is, thus,
    // 11025 Hz, that is exactly 1.5 times the upper cut-off
    // frequency of the filter.
    // EXAMPLE 2: Same sampling rate. If the upper cut-off
    // frequency is 400 Hz, it would be possible to downsample 
    // to 1200 Hz. A downsamplig factor of 36 would provide a
    // sampling rate of 1225 Hz, but it is not a power of 2. 
    // The chosen factor, 32, yields a sampling rate of 
    // 1378.125 Hz, which means that the Nyquist frequency is 
    // 1.72 times above the upper cut-off frequency, which is 
    // above the minimum acceptable factor of 1.5

    // FFT size
    // For Fs = 44100, an FFT size of 4096 is adequate to 
    // represent both the audible spectral resolution and a
    // reasonable time resolution. Here the FFT size is 
    // chosen to have an approximately equivalent resolution.
    // Consequently, the effect of reducing Fs is, mainly, to 
    // cut the high frequencies out without affecting too much
    // the spectral resolution.
    N = 2^ceil(log2(4096*Fs/44100));

    // The downsampling process can be skipped setting nd = 1,
    // gaining filtering quality at the expense of computational 
    // speed
    if nd ~= 1 
        if fsup < Fs/6   
            p = 1;
            // The downsampling factor is chosen as a power of 2
            q = 2^floor(log2(Fs/3/fsup));
            x = intdec(x, p/q);
            Fs = Fs/q;
            // FFT size is reduced so that the absolute duration 
            // of the FFT window is not much larger than the 
            // original one
            N = N/q;
        else
            // As the final upsampling process needed to get back
            // to the original Fs is performed only if q > 1, 
            // q needs to be defined anyway. 
            q = 1;
        end
    else
        q = 1;
    end

    // For very low frequency bands increase the FFT size
    // to increase frequency resolution, at the expense of
    // time resolution (time resolution is not much of a 
    // requirement at very low frequecies)
    if fsup < 20
        N = 2*N;
    end

    // Design of the bandpass filtering window, which is 1
    // for the samples of the spectrum between finf and fsup
    // and 0 elsewhere (note that "elsewhere" may be much 
    // less than 20000 Hz in the case of low frequency bands 
    // since the signal has been downsampled)
    // NOTE: If the original band exceeds the Nyquist 
    // frequency, the general formula for the upper frequency
    // index must be adjusted to N/2 
    w = zeros(1, N/2);
    //n1 = round(finf/Fs*N) + 1;
    //n2 = min(round(fsup/Fs*N) + 1, N/2);
    n1 = ceil(finf/Fs*N) + 1;
    n2 = min(floor(fsup/Fs*N) + 1, N/2);

    // disp ([N,n1,n2])

    w(n1:n2) = 1;

    if notch==1 
        // If it is a reject band (notch) filter
        w = 1 - w;
    end

    // [(0:N/2-1)'*Fs/N, w']
    // plot([0:N/2-1]*Fs/N, w)

    // A continuous FFT filtering is performed 
    y = fftfilter(x, w, N);

    // If the signal has been downsampled, the filtered signal
    // must be oversampled in order to recover the original
    // sampling rate

    if q > 1
        // Upsample output signal;
        y = intdec(y, q/p);
        // The result of upsampling a downsampled signal by the 
        // same factor may have a slightly different length.
        // It is necessary to adjust the difference.
        // NOTE: Scilab 6.1 has solved this bug so this may be 
        //       unnecessary
        ylen = length(y);
        if ylen > xlen
            // If longer than the original, truncate
            y = y(1:xlen);
        else
            // If shorter, it must be zero-padded, but care
            // must be taken since the signal may be a row
            // or column vector 
            ysiz = size(y);
            if ysiz(1) == 1
                // It is a row vector
                y = [y, zeros(1, xlen-ylen)];
            else
                // It is a column vector
                y = [y; zeros(xlen-ylen,1)];               
            end
        end
    end

endfunction


