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..e62cc75 --- /dev/null +++ b/pytta/loudspeaker.py @@ -0,0 +1,455 @@ +""" +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) + >>> pytta.loudspeaker.response(params, freq, area) + >>> pytta.loudspeaker.max_spl(params, freq, area, xmax, power_max) + >>> pytta.loudspeaker.radiation_impedance(freq, radius) + >>> pytta.loudspeaker.closed_box(params, Vb) + >>> pytta.loudspeaker.horn_profile(throat_area, cutoff, length) + >>> pytta.loudspeaker.fit_with_velocity(freq, impedance, velocity) + >>> pytta.loudspeaker.fit_voice_coil(freq, blocked, 'L2R2') + +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))) + + +def radiation_impedance(freq, radius, c=343.2, rho=1.204): + """ + Mechanical radiation impedance in N s/m on one side of a rigid piston of + the given radius (m) in an infinite baffle, S rho c (1 - J1(2ka) / (ka) + + j H1(2ka) / (ka)), with J1 and H1 the Bessel and Struve functions + (Beranek and Mellow 2012, eqs. 13.117 and 13.118). + """ + from scipy.special import j1, struve + x = 2 * 2 * np.pi * np.asarray(freq, dtype=float) / c * radius + with np.errstate(invalid='ignore', divide='ignore'): + z = 1 - 2 * j1(x) / x + 2j * struve(1, x) / x + z = np.where(x == 0, 0, z) + return np.pi * radius**2 * rho * c * z + + +def response(params, freq, area, voltage=2.83, distance=1., c=343.2, + rho=1.204): + """ + Small signal response of a driver in an infinite baffle to an rms + voltage, from its Thiele-Small parameters ('Re', 'Bl', 'Mms', 'Cms', + 'Rms' and optionally 'Le', as thiele_small returns them with a second + curve) and its effective area Sd (m^2), as ita_thielesmall2velocity. + + The diaphragm velocity is v = Bl U / ((Re + j w Le) Zm + Bl^2) with + Zm = Rms + j w Mms + 1 / (j w Cms) (Mms includes the air load), the + pressure on axis at distance r in the far field of the baffled piston is + j w rho Sd v e^(-jkr) / (2 pi r), and the radiated power is + |v|^2 Re(Zr) with Zr from radiation_impedance. + + Returns a dict with 'velocity' (m/s), 'displacement' (m), 'current' + (A), 'pressure' (Pa), all rms and complex, 'SPL' (dB re 20 uPa), + 'power' (W) and 'efficiency', the power over voltage^2 / Re (Beranek + and Mellow 2012, eq. 6.47 in the mid band). + """ + f = np.asarray(freq, dtype=float) + jw = 2j * np.pi * f + Ze = params['Re'] + jw * params.get('Le', 0.0) + Zm = params['Rms'] + jw * params['Mms'] + 1 / (jw * params['Cms']) + Bl = params['Bl'] + v = Bl * voltage / (Ze * Zm + Bl**2) + p = jw * rho * area * v * np.exp(-jw / c * distance) / (2 * np.pi * distance) + power = np.abs(v)**2 * radiation_impedance(f, np.sqrt(area / np.pi), c, rho).real + return {'velocity': v, 'displacement': v / jw, + 'current': voltage * Zm / (Ze * Zm + Bl**2), 'pressure': p, + 'SPL': 20 * np.log10(np.abs(p) / 2e-5), 'power': power, + 'efficiency': power / (voltage**2 / params['Re'])} + + +def max_spl(params, freq, area, xmax, power_max, distance=1., c=343.2, + rho=1.204): + """ + Largest SPL at each frequency limited by the peak excursion xmax (m) and + by the electrical power power_max (W) in Re, as + ita_loudspeakertools_maxSPL. Returns (SPL, rms voltage, limit), limit + being 'excursion' or 'thermal' at each frequency. + """ + one = response(params, freq, area, 1., distance, c, rho) + u_x = xmax / (np.sqrt(2) * np.abs(one['displacement'])) + u_t = np.sqrt(power_max / (np.abs(one['current'])**2 * params['Re'])) + u = np.minimum(u_x, u_t) + return one['SPL'] + 20 * np.log10(u), u, np.where(u_x < u_t, 'excursion', 'thermal') + + +def closed_box(params, Vb): + """ + Resonance frequency and total Q of the driver in a closed box of net + volume Vb (m^3), fc = fs sqrt(1 + Vas / Vb) and Qtc = Qts sqrt(1 + + Vas / Vb), without box losses (Small 1972, closed box; Beranek and + Mellow 2012, eq. 7.28). + """ + k = np.sqrt(1 + params['Vas'] / Vb) + return params['fs'] * k, params['Qts'] * k + + +def horn_profile(throat_area, cutoff, length, a=1., c=343.2, num=200): + """ + Cross-sectional area along a hyperbolic-exponential (hypex) horn, + S = St (cosh(x / xT) + a sinh(x / xT))^2 with xT = c / (2 pi fc), from + hyperbolic (a = 0) to exponential (a = 1, S = St e^(mx), m = 4 pi fc / + c) (Beranek and Mellow 2012, eqs. 9.39, 9.51 and 9.58). Returns + (x, S) at num points from the throat to length (m). + """ + x = np.linspace(0, length, num) + xT = c / (2 * np.pi * cutoff) + return x, throat_area * (np.cosh(x / xT) + a * np.sinh(x / xT))**2 + + +def fit_with_velocity(freq, impedance, velocity, weights=None): + """ + Thiele-Small parameters from the electrical impedance Z = U / i and the + diaphragm velocity per volt H = v / U measured together (e.g. with a + laser vibrometer), as ita_thiele_small_fit and, with weights per + frequency, ita_thiele_small_fit_weighted. With G = v / i = H Z, the + model gives two linear least squares problems: + + Z = Re + j w Le + Bl G, + Bl / G = Rms + j w Mms + 1 / (j w Cms). + + Returns a dict with 'Re', 'Le', 'Bl', 'Rms', 'Mms', 'Cms', 'fs', 'Qms', + 'Qes' and 'Qts'. + """ + w = 2 * np.pi * np.asarray(freq, dtype=float) + Z = np.asarray(impedance) + G = np.asarray(velocity) * Z + sw = np.sqrt(np.ones(w.size) if weights is None else np.asarray(weights)) + + def solve(columns, target): + M = np.column_stack(columns) * sw[:, None] + rhs = target * sw + return np.linalg.lstsq(np.vstack([M.real, M.imag]), + np.r_[rhs.real, rhs.imag], rcond=None)[0] + + Re, Le, Bl = solve([np.ones_like(w), 1j * w, G], Z) + Rms, Mms, K = solve([np.ones_like(w), 1j * w, 1 / (1j * w)], Bl / G) + ws = np.sqrt(K / Mms) + Qms, Qes = ws * Mms / Rms, ws * Mms * Re / Bl**2 + return {'Re': Re, 'Le': Le, 'Bl': Bl, 'Rms': Rms, 'Mms': Mms, + 'Cms': 1 / K, 'fs': ws / (2 * np.pi), 'Qms': Qms, 'Qes': Qes, + 'Qts': Qms * Qes / (Qms + Qes)} + + +def voice_coil_model(params, freq): + """ + Blocked electrical impedance of the voice coil: Re + K (j w)^n when + params has 'K' and 'n' (Leach, Loudspeaker voice-coil inductance + losses, J. Audio Eng. Soc. 50(6), 2002), otherwise Re + j w Le plus, + with 'L2' and 'R2', the eddy current branch j w L2 R2 / (R2 + j w L2) + (the L2R2 model). + """ + jw = 2j * np.pi * np.asarray(freq, dtype=float) + if 'K' in params: + return params['Re'] + params['K'] * jw ** params['n'] + z = params['Re'] + jw * params.get('Le', 0.) + if 'L2' in params: + z = z + jw * params['L2'] * params['R2'] / (params['R2'] + jw * params['L2']) + return z + + +def fit_voice_coil(freq, blocked, model='L2R2', Re=None): + """ + Parameters of voice_coil_model fitted by least squares on the relative + complex error to a blocked impedance (for instance Z - Bl G from + fit_with_velocity, or a measurement with the diaphragm held), as + ita_inductance_fit and ita_curvefit_inductance. model is 'Le', 'L2R2' + or 'Leach'; Re is fitted unless given. + """ + f = np.asarray(freq, dtype=float) + z = np.asarray(blocked) + names = {'Le': ['Le'], 'L2R2': ['Le', 'L2', 'R2'], 'Leach': ['K', 'n']}[model] + w = 2 * np.pi * f + Le0 = np.median(z.imag / w) + x0 = {'Le': Le0, 'L2': Le0, 'R2': z.real.max(), 'K': Le0, 'n': 0.9} + fit_re = Re is None + start = [x0[n] for n in names] + ([z.real.min()] if fit_re else []) + + def params(x): + p = dict(zip(names, x)) + p['Re'] = x[-1] if fit_re else Re + return p + + def residual(x): + err = voice_coil_model(params(x), f) / z - 1 + return np.concatenate([err.real, err.imag]) + + sol = least_squares(residual, start, bounds=(0, np.inf), x_scale='jac', + xtol=1e-14, ftol=1e-14, gtol=1e-14) + return params(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..9c884d1 --- /dev/null +++ b/tests/test_loudspeaker.py @@ -0,0 +1,258 @@ +""" +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() + + +DRIVER = {'Re': 6.0, 'Le': 0.0, 'Bl': 7.5, 'Mms': 0.02, 'Cms': 8e-4, + 'Rms': 1.5, 'fs': 1 / (2 * np.pi * np.sqrt(0.02 * 8e-4))} +SD = 0.022 + + +class TestLoudspeakerTools(unittest.TestCase): + + def test_reference_efficiency(self): + # In the mid band, where the mass dominates the mechanical + # impedance and ka is small (here 0.15), the efficiency is + # rho Bl^2 Sd^2 / (2 pi c Re Mms^2) (Beranek and Mellow 2012, + # eq. 6.47). + small = dict(DRIVER, Bl=3., Mms=0.05) + sd = np.pi * 0.02**2 + eta = loudspeaker.response(small, [400.], sd)['efficiency'][0] + expected = 1.204 * 3.**2 * sd**2 / (2 * np.pi * 343.2 * 6.0 * 0.05**2) + print(f'\n efficiency {eta:.4e} for {expected:.4e}', end='') + self.assertLess(abs(eta / expected - 1), 0.01) + + def test_radiation_impedance_limits(self): + # Low frequencies: R = S rho c (ka)^2 / 2 and X = S rho c 8 ka / + # (3 pi); high frequencies: R -> S rho c. + a, c, rho = 0.1, 343.2, 1.204 + S = np.pi * a**2 + ka = 0.01 + z = loudspeaker.radiation_impedance(ka * c / (2 * np.pi * a), a) / (S * rho * c) + self.assertLess(abs(z.real / (ka**2 / 2) - 1), 1e-4) + self.assertLess(abs(z.imag / (8 * ka / (3 * np.pi)) - 1), 1e-4) + high = loudspeaker.radiation_impedance(200 * c / (2 * np.pi * a), a) / (S * rho * c) + self.assertLess(abs(high.real - 1), 0.01) + + def test_max_spl_meets_both_limits(self): + f = np.array([20., 40., 1000.]) + spl, u, limit = loudspeaker.max_spl(DRIVER, f, SD, 0.004, 50.) + r = loudspeaker.response(DRIVER, f, SD, 1.) + x = np.sqrt(2) * np.abs(r['displacement']) * u + power = np.abs(r['current'] * u)**2 * DRIVER['Re'] + self.assertEqual(list(limit), ['excursion', 'excursion', 'thermal']) + np.testing.assert_allclose(x[:2], 0.004) + np.testing.assert_allclose(power[2], 50.) + + def test_closed_box_and_horn(self): + fc, Qtc = loudspeaker.closed_box({'fs': 30., 'Qts': 0.35, 'Vas': 0.06}, 0.02) + self.assertAlmostEqual(fc, 60.) + self.assertAlmostEqual(Qtc, 0.7) + # a = 1 is the exponential horn S = St e^(mx), m = 4 pi fc / c. + x, S = loudspeaker.horn_profile(0.001, 300., 0.8) + np.testing.assert_allclose(S, 0.001 * np.exp(4 * np.pi * 300. / 343.2 * x)) + + def test_fit_with_velocity_recovers_the_driver(self): + # Impedance and velocity of the lumped model with Le are fitted + # back exactly. + drv = dict(DRIVER, Le=0.5e-3) + f = np.geomspace(10, 2000, 200) + jw = 2j * np.pi * f + Ze = drv['Re'] + jw * drv['Le'] + Zm = drv['Rms'] + jw * drv['Mms'] + 1 / (jw * drv['Cms']) + Z = Ze + drv['Bl']**2 / Zm + H = drv['Bl'] / (Ze * Zm + drv['Bl']**2) + fit = loudspeaker.fit_with_velocity(f, Z, H) + e = max(abs(fit[k] / drv[k] - 1) for k in ('Re', 'Le', 'Bl', 'Rms', 'Mms', 'Cms', 'fs')) + print(f'\n fit with velocity: {e:.1e}', end='') + self.assertLess(e, 1e-9) + + def test_voice_coil_models_are_fitted_back(self): + f = np.geomspace(100, 20000, 120) + for model, true in (('L2R2', {'Re': 6., 'Le': 0.3e-3, 'L2': 0.6e-3, 'R2': 4.}), + ('Leach', {'Re': 6., 'K': 0.02, 'n': 0.7})): + fit = loudspeaker.fit_voice_coil(f, loudspeaker.voice_coil_model(true, f), model) + e = max(abs(fit[k] / true[k] - 1) for k in true) + print(f'\n voice coil {model}: {e:.1e}', end='') + self.assertLess(e, 1e-6)