diff --git a/Docs/Aeroacoustics.md b/Docs/Aeroacoustics.md
new file mode 100644
index 00000000000..a10364dce6e
--- /dev/null
+++ b/Docs/Aeroacoustics.md
@@ -0,0 +1,87 @@
+# Aeroacoustic analysis of unsteady pressure
+
+`aeroacoustics.py` converts pressure histories from an unsteady SU2
+simulation into frequency-domain acoustic quantities. It is intended for
+pressure probes from LES, DES, DDES, ZDES, or EDDES cases, but it can analyze
+any uniformly sampled pressure history.
+
+The tool estimates the one-sided power spectral density (PSD) with Welch's
+method and reports:
+
+- pressure PSD in Pa²/Hz;
+- PSD level in dB/Hz;
+- bin-integrated narrowband sound pressure level (SPL) in dB;
+- overall sound pressure level (OASPL) in dB; and
+- the dominant non-zero frequency.
+
+Only NumPy is required.
+
+## Recording pressure probes
+
+SU2 custom outputs can record pressure at one or more points. For example:
+
+```ini
+CUSTOM_OUTPUTS= 'mic1 : Probe{PRESSURE}[1.0, 0.0, 0.0]; \
+                 mic2 : Probe{PRESSURE}[0.0, 1.0, 0.0]'
+HISTORY_OUTPUT= TIME_ITER, CUR_TIME, mic1, mic2
+```
+
+Use a fixed physical time step and write every desired acoustic sample. The
+history column containing physical time is commonly named `Cur_Time`; inspect
+the CSV header and pass its exact name with `--time`.
+
+## Usage
+
+Analyze two pressure columns that are already dimensional:
+
+```bash
+python SU2_PY/aeroacoustics.py history.csv \
+  --pressure mic1 \
+  --pressure mic2 \
+  --skip-samples 10000 \
+  --segment-length 2048 \
+  --summary acoustic_summary.json \
+  --output acoustic_spectrum.csv
+```
+
+If the file contains an iteration rather than physical time, provide the
+sampling rate explicitly:
+
+```bash
+python SU2_PY/aeroacoustics.py history.csv \
+  --sample-rate 50000 \
+  --pressure mic1
+```
+
+Pressure must be dimensional in pascals for physical SPL. If an SU2 case uses
+nondimensional pressure, apply its dimensionalization factor in pascals per
+input unit:
+
+```bash
+python SU2_PY/aeroacoustics.py history.csv \
+  --time Cur_Time \
+  --pressure mic1 \
+  --pressure-scale 101325
+```
+
+The default acoustic reference is 20 µPa, appropriate for airborne sound. Use
+`--reference-pressure` for another medium or convention.
+
+## Sampling guidance
+
+- Resolve the highest frequency of interest with adequate samples per period;
+  the absolute Nyquist limit is half the sample rate.
+- Start collecting samples after initial transients have left the domain, or
+  remove them with `--skip-samples`.
+- Increase `--segment-length` for finer frequency resolution. Frequency-bin
+  width is the sample rate divided by segment length.
+- Longer signals give more Welch segments and a smoother PSD estimate.
+- The default Hann window, 50% overlap, and mean removal are suitable general
+  choices. `--detrend linear` can remove slow numerical drift.
+- Avoid interpreting probe pressure as a far-field acoustic result when the
+  probe is inside a vortical or hydrodynamic near field. Permeable-surface
+  Ffowcs Williams–Hawkings propagation is a larger, separate capability.
+
+The spectrum CSV contains one frequency column followed by PSD, PSD-level, and
+narrowband-SPL columns for every requested pressure signal. The optional JSON
+summary contains scalar metrics suitable for automated comparisons.
diff --git a/SU2_PY/aeroacoustics.py b/SU2_PY/aeroacoustics.py
new file mode 100755
index 00000000000..e7b8be01885
--- /dev/null
+++ b/SU2_PY/aeroacoustics.py
@@ -0,0 +1,420 @@
+#!/usr/bin/env python3
+
+## \file aeroacoustics.py
+#  \brief Spectral analysis of unsteady pressure histories.
+#  \version 8.5.0 "Harrier"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+#
+# SU2 is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# SU2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+"""Convert SU2 pressure histories into acoustic spectra and SPL metrics."""
+
+import argparse
+import csv
+import json
+import math
+import re
+from pathlib import Path
+
+import numpy as np
+
+
+WINDOWS = {
+    "hann": np.hanning,
+    "hamming": np.hamming,
+    "blackman": np.blackman,
+    "rectangular": np.ones,
+}
+
+
+def read_history(filename):
+    """Read a numeric SU2 CSV history file and return its columns."""
+    with open(filename, newline="", encoding="utf-8-sig") as history_file:
+        rows = (
+            row
+            for row in history_file
+            if row.strip() and not row.lstrip().startswith(("%", "#"))
+        )
+        reader = csv.reader(rows)
+        try:
+            headers = [field.strip().strip('"') for field in next(reader)]
+        except StopIteration as error:
+            raise ValueError("history file is empty") from error
+
+        if not headers or any(not field for field in headers):
+            raise ValueError("history file contains an empty column name")
+        if len(set(headers)) != len(headers):
+            raise ValueError("history file contains duplicate column names")
+
+        columns = {header: [] for header in headers}
+        for line_number, row in enumerate(reader, start=2):
+            if len(row) != len(headers):
+                raise ValueError(
+                    "line {} has {} values; expected {}".format(
+                        line_number, len(row), len(headers)
+                    )
+                )
+            try:
+                for header, value in zip(headers, row):
+                    columns[header].append(float(value))
+            except ValueError as error:
+                raise ValueError(
+                    "line {} contains a non-numeric value".format(line_number)
+                ) from error
+
+    if not next(iter(columns.values())):
+        raise ValueError("history file has a header but no data")
+    return {name: np.asarray(values, dtype=float) for name, values in columns.items()}
+
+
+def resolve_column(columns, requested):
+    """Resolve a column exactly, or uniquely without case sensitivity."""
+    if requested in columns:
+        return requested
+    matches = [name for name in columns if name.casefold() == requested.casefold()]
+    if len(matches) == 1:
+        return matches[0]
+    raise ValueError(
+        "column {!r} was not found; available columns: {}".format(
+            requested, ", ".join(columns)
+        )
+    )
+
+
+def sampling_rate(time, relative_tolerance=1.0e-5):
+    """Return the sample rate after verifying monotonically uniform time data."""
+    if time.size < 2:
+        raise ValueError("at least two time samples are required")
+    if not np.all(np.isfinite(time)):
+        raise ValueError("time values must be finite")
+
+    intervals = np.diff(time)
+    if np.any(intervals <= 0.0):
+        raise ValueError("time values must be strictly increasing")
+    mean_interval = float(np.mean(intervals))
+    maximum_error = float(np.max(np.abs(intervals - mean_interval)))
+    if maximum_error > relative_tolerance * mean_interval:
+        message = (
+            "time samples are not uniform " "(maximum relative interval error {:.3g})"
+        )
+        raise ValueError(message.format(maximum_error / mean_interval))
+    return 1.0 / mean_interval
+
+
+def _detrend(signal, mode):
+    if mode == "none":
+        return signal.copy()
+    if mode == "mean":
+        return signal - np.mean(signal)
+    positions = np.arange(signal.size, dtype=float)
+    slope, intercept = np.polyfit(positions, signal, 1)
+    return signal - (slope * positions + intercept)
+
+
+def welch_psd(
+    pressure,
+    sample_rate,
+    segment_length=None,
+    overlap=0.5,
+    window="hann",
+    detrend="mean",
+):
+    """Estimate a one-sided pressure PSD using Welch's averaged periodogram."""
+    pressure = np.asarray(pressure, dtype=float)
+    if pressure.ndim != 1:
+        raise ValueError("pressure data must be one-dimensional")
+    if pressure.size < 2:
+        raise ValueError("at least two pressure samples are required")
+    if not np.all(np.isfinite(pressure)):
+        raise ValueError("pressure values must be finite")
+    if not math.isfinite(sample_rate) or sample_rate <= 0.0:
+        raise ValueError("sample rate must be positive and finite")
+    if not 0.0 <= overlap < 1.0:
+        raise ValueError("overlap must be in the range [0, 1)")
+    if window not in WINDOWS:
+        raise ValueError("unsupported window {!r}".format(window))
+    if detrend not in ("none", "mean", "linear"):
+        raise ValueError("unsupported detrending mode {!r}".format(detrend))
+
+    if segment_length is None:
+        segment_length = min(1024, pressure.size)
+    if segment_length < 2 or segment_length > pressure.size:
+        raise ValueError("segment length must be between 2 and the number of samples")
+
+    step = max(1, int(round(segment_length * (1.0 - overlap))))
+    starts = range(0, pressure.size - segment_length + 1, step)
+    weights = WINDOWS[window](segment_length)
+    window_energy = float(np.sum(weights * weights))
+    if window_energy <= 0.0:
+        raise ValueError(
+            "segment length {} is too short for the {} window".format(
+                segment_length, window
+            )
+        )
+    spectrum = np.zeros(segment_length // 2 + 1, dtype=float)
+    segment_count = 0
+
+    for start in starts:
+        segment = _detrend(pressure[start : start + segment_length], detrend)
+        transform = np.fft.rfft(segment * weights)
+        periodogram = np.abs(transform) ** 2 / (sample_rate * window_energy)
+        if segment_length % 2 == 0:
+            periodogram[1:-1] *= 2.0
+        else:
+            periodogram[1:] *= 2.0
+        spectrum += periodogram
+        segment_count += 1
+
+    spectrum /= segment_count
+    frequencies = np.fft.rfftfreq(segment_length, d=1.0 / sample_rate)
+    return frequencies, spectrum, segment_count
+
+
+def acoustic_metrics(frequencies, psd, reference_pressure=20.0e-6):
+    """Calculate PSD level, bin SPL, OASPL, and the dominant frequency."""
+    if reference_pressure <= 0.0 or not math.isfinite(reference_pressure):
+        raise ValueError("reference pressure must be positive and finite")
+    if frequencies.size < 2:
+        raise ValueError("an acoustic spectrum needs at least two frequency bins")
+
+    frequency_step = float(frequencies[1] - frequencies[0])
+    tiny = np.finfo(float).tiny
+    psd_level = 10.0 * np.log10(np.maximum(psd, tiny) / reference_pressure**2)
+    narrowband_spl = 10.0 * np.log10(
+        np.maximum(psd * frequency_step, tiny) / reference_pressure**2
+    )
+    oaspl = 10.0 * math.log10(
+        max(float(np.sum(psd) * frequency_step), tiny) / reference_pressure**2
+    )
+    dominant_index = 1 + int(np.argmax(psd[1:]))
+    return {
+        "psd_level_db_per_hz": psd_level,
+        "narrowband_spl_db": narrowband_spl,
+        "oaspl_db": oaspl,
+        "dominant_frequency_hz": float(frequencies[dominant_index]),
+    }
+
+
+def _safe_column_name(name):
+    normalized = re.sub(r"[^A-Za-z0-9_.-]+", "_", name).strip("_")
+    return normalized or "pressure"
+
+
+def write_spectrum(filename, frequencies, spectra):
+    """Write frequency-domain results for one or more pressure signals."""
+    fields = ["Frequency_Hz"]
+    used_prefixes = set()
+    for name in spectra:
+        base = _safe_column_name(name)
+        prefix = base
+        suffix = 2
+        while prefix in used_prefixes:
+            prefix = "{}_{}".format(base, suffix)
+            suffix += 1
+        used_prefixes.add(prefix)
+        fields.extend(
+            [
+                prefix + "_PSD_Pa2_per_Hz",
+                prefix + "_PSD_level_dB_per_Hz",
+                prefix + "_SPL_dB",
+            ]
+        )
+
+    with open(filename, "w", newline="", encoding="utf-8") as spectrum_file:
+        writer = csv.writer(spectrum_file)
+        writer.writerow(fields)
+        for index, frequency in enumerate(frequencies):
+            row = [frequency]
+            for result in spectra.values():
+                row.extend(
+                    [
+                        result["psd"][index],
+                        result["psd_level_db_per_hz"][index],
+                        result["narrowband_spl_db"][index],
+                    ]
+                )
+            writer.writerow(row)
+
+
+def analyze(args):
+    """Run the command-line analysis and return its summary."""
+    columns = read_history(args.input)
+    pressure_names = [resolve_column(columns, name) for name in args.pressure]
+    if len(set(pressure_names)) != len(pressure_names):
+        raise ValueError("each pressure column may only be requested once")
+    if not math.isfinite(args.pressure_scale) or args.pressure_scale <= 0.0:
+        raise ValueError("pressure scale must be positive and finite")
+    if args.skip_samples < 0:
+        raise ValueError("skip samples must not be negative")
+    sample_count = len(next(iter(columns.values())))
+    if args.skip_samples > sample_count - 2:
+        raise ValueError("skip samples must leave at least two samples to analyze")
+
+    if args.sample_rate is not None:
+        sample_rate_hz = args.sample_rate
+    else:
+        time_name = resolve_column(columns, args.time)
+        sample_rate_hz = sampling_rate(columns[time_name][args.skip_samples :])
+
+    spectra = {}
+    summary = {
+        "input": str(args.input),
+        "sample_rate_hz": sample_rate_hz,
+        "reference_pressure_pa": args.reference_pressure,
+        "samples": sample_count - args.skip_samples,
+        "signals": {},
+    }
+    frequencies = None
+    for name in pressure_names:
+        dimensional_pressure = columns[name][args.skip_samples :] * args.pressure_scale
+        frequencies, psd, segment_count = welch_psd(
+            dimensional_pressure,
+            sample_rate_hz,
+            segment_length=args.segment_length,
+            overlap=args.overlap,
+            window=args.window,
+            detrend=args.detrend,
+        )
+        metrics = acoustic_metrics(frequencies, psd, args.reference_pressure)
+        spectra[name] = {"psd": psd, **metrics}
+        summary["signals"][name] = {
+            "oaspl_db": metrics["oaspl_db"],
+            "dominant_frequency_hz": metrics["dominant_frequency_hz"],
+            "segments": segment_count,
+        }
+
+    summary["frequency_resolution_hz"] = float(frequencies[1] - frequencies[0])
+    summary["nyquist_frequency_hz"] = float(frequencies[-1])
+    write_spectrum(args.output, frequencies, spectra)
+    if args.summary:
+        with open(args.summary, "w", encoding="utf-8") as summary_file:
+            json.dump(summary, summary_file, indent=2)
+            summary_file.write("\n")
+    return summary
+
+
+def build_parser():
+    parser = argparse.ArgumentParser(
+        description=(
+            "Estimate acoustic spectra from pressure columns in an SU2 history CSV."
+        )
+    )
+    parser.add_argument("input", type=Path, help="SU2 history CSV file")
+    parser.add_argument(
+        "-p",
+        "--pressure",
+        action="append",
+        required=True,
+        metavar="COLUMN",
+        help="pressure column to analyze; repeat for multiple probes",
+    )
+    sampling = parser.add_mutually_exclusive_group()
+    sampling.add_argument(
+        "-t",
+        "--time",
+        default="Cur_Time",
+        metavar="COLUMN",
+        help="physical-time column (default: Cur_Time)",
+    )
+    sampling.add_argument(
+        "--sample-rate",
+        type=float,
+        metavar="HZ",
+        help="sampling rate when the file has no physical-time column",
+    )
+    parser.add_argument(
+        "--pressure-scale",
+        type=float,
+        default=1.0,
+        metavar="PA_PER_UNIT",
+        help="multiply input pressures by this dimensionalization factor",
+    )
+    parser.add_argument(
+        "--reference-pressure",
+        type=float,
+        default=20.0e-6,
+        metavar="PA",
+        help="SPL reference pressure (default: 20e-6 Pa)",
+    )
+    parser.add_argument(
+        "-n",
+        "--segment-length",
+        type=int,
+        metavar="SAMPLES",
+        help="Welch segment length (default: min(1024, number of samples))",
+    )
+    parser.add_argument(
+        "--overlap",
+        type=float,
+        default=0.5,
+        metavar="FRACTION",
+        help="Welch segment overlap as a fraction (default: 0.5)",
+    )
+    parser.add_argument(
+        "--window",
+        choices=tuple(WINDOWS),
+        default="hann",
+        help="spectral window (default: hann)",
+    )
+    parser.add_argument(
+        "--detrend",
+        choices=("none", "mean", "linear"),
+        default="mean",
+        help="detrending applied to each segment (default: mean)",
+    )
+    parser.add_argument(
+        "--skip-samples",
+        type=int,
+        default=0,
+        metavar="COUNT",
+        help="discard initial transient samples before analysis (default: 0)",
+    )
+    parser.add_argument(
+        "-o",
+        "--output",
+        type=Path,
+        default=Path("acoustic_spectrum.csv"),
+        help="output spectrum CSV (default: acoustic_spectrum.csv)",
+    )
+    parser.add_argument(
+        "--summary", type=Path, help="optional JSON file for scalar acoustic metrics"
+    )
+    return parser
+
+
+def main():
+    args = build_parser().parse_args()
+    try:
+        summary = analyze(args)
+    except (OSError, ValueError) as error:
+        raise SystemExit("aeroacoustics: error: {}".format(error)) from error
+
+    print("Sample rate: {:.8g} Hz".format(summary["sample_rate_hz"]))
+    for name, result in summary["signals"].items():
+        print(
+            "{}: OASPL = {:.3f} dB, dominant frequency = {:.8g} Hz".format(
+                name, result["oaspl_db"], result["dominant_frequency_hz"]
+            )
+        )
+    print("Spectrum written to {}".format(args.output))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/SU2_PY/meson.build b/SU2_PY/meson.build
index ced81290346..469d0c16486 100644
--- a/SU2_PY/meson.build
+++ b/SU2_PY/meson.build
@@ -1,4 +1,5 @@
 install_data(['continuous_adjoint.py',
+             'aeroacoustics.py',
              'compute_uncertainty.py',
 	     'finite_differences.py',
 	     'mesh_deformation.py',
diff --git a/SU2_PY/tests/test_aeroacoustics.py b/SU2_PY/tests/test_aeroacoustics.py
new file mode 100644
index 00000000000..f05c713d221
--- /dev/null
+++ b/SU2_PY/tests/test_aeroacoustics.py
@@ -0,0 +1,136 @@
+#!/usr/bin/env python3
+
+import csv
+import importlib.util
+import math
+import tempfile
+import unittest
+from pathlib import Path
+from types import SimpleNamespace
+
+import numpy as np
+
+
+MODULE_PATH = Path(__file__).parents[1] / "aeroacoustics.py"
+SPEC = importlib.util.spec_from_file_location("aeroacoustics", MODULE_PATH)
+AEROACOUSTICS = importlib.util.module_from_spec(SPEC)
+SPEC.loader.exec_module(AEROACOUSTICS)
+
+
+class AeroacousticsTests(unittest.TestCase):
+    def test_welch_finds_sine_frequency_and_level(self):
+        sample_rate = 4096.0
+        frequency = 256.0
+        rms_pressure = 1.0
+        time = np.arange(4096) / sample_rate
+        pressure = (
+            math.sqrt(2.0) * rms_pressure * np.sin(2.0 * np.pi * frequency * time)
+        )
+
+        frequencies, psd, segments = AEROACOUSTICS.welch_psd(
+            pressure,
+            sample_rate,
+            segment_length=1024,
+            overlap=0.5,
+            window="hann",
+        )
+        metrics = AEROACOUSTICS.acoustic_metrics(frequencies, psd)
+
+        self.assertEqual(segments, 7)
+        self.assertEqual(metrics["dominant_frequency_hz"], frequency)
+        expected_oaspl = 20.0 * math.log10(rms_pressure / 20.0e-6)
+        self.assertAlmostEqual(metrics["oaspl_db"], expected_oaspl, places=3)
+
+    def test_sampling_rate_rejects_nonuniform_time(self):
+        with self.assertRaisesRegex(ValueError, "not uniform"):
+            AEROACOUSTICS.sampling_rate(np.array([0.0, 0.1, 0.21, 0.3]))
+
+    def test_read_history_accepts_su2_quoted_headers_and_comments(self):
+        with tempfile.TemporaryDirectory() as directory:
+            history = Path(directory) / "history.csv"
+            history.write_text(
+                '% SU2 history\n"Time","probe 1"\n0.0,101325.0\n0.1,101326.0\n',
+                encoding="utf-8",
+            )
+            columns = AEROACOUSTICS.read_history(history)
+
+        self.assertEqual(list(columns), ["Time", "probe 1"])
+        np.testing.assert_allclose(columns["probe 1"], [101325.0, 101326.0])
+
+    def test_write_spectrum_supports_multiple_signals(self):
+        frequencies = np.array([0.0, 1.0])
+        result = {
+            "psd": np.array([1.0, 2.0]),
+            "psd_level_db_per_hz": np.array([3.0, 4.0]),
+            "narrowband_spl_db": np.array([5.0, 6.0]),
+        }
+        with tempfile.TemporaryDirectory() as directory:
+            output = Path(directory) / "spectrum.csv"
+            AEROACOUSTICS.write_spectrum(
+                output, frequencies, {"probe 1": result, "probe/2": result}
+            )
+            with output.open(newline="", encoding="utf-8") as output_file:
+                rows = list(csv.reader(output_file))
+
+        self.assertEqual(len(rows), 3)
+        self.assertIn("probe_1_PSD_Pa2_per_Hz", rows[0])
+        self.assertIn("probe_2_SPL_dB", rows[0])
+
+    def test_write_spectrum_disambiguates_normalized_names(self):
+        frequencies = np.array([0.0, 1.0])
+        result = {
+            "psd": np.ones(2),
+            "psd_level_db_per_hz": np.ones(2),
+            "narrowband_spl_db": np.ones(2),
+        }
+        with tempfile.TemporaryDirectory() as directory:
+            output = Path(directory) / "spectrum.csv"
+            AEROACOUSTICS.write_spectrum(
+                output, frequencies, {"probe 1": result, "probe/1": result}
+            )
+            header = output.read_text(encoding="utf-8").splitlines()[0]
+
+        self.assertIn("probe_1_PSD_Pa2_per_Hz", header)
+        self.assertIn("probe_1_2_PSD_Pa2_per_Hz", header)
+
+    def test_analyze_writes_spectrum_and_summary(self):
+        with tempfile.TemporaryDirectory() as directory:
+            directory = Path(directory)
+            history = directory / "history.csv"
+            output = directory / "spectrum.csv"
+            summary_file = directory / "summary.json"
+            time = np.arange(256) / 128.0
+            pressure = np.sin(2.0 * np.pi * 16.0 * time)
+            with history.open("w", newline="", encoding="utf-8") as history_file:
+                writer = csv.writer(history_file)
+                writer.writerow(["Cur_Time", "mic"])
+                writer.writerows(zip(time, pressure))
+
+            summary = AEROACOUSTICS.analyze(
+                SimpleNamespace(
+                    input=history,
+                    pressure=["mic"],
+                    sample_rate=None,
+                    time="Cur_Time",
+                    pressure_scale=1.0,
+                    reference_pressure=20.0e-6,
+                    skip_samples=64,
+                    segment_length=64,
+                    overlap=0.5,
+                    window="hann",
+                    detrend="mean",
+                    output=output,
+                    summary=summary_file,
+                )
+            )
+
+            self.assertTrue(output.is_file())
+            self.assertTrue(summary_file.is_file())
+            self.assertEqual(summary["samples"], 192)
+            self.assertEqual(summary["frequency_resolution_hz"], 2.0)
+            self.assertEqual(summary["nyquist_frequency_hz"], 64.0)
+            self.assertEqual(summary["signals"]["mic"]["dominant_frequency_hz"], 16.0)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 7f10753db6e..73ae6d5d668 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1041,6 +1041,19 @@ def main():
     square_cylinder.unsteady  = True
     test_list.append(square_cylinder)
 
+    # Aeroacoustics postprocessing of a converged square-cylinder pressure history
+    aeroacoustics_square_cylinder = TestCase('aeroacoustics_square_cylinder')
+    aeroacoustics_square_cylinder.cfg_dir = "unsteady/square_cylinder"
+    aeroacoustics_square_cylinder.cfg_file = "history_aeroacoustics.csv"
+    aeroacoustics_square_cylinder.test_iter = 0
+    aeroacoustics_square_cylinder.test_vals = [103.1512634765, 1.9531250000]
+    aeroacoustics_square_cylinder.command = TestCase.Command(
+        exec = "python", param = "run_aeroacoustics.py"
+    )
+    aeroacoustics_square_cylinder.timeout = 60
+    aeroacoustics_square_cylinder.tol = 1e-6
+    test_list.append(aeroacoustics_square_cylinder)
+
     # Gust
     sine_gust           = TestCase('sine_gust')
     sine_gust.cfg_dir   = "gust"
diff --git a/TestCases/unsteady/square_cylinder/run_aeroacoustics.py b/TestCases/unsteady/square_cylinder/run_aeroacoustics.py
new file mode 100644
index 00000000000..48a6c58d91e
--- /dev/null
+++ b/TestCases/unsteady/square_cylinder/run_aeroacoustics.py
@@ -0,0 +1,86 @@
+#!/usr/bin/env python3
+
+## \file run_aeroacoustics.py
+#  \brief Regression driver for aeroacoustic pressure-history analysis.
+#  \version 8.5.0 "Harrier"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+#
+# SU2 is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# SU2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+import argparse
+import os
+import sys
+import tempfile
+from pathlib import Path
+from types import SimpleNamespace
+
+
+su2_run = os.environ.get("SU2_RUN")
+if su2_run:
+    sys.path.insert(0, su2_run)
+
+from aeroacoustics import analyze  # noqa: E402
+
+
+SEGMENT_LENGTH = 2048
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Run the square-cylinder aeroacoustics regression analysis."
+    )
+    parser.add_argument("history", help="SU2 pressure-history CSV")
+    args = parser.parse_args()
+
+    with tempfile.TemporaryDirectory() as directory:
+        summary = analyze(
+            SimpleNamespace(
+                input=Path(args.history),
+                pressure=["mic1"],
+                sample_rate=None,
+                time="Cur_Time",
+                pressure_scale=1.0,
+                reference_pressure=20.0e-6,
+                skip_samples=0,
+                segment_length=SEGMENT_LENGTH,
+                overlap=0.5,
+                window="hann",
+                detrend="mean",
+                output=Path(directory) / "spectrum.csv",
+                summary=None,
+            )
+        )
+
+    metrics = summary["signals"]["mic1"]
+    print("\n------------------------------ Begin Solver -----------------------------\n")
+    print(
+        "| 0 | {:.10f} | {:.10f} |".format(
+            metrics["oaspl_db"], metrics["dominant_frequency_hz"]
+        )
+    )
+    print(
+        "Analyzed {} pressure samples using {} Welch segments.".format(
+            summary["samples"], metrics["segments"]
+        )
+    )
+
+
+if __name__ == "__main__":
+    main()