function [y, y0, yy] = band_imp_resp(t, w, Q, G)
    // Analytic impulse response of a general band-pass filter  
    // of order n
    //
    // Usage: 
    //        [y, y0, yy] = band_imp_resp(t, w, Q, G)
    //
    // Where
    //        t:    Time vector in s
    //        w:    Vector containing the center frequencies
    //              of the n second-order sections
    //        Q:    Vector containing the values of Q of the n 
    //              second-orden sections
    //        G:    Filter gain at the center frequency 
    //        y:    Column vector containing the time reponse
    //        y0:   Column vector containing the forced response
    //        yy:   Matrix with n columns containing the 
    //              components of the natural response
    //
    // NOTES: 1) Normally vectors w and Q are obtained from an 
    //           analog filter designed, for instance, using the 
    //           function analpf() and a PB --> PBn conversion
    //        2) The n second-order factors have unit coefficients
    //           for the quadratic term. 
    //        3) The input impulse is considered of unit area.
    //        4) The center frequency is computed as the geometric
    //           mean of the center frequencies of the second-order 
    //           sections
    //        5) If the order n is high, numerical errors are 
    //           expected since the algorithm requires computing
    //           the inverse of a square matrix of order 2*n
    //        6) The unit impulse has infinite energy, so a 
    //           high-frequency filter, whose absolute bandwidth
    //           is large, will present a rather high spike 
    //           shortly after 0
    //       
    // Some theory
    // The input signal is
    //       x(t) = delta(t)
    // Its Laplace transform is
    //       X(s) = 1
    // The transfer function is 
    //       H(s) = K*prod(s / (w^2 + (w/Qk)*s + s^2))
    // where K is selected so that at the center frequency
    // the gain is 1. w is the vector of center frequencies. 
    // The Laplace transform of the output signal is
    //       Y(s) = X(s)*H(s) = H(s) = sum(Pk(s))
    // where Pk(s) are the partial fractions corresponding
    // to the natural modes, given by 
    //                       ak + bk*s  
    //       Pk(s) = ------------------------ 
    //                wk^2 + (wk/Qk)*s + s^2 
    // For each one of these the inverse Laplace transform
    // is obtained Laplace, which is the sum of a sine and
    // a cosine modulated by a decreasing exponential.
    // In order to minimize the arithmetic errors when 
    // solving the system of equations to find ak y bk the
    // filter is normalized so that its center angular
    // frequency is approximately 1, which makes the
    // magnitude order of the matrix coefficients more 
    // homogeneous
    //  
    // ------------------------------------------------
    // Author: Federico Miyara
    // Date:   2020-02-28
    //         2021-07-04
    //         2021-08-26
    
    // Input argument error handling
    rhs = argn(2);
    if rhs<4
        G = 1;
    end
    if rhs<3 | argn(2)>4
        msg = gettext("%s: Wrong number of input arguments: %d or %d expected.\n")
        error(msprintf(msg, "band_imp_resp", 3,4))        
    end
    if length(w)~=length(Q)
        msg = gettext("%s: Arguments #%d and #%d: Same sizes expected.\n")
        error(msprintf(msg, "band_imp_resp", 2,3))
    end
    if ~isreal(w)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "band_tone_resp", 2))
    elseif ~and(w>0)
       msg = gettext("%s: Argument #%d: Must be > %d.\n")
       error(msprintf(msg, "band_tone_resp", 2, 0))
    end
    if ~isreal(Q)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "band_tone_resp", 3))
    elseif ~and(Q>0)
       msg = gettext("%s: Argument #%d: Must be > %d.\n")
       error(msprintf(msg, "band_tone_resp", 3, 0))
    end
    if ~isreal(G)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "band_tone_resp", 4))
    end
            
    // Datos de prueba --BORRAR--
    if 1==2
        w0 = 2*%pi*1000;
        w = [1.357783939141098 1 0.73649202923460]*w0;
        Q = [2.961756551889214 1.414213562373100 2.961756551889216];
        G = 1; 
        Fs = 44100;
        t = (0:1/Fs:1)';
        x = [Fs; zeros(t)(1:$-1)];
        ww = 2*logspace(1,4,200)*2*%pi;
    end
    
    // Ensure t is a column vector
    t = t(:);
    
    // Order of the filter (actually, of its skirts since
    // as a dynamic system it has order 2*n)
    n = length(w);
    
    // Estimated center angular frequency
    Wo = prod(w)^(1/n);
        
    // Normalized Laplace variable
    p = poly(0,"p");
    
    // Normalized denominators of the filter's second-order
    // sections
    PH = w.^2/Wo^2 + w./Q/Wo*p + p^2; 
    // NOTE: The normalization respect to Wo will guarantee
    //       that the matrix of the equation system to get 
    //       the components of the partial fractions have 
    //       values with more similar magnitude orders, 
    //       which improves accuracy of the inversa matrix.  
    
    // Normalized transfer function of the filter
    H = p^n/Wo^n/prod(PH);
    // Normalization constant to have unit gain at the 
    // normalized center frequency
    K = 1/abs(horner(H,%i*Wo/Wo));
    H = K*H;
      
    // Factors of the denominator of the transform of
    // the output
    P = PH;
    
    // (n-1)-factor products 
    for k=1:n
        R(k) = prod([P(1:k-1),P(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 = [];
    for k=1:n
        CC = [CC, [C(k,:),0]', [0, C(k,:)]'];
    end
    
    // Independent vector corresponding to the coefficients
    // of the numerator of the response presented as column
    // vector. All frequencies must be unnormalized
    D = zeros(2*n,1);
    D(n+1) = G*K/Wo^n;
    
    // Calculation of the coefficient vector
    AB = inv(CC)*D;
        
    // Coefficients of the independent term of the numerator 
    // of the partial fractions
    a = AB(1:2:$)*Wo^2;
    // Coefficientes of the term in s of the numerator of 
    // the partial fractions
    b = AB(2:2:$)*Wo;
    
    if 1==2
        // Verification
        // Compare normalized partial fraction expansion
        // HHnor with normalized transfer function HHHnor
        apri = AB(1:2:$)
        bpri = AB(2:2:$)
        HHnor = sum((apri + bpri*p)./PH')
        HHHnor = G*K/Wo^n*p^n / prod(PH)
        // Compare unnormalized partial fraction expansion
        // HH with original transfer function HHH
        HH = sum((a + b*%s)./(horner(PH, %s/Wo)'*Wo^2))
        HHH = G*K*%s^n / prod(horner(PH, %s/Wo)*Wo^2)
    end
    
    
    // Forced response
    y0 = zeros(t); 
    
    // In order to denormalize the response we use the
    // natural angular frequencies wo obtained from the
    // orginal w
    wo = w.*sqrt(1 - 1/4./Q.^2);

    // Natural modes
    yy = zeros(length(t), n); 
    for k=1:n
        yy(:, k) = exp(-w(k)/2/Q(k)*t) .* ...
                   (b(k).* ...
                   cos(wo(k)*t) + ...
                   (a(k)/wo(k) - b(k)/2/Q(k)*w(k)/wo(k)) * ...
                   sin(wo(k)*t));
    end
             
    // Complete response
    y = sum(yy, 2);
    
    if 1==2
        // Impulse response energy
        Y2 = sum(y.^2)/Fs;
        N = 10000
        Wmax = 100*w0
        www = (0:N-1)*Wmax/N;
        // Fourier transform of impulse respose
        // which equals the transfer function
        // computed on the imaginary axis
        Y = abs(horner(H, %i*www/w0));
        // Plot absolute valeue of Fourier transform
        figure(101);
        gcf().background = 8;
        plot(www, Y)
        // Energy computed in the frequency domain
        // using Plancherel's identity. Integration
        // is performed numerically over a finite
        // but sufficiently large interval.
        Y2a = 1/%pi*sum(Y.^2)*Wmax/N;
    end
        
    if 1==2
        figure(100);
        scf(100);
        clf();
        gcf().background = 8;
        plot(t,[x/Fs, yy, y])
    end
    
endfunction

