feat: reactive silencers, HVAC attenuation, machine enclosures and the radiating piston

CHANGELOG.md

@@ -67,6 +67,39 @@ and this project adheres to [Semantic Versioning](https://semver.org/).
   `W = 0.5 |F|^2 Re{Y}` (Cremer, Heckl & Petersson 2005, *Structure-Borne Sound*
   3e, Table 5.1), the theoretical companions of the measured ISO 7626 mobilities.
 
+- Industrial noise control: silencers, HVAC and machine enclosures, a new
+  `phonometry.noise_control` domain (Bies, Hansen & Howard, *Engineering Noise
+  Control* 5th ed.; Munjal, *Acoustics of Ducts and Mufflers*). `expansion_chamber`,
+  `helmholtz_resonator`, `quarter_wave_resonator` and `extended_tube_chamber`
+  build **reactive silencers** by the one-dimensional four-pole
+  (transmission-matrix) method, returning a `ReactiveSilencerResult` with the
+  transmission loss, the insertion loss for configurable source/radiation
+  impedances, the compound transfer matrix and `.plot()`; the expansion-chamber
+  transmission loss matches the closed form `10 lg[1 + (1/4)(m - 1/m)^2 sin^2 kL]`
+  (Bies Eq. (8.111)) exactly, and the four-pole machinery is available under
+  `phonometry.noise_control` (`duct_matrix`, `shunt_matrix`, `cascade`,
+  `transmission_loss`, `insertion_loss`, `helmholtz_impedance`,
+  `quarter_wave_impedance`). The
+  `noise_control.hvac` module adds the duct end reflection (`end_reflection_loss`,
+  ASHRAE Table 8.14), bend insertion loss (`elbow_insertion_loss`, Table 8.11),
+  the plenum transmission loss by Wells' method (`plenum_attenuation`) and the
+  flow-generated (self) noise of straight ducts and mitred bends
+  (`flow_noise_straight_duct`, `flow_noise_bend`, VDI 2081). `enclosure_insertion_loss`
+  gives the machine-enclosure insertion loss `IL = R - C` (Bies Eqs. (7.103),
+  (7.111)) from a **caller-supplied** panel transmission loss `R` (array or
+  callable; never predicted here) and the interior room constant reused from
+  `phonometry.room.room_constant`. The four-pole expansion chamber is
+  cross-checked against the independent 2D FDTD solver.
+- Radiation of a rigid circular piston in an infinite baffle
+  (`phonometry.electroacoustics.piston`, re-exported at the top level). `radiating_piston`
+  computes the radiation impedance `rho c S (R1 + j X1)` from the piston
+  resistance `R1 = 1 - 2 J1(2ka)/2ka` and reactance `X1 = 2 H1(2ka)/2ka`
+  (Beranek & Mellow, *Acoustics* 2nd ed., Eqs. (13.117), (13.118)), the
+  low-frequency radiation mass `8 rho a^3 / 3` (Eq. (4.151)), the far-field
+  directivity `2 J1(ka sin theta)/(ka sin theta)` and the directivity index,
+  returning a `RadiatingPistonResult` with `.plot()`; the building blocks
+  `piston_resistance`, `piston_reactance` and `piston_directivity` are exposed
+  directly.
 - Objective speech intelligibility STOI and ESTOI
   (`phonometry.hearing.objective_intelligibility.stoi`, also `phonometry.stoi`).
   `stoi(clean, degraded, fs)` computes the short-time objective intelligibility

docs/CONFORMANCE.md

@@ -15,7 +15,7 @@
 
 ## Numerical conformance report
 
-✅ **338/338 conformance checks pass** across 42 domains and 210 standards - filters class 1 - weightings within IEC 61672-1 class 1.
+✅ **353/353 conformance checks pass** across 44 domains and 224 standards - filters class 1 - weightings within IEC 61672-1 class 1.
 
 <sub>Each row pins a standard clause to its expected normative value and the value the library computes. Every section below is collapsible and stays collapsed while all of its rows pass; a section with any failing row opens automatically.</sub>
 
@@ -719,3 +719,34 @@ Only **Butterworth** (the library default) and **Chebyshev-II** are class-compli
 
 </details>
 
+<details>
+<summary>&#9989; <b>Electroacoustics</b>: 100% (6/6)</summary>
+
+| Standard | Quantity | Expected (norm) | Computed | &#916; | Status |
+|:---|:---|:---|:---|:---|:---:|
+| Beranek & Mellow 2e Eq. (13.117) | Piston resistance R1(x) = 1 - 2 J1(x)/x at x = 2ka = 2 | 0.423275 (+/-0.00001) | 0.423275 | 0 | &#9989; |
+| Beranek & Mellow 2e Eq. (13.118) | Piston reactance X1(x) = 2 H1(x)/x at x = 2ka = 2 | 0.646764 (+/-0.00001) | 0.646764 | 0 | &#9989; |
+| Beranek & Mellow 2e Eq. (13.117) (low-frequency limit) | R1 -> (ka)^2/2 as ka -> 0 (x = 0.02, ka = 0.01) | 0.00005 (+/-0.01%) | 0.00005 | 0 | &#9989; |
+| Beranek & Mellow 2e Eq. (4.151) | Radiation mass M = 8 rho a^3 / 3  (a = 0.1 m, rho = 1.206) | 0.003216 kg (+/-0 kg) | 0.003216 kg | 0 kg | &#9989; |
+| Beranek & Mellow 2e Eq. (13.102), Table 14.1 | First directivity null at ka sin(theta) = 3.8317 (first zero of J1) | 0 (+/-0.000001) | 0 | 0 | &#9989; |
+| Beranek & Mellow 2e §4.19 (half-space baffle) | Directivity index DI -> 10 lg 2 = 3.01 dB as ka -> 0 | 3.0103 dB (+/-0.001 dB) | 3.0103 dB | 0 dB | &#9989; |
+
+</details>
+
+<details>
+<summary>&#9989; <b>Industrial noise control</b>: 100% (9/9)</summary>
+
+| Standard | Quantity | Expected (norm) | Computed | &#916; | Status |
+|:---|:---|:---|:---|:---|:---:|
+| Bies 5e Eq. (8.111) | Expansion-chamber peak TL = 10 lg[1 + (1/4)(m - 1/m)^2], m = 4 at kL = pi/2 | 6.5472 dB (+/-0 dB) | 6.5472 dB | 0 dB | &#9989; |
+| Bies 5e Eq. (8.111) | Expansion-chamber trough TL = 0 at kL = pi (chamber transparent) | 0 dB (+/-0 dB) | 0 dB | 0 dB | &#9989; |
+| Bies 5e Eq. (8.44) / Example 8.1 | Quarter-wave tube tuning f = c/(4 l_e), l_e = 1.516 m -> 56.6 Hz | 56.6 Hz (+/-0.1 Hz) | 56.6 Hz | 0.003 Hz | &#9989; |
+| Bies 5e Eq. (8.46) | Helmholtz resonance f0 = (c/2pi) sqrt(S/(l_e V))  (S=1e-4, l_e=0.02, V=1e-3) | 122.067 Hz (+/-0 Hz) | 122.067 Hz | 0 Hz | &#9989; |
+| Bies 5e Eq. (8.73) | Side-branch TL = 20 lg abs(1 + rho c/(2 Sd Zb)) (QWT branch, closed form) | 0.1638 dB (+/-0 dB) | 0.1638 dB | 0 dB | &#9989; |
+| Bies 5e Eqs. (8.141)/(8.148) (four-pole insertion loss) | Insertion loss = transmission loss for the anechoic reference Zs=Zr=rho c/S | 6.2498 dB (= TL) | 6.2498 dB | 0 dB | &#9989; |
+| Bies 5e Eq. (8.275) (Wells' plenum method) | Plenum TL = -10 lg[S_out(cos0/pi r^2 + (1-a)/(Sw a))] (S_out=.1,r=1,Sw=20,a=.2) | 12.8541 dB (+/-0 dB) | 12.8541 dB | 0 dB | &#9989; |
+| Bies 5e Table 8.14 (ASHRAE end reflection, flush) | Duct end reflection D = 200 mm at 125 Hz = 10 dB (table node) | 10 dB (+/-0 dB) | 10 dB | 0 dB | &#9989; |
+| Bies 5e Eqs. (7.103), (7.111) (enclosure, fully absorbing limit) | Enclosure correction C -> 10 lg 0.3 = -5.23 dB as alpha_i -> 1 | -5.2288 dB (+/-0.001 dB) | -5.2288 dB | 0 dB | &#9989; |
+
+</details>
+

docs/README.md

@@ -23,6 +23,7 @@ Full documentation for phonometry. Also available as a website:
 - [Objective intelligibility (STOI & ESTOI)](objective-intelligibility.md): the correlation-based intelligibility measures for time-frequency weighted noisy speech from a clean/degraded pair: STOI (Taal et al. 2011), the clipped per-band envelope correlation, and ESTOI (Jensen & Taal 2016), the row- and column-normalised spectral correlation that tracks modulated maskers
 - [Electroacoustics: distortion & frequency response](electroacoustics.md): the IEC 60268-3 distortion set (THD, nth-order harmonic, THD+N and SINAD via AES17, SMPTE and CCIF intermodulation, dynamic intermodulation and weighted THD), the AES17 dynamic range and idle channel noise, and the Bendat & Piersol H1/H2 frequency-response estimators with the ordinary coherence γ²
 - [Swept-sine distortion and phase utilities](swept-sine-distortion.md): harmonic separation from one exponential sweep (Farina 2000) with the synchronized swept-sine of Novak et al. 2015 for coherent harmonic phases, THD as a function of the excitation frequency, and minimum phase from |H| (real cepstrum), group delay and excess phase
+- [Industrial noise control: silencers, HVAC & enclosures](noise-control.md): reactive silencers by the four-pole transmission-matrix method (expansion chambers, Helmholtz and quarter-wave resonators, extended tubes) with transmission and insertion loss, the HVAC duct methods (end reflection, elbows, plenums and flow-generated noise) and machine-enclosure insertion loss from a supplied panel R and the interior room constant (Bies, Hansen & Howard)
 - [Programme loudness & true peak](program-loudness.md): the ITU-R BS.1770-5 programme loudness (K-weighting, gated 400 ms blocks, channel weights including the Annex 3 positions) and the oversampled true-peak level in dBTP, with the EBU R 128 −23 LUFS practice, the Tech 3341 EBU Mode momentary/short-term/integrated meters and the Tech 3342 loudness range, validated against the official EBU test signals
 - [Underwater acoustics: radiated noise & pile driving](underwater-acoustics.md): the ISO 18405 reference levels (SPL, SEL, peak re 1 µPa), the ISO 17208 ship radiated noise level and equivalent monopole source level via the Lloyd's-mirror correction, and the ISO 18406 single-strike, peak and cumulative pile-driving sound exposure
 - [Underwater sound propagation](underwater-propagation.md): closed-form transmission loss (geometrical spreading plus volume absorption by Francois-Garrison, Ainslie-McColm or Thorp), the speed of sound in sea water (UNESCO/Chen-Millero, Del Grosso, Mackenzie) with the sound-speed profile, the passive/active sonar equation, seabed reflection loss (Rayleigh), the ocean ambient-noise spectrum (Wenz wind/thermal plus JOMOPANS-ECHO ship-traffic source levels) and numerical solvers (normal modes, ray tracing, parabolic equation)

docs/api-reference.md

@@ -356,6 +356,20 @@ cycle and warn on use; they are removed in 4.0.
 | `synchronized_sweep_signal` | `function` | **Synchronized exponential sweep (Novak et al. 2015, Eqs. 47/49).**<br>• `fs`, `f1` < `f2` [Hz], `seconds` (quantized by the rounding)<br>• `amplitude` (Default: 1.0), `fade` (Default: 0.0)<br>Rate L = round(f1·T̃/ln(f2/f1))/f1 → f1·L integer, so a L·ln(n) shift equals the nth harmonic and harmonic phases become system properties | `x = synchronized_sweep_signal(48000, 20.0, 6000.0, 4.0)`<br><br>• sweep samples (duration L·ln(f2/f1)) |
 | `swept_sine_distortion` | `function` | **Harmonic separation and THD(f) from one sweep (Farina 2000 / Novak et al. 2015).**<br>• `recorded`, `fs`, plus the generator parameters `f1`, `f2`, `seconds`<br>• `method`: 'synchronized' (Default; analytic deconvolution, meaningful phases) / 'farina' (classical ESS, magnitudes only)<br>• `n_harmonics` (Default: 5), `ir_length` (Default: largest power of two fitting the closest arrivals)<br>• `amplitude` (Default: 1.0), `fade` (Default: generator default), `remove_dc` (Default: True) | `res = swept_sine_distortion(y, fs, 20.0, 6000.0, 4.0)`<br><br>• `SweptSineDistortionResult` |
 | `SweptSineDistortionResult` | `dataclass` | **Swept-sine harmonic separation.**<br>• `frequencies` [Hz], `harmonic_responses`: complex H₁..H_N, `harmonic_irs` (windowed, centred), `delays` L·ln(n) [s]<br>• `thd_frequencies` [Hz], `thd`, `distortion_ratios`: \|Hₙ(n·f)\|/\|H₁(f)\| per order<br>• `fs`, `f1`, `f2`, `duration`, `rate` L, `method`, `n_harmonics`<br>• `.plot()`: \|Hₙ\| magnitudes + THD(f) | `res.thd, res.harmonic_responses` |
+| `radiating_piston` | `function` | **Radiation of a rigid baffled circular piston (Beranek & Mellow §4.19/§13.7).**<br>• `radius` a [m], `frequencies` [Hz]<br>• `speed_of_sound` (Default: 343), `density` (Default: 1.206)<br>• `angles`: directivity polar angles [rad] (Default: None) | `res = radiating_piston(0.1, freqs)`<br><br>• `RadiatingPistonResult` |
+| `piston_resistance` / `piston_reactance` | `function` | **Normalized piston radiation impedance R1 = 1−2J1(x)/x, X1 = 2H1(x)/x, x = 2ka (Beranek & Mellow Eqs. 13.117/13.118).** | `piston_resistance(2.0)  # 0.4233` |
+| `piston_directivity` | `function` | **Far-field piston directivity 2·J1(ka·sinθ)/(ka·sinθ).**<br>• `ka`, `theta` [rad] | `piston_directivity(5.0, 0.3)` |
+| `RadiatingPistonResult` | `dataclass` | **Piston radiation result.**<br>• `ka`, `resistance`/`reactance`, `radiation_resistance`/`radiation_reactance` [N·s/m], `radiation_mass` = 8ρa³/3 [kg], `directivity_index` [dB], `directivity`<br>• `.plot()`: R1/X1 vs ka | `res.radiation_mass` |
+| `expansion_chamber` | `function` | **Expansion-chamber silencer TL/IL (Bies Eq. 8.111, four-pole).**<br>• `frequencies`, `length` L, `chamber_area`, `pipe_area`<br>• `source_impedance`/`radiation_impedance` (Default: None → no IL) | `res = expansion_chamber(f, 0.3, 0.04, 0.01)`<br><br>• `ReactiveSilencerResult` |
+| `helmholtz_resonator` / `quarter_wave_resonator` | `function` | **Side-branch resonator silencers (Bies Eqs. 8.46/8.44).**<br>• Helmholtz: `duct_area`, `neck_area`, `neck_length`, `cavity_volume`<br>• Quarter-wave: `duct_area`, `length`, `branch_area` | `quarter_wave_resonator(f, 0.01, 1.5, 2e-3)` |
+| `extended_tube_chamber` | `function` | **Extended-inlet/outlet expansion chamber (Bies §8.9.7); reduces to `expansion_chamber` at zero extension.**<br>• `length`, `chamber_area`, `pipe_area`, `inlet_extension`, `outlet_extension` | `extended_tube_chamber(f, 0.4, 0.04, 0.01, inlet_extension=0.1)` |
+| `ReactiveSilencerResult` | `dataclass` | **Reactive-silencer result.**<br>• `frequencies`, `transmission_loss`, `insertion_loss` (or None), `transfer_matrix`, `kind`, `resonances`<br>• `.plot()`: TL/IL vs frequency | `res.transmission_loss` |
+| `end_reflection_loss` / `elbow_insertion_loss` | `function` | **HVAC duct end reflection & bend insertion loss (Bies Tables 8.14/8.11, ASHRAE).**<br>• end: `frequencies`, `diameter`, `termination` 'flush'/'free'<br>• elbow: `frequencies`, `width`, `bend_type`, `vanes`, `lined` | `end_reflection_loss(bands, 0.3)`<br><br>• `HvacSpectrumResult` |
+| `plenum_attenuation` | `function` | **Plenum-chamber TL by Wells' method (Bies Eq. 8.275).**<br>• `exit_area`, `line_of_sight`, `wall_area`, `mean_absorption`, `angle` | `plenum_attenuation(0.1, 1.0, 20.0, 0.2)  # dB` |
+| `flow_noise_straight_duct` / `flow_noise_bend` | `function` | **HVAC flow-generated (self) noise sound power (VDI 2081, Bies Eqs. 8.251/8.254).**<br>• `frequencies`, `flow_velocity`, `area` (+ `height` for the bend) | `flow_noise_straight_duct(bands, 10.0, 0.04)`<br><br>• `HvacSpectrumResult` |
+| `HvacSpectrumResult` | `dataclass` | **Per-frequency HVAC attenuation or regenerated Lw.**<br>• `frequencies`, `values`, `quantity`, `label`<br>• `.plot()` | `res.values` |
+| `enclosure_insertion_loss` | `function` | **Machine-enclosure insertion loss IL = R − C (Bies Eqs. 7.103/7.111); panel R supplied by the caller, never predicted.**<br>• `panel_transmission_loss` (per-band array or callable of f)<br>• `external_area`, `internal_area`, `internal_absorption`<br>• `frequencies` (required for a callable R) | `enclosure_insertion_loss(R, 6.0, 5.0, 0.3)`<br><br>• `EnclosureResult` |
+| `EnclosureResult` | `dataclass` | **Enclosure insertion-loss result.**<br>• `panel_transmission_loss`, `correction` C, `insertion_loss`, `room_constant`, `external_area`, `internal_area`<br>• `.plot()` | `res.insertion_loss` |
 | `program_loudness` | `function` | **EBU R 128 measurement set of a programme (BS.1770-5 + Tech 3341/3342).**<br>• `x`: signal (1D mono or 2D [channels, samples], 1.0 = 0 dBFS)<br>• `fs` [Hz]<br>• `weights`: per-channel Gi (Default: Table 3 by channel count)<br>• `momentary_step` [s] (Default: 0.01), `short_term_step` [s] (Default: 0.1)<br>• `oversample`: true-peak factor (Default: None → reaches 192 kHz) | `res = program_loudness(x, fs)`<br><br>• `ProgramLoudnessResult` |
 | `integrated_loudness` | `function` | **Programme (integrated) loudness with the two-stage gate (BS.1770-5 Annex 1).**<br>• `x`, `fs` [Hz]<br>• `weights` (Default: Table 3 by channel count)<br>Absolute gate −70 LKFS, relative −10 LU | `integrated_loudness(x, fs)  # LUFS` |
 | `loudness_range` | `function` | **Loudness range LRA (EBU Tech 3342).**<br>• `short_term_loudness`: S readings [LUFS] at ≥ 10 Hz<br>Cascaded gate (−70 LUFS abs, −20 LU rel), 10th-95th percentile spread | `loudness_range(res.short_term)  # LU` |

docs/electroacoustics.md

@@ -289,6 +289,44 @@ hide in that sentence:
   "1 m" figure is then a referred quantity, not what a microphone placed at
   1 m would read.
 
+## 6. Radiating piston: radiation impedance and directivity
+
+The rigid circular **piston in an infinite baffle** is the canonical radiator
+behind a loudspeaker cone, the open end of a duct and the radiation efficiency
+of any finite vibrating surface (Beranek & Mellow §4.19, §13.7). Its mechanical
+radiation impedance is `Z_r = rho c S (R1 + j X1)` with `S = pi a^2` and the
+dimensionless resistance and reactance functions (Eqs. (13.117), (13.118))
+
+```
+R1(x) = 1 - 2 J1(x) / x ,   X1(x) = 2 H1(x) / x ,   x = 2ka,
+```
+
+with `J1` the Bessel function and `H1` the Struve function of order one. At low
+frequency `R1 -> (ka)^2 / 2` and the reactance is mass-like with the radiation
+mass `M_r = 8 rho a^3 / 3` (Eq. (4.151)); at high frequency `R1 -> 1`, `X1 -> 0`
+and the piston radiates as into an infinite tube. The far field follows the
+directivity `D(theta) = 2 J1(ka sin theta) / (ka sin theta)`, whose first null
+is at `ka sin theta = 3.8317`.
+
+```python
+import numpy as np
+from phonometry import radiating_piston
+
+res = radiating_piston(radius=0.1, frequencies=np.geomspace(20, 20000, 200),
+                       angles=np.linspace(0.0, np.pi / 2, 91))
+print(round(res.radiation_mass, 4))               # 8 rho a^3 / 3, kg
+print(round(float(res.directivity_index[0]), 2))  # 3.01 dB half-space limit
+res.plot()                                         # R1 and X1 vs ka
+```
+
+`radiating_piston` returns a `RadiatingPistonResult` with the normalized
+`resistance`/`reactance`, the mechanical `radiation_resistance`/`radiation_reactance`,
+the `radiation_mass`, the `directivity_index`, the far-field `directivity`
+pattern (when `angles` are given) and `.plot()`. The building blocks
+`piston_resistance`, `piston_reactance` and `piston_directivity` are also
+callable directly. The piston is the companion radiator of the
+[industrial noise-control silencers](noise-control.md).
+
 ## References
 
 - Beranek, L. L., & Mellow, T. J. (2012). *Acoustics: Sound fields and