function y = timevarconv(x, h, d)
   // Convolution with a time-variable impulse response
   //
   // Usage:
   //       y = timevarconv(x, h, d)
   // where
   //       x: Input signal vector
   //       h: Variable finite impulse response matrix
   //          or cell array
   //       d: Delay vector of impulse responses, in
   //          samples
   //       y: Output signal vector
   //
   // Impulse responses may be given either as an Lx*H 
   // matrix, where Lx is the input length and H the 
   // common length of all impulse responses or a 
   // cell-array with Lx components, each of which  
   // contains the impulse response corresponding to
   // each discrete time instant. In this case the
   // length of the impulse responses may differ.
   //
   // The delay vector is optional. If present its
   // size must be Lx. Positive entries mean that the
   // corresponding impulse responses are delayed, 
   // negative entries mean they are anticipated. If
   // the signal is a time signal, this means the 
   // response is non-causal. For other independent 
   // variables it just means a sort of blurring of
   // the signal.
   //
   // NOTE: Since the computation is done directly,
   //       for long impulse responses it is rather
   //       slow but in exchange, the function 
   //       allows changing the impulse response 
   //       at each sample  
   //
   // ------------------------
   // Author: Federico Miyara
   // Date:   2025-06-05
   //         2025-06-17
   //         2025-06-27
      
   // Length of signal
   Lx = length(x);
   
   // Missing delay
   if argn(2)<3
      d = zeros(1,Lx);
   end
   
   // Convert x to row
   if size(x,2)==1
      wascolumn = 1;
      x = x.';
   else
      wascolumn = 0;
   end
   
   // Compute length of impulse responses
   if iscell(h)
      // If h is a cell array get the length
      // of each cell
      Jo = zeros(1,Lx);
      for k=1:Lx
         Jo(k) = length(h{k});
      end
   else
      // If h is a matrix get the common length
      // to all impulse responses   
      Jo = size(h,2)*ones(1,Lx);
   end
   
   // Total length of convolution
   // Corresponds to the longest combination of
   // k + Jo(k)
   Ly = max((1:Lx) + Jo) - 1;
   
   // Example --delete--
   // d = -round(Jo/4);
   
   // Maximum anticipation
   // Corresponds to the maximum displacement
   // to the left of k + d(k) 
   A = 1 - min((1:Lx) + d);
   A = max(A,0);
   
   // Prepend zeros to x in order to prevent
   // indices smaller than 1  
   xx = [zeros(1, A), x];
   
   // Initialize y including the prepended
   // zeros
   y = zeros(1, Ly+A);
   
   // Index from which the impulse response
   // will be applied
   n = A + (1:Lx) + d;
   // Update y by adding only the Jo(m) values 
   // corresponding to h{m} or h(m,1)
   if iscell(h)
      for m=1:Lx
         y(n(m):n(m)+Jo(m)-1) = y(n(m):n(m)+Jo(m)-1) + xx(A+m)*h{m};
      end
   else
      for m=1:Lx
         y(n(m):n(m)+Jo(m)-1) = y(n(m):n(m)+Jo(m)-1) + xx(A+m)*h(m,1);
      end
   end
   
   // Remove the first A components of y   
   y = y(A+1:$);
   
   // If x was a column vector ensure y is also
   // a column vector
   if wascolumn==1
      y = y(:);
   end
  
endfunction
