jmrp.io / PyOctaveBand
[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
Fork 0

Configure Feed

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

feat: complete the .plot() convention and tidy the plotting layer

Audit batch 10:

Every result dataclass with a plottable series now exposes .plot()
(eleven new): open-plan spatial decay with the Clause 6.2 regression
and rD/rP markers; outdoor attenuation with signed pos/neg stacking (a
net ground gain stacks below zero); impedance-tube alpha(f) with the
muted |r| companion; Monte Carlo histogram with the coverage interval -
via a new opt-in monte_carlo(keep_samples=True) that stores the sample
on the result (the 8 MB/1M-trial cost stays off the default path, and
.plot() without samples raises with the hint); occupational-exposure
per-task contributions with the LEX and LEX+U lines; plus simple
renderers for static airflow, airborne/impact prediction,
airborne/impact insulation bands and band uncertainty. TransferMatrix
and ScatteringUncertainty are documented skips (value objects with no
stored frequency axis). The five priority plots are mentioned in their
guides on all three surfaces.

Plotting hygiene: _freq_axis in the enclosed-space renderer (minor-tick
suppression restored); shared _band_axis/_fractile_band/_hatch_invalid
helpers and plot_weighted_absorption folded into _plot_rating; named
color constants with the three neutral greys collapsed into one _C_MUTED
(the per-metric tonality red and roughness brown stay, commented);
every title now leads with its standard designation (17 aligned, with
per-type designations where one renderer serves several standards);
band x-labels unified to the dominant wording; the 11 result: Any
renderers typed; kwargs documented; the missing small legend fixed;
user color propagates into the companion fill_between in the five
loudness/tonality/roughness renderers.

Tests: 11 new kwargs-forwarding cases, 13 content/raise tests, and the
external-ax contract extended to all single-axes results (64/64).

José M. Requena Plens (Jul 10, 2026, 12:58 PM +0200) 49b1a2b0 0eed2cd6

+1459 -239
+7 -3
docs/gum-uncertainty.md
··· 113 113 ] 114 114 model = lambda a, b, c, d: a + b + c + d 115 115 result = ph.combine_uncertainty(model, quantities) 116 - mc = ph.monte_carlo(model, quantities, trials=1_000_000, coverage=0.95, seed=1) 116 + mc = ph.monte_carlo(model, quantities, trials=1_000_000, coverage=0.95, seed=1, 117 + keep_samples=True) 117 118 k, U = result.expanded(0.95) 118 119 119 - # One line for the budget: 120 + # One line per panel — the budget and the Monte Carlo output distribution: 120 121 result.plot() 122 + mc.plot() 121 123 plt.show() 122 124 123 125 # By hand, both panels — budget bars and the Monte Carlo output distribution: ··· 145 147 `sensitivities`, the per-input `contributions` and the `effective_dof`; its 146 148 `.plot()` draws the budget and `.expanded(coverage)` returns the pair 147 149 $(k, U)$. The `MonteCarloResult` carries the `value`, `standard_uncertainty`, 148 - the coverage `interval` and its `coverage`. The building-acoustics uncertainty 150 + the coverage `interval` and its `coverage`; with `keep_samples=True` it also 151 + retains the output `samples`, and its `.plot()` draws the output histogram 152 + with the coverage interval marked (the right panel above). The building-acoustics uncertainty 149 153 of ISO 12999-1 — which combines reproducibility terms for a single-number 150 154 rating — is a separate, domain-specific budget. 151 155
+2 -1
docs/materials.md
··· 307 307 `ImpedanceTubeResult` with absorption, reflection factor, surface impedance and 308 308 normalised impedance, applying the plane-wave frequency-range check and optional 309 309 tube attenuation; correct any microphone mismatch beforehand with 310 - `apply_mic_calibration`. 310 + `apply_mic_calibration`. Its `.plot()` draws the absorption spectrum 311 + $\alpha(f)$ with the reflection-factor magnitude $|r|$ overlaid. 311 312 312 313 **Transmission loss (ASTM E2611).** With four microphones — two upstream, two 313 314 downstream of the sample — a two-load (or one-load) measurement recovers the
+19 -1
docs/occupational-exposure.md
··· 29 29 30 30 <picture><source media="(prefers-color-scheme: dark)" srcset="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/exposure_uncertainty_dark.png"><img src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/exposure_uncertainty.png" alt="ISO 9612 Annex D task-based exposure: the three task LEX,8h contributions as bars, the energy-summed daily LEX,8h line and the one-sided 95 % upper limit LEX,8h + U band above it" width="80%"></picture> 31 31 32 + <details> 33 + <summary>Show the code for this figure</summary> 34 + 35 + ```python 36 + import matplotlib.pyplot as plt 37 + 38 + # One line, with the Annex D `res` computed below: task contribution bars 39 + # plus the LEX,8h and LEX,8h + U lines. 40 + res.plot() 41 + plt.show() 42 + ``` 43 + 44 + </details> 45 + 46 + 32 47 ```python 33 48 from phonometry.occupational_exposure import ( 34 49 Task, task_based_exposure, job_based_exposure, full_day_exposure, ··· 106 121 107 122 All three return an `ExposureResult` with `lex_8h`, `combined_standard_uncertainty` 108 123 $u$, `expanded_uncertainty` $U = 1.65\ u$, `upper_limit` = $L_{EX,8h} + U$, 109 - `sampling_advisory`, and (task-based) the per-task `tasks` breakdown. 124 + `sampling_advisory`, and (task-based) the per-task `tasks` breakdown; the 125 + result's `.plot()` draws the per-task contribution bars with the $L_{EX,8h}$ 126 + and upper-limit lines (task-based results only, since the other strategies 127 + carry no per-task breakdown). 110 128 111 129 ## See also 112 130
+6 -1
docs/outdoor-propagation.md
··· 184 184 relative_humidity=70.0, 185 185 ) 186 186 187 + # One line — the same stacked breakdown with the total overlaid: 188 + att.plot() 189 + 190 + # By hand: 187 191 x = np.arange(len(bands)) 188 192 fig, ax = plt.subplots() 189 193 # Separate positive and negative baselines: a negative term (Agr is a net ··· 294 298 | `projected_distance` | float or None | m | default $\sqrt{d^2-(h_s-h_r)^2}$ | Ground-plane $d_p$ | 295 299 296 300 Returns an `OutdoorAttenuation` with `a_div`, `a_atm`, `a_gr`, `a_bar`, 297 - `a_total` and `d_omega`, all one value per band. 301 + `a_total` and `d_omega`, all one value per band; its `.plot()` draws the 302 + stacked per-band breakdown with the total overlaid (the figure above). 298 303 299 304 ### `Barrier` fields 300 305
+10 -1
docs/room-acoustics.md
··· 352 352 <summary>Show the code for this figure</summary> 353 353 354 354 ```python 355 + # One line: the D2,S regression rebuilt from the result fields, with the 356 + # rD / rP crossings marked (the figure above adds the measured points and 357 + # the STI axis on top of it): 358 + m.plot() 359 + ``` 360 + 361 + ```python 355 362 import matplotlib.pyplot as plt 356 363 357 364 # Spatial decay: measured Lp,A,S vs distance on a log axis, the D2,S ··· 394 401 | `spl_a_speech` | 1D array | dB | same length | A-weighted speech level `Lp,A,S,n` per position | 395 402 | `sti_values` | 1D array | — | same length | STI per position (full IEC 60268-16 method) | 396 403 397 - Returns an `OpenPlanResult` with `d2s`, `lp_as_4m`, `rd` and `rp`. 404 + Returns an `OpenPlanResult` with `d2s`, `lp_as_4m`, `rd` and `rp`; its 405 + `.plot()` redraws the Clause 6.2 spatial-decay regression from those four 406 + fields and marks `rd` / `rp`. 398 407 `d2s`/`lp_as_4m` are `nan` if fewer than two positions fall in 2–16 m; 399 408 `rd`/`rp` are `nan` when STI does not decrease with distance. The per-position 400 409 STI can itself be measured with the STIPA tools in the
+7 -3
site/src/content/docs/es/guides/gum-uncertainty.md
··· 116 116 ] 117 117 model = lambda a, b, c, d: a + b + c + d 118 118 result = ph.combine_uncertainty(model, quantities) 119 - mc = ph.monte_carlo(model, quantities, trials=1_000_000, coverage=0.95, seed=1) 119 + mc = ph.monte_carlo(model, quantities, trials=1_000_000, coverage=0.95, seed=1, 120 + keep_samples=True) 120 121 k, U = result.expanded(0.95) 121 122 122 - # En una línea, el balance: 123 + # En una línea por panel — el balance y la distribución de salida de Monte Carlo: 123 124 result.plot() 125 + mc.plot() 124 126 plt.show() 125 127 126 128 # A mano, ambos paneles — las barras del balance y la distribución de salida de Monte Carlo: ··· 148 150 `sensitivities`, las `contributions` por entrada y los `effective_dof`; su 149 151 `.plot()` dibuja el balance y `.expanded(coverage)` devuelve el par $(k, U)$. El 150 152 `MonteCarloResult` lleva el `value`, la `standard_uncertainty`, el `interval` de 151 - cobertura y su `coverage`. La incertidumbre de acústica de la edificación de 153 + cobertura y su `coverage`; con `keep_samples=True` conserva además las 154 + muestras de salida (`samples`), y su `.plot()` dibuja el histograma de salida 155 + con el intervalo de cobertura marcado (el panel derecho de arriba). La incertidumbre de acústica de la edificación de 152 156 ISO 12999-1 — que combina términos de reproducibilidad para una magnitud de 153 157 número único — es un balance aparte, específico de ese dominio. 154 158
+3 -1
site/src/content/docs/es/guides/materials.md
··· 316 316 `ImpedanceTubeResult` con la absorción, el factor de reflexión, la impedancia 317 317 superficial y la impedancia normalizada, aplicando la comprobación del rango de 318 318 frecuencias de onda plana y la atenuación opcional del tubo; corrige antes 319 - cualquier desajuste de micrófonos con `apply_mic_calibration`. 319 + cualquier desajuste de micrófonos con `apply_mic_calibration`. Su `.plot()` 320 + dibuja el espectro de absorción $\alpha(f)$ con el módulo del factor de 321 + reflexión $|r|$ superpuesto. 320 322 321 323 **Pérdida por transmisión (ASTM E2611).** Con cuatro micrófonos —dos aguas 322 324 arriba, dos aguas abajo de la muestra— una medición de dos cargas (o de una
+19 -1
site/src/content/docs/es/guides/occupational-exposure.md
··· 35 35 36 36 <img class="light-only" src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/exposure_uncertainty_es.png" alt="Exposición por tareas del anexo D de la ISO 9612: las tres contribuciones de tarea a LEX,8h como barras, la línea del LEX,8h diario sumado en energía y la banda del límite superior unilateral al 95 % LEX,8h + U por encima" style="width:80%"><img class="dark-only" src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/exposure_uncertainty_es_dark.png" alt="Exposición por tareas del anexo D de la ISO 9612: las tres contribuciones de tarea a LEX,8h como barras, la línea del LEX,8h diario sumado en energía y la banda del límite superior unilateral al 95 % LEX,8h + U por encima" style="width:80%"> 37 37 38 + <details> 39 + <summary>Mostrar el código de esta figura</summary> 40 + 41 + ```python 42 + import matplotlib.pyplot as plt 43 + 44 + # En una línea, con el `res` del anexo D calculado más abajo: barras de 45 + # contribución por tarea más las líneas de LEX,8h y LEX,8h + U. 46 + res.plot() 47 + plt.show() 48 + ``` 49 + 50 + </details> 51 + 52 + 38 53 ```python 39 54 from phonometry.occupational_exposure import ( 40 55 Task, task_based_exposure, job_based_exposure, full_day_exposure, ··· 116 131 117 132 Las tres devuelven un `ExposureResult` con `lex_8h`, `combined_standard_uncertainty` 118 133 $u$, `expanded_uncertainty` $U = 1{,}65\ u$, `upper_limit` = $L_{EX,8h} + U$, 119 - `sampling_advisory` y (basada en tareas) el desglose por tarea en `tasks`. 134 + `sampling_advisory` y (basada en tareas) el desglose por tarea en `tasks`; el 135 + `.plot()` del resultado dibuja las barras de contribución por tarea con las 136 + líneas de $L_{EX,8h}$ y del límite superior (solo resultados por tareas: las 137 + otras estrategias no llevan desglose por tarea). 120 138 121 139 ## Véase también 122 140
+6 -1
site/src/content/docs/es/guides/outdoor-propagation.md
··· 193 193 relative_humidity=70.0, 194 194 ) 195 195 196 + # En una línea — el mismo desglose apilado con el total superpuesto: 197 + att.plot() 198 + 199 + # A mano: 196 200 x = np.arange(len(bands)) 197 201 fig, ax = plt.subplots() 198 202 # Líneas base positiva y negativa separadas: un término negativo (Agr es una ··· 307 311 | `projected_distance` | float o None | m | def. $\sqrt{d^2-(h_s-h_r)^2}$ | $d_p$ proyectada sobre el suelo | 308 312 309 313 Devuelve un `OutdoorAttenuation` con `a_div`, `a_atm`, `a_gr`, `a_bar`, 310 - `a_total` y `d_omega`, todos un valor por banda. 314 + `a_total` y `d_omega`, todos un valor por banda; su `.plot()` dibuja el 315 + desglose apilado por bandas con el total superpuesto (la figura anterior). 311 316 312 317 ### Campos de `Barrier` 313 318
+10 -1
site/src/content/docs/es/guides/room-acoustics.md
··· 366 366 <summary>Ver el código de esta figura</summary> 367 367 368 368 ```python 369 + # En una línea: la regresión D2,S reconstruida con los campos del resultado, 370 + # con los cruces rD / rP marcados (la figura superior añade además los puntos 371 + # medidos y el eje de STI): 372 + m.plot() 373 + ``` 374 + 375 + ```python 369 376 import matplotlib.pyplot as plt 370 377 371 378 # Decaimiento espacial: Lp,A,S medido frente a la distancia en eje logarítmico, ··· 408 415 | `spl_a_speech` | array 1D | dB | misma longitud | Nivel de habla ponderado A `Lp,A,S,n` por posición | 409 416 | `sti_values` | array 1D | — | misma longitud | STI por posición (método completo de IEC 60268-16) | 410 417 411 - Devuelve un `OpenPlanResult` con `d2s`, `lp_as_4m`, `rd` y `rp`. 418 + Devuelve un `OpenPlanResult` con `d2s`, `lp_as_4m`, `rd` y `rp`; su 419 + `.plot()` redibuja la regresión de decaimiento espacial del apartado 6.2 a 420 + partir de esos cuatro campos y marca `rd` / `rp`. 412 421 `d2s`/`lp_as_4m` son `nan` si menos de dos posiciones caen en 2–16 m; 413 422 `rd`/`rp` son `nan` cuando el STI no decrece con la distancia. El STI por 414 423 posición puede medirse a su vez con las herramientas STIPA de la
+7 -3
site/src/content/docs/guides/gum-uncertainty.md
··· 114 114 ] 115 115 model = lambda a, b, c, d: a + b + c + d 116 116 result = ph.combine_uncertainty(model, quantities) 117 - mc = ph.monte_carlo(model, quantities, trials=1_000_000, coverage=0.95, seed=1) 117 + mc = ph.monte_carlo(model, quantities, trials=1_000_000, coverage=0.95, seed=1, 118 + keep_samples=True) 118 119 k, U = result.expanded(0.95) 119 120 120 - # One line for the budget: 121 + # One line per panel — the budget and the Monte Carlo output distribution: 121 122 result.plot() 123 + mc.plot() 122 124 plt.show() 123 125 124 126 # By hand, both panels — budget bars and the Monte Carlo output distribution: ··· 146 148 `sensitivities`, the per-input `contributions` and the `effective_dof`; its 147 149 `.plot()` draws the budget and `.expanded(coverage)` returns the pair 148 150 $(k, U)$. The `MonteCarloResult` carries the `value`, `standard_uncertainty`, 149 - the coverage `interval` and its `coverage`. The building-acoustics uncertainty 151 + the coverage `interval` and its `coverage`; with `keep_samples=True` it also 152 + retains the output `samples`, and its `.plot()` draws the output histogram 153 + with the coverage interval marked (the right panel above). The building-acoustics uncertainty 150 154 of ISO 12999-1 — which combines reproducibility terms for a single-number 151 155 rating — is a separate, domain-specific budget. 152 156
+2 -1
site/src/content/docs/guides/materials.md
··· 308 308 `ImpedanceTubeResult` with absorption, reflection factor, surface impedance and 309 309 normalised impedance, applying the plane-wave frequency-range check and optional 310 310 tube attenuation; correct any microphone mismatch beforehand with 311 - `apply_mic_calibration`. 311 + `apply_mic_calibration`. Its `.plot()` draws the absorption spectrum 312 + $\alpha(f)$ with the reflection-factor magnitude $|r|$ overlaid. 312 313 313 314 **Transmission loss (ASTM E2611).** With four microphones — two upstream, two 314 315 downstream of the sample — a two-load (or one-load) measurement recovers the
+19 -1
site/src/content/docs/guides/occupational-exposure.md
··· 30 30 31 31 <img class="light-only" src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/exposure_uncertainty.png" alt="ISO 9612 Annex D task-based exposure: the three task LEX,8h contributions as bars, the energy-summed daily LEX,8h line and the one-sided 95 % upper limit LEX,8h + U band above it" style="width:80%"><img class="dark-only" src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/exposure_uncertainty_dark.png" alt="ISO 9612 Annex D task-based exposure: the three task LEX,8h contributions as bars, the energy-summed daily LEX,8h line and the one-sided 95 % upper limit LEX,8h + U band above it" style="width:80%"> 32 32 33 + <details> 34 + <summary>Show the code for this figure</summary> 35 + 36 + ```python 37 + import matplotlib.pyplot as plt 38 + 39 + # One line, with the Annex D `res` computed below: task contribution bars 40 + # plus the LEX,8h and LEX,8h + U lines. 41 + res.plot() 42 + plt.show() 43 + ``` 44 + 45 + </details> 46 + 47 + 33 48 ```python 34 49 from phonometry.occupational_exposure import ( 35 50 Task, task_based_exposure, job_based_exposure, full_day_exposure, ··· 107 122 108 123 All three return an `ExposureResult` with `lex_8h`, `combined_standard_uncertainty` 109 124 $u$, `expanded_uncertainty` $U = 1.65\ u$, `upper_limit` = $L_{EX,8h} + U$, 110 - `sampling_advisory`, and (task-based) the per-task `tasks` breakdown. 125 + `sampling_advisory`, and (task-based) the per-task `tasks` breakdown; the 126 + result's `.plot()` draws the per-task contribution bars with the $L_{EX,8h}$ 127 + and upper-limit lines (task-based results only, since the other strategies 128 + carry no per-task breakdown). 111 129 112 130 ## See also 113 131
+6 -1
site/src/content/docs/guides/outdoor-propagation.md
··· 185 185 relative_humidity=70.0, 186 186 ) 187 187 188 + # One line — the same stacked breakdown with the total overlaid: 189 + att.plot() 190 + 191 + # By hand: 188 192 x = np.arange(len(bands)) 189 193 fig, ax = plt.subplots() 190 194 # Separate positive and negative baselines: a negative term (Agr is a net ··· 295 299 | `projected_distance` | float or None | m | default $\sqrt{d^2-(h_s-h_r)^2}$ | Ground-plane $d_p$ | 296 300 297 301 Returns an `OutdoorAttenuation` with `a_div`, `a_atm`, `a_gr`, `a_bar`, 298 - `a_total` and `d_omega`, all one value per band. 302 + `a_total` and `d_omega`, all one value per band; its `.plot()` draws the 303 + stacked per-band breakdown with the total overlaid (the figure above). 299 304 300 305 ### `Barrier` fields 301 306
+10 -1
site/src/content/docs/guides/room-acoustics.md
··· 353 353 <summary>Show the code for this figure</summary> 354 354 355 355 ```python 356 + # One line: the D2,S regression rebuilt from the result fields, with the 357 + # rD / rP crossings marked (the figure above adds the measured points and 358 + # the STI axis on top of it): 359 + m.plot() 360 + ``` 361 + 362 + ```python 356 363 import matplotlib.pyplot as plt 357 364 358 365 # Spatial decay: measured Lp,A,S vs distance on a log axis, the D2,S ··· 395 402 | `spl_a_speech` | 1D array | dB | same length | A-weighted speech level `Lp,A,S,n` per position | 396 403 | `sti_values` | 1D array | — | same length | STI per position (full IEC 60268-16 method) | 397 404 398 - Returns an `OpenPlanResult` with `d2s`, `lp_as_4m`, `rd` and `rp`. 405 + Returns an `OpenPlanResult` with `d2s`, `lp_as_4m`, `rd` and `rp`; its 406 + `.plot()` redraws the Clause 6.2 spatial-decay regression from those four 407 + fields and marks `rd` / `rp`. 399 408 `d2s`/`lp_as_4m` are `nan` if fewer than two positions fall in 2–16 m; 400 409 `rd`/`rp` are `nan` when STI does not decrease with distance. The per-position 401 410 STI can itself be measured with the STIPA tools in the
+848 -212
src/phonometry/_plotting.py
··· 23 23 24 24 from __future__ import annotations 25 25 26 - from typing import TYPE_CHECKING, Any, Final, cast 26 + from typing import TYPE_CHECKING, Any, Final, Sequence, cast 27 27 28 28 import numpy as np 29 29 30 30 if TYPE_CHECKING: 31 31 from matplotlib.axes import Axes 32 + from matplotlib.container import BarContainer 32 33 34 + from .absorption_rating import AbsorptionRatingResult 35 + from .airflow_resistance import StaticAirflowResult 36 + from .building_prediction import AirbornePredictionResult, ImpactPredictionResult 37 + from .building_uncertainty import BandUncertainty 38 + from .impedance_tube import ImpedanceTubeResult 39 + from .insulation import ( 40 + AirborneInsulationResult, 41 + FacadeInsulationResult, 42 + ImpactInsulationResult, 43 + ImpactRatingResult, 44 + WeightedRatingResult, 45 + ) 33 46 from .intensity import IntensityResult 34 47 from .loudness_zwicker import ZwickerLoudness 35 48 from .loudness_ecma import EcmaLoudness 36 49 from .loudness_moore_glasberg import MooreGlasbergLoudness 37 50 from .loudness_moore_glasberg_time import MooreGlasbergTimeVaryingLoudness 51 + from .occupational_exposure import ExposureResult 52 + from .open_plan import OpenPlanResult 53 + from .outdoor_propagation import OutdoorAttenuation 54 + from .road_absorption import InsituAbsorptionResult 38 55 from .room_acoustics import DecayCurve, RoomAcousticsResult 39 56 from .room_ir import ImpulseResponseResult 40 57 from .hearing import AgeThresholdResult 41 58 from .enclosed_space_absorption import ReverberationResult 59 + from .human_vibration import ( 60 + DailyVibrationExposure, 61 + WeightedSpectrum, 62 + WeightingResponse, 63 + ) 42 64 from .noise_induced_hearing_loss import HtlanResult, NiptsResult 43 65 from .multiple_shock_vibration import MultipleShockResult 44 66 from .impulse_prominence import ImpulseProminenceResult 45 67 from .room_noise import NCResult, RCResult 68 + from .scattering_diffusion import DiffusionResult, ScatteringResult 69 + from .sound_power import SoundPowerResult 70 + from .sound_power_intensity import SoundPowerIntensityResult 71 + from .sound_power_reverberation import ReverberationSoundPowerResult 46 72 from .tonality_ecma import EcmaTonality 47 73 from .roughness_ecma import EcmaRoughness 48 74 from .sii import SIIResult 49 75 from .sti import STIResult 50 - from .uncertainty import UncertaintyResult 76 + from .uncertainty import MonteCarloResult, UncertaintyResult 51 77 52 78 _INSTALL_HINT = ( 53 79 "Plotting requires matplotlib. Install it with: pip install phonometry[plot]" ··· 59 85 #: Default legend placement shared by the band-spectrum figures. 60 86 _LEGEND_UPPER_RIGHT: Final = "upper right" 61 87 88 + # --------------------------------------------------------------------------- 89 + # Shared artist colors (matplotlib "tab10" hues plus neutral greys). 90 + # 91 + # The measured/primary series is drawn in _C_PRIMARY, reference curves and 92 + # limit lines in _C_REFERENCE, unfavourable/secondary annotations in 93 + # _C_SECONDARY, and de-emphasised context (invalid bands, threshold guides, 94 + # curve families, companion series) in the single neutral _C_MUTED. 95 + # 96 + # Two renderers keep a fixed per-metric identity color instead: ECMA-418-2 97 + # tonality is red and roughness is brown across the documentation, so those 98 + # literals live in their renderers with a comment, not here. 99 + # --------------------------------------------------------------------------- 100 + _C_PRIMARY: Final = "#1f77b4" 101 + _C_PRIMARY_LIGHT: Final = "#aec7e8" 102 + _C_REFERENCE: Final = "#d62728" 103 + _C_SECONDARY: Final = "#ff7f0e" 104 + _C_SECONDARY_LIGHT: Final = "#ffbb78" 105 + _C_TERTIARY: Final = "#2ca02c" 106 + _C_QUATERNARY: Final = "#9467bd" 107 + _C_MUTED: Final = "#9e9e9e" 108 + _C_EDGE: Final = "#555555" 109 + 110 + #: Standard-normal quantile of 0.9 (the 10 % / 90 % fractile offset). 111 + _Z90: Final = 1.2816 112 + 62 113 63 114 def _import_pyplot() -> Any: 64 115 """Import :mod:`matplotlib.pyplot` lazily with an actionable error.""" ··· 106 157 return text 107 158 108 159 160 + def _band_axis( 161 + ax: Axes, 162 + labels_or_freqs: "np.ndarray | Sequence[str] | Sequence[float]", 163 + *, 164 + xlabel: str | None = "Frequency [Hz]", 165 + ) -> np.ndarray: 166 + """Categorical band x-axis: evenly spaced positions labelled with centres. 167 + 168 + Band bars/curves are drawn on evenly spaced positions (a *linear* axis) 169 + so they stay legible; the tick labels carry the band centres (numeric 170 + input is shortened via :func:`_format_freq`) or the given strings. 171 + Returns the positions. ``xlabel=None`` leaves the axis label untouched 172 + (e.g. a shared-x upper panel). 173 + """ 174 + labels = [ 175 + item if isinstance(item, str) else _format_freq(float(item)) 176 + for item in list(labels_or_freqs) 177 + ] 178 + positions = np.arange(len(labels), dtype=np.float64) 179 + ax.set_xticks(positions) 180 + ax.set_xticklabels(labels, rotation=45, ha="right") 181 + if xlabel is not None: 182 + ax.set_xlabel(xlabel) 183 + return positions 184 + 185 + 186 + def _fractile_band( 187 + ax: Axes, 188 + freqs: np.ndarray, 189 + median: np.ndarray, 190 + spread_lower: np.ndarray, 191 + spread_upper: np.ndarray, 192 + *, 193 + color: str, 194 + floor: float | None = None, 195 + ) -> None: 196 + """Shade the 10-90 % fractile band around a median spectrum. 197 + 198 + The band spans ``median - z90*spread_lower`` to ``median + 199 + z90*spread_upper`` with the standard-normal :data:`_Z90` quantile; 200 + ``floor`` clamps the lower edge (e.g. NIPTS cannot be negative). 201 + """ 202 + lower = median - _Z90 * spread_lower 203 + if floor is not None: 204 + lower = np.maximum(lower, floor) 205 + ax.fill_between(freqs, lower, median + _Z90 * spread_upper, 206 + color=color, alpha=0.5, label="10-90 % fractile band") 207 + 208 + 209 + def _hatch_invalid(bars: "BarContainer", mask: np.ndarray) -> None: 210 + """Hatch (and outline) the bars flagged invalid/unusable by ``mask``.""" 211 + for bar, bad in zip(bars, np.asarray(mask, dtype=bool), strict=True): 212 + if bad: 213 + bar.set_hatch("//") 214 + bar.set_edgecolor(_C_EDGE) 215 + 216 + 109 217 # --------------------------------------------------------------------------- 110 218 # Zwicker loudness (ISO 532-1) 111 219 # --------------------------------------------------------------------------- ··· 139 247 else: 140 248 ax_specific = ax if ax is not None else _new_axes() 141 249 142 - kwargs.setdefault("color", "#1f77b4") 250 + kwargs.setdefault("color", _C_PRIMARY) 143 251 ax_specific.plot(bark, specific, **kwargs) 144 - ax_specific.fill_between(bark, specific, color="#1f77b4", alpha=0.25) 252 + ax_specific.fill_between(bark, specific, color=kwargs["color"], alpha=0.25) 145 253 ax_specific.set_xlabel("Critical-band rate z [Bark]") 146 254 ax_specific.set_ylabel("Specific loudness N' [sone/Bark]") 147 255 ax_specific.set_xlim(0.0, bark[-1]) 148 256 ax_specific.set_ylim(bottom=0.0) 149 257 ax_specific.set_title( 150 - f"Loudness N = {result.loudness:.2f} sone ({result.loudness_level:.1f} phon)" 258 + f"ISO 532-1 loudness N = {result.loudness:.2f} sone " 259 + f"({result.loudness_level:.1f} phon)" 151 260 ) 152 261 ax_specific.grid(True, alpha=0.3) 153 262 ··· 157 266 ax_time = cast("Axes", axes[1]) 158 267 time = np.asarray(result.time, dtype=np.float64) 159 268 lvt = np.asarray(result.loudness_vs_time, dtype=np.float64) 160 - ax_time.plot(time, lvt, color="#2ca02c", label="N(t)") 269 + ax_time.plot(time, lvt, color=_C_TERTIARY, label="N(t)") 161 270 if result.n5 is not None: 162 271 ax_time.axhline( 163 - result.n5, color="#d62728", ls="--", lw=1, label=f"N5={result.n5:.2f}" 272 + result.n5, color=_C_REFERENCE, ls="--", lw=1, label=f"N5={result.n5:.2f}" 164 273 ) 165 274 if result.n10 is not None: 166 275 ax_time.axhline( 167 - result.n10, color="#ff7f0e", ls=":", lw=1, label=f"N10={result.n10:.2f}" 276 + result.n10, color=_C_SECONDARY, ls=":", lw=1, label=f"N10={result.n10:.2f}" 168 277 ) 169 278 ax_time.set_xlabel("Time [s]") 170 279 ax_time.set_ylabel("Loudness N [sone]") ··· 204 313 else: 205 314 ax_specific = cast("Axes", ax) 206 315 207 - kwargs.setdefault("color", "#1f77b4") 316 + kwargs.setdefault("color", _C_PRIMARY) 208 317 ax_specific.plot(bark, specific, **kwargs) 209 - ax_specific.fill_between(bark, specific, color="#1f77b4", alpha=0.25) 318 + ax_specific.fill_between(bark, specific, color=kwargs["color"], alpha=0.25) 210 319 ax_specific.set_xlabel("Critical-band rate z [Bark_HMS]") 211 320 ax_specific.set_ylabel("Specific loudness N' [sone_HMS/Bark_HMS]") 212 321 ax_specific.set_xlim(0.0, bark[-1]) 213 322 ax_specific.set_ylim(bottom=0.0) 214 - ax_specific.set_title(f"Loudness N = {result.loudness:.2f} sone_HMS") 323 + ax_specific.set_title(f"ECMA-418-2 loudness N = {result.loudness:.2f} sone_HMS") 215 324 ax_specific.grid(True, alpha=0.3) 216 325 217 326 if not two_panel: ··· 220 329 ax_time = cast("Axes", axes[1]) 221 330 time = np.asarray(result.time, dtype=np.float64) 222 331 lvt = np.asarray(result.loudness_vs_time, dtype=np.float64) 223 - ax_time.plot(time, lvt, color="#2ca02c", label="N(l)") 332 + ax_time.plot(time, lvt, color=_C_TERTIARY, label="N(l)") 224 333 ax_time.set_xlabel("Time [s]") 225 334 ax_time.set_ylabel("Loudness N [sone_HMS]") 226 335 ax_time.set_ylim(bottom=0.0) ··· 249 358 erb_number = np.asarray(result.erb_number, dtype=np.float64) 250 359 ax = ax if ax is not None else _new_axes() 251 360 252 - kwargs.setdefault("color", "#1f77b4") 361 + kwargs.setdefault("color", _C_PRIMARY) 253 362 ax.plot(erb_number, specific, **kwargs) 254 - ax.fill_between(erb_number, specific, color="#1f77b4", alpha=0.25) 363 + ax.fill_between(erb_number, specific, color=kwargs["color"], alpha=0.25) 255 364 ax.set_xlabel("ERB number [Cam]") 256 365 ax.set_ylabel("Specific loudness N' [sone/Cam]") 257 366 ax.set_xlim(erb_number[0], erb_number[-1]) 258 367 ax.set_ylim(bottom=0.0) 259 368 ax.set_title( 260 - f"Loudness N = {result.loudness:.2f} sone ({result.loudness_level:.1f} phon)" 369 + f"ISO 532-2 loudness N = {result.loudness:.2f} sone " 370 + f"({result.loudness_level:.1f} phon)" 261 371 ) 262 372 ax.grid(True, alpha=0.3) 263 373 return ax ··· 279 389 ltl = np.asarray(result.long_term_loudness, dtype=np.float64) 280 390 ax = ax if ax is not None else _new_axes() 281 391 282 - ax.plot(time, stl, color="#aec7e8", lw=1.0, label="Short-term loudness") 283 - kwargs.setdefault("color", "#1f77b4") 392 + ax.plot(time, stl, color=_C_PRIMARY_LIGHT, lw=1.0, label="Short-term loudness") 393 + kwargs.setdefault("color", _C_PRIMARY) 284 394 kwargs.setdefault("lw", 1.8) 285 395 ax.plot(time, ltl, label="Long-term loudness", **kwargs) 286 - ax.axhline(result.n_max, color="#d62728", ls="--", lw=1.0, alpha=0.7) 396 + ax.axhline(result.n_max, color=_C_REFERENCE, ls="--", lw=1.0, alpha=0.7) 287 397 ax.set_xlabel("Time [s]") 288 398 ax.set_ylabel("Loudness [sone]") 289 399 if time.size: 290 400 ax.set_xlim(time[0], time[-1]) 291 401 ax.set_ylim(bottom=0.0) 292 402 ax.set_title( 293 - f"Peak long-term loudness N = {result.n_max:.2f} sone " 403 + f"ISO 532-3 peak long-term loudness N = {result.n_max:.2f} sone " 294 404 f"({result.loudness_level_max:.1f} phon)" 295 405 ) 296 - ax.legend(loc="best") 406 + ax.legend(loc="best", fontsize="small") 297 407 ax.grid(True, alpha=0.3) 298 408 return ax 299 409 ··· 323 433 else: 324 434 ax_specific = cast("Axes", ax) 325 435 436 + # Tonality's per-metric identity color is red across the documentation 437 + # figures (roughness is brown); kept literal on purpose, see the module 438 + # color-constant note. 326 439 kwargs.setdefault("color", "#d62728") 327 440 ax_specific.plot(bark, specific, **kwargs) 328 - ax_specific.fill_between(bark, specific, color="#d62728", alpha=0.25) 441 + ax_specific.fill_between(bark, specific, color=kwargs["color"], alpha=0.25) 329 442 ax_specific.set_xlabel("Critical-band rate z [Bark_HMS]") 330 443 ax_specific.set_ylabel("Specific tonality T' [tu_HMS]") 331 444 ax_specific.set_xlim(0.0, bark[-1]) 332 445 ax_specific.set_ylim(bottom=0.0) 333 - ax_specific.set_title(f"Tonality T = {result.tonality:.2f} tu_HMS") 446 + ax_specific.set_title(f"ECMA-418-2 tonality T = {result.tonality:.2f} tu_HMS") 334 447 ax_specific.grid(True, alpha=0.3) 335 448 336 449 if not two_panel: ··· 339 452 ax_time = cast("Axes", axes[1]) 340 453 time = np.asarray(result.time, dtype=np.float64) 341 454 tvt = np.asarray(result.tonality_vs_time, dtype=np.float64) 342 - ax_time.plot(time, tvt, color="#9467bd", label="T(l)") 455 + ax_time.plot(time, tvt, color=_C_QUATERNARY, label="T(l)") 343 456 ax_time.set_xlabel("Time [s]") 344 457 ax_time.set_ylabel("Tonality T [tu_HMS]") 345 458 ax_time.set_ylim(bottom=0.0) ··· 373 486 else: 374 487 ax_time = cast("Axes", ax) 375 488 489 + # Roughness's per-metric identity color is brown across the documentation 490 + # figures (tonality is red); kept literal on purpose, see the module 491 + # color-constant note. 376 492 kwargs.setdefault("color", "#8c564b") 377 493 ax_time.plot(time, rvt, **kwargs) 378 - ax_time.fill_between(time, rvt, color="#8c564b", alpha=0.25) 494 + ax_time.fill_between(time, rvt, color=kwargs["color"], alpha=0.25) 379 495 ax_time.set_xlabel("Time [s]") 380 496 ax_time.set_ylabel("Roughness R [asper]") 381 497 ax_time.set_ylim(bottom=0.0) 382 - ax_time.set_title(f"Roughness R = {result.roughness:.2f} asper") 498 + ax_time.set_title(f"ECMA-418-2 roughness R = {result.roughness:.2f} asper") 383 499 ax_time.grid(True, alpha=0.3) 384 500 385 501 if not two_panel: ··· 412 528 ax = ax if ax is not None else _new_axes() 413 529 mti = np.asarray(result.mti, dtype=np.float64) 414 530 positions = np.arange(mti.size) 415 - kwargs.setdefault("color", "#1f77b4") 531 + kwargs.setdefault("color", _C_PRIMARY) 416 532 ax.bar(positions, mti, **kwargs) 417 - ax.set_xticks(positions) 418 533 if mti.size == len(_STI_BAND_CENTERS): 419 - ax.set_xticklabels([_format_freq(f) for f in _STI_BAND_CENTERS]) 420 - ax.set_xlabel("Octave band [Hz]") 534 + _band_axis(ax, np.asarray(_STI_BAND_CENTERS)) 421 535 else: 422 - ax.set_xlabel("Octave band") 536 + ax.set_xticks(positions) 537 + ax.set_xlabel("Band") 423 538 ax.set_ylabel("Modulation transfer index MTI") 424 539 ax.set_ylim(0.0, 1.0) 425 - ax.set_title(f"STI = {result.sti:.2f} (rating {result.rating})") 540 + ax.set_title(f"IEC 60268-16 STI = {result.sti:.2f} (rating {result.rating})") 426 541 ax.grid(True, axis="y", alpha=0.3) 427 542 return ax 428 543 ··· 439 554 freqs = np.asarray(result.frequencies, dtype=np.float64) 440 555 audibility = np.asarray(result.band_audibility, dtype=np.float64) 441 556 contribution = audibility * np.asarray(result.band_importance, dtype=np.float64) 442 - positions = np.arange(freqs.size) 443 - ax.bar(positions, audibility, color="#c6dbef", label="Band audibility $A_i$") 444 - kwargs.setdefault("color", "#1f77b4") 557 + positions = _band_axis(ax, freqs) 558 + ax.bar(positions, audibility, color=_C_PRIMARY_LIGHT, 559 + label="Band audibility $A_i$") 560 + kwargs.setdefault("color", _C_PRIMARY) 445 561 # A fully masked speech signal (SII = 0) has an all-zero contribution; 446 562 # keep the zero bars rather than dividing 0/0 into NaN. 447 563 peak = float(contribution.max()) if contribution.size else 0.0 448 564 scaled = contribution / peak if peak > 0.0 else contribution 449 565 ax.bar(positions, scaled, width=0.5, 450 566 label=r"Importance-weighted $I_i A_i$ (scaled)", **kwargs) 451 - ax.set_xticks(positions) 452 - ax.set_xticklabels([_format_freq(f) for f in freqs], rotation=45, ha="right") 453 - ax.set_xlabel("One-third-octave band [Hz]") 454 567 ax.set_ylabel("Band audibility") 455 568 ax.set_ylim(0.0, 1.0) 456 - ax.set_title(f"SII = {result.sii:.3f}") 569 + ax.set_title(f"ANSI S3.5 SII = {result.sii:.3f}") 457 570 ax.legend(loc=_LEGEND_UPPER_RIGHT, fontsize="small") 458 571 ax.grid(True, axis="y", alpha=0.3) 459 572 return ax ··· 470 583 reference: np.ndarray, 471 584 *, 472 585 impact: bool, 473 - rating: int, 474 - unfavourable_sum: float, 475 586 title: str, 476 587 ylabel: str, 588 + measured_label: str = "Measured", 589 + ylim: tuple[float, float] | None = None, 477 590 ax: Axes | None, 478 591 **kwargs: Any, 479 592 ) -> Axes: 480 - """Shared renderer for airborne and impact rating figures. 593 + """Shared renderer for the shifted-reference rating figures. 481 594 482 595 Draws the measured curve against the shifted reference and shades the 483 596 unfavourable deviations: where the reference exceeds the measurement 484 - for airborne insulation (higher is better) and where the measurement 485 - exceeds the reference for impact sound (lower is better). Extra 486 - keyword arguments style the measured curve (its primary artist). 597 + for airborne insulation and absorption (higher is better) and where the 598 + measurement exceeds the reference for impact sound (lower is better). 599 + Extra keyword arguments style the measured curve (its primary artist). 487 600 """ 488 601 ax = ax if ax is not None else _new_axes() 489 - kwargs.setdefault("color", "#1f77b4") 490 - kwargs.setdefault("label", "Measured") 602 + kwargs.setdefault("color", _C_PRIMARY) 603 + kwargs.setdefault("label", measured_label) 491 604 ax.plot(band_centers, measured, "o-", **kwargs) 492 - ax.plot(band_centers, reference, "s--", color="#d62728", label="Shifted reference") 605 + ax.plot(band_centers, reference, "s--", color=_C_REFERENCE, 606 + label="Shifted reference") 493 607 unfavourable = _unfavourable_mask(measured, reference, impact) 494 608 ax.fill_between( 495 609 band_centers, 496 610 measured, 497 611 reference, 498 612 where=unfavourable.tolist(), 499 - color="#ff7f0e", 613 + color=_C_SECONDARY, 500 614 alpha=0.4, 501 615 label="Unfavourable deviations", 502 616 interpolate=True, 503 617 ) 504 618 _freq_axis(ax, band_centers) 505 619 ax.set_ylabel(ylabel) 506 - ax.set_title(f"{title} = {rating} dB (Sigma unfav. = {unfavourable_sum:.1f} dB)") 620 + if ylim is not None: 621 + ax.set_ylim(*ylim) 622 + ax.set_title(title) 507 623 ax.grid(True, which="both", alpha=0.3) 508 624 ax.legend(loc="best", fontsize="small") 509 625 return ax ··· 523 639 return np.asarray(measured < reference) 524 640 525 641 526 - def plot_weighted_rating(result: Any, ax: Axes | None = None, **kwargs: Any) -> Axes: 527 - """Airborne rating curve vs shifted reference (ISO 717-1).""" 642 + def plot_weighted_rating( 643 + result: "WeightedRatingResult", ax: Axes | None = None, **kwargs: Any 644 + ) -> Axes: 645 + """Airborne rating curve vs shifted reference (ISO 717-1). 646 + 647 + :param result: A :class:`~phonometry.insulation.WeightedRatingResult`. 648 + :param ax: Existing axes, or ``None`` to create a figure. 649 + :param kwargs: Forwarded to the measured-curve ``plot`` call. 650 + :return: The axes. 651 + """ 528 652 _require_rating_curve(result) 529 653 return _plot_rating( 530 654 np.asarray(result.band_centers, dtype=np.float64), 531 655 np.asarray(result.measured, dtype=np.float64), 532 656 np.asarray(result.shifted_reference, dtype=np.float64), 533 657 impact=False, 534 - rating=result.rating, 535 - unfavourable_sum=result.unfavourable_sum, 536 - title=f"Rw (C={result.c:+d}; Ctr={result.ctr:+d})", 658 + title=( 659 + f"ISO 717-1 Rw (C={result.c:+d}; Ctr={result.ctr:+d}) = " 660 + f"{result.rating} dB (Sigma unfav. = {result.unfavourable_sum:.1f} dB)" 661 + ), 537 662 ylabel="Sound reduction index [dB]", 538 663 ax=ax, 539 664 **kwargs, 540 665 ) 541 666 542 667 543 - def plot_impact_rating(result: Any, ax: Axes | None = None, **kwargs: Any) -> Axes: 668 + def plot_impact_rating( 669 + result: "ImpactRatingResult", ax: Axes | None = None, **kwargs: Any 670 + ) -> Axes: 544 671 """Impact rating curve vs shifted reference (ISO 717-2). 545 672 546 673 The drawn shifted-reference curve is the normatively honest ``ref - ··· 548 675 *minus 5 dB* (Clause 4.3.2), so the plot marks the 500 Hz read value on 549 676 the (undistorted) curve and annotates the -5 dB reduction rather than 550 677 pulling the curve down to the rating. 678 + 679 + :param result: An :class:`~phonometry.insulation.ImpactRatingResult`. 680 + :param ax: Existing axes, or ``None`` to create a figure. 681 + :param kwargs: Forwarded to the measured-curve ``plot`` call. 682 + :return: The axes. 551 683 """ 552 684 _require_rating_curve(result) 553 685 band_centers = np.asarray(result.band_centers, dtype=np.float64) ··· 557 689 np.asarray(result.measured, dtype=np.float64), 558 690 reference, 559 691 impact=True, 560 - rating=result.rating, 561 - unfavourable_sum=result.unfavourable_sum, 562 692 # The rated quantity depends on the input (Ln,w, L'n,w or L'nT,w); 563 693 # the dataclass does not carry which, so the figure uses the neutral 564 - # "Impact rating" label rather than hard-coding one specific symbol. 565 - title=f"Impact rating (CI={result.ci:+d})", 694 + # "impact rating" label rather than hard-coding one specific symbol. 695 + title=( 696 + f"ISO 717-2 impact rating (CI={result.ci:+d}) = {result.rating} dB" 697 + f" (Sigma unfav. = {result.unfavourable_sum:.1f} dB)" 698 + ), 566 699 ylabel="Impact sound pressure level [dB]", 567 700 ax=ax, 568 701 **kwargs, ··· 571 704 return ax 572 705 573 706 574 - def plot_weighted_absorption(result: Any, ax: Axes | None = None, **kwargs: Any) -> Axes: 707 + def plot_weighted_absorption( 708 + result: "AbsorptionRatingResult", ax: Axes | None = None, **kwargs: Any 709 + ) -> Axes: 575 710 """Practical absorption curve vs the shifted reference (ISO 11654:1997). 576 711 577 712 Draws the practical coefficients ``alpha_p`` against the shifted reference 578 713 curve and shades the unfavourable deviations (measured below the shifted 579 - reference, Clause 4.2), reusing :func:`_unfavourable_mask` and the shared 580 - frequency axis. Extra keyword arguments style the measured curve. 714 + reference, Clause 4.2) through the shared rating renderer. 715 + 716 + :param result: An 717 + :class:`~phonometry.absorption_rating.AbsorptionRatingResult`. 718 + :param ax: Existing axes, or ``None`` to create a figure. 719 + :param kwargs: Forwarded to the measured-curve ``plot`` call. 720 + :return: The axes. 581 721 """ 582 - ax = ax if ax is not None else _new_axes() 583 - band_centers = np.asarray(result.band_centers, dtype=np.float64) 584 - measured = np.asarray(result.measured, dtype=np.float64) 585 - reference = np.asarray(result.shifted_reference, dtype=np.float64) 586 - kwargs.setdefault("color", "#1f77b4") 587 - kwargs.setdefault("label", "Practical alpha_p") 588 - ax.plot(band_centers, measured, "o-", **kwargs) 589 - ax.plot(band_centers, reference, "s--", color="#d62728", label="Shifted reference") 590 - unfavourable = _unfavourable_mask(measured, reference, impact=False) 591 - ax.fill_between( 592 - band_centers, 593 - measured, 594 - reference, 595 - where=unfavourable.tolist(), 596 - color="#ff7f0e", 597 - alpha=0.4, 598 - label="Unfavourable deviations", 599 - interpolate=True, 722 + return _plot_rating( 723 + np.asarray(result.band_centers, dtype=np.float64), 724 + np.asarray(result.measured, dtype=np.float64), 725 + np.asarray(result.shifted_reference, dtype=np.float64), 726 + impact=False, 727 + title=( 728 + f"ISO 11654 alpha_w = {result.rating_label} " 729 + f"(class {result.absorption_class}, " 730 + f"Sigma unfav. = {result.unfavourable_sum:.2f})" 731 + ), 732 + ylabel="Sound absorption coefficient", 733 + measured_label="Practical alpha_p", 734 + ylim=(0.0, 1.05), 735 + ax=ax, 736 + **kwargs, 600 737 ) 601 - _freq_axis(ax, band_centers) 602 - ax.set_ylabel("Sound absorption coefficient") 603 - ax.set_ylim(0.0, 1.05) 604 - ax.set_title( 605 - f"alpha_w = {result.rating_label} (class {result.absorption_class}, " 606 - f"Sigma unfav. = {result.unfavourable_sum:.2f})" 607 - ) 608 - ax.grid(True, which="both", alpha=0.3) 609 - ax.legend(loc="best", fontsize="small") 610 - return ax 611 738 612 739 613 740 def _annotate_impact_500( ··· 626 753 [read_value], 627 754 marker="D", 628 755 ls="", 629 - color="#d62728", 756 + color=_C_REFERENCE, 630 757 ms=9, 631 758 mfc="none", 632 759 mew=1.6, ··· 642 769 textcoords="offset points", 643 770 ha="center", 644 771 fontsize="small", 645 - arrowprops={"arrowstyle": "->", "color": "#555555"}, 772 + arrowprops={"arrowstyle": "->", "color": _C_EDGE}, 646 773 ) 647 774 ax.legend(loc="best", fontsize="small") 648 775 649 776 650 - def _require_rating_curve(result: Any) -> None: 777 + def _require_rating_curve( 778 + result: "WeightedRatingResult | ImpactRatingResult", 779 + ) -> None: 651 780 if ( 652 781 result.band_centers is None 653 782 or result.measured is None ··· 660 789 661 790 662 791 def plot_facade_insulation( 663 - result: Any, ax: Axes | None = None, **kwargs: Any 792 + result: "FacadeInsulationResult", ax: Axes | None = None, **kwargs: Any 664 793 ) -> Axes: 665 794 """Per-band façade sound-insulation profile (ISO 16283-3). 666 795 ··· 672 801 :param result: A façade result exposing ``d_2m``, ``d_2m_nt``, 673 802 ``d_2m_n``, ``r_prime`` and (optionally) ``frequencies``. 674 803 :param ax: Existing axes, or ``None`` to create a figure. 804 + :param kwargs: Forwarded to the primary ``D2m,nT`` curve ``plot`` call. 675 805 :return: The axes. 676 806 """ 677 807 ax = ax if ax is not None else _new_axes() ··· 753 883 754 884 _draw_decay_times(ax_times, positions, result, **kwargs) 755 885 ax_times.set_ylabel("Reverberation time [s]") 756 - ax_times.set_title("Decay times and clarity") 757 - ax_times.set_xticks(positions) 758 - ax_times.set_xticklabels(labels, rotation=45, ha="right") 886 + ax_times.set_title("ISO 3382 decay times and clarity") 887 + _band_axis(ax_times, labels, xlabel=None) 759 888 ax_times.grid(True, axis="y", alpha=0.3) 760 889 ax_times.legend(loc="best", fontsize="small") 761 890 ··· 769 898 positions, 770 899 np.asarray(result.c50, dtype=np.float64), 771 900 "o-", 772 - color="#2ca02c", 901 + color=_C_TERTIARY, 773 902 label="C50", 774 903 ) 775 904 ax_clarity.plot( 776 905 positions, 777 906 np.asarray(result.c80, dtype=np.float64), 778 907 "s--", 779 - color="#9467bd", 908 + color=_C_QUATERNARY, 780 909 label="C80", 781 910 ) 782 911 ax_clarity.set_ylabel("Clarity [dB]") 783 - ax_clarity.set_xticks(positions) 784 - ax_clarity.set_xticklabels(labels, rotation=45, ha="right") 785 - ax_clarity.set_xlabel("Frequency [Hz]" if use_freq_axis else "Band") 912 + _band_axis( 913 + ax_clarity, labels, xlabel="Frequency [Hz]" if use_freq_axis else "Band" 914 + ) 786 915 ax_clarity.grid(True, alpha=0.3) 787 916 ax_clarity.legend(loc="best", fontsize="small") 788 917 return axes ··· 794 923 """Grouped EDT/T20/T30 bars, invalid bands hatched and greyed.""" 795 924 width = 0.27 796 925 series = ( 797 - ("EDT", result.edt, result.edt_valid, -width, "#1f77b4"), 798 - ("T20", result.t20, result.t20_valid, 0.0, "#ff7f0e"), 799 - ("T30", result.t30, result.t30_valid, width, "#2ca02c"), 926 + ("EDT", result.edt, result.edt_valid, -width, _C_PRIMARY), 927 + ("T20", result.t20, result.t20_valid, 0.0, _C_SECONDARY), 928 + ("T30", result.t30, result.t30_valid, width, _C_TERTIARY), 800 929 ) 801 930 for label, values, valid, offset, color in series: 802 931 vals = np.asarray(values, dtype=np.float64) 803 932 valid_arr = np.asarray(valid, dtype=bool) 804 - colors = [color if v else "#bbbbbb" for v in valid_arr] 805 - hatches = np.where(valid_arr, "", "//") 933 + colors = [color if v else _C_MUTED for v in valid_arr] 806 934 # Merge per-series defaults with the user kwargs (user wins) freshly 807 935 # each iteration so an overriding label/color is not frozen by the 808 936 # first band group. ··· 813 941 width=width, 814 942 **bar_kwargs, 815 943 ) 816 - for bar, hatch in zip(bars, hatches, strict=True): 817 - if hatch: 818 - bar.set_hatch(hatch) 819 - bar.set_edgecolor("#555555") 944 + _hatch_invalid(bars, ~valid_arr) 820 945 821 946 822 947 # --------------------------------------------------------------------------- ··· 824 949 # --------------------------------------------------------------------------- 825 950 826 951 827 - def plot_sound_power(result: Any, ax: Axes | None = None, **kwargs: Any) -> Axes: 952 + def plot_sound_power( 953 + result: ( 954 + "SoundPowerResult | ReverberationSoundPowerResult" 955 + " | SoundPowerIntensityResult | Any" 956 + ), 957 + ax: Axes | None = None, 958 + **kwargs: Any, 959 + ) -> Axes: 828 960 """Sound power level spectrum with the A-weighted total annotated. 829 961 830 962 Works for :class:`~phonometry.sound_power.SoundPowerResult`, ··· 837 969 ``sound_power_level``, ``sound_power_level_a`` and (optionally) 838 970 ``frequencies`` and ``negative_band``. 839 971 :param ax: Existing axes, or ``None`` to create a figure. 972 + :param kwargs: Forwarded to the band :meth:`~matplotlib.axes.Axes.bar`. 840 973 :return: The axes. 841 974 """ 842 975 ax = ax if ax is not None else _new_axes() 843 976 lw = np.asarray(result.sound_power_level, dtype=np.float64) 844 977 n = lw.size 845 - # Band bars are drawn on evenly spaced categorical positions (a linear 846 - # frequency-band axis) so they stay legible; the tick labels carry the 847 - # band centres when available. 848 - positions = np.arange(n, dtype=np.float64) 849 978 freqs = getattr(result, "frequencies", None) 850 979 if freqs is None: 851 - labels = [f"Band {i + 1}" for i in range(n)] 852 - xlabel = "Band" 980 + positions = _band_axis( 981 + ax, [f"Band {i + 1}" for i in range(n)], xlabel="Band" 982 + ) 853 983 else: 854 - labels = [_format_freq(f) for f in np.asarray(freqs, dtype=np.float64)] 855 - xlabel = "Frequency [Hz]" 984 + positions = _band_axis(ax, np.asarray(freqs, dtype=np.float64)) 856 985 857 986 # ``negative_band`` (ISO 9614-2) and ``not_applicable_band`` (ISO 9614-3) 858 987 # both flag bands whose net power is non-positive and therefore unusable. ··· 864 993 if negative is not None 865 994 else np.zeros(n, dtype=bool) 866 995 ) 867 - colors = ["#bbbbbb" if b else "#1f77b4" for b in neg] 996 + colors = [_C_MUTED if b else _C_PRIMARY for b in neg] 868 997 kwargs.setdefault("color", colors) 869 998 bars = ax.bar(positions, np.nan_to_num(lw), **kwargs) 870 - for bar, is_neg in zip(bars, neg, strict=True): 871 - if is_neg: 872 - bar.set_hatch("//") 873 - bar.set_edgecolor("#555555") 999 + _hatch_invalid(bars, neg) 874 1000 875 - ax.set_xticks(positions) 876 - ax.set_xticklabels(labels, rotation=45, ha="right") 877 - ax.set_xlabel(xlabel) 878 1001 ax.set_ylabel("Sound power level LW [dB]") 1002 + designation = _sound_power_designation(result) 879 1003 lwa = float(result.sound_power_level_a) 880 1004 if np.isfinite(lwa): 881 - ax.set_title(f"Sound power spectrum (LWA = {lwa:.1f} dB(A))") 1005 + ax.set_title(f"{designation} sound power spectrum (LWA = {lwa:.1f} dB(A))") 882 1006 else: 883 - ax.set_title("Sound power spectrum") 1007 + ax.set_title(f"{designation} sound power spectrum") 884 1008 if np.any(neg): 885 - ax.plot([], [], color="#bbbbbb", marker="s", ls="", label="Non-positive band") 1009 + ax.plot([], [], color=_C_MUTED, marker="s", ls="", label="Non-positive band") 886 1010 ax.legend(loc="best", fontsize="small") 887 1011 ax.grid(True, axis="y", alpha=0.3) 888 1012 return ax 889 1013 890 1014 1015 + def _sound_power_designation(result: Any) -> str: 1016 + """The standard designation matching a sound-power result's method. 1017 + 1018 + Distinguishes the reverberation-room (ISO 3741) and intensity (ISO 9614) 1019 + determinations by their result types; the enveloping-surface pressure 1020 + methods (:class:`~phonometry.sound_power.SoundPowerResult` and any other 1021 + duck-typed result) fall back to ISO 3744/3746. 1022 + """ 1023 + from .sound_power_intensity import SoundPowerIntensityResult 1024 + from .sound_power_reverberation import ReverberationSoundPowerResult 1025 + 1026 + if isinstance(result, ReverberationSoundPowerResult): 1027 + return "ISO 3741" 1028 + if isinstance(result, SoundPowerIntensityResult): 1029 + return "ISO 9614" 1030 + return "ISO 3744/3746" 1031 + 1032 + 891 1033 # --------------------------------------------------------------------------- 892 1034 # Sound intensity (ISO 9614 / IEC 61043) 893 1035 # --------------------------------------------------------------------------- ··· 918 1060 :param result: An :class:`~phonometry.intensity.IntensityResult` with 919 1061 per-band data (obtained by requesting a band ``fraction``). 920 1062 :param ax: Existing axes, or ``None`` to create a figure. 1063 + :param kwargs: Forwarded to the pressure-level curve ``plot`` call. 921 1064 :return: The axes. 922 1065 :raises ValueError: If the result carries no per-band data. 923 1066 """ ··· 932 1075 li = np.asarray(result.intensity_level, dtype=np.float64) 933 1076 index = np.asarray(result.pressure_intensity_index, dtype=np.float64) 934 1077 935 - kwargs.setdefault("color", "#1f77b4") 1078 + kwargs.setdefault("color", _C_PRIMARY) 936 1079 kwargs.setdefault("label", "Pressure level Lp") 937 1080 ax.plot(freqs, lp, "o-", **kwargs) 938 - ax.plot(freqs, li, "s--", color="#d62728", label="Intensity level LI") 1081 + ax.plot(freqs, li, "s--", color=_C_REFERENCE, label="Intensity level LI") 939 1082 _freq_axis(ax, freqs) 940 1083 ax.set_ylabel("Level [dB]") 941 1084 ax.grid(True, which="both", alpha=0.3) ··· 945 1088 freqs, 946 1089 index, 947 1090 width=_bar_width(freqs), 948 - color="#2ca02c", 1091 + color=_C_TERTIARY, 949 1092 alpha=0.25, 950 1093 label="δpI = Lp - LI", 951 1094 ) ··· 955 1098 tlines, tlabels = twin.get_legend_handles_labels() 956 1099 ax.legend(lines + tlines, labels + tlabels, loc="best", fontsize="small") 957 1100 ax.set_title( 958 - f"Lp vs LI (total δpI = {result.total_pressure_intensity_index:.1f} dB)" 1101 + "ISO 9614 Lp vs LI " 1102 + f"(total δpI = {result.total_pressure_intensity_index:.1f} dB)" 959 1103 ) 960 1104 return ax 961 1105 ··· 980 1124 ax = ax if ax is not None else _new_axes() 981 1125 time = np.asarray(result.time, dtype=np.float64) 982 1126 level = np.asarray(result.level, dtype=np.float64) 983 - kwargs.setdefault("color", "#1f77b4") 1127 + kwargs.setdefault("color", _C_PRIMARY) 984 1128 kwargs.setdefault("label", "Schroeder decay") 985 1129 ax.plot(time, level, **kwargs) 986 1130 ··· 999 1143 ax.set_ylim(top=3.0) 1000 1144 ax.set_xlim(left=0.0, right=float(time[-1]) if time.size else None) 1001 1145 band = result.band 1002 - title = "Schroeder decay curve" 1146 + title = "ISO 3382 Schroeder decay curve" 1003 1147 if band is not None: 1004 1148 title += f" ({_format_freq(float(band))} Hz band)" 1005 1149 ax.set_title(title) ··· 1061 1205 total = float(energy[0]) if energy.size else 0.0 1062 1206 edc_db = 10.0 * np.log10(np.maximum(energy, tiny) / (total if total > 0.0 else 1.0)) 1063 1207 1064 - color = kwargs.pop("color", "#1f77b4") 1208 + color = kwargs.pop("color", _C_PRIMARY) 1065 1209 1066 1210 def _decay(axd: Axes) -> None: 1067 - axd.plot(time, env_db, color="#9ecae1", lw=0.8, label="Log-magnitude envelope") 1068 - axd.plot(time, edc_db, color="#d62728", lw=1.8, label="Schroeder decay") 1211 + axd.plot(time, env_db, color=_C_PRIMARY_LIGHT, lw=0.8, 1212 + label="Log-magnitude envelope") 1213 + axd.plot(time, edc_db, color=_C_REFERENCE, lw=1.8, label="Schroeder decay") 1069 1214 axd.set_xlabel(xlabel) 1070 1215 axd.set_ylabel("Level re peak [dB]") 1071 1216 axd.set_ylim(bottom=-80.0, top=5.0) ··· 1076 1221 1077 1222 if ax is not None: 1078 1223 _decay(ax) 1079 - ax.set_title(f"Impulse response ({result.method})") 1224 + ax.set_title(f"ISO 18233 impulse response ({result.method})") 1080 1225 return ax 1081 1226 1082 1227 axes = _new_axes_column(2, sharex=True, figsize=(8.0, 6.0)) 1083 1228 axes[0].plot(time, h / norm, color=color, lw=0.8, **kwargs) 1084 1229 axes[0].set_ylabel("Amplitude (norm.)") 1085 - axes[0].set_title(f"Impulse response ({result.method})") 1230 + axes[0].set_title(f"ISO 18233 impulse response ({result.method})") 1086 1231 axes[0].grid(True, alpha=0.3) 1087 1232 _decay(axes[1]) 1088 1233 return axes ··· 1117 1262 if n == 0: 1118 1263 raise ValueError("excitation signal is empty; nothing to plot.") 1119 1264 t = np.arange(n) / float(fs) 1120 - color = kwargs.pop("color", "#1f77b4") 1265 + color = kwargs.pop("color", _C_PRIMARY) 1121 1266 1122 1267 two_panel = ax is None 1123 1268 if two_panel: ··· 1132 1277 ax_time.set_xlabel("Sample") 1133 1278 ax_time.set_ylabel("Amplitude") 1134 1279 ax_time.set_ylim(-1.4, 1.4) 1135 - ax_time.set_title(f"MLS excitation (first {show} of {n} samples)") 1280 + ax_time.set_title( 1281 + f"ISO 18233 MLS excitation (first {show} of {n} samples)" 1282 + ) 1136 1283 ax_time.grid(True, alpha=0.3) 1137 1284 if not two_panel: 1138 1285 return ax_time ··· 1143 1290 denom = float(np.median(ac)) if ac.size else 1.0 1144 1291 ax_f.semilogx(freqs[1:], 20.0 * np.log10( 1145 1292 np.maximum(ac, 1e-10) / (denom if denom > 0.0 else 1.0)), 1146 - color="#d62728", lw=0.8) 1293 + color=_C_REFERENCE, lw=0.8) 1147 1294 ax_f.set_xlabel("Frequency [Hz]") 1148 1295 ax_f.set_ylabel("Magnitude [dB]") 1149 1296 ax_f.set_title("Magnitude spectrum (flat)") ··· 1154 1301 ax_time.plot(t, x, color=color, lw=0.6, **kwargs) 1155 1302 ax_time.set_xlabel("Time [s]") 1156 1303 ax_time.set_ylabel("Amplitude") 1157 - ax_time.set_title("Exponential sine sweep") 1304 + ax_time.set_title("ISO 18233 exponential sine sweep") 1158 1305 ax_time.grid(True, alpha=0.3) 1159 1306 if not two_panel: 1160 1307 return ax_time ··· 1173 1320 1174 1321 1175 1322 def plot_scattering_coefficient( 1176 - result: Any, ax: Axes | None = None, **kwargs: Any 1323 + result: "ScatteringResult", ax: Axes | None = None, **kwargs: Any 1177 1324 ) -> Axes: 1178 1325 """Random-incidence scattering coefficient ``s`` versus frequency. 1179 1326 1180 1327 :param result: A :class:`~phonometry.scattering_diffusion.ScatteringResult` 1181 1328 exposing ``frequencies`` and ``scattering``. 1182 1329 :param ax: Existing axes, or ``None`` to create a figure. 1330 + :param kwargs: Forwarded to the coefficient curve ``plot`` call. 1183 1331 :return: The axes. 1184 1332 """ 1185 1333 ax = ax if ax is not None else _new_axes() 1186 1334 freqs = np.asarray(result.frequencies, dtype=np.float64) 1187 1335 s = np.asarray(result.scattering, dtype=np.float64) 1188 1336 kwargs.setdefault("marker", "o") 1189 - kwargs.setdefault("color", "#1f77b4") 1337 + kwargs.setdefault("color", _C_PRIMARY) 1190 1338 ax.plot(freqs, s, **kwargs) 1191 1339 _freq_axis(ax, freqs) 1192 1340 ax.set_ylabel("Scattering coefficient s") ··· 1200 1348 1201 1349 1202 1350 def plot_diffusion_polar( 1203 - result: Any, ax: Axes | None = None, **kwargs: Any 1351 + result: "DiffusionResult", ax: Axes | None = None, **kwargs: Any 1204 1352 ) -> Axes: 1205 1353 """Polar reflected-level response with the diffusion coefficient annotated. 1206 1354 1207 1355 :param result: A :class:`~phonometry.scattering_diffusion.DiffusionResult` 1208 1356 exposing ``angles`` (degrees), ``levels`` (dB) and ``coefficient``. 1209 1357 :param ax: Existing (ideally polar) axes, or ``None`` to create a polar one. 1358 + :param kwargs: Forwarded to the reflected-level curve ``plot`` call. 1210 1359 :return: The polar axes. 1211 1360 """ 1212 1361 if ax is None: ··· 1215 1364 angles = np.radians(np.asarray(result.angles, dtype=np.float64)) 1216 1365 levels = np.asarray(result.levels, dtype=np.float64) 1217 1366 kwargs.setdefault("marker", "o") 1218 - kwargs.setdefault("color", "#1f77b4") 1367 + kwargs.setdefault("color", _C_PRIMARY) 1219 1368 ax.plot(angles, levels, **kwargs) 1220 1369 ax.fill(angles, levels, alpha=0.15, color=kwargs["color"]) 1221 1370 ax.set_title( ··· 1226 1375 1227 1376 1228 1377 def plot_insitu_absorption( 1229 - result: Any, ax: Axes | None = None, **kwargs: Any 1378 + result: "InsituAbsorptionResult", ax: Axes | None = None, **kwargs: Any 1230 1379 ) -> Axes: 1231 1380 """In-situ one-third-octave absorption spectrum ``alpha(f)``. 1232 1381 ··· 1234 1383 :class:`~phonometry.road_absorption.InsituAbsorptionResult` exposing 1235 1384 ``frequencies`` and ``absorption``. 1236 1385 :param ax: Existing axes, or ``None`` to create a figure. 1386 + :param kwargs: Forwarded to the absorption :meth:`~matplotlib.axes.Axes.bar`. 1237 1387 :return: The axes. 1238 1388 """ 1239 1389 ax = ax if ax is not None else _new_axes() 1240 1390 freqs = np.asarray(result.frequencies, dtype=np.float64) 1241 1391 alpha = np.asarray(result.absorption, dtype=np.float64) 1242 - positions = np.arange(freqs.size, dtype=np.float64) 1243 - kwargs.setdefault("color", "#1f77b4") 1392 + positions = _band_axis(ax, freqs) 1393 + kwargs.setdefault("color", _C_PRIMARY) 1244 1394 ax.bar(positions, np.nan_to_num(alpha), **kwargs) 1245 - ax.set_xticks(positions) 1246 - ax.set_xticklabels([_format_freq(f) for f in freqs], rotation=45, ha="right") 1247 - ax.set_xlabel("Frequency [Hz]") 1248 1395 ax.set_ylabel("Absorption coefficient") 1249 1396 ax.set_ylim(0.0, 1.0) 1250 1397 ax.set_title("In-situ road-surface absorption (ISO 13472-1)") ··· 1258 1405 1259 1406 1260 1407 def plot_vibration_weighting( 1261 - result: Any, ax: Axes | None = None, **kwargs: Any 1408 + result: "WeightingResponse", ax: Axes | None = None, **kwargs: Any 1262 1409 ) -> Axes: 1263 1410 """Frequency-weighting factor (dB) versus frequency (ISO 8041-1). 1264 1411 ··· 1266 1413 :class:`~phonometry.human_vibration.WeightingResponse` exposing 1267 1414 ``name``, ``frequencies`` and ``magnitude_db``. 1268 1415 :param ax: Existing axes, or ``None`` to create a figure. 1416 + :param kwargs: Forwarded to the weighting curve ``plot`` call. 1269 1417 :return: The axes. 1270 1418 """ 1271 1419 ax = ax if ax is not None else _new_axes() 1272 1420 freqs = np.asarray(result.frequencies, dtype=np.float64) 1273 1421 mag_db = np.asarray(result.magnitude_db, dtype=np.float64) 1274 - kwargs.setdefault("color", "#1f77b4") 1422 + kwargs.setdefault("color", _C_PRIMARY) 1275 1423 ax.semilogx(freqs, mag_db, **kwargs) 1276 1424 ax.set_xlabel("Frequency [Hz]") 1277 1425 ax.set_ylabel("Weighting factor [dB]") ··· 1281 1429 1282 1430 1283 1431 def plot_weighted_spectrum( 1284 - result: Any, ax: Axes | None = None, **kwargs: Any 1432 + result: "WeightedSpectrum", ax: Axes | None = None, **kwargs: Any 1285 1433 ) -> Axes: 1286 1434 """Unweighted vs weighted one-third-octave acceleration spectrum. 1287 1435 ··· 1293 1441 ``frequencies``, ``band_accelerations``, ``weighted``, ``overall`` and 1294 1442 ``weighting_name``. 1295 1443 :param ax: Existing axes, or ``None`` to create a figure. 1444 + :param kwargs: Forwarded to the weighted (primary) bars. 1296 1445 :return: The axes. 1297 1446 """ 1298 1447 ax = ax if ax is not None else _new_axes() 1299 1448 freqs = np.asarray(result.frequencies, dtype=np.float64) 1300 1449 raw = np.asarray(result.band_accelerations, dtype=np.float64) 1301 1450 weighted = np.asarray(result.weighted, dtype=np.float64) 1302 - positions = np.arange(freqs.size, dtype=np.float64) 1451 + positions = _band_axis(ax, freqs) 1303 1452 width = 0.4 1304 1453 # The weighted bars are the primary artist; forward user kwargs there. 1305 - kwargs.setdefault("color", "#1f77b4") 1454 + kwargs.setdefault("color", _C_PRIMARY) 1306 1455 ax.bar( 1307 - positions - width / 2, raw, width, color="#bbbbbb", label="Unweighted $a_i$" 1456 + positions - width / 2, raw, width, color=_C_MUTED, label="Unweighted $a_i$" 1308 1457 ) 1309 1458 ax.bar( 1310 1459 positions + width / 2, ··· 1313 1462 label=f"Weighted $W_i a_i$ ({result.weighting_name})", 1314 1463 **kwargs, 1315 1464 ) 1316 - ax.set_xticks(positions) 1317 - ax.set_xticklabels([_format_freq(f) for f in freqs], rotation=45, ha="right") 1318 - ax.set_xlabel("Frequency [Hz]") 1319 1465 ax.set_ylabel(r"r.m.s. acceleration [m/s$^2$]") 1466 + # Wh is the hand-arm weighting of ISO 5349-1; the others (Wk, Wd, Wm...) 1467 + # are the whole-body weightings of ISO 2631. 1468 + designation = "ISO 5349-1" if str(result.weighting_name) == "Wh" else "ISO 2631" 1320 1469 ax.set_title( 1321 - f"Weighted acceleration spectrum ($a_w$ = {float(result.overall):.3f} " 1322 - r"m/s$^2$)" 1470 + f"{designation} weighted acceleration spectrum " 1471 + f"($a_w$ = {float(result.overall):.3f} " r"m/s$^2$)" 1323 1472 ) 1324 1473 ax.legend(loc="best", fontsize="small") 1325 1474 ax.grid(True, axis="y", alpha=0.3) ··· 1327 1476 1328 1477 1329 1478 def plot_daily_exposure( 1330 - result: Any, ax: Axes | None = None, **kwargs: Any 1479 + result: "DailyVibrationExposure", ax: Axes | None = None, **kwargs: Any 1331 1480 ) -> Axes: 1332 1481 """Partial daily exposures against the EAV / ELV (Directive 2002/44/EC). 1333 1482 ··· 1338 1487 :class:`~phonometry.human_vibration.DailyVibrationExposure` exposing 1339 1488 ``labels``, ``partials``, ``a8`` and ``assessment``. 1340 1489 :param ax: Existing axes, or ``None`` to create a figure. 1490 + :param kwargs: Forwarded to the exposure :meth:`~matplotlib.axes.Axes.bar`. 1341 1491 :return: The axes. 1342 1492 """ 1343 1493 ax = ax if ax is not None else _new_axes() ··· 1345 1495 labels = [*result.labels, "A(8)"] 1346 1496 values = [*partials.tolist(), float(result.a8)] 1347 1497 positions = np.arange(len(values), dtype=np.float64) 1348 - colors = ["#bbbbbb"] * partials.size + ["#1f77b4"] 1498 + colors = [_C_MUTED] * partials.size + [_C_PRIMARY] 1349 1499 kwargs.setdefault("color", colors) 1350 1500 ax.bar(positions, values, **kwargs) 1351 1501 ax.set_xticks(positions) ··· 1355 1505 assessment = result.assessment 1356 1506 eav = float(assessment.action_value) 1357 1507 elv = float(assessment.limit_value) 1358 - ax.axhline(eav, color="#ff7f0e", ls="--", label=f"EAV = {eav:g}") 1359 - ax.axhline(elv, color="#d62728", ls="--", label=f"ELV = {elv:g}") 1508 + ax.axhline(eav, color=_C_SECONDARY, ls="--", label=f"EAV = {eav:g}") 1509 + ax.axhline(elv, color=_C_REFERENCE, ls="--", label=f"ELV = {elv:g}") 1360 1510 top = max(elv, float(np.max(values))) * 1.15 1361 1511 ax.set_ylim(0.0, top) 1362 1512 kind = str(assessment.kind).upper() 1363 1513 ax.set_title( 1364 - f"Daily {kind} exposure (A(8) = {float(result.a8):.2f} " 1365 - rf"m/s$^2$, {assessment.zone})" 1514 + f"Directive 2002/44/EC daily {kind} exposure " 1515 + f"(A(8) = {float(result.a8):.2f} " rf"m/s$^2$, {assessment.zone})" 1366 1516 ) 1367 1517 ax.legend(loc="best", fontsize="small") 1368 1518 ax.grid(True, axis="y", alpha=0.3) ··· 1390 1540 freqs = np.asarray(result.frequencies, dtype=np.float64) 1391 1541 levels = np.asarray(result.levels, dtype=np.float64) 1392 1542 for row, idx in zip(NC_CURVES, NC_INDICES): 1393 - ax.plot(OCTAVE_BANDS, row, color="#bbbbbb", lw=0.8, zorder=1) 1543 + ax.plot(OCTAVE_BANDS, row, color=_C_MUTED, lw=0.8, zorder=1) 1394 1544 ax.annotate( 1395 1545 f"{idx:.0f}", (OCTAVE_BANDS[-1], row[-1]), 1396 - fontsize="x-small", color="#888888", va="center", 1546 + fontsize="x-small", color=_C_MUTED, va="center", 1397 1547 ) 1398 1548 valid = ~np.isnan(levels) 1399 - kwargs.setdefault("color", "#1f77b4") 1549 + kwargs.setdefault("color", _C_PRIMARY) 1400 1550 kwargs.setdefault("label", "Measured") 1401 1551 ax.plot(freqs[valid], levels[valid], "o-", zorder=3, **kwargs) 1402 1552 # Nearest *valid* band rather than float equality against the stored ··· 1411 1561 ax.plot( 1412 1562 [freqs[governing]], 1413 1563 [levels[governing]], 1414 - "D", color="#d62728", zorder=4, 1564 + "D", color=_C_REFERENCE, zorder=4, 1415 1565 label=f"Governing band ({_format_freq(result.governing_frequency)})", 1416 1566 ) 1417 1567 _freq_axis(ax, OCTAVE_BANDS) 1418 1568 ax.set_ylabel("Octave-band SPL [dB]") 1419 1569 ax.set_title( 1420 - f"NC-{result.rating:g} " 1570 + f"ANSI/ASA S12.2 NC-{result.rating:g} " 1421 1571 f"({_format_freq(result.governing_frequency)})" 1422 1572 ) 1423 1573 ax.legend(loc=_LEGEND_UPPER_RIGHT, fontsize="small") ··· 1444 1594 reference = np.asarray(result.reference_curve, dtype=np.float64) 1445 1595 valid = ~np.isnan(levels) 1446 1596 1447 - ax.plot(freqs, reference, "s--", color="#7f7f7f", 1597 + ax.plot(freqs, reference, "s--", color=_C_MUTED, 1448 1598 label=f"Reference RC-{result.rating}") 1449 1599 low = freqs <= 500.0 1450 1600 high = freqs >= 1000.0 1451 1601 ax.fill_between(freqs[low], reference[low], reference[low] + 5.0, 1452 - color="#ffbb78", alpha=0.35, label="Rumble tolerance (+5 dB)") 1602 + color=_C_SECONDARY_LIGHT, alpha=0.35, 1603 + label="Rumble tolerance (+5 dB)") 1453 1604 ax.fill_between(freqs[high], reference[high], reference[high] + 3.0, 1454 - color="#aec7e8", alpha=0.45, label="Hiss tolerance (+3 dB)") 1455 - kwargs.setdefault("color", "#1f77b4") 1605 + color=_C_PRIMARY_LIGHT, alpha=0.45, 1606 + label="Hiss tolerance (+3 dB)") 1607 + kwargs.setdefault("color", _C_PRIMARY) 1456 1608 kwargs.setdefault("label", "Measured") 1457 1609 ax.plot(freqs[valid], levels[valid], "o-", zorder=3, **kwargs) 1458 1610 _freq_axis(ax, freqs) 1459 1611 ax.set_ylabel("Octave-band SPL [dB]") 1460 - ax.set_title(result.label) 1612 + ax.set_title(f"ANSI/ASA S12.2 {result.label}") 1461 1613 ax.legend(loc=_LEGEND_UPPER_RIGHT, fontsize="small") 1462 1614 ax.grid(True, which="both", alpha=0.3) 1463 1615 return ax ··· 1483 1635 median = np.asarray(result.median, dtype=np.float64) 1484 1636 su = np.asarray(result.spread_upper, dtype=np.float64) 1485 1637 sl = np.asarray(result.spread_lower, dtype=np.float64) 1486 - z90 = 1.2816 # standard-normal quantile of 0.9 1487 1638 1488 - ax.fill_between(freqs, median - z90 * sl, median + z90 * su, 1489 - color="#aec7e8", alpha=0.5, label="10-90 % fractile band") 1490 - kwargs.setdefault("color", "#1f77b4") 1639 + _fractile_band(ax, freqs, median, sl, su, color=_C_PRIMARY_LIGHT) 1640 + kwargs.setdefault("color", _C_PRIMARY) 1491 1641 ax.plot(freqs, median, "o-", label="Median", **kwargs) 1492 1642 if abs(result.fractile - 0.5) > 1e-9: 1493 1643 ax.plot(freqs, np.asarray(result.threshold, dtype=np.float64), "s--", 1494 - color="#d62728", label=f"Fractile {result.fractile:g}") 1644 + color=_C_REFERENCE, label=f"Fractile {result.fractile:g}") 1495 1645 _freq_axis(ax, freqs) 1496 1646 ax.set_ylabel("Threshold deviation from age 18 [dB]") 1497 1647 ax.invert_yaxis() # audiogram convention: worse hearing downward ··· 1521 1671 median = np.asarray(result.median, dtype=np.float64) 1522 1672 du = np.asarray(result.spread_upper, dtype=np.float64) 1523 1673 dl = np.asarray(result.spread_lower, dtype=np.float64) 1524 - z90 = 1.2816 # standard-normal quantile of 0.9 1525 1674 1526 - ax.fill_between(freqs, np.maximum(median - z90 * dl, 0.0), median + z90 * du, 1527 - color="#ffbb78", alpha=0.5, label="10-90 % fractile band") 1528 - kwargs.setdefault("color", "#ff7f0e") 1675 + _fractile_band(ax, freqs, median, dl, du, color=_C_SECONDARY_LIGHT, floor=0.0) 1676 + kwargs.setdefault("color", _C_SECONDARY) 1529 1677 ax.plot(freqs, median, "o-", label="Median $N_{50}$", **kwargs) 1530 1678 if abs(result.fractile - 0.5) > 1e-9: 1531 1679 ax.plot(freqs, np.asarray(result.value, dtype=np.float64), "s--", 1532 - color="#d62728", label=f"Fractile {result.fractile:g}") 1680 + color=_C_REFERENCE, label=f"Fractile {result.fractile:g}") 1533 1681 _freq_axis(ax, freqs) 1534 1682 ax.set_ylabel("NIPTS [dB]") 1535 1683 ax.invert_yaxis() # audiogram convention: worse hearing downward ··· 1555 1703 ax = ax if ax is not None else _new_axes() 1556 1704 freqs = np.asarray(result.frequencies, dtype=np.float64) 1557 1705 ax.plot(freqs, np.asarray(result.htla, dtype=np.float64), "o-", 1558 - color="#1f77b4", label="Age (HTLA, ISO 7029)") 1706 + color=_C_PRIMARY, label="Age (HTLA, ISO 7029)") 1559 1707 ax.plot(freqs, np.asarray(result.nipts, dtype=np.float64), "^-", 1560 - color="#ff7f0e", label="Noise (NIPTS)") 1561 - kwargs.setdefault("color", "#d62728") 1708 + color=_C_SECONDARY, label="Noise (NIPTS)") 1709 + kwargs.setdefault("color", _C_REFERENCE) 1562 1710 ax.plot(freqs, np.asarray(result.threshold, dtype=np.float64), "s--", 1563 1711 label="Age + noise (HTLAN)", **kwargs) 1564 1712 _freq_axis(ax, freqs) ··· 1595 1743 per_max = float(per.max()) if per.size else 0.0 1596 1744 p_max = max(per_max, result.prominence, 15.0) + 1.0 1597 1745 grid = np.linspace(0.0, p_max, 200) 1598 - ax.plot(grid, impulse_adjustment(grid), color="#1f77b4", 1746 + ax.plot(grid, impulse_adjustment(grid), color=_C_PRIMARY, 1599 1747 label=r"$K_I = 1.8\,(P-5)$") 1600 - ax.axvline(ADJUSTMENT_THRESHOLD, color="#7f7f7f", ls=":", 1748 + ax.axvline(ADJUSTMENT_THRESHOLD, color=_C_MUTED, ls=":", 1601 1749 label=f"threshold $P = {ADJUSTMENT_THRESHOLD:g}$") 1602 1750 1603 - kwargs.setdefault("color", "#aec7e8") 1751 + kwargs.setdefault("color", _C_PRIMARY_LIGHT) 1604 1752 kwargs.setdefault("zorder", 3) 1605 1753 ax.scatter(per, impulse_adjustment(per), label="Impulses", **kwargs) 1606 - ax.scatter([result.prominence], [result.adjustment], color="#d62728", 1754 + ax.scatter([result.prominence], [result.adjustment], color=_C_REFERENCE, 1607 1755 zorder=4, s=90, marker="*", 1608 1756 label=f"Governing P = {result.prominence:.2f}, " 1609 1757 f"$K_I$ = {result.adjustment:.1f} dB") ··· 1629 1777 ax = ax if ax is not None else _new_axes() 1630 1778 freq = np.asarray(result.frequencies, dtype=np.float64) 1631 1779 rt = np.asarray(result.reverberation_time, dtype=np.float64) 1632 - kwargs.setdefault("color", "#1f77b4") 1780 + kwargs.setdefault("color", _C_PRIMARY) 1633 1781 kwargs.setdefault("marker", "o") 1634 - ax.semilogx(freq, rt, **kwargs) 1635 - ax.set_xticks(freq) 1636 - ax.set_xticklabels([f"{f:g}" if f < 1000 else f"{f / 1000:g}k" for f in freq]) 1637 - ax.set_xlabel("Octave-band centre frequency [Hz]") 1782 + ax.plot(freq, rt, **kwargs) 1783 + _freq_axis(ax, freq) 1638 1784 ax.set_ylabel("Reverberation time $T$ [s]") 1639 1785 ax.set_title("EN 12354-6 reverberation time") 1640 1786 ax.set_ylim(bottom=0.0) ··· 1659 1805 r_max = max(result.risk, r90) * 1.3 1660 1806 grid = np.linspace(0.0, r_max, 240) 1661 1807 prob = np.asarray(injury_probability(grid, sex=result.sex), dtype=np.float64) 1662 - ax.plot(grid, 100.0 * prob, color="#1f77b4", 1808 + ax.plot(grid, 100.0 * prob, color=_C_PRIMARY, 1663 1809 label=r"$\Pi(R) = 1 - e^{-(R/\alpha)^{\beta}}$") 1664 1810 for level, r_val in zip((10, 50, 90), (r10, r50, r90)): 1665 - ax.axhline(level, color="#7f7f7f", ls=":", lw=0.8) 1666 - ax.plot([r_val, r_val], [0.0, level], color="#7f7f7f", ls=":", lw=0.8) 1811 + ax.axhline(level, color=_C_MUTED, ls=":", lw=0.8) 1812 + ax.plot([r_val, r_val], [0.0, level], color=_C_MUTED, ls=":", lw=0.8) 1667 1813 1668 - kwargs.setdefault("color", "#d62728") 1814 + kwargs.setdefault("color", _C_REFERENCE) 1669 1815 kwargs.setdefault("zorder", 4) 1670 1816 kwargs.setdefault("s", 90) 1671 1817 ax.scatter([result.risk], [100.0 * result.probability], ··· 1700 1846 contributions = np.asarray(result.contributions, dtype=np.float64) 1701 1847 names = list(result.names) or [f"x{i + 1}" for i in range(contributions.size)] 1702 1848 positions = np.arange(contributions.size) 1703 - kwargs.setdefault("color", "#1f77b4") 1849 + kwargs.setdefault("color", _C_PRIMARY) 1704 1850 ax.barh(positions, contributions, **kwargs) 1705 - ax.axvline(result.combined_uncertainty, color="#d62728", ls="--", 1851 + ax.axvline(result.combined_uncertainty, color=_C_REFERENCE, ls="--", 1706 1852 label=f"$u_c$ = {result.combined_uncertainty:.3g}") 1707 1853 ax.set_yticks(positions) 1708 1854 ax.set_yticklabels(names) 1709 1855 ax.invert_yaxis() 1710 1856 ax.set_xlabel("Contribution to combined uncertainty $|c_i|\\,u(x_i)$") 1711 - ax.set_title(f"Uncertainty budget — y = {result.value:.4g}") 1857 + ax.set_title(f"GUM uncertainty budget — y = {result.value:.4g}") 1712 1858 ax.legend(loc="lower right", fontsize="small") 1713 1859 ax.grid(True, axis="x", alpha=0.3) 1714 1860 return ax 1861 + 1862 + 1863 + def plot_monte_carlo( 1864 + result: "MonteCarloResult", ax: Axes | None = None, **kwargs: Any 1865 + ) -> Axes: 1866 + """Histogram of the Monte Carlo output with the coverage interval marked. 1867 + 1868 + :param result: A :class:`~phonometry.uncertainty.MonteCarloResult` 1869 + obtained with ``keep_samples=True`` (the histogram needs the raw 1870 + output sample). 1871 + :param ax: Existing axes, or ``None`` to create a figure. 1872 + :param kwargs: Forwarded to :meth:`~matplotlib.axes.Axes.hist`. 1873 + :return: The axes. 1874 + :raises ValueError: If the result carries no output samples. 1875 + """ 1876 + if result.samples is None: 1877 + raise ValueError( 1878 + "plot() needs the Monte Carlo output samples; call " 1879 + "monte_carlo(..., keep_samples=True) to retain them." 1880 + ) 1881 + ax = ax if ax is not None else _new_axes() 1882 + samples = np.asarray(result.samples, dtype=np.float64) 1883 + kwargs.setdefault("color", _C_PRIMARY_LIGHT) 1884 + kwargs.setdefault("bins", 120) 1885 + ax.hist(samples, density=True, **kwargs) 1886 + low, high = result.interval 1887 + ax.axvspan(low, high, color=_C_PRIMARY, alpha=0.12, 1888 + label=f"{100.0 * result.coverage:g} % coverage interval") 1889 + ax.axvline(result.value, color=_C_REFERENCE, ls="--", 1890 + label=f"y = {result.value:.4g}") 1891 + ax.set_xlabel("Output quantity y") 1892 + ax.set_ylabel("Probability density") 1893 + ax.set_title( 1894 + "Monte Carlo distribution (GUM Supplement 1) — " 1895 + f"u(y) = {result.standard_uncertainty:.3g}" 1896 + ) 1897 + ax.legend(loc=_LEGEND_UPPER_RIGHT, fontsize="small") 1898 + ax.grid(True, axis="y", alpha=0.3) 1899 + return ax 1900 + 1901 + 1902 + # --------------------------------------------------------------------------- 1903 + # Open-plan offices (ISO 3382-3) 1904 + # --------------------------------------------------------------------------- 1905 + 1906 + 1907 + def plot_open_plan( 1908 + result: "OpenPlanResult", ax: Axes | None = None, **kwargs: Any 1909 + ) -> Axes: 1910 + """Spatial decay of speech with the distraction/privacy distances marked. 1911 + 1912 + Redraws the Clause 6.2 regression line ``Lp,A,S(r) = Lp,A,S,4m - 1913 + D2,S lg(r/4)/lg 2`` over the 2 m to 16 m fitting range (extended to 1914 + reach ``rP`` when it lies further out) and marks the distraction 1915 + distance ``rD`` and the privacy distance ``rP`` from the STI regression. 1916 + 1917 + :param result: An :class:`~phonometry.open_plan.OpenPlanResult`. 1918 + :param ax: Existing axes, or ``None`` to create a figure. 1919 + :param kwargs: Forwarded to the decay-line ``plot`` call. 1920 + :return: The axes. 1921 + :raises ValueError: If the spatial-decay regression is undefined 1922 + (``d2s`` / ``lp_as_4m`` are NaN). 1923 + """ 1924 + import matplotlib.ticker as mticker 1925 + 1926 + if not (np.isfinite(result.d2s) and np.isfinite(result.lp_as_4m)): 1927 + raise ValueError( 1928 + "plot() needs the spatial-decay regression; this result's d2s / " 1929 + "lp_as_4m are NaN (fewer than two positions in the 2 m to 16 m " 1930 + "range)." 1931 + ) 1932 + ax = ax if ax is not None else _new_axes() 1933 + r_max = 16.0 1934 + for marker in (result.rd, result.rp): 1935 + if np.isfinite(marker): 1936 + r_max = max(r_max, 1.15 * marker) 1937 + r = np.geomspace(2.0, r_max, 200) 1938 + level = result.lp_as_4m - result.d2s * np.log2(r / 4.0) 1939 + 1940 + kwargs.setdefault("color", _C_PRIMARY) 1941 + kwargs.setdefault( 1942 + "label", rf"$D_{{2,S}}$ = {result.d2s:.1f} dB per doubling" 1943 + ) 1944 + ax.plot(r, level, **kwargs) 1945 + ax.plot([4.0], [result.lp_as_4m], "o", color=_C_PRIMARY, ms=7, 1946 + label=rf"$L_{{p,A,S,4m}}$ = {result.lp_as_4m:.1f} dB") 1947 + if np.isfinite(result.rd): 1948 + ax.axvline(result.rd, color=_C_SECONDARY, ls="--", 1949 + label=rf"$r_D$ = {result.rd:.1f} m (STI 0.50)") 1950 + if np.isfinite(result.rp): 1951 + ax.axvline(result.rp, color=_C_REFERENCE, ls=":", 1952 + label=rf"$r_P$ = {result.rp:.1f} m (STI 0.20)") 1953 + 1954 + ax.set_xscale("log", base=2) 1955 + ticks = [float(2**k) for k in range(1, int(np.ceil(np.log2(r_max))) + 1)] 1956 + ax.set_xticks(ticks) 1957 + ax.set_xticklabels([f"{t:g}" for t in ticks]) 1958 + ax.xaxis.set_minor_formatter(mticker.NullFormatter()) 1959 + ax.set_xlabel("Distance from the sound source [m]") 1960 + ax.set_ylabel("A-weighted SPL of speech [dB]") 1961 + ax.set_title("ISO 3382-3 spatial decay of speech") 1962 + ax.legend(loc=_LEGEND_UPPER_RIGHT, fontsize="small") 1963 + ax.grid(True, which="both", alpha=0.3) 1964 + return ax 1965 + 1966 + 1967 + # --------------------------------------------------------------------------- 1968 + # Outdoor sound propagation (ISO 9613-2) 1969 + # --------------------------------------------------------------------------- 1970 + 1971 + 1972 + def plot_outdoor_attenuation( 1973 + result: "OutdoorAttenuation", ax: Axes | None = None, **kwargs: Any 1974 + ) -> Axes: 1975 + """Stacked per-band attenuation terms with the total overlaid (ISO 9613-2). 1976 + 1977 + The divergence, atmospheric, ground and barrier terms are stacked per 1978 + octave band on separate positive and negative baselines (the ground 1979 + effect can be a net *gain*), and the total attenuation ``A`` is drawn 1980 + as the primary marker line on top. 1981 + 1982 + :param result: An 1983 + :class:`~phonometry.outdoor_propagation.OutdoorAttenuation`. 1984 + :param ax: Existing axes, or ``None`` to create a figure. 1985 + :param kwargs: Forwarded to the total-attenuation ``plot`` call. 1986 + :return: The axes. 1987 + """ 1988 + ax = ax if ax is not None else _new_axes() 1989 + freqs = np.asarray(result.frequencies, dtype=np.float64) 1990 + positions = _band_axis(ax, freqs) 1991 + n = freqs.size 1992 + 1993 + # Separate positive and negative cumulative baselines so a negative term 1994 + # stacks below zero instead of being drawn on top of the previous bars; 1995 + # the signed heights sum to a_total. 1996 + pos_bottom = np.zeros(n) 1997 + neg_bottom = np.zeros(n) 1998 + terms = ( 1999 + (result.a_div, _C_PRIMARY, "$A_{div}$ — divergence"), 2000 + (result.a_atm, _C_TERTIARY, "$A_{atm}$ — atmospheric"), 2001 + (result.a_gr, _C_QUATERNARY, "$A_{gr}$ — ground"), 2002 + (result.a_bar, _C_SECONDARY, "$A_{bar}$ — barrier"), 2003 + ) 2004 + for values, color, label in terms: 2005 + term = np.asarray(values, dtype=np.float64) 2006 + bottom = np.where(term >= 0.0, pos_bottom, neg_bottom) 2007 + ax.bar(positions, term, bottom=bottom, color=color, label=label) 2008 + pos_bottom += np.maximum(term, 0.0) 2009 + neg_bottom += np.minimum(term, 0.0) 2010 + 2011 + kwargs.setdefault("color", _C_REFERENCE) 2012 + kwargs.setdefault("marker", "D") 2013 + kwargs.setdefault("label", "$A$ — total") 2014 + ax.plot(positions, np.asarray(result.a_total, dtype=np.float64), 2015 + zorder=4, **kwargs) 2016 + ax.axhline(0.0, color=_C_MUTED, lw=0.8) 2017 + ax.set_ylabel("Attenuation A [dB]") 2018 + ax.set_title("ISO 9613-2 attenuation breakdown") 2019 + ax.legend(loc="best", fontsize="small") 2020 + ax.grid(True, axis="y", alpha=0.3) 2021 + return ax 2022 + 2023 + 2024 + # --------------------------------------------------------------------------- 2025 + # Impedance tube (ISO 10534-2) 2026 + # --------------------------------------------------------------------------- 2027 + 2028 + 2029 + def plot_impedance_tube( 2030 + result: "ImpedanceTubeResult", ax: Axes | None = None, **kwargs: Any 2031 + ) -> Axes: 2032 + """Normal-incidence absorption spectrum with |r| overlaid (ISO 10534-2). 2033 + 2034 + Draws the absorption coefficient ``alpha(f)`` as the primary curve and 2035 + the magnitude of the reflection factor ``|r|(f)`` as a muted companion 2036 + (both are dimensionless and share the 0..1 axis). 2037 + 2038 + :param result: An :class:`~phonometry.impedance_tube.ImpedanceTubeResult`. 2039 + :param ax: Existing axes, or ``None`` to create a figure. 2040 + :param kwargs: Forwarded to the absorption-curve ``plot`` call. 2041 + :return: The axes. 2042 + """ 2043 + ax = ax if ax is not None else _new_axes() 2044 + freqs = np.asarray(result.frequency, dtype=np.float64) 2045 + alpha = np.asarray(result.absorption, dtype=np.float64) 2046 + kwargs.setdefault("color", _C_PRIMARY) 2047 + kwargs.setdefault("label", r"Absorption $\alpha$") 2048 + ax.plot(freqs, alpha, **kwargs) 2049 + ax.plot(freqs, np.abs(np.asarray(result.reflection, dtype=np.complex128)), 2050 + ls="--", color=_C_MUTED, label="Reflection factor $|r|$") 2051 + ax.set_xlabel("Frequency [Hz]") 2052 + ax.set_ylabel("Coefficient") 2053 + ax.set_ylim(0.0, 1.05) 2054 + ax.set_title("ISO 10534-2 normal-incidence absorption") 2055 + ax.legend(loc="best", fontsize="small") 2056 + ax.grid(True, alpha=0.3) 2057 + return ax 2058 + 2059 + 2060 + # --------------------------------------------------------------------------- 2061 + # Occupational noise exposure (ISO 9612) 2062 + # --------------------------------------------------------------------------- 2063 + 2064 + 2065 + def plot_occupational_exposure( 2066 + result: "ExposureResult", ax: Axes | None = None, **kwargs: Any 2067 + ) -> Axes: 2068 + """Per-task contributions to the daily exposure level (ISO 9612). 2069 + 2070 + One bar per task (its contribution to ``LEX,8h``), with the combined 2071 + ``LEX,8h`` and the one-sided upper limit ``LEX,8h + U`` as horizontal 2072 + lines. 2073 + 2074 + :param result: An 2075 + :class:`~phonometry.occupational_exposure.ExposureResult` from the 2076 + task-based strategy (the one that carries per-task contributions). 2077 + :param ax: Existing axes, or ``None`` to create a figure. 2078 + :param kwargs: Forwarded to the task :meth:`~matplotlib.axes.Axes.bar`. 2079 + :return: The axes. 2080 + :raises ValueError: If the result carries no per-task contributions. 2081 + """ 2082 + if not result.tasks: 2083 + raise ValueError( 2084 + "plot() needs per-task contributions; only task_based_exposure() " 2085 + "results carry them (the job/full-day strategies do not)." 2086 + ) 2087 + ax = ax if ax is not None else _new_axes() 2088 + contributions = [t.lex_8h_contribution for t in result.tasks] 2089 + labels = [t.label for t in result.tasks] 2090 + positions = np.arange(len(contributions), dtype=np.float64) 2091 + kwargs.setdefault("color", _C_PRIMARY) 2092 + ax.bar(positions, contributions, **kwargs) 2093 + ax.set_xticks(positions) 2094 + ax.set_xticklabels(labels, rotation=45, ha="right") 2095 + 2096 + ax.axhline(result.lex_8h, color=_C_REFERENCE, ls="--", 2097 + label=f"$L_{{EX,8h}}$ = {result.lex_8h:.1f} dB") 2098 + ax.axhline(result.upper_limit, color=_C_MUTED, ls=":", 2099 + label=f"$L_{{EX,8h}} + U$ = {result.upper_limit:.1f} dB") 2100 + ax.set_ylim(0.0, max(result.upper_limit, max(contributions)) * 1.12) 2101 + ax.set_ylabel("A-weighted level [dB]") 2102 + ax.set_title( 2103 + f"ISO 9612 daily noise exposure — $L_{{EX,8h}}$ = " 2104 + f"{result.lex_8h:.1f} dB (U = {result.expanded_uncertainty:.1f} dB)" 2105 + ) 2106 + ax.legend(loc="lower right", fontsize="small") 2107 + ax.grid(True, axis="y", alpha=0.3) 2108 + return ax 2109 + 2110 + 2111 + # --------------------------------------------------------------------------- 2112 + # Static airflow resistance (ISO 9053-1) 2113 + # --------------------------------------------------------------------------- 2114 + 2115 + 2116 + def plot_static_airflow( 2117 + result: "StaticAirflowResult", ax: Axes | None = None, **kwargs: Any 2118 + ) -> Axes: 2119 + """Fitted pressure-drop curve with the evaluation point (ISO 9053-1). 2120 + 2121 + Draws the clause 7.5 through-origin fit ``dp = a u + b u**2`` over twice 2122 + the evaluation range and marks the reference evaluation point. 2123 + 2124 + :param result: A 2125 + :class:`~phonometry.airflow_resistance.StaticAirflowResult`. 2126 + :param ax: Existing axes, or ``None`` to create a figure. 2127 + :param kwargs: Forwarded to the fitted-curve ``plot`` call. 2128 + :return: The axes. 2129 + """ 2130 + ax = ax if ax is not None else _new_axes() 2131 + u_eval = float(result.evaluation_velocity) 2132 + u = np.linspace(0.0, 2.0 * u_eval, 200) 2133 + dp = result.linear_coefficient * u + result.quadratic_coefficient * u**2 2134 + 2135 + kwargs.setdefault("color", _C_PRIMARY) 2136 + kwargs.setdefault("label", r"Fit $\Delta p = a\,u + b\,u^2$") 2137 + # Millimetres per second keep the clause 7.5 reference (0.5 mm/s) legible. 2138 + ax.plot(u * 1e3, dp, **kwargs) 2139 + ax.plot([u_eval * 1e3], [result.pressure_drop], "D", color=_C_REFERENCE, 2140 + ms=7, label=f"Evaluation point (u = {u_eval * 1e3:g} mm/s)") 2141 + ax.set_xlabel("Linear airflow velocity u [mm/s]") 2142 + ax.set_ylabel(r"Pressure difference $\Delta p$ [Pa]") 2143 + ax.set_title( 2144 + "ISO 9053-1 static airflow resistance — " 2145 + f"$R_s$ = {result.specific_resistance:.3g} Pa s/m" 2146 + ) 2147 + ax.legend(loc="upper left", fontsize="small") 2148 + ax.grid(True, alpha=0.3) 2149 + return ax 2150 + 2151 + 2152 + # --------------------------------------------------------------------------- 2153 + # Building performance prediction (EN 12354-1 / EN 12354-2) 2154 + # --------------------------------------------------------------------------- 2155 + 2156 + 2157 + def plot_airborne_prediction( 2158 + result: "AirbornePredictionResult", ax: Axes | None = None, **kwargs: Any 2159 + ) -> Axes: 2160 + """Per-path shares of the transmitted energy (EN 12354-1). 2161 + 2162 + One bar per transmission path (direct plus flanking), sorted by its 2163 + share of the total transmitted sound energy, largest first. 2164 + 2165 + :param result: An 2166 + :class:`~phonometry.building_prediction.AirbornePredictionResult`. 2167 + :param ax: Existing axes, or ``None`` to create a figure. 2168 + :param kwargs: Forwarded to the path :meth:`~matplotlib.axes.Axes.bar`. 2169 + :return: The axes. 2170 + """ 2171 + ax = ax if ax is not None else _new_axes() 2172 + contribs = sorted(result.paths, key=lambda c: c.fraction, reverse=True) 2173 + shares = [100.0 * c.fraction for c in contribs] 2174 + positions = np.arange(len(shares), dtype=np.float64) 2175 + kwargs.setdefault( 2176 + "color", [_C_PRIMARY if c.kind == "Dd" else _C_MUTED for c in contribs] 2177 + ) 2178 + ax.bar(positions, shares, **kwargs) 2179 + ax.set_xticks(positions) 2180 + ax.set_xticklabels([c.label for c in contribs], rotation=45, ha="right") 2181 + ax.set_xlabel("Transmission path") 2182 + ax.set_ylabel("Share of transmitted energy [%]") 2183 + ax.set_title( 2184 + f"EN 12354-1 flanking prediction — R'w = {result.r_prime_w:.1f} dB " 2185 + f"(RDd,w = {result.r_direct_w:.1f} dB)" 2186 + ) 2187 + ax.grid(True, axis="y", alpha=0.3) 2188 + return ax 2189 + 2190 + 2191 + def plot_impact_prediction( 2192 + result: "ImpactPredictionResult", ax: Axes | None = None, **kwargs: Any 2193 + ) -> Axes: 2194 + """Terms of the apparent impact-level prediction (EN 12354-2). 2195 + 2196 + Bars for the Formula 21 terms — the bare-floor equivalent level, the 2197 + covering improvement, the flanking correction — and the resulting 2198 + apparent level ``L'n,w = Ln,w,eq - DLw + K``. 2199 + 2200 + :param result: An 2201 + :class:`~phonometry.building_prediction.ImpactPredictionResult`. 2202 + :param ax: Existing axes, or ``None`` to create a figure. 2203 + :param kwargs: Forwarded to the term :meth:`~matplotlib.axes.Axes.bar`. 2204 + :return: The axes. 2205 + """ 2206 + ax = ax if ax is not None else _new_axes() 2207 + labels = ("$L_{n,w,eq}$", r"$-\Delta L_w$", "$+K$", "$L'_{n,w}$") 2208 + values = ( 2209 + result.ln_w_eq, 2210 + -result.delta_l_w, 2211 + result.k_correction, 2212 + result.l_prime_n_w, 2213 + ) 2214 + positions = np.arange(len(values), dtype=np.float64) 2215 + kwargs.setdefault("color", [_C_MUTED, _C_TERTIARY, _C_SECONDARY, _C_PRIMARY]) 2216 + ax.bar(positions, values, **kwargs) 2217 + ax.axhline(0.0, color=_C_MUTED, lw=0.8) 2218 + ax.set_xticks(positions) 2219 + ax.set_xticklabels(labels) 2220 + ax.set_ylabel("Level / correction [dB]") 2221 + ax.set_title( 2222 + f"EN 12354-2 impact prediction — L'n,w = {result.l_prime_n_w:.1f} dB" 2223 + ) 2224 + ax.grid(True, axis="y", alpha=0.3) 2225 + return ax 2226 + 2227 + 2228 + # --------------------------------------------------------------------------- 2229 + # Field sound insulation spectra (ISO 16283-1 / ISO 16283-2) 2230 + # --------------------------------------------------------------------------- 2231 + 2232 + 2233 + def plot_airborne_insulation( 2234 + result: "AirborneInsulationResult", ax: Axes | None = None, **kwargs: Any 2235 + ) -> Axes: 2236 + """Per-band airborne insulation quantities (ISO 16283-1). 2237 + 2238 + Draws the standardized level difference ``DnT`` first (the primary 2239 + curve), then the level difference ``D`` and, when available, the 2240 + apparent sound reduction index ``R'``. 2241 + 2242 + :param result: An :class:`~phonometry.insulation.AirborneInsulationResult`. 2243 + :param ax: Existing axes, or ``None`` to create a figure. 2244 + :param kwargs: Forwarded to the primary ``DnT`` curve ``plot`` call. 2245 + :return: The axes. 2246 + """ 2247 + curves = [ 2248 + ("$D_{nT}$", np.asarray(result.dnt, dtype=np.float64)), 2249 + ("$D$", np.asarray(result.d, dtype=np.float64)), 2250 + ] 2251 + if result.r_prime is not None: 2252 + curves.append(("$R'$", np.asarray(result.r_prime, dtype=np.float64))) 2253 + return _plot_insulation_bands( 2254 + curves, 2255 + ylabel="Level difference / reduction index [dB]", 2256 + title="Airborne sound insulation (ISO 16283-1)", 2257 + ax=ax, 2258 + **kwargs, 2259 + ) 2260 + 2261 + 2262 + def plot_impact_insulation( 2263 + result: "ImpactInsulationResult", ax: Axes | None = None, **kwargs: Any 2264 + ) -> Axes: 2265 + """Per-band impact sound pressure levels (ISO 16283-2). 2266 + 2267 + Draws the standardized level ``L'nT`` first (the primary curve) and, 2268 + when available, the normalized level ``L'n``. 2269 + 2270 + :param result: An :class:`~phonometry.insulation.ImpactInsulationResult`. 2271 + :param ax: Existing axes, or ``None`` to create a figure. 2272 + :param kwargs: Forwarded to the primary ``L'nT`` curve ``plot`` call. 2273 + :return: The axes. 2274 + """ 2275 + curves = [("$L'_{nT}$", np.asarray(result.l_n_t, dtype=np.float64))] 2276 + if result.l_n is not None: 2277 + curves.append(("$L'_n$", np.asarray(result.l_n, dtype=np.float64))) 2278 + return _plot_insulation_bands( 2279 + curves, 2280 + ylabel="Impact sound pressure level [dB]", 2281 + title="Impact sound insulation (ISO 16283-2)", 2282 + ax=ax, 2283 + **kwargs, 2284 + ) 2285 + 2286 + 2287 + def _plot_insulation_bands( 2288 + curves: "Sequence[tuple[str, np.ndarray]]", 2289 + *, 2290 + ylabel: str, 2291 + title: str, 2292 + ax: Axes | None, 2293 + **kwargs: Any, 2294 + ) -> Axes: 2295 + """Shared per-band insulation renderer (measurement bands are index-only). 2296 + 2297 + The ISO 16283 results do not carry their band centres, so the curves are 2298 + drawn over band indices; user kwargs style the first (primary) curve 2299 + only, mirroring :func:`plot_facade_insulation`. 2300 + """ 2301 + ax = ax if ax is not None else _new_axes() 2302 + n = curves[0][1].size 2303 + x = np.arange(n, dtype=np.float64) 2304 + ax.set_xticks(x) 2305 + ax.set_xticklabels([f"Band {i + 1}" for i in range(n)], 2306 + rotation=45, ha="right") 2307 + ax.set_xlabel("Band") 2308 + for index, (label, y) in enumerate(curves): 2309 + opts: dict[str, Any] = {"label": label} 2310 + if index == 0: 2311 + opts.update(kwargs) 2312 + ax.plot(x, y, "o-", **opts) 2313 + ax.set_ylabel(ylabel) 2314 + ax.set_title(title) 2315 + ax.legend(loc="best", fontsize="small") 2316 + ax.grid(True, alpha=0.3) 2317 + return ax 2318 + 2319 + 2320 + # --------------------------------------------------------------------------- 2321 + # Building-acoustics measurement uncertainty (ISO 12999-1) 2322 + # --------------------------------------------------------------------------- 2323 + 2324 + 2325 + def plot_band_uncertainty( 2326 + result: "BandUncertainty", ax: Axes | None = None, **kwargs: Any 2327 + ) -> Axes: 2328 + """Per-band standard uncertainty of an insulation quantity (ISO 12999-1). 2329 + 2330 + :param result: A 2331 + :class:`~phonometry.building_uncertainty.BandUncertainty`. 2332 + :param ax: Existing axes, or ``None`` to create a figure. 2333 + :param kwargs: Forwarded to the uncertainty curve ``plot`` call. 2334 + :return: The axes. 2335 + """ 2336 + ax = ax if ax is not None else _new_axes() 2337 + freqs, u = result.to_arrays() 2338 + kwargs.setdefault("color", _C_PRIMARY) 2339 + kwargs.setdefault("marker", "o") 2340 + ax.plot(freqs, u, **kwargs) 2341 + _freq_axis(ax, freqs) 2342 + ax.set_ylabel("Standard uncertainty u [dB]") 2343 + ax.set_ylim(bottom=0.0) 2344 + quantity = "sigma_R95 upper limit" if result.upper_limit else "u" 2345 + ax.set_title( 2346 + f"ISO 12999-1 band uncertainty ({quantity}) — " 2347 + f"{result.measurand}, situation {result.situation}" 2348 + ) 2349 + ax.grid(True, which="both", alpha=0.3) 2350 + return ax
+14
src/phonometry/airflow_resistance.py
··· 64 64 import math 65 65 import warnings 66 66 from dataclasses import dataclass 67 + from typing import TYPE_CHECKING, Any 67 68 68 69 import numpy as np 69 70 70 71 from numpy.typing import ArrayLike, NDArray 71 72 72 73 from ._warnings import PhonometryWarning 74 + 75 + if TYPE_CHECKING: 76 + from matplotlib.axes import Axes 73 77 74 78 __all__ = [ 75 79 "AirflowResistanceWarning", ··· 137 141 pressure_drop: float 138 142 linear_coefficient: float 139 143 quadratic_coefficient: float 144 + 145 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 146 + """Plot the fitted ``dp(u)`` curve with the evaluation point. 147 + 148 + Requires matplotlib (``pip install phonometry[plot]``); returns the 149 + :class:`~matplotlib.axes.Axes`. 150 + """ 151 + from ._plotting import plot_static_airflow 152 + 153 + return plot_static_airflow(self, ax=ax, **kwargs) 140 154 141 155 142 156 def linear_airflow_velocity(volume_flow_rate: float, area: float) -> float:
+24 -1
src/phonometry/building_prediction.py
··· 45 45 46 46 from dataclasses import dataclass 47 47 from math import isfinite, log10 48 - from typing import Literal, Sequence 48 + from typing import TYPE_CHECKING, Any, Literal, Sequence 49 + 50 + if TYPE_CHECKING: 51 + from matplotlib.axes import Axes 49 52 50 53 #: Reference coupling length ``l0`` in Formula (28a), in metres (Clause 4.4.1). 51 54 _L0 = 1.0 ··· 145 148 paths: tuple[PathContribution, ...] 146 149 dominant: PathContribution 147 150 151 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 152 + """Plot the per-path shares of the transmitted energy. 153 + 154 + Requires matplotlib (``pip install phonometry[plot]``); returns the 155 + :class:`~matplotlib.axes.Axes`. 156 + """ 157 + from ._plotting import plot_airborne_prediction 158 + 159 + return plot_airborne_prediction(self, ax=ax, **kwargs) 160 + 148 161 149 162 @dataclass(frozen=True) 150 163 class ImpactPredictionResult: ··· 161 174 ln_w_eq: float 162 175 delta_l_w: float 163 176 k_correction: float 177 + 178 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 179 + """Plot the Formula 21 terms and the resulting ``L'n,w``. 180 + 181 + Requires matplotlib (``pip install phonometry[plot]``); returns the 182 + :class:`~matplotlib.axes.Axes`. 183 + """ 184 + from ._plotting import plot_impact_prediction 185 + 186 + return plot_impact_prediction(self, ax=ax, **kwargs) 164 187 165 188 166 189 def _check_finite(value: float, name: str) -> float:
+14 -1
src/phonometry/building_uncertainty.py
··· 46 46 from dataclasses import dataclass 47 47 from math import sqrt 48 48 from types import MappingProxyType 49 - from typing import Dict, Literal, Mapping, Sequence, Tuple 49 + from typing import TYPE_CHECKING, Any, Dict, Literal, Mapping, Sequence, Tuple 50 50 51 51 import numpy as np 52 + 53 + if TYPE_CHECKING: 54 + from matplotlib.axes import Axes 52 55 53 56 Situation = Literal["A", "B", "C"] 54 57 Measurand = Literal["airborne", "impact", "impact_reduction"] ··· 211 214 np.asarray(self.frequencies, dtype=float), 212 215 np.asarray(self.uncertainties, dtype=float), 213 216 ) 217 + 218 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 219 + """Plot the per-band standard uncertainty spectrum. 220 + 221 + Requires matplotlib (``pip install phonometry[plot]``); returns the 222 + :class:`~matplotlib.axes.Axes`. 223 + """ 224 + from ._plotting import plot_band_uncertainty 225 + 226 + return plot_band_uncertainty(self, ax=ax, **kwargs) 214 227 215 228 216 229 @dataclass(frozen=True)
+14
src/phonometry/impedance_tube.py
··· 40 40 41 41 import warnings 42 42 from dataclasses import dataclass 43 + from typing import TYPE_CHECKING, Any 43 44 44 45 import numpy as np 45 46 from numpy.typing import ArrayLike, NDArray 46 47 47 48 from ._types import Real 48 49 from ._warnings import PhonometryWarning 50 + 51 + if TYPE_CHECKING: 52 + from matplotlib.axes import Axes 49 53 50 54 Complex = NDArray[np.complex128] 51 55 ··· 465 469 surface_impedance: Complex 466 470 normalized_impedance: Complex 467 471 absorption: Real 472 + 473 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 474 + """Plot the absorption spectrum ``alpha(f)`` with ``|r|`` overlaid. 475 + 476 + Requires matplotlib (``pip install phonometry[plot]``); returns the 477 + :class:`~matplotlib.axes.Axes`. 478 + """ 479 + from ._plotting import plot_impedance_tube 480 + 481 + return plot_impedance_tube(self, ax=ax, **kwargs) 468 482 469 483 470 484 def two_microphone_impedance(
+20
src/phonometry/insulation.py
··· 191 191 dnt: np.ndarray 192 192 r_prime: np.ndarray | None 193 193 194 + def plot(self, ax: Axes | None = None, **kwargs: Any) -> Axes: 195 + """Plot the per-band insulation quantities (``DnT``, ``D``, ``R'``). 196 + 197 + Requires matplotlib (``pip install phonometry[plot]``); returns the 198 + :class:`~matplotlib.axes.Axes`. 199 + """ 200 + from ._plotting import plot_airborne_insulation 201 + 202 + return plot_airborne_insulation(self, ax=ax, **kwargs) 203 + 194 204 195 205 @dataclass(frozen=True) 196 206 class WeightedRatingResult: ··· 250 260 251 261 l_n_t: np.ndarray 252 262 l_n: np.ndarray | None 263 + 264 + def plot(self, ax: Axes | None = None, **kwargs: Any) -> Axes: 265 + """Plot the per-band impact levels (``L'nT`` and, if present, ``L'n``). 266 + 267 + Requires matplotlib (``pip install phonometry[plot]``); returns the 268 + :class:`~matplotlib.axes.Axes`. 269 + """ 270 + from ._plotting import plot_impact_insulation 271 + 272 + return plot_impact_insulation(self, ax=ax, **kwargs) 253 273 254 274 255 275 @dataclass(frozen=True)
+16 -1
src/phonometry/occupational_exposure.py
··· 44 44 import warnings 45 45 from dataclasses import dataclass, field, replace 46 46 from math import log10, sqrt 47 - from typing import Any, Dict, Literal, Sequence, Tuple 47 + from typing import TYPE_CHECKING, Any, Dict, Literal, Sequence, Tuple 48 48 49 49 import numpy as np 50 50 51 51 from ._levels_math import energy_mean 52 52 from ._warnings import PhonometryWarning, _warn_renamed 53 + 54 + if TYPE_CHECKING: 55 + from matplotlib.axes import Axes 53 56 54 57 #: Reference duration T0 = 8 h (Clause 4). 55 58 _T0: float = 8.0 ··· 277 280 def upper_limit(self) -> float: 278 281 """Upper limit ``LEX,8h + U`` of the one-sided 95 % interval, dB.""" 279 282 return self.lex_8h + self.expanded_uncertainty 283 + 284 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 285 + """Plot the per-task contributions with the ``LEX,8h`` line. 286 + 287 + Only task-based results carry per-task contributions (the job and 288 + full-day strategies raise :class:`ValueError`). Requires matplotlib 289 + (``pip install phonometry[plot]``); returns the 290 + :class:`~matplotlib.axes.Axes`. 291 + """ 292 + from ._plotting import plot_occupational_exposure 293 + 294 + return plot_occupational_exposure(self, ax=ax, **kwargs) 280 295 281 296 282 297 # --------------------------------------------------------------------------- #
+16 -1
src/phonometry/open_plan.py
··· 32 32 from __future__ import annotations 33 33 34 34 from dataclasses import dataclass 35 - from typing import List, Tuple 35 + from typing import TYPE_CHECKING, Any, List, Tuple 36 36 37 37 import numpy as np 38 + 39 + if TYPE_CHECKING: 40 + from matplotlib.axes import Axes 38 41 39 42 #: Reference distance for the logarithmic distance axis of D2,S 40 43 #: (ISO 3382-3:2012, 6.2, Equation (5)): r0 = 1 m. ··· 80 83 lp_as_4m: float 81 84 rd: float 82 85 rp: float 86 + 87 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 88 + """Plot the spatial decay of speech with ``rD``/``rP`` marked. 89 + 90 + Redraws the Clause 6.2 regression line from ``d2s`` and 91 + ``lp_as_4m`` and marks the distraction and privacy distances. 92 + Requires matplotlib (``pip install phonometry[plot]``); returns the 93 + :class:`~matplotlib.axes.Axes`. 94 + """ 95 + from ._plotting import plot_open_plan 96 + 97 + return plot_open_plan(self, ax=ax, **kwargs) 83 98 84 99 85 100 def _linear_fit(x: np.ndarray, y: np.ndarray) -> Tuple[float, float]:
+14
src/phonometry/outdoor_propagation.py
··· 39 39 from __future__ import annotations 40 40 41 41 from dataclasses import dataclass 42 + from typing import TYPE_CHECKING, Any 42 43 43 44 import numpy as np 44 45 from numpy.typing import ArrayLike, NDArray 45 46 46 47 from .air_absorption import air_attenuation 47 48 from ._warnings import _warn_renamed 49 + 50 + if TYPE_CHECKING: 51 + from matplotlib.axes import Axes 48 52 49 53 #: Reference distance ``d0`` in the divergence term (ISO 9613-2:1996, Eq. (7)), m. 50 54 _D0 = 1.0 ··· 139 143 a_bar: NDArray[np.float64] 140 144 a_total: NDArray[np.float64] 141 145 d_omega: NDArray[np.float64] 146 + 147 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 148 + """Plot the stacked per-band attenuation terms with the total. 149 + 150 + Requires matplotlib (``pip install phonometry[plot]``); returns the 151 + :class:`~matplotlib.axes.Axes`. 152 + """ 153 + from ._plotting import plot_outdoor_attenuation 154 + 155 + return plot_outdoor_attenuation(self, ax=ax, **kwargs) 142 156 143 157 144 158 # --------------------------------------------------------------------------- #
+20
src/phonometry/uncertainty.py
··· 144 144 ``(low, high)`` (clause 7.7). 145 145 :ivar coverage: The coverage probability of ``interval``. 146 146 :ivar trials: Number of Monte Carlo trials. 147 + :ivar samples: The raw model-output sample (one value per trial), kept 148 + only when :func:`monte_carlo` is called with ``keep_samples=True`` 149 + (it feeds the output-distribution histogram of :meth:`plot`). 147 150 """ 148 151 149 152 value: float ··· 151 154 interval: tuple[float, float] 152 155 coverage: float 153 156 trials: int 157 + samples: np.ndarray | None = field(default=None, repr=False) 158 + 159 + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": 160 + """Plot the output histogram with the coverage interval marked. 161 + 162 + Needs the raw output sample, so call ``monte_carlo(..., 163 + keep_samples=True)``. Requires matplotlib (``pip install 164 + phonometry[plot]``); returns the :class:`~matplotlib.axes.Axes`. 165 + """ 166 + from ._plotting import plot_monte_carlo 167 + 168 + return plot_monte_carlo(self, ax=ax, **kwargs) 154 169 155 170 156 171 def _sensitivity(model: Model, values: np.ndarray, uncertainties: np.ndarray) -> np.ndarray: ··· 296 311 trials: int = 1_000_000, 297 312 coverage: float = 0.95, 298 313 seed: int | None = None, 314 + keep_samples: bool = False, 299 315 ) -> MonteCarloResult: 300 316 """Propagate uncertainty by the Monte Carlo method (Supplement 1). 301 317 ··· 309 325 :param trials: Number of Monte Carlo trials ``M``. 310 326 :param coverage: Coverage probability of the reported interval. 311 327 :param seed: Optional seed for the random generator (reproducibility). 328 + :param keep_samples: Retain the raw output sample on the result (one 329 + float per trial) so :meth:`MonteCarloResult.plot` can draw the 330 + output-distribution histogram. 312 331 :return: A :class:`MonteCarloResult`. 313 332 :raises ValueError: for no inputs, non-positive trials or bad coverage. 314 333 """ ··· 332 351 interval=(float(low), float(high)), 333 352 coverage=coverage, 334 353 trials=trials, 354 + samples=output if keep_samples else None, 335 355 )
+326 -2
tests/test_result_plots.py
··· 101 101 return ph.sound_intensity(p1, p2, FS, spacing=0.012, fraction=3, limits=[125, 4000]) 102 102 103 103 104 + def _open_plan() -> ph.OpenPlanResult: 105 + positions = np.array([2.0, 3.0, 4.0, 6.0, 8.0, 12.0, 16.0]) 106 + spl = 57.0 - 7.0 * np.log2(positions / 4.0) 107 + sti = np.clip(0.9 - 0.055 * positions, 0.0, 1.0) 108 + return ph.open_plan_metrics(positions, spl, sti) 109 + 110 + 111 + def _outdoor() -> ph.OutdoorAttenuation: 112 + bands = np.array([63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0]) 113 + barrier = ph.Barrier(source_to_edge=101.0, edge_to_receiver=101.0) 114 + return ph.outdoor_propagation_attenuation( 115 + 200.0, 1.5, 1.5, bands, ground_source=1.0, ground_middle=1.0, 116 + ground_receiver=1.0, barrier=barrier, temperature=15.0, 117 + relative_humidity=70.0, 118 + ) 119 + 120 + 121 + def _impedance_tube() -> ph.ImpedanceTubeResult: 122 + f = np.linspace(200.0, 1600.0, 60) 123 + r_true = 0.6 * np.exp(-f / 1200.0) * np.exp(0.8j) 124 + k = 2.0 * np.pi * f / 343.2 125 + s, x1 = 0.05, 0.12 126 + phase = np.exp(2j * k * x1) 127 + h12 = (np.exp(-1j * k * s) + r_true * phase * np.exp(1j * k * s)) / ( 128 + 1.0 + r_true * phase 129 + ) 130 + return ph.two_microphone_impedance( 131 + h12, frequency=f, spacing=s, x1=x1, speed_of_sound=343.2, 132 + characteristic_impedance=407.0, 133 + ) 134 + 135 + 136 + _MC_QUANTITIES = ( 137 + ph.Quantity(74.0, 0.0, name="Reading"), 138 + ph.rectangular(0.0, 0.20, name="Calibration"), 139 + ph.Quantity(0.0, 0.35, dof=9, name="Position"), 140 + ) 141 + 142 + 143 + def _monte_carlo() -> ph.MonteCarloResult: 144 + return ph.monte_carlo( 145 + lambda a, b, c: a + b + c, _MC_QUANTITIES, trials=2000, seed=7, 146 + keep_samples=True, 147 + ) 148 + 149 + 150 + def _exposure() -> ph.ExposureResult: 151 + tasks = [ 152 + ph.Task((86.4, 86.7, 87.0), 2.0, label="grinding"), 153 + ph.Task((80.1, 80.9, 80.5), 3.0, label="welding"), 154 + ph.Task((75.0, 74.6, 74.9), 3.0, label="assembly"), 155 + ] 156 + return ph.task_based_exposure(tasks) 157 + 158 + 159 + def _static_airflow() -> ph.StaticAirflowResult: 160 + u = np.array([0.2e-3, 0.4e-3, 0.6e-3, 0.8e-3, 1.0e-3]) 161 + dp = 30000.0 * u + 4.0e6 * u**2 162 + return ph.static_airflow_resistance(u, dp, area=0.01, thickness=0.05) 163 + 164 + 165 + def _airborne_prediction() -> ph.AirbornePredictionResult: 166 + paths = [] 167 + for name, rw, k_ff, k_side, lf in ( 168 + ("floor", 49.0, 12.4, 8.9, 4.5), 169 + ("facade", 42.0, 12.6, 6.7, 2.55), 170 + ): 171 + ff, df, fd = ph.flanking_element( 172 + label=name, r_flanking=rw, r_separating=57.0, k_ff=k_ff, 173 + k_fd=k_side, k_df=k_side, separating_area=11.5, coupling_length=lf, 174 + ) 175 + paths.extend((ff, df, fd)) 176 + return ph.predicted_airborne_insulation(r_direct=57.0, flanking_paths=paths) 177 + 178 + 179 + def _impact_prediction() -> ph.ImpactPredictionResult: 180 + return ph.predicted_impact_insulation( 181 + ln_w_eq=78.0, delta_l_w=17.0, k_correction=2.0 182 + ) 183 + 184 + 185 + def _airborne_insulation() -> ph.AirborneInsulationResult: 186 + return ph.airborne_insulation( 187 + [70.0, 72.0, 74.0], [40.0, 41.0, 42.0], [0.5, 0.5, 0.5], 188 + area=10.0, volume=50.0, 189 + ) 190 + 191 + 192 + def _impact_insulation() -> ph.ImpactInsulationResult: 193 + return ph.impact_insulation([60.0, 61.0, 62.0], [0.5, 0.5, 0.5], volume=50.0) 194 + 195 + 196 + def _band_uncertainty() -> ph.BandUncertainty: 197 + return ph.band_uncertainty("airborne", "B") 198 + 199 + 104 200 def _intensity_wide() -> ph.IntensityResult: 105 201 """IntensityResult with band centres spanning two decades (100 Hz-10 kHz) 106 202 so the log-axis bar-width scaling can be checked at the extremes.""" ··· 337 433 ("decay_curve", lambda: ph.decay_curve(_exp_ir(seconds=1.0, t60=0.6), FS), "line"), 338 434 ("facade", lambda: ph.facade_insulation( 339 435 [70.0, 72.0, 74.0], [40.0, 41.0, 42.0], [0.5, 0.5, 0.5]), "line"), 436 + ("open_plan", _open_plan, "line"), 437 + ("outdoor", _outdoor, "line"), 438 + ("impedance_tube", _impedance_tube, "line"), 439 + ("monte_carlo", _monte_carlo, "bar"), 440 + ("exposure", _exposure, "bar"), 441 + ("static_airflow", _static_airflow, "line"), 442 + ("airborne_prediction", _airborne_prediction, "bar"), 443 + ("impact_prediction", _impact_prediction, "bar"), 444 + ("airborne_insulation", _airborne_insulation, "line"), 445 + ("impact_insulation", _impact_insulation, "line"), 446 + ("band_uncertainty", _band_uncertainty, "line"), 340 447 ] 341 448 342 449 ··· 436 543 patches = axes[0].patches 437 544 assert len(patches) == 3 * n # EDT/T20/T30 grouped bars 438 545 hatched = [p for p in patches if p.get_hatch()] 439 - greyed = [p for p in patches if p.get_facecolor()[:3] == plt.matplotlib.colors.to_rgb("#bbbbbb")] 546 + greyed = [ 547 + p for p in patches 548 + if p.get_facecolor()[:3] == plt.matplotlib.colors.to_rgb(_plotting._C_MUTED) 549 + ] 440 550 # Exactly one bar per series (EDT/T20/T30) is invalid -> 3 hatched/greyed. 441 551 assert len(hatched) == 3 442 552 assert len(greyed) == 3 ··· 559 669 560 670 561 671 # -------------------------------------------------------------------------- 672 + # Open-plan spatial decay (ISO 3382-3) 673 + # -------------------------------------------------------------------------- 674 + def test_open_plan_plot_line_and_markers() -> None: 675 + res = _open_plan() 676 + ax = res.plot() 677 + # the regression line passes through Lp,A,S,4m at 4 m. 678 + line = ax.lines[0] 679 + x, y = np.asarray(line.get_xdata()), np.asarray(line.get_ydata()) 680 + at4 = float(np.interp(4.0, x, y)) 681 + assert at4 == pytest.approx(res.lp_as_4m, abs=0.05) 682 + # slope over one doubling equals -D2,S. 683 + at8 = float(np.interp(8.0, x, y)) 684 + assert at4 - at8 == pytest.approx(res.d2s, abs=0.05) 685 + # rD / rP are marked as vertical lines at their distances. 686 + vlines = [ 687 + np.asarray(ln.get_xdata())[0] for ln in ax.lines 688 + if np.asarray(ln.get_xdata()).size == 2 689 + and np.asarray(ln.get_xdata())[0] == np.asarray(ln.get_xdata())[1] 690 + ] 691 + assert any(v == pytest.approx(res.rd) for v in vlines) 692 + assert any(v == pytest.approx(res.rp) for v in vlines) 693 + plt.close("all") 694 + 695 + 696 + def test_open_plan_plot_without_regression_raises() -> None: 697 + bare = ph.OpenPlanResult( 698 + d2s=float("nan"), lp_as_4m=float("nan"), rd=float("nan"), rp=float("nan") 699 + ) 700 + with pytest.raises(ValueError, match="regression"): 701 + bare.plot() 702 + 703 + 704 + # -------------------------------------------------------------------------- 705 + # Outdoor attenuation breakdown (ISO 9613-2) 706 + # -------------------------------------------------------------------------- 707 + def test_outdoor_plot_stacks_terms_to_total() -> None: 708 + res = _outdoor() 709 + ax = res.plot() 710 + n = res.frequencies.size 711 + # four stacked terms -> 4 bars per band; signed heights sum to a_total. 712 + assert len(ax.patches) == 4 * n 713 + heights = np.array([p.get_height() for p in ax.patches]).reshape(4, n) 714 + np.testing.assert_allclose(heights.sum(axis=0), res.a_total, atol=1e-9) 715 + # the ground term is a net gain (negative) at 63 Hz in this scenario. 716 + assert res.a_gr[0] < 0.0 717 + # the total line echoes a_total. 718 + np.testing.assert_allclose(ax.lines[0].get_ydata(), res.a_total) 719 + plt.close("all") 720 + 721 + 722 + # -------------------------------------------------------------------------- 723 + # Impedance tube (ISO 10534-2) 724 + # -------------------------------------------------------------------------- 725 + def test_impedance_tube_plot_alpha_and_reflection() -> None: 726 + res = _impedance_tube() 727 + ax = res.plot() 728 + np.testing.assert_allclose(ax.lines[0].get_ydata(), res.absorption) 729 + np.testing.assert_allclose(ax.lines[1].get_ydata(), np.abs(res.reflection)) 730 + assert ax.get_ylim() == (0.0, 1.05) 731 + plt.close("all") 732 + 733 + 734 + # -------------------------------------------------------------------------- 735 + # Monte Carlo output distribution (GUM Supplement 1) 736 + # -------------------------------------------------------------------------- 737 + def test_monte_carlo_plot_histogram_and_interval() -> None: 738 + res = _monte_carlo() 739 + assert res.samples is not None and res.samples.size == res.trials 740 + ax = res.plot() 741 + bars = [ 742 + p for p in ax.patches 743 + if "coverage interval" not in str(p.get_label()) 744 + ] 745 + assert bars, "expected histogram bars" 746 + # the coverage-interval axvspan matches the result's interval. 747 + spans = [ 748 + p for p in ax.patches 749 + if "coverage interval" in str(p.get_label()) 750 + ] 751 + assert spans, "expected the coverage-interval axvspan" 752 + low, high = res.interval 753 + assert spans[0].get_x() == pytest.approx(low) 754 + assert spans[0].get_x() + spans[0].get_width() == pytest.approx(high) 755 + plt.close("all") 756 + 757 + 758 + def test_monte_carlo_plot_without_samples_raises() -> None: 759 + res = ph.monte_carlo( 760 + lambda a, b, c: a + b + c, _MC_QUANTITIES, trials=200, seed=7 761 + ) 762 + assert res.samples is None 763 + with pytest.raises(ValueError, match="keep_samples"): 764 + res.plot() 765 + 766 + 767 + # -------------------------------------------------------------------------- 768 + # Occupational exposure (ISO 9612) 769 + # -------------------------------------------------------------------------- 770 + def test_exposure_plot_task_bars_and_lex_line() -> None: 771 + res = _exposure() 772 + ax = res.plot() 773 + heights = [p.get_height() for p in ax.patches] 774 + np.testing.assert_allclose( 775 + heights, [t.lex_8h_contribution for t in res.tasks] 776 + ) 777 + hlines = [ 778 + np.asarray(ln.get_ydata())[0] for ln in ax.lines 779 + if np.asarray(ln.get_ydata()).size == 2 780 + and np.asarray(ln.get_ydata())[0] == np.asarray(ln.get_ydata())[1] 781 + ] 782 + assert any(v == pytest.approx(res.lex_8h) for v in hlines) 783 + assert any(v == pytest.approx(res.upper_limit) for v in hlines) 784 + plt.close("all") 785 + 786 + 787 + def test_exposure_plot_without_tasks_raises() -> None: 788 + levels = np.full(5, 80.0) 789 + res = ph.job_based_exposure(levels, 6.0) 790 + assert not res.tasks 791 + with pytest.raises(ValueError, match="per-task"): 792 + res.plot() 793 + 794 + 795 + # -------------------------------------------------------------------------- 796 + # Static airflow resistance (ISO 9053-1) 797 + # -------------------------------------------------------------------------- 798 + def test_static_airflow_plot_curve_through_evaluation_point() -> None: 799 + res = _static_airflow() 800 + ax = res.plot() 801 + x, y = ax.lines[0].get_xdata(), ax.lines[0].get_ydata() 802 + # x is in mm/s; the fitted curve passes through the evaluation point. 803 + at_eval = float(np.interp(res.evaluation_velocity * 1e3, x, y)) 804 + assert at_eval == pytest.approx(res.pressure_drop, rel=1e-3) 805 + plt.close("all") 806 + 807 + 808 + # -------------------------------------------------------------------------- 809 + # EN 12354 predictions 810 + # -------------------------------------------------------------------------- 811 + def test_airborne_prediction_plot_sorted_shares() -> None: 812 + res = _airborne_prediction() 813 + ax = res.plot() 814 + heights = [p.get_height() for p in ax.patches] 815 + assert heights == sorted(heights, reverse=True) 816 + assert sum(heights) == pytest.approx(100.0) 817 + assert f"{res.r_prime_w:.1f}" in ax.get_title() 818 + plt.close("all") 819 + 820 + 821 + def test_impact_prediction_plot_terms() -> None: 822 + res = _impact_prediction() 823 + ax = res.plot() 824 + heights = [p.get_height() for p in ax.patches] 825 + np.testing.assert_allclose( 826 + heights, 827 + [res.ln_w_eq, -res.delta_l_w, res.k_correction, res.l_prime_n_w], 828 + ) 829 + plt.close("all") 830 + 831 + 832 + # -------------------------------------------------------------------------- 833 + # ISO 16283 field insulation spectra 834 + # -------------------------------------------------------------------------- 835 + def test_airborne_insulation_plot_curves() -> None: 836 + res = _airborne_insulation() 837 + ax = res.plot() 838 + np.testing.assert_allclose(ax.lines[0].get_ydata(), res.dnt) 839 + np.testing.assert_allclose(ax.lines[1].get_ydata(), res.d) 840 + assert res.r_prime is not None 841 + np.testing.assert_allclose(ax.lines[2].get_ydata(), res.r_prime) 842 + plt.close("all") 843 + 844 + 845 + def test_impact_insulation_plot_curves_and_label_kwarg() -> None: 846 + res = _impact_insulation() 847 + ax = res.plot(label="my measurement") 848 + np.testing.assert_allclose(ax.lines[0].get_ydata(), res.l_n_t) 849 + labels = [str(ln.get_label()) for ln in ax.lines] 850 + # user label styles only the primary curve; companions keep theirs. 851 + assert "my measurement" in labels 852 + assert any("L'_n" in lbl for lbl in labels) 853 + plt.close("all") 854 + 855 + 856 + # -------------------------------------------------------------------------- 857 + # ISO 12999-1 band uncertainty 858 + # -------------------------------------------------------------------------- 859 + def test_band_uncertainty_plot_spectrum() -> None: 860 + res = _band_uncertainty() 861 + ax = res.plot() 862 + freqs, u = res.to_arrays() 863 + np.testing.assert_allclose(ax.lines[0].get_xdata(), freqs) 864 + np.testing.assert_allclose(ax.lines[0].get_ydata(), u) 865 + assert "12999" in ax.get_title() 866 + plt.close("all") 867 + 868 + 869 + # -------------------------------------------------------------------------- 562 870 # Common contract: ax=None creates a figure; passing ax composes 563 871 # -------------------------------------------------------------------------- 564 872 def test_single_axes_plots_accept_external_ax() -> None: 565 - for res in (_zwicker_stationary(), _sti(), _airborne_rating(), _sound_power()): 873 + for res in ( 874 + _zwicker_stationary(), 875 + _sti(), 876 + _airborne_rating(), 877 + _sound_power(), 878 + _open_plan(), 879 + _outdoor(), 880 + _impedance_tube(), 881 + _monte_carlo(), 882 + _exposure(), 883 + _static_airflow(), 884 + _airborne_prediction(), 885 + _impact_prediction(), 886 + _airborne_insulation(), 887 + _impact_insulation(), 888 + _band_uncertainty(), 889 + ): 566 890 fig, ax = plt.subplots() 567 891 out = res.plot(ax=ax) 568 892 assert out is ax