diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..baed105 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -57,6 +57,9 @@ Modules: * pytta.iso3741: calculations according to the standard; + * pytta.loudspeaker: + Thiele-Small parameters from the electrical impedance; + .. toctree:: classes/index utils/index @@ -64,4 +67,5 @@ Modules: functions generate rooms + loudspeaker properties diff --git a/docs/pyttadoc/loudspeaker.rst b/docs/pyttadoc/loudspeaker.rst new file mode 100644 index 0000000..e5f497b --- /dev/null +++ b/docs/pyttadoc/loudspeaker.rst @@ -0,0 +1,7 @@ +.. currentmodule:: pytta + +Loudspeaker +=========== + +.. automodule:: pytta.loudspeaker + :members: diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..9c15212 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -72,6 +72,9 @@ * pytta.iso3741: calculations according to the standard; + * pytta.loudspeaker: + Thiele-Small parameters from the electrical impedance; + You can find usage examples in the toolbox's examples folder. @@ -108,6 +111,7 @@ from . import utils from . import rooms from . import iso3741 +from . import loudspeaker from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +134,7 @@ 'generate', 'utils', 'iso3741', + 'loudspeaker', 'rooms', # Functions diff --git a/pytta/loudspeaker.py b/pytta/loudspeaker.py new file mode 100644 index 0000000..092a9f7 --- /dev/null +++ b/pytta/loudspeaker.py @@ -0,0 +1,279 @@ +""" +PyTTa loudspeaker: +------------------ + +Small signal (Thiele-Small) parameters of an electrodynamic loudspeaker from +its measured electrical impedance, with the lumped model of Thiele (1971) +and Small (1972): + + Z = Re + j w Le + Bl^2 / (Rms + j w Mms + 1 / (j w Cms)). + +The free air curve gives Re, fs, Qms, Qes and Qts. A second curve, with a +known mass added to the diaphragm or with the driver on a closed box of +known volume, gives Mms, Cms, Rms, Bl and, with the effective diaphragm +area, Vas. + +Available functions: + + >>> pytta.loudspeaker.thiele_small(impedance) + >>> pytta.loudspeaker.impedance_model(params, freq) + +References: + + A. N. Thiele, "Loudspeakers in vented boxes: Parts I and II", + J. Audio Eng. Soc. 19 (1971) 382-392 and 471-483. + + R. H. Small, "Direct-radiator loudspeaker system analysis", + J. Audio Eng. Soc. 20 (1972) 383-395. + + R. H. Small, "Closed-box loudspeaker systems, Part I: Analysis", + J. Audio Eng. Soc. 20 (1972) 798-808. + + L. L. Beranek and T. J. Mellow, "Acoustics: Sound Fields and + Transducers", Academic Press, 2012, section 6.10. + +For further information, check the function specific documentation. + +""" + +import numpy as np +from scipy.optimize import least_squares + +from pytta.classes import SignalObj + + +def impedance_model(params, freq): + """ + Electrical impedance of the lumped loudspeaker model. + + Z = Re + j 2 pi f Le + Re (Qms / Qes) / (1 + j Qms (f / fs - fs / f)) + + The last term is the motional impedance Bl^2 / (Rms + j w Mms + + 1 / (j w Cms)) written with fs = 1 / (2 pi sqrt(Mms Cms)), + Qms = 2 pi fs Mms / Rms and Qes = 2 pi fs Mms Re / Bl^2. At f = 0 it is + zero and Z = Re. + + Parameters + ---------- + params : dict + 'Re' (ohm), 'fs' (Hz), 'Qms', 'Qes' and optionally 'Le' (H, + default 0). The dict returned by thiele_small fits here. + freq : array_like + Frequencies in Hz. + + Returns + ------- + numpy.ndarray + Complex impedance in ohm, same shape as freq. + + """ + f = np.asarray(freq, dtype=float) + Re, fs, Qms = params['Re'], params['fs'], params['Qms'] + res = Re * Qms / params['Qes'] + zmot = res * f * fs / (f * fs + 1j * Qms * (f**2 - fs**2)) + return Re + 2j * np.pi * f * params.get('Le', 0.0) + zmot + + +def thiele_small(impedance, Re=None, area=None, added_mass=None, + impedance_with_mass=None, added_volume=None, + impedance_in_box=None, fit=False, c=343.2, rho=1.204): + """ + Thiele-Small parameters from the electrical impedance of a loudspeaker. + + Impedance method (Small 1972; Beranek and Mellow 2012, eqs. 6.58 and + 6.59): fs and Rmax at the impedance peak, found by a parabola through + its three bins (the peak is the local maximum of |Z| next to the + positive to negative sign change of the reactance, so the rise of |Z| + due to Le at high frequencies is skipped), and f1 < fs < f2 where |Z| = sqrt(Rmax Re), + interpolated linearly between bins. With r0 = Rmax / Re, + + Qms = fs sqrt(r0) / (f2 - f1), Qes = Qms / (r0 - 1), + Qts = Qms Qes / (Qms + Qes). + + Added mass dm (eq. 6.66), with fs' the resonance of the loaded curve: + + Mms = dm / ((fs / fs')^2 - 1), Cms = 1 / ((2 pi fs)^2 Mms). + + Added volume Vb (closed box, eq. 6.72), with fc and Qec from the curve + in the box: + + Vas = Vb (fc Qec / (fs Qes) - 1), Cms = Vas / (rho c^2 Sd^2). + + Then Rms = 2 pi fs Mms / Qms, Bl = sqrt(2 pi fs Mms Re / Qes) and + Vas = rho c^2 Sd^2 Cms. + + The closed form neglects the voice coil inductance; with fit=True every + curve is also fitted by least squares to impedance_model over the whole + band (Re, Le, fs, Qms, Qes; Re only when it is not given), starting + from the closed form. The relative complex error is minimized. + + Differences from ita_thiele_small (ITA-Toolbox): Re defaults to the + minimum of the real part below the resonance instead of the real part + at the lowest frequency; the peak and the side frequencies are interpolated instead + of read at whole bins; the bandwidth is taken on |Z| at sqrt(Rmax Re) + instead of on |Z - Re| at (Rmax - Re) / sqrt(2), which gives the same + Qms when Le = 0. + + Parameters + ---------- + impedance : SignalObj or tuple (freq, Z) + Free air impedance in ohm, complex. A SignalObj must have one + channel in freqSignal; bins at f <= 0 are ignored. + Re : float, optional + DC resistance of the voice coil in ohm, as measured with an + ohmmeter. Default: minimum of the real part of the impedance + below the resonance (in the model the real part is Re plus a + positive term that vanishes at low frequencies). + area : float, optional + Effective diaphragm area Sd in m^2. Needed for Vas (added mass) or + for the mechanical parameters (added volume). + added_mass : float, optional + Mass in kg added to the diaphragm for impedance_with_mass. + impedance_with_mass : SignalObj or tuple (freq, Z), optional + Impedance with added_mass on the diaphragm. + added_volume : float, optional + Net volume in m^3 of the closed box for impedance_in_box. + impedance_in_box : SignalObj or tuple (freq, Z), optional + Impedance of the driver on the closed box of added_volume. + fit : bool, optional + Refine the closed form by least squares and return Le. Default + False. + c : float, optional + Speed of sound in m/s. Default 343.2 (air at 20 degC). + rho : float, optional + Density of air in kg/m^3. Default 1.204 (air at 20 degC). + + Returns + ------- + dict + 'Re' (ohm), 'Le' (H, only with fit), 'fs' (Hz), 'Qms', 'Qes', + 'Qts'; with a second curve also 'Rms' (kg/s), 'Mms' (kg), + 'Cms' (m/N) and 'Bl' (T m), and 'Vas' (m^3) when it can be found. + With added_volume and no area only 'Vas' is added. + + Raises + ------ + ValueError + If the reactance does not turn from positive to negative inside + the band, the peak lies at the band edge, |Z| does not fall to + sqrt(Rmax Re) on both sides of it, if a second curve comes without + its mass or volume (or both methods are given), or if the added + mass does not lower the resonance or the box does not raise it. + + """ + if (added_mass is None) != (impedance_with_mass is None): + raise ValueError("added_mass and impedance_with_mass go together.") + if (added_volume is None) != (impedance_in_box is None): + raise ValueError("added_volume and impedance_in_box go together.") + if added_mass is not None and added_volume is not None: + raise ValueError("Give either the added mass or the added volume.") + for name, value in (('area', area), ('added_mass', added_mass), + ('added_volume', added_volume), ('Re', Re)): + if value is not None and not value > 0: + raise ValueError(f"{name} must be positive.") + + freq, z = _curve(impedance) + ts = _estimate(freq, z, Re, fit) + Re = ts['Re'] + Qms, Qes = ts['Qms'], ts['Qes'] + ws = 2 * np.pi * ts['fs'] + out = dict(ts, Qts=Qms * Qes / (Qms + Qes)) + + Mms = Cms = None + if added_mass is not None: + loaded = _estimate(*_curve(impedance_with_mass), Re, fit) + if not loaded['fs'] < ts['fs']: + raise ValueError("The added mass must lower the resonance.") + Mms = added_mass / ((ts['fs'] / loaded['fs'])**2 - 1) + Cms = 1 / (ws**2 * Mms) + elif added_volume is not None: + boxed = _estimate(*_curve(impedance_in_box), Re, fit) + if not boxed['fs'] > ts['fs']: + raise ValueError("The closed box must raise the resonance.") + out['Vas'] = added_volume * (boxed['fs'] * boxed['Qes'] + / (ts['fs'] * Qes) - 1) + if area is not None: + Cms = out['Vas'] / (rho * c**2 * area**2) + Mms = 1 / (ws**2 * Cms) + + if Mms is not None: + out.update(Rms=ws * Mms / Qms, Mms=Mms, Cms=Cms, + Bl=np.sqrt(ws * Mms * Re / Qes)) + if area is not None: + out['Vas'] = rho * c**2 * area**2 * Cms + return out + + +def _curve(impedance): + """Frequencies above 0 Hz and complex impedance from either input.""" + if isinstance(impedance, SignalObj): + if impedance.numChannels != 1: + raise ValueError("The impedance SignalObj must have one channel.") + freq, z = impedance.freqVector, impedance.freqSignal[:, 0] + else: + freq, z = impedance + freq = np.asarray(freq, dtype=float) + z = np.asarray(z, dtype=complex) + if freq.ndim != 1 or freq.shape != z.shape: + raise ValueError("freq and Z must be 1-D arrays of the same length.") + keep = freq > 0 + freq, z = freq[keep], z[keep] + if freq.size < 3 or np.any(np.diff(freq) <= 0): + raise ValueError("freq must increase and hold at least 3 bins > 0 Hz.") + return freq, z + + +def _estimate(freq, z, Re, fit): + """Re, fs, Qms, Qes (and Le with fit) of one curve.""" + mag = np.abs(z) + # Resonance: the + to - sign change of the reactance with the largest + # |Z|, then uphill to the local maximum of |Z|. + down = np.flatnonzero((z.imag[:-1] > 0) & (z.imag[1:] <= 0)) + if down.size == 0: + raise ValueError("The reactance never turns from positive to " + "negative: no resonance inside the band.") + down = down + (mag[down + 1] > mag[down]) + k = down[np.argmax(mag[down])] + while 0 < k < mag.size - 1 and mag[k + 1] > mag[k]: + k += 1 + while 0 < k < mag.size - 1 and mag[k - 1] > mag[k]: + k -= 1 + if k == 0 or k == mag.size - 1: + raise ValueError("The impedance peak lies at the band edge.") + fit_re = Re is None + if Re is None: + Re = np.min(z.real[:k]) + a, b, c = np.polyfit(freq[k-1:k+2] - freq[k], mag[k-1:k+2], 2) + fs = freq[k] - b / (2 * a) + zmax = c - b**2 / (4 * a) + level = np.sqrt(zmax * Re) + below = np.flatnonzero(mag[:k] < level) + above = k + np.flatnonzero(mag[k:] < level) + if below.size == 0 or above.size == 0: + raise ValueError("|Z| does not fall to sqrt(Rmax Re) on both sides " + "of the peak.") + i, j = below[-1], above[0] + f1 = np.interp(level, mag[i:i+2], freq[i:i+2]) + f2 = np.interp(level, mag[j-1:j+1][::-1], freq[j-1:j+1][::-1]) + r0 = zmax / Re + Qms = fs * np.sqrt(r0) / (f2 - f1) + ts = dict(Re=Re, fs=fs, Qms=Qms, Qes=Qms / (r0 - 1)) + if fit: + ts = _fit(freq, z, ts, fit_re) + return ts + + +def _fit(freq, z, ts, fit_re): + """Least squares fit of impedance_model, Le included.""" + names = ['Re'] * fit_re + ['fs', 'Qms', 'Qes', 'Le'] + le0 = max(z[-1].imag / (2 * np.pi * freq[-1]), 0.0) + x0 = [ts[n] for n in names[:-1]] + [le0] + + def residual(x): + err = (impedance_model(dict(ts, **dict(zip(names, x))), freq) - z) + err /= np.abs(z) + return np.concatenate([err.real, err.imag]) + + sol = least_squares(residual, x0, bounds=(0, np.inf), x_scale='jac', + xtol=1e-12, ftol=1e-12, gtol=1e-12) + return dict(ts, **dict(zip(names, sol.x))) diff --git a/tests/golden/make_thiele_small_ita.m b/tests/golden/make_thiele_small_ita.m new file mode 100644 index 0000000..0d458cf --- /dev/null +++ b/tests/golden/make_thiele_small_ita.m @@ -0,0 +1,54 @@ +function make_thiele_small_ita(itaRoot) +% MAKE_THIELE_SMALL_ITA Golden data for tests/test_loudspeaker.py. +% +% make_thiele_small_ita(itaRoot) builds the electrical impedance of a +% loudspeaker with known lumped parameters with ita_linear_loudspeaker +% from the ITA-Toolbox (RWTH Aachen, https://git.rwth-aachen.de/ita/toolbox), +% once in free air and once with an added mass, runs ita_thiele_small on +% the pair and writes thiele_small_ita.mat next to this file. +% +% The file holds the input parameters, the frequency grid (sampling rate +% and FFT degree), the free air impedance at every 20th bin up to 400 Hz +% (to check the model) and the parameters estimated by ita_thiele_small. +% The voice coil inductance is zero: ita_thiele_small reads the +% impedance peak and bandwidth without an inductance correction. +% +% Usage: matlab -batch "make_thiele_small_ita('/path/to/ITA-Toolbox')" + +addpath(genpath(itaRoot)); + +sr = 44100; +fftDegree = 18; +Re = 6.0; % Ohm +Mms = 0.012; % kg +Cms = 7.0e-4; % m/N +Rms = 1.2; % kg/s +Bl = 6.0; % T m +dm = 0.010; % kg, added mass +d = 0.13; % m, membrane diameter + +TS = struct('R_e', Re, 'L_e', 0, 'm', Mms, 'n', Cms, 'w', Rms, 'M', Bl); +Z = ita_linear_loudspeaker(TS, 'samplingRate', sr, 'fftDegree', fftDegree); +TSm = TS; +TSm.m = Mms + dm; +Zm = ita_linear_loudspeaker(TSm, 'samplingRate', sr, 'fftDegree', fftDegree); + +est = ita_thiele_small(Z, Zm, dm, d, 'frequency_limits', [1 7000]); + +idx = (1:20:Z.freq2index(400))'; +freqSample = Z.freqVector(idx); +zSample = Z.freqData(idx); +c = double(est.c); +rho0 = double(est.rho0); + +ita = struct('fs', double(est.f_s), 'Re', double(est.R_e), ... + 'Qms', double(est.Q_m), 'Qes', double(est.Q_e), 'Qts', double(est.Q_tot), ... + 'Mms', double(est.m), 'Cms', double(est.n), 'Rms', double(est.w), ... + 'Bl', double(est.M), 'Vas', double(est.V_equi)); + +out = fullfile(fileparts(mfilename('fullpath')), 'thiele_small_ita.mat'); +save(out, 'sr', 'fftDegree', 'Re', 'Mms', 'Cms', 'Rms', 'Bl', 'dm', 'd', ... + 'c', 'rho0', 'freqSample', 'zSample', 'ita', '-v7'); +fprintf('wrote %s\n', out); +disp(ita); +end diff --git a/tests/golden/thiele_small_ita.mat b/tests/golden/thiele_small_ita.mat new file mode 100644 index 0000000..70a9f0e Binary files /dev/null and b/tests/golden/thiele_small_ita.mat differ diff --git a/tests/test_loudspeaker.py b/tests/test_loudspeaker.py new file mode 100644 index 0000000..4380221 --- /dev/null +++ b/tests/test_loudspeaker.py @@ -0,0 +1,184 @@ +""" +Tests for pytta.loudspeaker. + +Each test pins one rule. The impedance curves are synthesized from known +parameters, so every estimate has an exact target; the only error left is +the interpolation between frequency bins. The last but one test compares +with the ITA-Toolbox (RWTH Aachen) on data generated by +golden/make_thiele_small_ita.m. + +The lumped model is the one of Thiele (1971) and Small (1972): + + Z = Re + j w Le + Bl^2 / (Rms + j w Mms + 1 / (j w Cms)), + +and its electrical form, the one impedance_model takes, + + Z = Re + j w Le + Re (Qms / Qes) / (1 + j Qms (f / fs - fs / f)), + +with fs = 1 / (2 pi sqrt(Mms Cms)), Qms = 2 pi fs Mms / Rms and +Qes = 2 pi fs Mms Re / Bl^2. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +import pytta +from pytta import loudspeaker + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden', + 'thiele_small_ita.mat') +FREQ = np.arange(1, 4001) * 0.05 # 0.05 Hz to 200 Hz, 0.05 Hz step +FREQ = np.concatenate([FREQ, np.arange(201, 20001, 1.0)]) +C, RHO = 343.2, 1.204 # the module defaults +TOL = 1e-5 # linear interpolation between 0.05 Hz bins + + +def mechanical(f, Re, Mms, Cms, Rms, Bl, Le=0.0): + """Electrical impedance of the mechanical circuit, written directly.""" + jw = 2j * np.pi * f + return Re + jw * Le + Bl**2 / (Rms + jw * Mms + 1 / (jw * Cms)) + + +def rel(a, b): + return abs(a / b - 1) + + +class TestThieleSmall(unittest.TestCase): + + def test_impedance_method_recovers_the_small_signal_parameters(self): + # fs, Qms and Qes from the peak and the two frequencies where + # |Z| = sqrt(Rmax Re) (Small 1972; Beranek and Mellow 2012, 6.58 and + # 6.59). With Le = 0 the method is exact, and Re is the minimum of + # the real part below fs, reached at the lowest frequency. + true = dict(Re=6.2, fs=48.0, Qms=4.5, Qes=0.45) + z = loudspeaker.impedance_model(true, FREQ) + est = loudspeaker.thiele_small((FREQ, z)) + true['Qts'] = true['Qms'] * true['Qes'] / (true['Qms'] + true['Qes']) + print() + for key in true: + print(f' {key}: true {true[key]:.6g}, estimated {est[key]:.6g}, ' + f'rel diff {rel(est[key], true[key]):.1e}') + self.assertLess(rel(est[key], true[key]), TOL, key) + for key in ('Le', 'Mms', 'Cms', 'Rms', 'Bl', 'Vas'): + self.assertNotIn(key, est) + + def test_added_mass_gives_the_mechanical_parameters(self): + # Mms = dm / ((fs / fs')^2 - 1), Cms = 1 / ((2 pi fs)^2 Mms), + # Rms = 2 pi fs Mms / Qms, Bl^2 = 2 pi fs Mms Re / Qes and + # Vas = rho c^2 Sd^2 Cms (Beranek and Mellow 2012, 6.66 and 6.27). + # Both curves come from the mechanical circuit, so this also pins + # the electrical form against it. + true = dict(Re=5.6, Mms=0.0215, Cms=9.0e-4, Rms=1.6, Bl=7.3) + dm, area = 0.015, 0.0201 + z = mechanical(FREQ, **true) + zm = mechanical(FREQ, **dict(true, Mms=true['Mms'] + dm)) + est = loudspeaker.thiele_small((FREQ, z), Re=true['Re'], area=area, + added_mass=dm, + impedance_with_mass=(FREQ, zm)) + true['Vas'] = RHO * C**2 * area**2 * true['Cms'] + print() + for key in ('Mms', 'Cms', 'Rms', 'Bl', 'Vas'): + print(f' {key}: true {true[key]:.6g}, estimated {est[key]:.6g}, ' + f'rel diff {rel(est[key], true[key]):.1e}') + self.assertLess(rel(est[key], true[key]), TOL, key) + + def test_added_volume_gives_vas_and_the_mechanical_parameters(self): + # In a closed box of volume Vb the compliance becomes Cms Cmb / + # (Cms + Cmb), with Cmb = Vb / (rho c^2 Sd^2), and + # Vas = Vb (fc Qec / (fs Qes) - 1) (Small 1972; Beranek and Mellow + # 2012, 6.69 and 6.72). With Sd the rest follows as above. + true = dict(Re=7.1, Mms=0.0098, Cms=1.1e-3, Rms=0.9, Bl=5.2) + vb, area = 0.008, 0.0133 + cmb = vb / (RHO * C**2 * area**2) + z = mechanical(FREQ, **true) + zb = mechanical(FREQ, **dict(true, Cms=true['Cms'] * cmb + / (true['Cms'] + cmb))) + est = loudspeaker.thiele_small((FREQ, z), Re=true['Re'], area=area, + added_volume=vb, + impedance_in_box=(FREQ, zb)) + true['Vas'] = RHO * C**2 * area**2 * true['Cms'] + print() + for key in ('Vas', 'Mms', 'Cms', 'Rms', 'Bl'): + print(f' {key}: true {true[key]:.6g}, estimated {est[key]:.6g}, ' + f'rel diff {rel(est[key], true[key]):.1e}') + self.assertLess(rel(est[key], true[key]), TOL, key) + + def test_fit_recovers_the_voice_coil_inductance(self): + # A voice coil inductance tilts the peak and moves the sqrt(Rmax Re) + # points, so the closed form is biased. The least squares fit of the + # model on the whole curve returns Le and removes that bias. At + # 20 kHz |Z| = 150 ohm, above the resonance peak of 67 ohm, so the + # peak search must skip the inductive rise. + true = dict(Re=3.4, Le=1.2e-3, fs=31.0, Qms=6.0, Qes=0.32) + z = loudspeaker.impedance_model(true, FREQ) + closed = loudspeaker.thiele_small((FREQ, z)) + est = loudspeaker.thiele_small((FREQ, z), fit=True) + print() + for key in true: + bias = rel(closed[key], true[key]) if key in closed else np.nan + print(f' {key}: true {true[key]:.6g}, fitted {est[key]:.6g}, ' + f'rel diff {rel(est[key], true[key]):.1e} ' + f'(closed form {bias:.1e})') + self.assertLess(rel(est[key], true[key]), 1e-6, key) + + def test_parity_with_ita_thiele_small(self): + # ita_linear_loudspeaker builds the curves on an FFT grid of + # 44100 / 2^18 = 0.168 Hz and ita_thiele_small reads them bin by bin + # (resonance at the bin where the phase changes sign, bandwidth at + # the first and last bins above Rp / sqrt(2)), so its estimates are + # off by up to a few bins. The model must match to rounding; the + # estimates within 1 %. + g = sio.loadmat(GOLDEN, squeeze_me=True) + Re, Mms, Cms, Rms, Bl = (float(g[k]) for k in + ('Re', 'Mms', 'Cms', 'Rms', 'Bl')) + ws = 1 / np.sqrt(Mms * Cms) + params = dict(Re=Re, fs=ws / (2 * np.pi), Qms=ws * Mms / Rms, + Qes=ws * Mms * Re / Bl**2) + Qts = params['Qms'] * params['Qes'] / (params['Qms'] + params['Qes']) + model = loudspeaker.impedance_model(params, g['freqSample']) + err = np.max(np.abs(model - g['zSample']) / np.abs(g['zSample'])) + print(f'\n impedance model vs ita_linear_loudspeaker: {err:.1e}') + self.assertLess(err, 1e-12) + + sr, n = int(g['sr']), 2 ** int(g['fftDegree']) + f = np.fft.rfftfreq(n, 1 / sr) + dm, area = float(g['dm']), np.pi * float(g['d'])**2 / 4 + z = mechanical(f[1:], Re, Mms, Cms, Rms, Bl) + zm = mechanical(f[1:], Re, Mms + dm, Cms, Rms, Bl) + z = pytta.SignalObj(np.r_[Re, z], 'freq', samplingRate=sr, + signalType='energy') + zm = pytta.SignalObj(np.r_[Re, zm], 'freq', samplingRate=sr, + signalType='energy') + est = loudspeaker.thiele_small(z, area=area, added_mass=dm, + impedance_with_mass=zm, + c=float(g['c']), rho=float(g['rho0'])) + ita = g['ita'] + true = dict(params, Mms=Mms, Cms=Cms, Rms=Rms, Bl=Bl, Qts=Qts, + Vas=float(g['rho0']) * float(g['c'])**2 * area**2 * Cms) + print() + for key in ('fs', 'Re', 'Qms', 'Qes', 'Qts', 'Mms', 'Cms', 'Rms', + 'Bl', 'Vas'): + d = rel(est[key], float(ita[key])) + t = rel(est[key], true[key]) if key in true else np.nan + print(f' {key}: ITA {float(ita[key]):.6g}, PyTTa {est[key]:.6g}, ' + f'rel diff {d:.1e} (PyTTa vs true {t:.1e})') + self.assertLess(d, 1e-2, key) + + def test_rejects_invalid_input(self): + true = dict(Re=6.2, fs=48.0, Qms=4.5, Qes=0.45) + z = loudspeaker.impedance_model(true, FREQ) + above = FREQ > 60 # the peak falls at the lower band edge + with self.assertRaises(ValueError): + loudspeaker.thiele_small((FREQ[above], z[above])) + with self.assertRaises(ValueError): # curve without its mass + loudspeaker.thiele_small((FREQ, z), impedance_with_mass=(FREQ, z)) + higher = loudspeaker.impedance_model(dict(true, fs=60.0), FREQ) + with self.assertRaises(ValueError): # added mass raised fs + loudspeaker.thiele_small((FREQ, z), added_mass=0.01, + impedance_with_mass=(FREQ, higher)) + + +if __name__ == '__main__': + unittest.main()