Resampling kernel

def kaiser_window_chebyshev(beta: float, tolerance: float = _WINDOW_FIT_TOLERANCE) -> tuple[float, ...] | None:
    """Cache-normalising wrapper; see :func:`_kaiser_window_chebyshev` for the fit itself.

    ``float()`` before the cache, on a reviewer's finding: ``32`` and ``32.0`` are distinct keys to
    ``lru_cache`` but the same Kaiser window, so an int caller silently paid for a second fit and
    filled a slot. numpy scalars have the same problem.
    """
    return _kaiser_window_chebyshev(float(beta), float(tolerance))

def oversample_fourier(samples: np.ndarray, factor: int = OVERSAMPLING_FACTOR) -> np.ndarray:
    """Return the band-limited interpolant of *samples* on a grid *factor* times finer.

    Sample ``k`` of *samples* is sample ``factor * k`` of the result, which has ``factor * (n - 1) + 1`` samples,
    so a position ``p`` on the input's grid is the position ``factor * p`` on the result's.

    Computed in the frequency domain: the input, padded with zeros to twice its length, is transformed, and its
    spectrum extended with zeros up to the finer grid's Nyquist frequency, so nothing is added or removed below
    the input's own. A component at exactly the input's Nyquist frequency has no unique interpolant (its samples
    are those of a cosine at that frequency whatever its phase); it is split evenly between the positive and the
    negative frequency, which keeps the result real and reproduces the input's samples.

    What this computes is the band-limited interpolant of the padded series repeated with period ``2 n``. For a
    series that falls smoothly to zero at both ends, such as a whole generated waveform, the repetitions
    contribute nothing and it is the band-limited interpolant of the series taken as zero outside its samples,
    accurate to roundoff up to the input's Nyquist frequency. For one that ends abruptly it is neither: the
    abrupt ends ring across the whole result -- 8e-5 of the signal mid-series for an abruptly ending tone at
    0.995 x Nyquist of 8192 samples, and more towards its ends -- where the windowed-sinc kernel alone is affected only within half its taps of them. Use it only on
    whole, tapered series. LAL's injection (``XLALSimAddInjectionREAL8TimeSeries``) shifts a series by a fraction
    of a sample the same way, in the frequency domain over the whole zero-padded series.

    Peak memory is about 110 bytes per input sample at a factor of 2 (6.1 GB at 5.2e7 samples, a 1.8-hour series
    at 8192 Hz), most of it the transforms' own working arrays.

    Args:
        samples: Uniformly sampled real series, shape ``(n,)``; taken as float, as by the kernel.
        factor: How many times finer the result's grid is; a positive integer.

    Returns:
        The interpolated series, shape ``(factor * (n - 1) + 1,)``; empty for an empty input.

    Raises:
        ValueError: If ``factor`` is not a positive integer, or ``samples`` is not one-dimensional.
    """
    if isinstance(factor, bool) or not isinstance(factor, (int, np.integer)) or factor < 1:
        raise ValueError(f"factor must be a positive integer; got {factor!r}.")
    samples = np.asarray(samples, dtype=float)
    if samples.ndim != 1:
        raise ValueError(f"samples must be one-dimensional; got shape {samples.shape}.")
    n = samples.shape[0]
    if n == 0 or factor == 1:
        return samples.copy()
    padded = 2 * n
    spectrum = np.fft.rfft(samples, padded)
    spectrum[-1] *= 0.5
    # `irfft` extends the spectrum with zeros to the requested length itself, so the finer spectrum is
    # never built as a separate array.
    return np.fft.irfft(spectrum, factor * padded)[: factor * (n - 1) + 1] * factor

def resample_uniform_sinc(
    samples: np.ndarray,
    index: np.ndarray,
    *,
    fraction: np.ndarray | float | None = None,
    taps: int = DEFAULT_SINC_TAPS,
    beta: float = DEFAULT_KAISER_BETA,
) -> np.ndarray:
    """Resample a uniformly sampled series at fractional positions (NumPy).

    Evaluates the Kaiser-windowed sinc series. Positions outside
    ``[0, len(samples) - 1]`` return zero, matching the ``fill_value=0.0`` behaviour of
    the interpolation this replaces.

    Weights are normalised to sum to one, which enforces unit gain at DC and removes
    the amplitude bias the window would otherwise introduce near the ends of the series.

    A position given as one float64 keeps about 16 significant digits, so its fractional part is
    resolved only to ``index * 2.2e-16`` of a sample: about 1e-9 of a sample at ten million samples,
    which the kernel turns into an error of that much time times the signal's slope. A caller whose
    positions are a large whole number plus a small offset should pass them as such: ``index`` the
    integer part and ``fraction`` the offset. The two are never added in floating point, so the
    offset keeps its own precision whatever the size of ``index``.

    Args:
        samples: Uniformly sampled series, shape ``(n,)``.
        index: Sample positions to evaluate at, any shape: fractional, or, with ``fraction``, whole
            numbers (an integer array).
        fraction: Offset added to each position of ``index``, in samples; a scalar or an array
            broadcastable to ``index``. Any size, though its precision is that of its own magnitude.
        taps: Number of kernel taps.
        beta: Kaiser window shape parameter.

    Returns:
        Interpolated values, the shape of ``index``.

    Raises:
        TypeError: If ``fraction`` is given and ``index`` is not an integer array castable to ``int64``
            (``uint64`` is not: its values above ``2**63`` would wrap).
    """
    taps, beta = validate_kernel(taps, beta)
    samples = np.asarray(samples, dtype=float)
    if fraction is None:
        positions = np.asarray(index, dtype=float)
        whole = np.floor(positions)
        part = positions - whole
        whole = whole.astype(np.int64)
    else:
        whole = np.asarray(index)
        if not np.issubdtype(whole.dtype, np.integer) or not np.can_cast(whole.dtype, np.int64):
            raise TypeError(
                f"index must be an integer array castable to int64 when fraction is given, not {whole.dtype}"
            )
        offset = np.broadcast_to(np.asarray(fraction, dtype=float), whole.shape)
        whole = whole.astype(np.int64)
        carried = np.floor(np.where(np.isfinite(offset), offset, 0.0))
        # A whole sample or a carried fraction this far out puts the position outside any series, and so
        # does a non-finite fraction (as a non-finite float position): such a position places nothing.
        # Excluding them first means the sum below cannot overflow.
        far = ~np.isfinite(offset) | (whole > _FAR) | (whole < -_FAR) | (np.abs(carried) > _FAR)
        whole = np.where(far, -1, whole + np.where(far, 0.0, carried).astype(np.int64))
        part = np.where(far, 0.0, offset - carried)
    shape = whole.shape
    whole, part = whole.ravel(), part.ravel()
    result = np.empty(whole.shape, dtype=float)
    for first in range(0, whole.size, _BLOCK_POSITIONS):
        block = slice(first, first + _BLOCK_POSITIONS)
        result[block] = _resample_block(samples, whole[block], part[block], taps=taps, beta=beta)
    return result.reshape(shape)

def validate_kernel(taps: int, beta: float) -> tuple[int, float]:
    """Validate a kernel specification and return it normalised.

    Args:
        taps: Number of kernel taps; must be an odd integer >= 3.
        beta: Kaiser window shape parameter; must be finite, non-negative and small enough
            that its transition band fits in ``taps``.

    Returns:
        ``(taps, beta)`` as an ``int`` and a ``float``.

    Raises:
        ValueError: If ``taps`` is not an odd integer >= 3, if ``beta`` is not finite or is
            negative, or if ``beta`` is too large for ``taps``. The last case is rejected
            rather than warned about because it fails *quietly*: the kernel still returns
            plausible numbers, several orders of magnitude less accurate than the tap count
            implies.
    """
    # Rejected rather than truncated: int(127.9) would silently accept a value the
    # integer-only contract forbids, and the caller would never learn which kernel ran.
    if isinstance(taps, float) and not taps.is_integer():
        raise ValueError(f"taps must be an odd integer >= 3; got {taps}.")
    taps = int(taps)
    if taps < _MINIMUM_TAPS or taps % 2 == 0:
        raise ValueError(f"taps must be an odd integer >= 3; got {taps}.")
    beta = float(beta)
    # NaN fails every comparison, so it would slip past a bare `beta < 0` check and then
    # propagate silently through every interpolated sample.
    if not np.isfinite(beta):
        raise ValueError(f"Kaiser beta must be finite; got {beta}.")
    if beta < 0.0:
        raise ValueError(f"Kaiser beta must be non-negative; got {beta}.")
    minimum_taps = _TAPS_PER_BETA * beta - 1.0
    if taps < minimum_taps:
        raise ValueError(
            f"Kaiser beta={beta} needs at least {minimum_taps:.0f} taps to hold its transition "
            f"band, but taps={taps}. Either raise taps or lower beta; a beta too large for the "
            f"kernel silently loses several orders of magnitude of accuracy."
        )
    return taps, beta