feat(aircraft): tone correction (encircling method) validated against ETM Table 3-7 · jmrp.io/PyOctaveBand@0b2cbef

+139

2 changed files

Expand all

src

phonometry

aircraft_noise.py

tests

test_aircraft_noise.py

+112

src/phonometry/aircraft_noise.py

··· 140 140 if total <= 0.0: 141 141 return 0.0 142 142 return float(40.0 + _PNL_FACTOR * np.log10(total)) 143 + 144 + 145 + #: First band of the tone-correction routine (band 3 = 80 Hz, 0-based index 2). 146 + _TONE_START = 2 147 + 148 + 149 + def _tone_background( 150 + spl: "NDArray[np.float64] | list[float]", 151 + ) -> "tuple[NDArray[np.float64], NDArray[np.float64]]": 152 + """Tone-correction background level ``SPL''`` and excess ``F`` (App. 2 §4.3). 153 + 154 + Implements Steps 1-8 of the slope ("encircling") method: slopes, the 155 + encircling of irregular slopes, the adjusted spectrum ``SPL'``, its smoothed 156 + slopes, the background ``SPL''`` and the tone excess ``F = SPL − SPL''``. 157 + Bands 1-2 (50, 63 Hz) do not participate; the routine starts at 80 Hz. 158 + 159 + :return: ``(spl_background, excess_f)``, each length 24; entries below the 160 + start band are zero. 161 + """ 162 + sig = _validate_spectrum(spl) 163 + n = 24 164 + start = _TONE_START 165 + 166 + # Step 1: slopes s(i) = SPL(i) − SPL(i−1), starting at band 4 (index 3). 167 + s = np.full(n, np.nan) 168 + for i in range(start + 1, n): 169 + s[i] = sig[i] - sig[i - 1] 170 + 171 + # Step 2: encircle a slope where |s(i) − s(i−1)| > 5. 172 + s_enc = np.zeros(n, dtype=bool) 173 + for i in range(start + 1, n): 174 + if not np.isnan(s[i]) and not np.isnan(s[i - 1]) and abs(s[i] - s[i - 1]) > 5.0: 175 + s_enc[i] = True 176 + 177 + # Step 3: mark which SPL value to adjust. 178 + spl_marked = np.zeros(n, dtype=bool) 179 + for i in range(start + 1, n): 180 + if not s_enc[i]: 181 + continue 182 + if s[i] > 0.0 and (np.isnan(s[i - 1]) or s[i] > s[i - 1]): 183 + spl_marked[i] = True 184 + elif s[i] <= 0.0 and not np.isnan(s[i - 1]) and s[i - 1] > 0.0: 185 + spl_marked[i - 1] = True 186 + 187 + # Step 4: adjusted spectrum SPL'. 188 + spl2 = sig.copy() 189 + for i in range(start, n): 190 + if spl_marked[i]: 191 + if i < n - 1: 192 + spl2[i] = 0.5 * (sig[i - 1] + sig[i + 1]) 193 + else: # last band: SPL'(24) = SPL(23) + s(23) 194 + spl2[i] = sig[i - 1] + s[i - 1] 195 + 196 + # Step 5: new slopes s'(i); s'(start) = s'(start+1); imaginary 25th slope. 197 + sp = np.full(n, np.nan) 198 + for i in range(start + 1, n): 199 + sp[i] = spl2[i] - spl2[i - 1] 200 + sp[start] = sp[start + 1] 201 + s_imaginary = sp[n - 1] 202 + 203 + # Step 6: 3-point running mean s̄(i), i = start … n−2. 204 + sbar = np.full(n, np.nan) 205 + for i in range(start, n - 1): 206 + third = sp[i + 2] if i + 2 < n else s_imaginary 207 + sbar[i] = (sp[i] + sp[i + 1] + third) / 3.0 208 + 209 + # Step 7: background SPL''(start) = SPL(start); accumulate the smoothed slopes. 210 + spl3 = np.full(n, np.nan) 211 + spl3[start] = sig[start] 212 + for i in range(start + 1, n): 213 + spl3[i] = spl3[i - 1] + sbar[i - 1] 214 + 215 + # Step 8: tone excess F = SPL − SPL''. 216 + excess = np.zeros(n) 217 + for i in range(start, n): 218 + excess[i] = sig[i] - spl3[i] 219 + 220 + spl3 = np.where(np.isnan(spl3), 0.0, spl3) 221 + return np.asarray(spl3, dtype=np.float64), np.asarray(excess, dtype=np.float64) 222 + 223 + 224 + def _tone_factor(excess: float, frequency: float) -> float: 225 + """Tone-correction factor ``C`` from the excess ``F`` (App. 2 Table A2-2).""" 226 + if not np.isfinite(excess) or excess < 1.5: 227 + return 0.0 228 + if 500.0 <= frequency <= 5000.0: 229 + if excess < 3.0: 230 + return 2.0 * excess / 3.0 - 1.0 231 + if excess < 20.0: 232 + return excess / 3.0 233 + return 20.0 / 3.0 234 + if excess < 3.0: 235 + return excess / 3.0 - 0.5 236 + if excess < 20.0: 237 + return excess / 6.0 238 + return 10.0 / 3.0 239 + 240 + 241 + def tone_correction(spl: "NDArray[np.float64] | list[float]") -> float: 242 + """Tone-correction factor ``C`` for a spectrum (ICAO Annex 16 App. 2 §4.3). 243 + 244 + The maximum over bands of the tone-correction factor derived from the tone 245 + excess ``F`` above the smoothed background, with the 1.5 dB threshold, the 246 + 500 Hz / 5000 Hz frequency split and the 6⅔ dB cap of Table A2-2. 247 + 248 + :param spl: The 24 one-third-octave-band sound pressure levels, in dB. 249 + :return: The tone-correction factor ``C``, in dB (0 if no tones qualify). 250 + :raises ValueError: If the spectrum is not 24 finite levels. 251 + """ 252 + _, excess = _tone_background(spl) 253 + factors = [_tone_factor(float(excess[i]), float(NOY_BANDS[i])) for i in range(24)] 254 + return float(max(factors))

+27

tests/test_aircraft_noise.py

··· 17 17 18 18 from phonometry.aircraft_noise import ( 19 19 NOY_BANDS, 20 + _tone_background, 20 21 perceived_noise_level, 21 22 perceived_noisiness, 23 + tone_correction, 22 24 ) 23 25 24 26 _B = list(NOY_BANDS) 27 + 28 + # ICAO Doc 9501 ETM Vol. I (2018) Table 3-7 turbofan spectrum (bands 1-2 blank). 29 + _SPL_37 = [ 30 + -999.0, -999.0, 70.0, 62.0, 70.0, 80.0, 82.0, 83.0, 76.0, 80.0, 31 + 80.0, 79.0, 78.0, 80.0, 78.0, 76.0, 79.0, 85.0, 79.0, 78.0, 71.0, 60.0, 54.0, 45.0, 32 + ] 25 33 26 34 27 35 def test_noy_breakpoints() -> None: ··· 67 75 spl[b11] = 40.0 68 76 expected = 40.0 + (10.0 / np.log10(2.0)) * np.log10(1.15) 69 77 assert perceived_noise_level(spl) == pytest.approx(expected, rel=1e-9) 78 + 79 + 80 + def test_tone_correction_etm_table_37() -> None: 81 + # ETM Vol. I Table 3-7: the max tone correction is C = 2.0 dB at 2500 Hz. 82 + assert tone_correction(_SPL_37) == pytest.approx(2.0, abs=1e-6) 83 + 84 + 85 + def test_tone_correction_background_column() -> None: 86 + # SPL'' (Step 7) column of Table 3-7 at 125, 160, 2500, 10000 Hz; F at 2500. 87 + spl_dd, excess = _tone_background(_SPL_37) 88 + assert spl_dd[4] == pytest.approx(71.0, abs=0.05) # 125 Hz 89 + assert spl_dd[5] == pytest.approx(77.0 + 2.0 / 3.0, abs=0.05) # 160 Hz 90 + assert spl_dd[17] == pytest.approx(79.0, abs=0.05) # 2500 Hz 91 + assert spl_dd[23] == pytest.approx(45.0, abs=0.05) # 10 kHz 92 + assert excess[17] == pytest.approx(6.0, abs=0.05) # F = 85 − 79 93 + 94 + 95 + def test_tone_correction_none_when_flat() -> None: 96 + assert tone_correction(np.full(24, 60.0)) == 0.0

Configure Feed

Configure Feed