// stepresponse.sce
//
// This script calculates the step response of a specified bass reflex system
// design based on the Weideman Contour method as described in the AES article:
// "A Contour Integral Method for Transient Response Calculation"
// J. Audio Eng. Soc., vol. 66, no. 5, pp. 360-368, (May 2018).
//
// The transducer data are from the following AES article (not the datasheet):
// "An Added-Mass Measurement Technique for Transducer Parameter Estimation"
// (because it includes visco-elastic parameters, classic Knudsen LOG model)
// and therefore this script is using the Knudsen LOG model for simulation.
//
// To "verify" the time domain response, an FFT is performed, which is compared
// to the exact solution (= the pre-defined frequency response function).
// Due to this extra step, the Weideman calculation of the time-response is
// "stupidly" long, but can of course be reduced significantly (as sketched in
// the AES paper) if the FFT is not needed.
//
// Script by Claus Futtrup and Jeff Candy
// Date created: 2018-08-05
// Last modified: 2018-10-30

clear; xdel(winsid()); clc; // Clear variables, graphics and console

// ====================
// ---- USER INPUT ----
// ====================

name = "SEAS L16RNX 8 ohm in 10 liter box"; // Included in PLOT title
saveplot = %t; // Save the plots? Yes = %t and No = %f

// --- DEFINE TRANSDUCER PARAMETERS
Re = 5.8; // ohm - DC resistance - world's simplest "inductance model"
Bl = 7.034; // Tm (or N/A)
Mms = 15.13; // gram
Sd = 104; // cm2
// Instead of the classical Thiele-Small parameters of Cms, Rms,
// we need three parameters, named R0, L0 and beta to describe the
// viscoelastic damping and compliance of the mechanical system.
R0 = 32.20; // ohm (damping) - for gbeta = 0, insert Res (T/S parameter)
L0 = 60.24; // mH (compliance) - for gbeta = 0, insert Lces (T/S parameter)
gbeta = 0.059; // Viscoelasticity, if gbeta = 0, then without creep (i.e. T/S)
// Above three parameters describe damping (R0) and compliance (L0, gbeta).
// We recommend that they are found using the dual-added-mass technique.
// Notice the zero ('0') in the names of the parameters, which show that these
// values are valid following Knudsen LOG model and _not_ classical T/S.
// The following model equations can degenerate to the Thiele-Small model just
// by setting gbeta = 0, R0 = Res = Bl^2 / Rms and L0 = Lces = Cms * Bl^2;

// --- DEFINE BOX PARAMETERS
Vb = 10; // liter
Fp = 40; // Hz - port tuning frequency, Fp = 0 for closed box simulation

// ===========================
// ---- END OF USER INPUT ----
// ===========================

// --- CONVERT TO SI-UNITS
Mms = Mms / 1e3; // convert to kg
Sd = Sd / 1e4; // convert to m2
L0 = L0 / 1e3; // convert to Henry
Vb = Vb / 1e3; // convert to m3

// --- SOME CONSTANTS
rho = 1.201; // Density of air, kg/m3
c = 343.7; // Speed of sound in air, m/s
p0 = 2e-5; // 20 micro-Pascal, Reference pressure level (SPL = 0 dB)
r = 1; // Distance 1 meter
Vin = sqrt(8); // Input voltage 2.83 Volt (step response signal input level)

// --- CALCULATE FURTHER PARAMETERS
Cm0 = L0 / Bl^2; // Compliance (mechanical)
ws = 1/sqrt(Cm0*Mms);
F0 = ws / (2 * %pi);
  // This is _not_ the resonance frequency of the complete speaker, but
  // only the frequency-independent part (= excluding viscoelasticity).
  // We do it this way because the viscoelastic compliance and damping
  // is added as a separate expression to the response function.
Rm0 = Bl^2 / R0; // Damping (mechanical), frequency-independent part
Qm0 = (Mms/Rm0) * ws;
Rme = Bl^2 / Re; // Mechanical equivalent of electrical (DC) resistance
Qe0 = (Mms/Rme) * ws;
Qt0 = (1/Qm0 + 1/Qe0)^-1;
Vas = Cm0 * Sd^2 * rho * c^2; // Volume of Air Equivalent to Cm0
glambda = gbeta * log(10) / (1-gbeta * log(ws)); // Knudsen LOG lambda
                                                 // Not used, just curious

// --- NORMALIZE DATA
alpha = Vas / Vb;
h = Fp / F0;
// NOTICE: When using the LOG model, these normalizations
// are in physical size quite different from Thiele/Small.

// Extract the filename and path of this script:
[units, typs, nams] = file(); 
nams(find(units==0 | units==5 | units==6)) = [];
scriptpath = fileparts(nams($),"path");
chdir(scriptpath);

// --- RESPONSE FUNCTION

function y=resp(s)
    cs = 1 - gbeta * log( s ); // Knudsen LOG viscoelastic compliance
    y = s.^4 ./ ( (s.^2 + h^2) .* (1.0./cs + s/Qt0 + s.^2) + alpha * s.^2);
endfunction
// NOTICE: This defines a 4th order system (bass reflex), utilizing the Knudsen
// LOG model with a viscoelastic compliance (i.e. the output is complex and
// gives both stiffness and damping to the final response). The parameters, h,
// Qt0, alpha and gbeta (greek-beta) are defined with a normalized frequency
// given by the LOG model, i.e. not the resonance frequency of the complete
// speaker. The only difference from the original Knudsen LOG article is that
// we use the natural logarithm (not base-10 logarithm).

// --- WEIDEMAN CONTOUR INTEGRAL

N0 = 12;     // resolution
steps = 8;   // steps per unit time,
t_max = 16;  // time (dimensionless)
N0 = 32;
steps = 128; // settings for length(t) = 4096 and MAX precision
t_max = 32;
t_step = 1/steps;
t = t_step:t_step:t_max;

mprintf("Phew! - Give me a minute.\n");

mu_c = max([1 h]);
t_c = %pi*N0/(12*mu_c);
nt = length(t);

for i=1:nt do
    ti = t(i);       // Selection defined in
    if ti < t_c then // Eqs (31) and (32)
        mu = %pi*N0/(12*ti);
        N = N0;
        d = 3/N;
    else
        mu = mu_c;
        N = ceil(N0*ti/t_c);
        d = 3/N;
    end

    xsk = 0.0; // Perform sum in Eq (24)
    for k=-N:N do
        u = k*d;
        s = mu * (%i*u + 1)^2;
        dsdu = 2.0*mu*(%i-u);
        xsk = xsk + imag(exp(s*ti) * resp(s)/s * dsdu);
    end
    xsk = xsk * d / (2*%pi);
    xs(i) = xsk; // Collect x-values (the normalized pressure response)
end

// --- UNDO NORMALIZATION
tau = t / ws; // time, in seconds (x-axis)
p1 = rho * Sd / (2 * %pi * r) * Vin * Bl / (Mms * Re); // Peak pressure (Pascal)
p = p1*xs; // pressure respones, in Pascal (y-axis)
tau = [0 tau]; // add startpoint (row vector)  -  data set is now 4097 long
p = [p1; p];   // add startpoint (column vector)

// --- PLOT STEP RESPONSE
g_tr = scf();
plot(tau.',p); // Transpose X to prevent WARNING
xgrid(color("grey70"));
title_string = "Step (air pressure) time response, " + name;
title(title_string);
xlabel("Time (seconds)");
ylabel("Pressure (Pascal)");
if saveplot then xs2png(gcf(),"step.png"); end

// --- FFT OF STEP RESPONSE

// Do an FFT of the step response, as a fact-check on above step-resp. calc.
// We do not calculate and plot the impulse response due to discontinuity issues
tau_step = t_step / ws; // not-normalized, convert to time in seconds
f = ws * linspace(0,1/t_step,nt); // Define f (not normalized), from 0 - 30 kHz
w = 2 * %pi * f;
iw = %i * w; // like s, but 's' variable name is already used (in Weideman loop)
g_w = tau_step * fft(p1*xs);
    // Above FFT is performed on the REAL (not normalized) pressure response.
g_mag = abs(iw.' .* g_w); // Transpose iw to prevent warning
g_phs = 180 / %pi * atan(imag(iw.' .* g_w) , real(iw.' .* g_w));
for i=1:length(g_phs) // unwrap phase, continuity
   if g_phs(i) < 0.0 then
      g_phs(i) = g_phs(i)+360.0 ;
   end
end

SPLw = 20 * log10(g_mag / p0 + 1e-20);
// Convert to deci-Bell - Note: g_mag(1) = 0, this equation returns -Inf for the
// first data point, so it's treated by adding a tiny 1e-20 to the expression

// --- FREQUENCY RESPONSE (the simple, normalized way)
fn = 0.1:0.01:2; // This is needed only for frequency response plot, normalized
wn = 2 * %pi * fn; // Omega, normalized
iwn = %i * wn; // like s, but 's' variable name is already used

Y = resp(iwn); // Use the pre-defined response function to plot the freq-resp.
MAG = 20 * log10( abs(Y) );
PHS = 180 / %pi * atan(imag(Y) , real(Y));

for i=1:length(Y) // unwrap phase, for continuity
   if PHS(i) < 0.0 then
      PHS(i) = PHS(i)+360.0 ;


   end
end

// Voltage SPL output in dB:
USPL = 20 * log10(Vin * rho * Sd * Bl / (2 * %pi * r * p0 * Re * Mms));

scf();
subplot(2,1,1);
plot(F0*wn,MAG + USPL,'-b',f(3:$/3),SPLw(3:$/3).','-r');
fm = gca(); // Get current axis
fm.log_flags = "lnn"; // Set x-axis = log
xgrid(color("grey70"));
title_string = "Magnitude Frequency response, " + name;
title(title_string);
xlabel("Frequency (Hz)");
ylabel("Magnitude (dB)");
legend("Direct","FFT",pos=4);

subplot(2,1,2);
plot(F0*wn,PHS,'-b',f(3:$/3),g_phs(3:$/3).','-r');
fp = gca(); // Get current axis
fp.log_flags = "lnn"; // Set x-axis = log
xgrid(color("grey70"));
title_string = "Phase Frequency response, " + name;
title(title_string);
xlabel("Frequency (Hz)");
ylabel("Phase (deg)");

if saveplot then xs2png(gcf(),"freqresp.png"); end
