#!/usr/bin/env catnip
# Sensibilité spectrale du bruit d'un profil NACA 0012.
#
# Quatre essais NASA gardent le même profil et le même angle d'attaque, mais
# changent la vitesse de l'écoulement. Une spline lissée reconstruit chaque
# spectre. Sa dérivée par rapport à log2(f) s'exprime directement en dB/octave.
#
# Source : Brooks, Pope & Marcolini, Airfoil Self-Noise (1989), NASA/UCI.
# DOI : 10.24432/C5VW2C — licence CC BY 4.0.
#
# DEPS: matplotlib numpy scipy
numpy = import('numpy')
differentiate = import('scipy.differentiate')
import('scipy.interpolate', 'make_smoothing_spline')
import('scipy.optimize', 'brentq')
mpl = import('matplotlib')
mpl.use('Agg')
plt = import('matplotlib.pyplot')
import('pathlib', 'Path')
script_dir = Path(META.file).parent
data_path = script_dir / 'data' / 'airfoil_self_noise.dat'
output_dir = script_dir / 'output'
output_dir.mkdir(exist_ok=True)
ANGLE = 7.3
CHORD = 0.2286
VELOCITIES = list(31.7, 39.6, 55.5, 71.3)
GRID_SIZE = 480
DERIVATIVE_TOLERANCE = 1.0e-5
PEAK_TOLERANCE_HZ = 0.05
# Colonnes : fréquence, angle, corde, vitesse, épaisseur de déplacement, dB.
data = numpy.loadtxt(str(data_path))
struct SpectrumAnalysis {
velocity: float; displacement: float; frequencies; pressure;
grid; smoothed; numeric_slope; raw_slope; peak_frequency: float; peak_pressure: float;
rmse: float; roughness_ratio: float; derivative_error: float; reported_error: float;
peak_error_hz: float; differentiation_ok: bool;
label(self): str => { f"{self.velocity:.1f} m/s" }
display(self): None => {
print(f" {self.label()} — {len(self.frequencies)} fréquences, δ*={self.displacement * 1000:.2f} mm")
print(
f" pic={self.peak_frequency:.1f} Hz, niveau={self.peak_pressure:.2f} dB, " +
f"écart spline={self.rmse:.3f} dB"
)
print(
f" dérivée : erreur oracle={self.derivative_error:.2e} dB/octave, " +
f"erreur SciPy≤{self.reported_error:.2e}"
)
print(f" pente brute/lissée : rugosité ×{self.roughness_ratio:.2f}")
}
}
select_spectrum = (velocity: float) => {
mask = numpy.logical_and(
numpy.isclose(data[:, 1], ANGLE),
numpy.logical_and(
numpy.isclose(data[:, 2], CHORD),
numpy.isclose(data[:, 3], velocity),
),
)
rows = data[mask]
rows[numpy.argsort(rows[:, 0])]
}
analyze = (velocity: float): SpectrumAnalysis => {
rows = select_spectrum(velocity)
frequencies = rows[:, 0]
pressure = rows[:, 5]
log_frequencies = numpy.log2(frequencies)
# Le lissage GCV évite que la dérivation amplifie chaque fluctuation locale.
spline = make_smoothing_spline(log_frequencies, pressure)
analytic_slope = spline.derivative()
grid = numpy.linspace(float(log_frequencies[0]), float(log_frequencies[-1]), GRID_SIZE)
smoothed = spline(grid)
# scipy.differentiate ne connaît pas la forme de la spline : il l'évalue
# comme une boîte noire. La dérivée de la BSpline sert ensuite d'oracle.
derivative_result = differentiate.derivative(spline, grid, initial_step=0.01)
numeric_slope = derivative_result.df
oracle_slope = analytic_slope(grid)
derivative_error = float(numpy.max(numpy.abs(numeric_slope - oracle_slope)))
reported_error = float(numpy.max(derivative_result.error))
# Les différences sur les mesures brutes montrent l'amplification du bruit.
raw_slope = numpy.gradient(pressure, log_frequencies)
smooth_on_samples = analytic_slope(log_frequencies)
raw_roughness = float(numpy.std(numpy.diff(raw_slope)))
smooth_roughness = float(numpy.std(numpy.diff(smooth_on_samples)))
roughness_ratio = raw_roughness / smooth_roughness
# Le premier passage de pente positive à négative encadre le pic spectral.
bracket = None
i = 0
while i < len(grid) - 1 {
if numeric_slope[i] >= 0.0 and numeric_slope[i + 1] < 0.0 {
bracket = list(float(grid[i]), float(grid[i + 1]))
break
}
i = i + 1
}
if bracket is None {
raise RuntimeError(f"aucun pic spectral encadré pour {velocity} m/s")
}
slope_at = (point) => {
float(differentiate.derivative(spline, point, initial_step=0.01).df)
}
# Robustesse : le bracket vient de la dérivée numérique mais borne aussi la racine
# analytique — valide tant que les deux coïncident (contrôlé par derivative_error).
peak_log = brentq(slope_at, bracket[0], bracket[1], xtol=1.0e-11)
oracle_peak_log = brentq(analytic_slope, bracket[0], bracket[1], xtol=1.0e-11)
peak_frequency = 2.0 ** peak_log
oracle_peak_frequency = 2.0 ** oracle_peak_log
SpectrumAnalysis(
velocity,
float(rows[0, 4]),
frequencies,
pressure,
grid,
smoothed,
numeric_slope,
raw_slope,
peak_frequency,
float(spline(peak_log)),
float(numpy.sqrt(numpy.mean((spline(log_frequencies) - pressure) ** 2))),
roughness_ratio,
derivative_error,
reported_error,
abs(peak_frequency - oracle_peak_frequency),
bool(numpy.all(derivative_result.success)),
)
}
print("⇒ Sensibilité spectrale du bruit d'un profil NACA 0012")
print(f" Source : NASA/UCI, {len(data)} mesures sous CC BY 4.0")
print(f" Configuration : angle={ANGLE}°, corde={CHORD} m")
print()
# Même analyse sur les quatre vitesses : le scénario change, le calcul reste.
analyses = VELOCITIES.[(velocity) => { analyze(velocity) }]
analyses.[(analysis) => { analysis.display() }]
# Deux oracles indépendants : dérivée analytique de la spline et position du pic.
oracle_ok = True
for analysis in analyses {
if len(analysis.frequencies) != 16 { oracle_ok = False }
if not analysis.differentiation_ok { oracle_ok = False }
if analysis.derivative_error > DERIVATIVE_TOLERANCE { oracle_ok = False }
if analysis.peak_error_hz > PEAK_TOLERANCE_HZ { oracle_ok = False }
if analysis.roughness_ratio <= 1.0 { oracle_ok = False }
}
print()
print(f"⇒ Oracle global (4 spectres + dérivée spline + pics) : {oracle_ok}")
# Mesures et splines en haut ; pentes brutes et dérivées numériques en bas.
fig, axes = plt.subplots(2, 1, figsize=tuple(12, 9), sharex=True)
colors = list('#4c78a8', '#f58518', '#54a24b', '#b279a2')
i = 0
while i < len(analyses) {
analysis = analyses[i]
color = colors[i]
smooth_frequency = 2.0 ** analysis.grid
axes[0].scatter(
analysis.frequencies,
analysis.pressure,
color=color,
marker='o',
s=22,
alpha=0.55,
)
axes[0].plot(smooth_frequency, analysis.smoothed, color=color, linewidth=2.0, label=analysis.label())
axes[0].scatter(
list(analysis.peak_frequency),
list(analysis.peak_pressure),
color=color,
edgecolor='#111111',
linewidth=0.7,
s=65,
zorder=4,
)
axes[1].scatter(
analysis.frequencies,
analysis.raw_slope,
color=color,
marker='x',
s=24,
alpha=0.35,
)
axes[1].plot(smooth_frequency, analysis.numeric_slope, color=color, linewidth=2.0, label=analysis.label())
axes[1].scatter(list(analysis.peak_frequency), list(0.0), color=color, s=45, zorder=4)
i = i + 1
}
axes[0].set_title("Spectres mesurés et splines lissées")
axes[0].set_ylabel("niveau acoustique (dB)")
axes[0].grid(True, alpha=0.2)
axes[0].legend(title="vitesse de l'écoulement")
axes[1].axhline(0.0, color='#111111', linewidth=1.0)
axes[1].set_xscale('log', base=2)
axes[1].set_xticks(list(250, 500, 1000, 2000, 4000, 8000))
axes[1].set_xticklabels(list('250', '500', '1 000', '2 000', '4 000', '8 000'))
axes[1].set_xlabel("fréquence (Hz, échelle logarithmique)")
axes[1].set_ylabel("pente spectrale (dB/octave)")
axes[1].set_title("Différences brutes (×) et dérivées des splines")
axes[1].grid(True, alpha=0.2)
fig.suptitle("NASA Airfoil Self-Noise — la pente révèle les pics", fontsize=15)
fig.tight_layout()
output_path = output_dir / 'airfoil_spectral_sensitivity.png'
fig.savefig(str(output_path), dpi=135, bbox_inches='tight')
print(f"⇒ Figure → {output_path}")
if not oracle_ok {
raise RuntimeError("les contrôles de différentiation spectrale ont échoué")
}
# Aperçu navigateur : sert la figure une fois puis rend la main.
# --no-browser garde un chemin headless (le PNG reste écrit ci-dessus).
if '--no-browser' not in import('sys').argv {
http = import('http')
b64 = import('base64').b64encode(output_path.read_bytes()).decode('ascii')
http.serve(f'<!doctype html><meta charset="utf-8"><title>Airfoil spectral sensitivity</title><body style="margin:0;background:#0d1117"><img style="max-width:100%;display:block;margin:0 auto" src="data:image/png;base64,{b64}">',
0, 'text/html; charset=utf-8', True)
}