def VEST_CoilCurrentNoiseReduction(data):  # noqa: N802
    """Deprecated compatibility wrapper for the snake-case function."""
    warnings.warn(
        "VEST_CoilCurrentNoiseReduction() is deprecated; use "
        "vest_coil_current_noise_reduction().",
        DeprecationWarning,
        stacklevel=2,
    )
    return vest_coil_current_noise_reduction(data)

def butterworth_bandpass(data, low: float, high: float, fs: float, order: int = 2,
                         *, zero_phase: bool = True) -> np.ndarray:
    """Butterworth band-pass filter along the last axis.

    Parameters
    ----------
    data : array_like
        The waveform or a stack of them; the time axis is the last [any].
    low : float
        Lower band edge [Hz].
    high : float
        Upper band edge, above ``low`` and below the Nyquist frequency [Hz].
    fs : float
        Sample rate [Hz].
    order : int, optional
        Filter order [-].
    zero_phase : bool, optional
        ``True`` for forward-backward ``filtfilt``; ``False`` for a causal
        ``lfilter`` [-].

    Returns
    -------
    np.ndarray
        The filtered waveform, same shape as ``data`` [any].

    Raises
    ------
    ValueError
        The band edges do not satisfy ``0 < low < high < fs / 2``.

    Defaults
    --------
    ``order = 2`` is a validated-workflow default shared with
    :func:`butterworth_lowpass`.  ``zero_phase = True`` is the opposite of the
    low-pass default because band-passed signals feed spectral and
    mode-number analysis, where phase across channels is the measurement.

    Convention
    ----------
    Zero-phase by default; see :func:`butterworth_lowpass` for the causal
    alternative and when each is right.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    As :func:`butterworth_lowpass`: a record too short for ``filtfilt`` raises
    from SciPy.
    """
    values = np.asarray(data, dtype=float)
    nyquist = 0.5 * float(fs)
    if not 0.0 < float(low) < float(high) < nyquist:
        raise ValueError(
            f"band edges must satisfy 0 < low < high < {nyquist:g} Hz for fs={fs:g}; "
            f"got ({low!r}, {high!r})"
        )
    b, a = scipy_signal.butter(
        int(order), [float(low) / nyquist, float(high) / nyquist], btype="band"
    )
    if zero_phase:
        return scipy_signal.filtfilt(b, a, values, axis=-1)
    return scipy_signal.lfilter(b, a, values, axis=-1)

def butterworth_lowpass(data, cutoff: float, fs: float, order: int = 2,
                        *, zero_phase: bool = False) -> np.ndarray:
    """Butterworth low-pass filter along the last axis.

    Parameters
    ----------
    data : array_like
        The waveform or a stack of them; the time axis is the last [any].
    cutoff : float
        Cutoff frequency, strictly between 0 and the Nyquist frequency [Hz].
    fs : float
        Sample rate [Hz].
    order : int, optional
        Filter order [-].
    zero_phase : bool, optional
        ``True`` for forward-backward ``filtfilt``; ``False`` for a causal
        ``lfilter`` [-].

    Returns
    -------
    np.ndarray
        The filtered waveform, same shape as ``data`` [any].

    Raises
    ------
    ValueError
        ``cutoff`` is not inside ``(0, fs / 2)``.

    Defaults
    --------
    ``order = 2`` and ``zero_phase = False`` are validated-workflow defaults:
    the values the VEST SXR viewer uses [SXR]_, kept so that
    :mod:`vaft.process.soft_x_rays` reproduces it.

    Convention
    ----------
    **Causal by default.**  With ``zero_phase=False`` every feature is delayed
    by the filter's group delay, identically on every channel, which is what
    ratio and reference arithmetic between SXR channels needs.  Pass
    ``zero_phase=True`` wherever a time is going to be read off the result --
    an onset, a peak -- because a causal filter moves it.

    Applicability
    -------------
    Machine-independent.  The defaults are VEST SXR practice but embed no VEST
    constant.

    Limitations
    -----------
    ``filtfilt`` needs a record longer than three times the filter's padding
    length and raises from SciPy otherwise; nothing here catches that.
    """
    values = np.asarray(data, dtype=float)
    nyquist = 0.5 * float(fs)
    if not 0.0 < float(cutoff) < nyquist:
        raise ValueError(
            f"cutoff must lie in (0, {nyquist:g}) Hz for fs={fs:g}; got {cutoff!r}"
        )
    b, a = scipy_signal.butter(int(order), float(cutoff) / nyquist, btype="low")
    if zero_phase:
        return scipy_signal.filtfilt(b, a, values, axis=-1)
    return scipy_signal.lfilter(b, a, values, axis=-1)

def define_baseline(time, onset_time, onset_window, offset_time=None, offset_window=None):
    """Select the baseline samples around a signal's onset and optional offset.

    Converts ``onset_time`` (and ``offset_time``) to indices with
    ``np.searchsorted`` and returns the ``onset_window`` samples immediately
    before the onset, plus the ``offset_window`` samples immediately after the
    offset when both are given.  The result is what :func:`subtract_baseline`
    fits through.

    Parameters
    ----------
    time : np.ndarray
        Sample times, sorted [s].
    onset_time : float
        Time at which the signal begins to deviate from baseline [s].
    onset_window : int
        Number of samples before the onset to include [-].
    offset_time : float or None, optional
        Time at which the signal returns to baseline [s].
    offset_window : int or None, optional
        Number of samples after the offset to include; both offset arguments
        must be given for the offset region to be used [-].

    Returns
    -------
    np.ndarray
        Indices into ``time`` of the baseline samples, onset region first [-].

    Assumptions
    -----------
    The signal is genuinely quiet in the selected windows.  A window that
    reaches back into a previous pulse, or forward into pickup from a coil
    ramp, is fitted as if it were baseline.

    Applicability
    -------------
    Machine-independent.  The onset and offset times, and the window widths,
    are the caller's; on VEST they come from the magnetics processing window
    in ``vest.yaml`` via :mod:`vaft.process.magnetics`.

    Limitations
    -----------
    Windows are clipped to the record silently, so a window wider than the
    pre-onset record yields fewer samples than asked, with no warning.

    Provenance
    ----------
    .. [1] The legacy VEST EFIT magnetics chain [VFIT]_, where the baseline is
       taken before the discharge and after it and removed by a fit.
    """
    baseline_indices = []

    # Convert onset_time -> onset_idx
    onset_idx = np.searchsorted(time, onset_time)

    # Add onset window
    if onset_window > 0:
        start_idx = max(0, onset_idx - onset_window)
        baseline_indices.extend(range(start_idx, onset_idx))

    # Convert offset_time -> offset_idx if provided
    if offset_time is not None and offset_window is not None:
        offset_idx = np.searchsorted(time, offset_time)
        end_idx = min(len(time), offset_idx + offset_window)
        baseline_indices.extend(range(offset_idx, end_idx))

    return np.array(baseline_indices)

def describe_time_grid(time, *, rtol: float = 1e-3) -> TimeGrid:
    """Measure a timebase's spacing, uniformity and monotonicity.

    Parameters
    ----------
    time : array_like
        Sample times, flattened [s].
    rtol : float, optional
        Largest ``|dt_i - median(dt)| / median(dt)`` for which the grid is
        still reported uniform [-].

    Returns
    -------
    TimeGrid
        ``n``, ``t0``, the median ``dt``, ``sample_rate``, whether the grid is
        ``uniform`` within ``rtol`` and ``strictly_increasing``, and the
        ``max_relative_jitter`` actually measured [any].

    Defaults
    --------
    ``rtol = 1e-3`` is a numerical convenience sized to an acquisition-era
    fact: VEST shots after 42190 store their timebase as ``linspace(0, span,
    n)``, so a nominally 4e-6 s grid is really ``span / (n - 1)`` =
    4.00016e-6 s.  That is uniform for every purpose here and a tight
    tolerance would only reject it.

    Convention
    ----------
    ``dt`` is the *median* spacing, robust to the single ragged interval a
    concatenated acquisition can leave behind.  ``dt`` and ``sample_rate``
    are ``nan`` for a single-sample grid, where neither is defined; an empty
    grid reports ``n = 0`` and ``nan`` throughout.

    Applicability
    -------------
    Machine-independent.
    """
    values = np.asarray(time, dtype=float).reshape(-1)
    n = int(values.size)
    if n == 0:
        return TimeGrid(0, float("nan"), float("nan"), float("nan"), False, False, float("nan"))
    t0 = float(values[0])
    if n == 1:
        return TimeGrid(1, t0, float("nan"), float("nan"), True, True, float("nan"))

    steps = np.diff(values)
    dt = float(np.median(steps))
    strictly_increasing = bool(np.all(steps > 0.0))
    if dt == 0.0 or not np.isfinite(dt):
        return TimeGrid(n, t0, dt, float("nan"), False, strictly_increasing, float("inf"))
    jitter = float(np.max(np.abs(steps - dt)) / abs(dt))
    return TimeGrid(
        n=n,
        t0=t0,
        dt=dt,
        sample_rate=1.0 / dt,
        uniform=bool(jitter <= float(rtol)),
        strictly_increasing=strictly_increasing,
        max_relative_jitter=jitter,
    )

def detect_clipped_samples(data, *, clip_values, tolerance: float) -> np.ndarray:
    """Return a boolean mask of samples sitting at an acquisition limit.

    ``clip_values`` is a single signed level or a sequence of them, and the
    mask is their union: a sample is saturated when it lies within
    ``tolerance`` of *any* supplied level.  Real acquisition hardware rails on
    both sides and rarely symmetrically -- VEST's diamagnetic Rogowski channel
    is a signed 16-bit ADC over +/-5 V, so its rails are exactly ``-5.0`` and
    ``5 * 32767 / 32768`` [1]_.

    Detecting every level in one pass is not a convenience.  Where a waveform
    oscillates hard enough to hit both rails within a few samples, repairing
    one rail at a time would fit the reconstruction through samples still
    pinned at the other one.

    Parameters
    ----------
    data : array_like
        The waveform, raw or calibrated; the rails are in its units [any].
    clip_values : float or sequence of float
        The acquisition limit or limits, signed, in the units of ``data`` [any].
    tolerance : float
        Half-width of the band around each limit inside which a sample counts
        as saturated [any].

    Returns
    -------
    np.ndarray
        Boolean mask, ``True`` where the sample is saturated; same shape as
        ``data`` [-].

    Raises
    ------
    SignalRepairError
        ``clip_values`` is empty or non-finite, or ``tolerance`` is not a
        positive finite width.

    Assumptions
    -----------
    Saturation shows as a sample *at* the rail, not as a distorted sample near
    it.  A digitizer that soft-clips will not be detected by a narrow
    tolerance.

    Applicability
    -------------
    Machine-independent.  The caller supplies the rails and tolerance; this
    function infers nothing about the hardware.  The VEST values are read from
    ``vest.yaml`` by :mod:`vaft.machine_mapping`, never here.

    Provenance
    ----------
    .. [1] ``vest.yaml``, the ``diamagnetic_flux`` processing block, where the
       asymmetric 16-bit rails are derived; issue #285.
    """
    values = np.asarray(data, dtype=float)
    levels = np.atleast_1d(np.asarray(clip_values, dtype=float)).reshape(-1)
    if levels.size == 0:
        raise SignalRepairError("`clip_values` must contain at least one acquisition limit.")
    if not np.all(np.isfinite(levels)):
        raise SignalRepairError(f"`clip_values` must be finite, got {clip_values!r}.")
    width = float(tolerance)
    if not np.isfinite(width) or width <= 0.0:
        raise SignalRepairError(f"`tolerance` must be a positive finite width, got {tolerance!r}.")

    saturated = np.zeros(values.shape, dtype=bool)
    for level in levels:
        saturated |= np.abs(values - level) < width
    return saturated

def detrend_moving_average(data, window_samples: int) -> np.ndarray:
    """Subtract a centered moving-average trend along the last axis.

    Parameters
    ----------
    data : array_like
        The waveform or a stack of them; the time axis is the last [any].
    window_samples : int
        Width of the trend window in samples; ``1`` or less removes nothing [-].

    Returns
    -------
    np.ndarray
        ``data`` minus its rolling mean, same shape [any].

    Convention
    ----------
    The trend is a centred rolling mean with ``min_periods=1`` semantics: the
    window at index ``i`` spans ``[i - w//2, i + (w-1)//2]`` clipped to the
    record, so edge windows shrink rather than yield NaN and the output has the
    input's length.  This is exactly ``pandas.Series.rolling(window,
    center=True, min_periods=1).mean()``, the convention of the validated VEST
    SXR viewer [SXR]_, reimplemented with cumulative sums so that
    :mod:`vaft.process.soft_x_rays` does not need pandas for it.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    A centred window is zero-phase but the shrinking edge windows bias the
    first and last ``w/2`` samples toward the local level.  The trend is a
    plain mean: one large excursion inside the window pulls the trend toward
    it and appears, inverted, in the neighbours.
    """
    values = np.asarray(data, dtype=float)
    window = int(window_samples)
    if window <= 1 or values.shape[-1] == 0:
        return values - values  # zero trend removal, preserving shape/dtype

    length = values.shape[-1]
    cumsum = np.zeros(values.shape[:-1] + (length + 1,), dtype=float)
    np.cumsum(values, axis=-1, out=cumsum[..., 1:])
    # A centered pandas window of size w at index i spans
    # [i - w//2, i + (w-1)//2], clipped to the record.
    index = np.arange(length)
    start = np.clip(index - window // 2, 0, length)
    stop = np.clip(index + (window - 1) // 2 + 1, 0, length)
    trend = (cumsum[..., stop] - cumsum[..., start]) / (stop - start)
    return values - trend

def exp_baseline(x, a, b, c):
    """Exponential model for baseline fitting: ``y = a * exp(b * x) + c``.

    Parameters
    ----------
    x : array_like
        Abscissa, normally time [any].
    a : float
        Amplitude [any].
    b : float
        Rate; the sign decides growth or decay [any].
    c : float
        Offset [any].

    Returns
    -------
    np.ndarray
        ``a * exp(b * x) + c`` [any].

    Applicability
    -------------
    Machine-independent.  A fit model for :func:`subtract_baseline`.
    """
    return a * np.exp(b * x) + c

def fir_filter(
    data,
    cutoff,
    fs: float,
    *,
    kind: str = "lowpass",
    numtaps: int | None = None,
    window: str = "hann",
    zero_phase: bool = True,
) -> np.ndarray:
    """Linear-phase FIR filter along the last axis.

    Parameters
    ----------
    data : array_like
        The waveform or a stack of them; the time axis is the last [any].
    cutoff : float or sequence of float
        Band edge, or the pair of edges for ``bandpass``/``bandstop`` [Hz].
    fs : float
        Sample rate [Hz].
    kind : str, optional
        ``lowpass``, ``highpass``, ``bandpass`` or ``bandstop`` [-].
    numtaps : int or None, optional
        Filter length; ``None`` derives it from ``fs``
        (:func:`fir_filter_coefficients`) [-].
    window : str, optional
        Window passed to :func:`scipy.signal.firwin` [-].
    zero_phase : bool, optional
        ``True`` for forward-backward ``filtfilt``; ``False`` for a causal
        ``lfilter`` that keeps the filter's own group delay [-].

    Returns
    -------
    np.ndarray
        The filtered waveform, same shape as ``data`` [any].

    Raises
    ------
    ValueError
        As :func:`fir_filter_coefficients`.

    Defaults
    --------
    ``zero_phase = True`` is a validated-workflow default matching
    :func:`butterworth_bandpass`: a filtered fluctuation signal feeds spectral
    and mode-number analysis, where phase across channels is the measurement.

    Convention
    ----------
    An FIR filter of this design is already linear-phase, so a causal
    application shifts every frequency by the same ``(numtaps - 1) / 2``
    samples and the delay can be removed exactly -- unlike an IIR filter, whose
    phase distortion is frequency-dependent and only ``filtfilt`` removes it.
    That is the reason to reach for FIR here rather than
    :func:`butterworth_bandpass`: the cost is a longer impulse response and
    more arithmetic per sample.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    ``filtfilt`` needs a record longer than roughly three times the tap count
    and raises from SciPy otherwise; a one-millisecond default response is long
    at a high sample rate.
    """
    values = np.asarray(data, dtype=float)
    taps = fir_filter_coefficients(
        cutoff, fs, kind=kind, numtaps=numtaps, window=window
    )
    if zero_phase:
        return scipy_signal.filtfilt(taps, 1.0, values, axis=-1)
    return scipy_signal.lfilter(taps, 1.0, values, axis=-1)

def is_signal_active(
    data,
    var_ratio_thresh=1e-2,
    change_ratio_thresh=1e-2,
    verbose=False,
):
    """Decide whether a trace carries a signal or is flat.

    Both tests are relative to the trace's own level, so the verdict does not
    depend on the signal's units or gain.

    Parameters
    ----------
    data : array_like
        The waveform; fewer than two samples is never active [any].
    var_ratio_thresh : float, optional
        Threshold on the variance divided by the squared mean absolute level,
        the squared coefficient of variation [-].
    change_ratio_thresh : float, optional
        Threshold on the mean absolute sample-to-sample change divided by the
        mean absolute level [-].
    verbose : bool, optional
        Print both ratios [-].

    Returns
    -------
    bool
        ``True`` unless *both* ratios fall below their thresholds [-].

    Defaults
    --------
    Both thresholds at ``1e-2`` are legacy compatibility values; no derivation
    is recorded.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    A trace held at a non-zero level with noise far below that level is
    inactive, whatever the level.  The converse is the caveat: the test is
    scale-free by design and has no absolute floor, so a zero-mean noise
    trace scores near 1 on both ratios and is reported active at any
    amplitude.  A caller that knows a noise floor must apply it before
    asking.  An identically zero trace is inactive.
    """
    data = np.asarray(data)

    if data.size < 2:
        return False

    variance = np.var(data)
    mean_abs_change = np.mean(np.abs(np.diff(data)))
    level = np.mean(np.abs(data))

    # Both ratios are relative to the magnitude of the trace (#463).  An
    # identically zero trace has no level to be relative to and is inactive.
    var_ratio = variance / level**2 if level > 0.0 else 0.0
    change_ratio = mean_abs_change / (level + 1e-12)

    if verbose:
        print(f"Variance ratio: {var_ratio:.3e} (thresh={var_ratio_thresh:.3e})")
        print(f"Mean |Δx| ratio: {change_ratio:.3e} (thresh={change_ratio_thresh:.3e})")

    if var_ratio < var_ratio_thresh and change_ratio < change_ratio_thresh:
        return False
    return True

def line_average_density(n_e_line, path_length_m: float) -> np.ndarray:
    """Return line-average electron density given an explicit chord length.

    Parameters
    ----------
    n_e_line : array_like
        Line-integrated electron density along the diagnostic chord [m^-2].
    path_length_m : float
        Length of the chord through the plasma [m].

    Returns
    -------
    np.ndarray
        Line-average electron density, ``n_e_line / path_length_m`` [m^-3].

    Raises
    ------
    ValueError
        ``path_length_m`` is not positive.

    Convention
    ----------
    Line-*integrated* in, line-*averaged* out.  An interferometer measures the
    former; the latter is what confinement scalings and Greenwald fractions
    want.  No calibration or geometry is inferred here -- the chord length is
    the caller's.

    Applicability
    -------------
    Machine-independent.  The chord length of a VEST interferometer channel
    comes from :mod:`vaft.machine_mapping`.
    """
    if path_length_m <= 0:
        raise ValueError(f"path_length_m must be positive, got {path_length_m}")
    return np.asarray(n_e_line, dtype=float) / float(path_length_m)

def linear_baseline(x, a, b):
    """Linear model for baseline fitting: ``y = a * x + b``.

    Parameters
    ----------
    x : array_like
        Abscissa, normally time [any].
    a : float
        Slope [any].
    b : float
        Intercept [any].

    Returns
    -------
    np.ndarray
        ``a * x + b`` [any].

    Applicability
    -------------
    Machine-independent.  A fit model for :func:`subtract_baseline`.
    """
    return a * x + b

def quadratic_baseline(x, a, b, c):
    """Quadratic model for baseline fitting: ``y = a * x^2 + b * x + c``.

    Parameters
    ----------
    x : array_like
        Abscissa, normally time [any].
    a : float
        Quadratic coefficient [any].
    b : float
        Linear coefficient [any].
    c : float
        Constant [any].

    Returns
    -------
    np.ndarray
        ``a * x**2 + b * x + c`` [any].

    Applicability
    -------------
    Machine-independent.  A fit model for :func:`subtract_baseline`.
    """
    return a * x**2 + b * x + c

def signal_on_offset(time, data, smooth_window=5, threshold=0.01, verbose=False):
    """Smooth a waveform, then return the window in which it is above a threshold.

    Parameters
    ----------
    time : array_like
        Sample times, one-dimensional [s].
    data : array_like
        The waveform, one-dimensional, same length as ``time`` [any].
    smooth_window : int, optional
        Savitzky-Golay window length in samples; must be odd and larger than
        the polynomial order, 3 [-].
    threshold : float, optional
        Level at or above which the smoothed signal counts as active, in the
        units of ``data`` [any].
    verbose : bool, optional
        Print the threshold [-].

    Returns
    -------
    tuple of (float, float)
        ``(onset, offset)`` times [s].

    Processing steps
    ----------------
    1. Smooth with ``scipy.signal.savgol_filter(data, smooth_window, 3)``.
    2. Hand the smoothed waveform to :func:`detect_active_window`.

    Defaults
    --------
    ``smooth_window = 5`` and the cubic polynomial are legacy compatibility
    values from the ported workflow, unrecorded beyond that.  ``threshold =
    0.01`` is inherited from :func:`detect_active_window`.

    Convention
    ----------
    Zero-phase: a Savitzky-Golay filter is a centred polynomial fit and moves
    nothing in time, so the returned onset is not delayed by the smoothing.

    Applicability
    -------------
    Machine-independent.  Callers in :mod:`vaft.omas.general` and
    :mod:`vaft.database` pass ``threshold=0.05`` for VEST H-alpha; the default
    is not a VEST policy.  For a plasma onset with evidence and a verdict, use
    :mod:`vaft.process.onset`.

    Limitations
    -----------
    Those of :func:`detect_active_window`: an absolute threshold and a walk
    that stops at the first sample below it.
    """

    if verbose:
        print("threshold for signal detection:", threshold)
    # Smooth the data
    smoothed = scipy_signal.savgol_filter(data, smooth_window, 3)
    return detect_active_window(time, smoothed, threshold=threshold)

def signal_onoffset(time, data, smooth_window=5, threshold=0.01, verbose=False):
    """Deprecated compatibility wrapper for :func:`signal_on_offset`."""
    warnings.warn(
        "signal_onoffset() is deprecated; use signal_on_offset().",
        DeprecationWarning,
        stacklevel=2,
    )
    return signal_on_offset(
        time,
        data,
        smooth_window=smooth_window,
        threshold=threshold,
        verbose=verbose,
    )

def smooth(array, span: int) -> np.ndarray:
    """Apply MATLAB-like moving-average smoothing with edge tapering.

    A centred moving average of odd width ``span``.  At each end the window
    shrinks symmetrically -- widths 1, 3, 5, ... -- so the output has the
    input's length and no sample is padded, which is what MATLAB's ``smooth``
    does and what the ported VEST workflows expect.

    Parameters
    ----------
    array : array_like
        The waveform, one-dimensional [any].
    span : int
        Window width in samples; an even value is reduced by one, and a value
        wider than the record is clipped to it [-].

    Returns
    -------
    np.ndarray
        The smoothed waveform, same length as ``array`` [any].

    Raises
    ------
    ValueError
        ``array`` is not one-dimensional.

    Convention
    ----------
    Zero-phase: the window is centred, so features are not shifted in time.
    Edge samples are averaged over progressively narrower windows rather than
    over a padded record, which biases the ends less than reflection padding
    would but leaves them noisier than the interior.

    Applicability
    -------------
    Machine-independent.

    Provenance
    ----------
    .. [1] MATLAB ``smooth(y, span)`` with the default ``'moving'`` method;
       ported from the legacy VEST workflow [VFIT]_ with its semantics
       preserved so that results match the donor sample for sample.
    """
    values = np.asarray(array, dtype=float)
    if values.ndim != 1:
        raise ValueError("`array` must be one-dimensional.")
    if values.size == 0 or span <= 1:
        return values.copy()
    if span % 2 == 0:
        span -= 1
    span = min(span, values.size if values.size % 2 == 1 else values.size - 1)
    if span < 1:
        return values.copy()

    count = values.size
    out = np.zeros(count, dtype=float)
    half = (span - 1) // 2

    # O(N) vectorized implementation using cumulative sums.
    cumsum = np.zeros(count + 1, dtype=float)
    cumsum[1:] = np.cumsum(values)

    # Interior: full-width windows of size `span`.
    interior = np.arange(half, count - half)
    if interior.size:
        out[interior] = (cumsum[interior + half + 1] - cumsum[interior - half]) / span

    # Edge tapering: left and right sides use progressively narrower windows.
    if half > 0:
        left_idx = np.arange(half)
        widths = 2 * left_idx + 1
        out[left_idx] = (cumsum[widths] - cumsum[0]) / widths
        right_idx = count - 1 - left_idx
        out[right_idx] = (cumsum[count] - cumsum[count - widths]) / widths

    return out

def vest_coil_current_noise_reduction(data) -> np.ndarray:
    """Suppress point spikes in coil current traces.

    Walks the record and, wherever a sample's magnitude exceeds the mean of
    its two neighbours' magnitudes by more than a fixed margin, replaces it
    with the previous sample.  A port of the donor routine, kept sample for
    sample so that the legacy PF-current preprocessing reproduces.

    Parameters
    ----------
    data : array_like
        Coil current, one-dimensional; records shorter than three samples are
        returned unchanged [A].

    Returns
    -------
    np.ndarray
        The de-spiked current, same length [A].

    Processing steps
    ----------------
    1. For each interior sample (the first two and the last are never
       touched), compute ``|x[i]| - |(x[i+1] + x[i-1]) / 2|``.
    2. Where that exceeds ``0.001``, overwrite ``x[i]`` with the already
       processed ``x[i-1]``.

    Defaults
    --------
    The ``0.001`` margin is a legacy compatibility value from the donor
    ``VEST_CoilCurrentNoiseReduction.m`` [1]_.  It is absolute, in the units of
    ``data``, and its derivation is not recorded.  On a trace in amperes it
    flags any excursion above a milliamp, so in practice every spike; on a
    trace in kiloamperes it would let real noise through.

    Convention
    ----------
    Causal: a spike is replaced by its predecessor, so a genuine step is
    delayed by one sample at its leading edge.

    Applicability
    -------------
    VEST-specific.  Written for the VEST PF coil current traces as digitized by
    the legacy DAQ; the margin has no meaning for another machine's units or
    noise floor.

    Limitations
    -----------
    Only a single-sample spike is caught: two consecutive bad samples pass,
    because the second is compared against the first.  Sign is ignored in the
    test but preserved in the replacement.

    Provenance
    ----------
    .. [1] ``Function/VEST_CoilCurrentNoiseReduction.m`` in the legacy VEST
       workflow [VFIT]_; loop bounds and margin preserved.
    """
    values = np.asarray(data, dtype=float)
    if values.size < 3:
        return values.copy()

    smoothed = values.copy()
    for index_i in range(2, values.size - 1):
        local_ref = abs((values[index_i + 1] + values[index_i - 1]) / 2.0)
        diff = abs(values[index_i]) - local_ref
        if diff > 0.001:
            smoothed[index_i] = smoothed[index_i - 1]
    return smoothed

def vfit_signal_start_end(time, data, threshold: float = 0.01) -> tuple[float, float]:
    """Deprecated compatibility wrapper for :func:`detect_active_window`."""
    warnings.warn(
        "vfit_signal_start_end() is deprecated; use detect_active_window().",
        DeprecationWarning,
        stacklevel=2,
    )
    return detect_active_window(time, data, threshold=threshold)

def vfit_signal_startend(time, data, threshold: float = 0.01):
    """Deprecated compatibility wrapper for :func:`detect_active_window`."""
    warnings.warn(
        "vfit_signal_startend() is deprecated; use detect_active_window().",
        DeprecationWarning,
        stacklevel=2,
    )
    return detect_active_window(time, data, threshold=threshold)