function [y, y0, yy] = biquad_imp_resp(t, w, Q, E, F, G)
    // Analytic impulse response of a general filter of order 2*n 
    //
    // Usage: 
    //        [y, y0, yy] = biquad_imp_resp(t, w, Q, E, F, G)
    //
    // Where
    //        t:    Time row vector in s
    //        w:    Vector containing the center frequencies of the n
    //              biquad sections
    //        Q:    Vector containing the values of Q of the n biquad 
    //              sections
    //        E:    Vector containing the n independent terms of the  
    //              biquad section numerators
    //        F:    Vector containing the first degre coefficients of
    //              the biquad section numerators
    //        G:    Vector contining the quadratic coefficients of the
    //              biquad section numerators  
    //        y:    Row ector containing the time reponse
    //        y0:   Row ector containing the forced response
    //        yy:   Matrix with n rows, each containing the components 
    //              of the natural response
    //
    // NOTES: 1) Normally vectors w, Q, E, F, G, are obtained from 
    //           a constant resistance network design, but they can
    //           also be obtained from an analog filter design, for 
    //           instance, using function analpf() and a LP --> BP 
    //           or LP --> HP conversion. They can also be obtained
    //           from published tables.
    //        2) The n biquad factors have unit coefficients for the
    //           quadratic term in the denominator. 
    //        3) The input impulse is considered of unit area. The
    //           output contains an impulse of area prod(G).
    //        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 requieres to compute the inverse
    //           of a square matrix of order 2*n
    //       
    // Some theory
    // The input signal is
    //       x(t) = delta(t)
    // Its Laplace transform is
    //       X(s) = 1
    // The transfer function is 
    //                    /   E + F*s + G*s^2)    \
    //       H(s) =  prod( ---------------------- )
    //                   \  w^2 + (w/Q)*s + s^2  /
    // where w is the vector of center frequencies and Q is the
    // vector of quality factors of the sections. 
    // The Laplace transformaofthe 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, 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 frecquency is approximately 1,
    // which makes the magnitude order of the matrix coefficients 
    // more homogeneous.
    //  
    // ------------------------------------------------
    // Author: Federico Miyara
    // Date:   2020-03-19
    //         2023-10-31
    
    // Input argument error handling
    rhs = argn(2)
    if rhs<6
        // Default quadratic coefficient
        G = 0;
    end
    if rhs<5
        // Default linear coefficient
        F = 0;
    end
    if rhs<4
        // Default independent term (low-pass filter)
        F = 1;
    end
    if rhs<3 | rhs>6
        msg = gettext("%s: Wrong number of input arguments: %d to %d expected.\n")
        error(msprintf(msg, "band_imp_resp", 3,6))        
    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

            
    // Datos de prueba --BORRAR--
    if 1==2
        // Parámetros del ecualizador para el TDH39
        R = 600;
        L1 = [2.65 0.0318];         L2 = [0.00669 0.00795];
        C1 = [18.4e-9 22.1e-9];     C2 = [7.3e-6 88.4e-9];
        R1 = [10070 2774];          R2 = [36 130];
        // Funciones de transferencia de cada biquad 
        H = 1 + R2.*C2*%s + L2.*C2*%s^2;
        H = H ./ (1 + (R2 + R).*C2*%s + L2.*C2*%s^2);
        // Parámetros delos denominadores
        w = 1 ./ sqrt(coeff(H.den)(5:6))
        Q = 1 ./ coeff(H.den)(3:4) ./ w
        // Coeficientes de los numeradores
        E = coeff(H.num)(1:2) .* w.^2
        F = coeff(H.num)(3:4) .* w.^2
        G = coeff(H.num)(5:6) .* w.^2

        Fs = 44100;
        t = 0:1/Fs:1;
        x = [Fs,zeros(t)(1:$-1)];
        ww = 2*logspace(1,4,200)*2*%pi;
        
        figure(1125);
        clf();
        gcf().background = 8;
        plot2d("ln",ww/2/%pi,-20*log10(abs(horner(prod(H),%i*ww))))
        
    end
    
    // Order of the filter (actually, of its skirts)
    n = length(w);
    
    // Estimated center frequency
    Wo = prod(w)^(1/n);
        
    // Normalized numerators of the filter's biquads
    NN =  E/Wo^2 + F/Wo*%s + G*%s^2;
    // Normalized denominators of the filter's biquads
    DD = w.^2/Wo^2 + w./Q/Wo*%s + %s^2; 
    // NOTE: The normalization respect to Wo will guarantee
    //       that the matrix of the equation system for geting 
    //       the components of the partial fractions have
    //       values with more similar magnitude orders, which
    //       improves accuracy of the inversa matrix.  
        
    // Transfer function of the filter
    N = prod(NN);
    D = prod(DD);
    Ho = N/D;
    
    // Normalization constant to have unit gain at the normalized 
    // center frequency (NOT USED--AND NOT USEFUL!)
    K = 1/abs(horner(Ho,%i*Wo/Wo));
  
    // (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 = [];
    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
    V = coeff(N - D*prod(G))';
    
    // Find the coefficient vector
    AB = inv(CC)*V;
        
    // Coefficients of the independent term of the numerator 
    // of the partial fractions for the normalized function
    A = AB(1:2:$);
    // Coefficientes of the term in s of the numerator of the 
    // partial fractions for the normalized function
    B = AB(2:2:$);
    
    // Unnormalized coefficients of the numertors of the partial
    // fractions
    a = A*Wo^2;
    b = B*Wo;
        
    // Forced response
    y0 = delta(t) * prod(G);
    
    // Natural modes (unnormalized)
    wo = w.*sqrt(1 - 1/4./Q.^2);
    
    // Free response
    yy = zeros(n, length(t)); 
    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 = y0 + sum(yy, 1);
        
    if 1==2
        figure(100);
        gcf().background = 8;
        plot(t',[x'/Fs,yy',y'])
    end
    
endfunction

function x = delta(t);
    // Unit pulse
    // t: Array of reals
    // x: Array of reals whose components are
    //    0 if the corresponding component of
    //    t is nonzero, and 1 if they are 0  
    // May be used to find the impulse response
    // of a filter. This response can then be 
    // convolved with any discrete signal to
    // find the corresponding response of the
    // filter    
    x = zeros(t);
    x(find(t==0)) = 1; 
endfunction
