// Testing filter_imp_resp()
//
// NOTE: The response of a single complex pole is
//       complex. The responses of a pair of conjugate
//       poles are conjugate, there sum being real.
//       However, when plotting each response, the
//       imaginary part is ignored, so both responses
//       seem to be equal. Their sum (twice either of
//       them) is the correct response to the complex
//       pair.
//
// ---------------------------------------------------
// Federico Miyara
// 2021-07-20
// 2025-10-08
    

caso = 2
select caso
case 1
    // Poles
    pp = [-4:-1]*10
    // Zeros
    zz = [0 0] 
    // Constant
    K = 1
case 2
    // Low-pass prototype
    [hs, pols, zers, gain] = analpf(3, "butt", [], 1)
    // Relative bandwidth for one third octave band
    B = 2^(1/6) - 2^(-1/6);
    wo = 2*%pi*1000;
    H = horner(hs, (wo/%s + %s/wo)/B);
    pp = roots(H.den)'
    zz = roots(H.num)'
    // Constant required to express numerator
    // and denominator as monic polynomials
    K = coeff(H.num)($)/coeff(H.den)($)
case 3
    // Low-pass prototype
    [hs, pols, zers, gain] = analpf(3, "butt", [], 1)
    // Relative bandwidth for octave band
    B = 2^(1/2) - 2^(-1/2);
    // Center angular frequency
    wo = 2*%pi*1000;
    // Conversion to band-pass
    H = horner(hs, (wo/%s + %s/wo)/B);
    // Poles and zeros
    pp = roots(H.den)'
    zz = roots(H.num)'
    // Constant required to express numerator
    // and denominator as monic polynomials
    K = coeff(H.num)($)/coeff(H.den)($)
end

// Sample rate
Fs = 44100;
// Time vector
t = 0:1/Fs:1;
// Input signal (a Dirac's delta with ideal 
// antialiasing filter with cutoff frequency
// at Fs/2)
x = [Fs, zeros(t)(1:$-1)];
// Frequency vector
ww = 2*logspace(1,4,1000)*2*%pi;    

// Load function
exec filter_imp_resp.sci;

// Impulse response
[y, y0, yy] = filter_imp_resp(t, pp, zz, K);

// Maximum time constant
// Used to plot the most relevant part of the signal
taumax = max(-1./real(pp))

figure(200);
scf(200);
clf();
gcf().background = -2;
subplot(2,1,1)
plot(t, y)
gca().data_bounds = [0, -max(abs(y)); 5*taumax, max(abs(y))];
gca().grid = [1 1];
xlabel("t  [s]")
ylabel("y") 
subplot(2,1,2)
plot(t', yy') 
gca().data_bounds = [0,-max(abs(y)); 5*taumax,max(abs(y))];
gca().grid = [1 1];
xlabel("t  [s]")
ylabel("yy") 

figure(201);
scf(201);
clf();
gcf().background = -2;
plot("ln", ww/2/%pi, 20*log10(abs(horner(H,%i*ww))))
gca().grid = [1 1];
xlabel("f  [Hz]")
ylabel("|H|   [dB]") 
xstring(100, 100,"")
