function [z] = fftfilter(x, w, N)
    // FFT-filtering procedure
    // This function takes an input vector signal x (either row
    // or column) of arbitrary length and filters it according
    // to the frequency response given by the filtering window w,
    // assuming an FFT size N.
    //
    // Usage
    //     z = fftfilter(x, w, N)
    //
    // where
    //     x: Input signal as a column or row vector 
    //     w: Filtering window
    //     N: Size of the FFT 
    //     z: Filtered signal       
    //
    // Output signal z is a row or column vector according to x.
    // If N is omitted, it is assumed to be 4096. No value is
    // assumed as sampling rate Fs but, of course, the actual
    // spectral meaning of the filtering window depends on Fs.
    // If w is longer than N it is truncated.
    // If w has a length N/2 or less, it is first zero-padded to
    // N/2 samples and then it is symmetrically completed so that
    // it works as a real filter.
    // If N/2 < length(w) < N, w is zero-padded to complete
    // N samples; in this case the filter response is complex.
    // Filter resolution is Fs/N. If, for instance, Fs = 44100 Hz
    // and N = 4096, Fs/N = 10.7666 Hz. Maximum non-aliased
    // frequency is Fs/2, which corresponds to the (N/2)-th
    // spectral line.
    //
    // Technical stuff
    //
    // The signal x is segmented into frames of size N with an
    // N/2-sample overlapping. Each one is FFT transformed, then
    // each FFT is multiplied by the filtering window and inversely
    // transformed into filtered frames. Each of them is multiplied
    // by a Hann window. Odd frames and even frames are independently
    // concatenated to form two segments. Finally, both segments are
    // added together after applying an offset of N/2 samples to the 
    // even segment. 
    // Due to memory issues it is not recommended to filter signals 
    // lasting in excess of one or two minutes.
    //
    // ----------------------------------------
    // Authors: Federico Miyara, David Giardini, Ernesto Accolti
    // Date:    2007-12-19
    //          2008-07-31
    //          2009-08-10
    //          2010-03-28
    //          2010-03-30 (EA)
    //          2010-04-01
    //          2020-02-28 

    // BUREAUCRACY
    // -----------

    // Output variables initialisation (not found in input variables)
    //z=[];
    // Display mode
    //mode(0);
    // Display warning for floating point exception
    //ieee(1);

    // Number of arguments in function call
    [lhs, rhs] = argn()

    // If not enough arguments, assume N = 4096
    if rhs == 2 then
        N = 4096;
    end;
    if rhs<2 | rhs>3
        msg = gettext("%s: Wrong number of input arguments: %d or %d expected.\n");
        error(msprintf(msg, "fftfilter", 2,3))
    end

    // If necessary, signal is converted into a column vector
    rowvect = 0;
    if size(x,1) == 1 then
        x = x.';
        rowvect = 1;
    end;

    // If necessary, filtering window is converted into a
    // column vector
    if size(w,1) == 1 then
        w = w.';
    end;

    // Signal length
    Lo = max(size(x));

    // Zero-padding of signal to have an integer number of frames
    Q = max([floor(Lo/(N/2)),1]);        // Number of frames
    x = [x; zeros((Q+1)*N/2 - Lo, 1)];   // Zero-padded signal
    L = (Q + 1)*N/2;                     // Length of padded signal

    // Dealing with filtering window mismatches
    Lw = max(size(w));
    // If w is longer than N it is truncated 
    if Lw > N then
        w = w(1:N);
    end
    if N/2 < Lw then    
        // Zero-padding of incomplete non-symmetric filtering windows
        // In this case the filter is complex
        w = [w; zeros(N-Lw, 1)];
        Lw = N;
        // Set real-filter flag to 0
        re = 0;
    else
        // Symmetrisation of filtering window for real filters
        w = [w; zeros(N/2 - length(w), 1)];
        w = [w; 0; conj(w($:-1:2))];
        // Set real-filter flag to 1
        re = 1;
    end;

    // ALGORITHM
    // ---------

    // Decomposition into column vector frames
    // NOTE: Calculation speed increases if the memory
    //       space required for the entire matrix xx
    //       is pre-allocated
    xx = zeros(N,Q);
    for k = 1:Q
        xx(:,k) = x((k-1)*N/2+1:(k-1)*N/2+N);
    end;

    // FFT spectra of frames
    //yy = mtlb_fft(xx,N);
    yy = fft(xx, -1, 1);

    // Replication of filtering window once for each
    // column vector using Kronecker product
    ww = ones(1,Q) .*. w;

    // Multiplication of frames by the filtering window
    // and application of the inverse FFT transform
    zz = fft(yy .*ww, 1, 1);

    // Concatenation of odd frames and even frames

    // Pre-allocate vector of odd frames for increased speed
    z1 = zeros(L,1);
    // Pre-allocate vector of even frames for increased speed
    z2 = z1;

    // Hann window for cross-fade between overlapping
    // frames
    hann = 0.5*(1 - cos(2*%pi*[0:N-1]/N))';

    // The first half of the first window is set equal
    // to 1. Similarly, the second half of the last
    // window is set to 1. This is to avoid fade-in and
    // fade-out effects.
    hann1 = [ones(N/2,1); hann(N/2+1:N)];
    hann2 = [hann(1:N/2); ones(N/2,1)];

    // First frame
    z1(1:N) = zz(:,1).*hann1;

    // Odd-order frames; offset = 0
    for k = 3:2:Q-1
        z1((k-1)*N/2+1:(k-1)*N/2+N) = zz(:,k).*hann;   
    end;

    // Even-order frames; offset = N/2
    for k = 2:2:Q-1
        z2((k-1)*N/2+1:(k-1)*N/2+N) = zz(:,k).*hann;   
    end;


    // Last frame (it belongs to z1 if Q is odd and to
    // z2 if Q is even)
    if Q/2==floor(Q/2) then
        z2((Q-1)*N/2+1:(Q-1)*N/2+N) = zz(:,Q).*hann2;
    else
        z1((Q-1)*N/2+1:(Q-1)*N/2+N) = zz(:,Q).*hann2;
    end   

    // Overlapping
    z = z1 + z2;

    // A FINAL BIT OF BUREAUCRACY...
    // -----------------------------

    // Restore original signal length
    z = z(1:Lo);

    // If x was a row vector, z is converted back into a row vector
    if rowvect==1 then
        z = z.';
    end;

    // If the filter must be real, keep only the real part
    // (hopefully the imaginary part will be vanishingly small,
    // but it will mess things up if not removed)  
    if re==1 then
        z = real(z);
    end;

endfunction
