function [y, y0, yy, x] = band_tone_resp(t, w0, phi0, w, Q, G)
    // Pure tone analytic response of a general band-pass filter  
    // of order n
    //
    // Usage: 
    //        [y, y0, yy, x] = band_tone_resp(t, w0, phi0, w, Q, G)
    //
    // Where
    //        t:    Time vector containing the instants in s
    //              at which the response is to be computed 
    //        w0:   Angular frequency of tone in 1/s. May be 
    //              a row vector 
    //        phi0: Phase of tone in rad. May be a row cector
    //              the same size as w0
    //        w:    Vector containing the center frequencies
    //              of the n second-order sections in 1/s
    //        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.
    //              May be an array with as many responses as 
    //              frequencies in w0
    //        y0:   Column vector containing the forced response.
    //              May be an array with as many responses as 
    //              frequencies in w0
    //        yy:   Matrix with n columns, each containing the 
    //              components of the natural response. May 
    //              be an array with as many responses as the 
    //              length of w times the number of frequencies 
    //              in w0
    //        x:    Column vector containing the input signal.  
    //              May be an array with as many signals as 
    //              frequencies in w0
    //
    // NOTES: 1) Normally vectors w and Q are obtained from 
    //           an analog filter designed, for instance, 
    //           using the function analpf() and a LP --> BP
    //           conversion
    //        2) The n second-order factors have unit
    //           coefficients for the quadratic term. 
    //        3) The input signal is considered of unit
    //           amplitude.
    //        4) The center frequency is computed as the
    //           geometric mean of the center frequencies
    //           of the second-order sections
    //        5) x, y0, and y are in general arrays whose
    //           first dimension (columns) is the same 
    //           length as t and whose last dimensions 
    //           match those of w0
    //        6) yy is in general an array whose first 
    //           dimension (columns) is the same length as
    //           t, whose second dimension is equal to the
    //           order n of the filter and whose last 
    //           dimensions match those of w0
    //        7) If the order n is high, numerical errors 
    //           are expected since the algorithm requieres
    //           the computation of the inverse of a square
    //           matrix of order 2*n
    //       
    // Some theory
    // The input signal is
    //       x(t) = A*sin(w0*t) + B*cos(w0*t)
    // where 
    //       A = cos(phi0)
    //       B = sin(phi0)
    // Its Laplace transform is
    //       X(s) = (A*w0 + B*s) / (w0^2 + s^2)
    // 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) = 
    //              (a0 + b0*s)/(w0^2 + s^2) + 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
    // frequency is approximately 1, which makes the
    // magnitude order of the matrix coefficients more 
    // homogeneous
    //  
    // ------------------------------------------------
    // Author: Federico Miyara
    // Date:   2020-02-27
    //         2020-08-31
    //         2021-07-04
    //         2021-08-25
    
    // ------------------------------------------------
    // Input argument handling
    // ------------------------------------------------
    
    rhs = argn(2);
    if rhs<6
       G = 1;
    end
    if rhs<5
       msg = gettext("%s: Wrong number of input arguments: %d or %d expected.\n")
       error(msprintf(msg, "band_tone_resp", 5,6))        
    end
    if size(phi0)~=size(w0)
       msg = gettext("%s: Arguments #%d and #%d: Same sizes expected.\n")
       error(msprintf(msg, "resp_tono_banda", 2, 3))
    end
    if ~isreal(w0)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "band_tone_resp", 2))
    elseif ~and(w0>0)
       msg = gettext("%s: Argument #%d: Must be > %d.\n")
       error(msprintf(msg, "band_tone_resp", 2, 0))
    end
    if ~isreal(phi0)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "resp_tono_banda", 2))
    end
    if length(w)~=length(Q)
       msg = gettext("%s: Arguments #%d and #%d: Same sizes expected.\n")
       error(msprintf(msg, "band_tone_resp", 4,5))
    end
    if ~isreal(w)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "band_tone_resp", 4))
    elseif ~and(w>0)
       msg = gettext("%s: Argument #%d: Must be > %d.\n")
       error(msprintf(msg, "band_tone_resp", 4, 0))
    end
    if ~isreal(Q)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "band_tone_resp", 5))
    elseif ~and(Q>0)
       msg = gettext("%s: Argument #%d: Must be > %d.\n")
       error(msprintf(msg, "band_tone_resp", 5, 0))
    end
    if ~isreal(G)
       msg = gettext("%s: Argument #%d: Decimal number expected.\n")
       error(msprintf(msg, "band_tone_resp", 6))
    end
            
    // Test data --DELETE--
    if 1==2
       // Centre frequency of band filter
       wc = 2*%pi*1000
       // Frequencies and Q's of the second-order sections
       // of filter
       w = [1.357783939141098 1 0.73649202923462]*wc
       Q = [2.961756551889217 1.414213562373096 2.961756551889228] 
       // Filter gain
       G = 1
       // Sample rate
       Fs = 44100
       // Time vector
       t = (0:1/Fs:0.2)';
       // Frequency and phase of input sine wave
       cas = 1
       select cas
       case 1
          w0 = 2*%pi*[400, 500; 800, 1000]
          phi0 = [0, 0; 0 0]
       case 2
          w0 = 2*%pi*[400, 500]
          phi0 = [0, 0]
       case 3
          w0 = 2*%pi*[400]
          phi0 = [0]
       end
    end
        
    // Ensure that t is a column vector
    t = t(:);
    
    
    // ------------------------------------------------
    // Compute response of the filter
    // ------------------------------------------------
    
    // Order of the filter (actually, of its skirts since
    // as a dynamic system it has order 2*n)
    n = length(w);

    // Estimated center frequency of filter
    Wo = prod(w)^(1/n);
        
    // Sine and cosine oefficients of the input signal
    A = cos(phi0);
    B = sin(phi0);
    
    // Normalized Laplace variable
    p = poly(0,"p");
    
    // Normalized denominator of the input signal transform
    //P0 = w0^2/Wo^2 + p^2;
    P0 = w0(:).^2/Wo^2 + p^2;
    
    // 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 inverse matrix.  
    
    // Normalized transfer function of the filter
    H = p^n/prod(PH)/Wo^n;
    
    // Normalization constant to have unit gain at the 
    // normalized center frequency
    K = 1/abs(horner(H,%i*Wo/Wo));
  
    // Factors of the denominator of the transform of the
    // output
    P = [P0, kron(PH, ones(P0))];
    
    // (n-1)-factor products 
    for k=1:size(P,2)
        R(:,k) = prod([P(:,1:k-1), P(:,k+1:$)], 2);
    end
    
    // Coefficient array of the (n-1)-factor products.
    // For each page (corresponding to a frequency in w0)
    // each row contains the coefficients of each polynomial
    // obtained as the product of n-1 second-order factors
    // in ascending power orden 
    for h=1:length(w0)
       C(:,:,h) = coeff(R(h,:)');
    end
   
    // System matrix to obtain the coefficients of the partial 
    // fraction decomposition. Each page corresponds to a
    // frequency of the input signal
    for h = 1:length(w0)
       CCC = [];
       for k = 1:size(P,2)
          CCC = [CCC, [C(k,:,h),0]', [0, C(k,:,h)]'];
       end
       CC(:,:,h) = CCC;
    end
    
    // Independent vector corresponding to the coefficients
    // of the numerator of the response presented as column
    // vector. All frequencies must be unnormalized
    for h = 1:length(w0)
       V(:,h) = [zeros(PH), G*K*A(h)*w0(h)/Wo^(n+2), ...
                 G*K*B(h)/Wo^(n+1), zeros(PH)]';
    end
    
    // Calculation of the coefficient vector
    for h = 1:length(w0)
       AB(:,h) = inv(CC(:,:,h))*V(:,h);
    end
        
    // Coefficients of the independent term of the numerator 
    // of the partial fractions. Each row corresponds to a
    // frequency
    a = AB(1:2:$,:)*Wo^2;
    
    // Coefficientes of the term in s of the numerator of 
    // the partial fractions
    b = AB(2:2:$,:)*Wo;
    
    // Forced response
    // The frequency has been denormalized by using the
    // original frequency w0 of the signal. a(1) must be
    // denormalized multiplying by Wo since it has been 
    // computed using frequencies normalized by dividing
    // by Wo
    y0 = kron(a(1,:)./w0(:)', ones(t)) .* ...
            sin(kron(w0(:)', t) + kron(phi0(:)', ones(t))) + ...
         kron(b(1,:),ones(t)) .* ...
            cos(kron(w0(:)', t) + kron(phi0(:)', ones(t)));
    
    // In order to denormalize the response we use the
    // wo obtained from the orginal w
    wo = w.*sqrt(1 - 1/4./Q.^2);
    
    // Natural modes
    yy = zeros(length(t), n, length(w0));
    for k=1:n
        yy(:,k,:) = ...
           kron(exp(-w(k)/2/Q(k)*t), ones(w0(:)')) .* ...
           (kron(ones(t), b(k+1,:)) .* ...
           kron(cos(wo(k)*t), ones(w0(:)')) + ...
           kron(ones(t), a(k+1,:)/wo(k)-b(k+1,:)/2/Q(k)*w(k)/wo(k)).* ...
           kron(sin(wo(k)*t), ones(w0(:)')));
    end
             
    // Complete response
    y = sum([matrix(y0, size(y0,1),1,size(y0,2)), yy], 2);
    
    // Input signal
    x = sin(kron(w0(:)', t) + kron(phi0(:)', ones(t)));    


    // ------------------------------------------------
    // Output arguments formatting
    // ------------------------------------------------
    
    // Size conversion of signals to match the array
    // structure of w0. The last dimensions are those
    // of w0
    y0 = matrix(y0, [size(y0,1), size(w0)]); 
    yy = matrix(yy, [size(yy,1), size(yy,2), size(w0)]); 
    y = matrix(y, [size(y,1), size(w0)]);
    x = matrix(x, [size(x,1), size(w0)]);    


    // For testing purposes only --DELETE--
    if 1==2
       for k=1:length(w0)
          figure(1000+k);
          clf(1000+k);
          gcf().background = -2
          plot(t,[x(:,k),y0(:,k),y(:,k)])
       end

       for k=1:length(w0)
          figure(1000+length(w0)+k);
          clf(1000+length(w0)+k);
          gcf().background = -2
          plot(t,yy(:,:,k))
       end
    end
    
endfunction

