fix: SII disturbance per ANSI S3.5, IEC-exact intermodulation and the AES17 measurement chain (#173)
* fix: SII equivalent disturbance is the clause 5.6 maximum, not an energy sum
ANSI S3.5-1997 clause 5.6 defines the equivalent disturbance spectrum
level Di as the larger of the equivalent masking level Zi and the
equivalent internal noise level Xi'. The module summed the two in energy,
underestimating the SII by up to 0.042 in noise-plus-hearing-loss
conditions (0.1841 vs 0.2185 for normal speech in flat 30 dB noise with a
flat 40 dB loss). The official Hornsby worksheet computes the maximum in
every band and the R CRAN worked example confirms it.
Oracles: quiet-condition SII pinned at the worksheet full precision
(0.99582516666667, abs 1e-6), the disturbance intermediate locked at
5000 Hz (-23.6 dB), a discriminating noise-plus-loss case (0.2185) and
the independent R package Example C.2 (0.8513749). A code comment records
that the level-distortion factor referencing the normal-effort spectrum
for every vocal effort is the standard's own definition.
* test: STIPA direct-method modulation-depth staircase and NIPTS reading notes
The IEC 60268-16:2020 C.3.2 verification signal (Formula C.1, sinusoidal
carriers at the A.6.1 male levels) is synthesized end to end and the
published STI staircase for m = 0.0 to 1.0 (0, 0.18, 0.30, 0.38, 0.44,
0.50, 0.56, 0.62, 0.70, 0.82, 1) is asserted within the published
tolerance of 0.05 through the full stipa() audio path; the m = 0.5 point
also joins the conformance report.
ISO 1999 notes: the sub-10-year NIPTS spread keeps the raw lg(t) reading
(Formulae 6/7 carry no Formula-3 extrapolation provision) and the Annex D
half-integer rounding ambiguity of Table D.4 is recorded next to the
pinned cells.
* fix: IEC-exact intermodulation distortion (14.12.7, 14.12.8, 14.12.10)
modulation_distortion now returns the per-order IEC 60268-3 values: d_m,2
and d_m,3 are the arithmetic sums of the sideband amplitudes referenced to
the output voltage at f2 (14.12.7.2 g-h), bundled with the SMPTE-analyzer
combined RMS as an explicitly labelled convention in a new
ModulationDistortionResult.
difference_frequency_distortion references U_2,ref = 2 U_2,f2 (the sum of
both tone amplitudes) and sums the third-order products arithmetically per
14.12.8.1; the former half-difference search window could swallow the
primary tones (octave-spaced clean tones read d2 = 1.0) or the DC bin, so
product windows are now a few bins wide, capped by the component spacing,
clamped to (0, Nyquist) and forced to zero when they cannot be separated
from a primary tone or DC.
total_difference_frequency_distortion implements 14.12.10.1 directly: only
the two in-band products at f0 -/+ delta, rms-summed over the arithmetic
sum of the tone output voltages, with the standard 8 kHz / 11.95 kHz tones
as defaults.
The conformance entries carry the corrected per-clause labels and exact
values (d_m,2 = 0.16, d_m,3 = 0.08, d_d,2 = 0.03, d_d,3 = 0.04, d_TDFD =
sqrt(0.0013)). The 14.12.9.2 f) print of the DIM denominator (U2 vs the
14.12.9.1 f_s amplitude) is recorded in docs/ERRATA.md with a code comment
at the reference measurement.
* fix: AES17 measurement bandwidth, effective notch Q and the ITU-R 468 weighted THD
thd_plus_noise and sinad now measure both the residual and the total
signal through the AES17 chain (5.2.5 / 6.3.1): a 20 Hz high-pass plus
the standard low-pass at a configurable bandwidth (default 20 kHz),
realised as a brick-wall frequency-domain band per 5.2.10. Out-of-band
noise at high sample rates (+7 dB at fs = 192 kHz) and DC offsets no
longer inflate the result; bandwidth=None keeps the full-Nyquist
measurement.
The standard notch is designed so the APPLIED zero-phase response has the
requested quality factor (AES17 5.2.8 defines Q on the combined response;
filtfilt squares the magnitude, so the single-pass design is sharpened by
sqrt(1 + sqrt(2))). Previously a nominal Q of 2.0 acted as 1.29 and
requests below ~1.87 silently fell outside the AES17 [1.2, 3] range.
weighted_thd defaults to the network IEC 60268-3 14.12.11 actually
requires: the IEC 60268-1 Appendix A curve, i.e. ITU-R BS.468-4 Table 1
(peak +12.2 dB at 6.3 kHz), interpolated per the recommendation's own
log-frequency rule and exposed as itu_r_468_weighting; A and C remain as
labelled options and the 31.5 Hz - 400 Hz fundamental validity range is
documented. THD now raises when no harmonic fits below Nyquist instead of
returning 0, captures shorter than 64 samples are rejected, and the THD/
SINAD docstrings state the 14.12.3.2 R-convention and the derived nature
of SINAD explicitly.
New oracles: clipped-sine Fourier constants, the full-signal DIM
regression recipe (1.536 MHz synthesis, weak polynomial nonlinearity,
FIR decimation to 192 kHz, exact-bin FFT oracle), the H1 input-noise
bias SNR/(1+SNR), the printed ITU-R 468 table values including the AES17
CCIR-RMS cross-check, and the effective-notch-Q validation.
* test: ISO 532-1 Annex B.4 level-ramp oracles and reference-source provenance
The Annex B.4 ramp signals 6-9 join the validation suite: the tone ramps
(250 Hz / 1 kHz / 4 kHz, 30 dB to 80 dB at 5 dB/s over 10 s inside a
10.6 s file) are synthesized from parameters measured on the official
electronic-attachment WAVs and pinned to the workbook Nmax/N5 with
tolerances calibrated against those WAVs (module-on-official-WAV agrees
with the workbook to < 0.01 % on Nmax; the synthesis lands within +0.3 %
/ -0.8 %, pinned at 0.6 % / 2 %). The pink-noise ramp (overall SPL
25.85 + 5t dB per the WAV's Annex B.1 calibration) is pinned at a
realization-honest 10 %.
The sub-1-sone phon conversion (0.0005 offset, 3.0 phon floor, exact
10/lg 2 factor) and the NX percentile ((k-1, k) mean at
k = int((1 - X/100) n)) are now confirmed line by line against the
attachment's Annex A.4 reference program source; the provenance
docstrings record the exact code they match.
* docs: IEC-exact IMD semantics, AES17 bandwidth and ITU-R 468 weighting across guides and reference
The electroacoustics guide (EN and ES) and the API reference describe the
per-order modulation and difference-frequency definitions with their
references and arithmetic sums, the 14.12.10 two-product TDFD with the
standard 8 kHz / 11.95 kHz defaults, the AES17 measurement-bandwidth chain
with its bandwidth escape hatch, the effective notch Q, the ITU-R BS.468-4
default weighting for weighted THD with the 31.5-400 Hz validity range,
the DIM denominator erratum link and the corrected THD/SINAD convention
wording. CHANGELOG records the added API and the fixes; the conformance
report and the distortion figure variants are regenerated.
* docs: the ISO 12354-1:2017 edition repeats the K24 clamp misprint and flips the double-leaf K24 sign
The 2017 edition repeats the impossible '0 <= K24 <= -4 dB' interval of the
2000 print verbatim, and its E.3.5 prints the double-leaf K24 with a plus
where the 2000 edition's formula and its own corroborating Figure E.9 curve
carry a minus; both are recorded in the errata registry and the sign choice
is noted at the implementation.
* docs: the theory pages and SII diagram teach the max-form disturbance
The equivalent disturbance prose in the theory pages and the diagram label
follow the corrected ANSI S3.5 clause 5.6 reading (the larger of masking
and internal noise, not their energy sum), in English and Spanish; the
ISO 532-1 provenance docstrings paraphrase the reference program instead of
quoting it; and the notch-Q compensation note states the bilinear-warp
caveat near Nyquist.
* fix: regenerate the SII figure for the max-form disturbance and analysis follow-ups
The speech-intelligibility figure consumes the corrected clause 5.6
disturbance, so its four variants are regenerated; the 468 weighting's
zero-frequency guard avoids floating-point equality and the exception
tests keep one throwing call per block.
feat: swept-sine harmonic separation, THD(f) by order and phase utilities (Farina / Novak) (#226)
* feat: swept-sine harmonic separation and phase utilities (Farina / Novak / Bendat & Piersol)
synchronized_sweep_signal() generates the synchronized exponential sweep
of Novak, Lotton & Simon (2015, Eqs. 47/49): the rounded rate L makes
f1*L an integer, so delaying the signal by L*ln(n) is exactly equivalent
to generating its nth harmonic (Eq. 18) and the harmonic phases measured
downstream are properties of the system under test, not of the
excitation.
swept_sine_distortion() in electroacoustics deconvolves a recorded sweep
response and windows the impulse response of every distortion order at
its exact L*ln(n) advance (integer part in the window position, the
fractional remainder removed with a frequency-domain phase ramp on the
one-sided spectrum, and the centring delay compensated). Two paths:
method="synchronized" (default) uses the closed-form analytic
inverse-filter spectrum X~(f) = 2*sqrt(f/L)*exp(-j2pifL(1-ln(f/f1)) +
jpi/4) (Novak Eqs. 50-51), which extends each Hn to [n*f1, n*f2];
method="farina" reuses the ISO 18233 inverse-filter machinery of
room_ir for recordings made with the classical ESS (magnitudes exact,
phases documented as excitation-dependent). The frozen
SweptSineDistortionResult carries H1..HN, the windowed harmonic IRs, the
arrival advances, THD(f) with per-order distortion ratios evaluated at
|Hn(n*f)|/|H1(f)| inside each order's valid band, and .plot(). The
recording mean is removed by default (Novak Sec. A.2.3) and the analysis
is referenced to the excitation amplitude. room_ir cross-references the
new analyzer for the distortion products impulse_response discards; the
electroacoustics -> room edge is whitelisted in the architecture test.
The metrology phase utilities implement the B&P Sec. 13.1.4 Hilbert
relation: minimum_phase() reconstructs the phase of a magnitude via the
real cepstrum, with exact trigonometric-interpolation oversampling
against cepstral aliasing and a documented floor for magnitude zeros;
group_delay() differentiates the unwrapped phase; excess_phase() returns
the all-pass part (measured minus minimum, DC-referenced); and
phase_decomposition() bundles everything into a frozen result with the
total and excess group delays and .plot().
Oracles: the Chebyshev identities of the memoryless cubic polynomial
(|H1| = 1+3a3/4, |H2| = a2/2 at -pi/2, |H3| = a3/4 at pi, closed-form
THD, all within a few 1e-3), the time-shift/harmonic identity of the
synchronized sweep at 1e-9, agreement with the tone-by-tone thd()
analyzer at 1e-3, the pure-linear THD floor below 1e-3, a delayed/gained
Hammerstein composition, the band extension beyond f2, a strictly
minimum-phase biquad reconstructed to 1e-12 rad, the exact first-order
allpass group delay (1-a^2)/(1+2a cos w+a^2), and the exact -wD excess
phase and constant excess group delay of a fractional latency.
* feat(conformance): swept-sine distortion and phase-utility checks
New domain with seven anchors: the Chebyshev |H3| = a3/4 magnitude and
the pi phase of the third harmonic separated from a synchronized sweep,
the closed-form THD(1 kHz) of the cubic polynomial, the THD floor of a
purely linear path, the minimum-phase reconstruction of a strictly
minimum-phase biquad at the 1e-9 level, the exact first-order-allpass
group delay at w = pi/2, and the excess group delay of a biquad delayed
7.25 samples. 306/306 checks pass; the report regenerates
byte-identically.
* feat(figures): swept-sine THD by order (Farina / Novak)
One figure for the concept the separator adds: the total THD(f) and the
per-order ratios d2(f)/d3(f) of a cubic polynomial followed by a 3 kHz
low-pass, measured from one synchronized sweep. Each order sits exactly
on its Chebyshev level at low frequency and rolls off where its own
product n*f crosses the filter corner. Four language/theme variants; ES
strings added to the translation tables. The remaining committed figures
regenerate within the check_figures tolerance.
* docs: swept-sine distortion and phase utilities guide (EN/ES), API reference and indices
New guide in the three documentation trees: the L*ln(n) harmonic
separation with the THD(f) definition, the synchronized sweep with its
analytic deconvolution and band extension, the practice notes for
classical ESS recordings, and the phase utilities with their numerical
contract and sampling precautions. Wired into the Sources and devices
section pages and the sidebar. Adds the new modules to the API taxonomy
(electroacoustics / spectral analysis), regenerates the Starlight API
reference, extends the curated API table, the README highlights and
guide table, the bibliography (Farina, Novak et al., Mueller &
Massarani), the llms files (regenerated from the current docs) and the
changelog.
* fix: review round for the swept-sine separator and its docs
- The Farina buffer is a linear convolution whose anticausal span is only
sweep_samples - 1: an order arriving beyond it (n > f2/f1) wrapped onto
the causal linear response and read back a spurious copy of H1 in
harmonic_responses. swept_sine_distortion now rejects the request with
a sizing error naming the sweep ratio; the synchronized path is
unaffected (its padded circular buffer keeps absent orders at the
deconvolution floor, pinned by a new regression test).
- The duration-quantization bound in the synchronized_sweep_signal
docstring was wrong by a factor ln(f2/f1): the rounding perturbs the
duration by up to ln(f2/f1)/2 periods of f1, not half a period.
- An explicitly tiny ir_length now gets its own message instead of the
arrival-spacing diagnosis meant for the default-window path.
- The guides no longer attribute the |H1|/|H2| Chebyshev anchors to the
conformance suite alone (they are pinned by the test suite; conformance
pins |H3|, its phase, the THD and the linear floor).
- _sweep_rate drops an unused parameter; the thd_frequencies field doc no
longer cross-links to the unrelated tone-by-tone thd() function.
- Regenerated API reference (including the room_ir cross-reference to the
new analyzer) and llms files.
* fix(electroacoustics): name the actual culprit in the Farina sizing error
The anticausal-span guard raised the same "orders beyond f2/f1" message
whether the order itself did not fit or only its window did. The two
causes now get their own diagnosis and remedy: order beyond the ratio
(lower n_harmonics or widen the band) versus a window that overruns the
span around an order that fits (shorten ir_length, lower n_harmonics or
lengthen the sweep). Both branches pinned by tests.
* fix: PR review round - Sonar smells and bot findings
Sonar:
- S3776 (swept_sine_distortion, complexity 26): the validation cascade
moves into _validated_recording, _sweep_timing, _harmonic_window and
_check_farina_span with identical semantics and messages; every pinned
oracle (Chebyshev magnitudes/phases, THD closed form, min-phase 1e-12)
is unchanged.
- S1192 (_plot/electroacoustics.py): the three "Frequency [Hz]" literals
become the module constant _FREQ_LABEL, matching _plot/metrology.py.
- S1172: _window_harmonics drops its unused fs parameter.
- S5778 (test_swept_sine.py, test_phase.py): the array constructions move
out of the pytest.raises blocks so each block holds one statement.
CodeRabbit:
- synchronized_sweep_signal now validates 'amplitude' like the analyzer
does (zero, negative or non-finite raise); regression test.
- _thd_curves raises an actionable ValueError when no excitation bin has
a measurable order-2 product (f1 above fs/4) instead of returning
silently empty THD arrays; the farina wording was dropped from the
message because that path is unreachable there (the anticausal guard
already caps farina orders at the sweep ratio, which excludes
sub-octave sweeps). Regression test and raises-doc updated.
- The guide example asks for n_harmonics=3 so the highest product stays
below Nyquist, consistent with its own sizing rule (EN/ES, three
trees).
- The intro no longer claims a sweep measures a nonlinear system
"completely": it characterises the linear response and every harmonic
order at once (guide three trees + module docstring).
- excess_group_delay is described as the all-pass part, not "the
pure-latency part", and the phase-utility input note now states that
only minimum_phase accepts a plain magnitude array (three trees).
- The figure asymptote labels name what is drawn, (a2/2)/H1 and
(a3/4)/H1, in the generator, the ES translation table and both guide
snippets; the four swept_sine_thd variants regenerate.
- api-reference group_delay example comment reads "# seconds".
- Rejected: exact-integer validation for n_harmonics/ir_length/oversample
(int() conversion is the library-wide convention, e.g.
envelope.decimation_factor and correlation upsample); the llms-full
loudness wording (pre-existing docs/loudness.md prose outside this PR,
already hedged and followed by the exact phon/sone definition).
Gemini (all four rejected with evidence): no English fragment exists in
the ES guide (the flagged line is fully Spanish); the phase sentence
already uses "queda totalmente determinada"; "paso-todo" follows the
repo-wide ES convention (paso bajo/paso banda/paso alto in the existing
guides), not "pasa-todo".
Regenerated: API reference, llms files, conformance byte-identical
(306/306), figures within check_figures tolerance.