[READ-ONLY] Mirror of https://github.com/jmrplens/PyOctaveBand. [Python3] Octave-Band and Fractional Octave-Band filter. For signal in time domain. jmrplens.github.io/PyOctaveBand/
acoustics audio filter frequency frequency-analysis frequency-domain octave python3 signal time-domain
0

Configure Feed

Select the types of activity you want to include in your feed.

Correct DSP estimator edge cases: TSA time axis, echo detection, EQ validation (#356)

* Report the true 1/fs time axis from time_synchronous_average

The averaged samples sit on the 1/fs sampling grid (each block is aligned
back to it by the band-limited fractional delay before averaging), but the
result reported a time axis with T/M spacing. For a non-integer number of
samples per period the two grids drift apart across the period: with
100.37 samples per period the noiseless recovery error read against the
reported axis was 9.3e-2, versus 2.8e-6 on the true m/fs grid.

Return times = arange(M)/fs, document that the M samples cover the period
only approximately when fs*T is not an integer (McFadden 1987 revised
model: a block average on the sampling grid, no angular resampling), and
assert the non-integer-period tests on the true grid, with a regression
test pinning the 100.37-sample case on both grids. Also soften the comb
response [0, 1] claim (float cancellation beside a tooth can exceed 1 by
a few 1e-9) and mark the ~14 dB N = 32 figure as derived from Eq. 8
rather than printed in the paper.

* Detect inverting echoes and off-sample delays in echo_detection

echo_detection picked the maximum of the signed power cepstrum, so an
inverting reflection (a < 0), whose first rahmonic is a itself and hence
negative, was silently missed: with a = -0.4 at 313 samples the true
rahmonic c[313] = -0.42 was ignored and an unrelated 4739-sample peak
returned. Peak-pick on |cepstrum| within the searched band instead and
report the signed cepstrum value at the peak as the reflection
coefficient (the Mercator series ln(1 + a e^{-j theta}) holds for either
sign of a).

Also refine the reported delay by three-point parabolic interpolation of
|cepstrum| through the peak, so an off-sample delay (e.g. 313.5 samples)
is estimated to a fraction of a sample, and document the bin-splitting
caveat: at off-sample delays the peak-sample coefficient underestimates
|a| (about 65 % of it midway between bins). One-line convention note
distinguishing Milner's IDFT(ln|X|^2) power cepstrum from Bogert's
squared original. Tests cover a < 0 (exact impulse-echo closed form and
a noise source) and the interpolated half-sample delay.

* Validate parametric-EQ sections at the numerical extremes

Three validation gaps in the RBJ cookbook designer:

- A wide 'bw' at an f0 near Nyquist drives the w0/sin(w0) warping factor
into math.sinh overflow, escaping as a raw OverflowError; slightly less
extreme pairs produced a large but finite alpha whose rounded a2 landed
exactly at -1.0: poles on the unit circle designed silently. Wrap the
sinh with an informative ValueError and check the realized coefficients
against the stability triangle (|a2| < 1, |a1| < 1 + a2) after
normalization.
- A shelf slope beyond the gain-dependent bound (A + 1/A)/(A + 1/A - 2)
makes the cookbook's alpha complex and surfaced as a bare 'math domain
error'. Validate the bound at section construction with the limit in
the message (any slope is fine at 0 dB, where A = 1).
- Non-finite f0/gain_db/q/bw/slope passed validation silently (NaN
compares false against every bound) and produced NaN/inf coefficient
rows. Require finiteness explicitly.

Tests cover all three paths plus extreme-but-valid sections that must
keep designing stably.

* Scope the envelope-spectrum amplitude claim and add the band-pass front end

The docstring promised that a sinusoidal modulation 'reads out as a line
at its exact amplitude'; that holds only when the modulation frequency
falls on an analysis bin. Off-bin lines read low by the taper's
scalloping loss (about 1.42 dB, ~15 %, for the default Hann midway
between bins), like any single-record amplitude spectrum. State the
on-bin condition on the closed forms and the caveat in the module,
function and result docstrings, pinned by a midway-between-bins test.

Also add the optional band=(low, high) pre-filter to envelope_spectrum:
a zero-phase 4th-order Butterworth band-pass ahead of the detector, the
band-pass front end of Bendat & Piersol Figure 13.11 and the classical
bearing-envelope chain (isolate the excited resonance band, then
envelope it). The band travels with the result; a test shows it
recovering the exact modulation line under a strong out-of-band
interferer that swamps the raw envelope.

* Warn when a tone burst cannot close on the tone's zero crossing

IEC 60268-1 Clause A2.1 asks for an integral number of full periods, but
the gate closes after round(fs*cycles/frequency) samples: sample-exact
only when f and fs are commensurate. For 10 cycles of 997 Hz at 48 kHz
the burst spans 481.44 samples, gated at 481, and the waveform carries a
residual step at the gate edge. Document the condition and emit a
PhonometryWarning quantifying the realized residual (bounded by
amplitude*sin(pi*f/fs)) for incommensurate pairs; commensurate bursts
stay silent.

Also tighten the resampler's 'equiripple-bounded' wording: the Kaiser
window method bounds both bands by the same delta = 10^(-A/20), but its
ripple decays away from the band edges, unlike a true equiripple design.

* Document estimator conventions verified against their sources

Documentation and comment corrections, no behavior changes:

- random_data: keep the reverse-arrangement acceptance bounds on the
normal approximation deliberately and say why: B&P Table A.6 itself
follows it, and no percentage-point convention on the exact Mahonian
distribution reproduces the published entries (N = 40 and N = 100 at
a = 0.975 pin it down), so exact bounds would disagree with the table
users check against. The p-value stays exact; the one-count boundary
disagreement is now documented on TrendTestResult.
- time_frequency: the spectrogram energy-recovery statement was wrong as
written; integrating 'density' columns over frequency and summing over
time recovers the record energy only after multiplying by the hop
duration, up to the taper roll-off at the record edges (verified
0.988 vs 1 for Hann at 75 % overlap on white noise).
- spectra: note that the adaptive multitaper sigma^2 is mean(x^2) with no
mean removal (consistent with the documented no-detrending
calibration), and mark the 'P&W Eq. 369a' citation as pending
verification against the book.
- compliance + tests/reference_data: transcription note for the bare
'+3' at 20 Hz, Type 2, in ANSI S1.4-1983 Table V (read as upper-only;
+/-3 plausible; either way the verdict cannot change, the realized
deviation is -0.05 dB).
- room_ir: the Mueller & Massarani 'below 4 dB' crest-factor figure
assumes their slight magnitude smoothing, which shaped_sweep_signal
does not implement; measured sweeps here sit at 4.2-4.3 dB.

* Translate EQ section labels, complete the bibliography and changelog

- The per-section legend labels of the parametric-EQ plot printed the raw
English filter-type names in Spanish figures; route them through the
plot string table (peaking -> campana, shelves, band-pass, notch,
all-pass; lowpass/highpass were already translated for the lifter
plot). English output stays byte-identical.
- docs/references.md: add the Harris 1978 window paper and IEC
60268-1:1985, both cited by guide pages but missing from the collected
bibliography; link the IEC entry (and the one in test-signals.md) to
the IEC webstore page.
- Changelog entries for the DSP fixes; regenerate the API reference,
llms-full.txt and the cepstrum_echo figures (the detected-echo marker
moved by the sub-sample delay interpolation); the conformance report is
unchanged (419/419).

* Propagate the echo and envelope semantics to the site guides

Review follow-ups:

- The EN and ES site guides (cepstrum-echoes.mdx) still described the
old echo reading and the unqualified 'exact amplitude' claim; bring
them in step with docs/cepstrum-echoes.md: |cepstrum| peak-picking,
signed coefficient, interpolated delay, bin-splitting caveat, the
on-bin/scalloping condition and the new band= front end (the docs page
gains the same envelope-section text).
- State in the changelog and the result docstring that delay is no
longer exactly delay_samples/fs when the interpolated offset is
nonzero (delay_samples remains the whole-sample index the coefficient
is read at), and merge the duplicated Unreleased headings.
- Clearer stability message (poles 'do not lie strictly inside the unit
circle', full-precision coefficients) and an informative short-record
error for the envelope band-pass pre-filter instead of scipy's padlen
message.
- Sharper tone-burst warning test (asserts the quantified span and
residual) and a tighter interpolated-delay tolerance.

* Address review notes: cepstrum terminology, band-pass claims, edge validation

- Guides (docs and site, EN and ES): call the signed quantity what it is,
Milner's cepstrum of the log-power spectrum, and spell out that
Bogert's original 1963 power cepstrum squares once more and is
non-negative, so the signed rahmonics are the library's documented
convention; qualify the band-pass front end (out-of-band content is
strongly attenuated, not eliminated: finite Butterworth rejection and
small zero-phase edge transients).
- cepstrum: render reflection_coefficient as inline code in the
delay_samples docstring so the generated API page stops linking the
unrelated underwater seabed function.
- envelope: validate that 'band' is exactly two numeric edges before
indexing (a one-item tuple raised IndexError, extra entries were
ignored); the validated edges now travel to the result.
- equalizer: reject finite but unrepresentable gains through a shared
_gain_amplitude helper (10^(gain_db/40) overflowing, or underflowing
to zero and then dividing by zero in the designers); document the 0 dB
exception to the shelf-slope bound on the slope field.
- tone_burst: validate the repetition configuration before emitting the
incommensurate-frequency warning, so invalid setups raise even under
a warnings-as-errors filter.
- comb_filter_response: reject phase products that overflow to infinity
(sin(inf) returned NaN instead of the documented bounded response).

Tests for every new rejection path; API reference and llms regenerated.

* Test the representable gain amplitude as a positive finite bound

* Scope the first-rahmonic height claim to the source-free identity

authored by

José M. Requena Plens and committed by
GitHub
(Jul 25, 2026, 12:36 PM +0200) 77a55e7e c7bfded3

+936 -153
+34
CHANGELOG.md
··· 74 74 impulse-response onset detector is ISO 3382-1:2009 A.3.4, not A.2.1. 75 75 - English fiche strings and docstrings for ISO 3382-3 use point decimals 76 76 ("STI = 0.50"/"0.20"); the Spanish strings keep the comma. 77 + - `metrology.time_synchronous_average`: the result's `times` axis now 78 + reports the true `1/fs` sampling grid the averaged samples sit on, 79 + instead of a `T/M` spacing that drifted across the period when 80 + `fs*period` is not an integer (noiseless recovery read against the old 81 + axis was wrong by up to ~9e-2 for a 100.37-sample period; on the true 82 + grid the error is the fractional-delay interpolation level, ~3e-6). 83 + - `metrology.echo_detection`: the peak is now picked on `|cepstrum|`, so 84 + an inverting reflection (negative reflection coefficient), whose first 85 + rahmonic is negative, is detected at its true delay instead of being 86 + silently missed; the signed cepstrum value at the peak is returned as 87 + the reflection coefficient, and the reported delay is refined by 88 + three-point parabolic interpolation so off-sample delays are estimated 89 + to a fraction of a sample (with the bin-splitting caveat on the 90 + coefficient documented). Note that `delay_samples` stays the whole 91 + peak-sample index, so `delay == delay_samples/fs` no longer holds 92 + exactly when the interpolated offset is nonzero. 93 + - `metrology.EQSection`/`ParametricEQ`: informative `ValueError`s replace 94 + a raw `OverflowError` (bandwidth-derived alpha overflowing near 95 + Nyquist), a bare `math domain error` (shelf slope beyond the 96 + gain-dependent bound `(A + 1/A)/(A + 1/A - 2)`) and two silent 97 + failures (sections whose rounded poles land exactly on the unit 98 + circle, and non-finite `f0`/`gain_db`/`q`/`bw`/`slope` values passing 99 + validation). 100 + - `metrology.tone_burst`: warns when `frequency` is incommensurate with 101 + `fs`, in which case the gate cannot close exactly at the tone's final 102 + zero crossing (the "integral number of full periods" of IEC 60268-1 103 + Clause A2.1 is then sample-exact only to the nearest sample) and the 104 + waveform carries a small residual step at the gate edge. 77 105 78 106 ### Fixed 79 107 ··· 248 276 each kind must be fed with for a Directive-conforming assessment: the Part A 249 277 vector total `a_hv` for hand-arm, the Part B dominant-axis value for 250 278 whole-body (not the ISO 2631-1 Eq. (10) vector total `a_v`). 279 + - `metrology.envelope_spectrum`: optional `band=(low, high)` zero-phase 280 + band-pass pre-filter ahead of the envelope detector -- the band-pass 281 + front end of Bendat & Piersol Figure 13.11 and the classical 282 + bearing-envelope chain; the band travels with the result. The 283 + amplitude-calibration docs now state the on-bin condition and the 284 + taper scalloping loss for off-bin modulation lines. 251 285 - `materials.diffuser_design`: far-field polar-response and diffusion-coefficient 252 286 prediction of a Schroeder phase-grating diffuser from its surface design, 253 287 complementing the measurement-only `materials.scattering_diffusion`.
+32 -8
llms-full.txt
··· 11359 11359 11360 11360 - International Electrotechnical Commission. (1985). *Sound system 11361 11361 equipment — Part 1: General* (IEC 60268-1:1985). 11362 + [IEC webstore](https://webstore.iec.ch/en/publication/1204). 11362 11363 Annex A, Clause A2: tone bursts starting at the zero crossing of the tone 11363 11364 with an integral number of full periods (A2.1), repetitive burst trains 11364 11365 at a stated repetition rate (A2.2), and the Table AII burst durations the ··· 11396 11397 11397 11398 - `'power'` (default): the inverse DFT of `ln|X|²` (Milner's Fig. 21). Real, 11398 11399 even and phase-blind - the workhorse for echo and harmonic-family 11399 - detection; 11400 + detection. This is Milner's *signed* cepstrum of the log-power spectrum; 11401 + Bogert's original 1963 "power cepstrum" squares once more and is 11402 + non-negative - the library follows Milner throughout, so negative 11403 + rahmonics keep their sign; 11400 11404 - `'real'`: the inverse DFT of `ln|X|` - exactly half the power cepstrum, 11401 11405 and the quantity whose causal folding is the minimum-phase reconstruction 11402 11406 (below); ··· 11439 11443 11440 11444 so the cepstrum carries a spike train at the **rahmonics** `n·t0` with 11441 11445 amplitudes `a, -a²/2, a³/3, ...` (their sum is `ln(1+a)`), regardless of the 11442 - spectrum of `s` itself, which concentrates at low quefrencies. On the power 11443 - cepstrum the first spike's height is exactly the reflection coefficient `a` - 11444 - a closed form the tests and the conformance suite pin to 1e-10. 11445 - `echo_detection` automates the reading: 11446 + spectrum of `s` itself, which concentrates at low quefrencies. On the signed 11447 + cepstrum of the log-power spectrum (`kind='power'`, Milner's convention) the 11448 + first spike's height is the reflection coefficient `a` - of either 11449 + sign - plus whatever the source cepstrum contributes at that quefrency 11450 + (negligible for broadband sources, whose cepstrum concentrates at low 11451 + quefrencies); on an ideal impulse-plus-echo the identity is a closed form 11452 + the tests and the conformance suite pin to 1e-10. `echo_detection` automates the reading: it picks the largest 11453 + `|cepstrum|` peak in the searched band (so an inverting reflection, `a < 0`, 11454 + is found at its true delay rather than missed), refines the delay by 11455 + quadratic interpolation through the peak and its neighbours, and reports the 11456 + signed peak value as the reflection coefficient. When the true delay falls 11457 + between samples the rahmonic splits across quefrency bins: the interpolated 11458 + delay still lands on it, but the reported coefficient underestimates `|a|` 11459 + (down to about 65 % of it midway between samples). 11446 11460 11447 11461 ```python 11448 11462 import numpy as np ··· 11572 11586 DC remover, and a spectral view of what remains. `envelope_spectrum` runs the 11573 11587 [Hilbert envelope](https://jmrplens.github.io/phonometry/guides/correlation-delay/) (`kind="magnitude"`, 11574 11588 the practical default) or the book's square-law detector (`kind="squared"`) 11575 - through exactly that chain, scaled so a sinusoidal modulation reads out as a 11576 - line at its exact amplitude. 11589 + through exactly that chain, scaled by the taper's coherent gain so a 11590 + sinusoidal modulation whose frequency falls on an analysis bin reads out as 11591 + a line at its exact amplitude (off-bin lines read low by the taper's 11592 + scalloping loss, up to about 1.4 dB for the default Hann). The optional 11593 + `band=(low, high)` argument reproduces the figure's band-pass front end - 11594 + the classical bearing-envelope chain: isolate the structural-resonance band 11595 + excited by the defect impacts with a zero-phase band-pass, then envelope 11596 + it - so an out-of-band interferer is strongly attenuated before the 11597 + detector (the roll-off of a fourth-order Butterworth applied forward and 11598 + backward; the rejection is finite, and the zero-phase pass leaves small 11599 + transients at the record edges). 11577 11600 11578 - For an AM tone `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` the closed forms are: 11601 + For an AM tone `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` with `fm` on an 11602 + analysis bin the closed forms are: 11579 11603 11580 11604 | `kind` | mean level | line at `fm` | line at `2fm` | 11581 11605 |---|---|---|---|
+31 -8
docs/cepstrum-echoes.md
··· 22 22 23 23 - `'power'` (default): the inverse DFT of `ln|X|²` (Milner's Fig. 21). Real, 24 24 even and phase-blind - the workhorse for echo and harmonic-family 25 - detection; 25 + detection. This is Milner's *signed* cepstrum of the log-power spectrum; 26 + Bogert's original 1963 "power cepstrum" squares once more and is 27 + non-negative - the library follows Milner throughout, so negative 28 + rahmonics keep their sign; 26 29 - `'real'`: the inverse DFT of `ln|X|` - exactly half the power cepstrum, 27 30 and the quantity whose causal folding is the minimum-phase reconstruction 28 31 (below); ··· 65 68 66 69 so the cepstrum carries a spike train at the **rahmonics** `n·t0` with 67 70 amplitudes `a, -a²/2, a³/3, ...` (their sum is `ln(1+a)`), regardless of the 68 - spectrum of `s` itself, which concentrates at low quefrencies. On the power 69 - cepstrum the first spike's height is exactly the reflection coefficient `a` - 70 - a closed form the tests and the conformance suite pin to 1e-10. 71 - `echo_detection` automates the reading: 71 + spectrum of `s` itself, which concentrates at low quefrencies. On the signed 72 + cepstrum of the log-power spectrum (`kind='power'`, Milner's convention) the 73 + first spike's height is the reflection coefficient `a` - of either 74 + sign - plus whatever the source cepstrum contributes at that quefrency 75 + (negligible for broadband sources, whose cepstrum concentrates at low 76 + quefrencies); on an ideal impulse-plus-echo the identity is a closed form 77 + the tests and the conformance suite pin to 1e-10. `echo_detection` automates the reading: it picks the largest 78 + `|cepstrum|` peak in the searched band (so an inverting reflection, `a < 0`, 79 + is found at its true delay rather than missed), refines the delay by 80 + quadratic interpolation through the peak and its neighbours, and reports the 81 + signed peak value as the reflection coefficient. When the true delay falls 82 + between samples the rahmonic splits across quefrency bins: the interpolated 83 + delay still lands on it, but the reported coefficient underestimates `|a|` 84 + (down to about 65 % of it midway between samples). 72 85 73 86 ```python 74 87 import numpy as np ··· 198 211 DC remover, and a spectral view of what remains. `envelope_spectrum` runs the 199 212 [Hilbert envelope](correlation-delay.md) (`kind="magnitude"`, 200 213 the practical default) or the book's square-law detector (`kind="squared"`) 201 - through exactly that chain, scaled so a sinusoidal modulation reads out as a 202 - line at its exact amplitude. 214 + through exactly that chain, scaled by the taper's coherent gain so a 215 + sinusoidal modulation whose frequency falls on an analysis bin reads out as 216 + a line at its exact amplitude (off-bin lines read low by the taper's 217 + scalloping loss, up to about 1.4 dB for the default Hann). The optional 218 + `band=(low, high)` argument reproduces the figure's band-pass front end - 219 + the classical bearing-envelope chain: isolate the structural-resonance band 220 + excited by the defect impacts with a zero-phase band-pass, then envelope 221 + it - so an out-of-band interferer is strongly attenuated before the 222 + detector (the roll-off of a fourth-order Butterworth applied forward and 223 + backward; the rejection is finite, and the zero-phase pass leaves small 224 + transients at the record edges). 203 225 204 - For an AM tone `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` the closed forms are: 226 + For an AM tone `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` with `fm` on an 227 + analysis bin the closed forms are: 205 228 206 229 | `kind` | mean level | line at `fm` | line at `2fm` | 207 230 |---|---|---|---|
+12
docs/references.md
··· 66 66 Cited by [Calibrated Spectral Analysis](spectral-analysis.md), 67 67 [Correlation, Time Delay and Envelope](correlation-delay.md) and 68 68 [Swept-sine distortion and phase utilities](swept-sine-distortion.md). 69 + - Harris, F. J. (1978). On the use of windows for harmonic analysis with the 70 + discrete Fourier transform. *Proceedings of the IEEE*, 66(1), 51-83. 71 + [doi:10.1109/PROC.1978.10837](https://doi.org/10.1109/PROC.1978.10837). 72 + The window figures of merit (Table 1) computed by `window_metrics`. 73 + Cited by [Calibrated Spectral Analysis](spectral-analysis.md). 69 74 - Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. 70 75 *Proceedings of the IEEE*, 70(9), 1055-1096. 71 76 [doi:10.1109/PROC.1982.12433](https://doi.org/10.1109/PROC.1982.12433). ··· 138 143 The calibrator classes, level tolerances and the short-term stability 139 144 criterion applied to calibration recordings. 140 145 Cited by [Calibration and dBFS](calibration.md). 146 + - International Electrotechnical Commission. (1985). *Sound system 147 + equipment — Part 1: General* (IEC 60268-1:1985). 148 + [IEC webstore](https://webstore.iec.ch/en/publication/1204). 149 + Annex A, Clause A2: tone bursts starting at the zero crossing of the tone 150 + with an integral number of full periods, repetitive burst trains at a 151 + stated repetition rate, and the Table AII burst durations. 152 + Cited by [Test signals and sample-rate tools](test-signals.md). 141 153 142 154 ## Broadcast programme loudness 143 155
+1
docs/test-signals.md
··· 139 139 140 140 - International Electrotechnical Commission. (1985). *Sound system 141 141 equipment — Part 1: General* (IEC 60268-1:1985). 142 + [IEC webstore](https://webstore.iec.ch/en/publication/1204). 142 143 Annex A, Clause A2: tone bursts starting at the zero crossing of the tone 143 144 with an integral number of full periods (A2.1), repetitive burst trains 144 145 at a stated repetition rate (A2.2), and the Table AII burst durations the
+4
tests/reference_data.py
··· 157 157 # response levels for Type 0 (laboratory), Type 1 (precision) and Type 2 158 158 # (general purpose) instruments; they apply to every weighting. A lower 159 159 # limit of -inf means only the upper limit applies. 160 + # Transcription note, 20 Hz Type 2: the standard prints a bare "+3"; read 161 + # as +3/upper-only (like the surrounding upper-only rows), with +/-3 a 162 + # plausible alternative. The realized B response there is only 0.05 dB 163 + # below nominal, so the reading cannot change any verdict. 160 164 # Row = (freq_Hz, t0_up, t0_lo, t1_up, t1_lo, t2_up, t2_lo). 161 165 ANSIS14_TABLE5: list[tuple[float, float, float, float, float, float, float]] = [ 162 166 (10, 2.0, -5.0, 4.0, -4.0, 5.0, -INF),
+1 -1
.github/images/cepstrum_echo.svg
··· 1652 1652 " style="stroke: #d62728; stroke-linejoin: miter"/> 1653 1653 </defs> 1654 1654 <g clip-path="url(#p11e2239aa4)"> 1655 - <use xlink:href="#me284da6c85" x="224.930813" y="89.260908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1655 + <use xlink:href="#me284da6c85" x="224.930825" y="89.260908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1656 1656 </g> 1657 1657 </g> 1658 1658 <g id="patch_4">
+1 -1
.github/images/cepstrum_echo_dark.svg
··· 1652 1652 " style="stroke: #d62728; stroke-linejoin: miter"/> 1653 1653 </defs> 1654 1654 <g clip-path="url(#p11e2239aa4)"> 1655 - <use xlink:href="#me284da6c85" x="224.930813" y="89.260908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1655 + <use xlink:href="#me284da6c85" x="224.930825" y="89.260908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1656 1656 </g> 1657 1657 </g> 1658 1658 <g id="patch_4">
+1 -1
.github/images/cepstrum_echo_es.svg
··· 1652 1652 " style="stroke: #d62728; stroke-linejoin: miter"/> 1653 1653 </defs> 1654 1654 <g clip-path="url(#pfc73f456aa)"> 1655 - <use xlink:href="#me284da6c85" x="220.158938" y="89.740908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1655 + <use xlink:href="#me284da6c85" x="220.15895" y="89.740908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1656 1656 </g> 1657 1657 </g> 1658 1658 <g id="patch_4">
+1 -1
.github/images/cepstrum_echo_es_dark.svg
··· 1652 1652 " style="stroke: #d62728; stroke-linejoin: miter"/> 1653 1653 </defs> 1654 1654 <g clip-path="url(#pfc73f456aa)"> 1655 - <use xlink:href="#me284da6c85" x="220.158938" y="89.740908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1655 + <use xlink:href="#me284da6c85" x="220.15895" y="89.740908" style="fill: #d62728; stroke: #d62728; stroke-linejoin: miter"/> 1656 1656 </g> 1657 1657 </g> 1658 1658 <g id="patch_4">
+44
tests/metrology/test_cepstrum.py
··· 173 173 assert res.reflection_coefficient == pytest.approx(0.25, abs=0.02) 174 174 175 175 176 + def test_echo_detection_negative_reflection_exact() -> None: 177 + """An inverting echo (a < 0) is a *negative* rahmonic of height a. 178 + 179 + The Mercator series holds for either sign of ``a``: the first rahmonic 180 + of ``delta[n] - 0.4*delta[n-d]`` is exactly ``-0.4`` at quefrency 181 + ``d``. Peak-picking on ``|cepstrum|`` must find the true delay and 182 + report the signed coefficient. 183 + """ 184 + res = ph.echo_detection(_impulse_echo(a=-ALPHA), FS) 185 + assert res.delay_samples == DELAY 186 + assert res.delay == pytest.approx(DELAY / FS) 187 + assert res.reflection_coefficient == pytest.approx(-ALPHA, abs=1e-12) 188 + 189 + 190 + def test_echo_detection_negative_reflection_on_noise_source() -> None: 191 + s = ph.noise_signal(FS, N / FS, color="white", seed=42) 192 + x = s - 0.25 * np.roll(s, 200) 193 + res = ph.echo_detection(x, FS, min_quefrency=0.005) 194 + assert res.delay_samples == 200 195 + assert res.reflection_coefficient == pytest.approx(-0.25, abs=0.02) 196 + 197 + 198 + def test_echo_detection_off_sample_delay_is_interpolated() -> None: 199 + """A delay midway between samples: quadratic interpolation locates it. 200 + 201 + The echo is synthesised circularly by an exact half-sample phase ramp, 202 + so the true delay is 313.5 samples. The rahmonic splits across the 203 + neighbouring bins: the interpolated ``delay`` recovers the fraction, 204 + while the coefficient at the peak sample underestimates ``a`` (the 205 + documented bin-splitting caveat). 206 + """ 207 + rng = np.random.default_rng(7) 208 + s = rng.standard_normal(N) 209 + true_delay = DELAY + 0.5 210 + ramp = np.exp(-2j * np.pi * np.fft.rfftfreq(N) * true_delay) 211 + x = s + ALPHA * np.fft.irfft(np.fft.rfft(s) * ramp, N) 212 + 213 + res = ph.echo_detection(x, FS) 214 + assert res.delay_samples in (DELAY, DELAY + 1) 215 + assert res.delay * FS == pytest.approx(true_delay, abs=0.15) 216 + # Bin splitting: the peak-sample coefficient reads low but keeps sign. 217 + assert 0.15 < res.reflection_coefficient < ALPHA 218 + 219 + 176 220 def test_echo_detection_search_band_is_respected() -> None: 177 221 # Restrict the band away from the true echo: the peak reported must 178 222 # come from inside the band.
+60
tests/metrology/test_envelope.py
··· 291 291 ) 292 292 293 293 294 + def test_envelope_spectrum_off_bin_line_shows_hann_scalloping() -> None: 295 + """An off-bin modulation frequency reads low by the scalloping loss. 296 + 297 + Midway between bins the Hann taper's scalloping loss is about 298 + 1.42 dB (a factor ~0.85): the documented limit of the on-bin 'exact 299 + amplitude' calibration. 300 + """ 301 + df = FS / ES_N 302 + fm = (round(ES_FM / df) + 0.5) * df # exactly midway between bins 303 + t = np.arange(ES_N) / FS 304 + x = ( 305 + ES_A0 * (1.0 + ES_M * np.cos(2.0 * np.pi * fm * t)) 306 + * np.cos(2.0 * np.pi * 1000.0 * t) 307 + ) 308 + res = ph.envelope_spectrum(x, FS) 309 + peak = float(np.max(res.amplitude[1:])) 310 + # Hann scalloping at the bin midpoint: |W(1/2)| / |W(0)| ~ 0.8488. 311 + assert peak == pytest.approx(0.8488 * ES_A0 * ES_M, rel=0.02) 312 + assert peak < 0.9 * ES_A0 * ES_M 313 + 314 + 315 + def test_envelope_spectrum_bandpass_prefilter_isolates_the_carrier() -> None: 316 + """The optional band isolates the modulated carrier before detection. 317 + 318 + With a strong out-of-band interfering tone the raw Hilbert envelope is 319 + dominated by the two-tone beat; band-passing around the carrier first 320 + (the classical bearing-envelope chain) recovers the clean modulation 321 + line at its exact amplitude. 322 + """ 323 + t = np.arange(ES_N) / FS 324 + x = _am_signal() + 3.0 * np.cos(2.0 * np.pi * 3000.0 * t) 325 + 326 + raw = ph.envelope_spectrum(x, FS) 327 + filtered = ph.envelope_spectrum(x, FS, band=(700.0, 1300.0)) 328 + 329 + assert filtered.band == (700.0, 1300.0) 330 + assert raw.band is None 331 + assert filtered.amplitude[_bin(ES_FM)] == pytest.approx( 332 + ES_A0 * ES_M, rel=1e-2 333 + ) 334 + assert filtered.mean_level == pytest.approx(ES_A0, rel=1e-2) 335 + # Without the pre-filter the beat envelope swamps the modulation line. 336 + beat_error = abs(raw.amplitude[_bin(ES_FM)] - ES_A0 * ES_M) 337 + assert beat_error > 0.05 * ES_A0 * ES_M 338 + 339 + 340 + def test_envelope_spectrum_band_validation() -> None: 341 + x = _am_signal() 342 + for band in ((0.0, 100.0), (200.0, 100.0), (100.0, FS / 2.0)): 343 + with pytest.raises(ValueError, match="band"): 344 + ph.envelope_spectrum(x, FS, band=band) 345 + # Malformed shapes: not exactly two numeric edges. 346 + for bad in ((700.0,), (700.0, 900.0, 1300.0), (700.0, None)): 347 + with pytest.raises(ValueError, match="pair"): 348 + ph.envelope_spectrum(x, FS, band=bad) # type: ignore[arg-type] 349 + # A record shorter than the zero-phase padding fails informatively. 350 + with pytest.raises(ValueError, match="too short"): 351 + ph.envelope_spectrum(x[:20], FS, band=(700.0, 1300.0)) 352 + 353 + 294 354 def test_envelope_spectrum_validates_inputs() -> None: 295 355 x = _am_signal() 296 356 with pytest.raises(ValueError, match="kind"):
+2
tests/metrology/test_metrology_plot_i18n.py
··· 331 331 axes = res.plot(language="es") 332 332 assert "Respuesta del EQ paramétrico (Audio EQ Cookbook)" in _titles(axes) 333 333 assert "Cascada" in _labels(axes) 334 + assert "shelf de graves 100 Hz" in _labels(axes) 335 + assert "campana 1000 Hz" in _labels(axes) 334 336 assert "Fase [grados]" in _labels(axes) 335 337 plt.close("all") 336 338 with pytest.raises(ValueError):
+73
tests/metrology/test_parametric_eq.py
··· 460 460 ph.ParametricEQ(FS, []) 461 461 462 462 463 + @pytest.mark.parametrize( 464 + "kwargs", 465 + [ 466 + {"f0": float("nan"), "gain_db": 3.0}, 467 + {"f0": float("inf"), "gain_db": 3.0}, 468 + {"f0": 1000.0, "gain_db": float("nan")}, 469 + {"f0": 1000.0, "gain_db": float("inf")}, 470 + {"f0": 1000.0, "gain_db": 3.0, "q": float("nan")}, 471 + {"f0": 1000.0, "gain_db": 3.0, "bw": float("inf")}, 472 + ], 473 + ) 474 + def test_rejects_non_finite_parameters(kwargs: dict) -> None: 475 + with pytest.raises(ValueError, match="finite"): 476 + ph.EQSection("peaking", **kwargs) 477 + 478 + 479 + @pytest.mark.parametrize("gain_db", [1e308, -1e308]) 480 + @pytest.mark.parametrize( 481 + ("filter_type", "kwargs"), 482 + [("peaking", {}), ("lowshelf", {"slope": 1.0})], 483 + ) 484 + def test_extreme_finite_gain_raises_value_error( 485 + gain_db: float, filter_type: str, kwargs: dict 486 + ) -> None: 487 + # 10^(gain_db/40) overflows for huge positive gains and underflows to 488 + # zero (then divides by zero in the designers) for huge negative ones; 489 + # both used to leak raw OverflowError/ZeroDivisionError. 490 + with pytest.raises(ValueError, match="representable"): 491 + ph.EQSection(filter_type, 1000.0, gain_db=gain_db, **kwargs) 492 + 493 + 494 + def test_shelf_slope_beyond_gain_bound_raises_informatively() -> None: 495 + # A = 10^(9/40): the alpha recipe is real only below 496 + # (A + 1/A)/(A + 1/A - 2) = 8.2869...; beyond it the raw math would 497 + # fail with a bare "math domain error". 498 + with pytest.raises(ValueError, match="too steep"): 499 + ph.EQSection("lowshelf", 1000.0, gain_db=9.0, slope=12.0) 500 + # Just below the bound the design succeeds and is strictly stable. 501 + section = ph.EQSection("lowshelf", 1000.0, gain_db=9.0, slope=8.0) 502 + sos = ph.ParametricEQ(FS, section).sos 503 + assert np.max(np.abs(np.roots(sos[0, 3:]))) < 1.0 504 + # Gain 0 dB (A = 1): every positive slope is admissible. 505 + ph.EQSection("lowshelf", 1000.0, gain_db=0.0, slope=50.0) 506 + 507 + 508 + def test_bw_overflow_near_nyquist_raises_value_error() -> None: 509 + # The w0/sin(w0) warping factor diverges towards Nyquist: this bw/f0 510 + # pair used to escape as a raw OverflowError from math.sinh. 511 + section = ph.EQSection("notch", 23999.0, bw=1.0) 512 + with pytest.raises(ValueError, match="too wide"): 513 + ph.ParametricEQ(FS, section) 514 + 515 + 516 + def test_bw_marginally_unstable_section_raises_value_error() -> None: 517 + # Large but finite alpha: the rounded a2 lands exactly at -1.0 (poles 518 + # on the unit circle) with no arithmetic error to flag it. 519 + for f0, bw in ((23800.0, 1.0), (23900.0, 6.0)): 520 + with pytest.raises(ValueError, match="not strictly stable"): 521 + ph.ParametricEQ(FS, ph.EQSection("notch", f0, bw=bw)) 522 + 523 + 524 + def test_extreme_but_valid_sections_still_design_stably() -> None: 525 + for filter_type, f0, kwargs in ( 526 + ("peaking", 23990.0, {"gain_db": 40.0, "q": 100.0}), 527 + ("lowpass", 23999.0, {"q": 100.0}), 528 + ("highpass", 0.01, {"q": 0.001}), 529 + ("allpass", 12000.0, {"q": 1e-6}), 530 + ): 531 + sos = ph.ParametricEQ(FS, ph.EQSection(filter_type, f0, **kwargs)).sos 532 + assert np.all(np.isfinite(sos)) 533 + assert np.max(np.abs(np.roots(sos[0, 3:]))) < 1.0 534 + 535 + 463 536 def test_response_grid_validation() -> None: 464 537 eq = ph.ParametricEQ(FS, ph.EQSection("peaking", 1000.0, gain_db=3.0)) 465 538 with pytest.raises(ValueError, match="n_points"):
+41
tests/metrology/test_signal_toolbox.py
··· 113 113 assert res.period_samples == round(FS / rate) 114 114 115 115 116 + def test_tone_burst_incommensurate_frequency_warns() -> None: 117 + """997 Hz at 48 kHz: 10 periods span 481.44 samples, gated at 481. 118 + 119 + The 'integral number of full periods' of Clause A2.1 is sample-exact 120 + only for commensurate f/fs; otherwise a warning quantifies the 121 + residual step at the gate edge. 122 + """ 123 + with pytest.warns(ph.PhonometryWarning, match="incommensurate") as record: 124 + res = ph.tone_burst(FS, 997.0, 10) 125 + assert res.burst_samples == 481 126 + # The warning quantifies the exact span and the residual step at the 127 + # gate edge, sin(2*pi*997*481/48000) ~ -0.058 (a commensurate burst 128 + # would close exactly at the zero crossing). 129 + message = str(record[0].message) 130 + assert "481.444" in message 131 + assert "-0.058" in message 132 + 133 + 134 + def test_tone_burst_invalid_config_raises_before_warning() -> None: 135 + """An invalid repetition setup raises even for incommensurate bursts. 136 + 137 + Under a warnings-as-errors filter the incommensurate warning must not 138 + pre-empt the configuration ValueError. 139 + """ 140 + import warnings as _warnings 141 + 142 + with _warnings.catch_warnings(): 143 + _warnings.simplefilter("error", ph.PhonometryWarning) 144 + with pytest.raises(ValueError, match="repetition_rate"): 145 + ph.tone_burst(FS, 997.0, 10, repetitions=2) 146 + 147 + 148 + def test_tone_burst_commensurate_frequency_does_not_warn() -> None: 149 + import warnings as _warnings 150 + 151 + with _warnings.catch_warnings(): 152 + _warnings.simplefilter("error", ph.PhonometryWarning) 153 + res = ph.tone_burst(FS, 5000.0, 25) # 240 samples exactly 154 + assert res.burst_samples == 240 155 + 156 + 116 157 def test_tone_burst_result_is_frozen() -> None: 117 158 res = ph.tone_burst(FS, 5000.0, 25) 118 159 with pytest.raises(dataclasses.FrozenInstanceError):
+62 -7
tests/metrology/test_synchronous_average.py
··· 147 147 assert result.residual_rms < 1e-12 148 148 149 149 150 - def test_times_span_one_period() -> None: 150 + def test_times_are_the_sampling_grid() -> None: 151 151 one = _periodic(PERIOD, M, (2.0,)) 152 152 result = ph.time_synchronous_average(_repeat(one, 10), FS, PERIOD) 153 153 assert result.times.size == M 154 154 assert result.times[0] == pytest.approx(0.0, abs=1e-15) 155 155 assert result.times[-1] < PERIOD 156 + assert np.allclose(result.times, np.arange(M) / FS, atol=0.0) 156 157 157 158 158 159 # --------------------------------------------------------------------------- ··· 161 162 162 163 163 164 def test_noninteger_period_recovered_within_bound() -> None: 164 - """A non-integer M is aligned by band-limited fractional delay.""" 165 + """A non-integer M is aligned by band-limited fractional delay. 166 + 167 + The averaged samples stay on the ``1/fs`` sampling grid, so the 168 + reference is the true waveform evaluated at ``m/fs`` -- exactly the 169 + grid :attr:`times` reports. 170 + """ 165 171 period = 1.0 / 31.7 # FS * period is not an integer 166 172 m_int = round(FS * period) 167 173 phase = np.arange(30 * m_int) / FS / period ··· 171 177 result = ph.time_synchronous_average(signal, FS, period) 172 178 173 179 assert result.interpolated is True 174 - reference = np.cos( 175 - 2.0 * np.pi * np.arange(m_int) * (period / m_int) / period 176 - ) + 0.4 * np.cos( 177 - 2.0 * np.pi * 2.0 * np.arange(m_int) * (period / m_int) / period + 0.3 180 + assert np.allclose(result.times, np.arange(m_int) / FS, atol=0.0) 181 + ref_phase = result.times / period 182 + reference = np.cos(2.0 * np.pi * ref_phase) + 0.4 * np.cos( 183 + 2.0 * np.pi * 2.0 * ref_phase + 0.3 178 184 ) 179 - assert np.max(np.abs(result.period_waveform - reference)) < 0.05 185 + assert np.max(np.abs(result.period_waveform - reference)) < 1e-5 186 + 187 + 188 + def test_noninteger_period_samples_on_fs_grid_not_angular_grid() -> None: 189 + """Regression: 100.37 samples per period, noiseless, 12 harmonics. 190 + 191 + Each output sample ``m`` averages the input at the times 192 + ``n·T + m/fs``: evaluated on the ``m/fs`` grid the noiseless recovery 193 + error is at the fractional-delay interpolation level (~3e-6 here), 194 + while evaluating the same samples on an ``m·T/M`` angular grid -- 195 + the axis reported before the fix -- misreads them by ~9e-2. 196 + """ 197 + fs = 1000.0 198 + period = 100.37 / fs 199 + m_int = round(fs * period) 200 + n_avg = 50 201 + n_samples = int(np.ceil((n_avg + 1) * fs * period)) 202 + t = np.arange(n_samples) / fs 203 + 204 + orders = np.arange(1, 13) 205 + amps = 1.0 / orders 206 + phases = 0.821 * orders**2 207 + 208 + def wave(tv: np.ndarray) -> np.ndarray: 209 + return sum( 210 + a * np.cos(2.0 * np.pi * k / period * tv + p) 211 + for k, a, p in zip(orders, amps, phases) 212 + ) 213 + 214 + result = ph.time_synchronous_average(wave(t), fs, period, n_averages=n_avg) 215 + 216 + assert result.interpolated is True 217 + assert result.samples_per_period == m_int 218 + assert np.allclose(result.times, np.arange(m_int) / fs, atol=0.0) 219 + 220 + err_fs_grid = np.max(np.abs(result.period_waveform - wave(result.times))) 221 + assert err_fs_grid < 1e-4 # measured ~2.8e-6 222 + 223 + angular_grid = np.arange(m_int) * (period / m_int) 224 + err_angular = np.max(np.abs(result.period_waveform - wave(angular_grid))) 225 + assert err_angular > 0.05 # the angular grid is the wrong axis 180 226 181 227 182 228 # --------------------------------------------------------------------------- ··· 257 303 comb_filter_response(one, PERIOD, 0) 258 304 with pytest.raises(ValueError, match="positive"): 259 305 comb_filter_response(one, -1.0, 2) 306 + 307 + 308 + def test_comb_filter_response_rejects_overflowing_order() -> None: 309 + # Finite f and T whose product overflows would give sin(inf) = NaN 310 + # instead of the documented bounded response. 311 + with pytest.raises(ValueError, match="overflows"): 312 + comb_filter_response(np.array([1e308]), 1e308, 3) 313 + with pytest.raises(ValueError, match="overflows"): 314 + comb_filter_response(np.array([1e307]), 100.0, 2**40) 260 315 261 316 262 317 def test_plot_returns_axes() -> None:
+12 -1
src/phonometry/_plot/metrology.py
··· 218 218 "Cascade": "Cascada", 219 219 "Parametric EQ response (Audio EQ Cookbook)": 220 220 "Respuesta del EQ paramétrico (Audio EQ Cookbook)", 221 + # RBJ biquad type names (per-section legend labels); "lowpass" and 222 + # "highpass" are already translated above for the lifter plot. 223 + "peaking": "campana", 224 + "lowshelf": "shelf de graves", 225 + "highshelf": "shelf de agudos", 226 + "bandpass": "paso banda", 227 + "bandpass_skirt": "paso banda (faldón)", 228 + "notch": "muesca", 229 + "allpass": "paso todo", 221 230 "Time synchronous average (McFadden 1987)": 222 231 "Promediado síncrono en el tiempo (McFadden 1987)", 223 232 "Averaged periodic waveform (N = {n})": ··· 1831 1840 if show_sections and result.section_magnitude_db.shape[0] > 1: 1832 1841 for idx, section in enumerate(result.sections): 1833 1842 label = decimal_comma( 1834 - f"{section.filter_type} {section.f0:g} Hz", language 1843 + f"{_t(section.filter_type, language)} " 1844 + f"{section.f0:g} Hz", 1845 + language, 1835 1846 ) 1836 1847 axm.semilogx( 1837 1848 freqs, result.section_magnitude_db[idx],
+65 -21
src/phonometry/metrology/cepstrum.py
··· 13 13 at the quefrency of its period. Three variants are standard: 14 14 15 15 * the **power cepstrum**, the inverse transform of the log *power* spectrum 16 - ``ln|X|^2`` (Milner Fig. 21) -- real, even, phase-blind; 16 + ``ln|X|^2`` (Milner Fig. 21) -- real, even, phase-blind. (Convention 17 + note: Bogert, Healy & Tukey's original 1963 power cepstrum squares once 18 + more, ``|IDFT(ln|X|^2)|^2``; this module follows Milner's linear 19 + ``IDFT(ln|X|^2)`` throughout.); 17 20 * the **real cepstrum**, the inverse transform of ``ln|X|`` -- exactly half 18 21 the power cepstrum, and the quantity whose causal folding yields the 19 22 minimum-phase reconstruction of :func:`phonometry.minimum_phase` ··· 397 400 ) 398 401 399 402 403 + def _parabolic_offset( 404 + magnitude: NDArray[np.float64], peak: int, lo: int, hi: int 405 + ) -> float: 406 + """Sub-sample offset of a peak by quadratic (parabolic) interpolation. 407 + 408 + Fits a parabola through ``magnitude`` at ``peak - 1, peak, peak + 1`` 409 + and returns the vertex offset in ``[-0.5, 0.5]`` samples: the standard 410 + three-point refinement of a delay that falls between quefrency bins. 411 + Zero when the peak sits at an edge of the searched band or the triple 412 + is not concave. 413 + """ 414 + if peak <= lo or peak >= hi: 415 + return 0.0 416 + y0, y1, y2 = magnitude[peak - 1 : peak + 2] 417 + curvature = y0 - 2.0 * y1 + y2 418 + if curvature >= 0.0: 419 + return 0.0 420 + return float(np.clip(0.5 * (y0 - y2) / curvature, -0.5, 0.5)) 421 + 422 + 400 423 @dataclass(frozen=True) 401 424 class EchoDetectionResult: 402 425 """An echo delay and reflection coefficient read off the power cepstrum. 403 426 404 427 :ivar quefrencies: Quefrency axis, in seconds. 405 428 :ivar cepstrum: Power cepstrum searched. 406 - :ivar delay: Quefrency of the highest cepstral peak in the searched 407 - band: the echo delay, in seconds. 408 - :ivar delay_samples: The same delay in samples. 409 - :ivar reflection_coefficient: Height of the peak. For a single in-record 410 - echo ``x(t) = s(t) + a s(t - t0)`` the power cepstrum's first 411 - rahmonic height is exactly ``a`` (the ``n = 1`` term of the 412 - ``ln(1 + a e^{-j theta})`` series), so the height estimates the 413 - reflection coefficient directly. 429 + :ivar delay: The echo delay, in seconds: the quefrency of the largest 430 + ``|cepstrum|`` peak in the searched band, refined by quadratic 431 + (parabolic) interpolation of ``|cepstrum|`` through the peak sample 432 + and its two neighbours, so an off-sample delay is estimated to a 433 + fraction of a sample. 434 + :ivar delay_samples: The peak position in whole samples (no 435 + interpolation): the index at which ``reflection_coefficient`` 436 + is read, so ``delay`` differs from ``delay_samples/fs`` by the 437 + interpolated sub-sample offset (at most half a sample). 438 + :ivar reflection_coefficient: Signed cepstrum value at the peak sample. 439 + For a single in-record echo ``x(t) = s(t) + a s(t - t0)`` the power 440 + cepstrum's first rahmonic height is exactly ``a`` -- of either sign 441 + -- (the ``n = 1`` term of the ``ln(1 + a e^{-j theta})`` series), 442 + so the value estimates the reflection coefficient directly, 443 + including its polarity. When the true delay falls between samples 444 + the rahmonic is split across neighbouring quefrency bins and the 445 + reported coefficient underestimates ``|a|`` (down to roughly 65 % 446 + of it for a delay midway between samples); the interpolated 447 + :attr:`delay` is unaffected. 414 448 :ivar search_range: The ``(min, max)`` quefrency band searched, s. 415 449 :ivar fs: Sample rate of the analysed record, in Hz. 416 450 :ivar nfft: FFT length used. ··· 448 482 nfft: int | None = None, 449 483 ) -> EchoDetectionResult: 450 484 """ 451 - Detect an echo as the highest power-cepstrum peak in a quefrency band. 485 + Detect an echo as the largest power-cepstrum peak in a quefrency band. 452 486 453 - A reflection ``x(t) = s(t) + a s(t - t0)`` leaves a positive spike of 454 - height ``a`` at quefrency ``t0`` in the power cepstrum (module note; 455 - the seismic reverberation spike trains of Neelamani Sec. 3.3 are the 456 - same signature), regardless of the spectrum of ``s`` itself, which 457 - concentrates at low quefrencies. The search band therefore starts 458 - above the low-quefrency region occupied by the source's spectral 459 - envelope: raise ``min_quefrency`` if the source is very reverberant 460 - or narrowband. 487 + A reflection ``x(t) = s(t) + a s(t - t0)`` leaves a spike of height 488 + ``a`` -- positive or negative with the sign of the reflection -- at 489 + quefrency ``t0`` in the power cepstrum (module note; the seismic 490 + reverberation spike trains of Neelamani Sec. 3.3 are the same 491 + signature), regardless of the spectrum of ``s`` itself, which 492 + concentrates at low quefrencies. The peak is therefore picked on 493 + ``|cepstrum|``, so an inverting reflection (``a < 0``, e.g. a 494 + pressure-release boundary) is found at its true delay, and the 495 + *signed* cepstrum value at the peak is returned as the reflection 496 + coefficient. The search band starts above the low-quefrency region 497 + occupied by the source's spectral envelope: raise ``min_quefrency`` 498 + if the source is very reverberant or narrowband. 499 + 500 + The delay is refined by quadratic interpolation of ``|cepstrum|`` 501 + around the peak; see :class:`EchoDetectionResult` for the bin-splitting 502 + caveat on the coefficient at off-sample delays. 461 503 462 504 :param x: Signal (an impulse response, or any record with an in-record 463 505 echo), 1-D. ··· 486 528 f"samples of nfft = {n}." 487 529 ) 488 530 489 - band = result.cepstrum[lo : hi + 1] 490 - peak = int(np.argmax(band)) + lo 531 + # Peak-pick on the magnitude so a negative (inverting) reflection is 532 + # found at its true delay; the signed value at the peak is reported. 533 + magnitude = np.abs(result.cepstrum) 534 + peak = int(np.argmax(magnitude[lo : hi + 1])) + lo 491 535 return EchoDetectionResult( 492 536 quefrencies=result.quefrencies, 493 537 cepstrum=result.cepstrum, 494 - delay=peak / fs_v, 538 + delay=(peak + _parabolic_offset(magnitude, peak, lo, hi)) / fs_v, 495 539 delay_samples=peak, 496 540 reflection_coefficient=float(result.cepstrum[peak]), 497 541 search_range=(lo / fs_v, hi / fs_v),
+5
src/phonometry/metrology/compliance.py
··· 614 614 # column lives in tests/reference_data.py and is pinned by the CI 615 615 # conformance report.) Row = (nominal Hz, type1 upper, type1 lower, 616 616 # type2 upper, type2 lower); a -inf lower limit means upper-only. 617 + # Transcription note, 20 Hz Type 2: the standard prints a bare "+3" there, 618 + # read here as +3/upper-only (matching the surrounding upper-only rows); a 619 + # plausible alternative reading is +/-3. The B-weighting response at 20 Hz 620 + # sits only 0.05 dB below nominal, so either reading yields the same 621 + # verdict. 617 622 _ANSI_S14_TABLE5_12: list[tuple[float, float, float, float, float]] = [ 618 623 (10.0, 4.0, -4.0, 5.0, -_INF), 619 624 (12.5, 3.5, -3.5, 5.0, -_INF),
+83 -16
src/phonometry/metrology/envelope.py
··· 30 30 narrowband. 31 31 32 32 The **envelope spectrum** (:func:`envelope_spectrum`) transforms the 33 - detected envelope itself: Section 13.3 of the book runs a square-law 34 - envelope detector into a DC remover before correlating (Figure 13.11), 35 - because the spectral content of the envelope - not of the signal - is 36 - where amplitude modulations show as discrete lines. An AM tone with 37 - modulation frequency ``f_m`` and depth ``m`` puts a line of closed-form 38 - amplitude at exactly ``f_m``, the anchor the tests pin. 33 + detected envelope itself: Section 13.3 of the book runs a band-pass 34 + filter and a square-law envelope detector into a DC remover before 35 + correlating (Figure 13.11), because the spectral content of the envelope 36 + - not of the signal - is where amplitude modulations show as discrete 37 + lines. The optional ``band`` argument reproduces the figure's band-pass 38 + front end (the classical bearing-envelope chain: isolate the resonance 39 + band, then envelope it). An AM tone with modulation frequency ``f_m`` 40 + (on an analysis bin) and depth ``m`` puts a line of closed-form 41 + amplitude at exactly ``f_m``, the anchor the tests pin; off-bin 42 + modulation lines read low by the taper's scalloping loss. 39 43 """ 40 44 41 45 from __future__ import annotations ··· 200 204 :ivar frequencies: Frequency axis of the spectrum, in Hz. 201 205 :ivar amplitude: One-sided amplitude spectrum of the (mean-removed) 202 206 envelope: the height of a discrete modulation line in the units of 203 - the envelope itself. The zero-frequency bin is not doubled. 207 + the envelope itself, exact when the modulation frequency falls on 208 + an analysis bin (off-bin lines read low by the taper's scalloping 209 + loss; see :func:`envelope_spectrum`). The zero-frequency bin is 210 + not doubled. 204 211 :ivar mean_level: Mean of the detected envelope (the DC the remover of 205 212 Figure 13.11 takes out): the carrier amplitude for 206 213 ``kind="magnitude"``, its mean square for ``kind="squared"``. ··· 213 220 :ivar remove_dc: Whether the envelope mean was removed first. 214 221 :ivar fs: Sample rate of the analysed record, in Hz. 215 222 :ivar nfft: FFT length used. 223 + :ivar band: ``(low, high)`` edges of the zero-phase band-pass 224 + pre-filter applied before envelope detection, in Hz, or ``None`` 225 + (no pre-filter). 216 226 """ 217 227 218 228 frequencies: NDArray[np.float64] ··· 225 235 remove_dc: bool 226 236 fs: float 227 237 nfft: int 238 + band: tuple[float, float] | None = None 228 239 229 240 def plot( 230 241 self, ax: Axes | None = None, *, language: str = "en", **kwargs: Any ··· 242 253 return plot_envelope_spectrum(self, ax=ax, language=language, **kwargs) 243 254 244 255 256 + def _bandpass_pre_filter( 257 + xa: NDArray[np.float64], fs: float, band: tuple[float, float] 258 + ) -> tuple[NDArray[np.float64], tuple[float, float]]: 259 + """Zero-phase 4th-order Butterworth band-pass ahead of the detector. 260 + 261 + Returns the filtered record and the validated ``(low, high)`` edges. 262 + """ 263 + from scipy import signal as sp_signal 264 + 265 + try: 266 + low, high = (float(edge) for edge in band) 267 + except (TypeError, ValueError) as exc: 268 + raise ValueError( 269 + "'band' must be a pair of numeric (low, high) edges in Hz, " 270 + f"got {band!r}." 271 + ) from exc 272 + if not 0.0 < low < high < fs / 2.0: 273 + raise ValueError( 274 + "'band' must satisfy 0 < low < high < fs/2; got " 275 + f"({low:g}, {high:g}) at fs = {fs:g} Hz." 276 + ) 277 + sos = sp_signal.butter(4, (low, high), btype="bandpass", fs=fs, 278 + output="sos") 279 + # sosfiltfilt's default edge padding needs 3*(2*n_sections + 1) 280 + # samples; fail with a clear message instead of scipy's padlen error. 281 + min_length = 3 * (2 * sos.shape[0] + 1) + 1 282 + if xa.size < min_length: 283 + raise ValueError( 284 + f"'x' is too short ({xa.size} samples) for the zero-phase " 285 + f"band-pass pre-filter, which needs at least {min_length} " 286 + "samples of padding; lengthen the record or omit 'band'." 287 + ) 288 + filtered = np.asarray(sp_signal.sosfiltfilt(sos, xa), dtype=np.float64) 289 + return filtered, (low, high) 290 + 291 + 245 292 def envelope_spectrum( 246 293 x: NDArray[np.float64] | list[float], 247 294 fs: float, ··· 250 297 window: str = "hann", 251 298 nfft: int | None = None, 252 299 remove_dc: bool = True, 300 + band: tuple[float, float] | None = None, 253 301 ) -> EnvelopeSpectrumResult: 254 302 """Amplitude spectrum of the envelope: where modulations become lines. 255 303 256 304 Follows the structure of Bendat & Piersol Section 13.3 (Figure 13.11): 257 - an envelope detector, a DC remover, and a spectral view of what is 258 - left. The detector is the Hilbert envelope ``A(t) = |z(t)|`` 259 - (``kind="magnitude"``, the practical default) or the book's square-law 260 - detector ``A^2(t) = x^2 + x_hat^2`` (``kind="squared"``); its mean is 261 - removed (kept in :attr:`EnvelopeSpectrumResult.mean_level`) and the 262 - remainder is tapered and transformed once, scaled so a sinusoidal 263 - modulation reads out as a line at its exact amplitude. 305 + a band-pass filter (optional here), an envelope detector, a DC 306 + remover, and a spectral view of what is left. The detector is the 307 + Hilbert envelope ``A(t) = |z(t)|`` (``kind="magnitude"``, the 308 + practical default) or the book's square-law detector 309 + ``A^2(t) = x^2 + x_hat^2`` (``kind="squared"``); its mean is removed 310 + (kept in :attr:`EnvelopeSpectrumResult.mean_level`) and the remainder 311 + is tapered and transformed once, scaled by the taper's coherent gain 312 + so a sinusoidal modulation whose frequency falls on an analysis bin 313 + reads out as a line at its exact amplitude. An off-bin modulation 314 + frequency reads low by the taper's scalloping loss -- up to about 315 + 1.4 dB (~15 %) for the default Hann midway between bins -- like any 316 + single-record amplitude spectrum. 264 317 265 318 Closed forms for an AM tone ``A0 (1 + m cos(2 pi f_m t)) cos(2 pi f_c t)`` 266 - with ``0 <= m < 1``: 319 + with ``0 <= m < 1`` and ``f_m`` on an analysis bin: 267 320 268 321 * ``kind="magnitude"``: a line of amplitude ``A0 m`` at ``f_m``; 269 322 mean level ``A0``. ··· 273 326 Amplitude modulation of rotating machinery (bearing and gear defect 274 327 frequencies), mains hum and wind-turbine amplitude modulation appear 275 328 the same way: lines at the modulation frequency and its harmonics, 276 - separated from the carrier's own spectrum. 329 + separated from the carrier's own spectrum. For the classical bearing 330 + chain -- isolate a structural-resonance band excited by the defect 331 + impacts, then envelope it -- pass the resonance band as ``band``: the 332 + record is band-pass filtered (zero-phase, so the modulation phase is 333 + untouched) before the detector, the Figure 13.11 front end. 277 334 278 335 :param x: Signal, 1-D. 279 336 :param fs: Sample rate, in Hz. ··· 285 342 :param remove_dc: Remove the envelope mean before the transform 286 343 (default ``True``, the Figure 13.11 DC remover); the mean is 287 344 reported either way. 345 + :param band: Optional ``(low, high)`` band-pass edges, in Hz 346 + (``0 < low < high < fs/2``), applied to the record before 347 + envelope detection as a zero-phase 4th-order Butterworth 348 + (:func:`scipy.signal.sosfiltfilt`, giving an 8th-order magnitude 349 + roll-off). Default ``None``: detect on the record as given. 288 350 :return: An :class:`EnvelopeSpectrumResult`. 289 351 :raises ValueError: If the inputs or parameters are invalid. 290 352 """ ··· 302 364 raise ValueError( 303 365 f"'nfft' must be at least the record length ({n} samples)." 304 366 ) 367 + 368 + band_v: tuple[float, float] | None = None 369 + if band is not None: 370 + xa, band_v = _bandpass_pre_filter(xa, fs_v, band) 305 371 306 372 env = np.asarray(np.abs(sp_signal.hilbert(xa)), dtype=np.float64) 307 373 detector = env**2 if kind == "squared" else env ··· 326 392 remove_dc=bool(remove_dc), 327 393 fs=fs_v, 328 394 nfft=nfft_v, 395 + band=band_v, 329 396 )
+103 -5
src/phonometry/metrology/equalizer.py
··· 114 114 :ivar bw: Bandwidth in octaves (peaking: between the midpoint-gain 115 115 frequencies; band-pass/notch: between the -3 dB frequencies). 116 116 :ivar slope: Shelf-slope parameter ``S`` (shelves only; ``S = 1`` is the 117 - steepest monotonic slope). 117 + steepest monotonic slope, and ``S`` must stay below the 118 + gain-dependent bound ``(A + 1/A)/(A + 1/A - 2)`` with 119 + ``A = 10^(gain_db/40)``, beyond which the cookbook's alpha turns 120 + complex; at 0 dB gain ``A = 1``, the bound is unbounded and every 121 + positive slope is admissible). 118 122 """ 119 123 120 124 filter_type: EQFilterType ··· 135 139 f"Unknown filter_type {section.filter_type!r}; expected one of " 136 140 f"{sorted(_DESIGNERS)}." 137 141 ) 142 + for name, value in ( 143 + ("f0", section.f0), 144 + ("gain_db", section.gain_db), 145 + ("q", section.q), 146 + ("bw", section.bw), 147 + ("slope", section.slope), 148 + ): 149 + if value is not None and not math.isfinite(value): 150 + raise ValueError(f"'{name}' must be finite, got {value!r}.") 138 151 if section.f0 <= 0: 139 152 raise ValueError("Centre frequency 'f0' must be positive.") 153 + _gain_amplitude(section.gain_db) 140 154 _validate_section_parameterization(section) 155 + 156 + 157 + def _gain_amplitude(gain_db: float) -> float: 158 + """The cookbook's amplitude ``A = 10^(gain_db/40)``, representable. 159 + 160 + A finite but extreme ``gain_db`` overflows the power (about 161 + +12320 dB) or underflows it to zero (about -12320 dB), where the 162 + designers would divide by zero; both are rejected with the documented 163 + :class:`ValueError` instead of leaking raw arithmetic errors. 164 + """ 165 + try: 166 + big_a = float(10.0 ** (gain_db / 40.0)) 167 + except OverflowError: 168 + big_a = math.inf 169 + if not (big_a > 0.0 and math.isfinite(big_a)): 170 + raise ValueError( 171 + f"'gain_db' = {gain_db:g} dB is outside the representable " 172 + "range: 10^(gain_db/40) must be a positive finite number." 173 + ) 174 + return big_a 175 + 176 + 177 + def _validate_shelf_slope(section: EQSection) -> None: 178 + """Bound the shelf slope so the cookbook's alpha stays real. 179 + 180 + The alpha recipe takes ``sqrt((A + 1/A)·(1/S - 1) + 2)``, which is 181 + real only for ``S < (A + 1/A)/(A + 1/A - 2)`` (any ``S`` when the 182 + gain is 0 dB, i.e. ``A = 1``); at the bound itself alpha is 0 and the 183 + poles land on the unit circle. 184 + """ 185 + if section.slope is None: 186 + return 187 + big_a = _gain_amplitude(section.gain_db) 188 + a_sum = big_a + 1.0 / big_a 189 + if a_sum <= 2.0: # A = 1 (gain 0 dB): every positive slope is fine 190 + return 191 + bound = a_sum / (a_sum - 2.0) 192 + if section.slope >= bound: 193 + raise ValueError( 194 + f"Shelf slope 'slope' = {section.slope:g} is too steep for " 195 + f"gain_db = {section.gain_db:g}: the cookbook's alpha is real " 196 + f"only for slopes below (A + 1/A)/(A + 1/A - 2) = {bound:.6g}." 197 + ) 141 198 142 199 143 200 def _validate_section_parameterization(section: EQSection) -> None: ··· 165 222 raise ValueError( 166 223 "'gain_db' applies to 'peaking', 'lowshelf' and 'highshelf' only." 167 224 ) 225 + _validate_shelf_slope(section) 168 226 169 227 170 228 # --------------------------------------------------------------------------- ··· 259 317 sin_w0 = math.sin(w0) 260 318 if section.bw is not None: 261 319 # Digital-domain bandwidth relation: the w0/sin(w0) factor 262 - # compensates the bilinear-transform frequency warping. 263 - return sin_w0 * math.sinh(math.log(2.0) / 2 * section.bw * w0 / sin_w0) 320 + # compensates the bilinear-transform frequency warping. Near the 321 + # Nyquist frequency that factor diverges, so a wide 'bw' can push 322 + # the sinh argument past the floating-point range. 323 + try: 324 + return sin_w0 * math.sinh( 325 + math.log(2.0) / 2 * section.bw * w0 / sin_w0 326 + ) 327 + except OverflowError: 328 + raise ValueError( 329 + f"Bandwidth 'bw' = {section.bw:g} octaves is too wide for " 330 + f"f0 = {section.f0:g} Hz this close to the Nyquist " 331 + "frequency: the cookbook's bandwidth-to-alpha relation " 332 + "overflows. Reduce 'bw', lower 'f0' or parameterize the " 333 + "section with 'q' instead." 334 + ) from None 264 335 if section.slope is not None: 265 336 return (sin_w0 / 2) * math.sqrt( 266 337 (big_a + 1 / big_a) * (1 / section.slope - 1) + 2 267 338 ) 268 339 q = section.q if section.q is not None else _DEFAULT_Q 269 340 return sin_w0 / (2 * q) 341 + 342 + 343 + def _check_section_stability( 344 + fs: float, section: EQSection, a1: float, a2: float 345 + ) -> None: 346 + """Require the realized poles strictly inside the unit circle. 347 + 348 + In exact arithmetic every cookbook section with ``alpha > 0`` is 349 + stable, but at the extremes (a very wide ``bw`` or ``f0`` within a few 350 + ppm of DC or the Nyquist frequency) the rounded coefficients can land 351 + exactly on the stability-triangle boundary -- e.g. ``a2 == -1.0``, a 352 + pole pair on the unit circle -- and the filter would ring forever 353 + without any arithmetic error to flag it. 354 + """ 355 + if abs(a2) < 1.0 and abs(a1) < 1.0 + a2: 356 + return 357 + raise ValueError( 358 + f"The designed {section.filter_type!r} section at f0 = " 359 + f"{section.f0:g} Hz (fs = {fs:g} Hz) is not strictly stable after " 360 + f"rounding (a1 = {a1:.17g}, a2 = {a2:.17g}): its poles do not lie " 361 + "strictly inside the unit circle. This happens when 'bw' (or " 362 + "1/'q') is very large or f0 sits within a few ppm of DC or the " 363 + "Nyquist frequency; bring the section parameters back from the " 364 + "extreme." 365 + ) 270 366 271 367 272 368 def _section_sos(fs: float, section: EQSection) -> NDArray[np.float64]: ··· 277 373 f"Nyquist frequency {fs / 2} Hz." 278 374 ) 279 375 w0 = 2 * math.pi * section.f0 / fs 280 - big_a = 10.0 ** (section.gain_db / 40.0) 376 + big_a = _gain_amplitude(section.gain_db) 281 377 alpha = _section_alpha(section, w0, big_a) 282 378 b, a = _DESIGNERS[section.filter_type](math.cos(w0), math.sin(w0), alpha, big_a) 283 379 a0 = a[0] 380 + a1, a2 = a[1] / a0, a[2] / a0 381 + _check_section_stability(fs, section, a1, a2) 284 382 return np.array( 285 - [b[0] / a0, b[1] / a0, b[2] / a0, 1.0, a[1] / a0, a[2] / a0], 383 + [b[0] / a0, b[1] / a0, b[2] / a0, 1.0, a1, a2], 286 384 dtype=np.float64, 287 385 ) 288 386
+21 -1
src/phonometry/metrology/random_data.py
··· 153 153 correction. This reproduces every ``alpha = 0.05`` entry of B&P 154 154 Table A.6 (``N`` = 10 to 100), e.g. ``(64, 125)`` for ``N = 20`` -- 155 155 the book's Examples 4.4 and 10.3. 156 + 157 + The normal approximation is kept *deliberately*, for all ``N``: the 158 + book's table itself follows it. No percentage-point convention on the 159 + exact Mahonian null distribution reproduces every tabulated entry 160 + (checked digit by digit: e.g. at ``N = 40``, ``a = 0.975`` the exact 161 + tail jumps from ``P(A > 304) = 0.9771`` to ``P(A > 305) = 0.9758``, 162 + neither equal to 0.975, and the nearest-tail point would be 306, yet 163 + the table prints 305 -- exactly where the continuity-corrected normal 164 + point 305.39 rounds to; ``N = 100``, ``a = 0.975`` behaves the same 165 + way). Exact-distribution bounds would therefore *disagree* with the 166 + published Table A.6 that users check against. The p-value, which has 167 + no tabulated convention to honour, does use the exact distribution 168 + (:func:`_reverse_arrangement_p_value`), so verdict and p-value can 169 + disagree by at most one count right at an acceptance boundary. 156 170 """ 157 171 mean, std = _reverse_arrangement_moments(n) 158 172 spread = float(special.ndtri(1.0 - alpha / 2.0)) * std ··· 290 304 :ivar std: Null standard deviation (B&P Eq. (4.55) for ``A``). 291 305 :ivar bounds: Acceptance region ``(lower, upper)``: percentage points 292 306 such that the no-trend hypothesis is accepted when 293 - ``lower < statistic <= upper``. 307 + ``lower < statistic <= upper``. For reverse arrangements these 308 + follow B&P Table A.6's own convention (normal approximation with 309 + continuity correction, which reproduces the book's tabulated 310 + ``alpha = 0.05`` entries exactly; the table is not derivable from 311 + the exact Mahonian distribution), so at an acceptance boundary 312 + the verdict and the exact :attr:`p_value` can disagree by one 313 + count. 294 314 :ivar p_value: Two-sided p-value from the exact null distribution 295 315 (normal approximation above ``n = 100`` for reverse arrangements). 296 316 :ivar trend_free: ``True`` when the statistic falls inside the
+36 -4
src/phonometry/metrology/signals.py
··· 49 49 50 50 from __future__ import annotations 51 51 52 + import math 53 + import warnings 52 54 from dataclasses import dataclass 53 55 from typing import TYPE_CHECKING, Any, Literal 54 56 55 57 import numpy as np 56 58 59 + from .._internal.warnings import PhonometryWarning 57 60 from .spectra import _positive 58 61 59 62 if TYPE_CHECKING: ··· 295 298 repetition period, as in the repetitive-burst test of Clause A2.2 296 299 (there: 5 ms bursts of 5 kHz tone at 2, 10 or 100 bursts per second). 297 300 301 + The gate closes after ``round(fs·cycles/frequency)`` samples, so the 302 + "integral number of full periods" is sample-exact only when ``fs`` 303 + and ``frequency`` are commensurate (``fs·cycles/frequency`` an 304 + integer). Otherwise the gate closes up to half a sample away from 305 + the tone's final zero crossing and the gated waveform carries a 306 + residual step of up to ``amplitude·sin(π·frequency/fs)`` there 307 + (e.g. 10 cycles of 997 Hz at 48 kHz span 481.44 samples, gated at 308 + 481); a :class:`~phonometry.PhonometryWarning` quantifies the 309 + realized residual. 310 + 298 311 :param fs: Sample rate, in Hz. 299 312 :param frequency: Tone frequency, in Hz (below the Nyquist rate). 300 313 :param cycles: Full tone periods per burst (positive integer). ··· 316 329 ) 317 330 318 331 burst_seconds = cycles_v / f_v 319 - n_on = round(fs_v * burst_seconds) 332 + exact_samples = fs_v * burst_seconds 333 + n_on = round(exact_samples) 320 334 if n_on < 2: 321 335 raise ValueError( 322 336 "The burst is shorter than 2 samples; increase 'cycles' or 'fs'." 323 337 ) 324 - 325 338 rate_v, period, duty = _tone_burst_period( 326 339 fs_v, n_on, repetitions_v, repetition_rate 327 340 ) 341 + 342 + # Warn only once every configuration check has passed, so an invalid 343 + # repetition setup raises instead of warning first (which would mask 344 + # the error under a warnings-as-errors filter). 345 + delta = n_on - exact_samples # gate-close offset from the zero crossing 346 + if abs(delta) > 1e-9: 347 + residual = amplitude_v * math.sin(2.0 * math.pi * f_v * delta / fs_v) 348 + warnings.warn( 349 + f"'frequency' = {f_v:g} Hz is incommensurate with fs = " 350 + f"{fs_v:g} Hz: {cycles_v} full periods span {exact_samples:.6g} " 351 + f"samples, gated at {n_on}. The gate closes {abs(delta):.3g} " 352 + "samples away from the tone's final zero crossing, leaving a " 353 + f"residual step of {residual:.3g} (peak amplitude " 354 + f"{amplitude_v:g}).", 355 + PhonometryWarning, 356 + stacklevel=2, 357 + ) 328 358 329 359 n_pre = round(float(pre_silence) * fs_v) 330 360 n_post = round(float(post_silence) * fs_v) ··· 388 418 :ivar stopband_edge_hz: Stopband edge of the design (the smaller of 389 419 the two Nyquist frequencies), in Hz. 390 420 :ivar stopband_attenuation_db: Designed stopband attenuation, in dB 391 - (also the passband ripple bound: the Kaiser method is 392 - equiripple-bounded by the same ``δ = 10^(-A/20)`` in both bands). 421 + (also the passband ripple bound: the Kaiser window method holds 422 + the ripple of both bands within the same ``δ = 10^(-A/20)``, 423 + though -- unlike a true equiripple design -- its ripple decays 424 + away from the band edges rather than staying at the bound). 393 425 :ivar transition_width: Transition-band width as a fraction of the 394 426 smaller Nyquist frequency. 395 427 """
+13 -4
src/phonometry/metrology/spectra.py
··· 1318 1318 convergence): wherever the local spectrum is weak relative to the 1319 1319 broad-band leakage each taper could carry, the leakier high-order 1320 1320 tapers are downweighted, trading degrees of freedom (Eq. 370b) for 1321 - leakage protection in high-dynamic-range spectra. For a locally white 1322 - spectrum the weights converge to uniform and nothing is lost. With 1323 - ``adaptive=False`` the eigenvalue-weighted average of P&W Eq. 369a is 1324 - returned. 1321 + leakage protection in high-dynamic-range spectra. The broadband 1322 + ``σ²`` driving the weights is ``mean(x²)`` with no mean removal, 1323 + consistent with the no-detrending calibration below. For a locally 1324 + white spectrum the weights converge to uniform and nothing is lost. 1325 + With ``adaptive=False`` the eigenvalue-weighted average of P&W 1326 + Eq. 369a is returned. 1325 1327 1326 1328 Calibration matches the Welch estimators of this module exactly: no 1327 1329 detrending, ``'density'`` scaling integrates to the signal power ··· 1353 1355 scaling_v = _validate_scaling(scaling) 1354 1356 conf = _validate_confidence(confidence) 1355 1357 nw, k = _validate_multitaper_params(xa.size, time_half_bandwidth, n_tapers) 1358 + # Broadband power sigma^2 for the adaptive weights: mean(x^2) without 1359 + # mean removal, consistent with the module-wide no-detrending 1360 + # calibration (a DC offset counts as signal power here too). 1356 1361 power = float(np.mean(xa * xa)) 1357 1362 if power <= 0.0: 1358 1363 raise ValueError("'x' must not be identically zero.") ··· 1361 1366 if adaptive: 1362 1367 d = _adaptive_multitaper_weights(sk, eigenvalues, power / fs_v) 1363 1368 else: 1369 + # Eigenvalue-weighted average, cited in the docstring as P&W 1370 + # Eq. 369a; the equation number is pending verification against 1371 + # the book (source copy not yet acquired). The formula itself is 1372 + # Thomson's standard lambda_k-weighted eigenspectrum average. 1364 1373 d = eigenvalues[:, np.newaxis] * np.ones_like(sk) 1365 1374 weights = d / np.sum(d, axis=0) 1366 1375 psd = np.sum(weights * sk, axis=0)
+30 -9
src/phonometry/metrology/synchronous_average.py
··· 43 43 not removed, its rejection is optimised by choosing ``N`` so that a comb 44 44 node lands exactly on it, i.e. the smallest ``N`` with ``N·q`` an integer. 45 45 McFadden's own example, a tone at 32.05 orders, is suppressed by more than 46 - 100 dB with ``N = 20`` (since ``20·32.05 = 641``) yet only ~14 dB with the 47 - common power-of-two choice ``N = 32`` (``32·32.05 = 1025.6``). Thus the 48 - habit of taking a power-of-two number of averages is not, in general, 49 - optimal. 46 + 100 dB with ``N = 20`` (since ``20·32.05 = 641``) yet -- evaluating Eq. 8 47 + at that order -- by only about 14 dB with the common power-of-two choice 48 + ``N = 32`` (``32·32.05 = 1025.6``; the paper makes the comparison but does 49 + not print this figure). Thus the habit of taking a power-of-two number of 50 + averages is not, in general, optimal. 50 51 51 52 **Non-integer samples per period.** When ``fs·T`` is not an integer the 52 53 period boundaries fall between samples. Each block is then aligned to a ··· 54 55 :func:`phonometry.metrology.signals.fractional_delay` before averaging, so 55 56 the periodic waveform is recovered within the interpolation error of that 56 57 band-limited shift. An integer ``fs·T`` needs no interpolation and the 57 - waveform is recovered to machine precision. 58 + waveform is recovered to machine precision. The averaged samples stay on 59 + the ``1/fs`` sampling grid throughout: output sample ``m`` is the average 60 + of the input at the times ``n·T + m/fs``, so the returned time axis is 61 + ``m/fs`` and the ``M = round(fs·T)`` samples cover one period exactly when 62 + ``fs·T`` is an integer and to within half a sample otherwise. No 63 + resampling onto an ``M``-point angular grid is performed. 58 64 """ 59 65 60 66 from __future__ import annotations ··· 97 103 :param frequencies: Frequencies at which to evaluate, in Hz. 98 104 :param period: Repetition period ``T``, in seconds. 99 105 :param n_averages: Number of averaged periods ``N`` (at least 1). 100 - :return: The filter magnitude at each frequency (unitless, in [0, 1]). 106 + :return: The filter magnitude at each frequency (unitless; bounded by 1, 107 + though floating-point cancellation immediately beside a tooth can 108 + return values above 1 by a few parts in 1e9). 101 109 :raises ValueError: If the parameters are invalid. 102 110 """ 103 111 period_v = _positive(period, "period") ··· 108 116 if not np.all(np.isfinite(freqs)): 109 117 raise ValueError("'frequencies' must be finite.") 110 118 order = freqs * period_v 119 + # Finite inputs can still overflow the phase products f*T and N*pi*f*T 120 + # (sin(inf) is NaN); reject them so the bounded-response contract holds. 121 + if not np.all(np.isfinite(n * np.pi * order)): 122 + raise ValueError( 123 + "'frequencies' * 'period' (times n_averages*pi) overflows the " 124 + "floating-point range; the comb filter cannot be evaluated at " 125 + "such orders." 126 + ) 111 127 lower = np.sin(np.pi * order) 112 128 upper = np.sin(n * np.pi * order) 113 129 with np.errstate(divide="ignore", invalid="ignore"): ··· 124 140 125 141 :ivar period_waveform: The averaged periodic waveform, one period of 126 142 :attr:`samples_per_period` samples. 127 - :ivar times: Time axis of :attr:`period_waveform`, in seconds, spanning 128 - one period ``[0, T)``. 143 + :ivar times: Time axis of :attr:`period_waveform`, in seconds: the 144 + sampling grid ``m/fs``, ``m = 0 .. M-1`` (the averaged samples stay 145 + on the ``1/fs`` grid; see the module note). The axis spans one 146 + period exactly when ``fs·T`` is an integer, and to within half a 147 + sample otherwise. 129 148 :ivar residual: Input minus the periodic reconstruction, over the 130 149 analysed span (``n_averages·samples_per_period`` samples, aligned 131 150 to the integer period grid): what is left after the synchronous ··· 296 315 noise_reduction_db = 10.0 * float(np.log10(n_avg)) 297 316 298 317 comb_freqs, comb_response = _comb_grid(period_v, n_avg, n_harmonics_v) 299 - times = np.arange(m_int, dtype=np.float64) * (period_v / m_int) 318 + # The averaged samples sit on the 1/fs sampling grid (each block is 319 + # aligned back to it before averaging), not on a T/M angular grid. 320 + times = np.arange(m_int, dtype=np.float64) / fs_v 300 321 301 322 return SynchronousAverageResult( 302 323 period_waveform=period_waveform,
+7 -4
src/phonometry/metrology/time_frequency.py
··· 94 94 exact calibration of 95 95 :func:`~phonometry.metrology.spectra.power_spectral_density`: 96 96 the column mean over time reproduces the Welch spectrum bin by 97 - bin, and integrating ``'density'`` columns over frequency and 98 - summing over time recovers the record energy exactly when the 99 - squared taper overlap-adds to a constant (e.g. Hann at 75 % 100 - overlap). 97 + bin. Integrating a ``'density'`` column over frequency gives that 98 + segment's taper-weighted mean square ``Σ(x·w)²/Σw²``; summing 99 + those over time *and multiplying by the hop duration* ``hop/fs`` 100 + recovers the record energy ``Σx²/fs`` when the squared taper 101 + overlap-adds to a constant (e.g. Hann at 75 % overlap), up to the 102 + taper roll-off at the record edges (the first and last segments 103 + are under-weighted: about 1-2 % low for typical records). 101 104 :ivar time_resolution: Segment duration ``T_B = nperseg/fs``, in 102 105 seconds - the time resolution of the display. 103 106 :ivar resolution_bandwidth: Effective noise bandwidth ``Bₑ`` of the
+3
src/phonometry/room/room_ir.py
··· 734 734 (constant-envelope) interval, in dB. A time-domain swept sine has 735 735 the ideal 3.02 dB; the frequency-domain synthesis stays close to 736 736 it (Mueller & Massarani 2001, Sec. 4.3: normally below 4 dB). 737 + Their figure assumes the slight magnitude smoothing the paper 738 + applies before synthesis, which this module does not implement: 739 + typical shaped sweeps here measure about 4.2-4.3 dB. 737 740 """ 738 741 739 742 signal: np.ndarray
+31 -8
site/src/content/docs/guides/cepstrum-echoes.mdx
··· 41 41 42 42 - `'power'` (default): the inverse DFT of `ln|X|²` (Milner's Fig. 21). Real, 43 43 even and phase-blind - the workhorse for echo and harmonic-family 44 - detection; 44 + detection. This is Milner's *signed* cepstrum of the log-power spectrum; 45 + Bogert's original 1963 "power cepstrum" squares once more and is 46 + non-negative - the library follows Milner throughout, so negative 47 + rahmonics keep their sign; 45 48 - `'real'`: the inverse DFT of `ln|X|` - exactly half the power cepstrum, 46 49 and the quantity whose causal folding is the minimum-phase reconstruction 47 50 (below); ··· 84 87 85 88 so the cepstrum carries a spike train at the **rahmonics** `n·t0` with 86 89 amplitudes `a, -a²/2, a³/3, ...` (their sum is `ln(1+a)`), regardless of the 87 - spectrum of `s` itself, which concentrates at low quefrencies. On the power 88 - cepstrum the first spike's height is exactly the reflection coefficient `a` - 89 - a closed form the tests and the conformance suite pin to 1e-10. 90 - `echo_detection` automates the reading: 90 + spectrum of `s` itself, which concentrates at low quefrencies. On the signed 91 + cepstrum of the log-power spectrum (`kind='power'`, Milner's convention) the 92 + first spike's height is the reflection coefficient `a` - of either 93 + sign - plus whatever the source cepstrum contributes at that quefrency 94 + (negligible for broadband sources, whose cepstrum concentrates at low 95 + quefrencies); on an ideal impulse-plus-echo the identity is a closed form 96 + the tests and the conformance suite pin to 1e-10. `echo_detection` automates the reading: it picks the largest 97 + `|cepstrum|` peak in the searched band (so an inverting reflection, `a < 0`, 98 + is found at its true delay rather than missed), refines the delay by 99 + quadratic interpolation through the peak and its neighbours, and reports the 100 + signed peak value as the reflection coefficient. When the true delay falls 101 + between samples the rahmonic splits across quefrency bins: the interpolated 102 + delay still lands on it, but the reported coefficient underestimates `|a|` 103 + (down to about 65 % of it midway between samples). 91 104 92 105 ```python 93 106 import numpy as np ··· 217 230 DC remover, and a spectral view of what remains. `envelope_spectrum` runs the 218 231 [Hilbert envelope](/phonometry/guides/correlation-delay/) (`kind="magnitude"`, 219 232 the practical default) or the book's square-law detector (`kind="squared"`) 220 - through exactly that chain, scaled so a sinusoidal modulation reads out as a 221 - line at its exact amplitude. 233 + through exactly that chain, scaled by the taper's coherent gain so a 234 + sinusoidal modulation whose frequency falls on an analysis bin reads out as 235 + a line at its exact amplitude (off-bin lines read low by the taper's 236 + scalloping loss, up to about 1.4 dB for the default Hann). The optional 237 + `band=(low, high)` argument reproduces the figure's band-pass front end - 238 + the classical bearing-envelope chain: isolate the structural-resonance band 239 + excited by the defect impacts with a zero-phase band-pass, then envelope 240 + it - so an out-of-band interferer is strongly attenuated before the 241 + detector (the roll-off of a fourth-order Butterworth applied forward and 242 + backward; the rejection is finite, and the zero-phase pass leaves small 243 + transients at the record edges). 222 244 223 - For an AM tone `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` the closed forms are: 245 + For an AM tone `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` with `fm` on an 246 + analysis bin the closed forms are: 224 247 225 248 | `kind` | mean level | line at `fm` | line at `2fm` | 226 249 |---|---|---|---|
+36 -8
site/src/content/docs/es/guides/cepstrum-echoes.mdx
··· 43 43 44 44 - `'power'` (por defecto): la DFT inversa de `ln|X|²` (Fig. 21 de Milner). 45 45 Real, par y ciega a la fase - el caballo de batalla para detectar ecos y 46 - familias de armónicos; 46 + familias de armónicos. Es el cepstro *con signo* del espectro logarítmico 47 + de potencia según Milner; el "cepstro de potencia" original de Bogert 48 + (1963) eleva al cuadrado una vez más y es no negativo - la biblioteca 49 + sigue a Milner en todo, así que los rahmónicos negativos conservan su 50 + signo; 47 51 - `'real'`: la DFT inversa de `ln|X|` - exactamente la mitad del cepstro de 48 52 potencia, y la cantidad cuyo plegado causal es la reconstrucción de fase 49 53 mínima (más abajo); ··· 88 92 así que el cepstro lleva un tren de picos en los **rahmónicos** `n·t0` con 89 93 amplitudes `a, -a²/2, a³/3, ...` (su suma es `ln(1+a)`), independientemente 90 94 del espectro de la propia `s`, que se concentra en las quefrencias bajas. En 91 - el cepstro de potencia la altura del primer pico es exactamente el 92 - coeficiente de reflexión `a` - una forma cerrada que los tests y la batería 93 - de conformidad fijan a 1e-10. `echo_detection` automatiza la lectura: 95 + el cepstro con signo del espectro logarítmico de potencia (`kind='power'`, 96 + la convención de Milner) la altura del primer pico es el coeficiente de 97 + reflexión `a` - con su signo - más lo que el cepstro de la fuente aporte en 98 + esa quefrencia (despreciable para fuentes de banda ancha, cuyo cepstro se 99 + concentra en las quefrencias bajas); sobre un impulso ideal con eco la 100 + identidad es una forma cerrada que los tests y la batería de conformidad 101 + fijan a 1e-10. `echo_detection` 102 + automatiza la lectura: busca el mayor pico de `|cepstrum|` en la banda de 103 + búsqueda (así una reflexión inversora, `a < 0`, se encuentra en su retardo 104 + verdadero en vez de perderse), refina el retardo por interpolación 105 + cuadrática a través del pico y sus vecinos, y devuelve como coeficiente de 106 + reflexión el valor del pico con su signo. Cuando el retardo verdadero cae 107 + entre muestras, el rahmónico se reparte entre bins de quefrencia vecinos: 108 + el retardo interpolado sigue acertando, pero el coeficiente devuelto 109 + subestima `|a|` (hasta cerca del 65 % de su valor a mitad de camino entre 110 + muestras). 94 111 95 112 ```python 96 113 import numpy as np ··· 225 242 `envelope_spectrum` pasa la 226 243 [envolvente de Hilbert](/phonometry/es/guides/correlation-delay/) 227 244 (`kind="magnitude"`, el valor práctico por defecto) o el detector cuadrático 228 - del libro (`kind="squared"`) por exactamente esa cadena, escalada para que 229 - una modulación sinusoidal se lea como una línea con su amplitud exacta. 245 + del libro (`kind="squared"`) por exactamente esa cadena, escalada por la 246 + ganancia coherente de la ventana para que una modulación sinusoidal cuya 247 + frecuencia cae en un bin de análisis se lea como una línea con su amplitud 248 + exacta (las líneas fuera de bin leen de menos por la pérdida de scalloping 249 + de la ventana, hasta cerca de 1,4 dB con la Hann por defecto). El argumento 250 + opcional `band=(low, high)` reproduce el filtro paso banda de entrada de la 251 + figura - la cadena clásica de envolvente para rodamientos: aislar con un 252 + paso banda de fase cero la banda de resonancia estructural que excitan los 253 + impactos del defecto y después obtener su envolvente - de modo que un 254 + interferente fuera de banda llega muy atenuado al detector (la caída de un 255 + Butterworth de cuarto orden aplicado hacia delante y hacia atrás; el 256 + rechazo es finito, y el paso de fase cero deja pequeños transitorios en 257 + los bordes del registro). 230 258 231 - Para un tono AM `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` las formas cerradas 232 - son: 259 + Para un tono AM `A0·(1 + m·cos(2πfm·t))·cos(2πfc·t)` con `fm` en un bin de 260 + análisis las formas cerradas son: 233 261 234 262 | `kind` | nivel medio | línea en `fm` | línea en `2fm` | 235 263 |---|---|---|---|
+33 -16
site/src/content/docs/reference/api/correlation/envelope.md
··· 37 37 narrowband. 38 38 39 39 The **envelope spectrum** ([`envelope_spectrum`](/phonometry/reference/api/correlation/envelope/#envelope_spectrum)) transforms the 40 - detected envelope itself: Section 13.3 of the book runs a square-law 41 - envelope detector into a DC remover before correlating (Figure 13.11), 42 - because the spectral content of the envelope - not of the signal - is 43 - where amplitude modulations show as discrete lines. An AM tone with 44 - modulation frequency `f_m` and depth `m` puts a line of closed-form 45 - amplitude at exactly `f_m`, the anchor the tests pin. 40 + detected envelope itself: Section 13.3 of the book runs a band-pass 41 + filter and a square-law envelope detector into a DC remover before 42 + correlating (Figure 13.11), because the spectral content of the envelope 43 + - not of the signal - is where amplitude modulations show as discrete 44 + lines. The optional `band` argument reproduces the figure's band-pass 45 + front end (the classical bearing-envelope chain: isolate the resonance 46 + band, then envelope it). An AM tone with modulation frequency `f_m` 47 + (on an analysis bin) and depth `m` puts a line of closed-form 48 + amplitude at exactly `f_m`, the anchor the tests pin; off-bin 49 + modulation lines read low by the taper's scalloping loss. 46 50 47 51 ## envelope 48 52 ··· 104 108 window: str = 'hann', 105 109 nfft: int | None = None, 106 110 remove_dc: bool = True, 111 + band: tuple[float, float] | None = None, 107 112 ) -> EnvelopeSpectrumResult 108 113 ``` 109 114 110 115 Amplitude spectrum of the envelope: where modulations become lines. 111 116 112 117 Follows the structure of Bendat & Piersol Section 13.3 (Figure 13.11): 113 - an envelope detector, a DC remover, and a spectral view of what is 114 - left. The detector is the Hilbert envelope `A(t) = |z(t)|` 115 - (`kind="magnitude"`, the practical default) or the book's square-law 116 - detector `A^2(t) = x^2 + x_hat^2` (`kind="squared"`); its mean is 117 - removed (kept in [`EnvelopeSpectrumResult.mean_level`](/phonometry/reference/api/correlation/envelope/#envelopespectrumresult)) and the 118 - remainder is tapered and transformed once, scaled so a sinusoidal 119 - modulation reads out as a line at its exact amplitude. 118 + a band-pass filter (optional here), an envelope detector, a DC 119 + remover, and a spectral view of what is left. The detector is the 120 + Hilbert envelope `A(t) = |z(t)|` (`kind="magnitude"`, the 121 + practical default) or the book's square-law detector 122 + `A^2(t) = x^2 + x_hat^2` (`kind="squared"`); its mean is removed 123 + (kept in [`EnvelopeSpectrumResult.mean_level`](/phonometry/reference/api/correlation/envelope/#envelopespectrumresult)) and the remainder 124 + is tapered and transformed once, scaled by the taper's coherent gain 125 + so a sinusoidal modulation whose frequency falls on an analysis bin 126 + reads out as a line at its exact amplitude. An off-bin modulation 127 + frequency reads low by the taper's scalloping loss -- up to about 128 + 1.4 dB (~15 %) for the default Hann midway between bins -- like any 129 + single-record amplitude spectrum. 120 130 121 131 Closed forms for an AM tone `A0 (1 + m cos(2 pi f_m t)) cos(2 pi f_c t)` 122 - with `0 <= m < 1`: 132 + with `0 <= m < 1` and `f_m` on an analysis bin: 123 133 124 134 * `kind="magnitude"`: a line of amplitude `A0 m` at `f_m`; 125 135 mean level `A0`. ··· 129 139 Amplitude modulation of rotating machinery (bearing and gear defect 130 140 frequencies), mains hum and wind-turbine amplitude modulation appear 131 141 the same way: lines at the modulation frequency and its harmonics, 132 - separated from the carrier's own spectrum. 142 + separated from the carrier's own spectrum. For the classical bearing 143 + chain -- isolate a structural-resonance band excited by the defect 144 + impacts, then envelope it -- pass the resonance band as `band`: the 145 + record is band-pass filtered (zero-phase, so the modulation phase is 146 + untouched) before the detector, the Figure 13.11 front end. 133 147 134 148 **Parameters** 135 149 ··· 141 155 | `window` | Taper (any scipy window name; default Hann). The amplitude is corrected for the taper's coherent gain. | 142 156 | `nfft` | FFT length, at least `x.size` (default: the record length). | 143 157 | `remove_dc` | Remove the envelope mean before the transform (default `True`, the Figure 13.11 DC remover); the mean is reported either way. | 158 + | `band` | Optional `(low, high)` band-pass edges, in Hz (`0 < low < high < fs/2`), applied to the record before envelope detection as a zero-phase 4th-order Butterworth (`scipy.signal.sosfiltfilt`, giving an 8th-order magnitude roll-off). Default `None`: detect on the record as given. | 144 159 145 160 **Returns:** An [`EnvelopeSpectrumResult`](/phonometry/reference/api/correlation/envelope/#envelopespectrumresult). 146 161 ··· 218 233 remove_dc: bool, 219 234 fs: float, 220 235 nfft: int, 236 + band: tuple[float, float] | None = None, 221 237 ) 222 238 ``` 223 239 ··· 228 244 | Name | Description | 229 245 | :--- | :--- | 230 246 | `frequencies` | Frequency axis of the spectrum, in Hz. | 231 - | `amplitude` | One-sided amplitude spectrum of the (mean-removed) envelope: the height of a discrete modulation line in the units of the envelope itself. The zero-frequency bin is not doubled. | 247 + | `amplitude` | One-sided amplitude spectrum of the (mean-removed) envelope: the height of a discrete modulation line in the units of the envelope itself, exact when the modulation frequency falls on an analysis bin (off-bin lines read low by the taper's scalloping loss; see [`envelope_spectrum`](/phonometry/reference/api/correlation/envelope/#envelope_spectrum)). The zero-frequency bin is not doubled. | 232 248 | `mean_level` | Mean of the detected envelope (the DC the remover of Figure 13.11 takes out): the carrier amplitude for `kind="magnitude"`, its mean square for `kind="squared"`. | 233 249 | `kind` | `"magnitude"` (Hilbert envelope `A(t)`) or `"squared"` (the book's square-law detector, `A^2(t)`). | 234 250 | `times` | Time axis of [`envelope`](/phonometry/reference/api/correlation/envelope/#envelope), in seconds. | ··· 237 253 | `remove_dc` | Whether the envelope mean was removed first. | 238 254 | `fs` | Sample rate of the analysed record, in Hz. | 239 255 | `nfft` | FFT length used. | 256 + | `band` | `(low, high)` edges of the zero-phase band-pass pre-filter applied before envelope detection, in Hz, or `None` (no pre-filter). | 240 257 241 258 ### EnvelopeSpectrumResult.plot() 242 259
+1 -1
site/src/content/docs/reference/api/filters/equalizer.md
··· 137 137 | `gain_db` | Gain `G` in dB - peaking and shelving types only. | 138 138 | `q` | Quality factor. Exactly one of `q`, `bw` and `slope` may be given; with none, `q = 1/sqrt(2)` (Butterworth alignment). | 139 139 | `bw` | Bandwidth in octaves (peaking: between the midpoint-gain frequencies; band-pass/notch: between the -3 dB frequencies). | 140 - | `slope` | Shelf-slope parameter `S` (shelves only; `S = 1` is the steepest monotonic slope). | 140 + | `slope` | Shelf-slope parameter `S` (shelves only; `S = 1` is the steepest monotonic slope, and `S` must stay below the gain-dependent bound `(A + 1/A)/(A + 1/A - 2)` with `A = 10^(gain_db/40)`, beyond which the cookbook's alpha turns complex; at 0 dB gain `A = 1`, the bound is unbounded and every positive slope is admissible). | 141 141 142 142 ## parametric_eq 143 143
+1 -1
site/src/content/docs/reference/api/metrology/random-data.md
··· 520 520 | `n` | Number of observations used. | 521 521 | `mean` | Null mean of the statistic (B&P Eq. (4.54) for `A`). | 522 522 | `std` | Null standard deviation (B&P Eq. (4.55) for `A`). | 523 - | `bounds` | Acceptance region `(lower, upper)`: percentage points such that the no-trend hypothesis is accepted when `lower < statistic <= upper`. | 523 + | `bounds` | Acceptance region `(lower, upper)`: percentage points such that the no-trend hypothesis is accepted when `lower < statistic <= upper`. For reverse arrangements these follow B&P Table A.6's own convention (normal approximation with continuity correction, which reproduces the book's tabulated `alpha = 0.05` entries exactly; the table is not derivable from the exact Mahonian distribution), so at an acceptance boundary the verdict and the exact `p_value` can disagree by one count. | 524 524 | `p_value` | Two-sided p-value from the exact null distribution (normal approximation above `n = 100` for reverse arrangements). | 525 525 | `trend_free` | `True` when the statistic falls inside the acceptance region. | 526 526 | `alpha` | Significance level of the region (default 0.05). |
+1 -1
site/src/content/docs/reference/api/rooms/room-ir.md
··· 475 475 | `magnitude` | Band-limited magnitude imposed on the synthesis spectrum, normalised to a peak of 1 (linear). | 476 476 | `group_delay` | Synthesized group delay `tau_G(f)` on `frequencies`, in seconds: the time at which each frequency is swept through. | 477 477 | `f_range` | `(f1, f2)` band covered by the sweep, in Hz. | 478 - | `crest_factor_db` | Peak-to-RMS ratio over the sweep's central (constant-envelope) interval, in dB. A time-domain swept sine has the ideal 3.02 dB; the frequency-domain synthesis stays close to it (Mueller & Massarani 2001, Sec. 4.3: normally below 4 dB). | 478 + | `crest_factor_db` | Peak-to-RMS ratio over the sweep's central (constant-envelope) interval, in dB. A time-domain swept sine has the ideal 3.02 dB; the frequency-domain synthesis stays close to it (Mueller & Massarani 2001, Sec. 4.3: normally below 4 dB). Their figure assumes the slight magnitude smoothing the paper applies before synthesis, which this module does not implement: typical shaped sweeps here measure about 4.2-4.3 dB. | 479 479 480 480 ### ShapedSweepResult.plot() 481 481
+24 -13
site/src/content/docs/reference/api/spectra/cepstrum.md
··· 20 20 at the quefrency of its period. Three variants are standard: 21 21 22 22 * the **power cepstrum**, the inverse transform of the log *power* spectrum 23 - `ln|X|^2` (Milner Fig. 21) -- real, even, phase-blind; 23 + `ln|X|^2` (Milner Fig. 21) -- real, even, phase-blind. (Convention 24 + note: Bogert, Healy & Tukey's original 1963 power cepstrum squares once 25 + more, `|IDFT(ln|X|^2)|^2`; this module follows Milner's linear 26 + `IDFT(ln|X|^2)` throughout.); 24 27 * the **real cepstrum**, the inverse transform of `ln|X|` -- exactly half 25 28 the power cepstrum, and the quantity whose causal folding yields the 26 29 minimum-phase reconstruction of [`phonometry.minimum_phase`](/phonometry/reference/api/spectra/phase/#minimum_phase) ··· 190 193 ) -> EchoDetectionResult 191 194 ``` 192 195 193 - Detect an echo as the highest power-cepstrum peak in a quefrency band. 196 + Detect an echo as the largest power-cepstrum peak in a quefrency band. 194 197 195 - A reflection `x(t) = s(t) + a s(t - t0)` leaves a positive spike of 196 - height `a` at quefrency `t0` in the power cepstrum (module note; 197 - the seismic reverberation spike trains of Neelamani Sec. 3.3 are the 198 - same signature), regardless of the spectrum of `s` itself, which 199 - concentrates at low quefrencies. The search band therefore starts 200 - above the low-quefrency region occupied by the source's spectral 201 - envelope: raise `min_quefrency` if the source is very reverberant 202 - or narrowband. 198 + A reflection `x(t) = s(t) + a s(t - t0)` leaves a spike of height 199 + `a` -- positive or negative with the sign of the reflection -- at 200 + quefrency `t0` in the power cepstrum (module note; the seismic 201 + reverberation spike trains of Neelamani Sec. 3.3 are the same 202 + signature), regardless of the spectrum of `s` itself, which 203 + concentrates at low quefrencies. The peak is therefore picked on 204 + `|cepstrum|`, so an inverting reflection (`a < 0`, e.g. a 205 + pressure-release boundary) is found at its true delay, and the 206 + *signed* cepstrum value at the peak is returned as the reflection 207 + coefficient. The search band starts above the low-quefrency region 208 + occupied by the source's spectral envelope: raise `min_quefrency` 209 + if the source is very reverberant or narrowband. 210 + 211 + The delay is refined by quadratic interpolation of `|cepstrum|` 212 + around the peak; see [`EchoDetectionResult`](/phonometry/reference/api/spectra/cepstrum/#echodetectionresult) for the bin-splitting 213 + caveat on the coefficient at off-sample delays. 203 214 204 215 **Parameters** 205 216 ··· 242 253 | :--- | :--- | 243 254 | `quefrencies` | Quefrency axis, in seconds. | 244 255 | `cepstrum` | Power cepstrum searched. | 245 - | `delay` | Quefrency of the highest cepstral peak in the searched band: the echo delay, in seconds. | 246 - | `delay_samples` | The same delay in samples. | 247 - | `reflection_coefficient` | Height of the peak. For a single in-record echo `x(t) = s(t) + a s(t - t0)` the power cepstrum's first rahmonic height is exactly `a` (the `n = 1` term of the `ln(1 + a e^{-j theta})` series), so the height estimates the reflection coefficient directly. | 256 + | `delay` | The echo delay, in seconds: the quefrency of the largest `\|cepstrum\|` peak in the searched band, refined by quadratic (parabolic) interpolation of `\|cepstrum\|` through the peak sample and its two neighbours, so an off-sample delay is estimated to a fraction of a sample. | 257 + | `delay_samples` | The peak position in whole samples (no interpolation): the index at which `reflection_coefficient` is read, so `delay` differs from `delay_samples/fs` by the interpolated sub-sample offset (at most half a sample). | 258 + | `reflection_coefficient` | Signed cepstrum value at the peak sample. For a single in-record echo `x(t) = s(t) + a s(t - t0)` the power cepstrum's first rahmonic height is exactly `a` -- of either sign -- (the `n = 1` term of the `ln(1 + a e^{-j theta})` series), so the value estimates the reflection coefficient directly, including its polarity. When the true delay falls between samples the rahmonic is split across neighbouring quefrency bins and the reported coefficient underestimates `\|a\|` (down to roughly 65 % of it for a delay midway between samples); the interpolated `delay` is unaffected. | 248 259 | `search_range` | The `(min, max)` quefrency band searched, s. | 249 260 | `fs` | Sample rate of the analysed record, in Hz. | 250 261 | `nfft` | FFT length used. |
+11 -1
site/src/content/docs/reference/api/spectra/signals.md
··· 223 223 | `filter_taps` | Anti-alias FIR taps (unit passband gain; the polyphase engine applies the `up` interpolation gain), designed at the intermediate rate `original_fs·up`. A single `1.0` tap when the ratio is 1 (no filtering). | 224 224 | `passband_edge_hz` | Passband edge of the design, in Hz. | 225 225 | `stopband_edge_hz` | Stopband edge of the design (the smaller of the two Nyquist frequencies), in Hz. | 226 - | `stopband_attenuation_db` | Designed stopband attenuation, in dB (also the passband ripple bound: the Kaiser method is equiripple-bounded by the same `δ = 10^(-A/20)` in both bands). | 226 + | `stopband_attenuation_db` | Designed stopband attenuation, in dB (also the passband ripple bound: the Kaiser window method holds the ripple of both bands within the same `δ = 10^(-A/20)`, though -- unlike a true equiripple design -- its ripple decays away from the band edges rather than staying at the bound). | 227 227 | `transition_width` | Transition-band width as a fraction of the smaller Nyquist frequency. | 228 228 229 229 ### ResampledSignalResult.n_taps ··· 257 257 train of `repetitions` identical bursts is produced, one per 258 258 repetition period, as in the repetitive-burst test of Clause A2.2 259 259 (there: 5 ms bursts of 5 kHz tone at 2, 10 or 100 bursts per second). 260 + 261 + The gate closes after `round(fs·cycles/frequency)` samples, so the 262 + "integral number of full periods" is sample-exact only when `fs` 263 + and `frequency` are commensurate (`fs·cycles/frequency` an 264 + integer). Otherwise the gate closes up to half a sample away from 265 + the tone's final zero crossing and the gated waveform carries a 266 + residual step of up to `amplitude·sin(π·frequency/fs)` there 267 + (e.g. 10 cycles of 997 Hz at 48 kHz span 481.44 samples, gated at 268 + 481); a [`PhonometryWarning`](/phonometry/reference/api/filters/phonometry/#phonometrywarning) quantifies the 269 + realized residual. 260 270 261 271 **Parameters** 262 272
+6 -4
site/src/content/docs/reference/api/spectra/spectra.md
··· 369 369 convergence): wherever the local spectrum is weak relative to the 370 370 broad-band leakage each taper could carry, the leakier high-order 371 371 tapers are downweighted, trading degrees of freedom (Eq. 370b) for 372 - leakage protection in high-dynamic-range spectra. For a locally white 373 - spectrum the weights converge to uniform and nothing is lost. With 374 - `adaptive=False` the eigenvalue-weighted average of P&W Eq. 369a is 375 - returned. 372 + leakage protection in high-dynamic-range spectra. The broadband 373 + `σ²` driving the weights is `mean(x²)` with no mean removal, 374 + consistent with the no-detrending calibration below. For a locally 375 + white spectrum the weights converge to uniform and nothing is lost. 376 + With `adaptive=False` the eigenvalue-weighted average of P&W 377 + Eq. 369a is returned. 376 378 377 379 Calibration matches the Welch estimators of this module exactly: no 378 380 detrending, `'density'` scaling integrates to the signal power
+13 -7
site/src/content/docs/reference/api/spectra/synchronous-average.md
··· 50 50 not removed, its rejection is optimised by choosing `N` so that a comb 51 51 node lands exactly on it, i.e. the smallest `N` with `N·q` an integer. 52 52 McFadden's own example, a tone at 32.05 orders, is suppressed by more than 53 - 100 dB with `N = 20` (since `20·32.05 = 641`) yet only ~14 dB with the 54 - common power-of-two choice `N = 32` (`32·32.05 = 1025.6`). Thus the 55 - habit of taking a power-of-two number of averages is not, in general, 56 - optimal. 53 + 100 dB with `N = 20` (since `20·32.05 = 641`) yet -- evaluating Eq. 8 54 + at that order -- by only about 14 dB with the common power-of-two choice 55 + `N = 32` (`32·32.05 = 1025.6`; the paper makes the comparison but does 56 + not print this figure). Thus the habit of taking a power-of-two number of 57 + averages is not, in general, optimal. 57 58 58 59 **Non-integer samples per period.** When `fs·T` is not an integer the 59 60 period boundaries fall between samples. Each block is then aligned to a ··· 61 62 [`phonometry.metrology.signals.fractional_delay`](/phonometry/reference/api/spectra/signals/#fractional_delay) before averaging, so 62 63 the periodic waveform is recovered within the interpolation error of that 63 64 band-limited shift. An integer `fs·T` needs no interpolation and the 64 - waveform is recovered to machine precision. 65 + waveform is recovered to machine precision. The averaged samples stay on 66 + the `1/fs` sampling grid throughout: output sample `m` is the average 67 + of the input at the times `n·T + m/fs`, so the returned time axis is 68 + `m/fs` and the `M = round(fs·T)` samples cover one period exactly when 69 + `fs·T` is an integer and to within half a sample otherwise. No 70 + resampling onto an `M`-point angular grid is performed. 65 71 66 72 ## comb_filter_response 67 73 ··· 88 94 | `period` | Repetition period `T`, in seconds. | 89 95 | `n_averages` | Number of averaged periods `N` (at least 1). | 90 96 91 - **Returns:** The filter magnitude at each frequency (unitless, in [0, 1]). 97 + **Returns:** The filter magnitude at each frequency (unitless; bounded by 1, though floating-point cancellation immediately beside a tooth can return values above 1 by a few parts in 1e9). 92 98 93 99 **Raises** 94 100 ··· 122 128 | Name | Description | 123 129 | :--- | :--- | 124 130 | `period_waveform` | The averaged periodic waveform, one period of `samples_per_period` samples. | 125 - | `times` | Time axis of `period_waveform`, in seconds, spanning one period `[0, T)`. | 131 + | `times` | Time axis of `period_waveform`, in seconds: the sampling grid `m/fs`, `m = 0 .. M-1` (the averaged samples stay on the `1/fs` grid; see the module note). The axis spans one period exactly when `fs·T` is an integer, and to within half a sample otherwise. | 126 132 | `residual` | Input minus the periodic reconstruction, over the analysed span (`n_averages·samples_per_period` samples, aligned to the integer period grid): what is left after the synchronous component is removed. | 127 133 | `n_averages` | Number of periods averaged, `N`. | 128 134 | `samples_per_period` | Integer samples per period `M` after any alignment. |
+1 -1
site/src/content/docs/reference/api/spectra/time-frequency.md
··· 130 130 | :--- | :--- | 131 131 | `times` | Segment-centre times, in seconds (one per column). | 132 132 | `frequencies` | One-sided frequency axis, in Hz (one per row). | 133 - | `power` | Power spectrogram, shape `(frequencies, times)` (units²/Hz for `'density'` scaling, units² for `'spectrum'`). Each column is the tapered periodogram of one segment, with the exact calibration of [`power_spectral_density`](/phonometry/reference/api/spectra/spectra/#power_spectral_density): the column mean over time reproduces the Welch spectrum bin by bin, and integrating `'density'` columns over frequency and summing over time recovers the record energy exactly when the squared taper overlap-adds to a constant (e.g. Hann at 75 % overlap). | 133 + | `power` | Power spectrogram, shape `(frequencies, times)` (units²/Hz for `'density'` scaling, units² for `'spectrum'`). Each column is the tapered periodogram of one segment, with the exact calibration of [`power_spectral_density`](/phonometry/reference/api/spectra/spectra/#power_spectral_density): the column mean over time reproduces the Welch spectrum bin by bin. Integrating a `'density'` column over frequency gives that segment's taper-weighted mean square `Σ(x·w)²/Σw²`; summing those over time *and multiplying by the hop duration* `hop/fs` recovers the record energy `Σx²/fs` when the squared taper overlap-adds to a constant (e.g. Hann at 75 % overlap), up to the taper roll-off at the record edges (the first and last segments are under-weighted: about 1-2 % low for typical records). | 134 134 | `time_resolution` | Segment duration `T_B = nperseg/fs`, in seconds - the time resolution of the display. | 135 135 | `resolution_bandwidth` | Effective noise bandwidth `Bₑ` of the tapered segment, in Hz - the frequency resolution (`≈ 1/T_B` for a light taper; the `BₑT_B` product per cell is close to 1). | 136 136 | `random_error` | Normalized random error of each (unaveraged) power cell for random data, `1/√nd = 1` with `nd = 1` (Eq. 8.158); Bendat & Piersol quote `√2/1.25 ≈ 1.13` for the magnitude display (Section 12.6.4.2). Deterministic components are unaffected. |