From 7e25f5977a5efa3b16bdd0d6606cf3815554888e Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 21:56:32 +0100 Subject: [PATCH] [FEATURE] add pytta.loudspeaker with Thiele-Small parameters New module pytta.loudspeaker estimates the small signal parameters of an electrodynamic loudspeaker from its electrical impedance, given as a SignalObj or as (freq, Z) arrays. thiele_small reads fs, Re, Qms, Qes and Qts from the impedance peak and the two frequencies where |Z| = sqrt(Rmax Re) (Small 1972; Beranek and Mellow 2012, section 6.10). A second curve with an added mass or in a closed box of known volume gives Mms, Cms, Rms, Bl and Vas. With fit=True every curve is refined by a least squares fit of the lumped model, which also returns Le and removes the bias that Le causes in the closed form. impedance_model returns Z = Re + j w Le + Zmot of the same model. Six tests: recovery of the parameters from synthetic curves for the impedance, added mass, added volume and fit methods, parity with ita_linear_loudspeaker and ita_thiele_small of the ITA-Toolbox (golden data and generator script in tests/golden), and rejection of invalid input. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/index.rst | 4 + docs/pyttadoc/loudspeaker.rst | 7 + pytta/__init__.py | 5 + pytta/loudspeaker.py | 279 +++++++++++++++++++++++++++ tests/golden/make_thiele_small_ita.m | 54 ++++++ tests/golden/thiele_small_ita.mat | Bin 0 -> 3281 bytes tests/test_loudspeaker.py | 184 ++++++++++++++++++ 7 files changed, 533 insertions(+) create mode 100644 docs/pyttadoc/loudspeaker.rst create mode 100644 pytta/loudspeaker.py create mode 100644 tests/golden/make_thiele_small_ita.m create mode 100644 tests/golden/thiele_small_ita.mat create mode 100644 tests/test_loudspeaker.py 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 0000000000000000000000000000000000000000..70a9f0ee2585ee6350bc58ad8e6ea078103d8dba GIT binary patch literal 3281 zcma);c{CJy8^;G@&lZKkj4?5mv9BXjlP&w$vop4Ytc9^_SS)R*YrFn|sf_=brPPd*Amt=lPxI{Lc4tp65BwKVKu1g(2z^SV>wIY=p9adwBV{ zOM%foF2Nptfxd`ef0`&|MJcdmpu0=3yBpZ=Is$C&5)3wR4F)U7f#u~8N^%H!Ww5;L zujl^_z@I^*XTEYw`Kqr9UW`>O`DE3001;e zd+}Hf|4Z)wcR9qx7;J1H8K7q*qi>W@?n>LsMB@c;BLM}EA?RPwfAv61(51TpD`)^` zDrwFgYbF2ER{y<0yPQIt-*q`UF^fZM2NM2%kKhhKYnXxWe;`4BAml2+S3aQWH5$zi z^oz#`>>s;mqJKxG`A#6G=m07ez*GOeL7X}JpBo@(1paj#ZcTk@x?dYK+qxAUbd;YC z2dHWSH;)mCe;P+blO_FAe6oyq^>R}S`IZJimIQo$Onjj^KAR;@&(PeIF9C#NcCn_!)14Ac z6}BP@?AvbF1TU;_M{S3(SI0XwazOK-Wso>W|GCnHbuFB$i5i8YNPc5$b!RV@nsLMk z6+d;jtH6uKMs?;5$_EET=c%15HnHY@$I4Sx0u7?;&*+D4({;Z5&JzSQ@)Yenj|$il z>$GPh(-U4tpcW1>vYyKz^^tdHc@j9hKC$PgZYypNwp}186w+P&6z!&(CKb`$A+Yv^7UrP`aJ829N$Q9R8##Fu~>VW`DiaQABv<+^pxb$bt2>Iot3)Ggzk1nFNS% zN?&|BnZTV-PCIYPgadIpmI_SPD#&{{@oAfdwirD)lm+hB@;cQSK<3$Gt~HNdLgPl2 zyIHh~WkZ(-bFo0BSCgV=Dsqd*5A5n)Ik?{?hXnLqC$^0{2X(q4&W_=Q7Fe{@(v^P1 z&jw!+OcZdKN+BEqQMAj^9GTviye1uA^YsA~D`#xwvkbx44er|c-9Dwon6F`Z?L|Tt zEIgrX7ch_`l2*5D7u&shuel&wRMr)CG6&j4AD&FA5_#SPPK=n^4}R|a)s-#BJ-U*% zqY*AL_2EpLZ;enY$bs4DzZ4-W+V;xNp%mRlBiNSi?fT_VeX! zLt=4gTuqR@Bx^RKzxs2cP4QQOS5sFRS)-yBN3eq_^n1LBb%bQ~e3tpO^?WZ;i=$1D z(F2s|3E9A@ExWwxm=}sd=+5}Vnl}ecYKD9+b~eo!~^MAQi%!t zQ$5#2rhj9sH&e6Xh*{PuXRbD*9@I&vdWU5t)`LFQ?=s|uws`Dj@)LiU7JBdX%)yoo zkXMIvbS>lz%6a=9bEA(iVO=eCIsm91co?BlChROpLF7gApi+`IqxRf@r>EE23+^!e z=GN_-)+OA4L+8|~pN>f;4A17mfGN?~5FQf&}Qk!u&{Jrbv7e07T>u|ulEAex2R3%(3Fs9iN;tKG-ClQ@|3uyV zDQ|jvZHf|x;byCQHj0Or54&htyd(u^v8SF@%H=0XwwLm6Ns}(B4*Jc2UVM5IlhYcH zt1M%H!q^o_k|)NTqifjGig5aC+)qdxe=qaRhCr82M)7eT|{^)jkJOS zMG}A9)|#~6udtyzSu`X>?zAM>Q#yTN^>#^PTK13SHpU@i4}}IuLl zP7H)I)b_^jQFONEV<8_LMb1idnnD@IJ5-`f?NFG-0ufe?Q7)ElG1fWTg@UUE0RRz%}EcNNjG=h=H)0UY^XX zyZKwPiynucHH&yEdq%6#BTg}{pOZc1=cHErq29w+nRhF)>(@JJr+ukcFZLMRGAT2} zM=4PSL#$gsu*;Oxy9RhqZ({4>*joMT1G^TNf;)TaEi`i*cj^=pd}FoOT*r9}DN_(q z>Z3HLK)W~kWYC4)1Sv+w?B@GBH5>iNPQ0u z5NNTT#Y+h(q?%O_6yh&LK##$#1`4`M!pd-;5NW5>GV?4Z&p&*h<0Tf{@4_PZ{*@u2 zqa67oGAlLc%64Ao*PD|9cnLzj%IAdzyVKY9I~vcwV{iB}x6sjnncSnq`3QL|)rHA+ zc?5V=FC|s)+eCnk6F158QakKS@amAs)0q>>g4O+p z_s_?%h0_vN8~f*Lwx@a>EIc<>)y4cKjE!AaW(9j|1xLw;^B*3pZjX{%FFQE!{tuHT ze{9kQBM(}{Dyqh404m$`vB$(edRbOYQ=Jx?UtTtgXNvTQf$)t+hn@cD7-{;If7tQR zwl{5U3}tG!_}&8PtY8GL)I5H9gOf2Ir7kjYjRHQ_VHgsNEf^BprZI_a7GP;{j}7|T zq8L@C+eLwVeC;Ri@{?~}gY$9jR4>#q(yHR#$asdUr&nC?#zjsT2tA}Z cr+lQz>Lb+;3o`xq%%c6Pzi9wM(-fxv1A(^RWB>pF literal 0 HcmV?d00001 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()