#  Copyright (c) 2026. Jose M. Requena-Plens
"""
Speech Transmission Index (STI) per IEC 60268-16:2020 (Edition 5).

Implements the full-STI indirect method from impulse responses (Schroeder
modulation transfer function), the direct STIPA method on recorded signals
(Annex B) and an Ed.5-conformant STIPA test-signal generator (clauses A.4
and A.6.1). Only the male speech option exists: Edition 5 removed the
female spectrum and weighting factors (foreword, item d).

The computation chain (octave-band MTF -> auditory masking and reception
threshold correction -> effective SNR clipped to +/-15 dB -> transmission
indices -> band MTI -> weighted STI) is numerically identical between
Ed.4 (2011) clauses A.5.2-A.5.6 and Ed.5; the only Ed.5 numeric change is
the male test-signal spectrum (A.6.1).
"""

from __future__ import annotations

import warnings
from dataclasses import dataclass
from typing import TYPE_CHECKING, Any, List, Sequence

import numpy as np

if TYPE_CHECKING:
    from matplotlib.axes import Axes
from scipy import signal

from ._warnings import PhonometryWarning
from .core import OctaveFilterBank
from .frequencies import nominal_frequencies
from .utils import _typesignal


class STIWarning(PhonometryWarning):
    """Warns about suspect STI/STIPA measurements or inputs."""


# Octave-band analysis range: 125 Hz - 8 kHz (Ed.4 A.5.1 = Ed.5), 7 bands.
_BAND_LIMITS = (125.0, 8000.0)
_NUM_BANDS = 7

# The 8 kHz octave band extends to 11.2 kHz (IEC 61260-1 exact upper edge
# 7943.3 * 2^0.5 = 11220 Hz), so fs/2 must exceed it.
_MIN_FS = 22500

# The 14 full-STI modulation frequencies, 0.63-12.5 Hz in nominal
# one-third-octave steps (Ed.4 A.2.2 = Ed.5).
_MOD_FREQS = np.array(
    [0.63, 0.80, 1.00, 1.25, 1.60, 2.00, 2.50, 3.15, 4.00, 5.00, 6.30, 8.00, 10.0, 12.5]
)

# Male octave-band weighting factors alpha and redundancy factors beta
# (Ed.5 Table A.1; identical to the male rows of Ed.4 Table A.3).
_ALPHA_MALE = np.array([0.085, 0.127, 0.230, 0.233, 0.309, 0.224, 0.173])
_BETA_MALE = np.array([0.085, 0.078, 0.065, 0.011, 0.047, 0.095])

# Absolute speech reception threshold ART_k in dB SPL per octave band
# (Ed.5 Table A.3; identical to Ed.4 Table A.2).
_ART_DB = np.array([46.0, 27.0, 12.0, 6.5, 7.5, 8.0, 12.0])

# STIPA modulation frequency pairs per octave band (Table B.1, unchanged
# between Ed.4 and Ed.5) and the source modulation index per component.
_STIPA_F1 = np.array([1.60, 1.00, 0.63, 2.00, 1.25, 0.80, 2.50])
_STIPA_F2 = np.array([8.00, 5.00, 3.15, 10.00, 6.25, 4.00, 12.50])
_STIPA_MOD_INDEX = 0.55  # Annex B: m = 0.55 per component, unchanged in Ed.5

#: Recommended minimum STIPA recording length (IEC 60268-16 STIPA practice
#: recommends 15 s to 25 s). Below this the recovered modulation depths are
#: biased low - and hence the STI too - because the slow modulation
#: components are averaged over too few periods: on an ideal loopback the
#: recovered STI is ~0.944 at 5 s, ~0.991 at 15 s and ~0.998 at 18 s.
_STIPA_MIN_SECONDS = 15.0

# Ed.5 male speech test-signal spectrum, clause A.6.1, in dB relative to
# the 500 Hz band (125 Hz and 250 Hz reduced by 6.2 dB and 3.2 dB vs Ed.4
# Table A.4 - the only numeric change of Edition 5).
_SPEECH_SPECTRUM_ED5 = np.array([-2.5, 0.5, 0.0, -6.0, -12.0, -18.0, -24.0])

# Intensity-envelope low-pass cut-off, ~100 Hz (Ed.4 A.5.2 = Ed.5).
_ENVELOPE_LPF_HZ = 100.0

# Annex F (informative) qualification bands: edges 0.36-0.76 in 0.04
# steps; categories U (< 0.36), J..A, A+ (>= 0.76). Unchanged in Ed.5.
_RATING_EDGES = (0.36, 0.40, 0.44, 0.48, 0.52, 0.56, 0.60, 0.64, 0.68, 0.72, 0.76)
_RATING_LETTERS = ("U", "J", "I", "H", "G", "F", "E", "D", "C", "B", "A", "A+")


@dataclass(frozen=True)
class STIResult:
    """Result of a Speech Transmission Index computation.

    ``mtf`` holds the modulation transfer values actually used for the
    transmission indices, i.e. after the optional SNR / masking /
    reception-threshold corrections and after clipping to [0, 1]; its
    shape is (7, 14) for full STI and (7, 2) for STIPA. ``mti`` is the
    per-band modulation transfer index (7,), ``band_levels`` echoes the
    speech octave-band levels used for the level-dependent corrections
    (None when they were skipped) and ``rating`` is the Annex F
    qualification letter (``A+`` .. ``U``).
    """

    sti: float
    mti: np.ndarray
    mtf: np.ndarray
    band_levels: np.ndarray | None
    rating: str

    def plot(self, ax: Axes | None = None, **kwargs: Any) -> Axes:
        """Plot the per-band MTI bars with the STI and rating letter.

        Requires matplotlib (``pip install phonometry[plot]``); returns the
        :class:`~matplotlib.axes.Axes`.
        """
        from ._plotting import plot_sti

        return plot_sti(self, ax=ax, **kwargs)


def _rating(sti: float) -> str:
    """Annex F qualification letter for an STI value (band edges 0.36-0.76).

    Note: the qualification bands are 0.04 STI wide, so every band edge flips a
    letter under +/-0.005 of STI jitter. Because typical STIPA measurement
    uncertainty (+/-0.02 to 0.03) exceeds the 0.02 half-width, a value near an
    edge can legitimately be graded one letter either way; treat the letter as
    indicative near a boundary and report the numeric STI alongside it.
    """
    idx = int(np.searchsorted(np.asarray(_RATING_EDGES), sti, side="right"))
    return _RATING_LETTERS[idx]


def _masking_amdb(level: np.ndarray | float) -> np.ndarray:
    """Level-dependent auditory masking ``amdB`` in dB (Ed.5 Table A.2).

    Four-segment continuous function of the total (signal + noise) level
    ``L`` of the next lower octave band, identical to Ed.4 Table A.1:
    ``0,5*L - 65`` (L < 63), ``1,8*L - 146,9`` (63 <= L < 67),
    ``0,5*L - 59,8`` (67 <= L < 100) and ``-10`` (L >= 100).
    """
    lvl = np.asarray(level, dtype=np.float64)
    return np.asarray(
        np.select(
            [lvl < 63.0, lvl < 67.0, lvl < 100.0],
            [0.5 * lvl - 65.0, 1.8 * lvl - 146.9, 0.5 * lvl - 59.8],
            default=-10.0,
        )
    )


def _octave_bank(fs: int) -> OctaveFilterBank:
    """Octave filter bank for the 7 STI bands (IEC 61260-1 class 0/1)."""
    bank = OctaveFilterBank(fs=fs, fraction=1, order=6, limits=list(_BAND_LIMITS))
    if bank.num_bands != _NUM_BANDS:
        raise ValueError(
            f"Expected {_NUM_BANDS} octave bands between 125 Hz and 8 kHz, "
            f"got {bank.num_bands} at fs={fs}."
        )
    return bank


def _validate_band_vector(values: Sequence[float] | np.ndarray, name: str) -> np.ndarray:
    arr = np.asarray(values, dtype=np.float64)
    if arr.shape != (_NUM_BANDS,):
        raise ValueError(
            f"'{name}' must contain exactly {_NUM_BANDS} octave-band values "
            f"(125 Hz - 8 kHz), got shape {arr.shape}."
        )
    return arr


def _sti_from_mtf(
    mtf: np.ndarray,
    snr: float | Sequence[float] | np.ndarray | None = None,
    level: Sequence[float] | np.ndarray | None = None,
    ambient: Sequence[float] | np.ndarray | None = None,
) -> STIResult:
    """STI computation chain from a matrix of modulation transfer values.

    Applies, in order: m > 1 truncation (Ed.4 A.5.3 NOTE 1), the optional
    signal-to-noise and level-dependent corrections (Ed.4 A.5.3 = Ed.5),
    the effective SNR with +/-15 dB limits (A.5.4), the transmission
    indices TI = (SNR_eff + 15)/30 (A.5.5), the band MTIs and the final
    weighted STI truncated to 1,0 (A.5.6, male factors of Ed.5 Table A.1).

    Noise handling: ``snr`` alone multiplies m by 1/(1 + 10^(-SNR/10)),
    which is exactly the I_k/(I_k + I_n,k) factor of the standard. When
    ``level`` is provided the full intensity-domain correction
    ``m' = m * I_k / (I_k + I_am,k + I_rt,k + I_n,k)`` is used instead,
    with ``ambient`` (or ``level - snr``) defining I_n,k, so the noise
    degradation is never applied twice.
    """
    m = np.array(mtf, dtype=np.float64)
    if m.ndim != 2 or m.shape[0] != _NUM_BANDS:
        raise ValueError(
            f"'mtf' must have shape ({_NUM_BANDS}, n_modulation_frequencies), "
            f"got {m.shape}."
        )
    if np.any(m < 0.0) or not np.all(np.isfinite(m)):
        raise ValueError("Modulation transfer values must be finite and >= 0.")
    if np.any(m > 1.3):
        # Ed.4 A.5.3 NOTE 1 (= Ed.5): m > 1,3 indicates an invalid measurement.
        warnings.warn(
            "Modulation transfer values above 1.3 detected: the measurement "
            "is likely invalid (IEC 60268-16 A.5.3). Values truncated to 1.0.",
            STIWarning,
            stacklevel=3,
        )
    m = np.minimum(m, 1.0)

    if snr is not None and ambient is not None:
        raise ValueError("Provide either 'snr' or 'ambient' noise levels, not both.")

    snr_arr: np.ndarray | None = None
    if snr is not None:
        snr_arr = np.asarray(snr, dtype=np.float64)
        if snr_arr.ndim == 0:
            snr_arr = np.full(_NUM_BANDS, float(snr_arr))
        elif snr_arr.shape != (_NUM_BANDS,):
            raise ValueError(
                f"'snr' must be a scalar or a vector of {_NUM_BANDS} "
                f"octave-band values, got shape {snr_arr.shape}."
            )

    band_levels: np.ndarray | None = None
    if level is None:
        if ambient is not None:
            raise ValueError(
                "'ambient' requires the speech octave-band levels 'level' to "
                "form the intensity-domain correction; pass 'snr' instead."
            )
        if snr_arr is not None:
            # I_k / (I_k + I_n,k) with I_n,k = I_k * 10^(-SNR/10).
            m = m / (1.0 + 10.0 ** (-snr_arr[:, np.newaxis] / 10.0))
        # No absolute level information: the auditory masking and absolute
        # reception threshold corrections are skipped.
    else:
        band_levels = _validate_band_vector(level, "level")
        i_signal = 10.0 ** (band_levels / 10.0)
        if ambient is not None:
            ambient_arr = _validate_band_vector(ambient, "ambient")
        elif snr_arr is not None:
            ambient_arr = band_levels - snr_arr
        else:
            ambient_arr = None
        i_noise = (
            10.0 ** (ambient_arr / 10.0) if ambient_arr is not None else np.zeros(_NUM_BANDS)
        )
        i_total = i_signal + i_noise
        level_total = 10.0 * np.log10(i_total)
        # Masking only acts on the next higher band; 125 Hz is unmasked.
        i_masking = np.zeros(_NUM_BANDS)
        i_masking[1:] = i_total[:-1] * 10.0 ** (_masking_amdb(level_total[:-1]) / 10.0)
        i_threshold = 10.0 ** (_ART_DB / 10.0)
        factor = i_signal / (i_signal + i_masking + i_threshold + i_noise)
        m = m * factor[:, np.newaxis]

    # Effective SNR clipped to +/-15 dB (A.5.4); m = 0 and m = 1 map to the
    # clip limits through the log divergences.
    with np.errstate(divide="ignore"):
        snr_eff = 10.0 * np.log10(m / (1.0 - m))
    snr_eff = np.clip(snr_eff, -15.0, 15.0)
    ti = (snr_eff + 15.0) / 30.0  # A.5.5
    mti = ti.mean(axis=1)

    # A.5.6 with the Ed.5 Table A.1 male factors; truncated to 1,0 (the
    # known alpha/beta artefact can push the raw sum above 1).
    sti = float(
        np.dot(_ALPHA_MALE, mti) - np.dot(_BETA_MALE, np.sqrt(mti[:-1] * mti[1:]))
    )
    sti = min(sti, 1.0)
    return STIResult(
        sti=sti, mti=mti, mtf=m, band_levels=band_levels, rating=_rating(sti)
    )


def sti_from_impulse_response(
    ir: List[float] | np.ndarray,
    fs: int,
    snr: float | Sequence[float] | np.ndarray | None = None,
    level: Sequence[float] | np.ndarray | None = None,
    ambient: Sequence[float] | np.ndarray | None = None,
) -> STIResult:
    """
    Full STI from a room/system impulse response (indirect method).

    The impulse response is filtered into the seven octave bands 125 Hz -
    8 kHz (IEC 61260-1 filters) and the modulation transfer function is
    obtained from the Schroeder integral (IEC 60268-16, indirect method):
    ``m_k(f_m) = |integral h_k^2(t) exp(-j 2 pi f_m t) dt| /
    integral h_k^2(t) dt`` at the 14 modulation frequencies 0,63-12,5 Hz
    (A.2.2). The result then follows the standard chain: optional noise
    degradation, optional auditory masking and absolute reception
    threshold correction, effective SNR clipped to +/-15 dB, transmission
    indices, band MTIs and the male-weighted STI (Ed.5 Table A.1).

    When neither ``level`` nor ``ambient`` is given the level-dependent
    auditory masking and the absolute reception threshold corrections are
    skipped (they require absolute band levels), matching the common
    "noise-free indirect measurement" use of the standard.

    .. note:: The indirect method has a small positive MTF bias that grows
       with reverberation time for a finite noise-carrier IR. It stays within
       the IEC 60268-16 A.5.1.2 systematic-error allowance (<= 0.01 STI) up to
       about T60 = 4 s, reaching ~+0.012 STI only in very reverberant rooms
       (T60 ~ 8 s, STI ~ 0.19). It is a property of the finite IR, not of the
       (exact) Schroeder integration, and does not depend on IR truncation.

    :param ir: Impulse response (1D).
    :param fs: Sample rate in Hz (>= 22,5 kHz so the 8 kHz band fits).
    :param snr: Optional signal-to-noise ratio in dB, scalar or one value
        per octave band. Degrades m by 1/(1 + 10^(-SNR/10)); combined
        with ``level`` it is interpreted as ambient levels
        ``level - snr`` so noise is not applied twice. Mutually
        exclusive with ``ambient``.
    :param level: Optional speech octave-band levels in dB SPL (7 values)
        at the listener position; enables the auditory masking (Ed.5
        Table A.2) and reception threshold (Ed.5 Table A.3) corrections.
    :param ambient: Optional ambient noise octave-band levels in dB SPL
        (7 values); requires ``level``.
    :return: :class:`STIResult` with ``mtf`` of shape (7, 14).
    """
    ir_proc = _typesignal(ir)
    if ir_proc.ndim != 1:
        raise ValueError("sti_from_impulse_response expects a 1D impulse response.")
    if fs <= 0:
        raise ValueError("Sample rate 'fs' must be positive.")
    if fs < _MIN_FS:
        raise ValueError(
            f"Sample rate 'fs' must be >= {_MIN_FS} Hz: the 8 kHz octave band "
            "extends to 11,2 kHz (IEC 61260-1)."
        )
    if not np.any(ir_proc):
        raise ValueError("Impulse response 'ir' is silent.")

    bank = _octave_bank(fs)
    _, _, bands = bank.filter(
        ir_proc, sigbands=True, detrend=False, calculate_level=False
    )
    p_bands = np.asarray(bands) ** 2  # h_k^2(t), shape (7, n)
    denom = p_bands.sum(axis=1)
    if np.any(denom <= 0.0):
        raise ValueError("An octave band of the impulse response has no energy.")

    t = np.arange(ir_proc.shape[-1]) / fs
    mtf = np.empty((_NUM_BANDS, _MOD_FREQS.size))
    for j, fm in enumerate(_MOD_FREQS):
        kernel = np.exp(-2j * np.pi * fm * t)
        mtf[:, j] = np.abs(p_bands @ kernel) / denom
    return _sti_from_mtf(mtf, snr=snr, level=level, ambient=ambient)


def _intensity_envelopes(x: np.ndarray, fs: int) -> np.ndarray:
    """Octave-band intensity envelopes I_k(t) = x_k^2(t) low-passed at
    ~100 Hz (Ed.4 A.5.2 = Ed.5), shape (7, n)."""
    bank = _octave_bank(fs)
    _, _, bands = bank.filter(x, sigbands=True, calculate_level=False)
    sos = signal.butter(4, _ENVELOPE_LPF_HZ, btype="lowpass", fs=fs, output="sos")
    env = np.empty((_NUM_BANDS, x.shape[-1]))
    for k, band in enumerate(bands):
        env[k] = signal.sosfiltfilt(sos, band**2)
    return env


def _stipa_modulation_depths(env: np.ndarray, fs: int) -> np.ndarray:
    """Modulation depths mdr_{k,fm} of the intensity envelopes at the two
    Table B.1 frequencies of each band, over an integer number of periods:
    ``mdr = 2 sqrt[(sum I sin)^2 + (sum I cos)^2] / sum I`` (Ed.4 A.5.2)."""
    n = env.shape[-1]
    duration = n / fs
    mdr = np.empty((_NUM_BANDS, 2))
    for k in range(_NUM_BANDS):
        for j, fm in enumerate((float(_STIPA_F1[k]), float(_STIPA_F2[k]))):
            periods = int(np.floor(duration * fm))
            if periods < 1:
                raise ValueError(
                    f"Signal too short for STIPA: it must contain at least one "
                    f"full period of the {fm} Hz modulation (>= {1.0 / fm:.2f} s; "
                    "the standard recommends 15 s to 25 s)."
                )
            n_use = int(round(periods / fm * fs))
            seg = env[k, :n_use]
            total = float(np.sum(seg))
            if total <= 0.0:
                raise ValueError(
                    "An octave band of the signal has no energy; STIPA needs "
                    "all seven carriers of the test signal."
                )
            t = np.arange(n_use) / fs
            s = float(np.dot(seg, np.sin(2.0 * np.pi * fm * t)))
            c = float(np.dot(seg, np.cos(2.0 * np.pi * fm * t)))
            mdr[k, j] = 2.0 * float(np.hypot(s, c)) / total
    return mdr


def stipa(
    x: List[float] | np.ndarray,
    fs: int,
    reference: List[float] | np.ndarray | None = None,
    level: Sequence[float] | np.ndarray | None = None,
    ambient: Sequence[float] | np.ndarray | None = None,
) -> STIResult:
    """
    STIPA on a recorded test signal (direct method, Annex B).

    The recording is filtered into the seven octave bands, squared and
    low-passed (~100 Hz) into intensity envelopes, and the modulation
    depths at the two Table B.1 modulation frequencies of each band are
    measured with the sine/cosine correlation over an integer number of
    periods (Ed.4 A.5.2 = Ed.5). The modulation transfer values are the
    measured depths normalized by the source modulation index 0,55
    (Annex B) - or by the depths measured on ``reference`` when the
    actually emitted signal is supplied - and feed the same masking /
    threshold / TI / STI chain as the full method.

    Physical background noise is already contained in the recording; use
    ``level`` (and optionally ``ambient``) only to enable the absolute
    level-dependent corrections, which are otherwise skipped.

    An :class:`STIWarning` is emitted when the recording is shorter than
    the recommended 15 s (IEC 60268-16 STIPA practice, 15 s to 25 s),
    because the slow modulation components are then averaged over too few
    periods and the recovered modulation depths - and hence the STI - are
    biased low (an ideal loopback gives STI ~0.944 at 5 s vs ~0.998 at 18 s).

    :param x: Recorded STIPA signal (1D), 15 s to 25 s recommended.
    :param fs: Sample rate in Hz (>= 22,5 kHz).
    :param reference: Optional reference recording of the undistorted
        test signal; its measured modulation depths replace the nominal
        0,55 as normalization (useful for non-conformant sources).
    :param level: Optional speech octave-band levels in dB SPL (7 values)
        enabling auditory masking and reception threshold corrections.
    :param ambient: Optional ambient noise octave-band levels in dB SPL
        (7 values); requires ``level``.
    :return: :class:`STIResult` with ``mtf`` of shape (7, 2).
    """
    x_proc = _typesignal(x)
    if x_proc.ndim != 1:
        raise ValueError("stipa expects a 1D signal.")
    if fs <= 0:
        raise ValueError("Sample rate 'fs' must be positive.")
    if fs < _MIN_FS:
        raise ValueError(
            f"Sample rate 'fs' must be >= {_MIN_FS} Hz: the 8 kHz octave band "
            "extends to 11,2 kHz (IEC 61260-1)."
        )

    if x_proc.size / fs < _STIPA_MIN_SECONDS:
        warnings.warn(
            f"STIPA recording is {x_proc.size / fs:.1f} s, shorter than the "
            f"recommended {_STIPA_MIN_SECONDS:.0f} s (IEC 60268-16 STIPA "
            "practice, 15 s to 25 s): the recovered modulation depths, and "
            "hence the STI, are biased low.",
            STIWarning,
            stacklevel=2,
        )

    mdr = _stipa_modulation_depths(_intensity_envelopes(x_proc, fs), fs)
    if reference is not None:
        ref_proc = _typesignal(reference)
        if ref_proc.ndim != 1:
            raise ValueError("'reference' must be a 1D signal.")
        mdt = _stipa_modulation_depths(_intensity_envelopes(ref_proc, fs), fs)
        if np.any(mdt <= 0.0):
            raise ValueError("'reference' has zero modulation depth in a band.")
        mtf = mdr / mdt
    else:
        mtf = mdr / _STIPA_MOD_INDEX
    return _sti_from_mtf(mtf, level=level, ambient=ambient)


def _pink_noise(n: int, rng: np.random.Generator) -> np.ndarray:
    """Gaussian pink (1/f power) noise of length ``n`` via FFT shaping."""
    spectrum = np.fft.rfft(rng.standard_normal(n))
    freq = np.fft.rfftfreq(n)
    freq[0] = freq[1]  # placeholder; DC is zeroed below
    spectrum /= np.sqrt(freq)
    spectrum[0] = 0.0
    return np.asarray(np.fft.irfft(spectrum, n))


def stipa_signal(
    fs: int,
    seconds: float = 18.0,
    level_db: float | None = None,
    seed: int | None = None,
) -> np.ndarray:
    """
    Generate an IEC 60268-16:2020 conformant STIPA test signal.

    Pink-noise carriers are band-limited to half-octave bands centred on
    the seven octave-band frequencies (clause A.4), set to the Ed.5 male
    speech spectrum of clause A.6.1 (``-2,5; 0,5; 0; -6; -12; -18; -24``
    dB re the 500 Hz band) and intensity-modulated with
    ``0,5 (1 + 0,55 (sin 2 pi f1 t - sin 2 pi f2 t))`` - the Table B.1
    frequency pair of each band, 180 degrees between components, applied
    in amplitude through its square root (Annex B).

    :param fs: Sample rate in Hz (>= 22,5 kHz).
    :param seconds: Duration in seconds (the standard recommends 15 s to
        25 s; default 18 s).
    :param level_db: Optional overall level in dB SPL: the output is
        scaled so its RMS, taken as pascals, sits at ``level_db`` re
        20 uPa. Default (None) normalizes the RMS to 0,1 (digital full
        scale headroom for the 12-14 dB crest factor).
    :param seed: Seed for the pink-noise generator (None: random).
    :return: Test signal, 1D array of ``round(seconds * fs)`` samples.
    """
    if fs <= 0:
        raise ValueError("Sample rate 'fs' must be positive.")
    if fs < _MIN_FS:
        raise ValueError(
            f"Sample rate 'fs' must be >= {_MIN_FS} Hz: the 8 kHz half-octave "
            "carrier extends beyond 9,4 kHz."
        )
    if seconds <= 0:
        raise ValueError("'seconds' must be positive.")

    n = int(round(seconds * fs))
    rng = np.random.default_rng(seed)
    pink = _pink_noise(n, rng)
    t = np.arange(n) / fs
    centers, _, _, _ = nominal_frequencies(1, list(_BAND_LIMITS))
    half_octave = 2.0**0.25
    out = np.zeros(n)
    for k in range(_NUM_BANDS):
        sos = signal.butter(
            5,
            [centers[k] / half_octave, centers[k] * half_octave],
            btype="bandpass",
            fs=fs,
            output="sos",
        )
        carrier = signal.sosfiltfilt(sos, pink)
        carrier = carrier / np.sqrt(np.mean(carrier**2))
        envelope = 0.5 * (
            1.0
            + _STIPA_MOD_INDEX
            * (np.sin(2.0 * np.pi * _STIPA_F1[k] * t) - np.sin(2.0 * np.pi * _STIPA_F2[k] * t))
        )
        # Intensity modulation: amplitude carries the square root. The
        # exact 1:5 frequency ratios keep the envelope non-negative; the
        # clip only guards against float round-off.
        out += (
            10.0 ** (_SPEECH_SPECTRUM_ED5[k] / 20.0)
            * carrier
            * np.sqrt(np.maximum(envelope, 0.0))
        )

    rms = float(np.sqrt(np.mean(out**2)))
    if level_db is None:
        return out * (0.1 / rms)
    return np.asarray(out * (2e-5 * 10.0 ** (level_db / 20.0) / rms))
