function z = morphfilter(x, W, T, N)
   // Filter with time-varying frequency response 
   // 
   // This function takes an input vector signal x 
   // (either row or column) of arbitrary length 
   // (with the sole limitation of memory size) and 
   // filters it according to a variable frequency 
   // response. 
   //
   // Usage:
   //       z = morphfilter(x, W, T, N)
   //
   // where
   //       x: signal vector
   //       W: filtering matrix; each column contains
   //          the filtering window which corresponds
   //          to the respective discrete time given
   //          by the row vector T
   //       T: vector containing the discrete instants
   //          at which the filtering windows are 
   //          specified
   //       N: FFT size
   //       z: output filtered signal
   // 
   // NOTES: 1) Output signal z is a row or column
   //           vector, according to x. 
   //        2) W must have the same number of columns 
   //           as T, and N/2 or less rows (if less, 
   //           it will be zero-padded to complete N/2). 
   //        3) N and T are optional parameters. If N 
   //           is omitted, it is assumed to be 4096. 
   //           If T is omitted, it is assumed to be
   //           uniformly distributed along the signal
   //           duration.
   //        4) If T does not start at the beginning or
   //           end at the ending of the signal, it is
   //           completed adding an initial and/or last
   //           time, and W is consistently extrapolated.
   //           Note: The timing of the frequency response
   //           must be associated to the signal duration.   
   //        5) Since filtering windows are N/2 points
   //           apart, some interpolation is applied
   //           internally on a frequency basis. To this  
   //           end, the time evolution of the response 
   //           for each frequency (specified at instants 
   //           T(1), T(2), etc.) will be resampled at
   //           equal intervals of size N/2. Attention
   //           should be payed to avoid fast variations
   //           of the frequency response, which could
   //           cause aliasing issues.
   //        6) Note that no value is assumed as sample
   //           rate Fs. Hoewever, the actual spectral
   //           meaning of the filtering windows depends
   //           on Fs. Filter resolution is Fs/N. If, 
   //           for instance, Fs = 44100 Hz and N = 4096, 
   //           Fs/N = 10.767 Hz. Maximum non-aliased 
   //           frequency is Fs/2, which corresponds to 
   //           the (N/2)-th spectral line.
   //
   // Some 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 corresponding to center time of
   // the frame and inversely transformed into filtered
   // frames. Each of them is multiplied by a Hann  
   // window. Even and odd frames are concatenated 
   // independently to form two segments. Finally, both
   // segments are added together after applying an 
   // offset of N/2 samples to the even segment. Note: 
   // this is the overlap-add FFT filtering method.
   //
   // See also fftfilter, trackfilter
   // timevarfilter is similar but deprecatd
   //
   // --------------------------
   // Author: Federico Miyara
   // Date:   2007-12-31
   //         2008-01-07
   //         2021-07-21
   //         2023-10-29
   //         2024-08-15

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

   // If necessary, signal is converted into a column
   // vector
   rowvect = 0;
   if size(x,1) == 1
      x = x';
      // Flag indicating that the original is a 
      // row vector
      rowvect = 1;
   end

   // If argument N is missing, assume N = 4096
   if argn(2) < 4
      N = 4096;
   end

   // Original signal length 
   Lo = size(x,1);

   // Zero-padding of signal so as to have an integer
   // number of frames (actually, semiinteger, but as 
   // frames are overlapped 50 %, it is equivalent to
   // an integer number of frames) 
   // Number of frames
   Q = max([floor(Lo/(N/2)), 1]);
   // Zero-padded signal
   x = [x; zeros((Q+1)*N/2 - Lo, 1)];           
   // Length of zero-padded signal
   L = (Q+1)*N/2;                               
   
   // Number of specified filtering windows
   LT = size(W,2);
   
   // If argument T is lacking, assume T equally
   // spaced along the duration of x. However,
   // as windows have length N, and the first one
   // starts at sample 1 and the last one ends at
   // sample L, then the last one starts at sample
   // L-N+1, so the distance between the starting
   // samples of first and last specified windows 
   // is L-N 
   if argn(2) < 3
      T = [0:LT-1] / (LT-1) * (L-N) + 1;                
   end

   // If frequency response specification is too 
   // short for susequent interpolation, complete
   // it by including a first and/or a last window
   // with a frequency response equal to the first
   // and/or last specified ones. 
   if T(1) > 1
      T = [1, T];
      W = [W(:,1), W];
   end
   if T($) < L - N + 1
      T = [T, L - N + 1];
      W = [W, W(:, $)];
   end

   // If W has less than N/2 rows, complete it with
   // zeros
   W = [W; zeros(N/2-size(W,1), size(W,2))];

   // Get an interpolated matrix with one column for 
   // each frame. Since for each frequency we need 
   // to interpolate along the time axis but interp1
   // interpolates along the rows of each column, we 
   // first need to transpose W and then interpolate
   // the result back.
   WQ = interp1(T, W', [1:Q]*N/2, 'linear')';

   // Frequency response is completed for post-Nyquist
   // frequencies by conjugate symmetrization in order 
   // to get a real response. Each window has now length 
   // N and it can be multiplied point-wise with the FFT 
   // of the column frames
   ww = [WQ; zeros(1,size(WQ,2)); conj(WQ($:-1:2,:))];


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

   // Decomposition into column vector frames 
   xx = [];
   for k = 1:Q
      xx = [xx, x((k-1)*N/2+1:(k-1)*N/2+N)];
   end

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

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

   // Replication of Hann window once for each column 
   // vector 
   // NOTE: The first half of the first window is kept 
   //       equal to 1 and the same holds for the second 
   //       half of the last window. This is to avoid 
   //       fade-in and fade-out effects.
   w = 0.5 - 0.5*cos(2*%pi*(0:N-1)/N);
   hannw = kron(ones(1,Q), w');
   hannw(1:N/2,1) = ones(N/2,1);
   hannw(N/2+1:N,Q) = ones(N/2,1);

   // Application of Hann window to each column
   zz = zz.*hannw;

   // Concatenation of even frames and odd frames  
   // Odd-order frames offset (= 0) 
   z1 = [];                           
   // Even-order frames offset (= N/2)
   z2 = zeros(N/2,1);                
   for k = 1:2:Q
      // odd-order frames segment
      z1 = [z1; zz(:,k)];            
   end

   for k = 2:2:Q
      // even-order frames segment
      z2 = [z2; zz(:,k)];            
   end

   // Zero padding of shorter segment
   if length(z1) < length(z2)
      z1 = [z1; zeros(N/2,1)];
   elseif length(z2) < length(z1)
      z2 = [z2; zeros(N/2,1)];
   end

   // Overlapping and adding
   z = z1 + z2;


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

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

endfunction
