function z = trackfilter(x, F, T, N, mode, B)
   // Filter whose frequency response tracks the specified 
   // evolution of a frequency component
   //
   // Usage:
   //       z = trackfilter(x, F, T, N, mode, B)
   // Where
   //       x:    signal vector 
   //       F:    frequency vector (as fractions of
   //             sampling rate)
   //       T:    discrete time vector (in samples) 
   //       N:    FFT size
   //       mode: type of filter action (see below)
   //       B:    bandwidth in % of sample rate or
   //             center frequency (See Note 5)
   //       z:    filtered signal vector
   // 
   // NOTES: 
   // 1) Signal vector x may be a row or column vector
   //    and it may have arbitrary length (compatible 
   //    with the available memory). 
   // 2) Frequency vector F may have several values of
   //    frequency distributed along the discrete time
   //    axis.
   // 3) Frequency values are given as fractions of the 
   //    sampling rate. For instance, if Fs = 44100 Hz 
   //    and f = 4410 Hz, then F = 0.1. 
   // 4) T must have the same size as F and it contains
   //    the discrete time values for which F values are
   //    specified. Usually it is uniformly distributed 
   //    along the signal duration, but this is not 
   //    mandatory. 
   // 5) N, mode and B are optional. If N is omitted, it
   //    is assumed to be 4096. It must be any integer 
   //    greater than 127 but the algorithm is much faster
   //    if N is a power of 2. As N  becomes larger, 
   //    frequency discrimination is better, but the 
   //    frequency response cannot vary very quickly. When
   //    N is missing B cannot be provided as argument if
   //    mode is also missing since B is a numeric argument
   //    and there is no way to discriminate if the fifth
   //    argument is N or B. If B is missing it is assumed
   //    to be 5 %. Argument mode must be a string and 
   //    determines how the filter behaves, according to 
   //    the following table:
   //
   //    "lp"       Low pass
   //    "hp"       High pass
   //    "bp"       Single harmonic band pass 
   //    "br"       Single harmonic band reject
   //    "bp%"      Single harmonic percentage band pass 
   //    "br%"      Single harmonic percentage band reject
   //    "bpnnn"    Multiple harmonic band pass
   //    "brnnn"    Multiple harmonic band reject
   //    "bp%nnn"   Multiple harmonic percentage band pass
   //    "br%nnn"   Multiple harmonic percentage band reject
   //
   //    This argument is case insensitive.
   //    In all multiple harmonic options, nnn stands for
   //    up to three decimal digits indicating the number
   //    of harmonics to be passed or rejected. For
   //    instance "pb17" means that 17 harmonics of the
   //    frequencies specified by F (the fundamental 
   //    included) will be allowed at the output, whereas 
   //    "rb8" means that 8 harmonics will be removed. 
   //    Note that "pb1" is equivalent to 'pb'.
   //    The percentage option indicates that the bandwidth
   //    as specified by B is a percentage of the center 
   //    frequency. Otherwise, it is a fixed percentage of
   //    the sampling rate. 
   //    If mode is missing or different from the legal 
   //    values according to the table above it is assumed 
   //    to be "pb%".
   //
   // ------------------------------------------------------
   // Author: Federico Miyara
   // Date: 2008-01-07
   //       2010-05-23
   //       2021-08-29

   // ***********
   // BUREAUCRACY
   // ***********

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

   // If N is missing assume 4096
   if nargin < 4
      N = 4096;
   end

   // If mode is missing or illegal assume 'BP%'
   if nargin < 5 
      mode = 'BP%';
   else 
      if length(mode) < 2
         mode = 'BP%';
      end
   end
   // Take upper case version to simplify parsing
   mode = convstr(mode, "u");

   // If B is missing assume 5 percent
   if nargin < 6
      B = 5;
   end
   // Convert percentage to fraction
   B = B/100;

   // 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 padded signal
   L = (Q+1)*N/2;                       

   // Get an interpolated vector with one frequency value 
   // for each frame. Frequency is normalised to the FFT
   // index scale (N corresponds to Fs). For instance, if 
   // Fs = 44100 Hz and the real frequency is 4410 Hz, 
   // then F = 0.1 and FF = 4096*0.1 = 409.6
   FF = N*interp1(T, F, [1:Q]*N/2, "spline")';

   // Get a filtering matrix according to the filtering
   // behaviour and parameters passed as arguments
   select part(mode, [1 2])
   case 'LP'                   // Low pass
      W = zeros(N/2, Q);
      for i=1:Q
         W(1:round(FF(i))+1, i) = 1;
      end
   case 'HP'                   // High pass
      W = zeros(N/2, Q);
      for i=1:Q
         W(round(FF(i))+1:N/2, i) = 1;
      end
   case 'BR'                   // Band reject
      W = ones(N/2, Q);
      M = harm(mode);
      for k=1:M
         [Fmin, Fmax] = cutoff(k*FF, mode, N, B);
         for i = 1:Q
            W(Fmin(i):Fmax(i), i) = 0; 
         end
      end
   case 'BP'                   // Band pass 
      W = zeros(N/2, Q);
      M = harm(mode);
      for k=1:M
         [Fmin, Fmax] = cutoff(k*FF, mode, N, B);
         for i = 1:Q
            W(Fmin(i):Fmax(i), i) = 1; 
         end
      end 
   end

// Symmetrization of filtering window for real filters
ww = [W; zeros(1,size(W,2)); W($:-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);

// 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)
   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...
// -----------------------------

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

endfunction

// **************************************


// ******************
// INTERNAL FUNCTIONS
// ******************

function M = harm(mode)
   // Internal function for getting the number of harmonics   
   if length(mode) > 2
      if part(mode, 3)=="%"
         M = strtod(part(mode, 4:$));
      else
         M = strtod(part(mode, 3:$));
      end
   end
   if isempty(M) | M==0
      M = 1;
   end
endfunction


function [Fmin, Fmax] = cutoff(FF, mode, N, B)
   // Internal function for getting band cutoff frecuencies
   // Frequencies are normalised so that they represent
   // indexes in the spectral domain
   if length(mode)>2
      if part(mode, 3) == "%"
         // H is the frequency multiplier and divider such
         // that B is the relative bandwidth
         H = sqrt(1 + (B/2)^2) + B/2;
         Fmin = FF / H;
         Fmax = FF * H;
      else
         H = N*B/2;
         Fmin = FF - H;
         Fmax = FF + H;
      end
   end
   // Indexes must be integers
   Fmin = round(Fmin);
   Fmax = round(Fmax);
   // Indexes cannot be smaller than 1 or greater than N/2
   Fmin = max(Fmin, ones(Fmin));
   Fmax = min(Fmax, N/2*ones(Fmax));
endfunction



