function [y, y0, yy] = filter_imp_resp(t, pp, zz, K)
    // Impulse response of a general filter
    //
    // Usage:     
    //        [y, y0, yy] = filter_imp_resp(t, pp, zz, K)
    //
    // Where
    //        t:    Time row vector in s
    //        pp:   Vector containing the poles  
    //        zz:   Vector containing the zeros
    //        K:    Constant when the coefficients
    //              of the maximum powers of the 
    //              numerator and denominator are 1   
    //        y:    Row vector containing the time
    //              reponse
    //        y0:   Row vector containing the forced
    //              response
    //        yy:   Matrix with n rows, each one 
    //              containing the components of 
    //              the natural response
    //
    // NOTES: 1) The number of zeros must be smaller 
    //           than the number of poles. If it is 
    //           needed to get the impulse response of
    //           a high-pass filter, which contains 
    //           the same number of poles and zeros,
    //           first divide both polynomials to get 
    //           a constant plus a filter complying
    //           with this condition. Then apply this 
    //           function. The constant will produce
    //           the constant times a Dirac delta,
    //           which can only be represented as a 
    //           discrete time signal if the sample 
    //           rate is known, in which case its 
    //           non-aliased version is equal to a
    //           unit pulse times Fs.
    //        2) If there are conjugated complex poles
    //           or zeros, the result should be real,
    //           but because of arithmetic errors, a
    //           small imaginary part could appear. 
    //           In such case the real part should be
    //           taken.
    //        3) The poles must be different. This
    //           implementation does not support
    //           multiple poles. 
    //      
    // -----------------------------------------------
    // Author: Federico Miyara
    // Date:   2021-07-18
    //         2025-10-08
    
    // Input arguments
    rhs = argn(2);
    
    if rhs<4
        K = 1;
    end
    if rhs<3
        zz = [];    
    end
    if rhs<2 | rhs>4
        msg = gettext("%s: Wrong number of input arguments: %d to %d expected.\n");
        error(msprintf(msg, "filter_imp_response", 2,4))
    end
    
    if prod(size(t))<>max(size(t)) 
       msg = gettext("%s: Argument #%d: Vector expected.\n");
       error(msprintf(msg, "filter_imp_response", 1))
    else
       if type(t)==1
          if isreal(t)
          else
             msg = gettext("%s: Argument #%d: Decimal number expected.\n");
             error(msprintf(msg, "filter_imp_response", 1))
          end            
       else
          msg = gettext("%s: Argument #%d: Decimal number expected.\n");
          error(msprintf(msg, "filter_imp_response", 1))
       end        
    end

    if prod(size(pp))<>max(size(pp))  
       msg = gettext("%s: Argument #%d: Vector expected.\n");
       error(msprintf(msg, "filter_imp_response", 2))
    else
       if type(pp)<>1    
          msg = gettext("%s: Argument #%d: Decimal or complex number expected.\n");
          error(msprintf(msg, "filter_imp_response", 2))
       end
    end
    
    if prod(size(zz))<>max(size(zz))  
       msg = gettext("%s: Argument #%d: Vector expected.\n");
       error(msprintf(msg, "filter_imp_response", 3))
    else
       if type(zz)<>1    
          msg = gettext("%s: Argument #%d: Decimal or complex number expected.\n");
          error(msprintf(msg, "filter_imp_response", 3))
       end
    end
    
    if length(K)<>1
       msg = gettext("%s: Argument #%d: Scalar (1 element) expected.\n");
       error(msprintf(msg, "filter_imp_response", 4))
    else
       if type(K)<>1
          msg = gettext("%s: Argument #%d: Decimal or complex number expected.\n");
          error(msprintf(msg, "filter_imp_response", 4))
       end
    end


    // For testing purposes --delete--
    if 1==2
        pp = [-4:-1]*10
        zz = [0 0] 
        K = 1
        Fs = 44100;
        t = 0:1/Fs:1;
        x = [Fs,zeros(t)(1:$-1)];
        ww = 2*logspace(1,4,200)*2*%pi;    
    end
       
    // Degree of numerator and denominator
    m = length(zz);
    n = length(pp);
    
    if m>=n
        msg = "%s: The degree of the numerator must be smaller"
        msg = msg + " than the degree of the denominator.\n";
        error(msprintf(msg, "filter_imp_response"))
    end
    
    // Normalization angular frequency
    Wo = abs(prod(pp)^(1/n));    
    
    // Normalized Laplace variable
    p = poly(0,"p");
        
    // Normalized numerators of the filter's first-order
    // sections
    NN = p - zz/Wo;     
    
    // Normalized denominators of the filter's first-order
    // sections
    DD = p - pp/Wo;     
    
    // Normalized transfer function of the filter
    H = K/Wo^(n-m)*prod(NN)/prod(DD);
    
    // (n-1)-factor products 
    for k=1:n
        R(k) = prod([DD(1:k-1), DD(k+1:$)]);
    end

    // Coefficient matrix of the (n-1)-factor products.
    // Each row contains the coefficients of each
    // polynomial obtained as the product of n-1
    // second-order factors in ascending power orden 
    C = coeff(R);
    
    // System matrix to obtain the coefficients of the 
    // partial fraction decomposition
    CC = C.';

    // Independent vector corresponding to the coefficients
    // of the numerator of the response presented as column
    // vector. All frequencies must be unnormalized
    D = zeros(1, n);
    D(1:m+1) = K/Wo^(n-m)*coeff(prod(p - zz/Wo));    
    
    // Partial fraction coefficients 
    A = inv(CC)*D';
    
    // Forced response
    y0 = zeros(t);
    
    // Components of the free response
    yy = zeros(n,length(t));
    for k = 1:n
        yy(k,1:length(t)) = A(k)*Wo*exp(pp(k)*t);
    end 
    
    // Total response
    y = y0 + sum(yy,1);
       
    // Testing
    // sum(A'./DD)
    
endfunction
