sgnts.transforms.resampler

Resampler dataclass

Bases: TSTransform

              flowchart TD
              sgnts.transforms.resampler.Resampler[Resampler]
              sgnts.base.base.TSTransform[TSTransform]
              sgnts.base.base.TimeSeriesMixin[TimeSeriesMixin]

                              sgnts.base.base.TSTransform --> sgnts.transforms.resampler.Resampler
                                sgnts.base.base.TimeSeriesMixin --> sgnts.base.base.TSTransform
                



              click sgnts.transforms.resampler.Resampler href "" "sgnts.transforms.resampler.Resampler"
              click sgnts.base.base.TSTransform href "" "sgnts.base.base.TSTransform"
              click sgnts.base.base.TimeSeriesMixin href "" "sgnts.base.base.TimeSeriesMixin"
            

Up/down samples time-series data

Parameters:

Name	Type	Description	Default
inrate	int	int, sample rate of the input frames	required
outrate	int	int, sample rate of the output frames	required
gstlal_norm	bool	boolean: If true it will normalize consistent with SGNL filter matching. If false it have a slightly more accurate normalization	True
use_gstlal_cpu_upsample	bool	boolean: If true, use the fast C-based gstlal implementation (sgnl_cpu_interp) for upsampling only; raises ImportError at configure time if the package is not installed	False
use_strided_downsample	bool | None	boolean: If None (the default), decide automatically whether scipy fft or strided correlation is used. Set to True/False do enable/ disable usage of this method for all kernel sizes.	None
up_half_length	int	int, half length (in input-rate samples) of the upsampling (anti-imaging) kernel. Defaults to UP_HALF_LENGTH (8). The full kernel is 2 * up_half_length * factor + 1 taps. Increase it when upsampling data with strong narrow-band content: a longer kernel deepens the stopband and suppresses the spectral images that otherwise appear at k * inrate +/- f. Only used on upsample paths (outrate > inrate); ignored when downsampling.	UP_HALF_LENGTH
use_simd_resample	bool	bool: opt in to dispatching numpy-backend resampling (both directions) to the sgnl_cpu_interp SIMD C extension. False (the default) preserves the historical behavior exactly: the pure numpy/scipy path, bit-for-bit, with the extension used only where the legacy use_gstlal_cpu_upsample flag engages it. True requires the package (ImportError at configure time if missing) and is mutually exclusive with use_gstlal_cpu_upsample (ValueError if both are set). The SIMD kernels are the same double-precision Lanczos-windowed sinc as the numpy path, so outputs agree at the ~1e-14 level but are not bit-identical. Inputs whose dtype the extension does not handle natively (anything other than float32/float64/complex64/complex128), or buffers shorter than the kernel, silently fall back to the numpy path. The torch backend is unaffected.	False

Notes

Thread safety: Marked thread_safe = True. With Pipeline.run(threaded=N) the pad callbacks for this element are dispatched onto worker threads.

Pad layout: 1 sink + 1 source pad
(enforced by ``@validator.one_to_one``). No same-element
``pull``/``new`` concurrency. ``internal`` runs alone.

Where the GIL-releasing work lives: ``internal()`` →
``process()`` → ``self.resample()``, which is one of
``scipy.signal.correlate``, ``torch.nn.functional.conv1d``,
or the C extension ``sgnl_cpu_interp.upsample_transposed``
— all release the GIL. Significant speedup expected with
multiple parallel resampling branches.

State touched per call:

- ``pull`` (inherited): per-pad-keyed dict writes.
- ``new`` (inherited): read-only lookup in ``self.outframes``.
- ``process``: reads ``self.thiskernel``, ``self.half_length``,
  ``self.resample`` (the bound method), and the
  ``adapter_config`` — all set in ``configure()`` and
  read-only afterwards. Writes only the local output buffer.

**Future editors MUST preserve thread safety**: kernel and
adapter state must remain post-init read-only. Do not
relax the one-to-one constraint without re-auditing.

Source code in src/sgnts/transforms/resampler.py

@dataclass(kw_only=True)
class Resampler(TSTransform):
    """Up/down samples time-series data

    Args:
        inrate:
            int, sample rate of the input frames
        outrate:
            int, sample rate of the output frames
        gstlal_norm:
            boolean: If true it will normalize consistent with SGNL
            filter matching. If false it have a slightly more accurate normalization
        use_gstlal_cpu_upsample:
            boolean: If true, use the fast C-based gstlal implementation
            (sgnl_cpu_interp) for upsampling only; raises ImportError at
            configure time if the package is not installed
        use_strided_downsample:
            boolean: If None (the default), decide automatically whether scipy
            fft or strided correlation is used. Set to True/False do enable/
            disable usage of this method for all kernel sizes.
        up_half_length:
            int, half length (in input-rate samples) of the upsampling
            (anti-imaging) kernel. Defaults to ``UP_HALF_LENGTH`` (8). The full
            kernel is ``2 * up_half_length * factor + 1`` taps. Increase it when
            upsampling data with strong narrow-band content: a
            longer kernel deepens the stopband and suppresses the
            spectral images that otherwise appear at ``k * inrate +/- f``. Only
            used on upsample paths (outrate > inrate); ignored when downsampling.
        use_simd_resample:
            bool: opt in to dispatching numpy-backend resampling (both
            directions) to the ``sgnl_cpu_interp`` SIMD C extension.
            ``False`` (the default) preserves the historical behavior
            exactly: the pure numpy/scipy path, bit-for-bit, with the
            extension used only where the legacy ``use_gstlal_cpu_upsample``
            flag engages it. ``True`` requires the package (ImportError at
            configure time if missing) and is mutually exclusive with
            ``use_gstlal_cpu_upsample`` (ValueError if both are set). The
            SIMD kernels are the same double-precision Lanczos-windowed sinc
            as the numpy path, so outputs agree at the ~1e-14 level but are
            not bit-identical. Inputs whose dtype the extension does not
            handle natively (anything other than
            float32/float64/complex64/complex128), or buffers shorter than
            the kernel, silently fall back to the numpy path. The torch
            backend is unaffected.

    Notes:
        Thread safety:
            Marked ``thread_safe = True``. With
            ``Pipeline.run(threaded=N)`` the pad callbacks for this
            element are dispatched onto worker threads.

            Pad layout: 1 sink + 1 source pad
            (enforced by ``@validator.one_to_one``). No same-element
            ``pull``/``new`` concurrency. ``internal`` runs alone.

            Where the GIL-releasing work lives: ``internal()`` →
            ``process()`` → ``self.resample()``, which is one of
            ``scipy.signal.correlate``, ``torch.nn.functional.conv1d``,
            or the C extension ``sgnl_cpu_interp.upsample_transposed``
            — all release the GIL. Significant speedup expected with
            multiple parallel resampling branches.

            State touched per call:

            - ``pull`` (inherited): per-pad-keyed dict writes.
            - ``new`` (inherited): read-only lookup in ``self.outframes``.
            - ``process``: reads ``self.thiskernel``, ``self.half_length``,
              ``self.resample`` (the bound method), and the
              ``adapter_config`` — all set in ``configure()`` and
              read-only afterwards. Writes only the local output buffer.

            **Future editors MUST preserve thread safety**: kernel and
            adapter state must remain post-init read-only. Do not
            relax the one-to-one constraint without re-auditing.
    """

    thread_safe = True

    # scipy.signal.correlate (numpy) and torch.nn.functional.conv1d (torch).
    backends = frozenset({"numpy", "torch"})

    inrate: int
    outrate: int
    gstlal_norm: bool = True
    use_gstlal_cpu_upsample: bool = False
    use_strided_downsample: bool | None = None
    up_half_length: int = UP_HALF_LENGTH
    use_simd_resample: bool = False

    def configure(self) -> None:
        self.next_out_offset = None

        # Correction from sgnl_cpu_interp's downsample normalization (unit DC
        # gain, i.e. the gstlal_norm=False convention) to this element's;
        # downkernel() overwrites it for gstlal_norm=True downsampling.
        self._simd_down_scale = 1.0
        if self.use_simd_resample:
            if self.use_gstlal_cpu_upsample:
                msg = (
                    "use_simd_resample and use_gstlal_cpu_upsample are "
                    "mutually exclusive: use_simd_resample covers both "
                    "directions; the legacy flag selects the upsample-only "
                    "path inside resample_numpy"
                )
                raise ValueError(msg)
            if not GSTLAL_AVAILABLE:
                msg = "use_simd_resample=True requires the sgnl_cpu_interp package"
                raise ImportError(msg)
            self._simd_enabled = True
        else:
            if self.use_gstlal_cpu_upsample and not GSTLAL_AVAILABLE:
                # Historically this fell back to the numpy path silently;
                # fail loudly instead so a requested acceleration (and its
                # float32 kernel numerics) can't quietly disappear when the
                # package is absent.
                msg = (
                    "use_gstlal_cpu_upsample=True requires the "
                    "sgnl_cpu_interp package"
                )
                raise ImportError(msg)
            self._simd_enabled = False

        if self.outrate < self.inrate:
            # downsample parameters
            factor = self.inrate // self.outrate
            self.half_length = int(DOWN_HALF_LENGTH * factor)
            self.kernel_length = self.half_length * 2 + 1
            self.thiskernel = self.downkernel(factor)
        elif self.outrate > self.inrate:
            # upsample parameters
            factor = self.outrate // self.inrate
            self.half_length = self.up_half_length
            self.kernel_length = self.half_length * 2 + 1
            self.thiskernel = self.upkernel(factor)
        else:
            # same rate
            raise ValueError("Inrate {self.inrate} is the same as outrate {outrate}")

        # ``self.thiskernel`` stays the numpy base kernel; its torch form is
        # built lazily per (dtype, device) from the data itself (see
        # ``_torch_kernel``). ``resample`` dispatches on the data's backend, so
        # this element adapts to whatever flows in. (A torch stream without torch
        # installed fails in ``resample_torch`` with a clear ImportError.)
        self._torch_kernels: dict = {}
        # Promoted-dtype forms of the numpy kernel, built lazily per data dtype
        # (see ``_numpy_kernel``).
        self._numpy_kernels: dict = {}

        self.adapter_config.overlap = (
            Offset.fromsamples(self.half_length, self.inrate),
            Offset.fromsamples(self.half_length, self.inrate),
        )
        self.adapter_config.on_startup(pad_zeros=True)

        self.pad_length = self.half_length

    @validator.one_to_one
    def validate(self) -> None:
        assert (
            self.inrate in Offset.ALLOWED_RATES
        ), f"Input rate {self.inrate} not in ALLOWED_RATES: {Offset.ALLOWED_RATES}"
        assert (
            self.outrate in Offset.ALLOWED_RATES
        ), f"Output rate {self.outrate} not in ALLOWED_RATES: {Offset.ALLOWED_RATES}"

    def downkernel(self, factor: int) -> Array:
        """Compute the kernel for downsampling. Modified from gstlal_interpolator.c

        This is a sinc windowed sinc function kernel
        The baseline kernel is defined as

        g[k] = sin(pi / f * (k-c)) / (pi / f * (k-c)) * (1 - (k-c)^2 / c / c)   k != c
        g[k] = 1                                                                k = c

        Where:

            f: downsample factor, must be power of 2, e.g., 2, 4, 8, ...
            c: defined as half the full kernel length

        You specify the half filter length at the target rate in samples,
        the kernel length is then given by:

            kernel_length = half_length_at_original_rate * 2 * f + 1


        Args:
            factor:
                int, factor = inrate/outrate

        Returns:
            Array, the downsampling kernel
        """
        kernel_length = int(2 * self.half_length + 1)

        # the domain should be the kernel_length divided by two
        c = kernel_length // 2
        x = np.arange(-c, c + 1)
        vecs = np.sinc(x / factor) * np.sinc(x / c)
        if self.gstlal_norm:
            norm = np.linalg.norm(vecs) * factor**0.5
            # sgnl_cpu_interp normalizes its downsample kernel to unit DC
            # gain (the gstlal_norm=False convention); rescale its output to
            # this element's convention.
            self._simd_down_scale = float(sum(vecs) / norm)
        else:
            norm = sum(vecs)
        vecs = vecs / norm
        return vecs.reshape(1, -1)

    def upkernel(self, factor: int) -> Array:
        """Compute the kernel for upsampling. Modified from gstlal_interpolator.c

        This is a sinc windowed sinc function kernel
        The baseline kernel is defined as

        $$\\begin{align}
        g(k) &= \\sin(\\pi / f * (k-c)) /
                (\\pi / f * (k-c)) * (1 - (k-c)^2 / c / c)  & k != c \\\\
        g(k) &= 1 & k = c
        \\end{align}$$

        Where `f` is the interpolation factor (must be power of 2, e.g., 2, 4, 8, ...)
        and `c` is defined as half the full kernel length.

        You specify the half filter length at the original rate in samples,
        the kernel length is then given by:

            kernel_length = half_length_at_original_rate * 2 * f + 1

        Interpolation is then defined as a two step process.  First the
        input data is zero filled to bring it up to the new sample rate,
        i.e., the input data, x, is transformed to x' such that:

        x'[i] = x[i/f]	if (i%f) == 0
              = 0       if (i%f) > 0

        y[i] = sum_{k=0}^{2c+1} x'[i-k] g[k]

        Since more than half the terms in this series would be zero, the
        convolution is implemented by breaking up the kernel into f separate
        kernels each 1/f as large as the originalcalled z, i.e.,:

        z[0][k/f] = g[k*f]
        z[1][k/f] = g[k*f+1]
        ...
        z[f-1][k/f] = g[k*f + f-1]

        Now the convolution can be written as:

        y[i] = sum_{k=0}^{2c/f+1} x[i/f] z[i%f][k]

        which avoids multiplying zeros.  Note also that by construction the
        sinc function has its zeros arranged such that z[0][:] had only one
        nonzero sample at its center. Therefore the actual convolution is:

        y[i] = x[i/f]					if i%f == 0
        y[i] = sum_{k=0}^{2c/f+1} x[i/f] z[i%f][k]	otherwise


        Args:
            factor:
                int, factor = outrate/inrate

        Returns:
            Array, the upsampling kernel
        """
        kernel_length = int(2 * self.half_length * factor + 1)
        sub_kernel_length = int(2 * self.half_length + 1)

        # the domain should be the kernel_length divided by two
        c = kernel_length // 2
        x = np.arange(-c, c + 1)
        out = np.sinc(x / factor) * np.sinc(x / c)
        out = np.pad(out, (0, factor - 1))
        # FIXME: check if interleave same as no interleave
        vecs = out.reshape(-1, factor).T[:, ::-1]

        return vecs.reshape(int(factor), 1, sub_kernel_length)

    def upsample_gstlal(self, data):
        """Upsample using gstlal implementation.

        Handles both numpy arrays and torch tensors.

        Args:
            data: Input data (numpy array or torch tensor), shape (-1, n_samples)

        Returns:
            Upsampled data (same type as input), not reshaped
        """
        # Check if input is torch tensor
        is_torch = TORCH_AVAILABLE and torch.is_tensor(data)
        if is_torch:
            # Convert torch -> numpy
            torch_device = data.device
            torch_dtype = data.dtype
            data_np = data.cpu().numpy()
        else:
            data_np = data

        # Call gstlal
        factor = self.outrate // self.inrate
        out_np = simd_upsample_transposed(
            data_np, factor=factor, half_length=self.half_length
        )

        # Convert back to torch if needed
        if is_torch:
            out = torch.from_numpy(out_np).to(torch_device).to(torch_dtype)
        else:
            out = out_np

        return out

    def resample(self, data: Array, outshape: tuple[int, ...]) -> Array:
        """Resample ``data``, dispatching on the backend of the data itself.

        The numpy/scipy and torch/conv1d kernels are the sanctioned
        backend-specific "escape hatch"; which one runs is decided by the array,
        not by a pre-declared backend.

        Args:
            data:
                Array, the data to be up/downsampled
            outshape:
                tuple[int, ...], the shape of the output array

        Returns:
            Array, the resampled data (same backend as ``data``)
        """
        if backend_name(data) == "torch":
            return self.resample_torch(data, outshape)
        if self._simd_enabled:
            out = self.resample_simd(data, outshape)
            if out is not None:
                return out
        return self.resample_numpy(data, outshape)

    def resample_simd(self, data0: Array, outshape: tuple[int, ...]) -> Array | None:
        """Resample via the ``sgnl_cpu_interp`` SIMD C extension.

        Same kernels and output-length contract as ``resample_numpy`` (the
        extension generates the identical double-precision Lanczos-windowed
        sinc; downsample output is rescaled to this element's normalization
        convention via ``self._simd_down_scale``).

        Returns None when the input cannot be dispatched -- dtype outside the
        extension's native set, or a buffer shorter than the kernel (where
        ``resample_numpy`` defines the empty-output contract) -- in which
        case the caller falls back to ``resample_numpy``.

        Args:
            data0:
                Array, the data to be up/downsampled
            outshape:
                tuple[int, ...], the shape of the output array

        Returns:
            Array | None, the resampled data, or None to request fallback
        """
        data = data0.reshape(-1, data0.shape[-1])
        if data.dtype not in _SIMD_NATIVE_DTYPES:
            return None

        if self.outrate > self.inrate:
            if data.shape[-1] < 2 * self.half_length + 1:
                return None
            out = simd_upsample_transposed(
                data,
                factor=self.outrate // self.inrate,
                half_length=self.half_length,
            )
        else:
            factor = self.inrate // self.outrate
            if data.shape[-1] < self.kernel_length:
                return None
            out = simd_downsample_transposed(
                data,
                factor=factor,
                half_length=self.half_length // factor,
            )
            if self._simd_down_scale != 1.0:
                out = out * self._simd_down_scale
        return out.reshape(outshape)

    def _torch_kernel(self, data: Array) -> Array:
        """Return the conv kernel as a tensor matching ``data``'s dtype/device.

        Derived from the numpy base kernel (``self.thiskernel``) and cached per
        (dtype, device). The kernel follows the data (no global dtype/device), so
        mixed precision and non-CPU devices work without reconfiguration.

        Args:
            data:
                Array, the incoming data whose dtype/device the kernel
                must match

        Returns:
            Array, the conv kernel
        """
        key = (data.dtype, data.device)
        kernel = self._torch_kernels.get(key)
        if kernel is None:
            base = self.thiskernel
            if self.outrate < self.inrate:
                # downsample: numpy (1, N) -> torch (1, 1, N)
                kernel = torch.from_numpy(base).view(1, 1, -1)
            else:
                # upsample: numpy (factor, 1, sub) -> torch (factor, 1, sub)
                kernel = torch.tensor(base.copy()).view(
                    self.outrate // self.inrate, 1, -1
                )
            kernel = kernel.to(device=data.device, dtype=data.dtype)
            self._torch_kernels[key] = kernel
        return kernel

    def _numpy_kernel(self, data: Array) -> Array:
        """Return the base kernel pre-cast to the promoted correlation dtype,
        cached per data dtype.

        numpy promotes the (real) base kernel against complex input on every
        matmul/einsum/correlate call, re-allocating and casting the kernel
        each time; caching the cast form makes that a one-time cost. The
        arithmetic is unchanged (same promotion, identical values), and for
        real data the base kernel is returned as-is.

        Args:
            data:
                Array, the incoming data whose dtype the kernel must promote
                against

        Returns:
            Array, the promoted kernel
        """
        dtype = np.result_type(data.dtype, self.thiskernel.dtype)
        kernel = self._numpy_kernels.get(dtype)
        if kernel is None:
            kernel = self.thiskernel.astype(dtype, copy=False)
            self._numpy_kernels[dtype] = kernel
        return kernel

    def resample_numpy(self, data0: Array, outshape: tuple[int, ...]) -> Array:
        """Correlate the data with the kernel.

        Args:
            data0:
                Array, the data to be up/downsampled
            outshape:
                tuple[int, ...], the shape of the output array

        Returns:
            Array, the resulting array of the up/downsamping
        """
        data = data0.reshape(-1, data0.shape[-1])
        thiskernel = self._numpy_kernel(data)

        if self.outrate > self.inrate:
            # upsample
            factor = self.outrate // self.inrate
            sub_len = thiskernel.shape[-1]
            if self.use_gstlal_cpu_upsample and GSTLAL_AVAILABLE:
                # Use fast C-based gstlal implementation
                out = self.upsample_gstlal(data)
            elif data.shape[-1] < sub_len:
                # Buffer shorter than a sub-kernel: no valid output window, so
                # the result is empty, in the promoted dtype the vectorized
                # contraction below would produce. (scipy's "valid" correlate
                # cannot express this: with in2 longer than in1 it swaps its
                # inputs and yields spurious samples, or raises for
                # multi-channel data.)
                out = np.empty((data.shape[0], 0), dtype=thiskernel.dtype)
            else:
                # Vectorized polyphase upsample: apply all ``factor`` sub-kernels
                # in one batched sliding-window contraction instead of ``factor``
                # separate scipy.correlate calls (a Python loop that dominated the
                # cost for large factors).
                #
                # ``corr`` is (nchan, factor, nwin) indexed [channel, phase,
                # time]. Each channel is upsampled independently, and its output
                # stream interleaves the phases within a time step: output sample
                # ``j*factor + i`` is phase ``i`` of input time ``j``. So the
                # per-channel layout is [time, phase] collapsed, i.e. transpose to
                # [channel, time, phase] then flatten the last two axes -- keeping
                # channels separate. (The previous ``(nwin, factor*nchan)`` layout
                # interleaved channels into the phase axis, which scrambled
                # multi-channel output; single-channel was coincidentally
                # unaffected. This layout matches ``resample_torch``.)
                kernels = thiskernel.reshape(factor, sub_len)
                windows = sliding_window_view(data, sub_len, axis=-1)
                corr = np.einsum("cjk,ik->cij", windows, kernels)
                out = corr.transpose(0, 2, 1).reshape(data.shape[0], -1)
        else:
            # downsample
            factor = self.inrate // self.outrate

            use_strided = self.use_strided_downsample
            if use_strided is None:
                use_strided = self._prefer_strided(data.shape[-1])

            if data.shape[-1] < self.kernel_length:
                # Buffer shorter than the kernel: no valid output window (the
                # SIMD path punts these here for the authoritative empty
                # shape; scipy's "valid" correlate would swap its inputs or
                # raise instead of yielding empty).
                out = np.empty((data.shape[0], 0), dtype=thiskernel.dtype)
            elif use_strided:
                # Strided correlation: compute only the kept outputs. For large
                # decimation factors the kernel (``64*factor+1`` taps) dwarfs
                # the buffer, so scipy's FFT path spends almost all its work
                # computing samples the ``::factor`` slice then discards. A
                # strided sliding-window dot product costs ``n_out * L`` MACs
                # (multiply-accumulate; do a * b and add this to a running sum)
                # regardless of factor and wins once the FFT block is
                # kernel-dominated. Bit-comparable to the FFT path.
                kernel_1d = thiskernel.reshape(-1)
                windows = sliding_window_view(data, kernel_1d.size, axis=-1)[
                    :, ::factor, :
                ]
                out = windows @ kernel_1d
            else:
                out = correlate(data, thiskernel, mode="valid")[..., ::factor]
        return out.reshape(outshape)

    def _prefer_strided(self, n_in: int) -> bool:
        """Whether strided correlation beats scipy's FFT path for this buffer.

        ``scipy.signal.correlate(method="auto")`` computes *every* output, then
        the ``::factor`` slice throws away all but every ``factor``-th. When the
        kernel is long relative to the buffer (large decimation factor), that
        FFT is dominated by the kernel and strided direct correlation -- which
        computes only the kept outputs -- is cheaper. We compare a MAC count
        (``n_out * L``, roughly constant at ``64 * buffer``) against an FFT
        butterfly estimate (``(n_in + L) * log2(n_in + L)``); empirically this
        switches over around factor 512 for strain-rate buffers. Only downsample
        paths call this (upsample uses a different kernel layout).

        Args:
            n_in:
                int, the number of input samples along the last axis

        Returns:
            bool, True to use the strided path, False to keep scipy's FFT path
        """
        factor = self.inrate // self.outrate
        length = self.kernel_length
        n_out = (n_in - length) // factor + 1
        if n_out <= 0:
            # buffer shorter than the kernel: no valid output either way, let
            # scipy handle the (empty) result to keep behaviour identical.
            return False
        return n_out * length < (n_in + length) * math.log2(n_in + length)

    def resample_torch(self, data0: Array, outshape: tuple[int, ...]) -> Array:
        """Correlate the data with the kernel.

        Args:
            data0:
                Array, the data to be up/downsampled
            outshape:
                tuple[int, ...], the shape of the output array

        Returns:
            Array, the resulting array of the up/downsamping
        """
        if not TORCH_AVAILABLE:
            raise ImportError(
                "PyTorch is not installed. Install it with 'pip install sgn-ts[torch]'"
            )

        if self.outrate > self.inrate:  # upsample
            if self.use_gstlal_cpu_upsample and GSTLAL_AVAILABLE:
                # Use gstlal (handles torch->numpy->torch conversion)
                data = data0.view(-1, data0.shape[-1])
                out = self.upsample_gstlal(data)
                return out.view(outshape)
            else:
                # Use PyTorch conv1d
                data = data0.view(-1, 1, data0.shape[-1])
                # The kernel is built to match the data's dtype/device.
                thiskernel = self._torch_kernel(data0)
                out = Fconv1d(data, thiskernel)
                out = out.mT.reshape(data.shape[0], -1)
                return out.view(outshape)
        else:  # downsample
            data = data0.view(-1, 1, data0.shape[-1])
            # The kernel is built to match the data's dtype/device.
            thiskernel = self._torch_kernel(data0)
            out = Fconv1d(data, thiskernel, stride=self.inrate // self.outrate)
            out = out.squeeze(1)

        return out.view(outshape)

    @transform.one_to_one
    def process(self, input_frame: TSFrame, output_frame: TSCollectFrame) -> None:
        """Resample input frame to output sample rate."""
        assert input_frame.sample_rate == self.inrate, (
            f"Frame sample rate {input_frame.sample_rate} doesn't match "
            f"resampler input rate {self.inrate}"
        )

        if input_frame.shape[-1] == 0:
            buf = SeriesBuffer(
                offset=output_frame.offset,
                sample_rate=self.outrate,
                data=None,
                shape=input_frame.shape,
            )
            output_frame.append(buf)
        else:
            for buf in input_frame:
                shape = input_frame.shape[:-1] + (
                    Offset.tosamples(output_frame.noffset, self.outrate),
                )
                if buf.is_gap:
                    data = None
                else:
                    assert buf.data is not None and not isinstance(buf.data, int)
                    data = self.resample(buf.data, shape)
                buf = buf.copy(
                    offset=output_frame.offset,
                    sample_rate=self.outrate,
                    data=data,
                    shape=shape,
                )
                output_frame.append(buf)

downkernel(factor)

Compute the kernel for downsampling. Modified from gstlal_interpolator.c

This is a sinc windowed sinc function kernel The baseline kernel is defined as

g[k] = sin(pi / f * (k-c)) / (pi / f * (k-c)) * (1 - (k-c)^2 / c / c) k != c g[k] = 1 k = c

Where:

f: downsample factor, must be power of 2, e.g., 2, 4, 8, ...
c: defined as half the full kernel length

You specify the half filter length at the target rate in samples, the kernel length is then given by:

kernel_length = half_length_at_original_rate * 2 * f + 1

Parameters:

Name	Type	Description	Default
factor	int	int, factor = inrate/outrate	required

Returns:

Type	Description
Array	Array, the downsampling kernel

Source code in src/sgnts/transforms/resampler.py

def downkernel(self, factor: int) -> Array:
    """Compute the kernel for downsampling. Modified from gstlal_interpolator.c

    This is a sinc windowed sinc function kernel
    The baseline kernel is defined as

    g[k] = sin(pi / f * (k-c)) / (pi / f * (k-c)) * (1 - (k-c)^2 / c / c)   k != c
    g[k] = 1                                                                k = c

    Where:

        f: downsample factor, must be power of 2, e.g., 2, 4, 8, ...
        c: defined as half the full kernel length

    You specify the half filter length at the target rate in samples,
    the kernel length is then given by:

        kernel_length = half_length_at_original_rate * 2 * f + 1


    Args:
        factor:
            int, factor = inrate/outrate

    Returns:
        Array, the downsampling kernel
    """
    kernel_length = int(2 * self.half_length + 1)

    # the domain should be the kernel_length divided by two
    c = kernel_length // 2
    x = np.arange(-c, c + 1)
    vecs = np.sinc(x / factor) * np.sinc(x / c)
    if self.gstlal_norm:
        norm = np.linalg.norm(vecs) * factor**0.5
        # sgnl_cpu_interp normalizes its downsample kernel to unit DC
        # gain (the gstlal_norm=False convention); rescale its output to
        # this element's convention.
        self._simd_down_scale = float(sum(vecs) / norm)
    else:
        norm = sum(vecs)
    vecs = vecs / norm
    return vecs.reshape(1, -1)

process(input_frame, output_frame)

Resample input frame to output sample rate.

Source code in src/sgnts/transforms/resampler.py

@transform.one_to_one
def process(self, input_frame: TSFrame, output_frame: TSCollectFrame) -> None:
    """Resample input frame to output sample rate."""
    assert input_frame.sample_rate == self.inrate, (
        f"Frame sample rate {input_frame.sample_rate} doesn't match "
        f"resampler input rate {self.inrate}"
    )

    if input_frame.shape[-1] == 0:
        buf = SeriesBuffer(
            offset=output_frame.offset,
            sample_rate=self.outrate,
            data=None,
            shape=input_frame.shape,
        )
        output_frame.append(buf)
    else:
        for buf in input_frame:
            shape = input_frame.shape[:-1] + (
                Offset.tosamples(output_frame.noffset, self.outrate),
            )
            if buf.is_gap:
                data = None
            else:
                assert buf.data is not None and not isinstance(buf.data, int)
                data = self.resample(buf.data, shape)
            buf = buf.copy(
                offset=output_frame.offset,
                sample_rate=self.outrate,
                data=data,
                shape=shape,
            )
            output_frame.append(buf)

resample(data, outshape)

Resample data, dispatching on the backend of the data itself.

The numpy/scipy and torch/conv1d kernels are the sanctioned backend-specific "escape hatch"; which one runs is decided by the array, not by a pre-declared backend.

Parameters:

Name	Type	Description	Default
data	Array	Array, the data to be up/downsampled	required
outshape	tuple[int, ...]	tuple[int, ...], the shape of the output array	required

Returns:

Type	Description
Array	Array, the resampled data (same backend as data)

Source code in src/sgnts/transforms/resampler.py

def resample(self, data: Array, outshape: tuple[int, ...]) -> Array:
    """Resample ``data``, dispatching on the backend of the data itself.

    The numpy/scipy and torch/conv1d kernels are the sanctioned
    backend-specific "escape hatch"; which one runs is decided by the array,
    not by a pre-declared backend.

    Args:
        data:
            Array, the data to be up/downsampled
        outshape:
            tuple[int, ...], the shape of the output array

    Returns:
        Array, the resampled data (same backend as ``data``)
    """
    if backend_name(data) == "torch":
        return self.resample_torch(data, outshape)
    if self._simd_enabled:
        out = self.resample_simd(data, outshape)
        if out is not None:
            return out
    return self.resample_numpy(data, outshape)

resample_numpy(data0, outshape)

Correlate the data with the kernel.

Parameters:

Name	Type	Description	Default
data0	Array	Array, the data to be up/downsampled	required
outshape	tuple[int, ...]	tuple[int, ...], the shape of the output array	required

Returns:

Type	Description
Array	Array, the resulting array of the up/downsamping

Source code in src/sgnts/transforms/resampler.py

def resample_numpy(self, data0: Array, outshape: tuple[int, ...]) -> Array:
    """Correlate the data with the kernel.

    Args:
        data0:
            Array, the data to be up/downsampled
        outshape:
            tuple[int, ...], the shape of the output array

    Returns:
        Array, the resulting array of the up/downsamping
    """
    data = data0.reshape(-1, data0.shape[-1])
    thiskernel = self._numpy_kernel(data)

    if self.outrate > self.inrate:
        # upsample
        factor = self.outrate // self.inrate
        sub_len = thiskernel.shape[-1]
        if self.use_gstlal_cpu_upsample and GSTLAL_AVAILABLE:
            # Use fast C-based gstlal implementation
            out = self.upsample_gstlal(data)
        elif data.shape[-1] < sub_len:
            # Buffer shorter than a sub-kernel: no valid output window, so
            # the result is empty, in the promoted dtype the vectorized
            # contraction below would produce. (scipy's "valid" correlate
            # cannot express this: with in2 longer than in1 it swaps its
            # inputs and yields spurious samples, or raises for
            # multi-channel data.)
            out = np.empty((data.shape[0], 0), dtype=thiskernel.dtype)
        else:
            # Vectorized polyphase upsample: apply all ``factor`` sub-kernels
            # in one batched sliding-window contraction instead of ``factor``
            # separate scipy.correlate calls (a Python loop that dominated the
            # cost for large factors).
            #
            # ``corr`` is (nchan, factor, nwin) indexed [channel, phase,
            # time]. Each channel is upsampled independently, and its output
            # stream interleaves the phases within a time step: output sample
            # ``j*factor + i`` is phase ``i`` of input time ``j``. So the
            # per-channel layout is [time, phase] collapsed, i.e. transpose to
            # [channel, time, phase] then flatten the last two axes -- keeping
            # channels separate. (The previous ``(nwin, factor*nchan)`` layout
            # interleaved channels into the phase axis, which scrambled
            # multi-channel output; single-channel was coincidentally
            # unaffected. This layout matches ``resample_torch``.)
            kernels = thiskernel.reshape(factor, sub_len)
            windows = sliding_window_view(data, sub_len, axis=-1)
            corr = np.einsum("cjk,ik->cij", windows, kernels)
            out = corr.transpose(0, 2, 1).reshape(data.shape[0], -1)
    else:
        # downsample
        factor = self.inrate // self.outrate

        use_strided = self.use_strided_downsample
        if use_strided is None:
            use_strided = self._prefer_strided(data.shape[-1])

        if data.shape[-1] < self.kernel_length:
            # Buffer shorter than the kernel: no valid output window (the
            # SIMD path punts these here for the authoritative empty
            # shape; scipy's "valid" correlate would swap its inputs or
            # raise instead of yielding empty).
            out = np.empty((data.shape[0], 0), dtype=thiskernel.dtype)
        elif use_strided:
            # Strided correlation: compute only the kept outputs. For large
            # decimation factors the kernel (``64*factor+1`` taps) dwarfs
            # the buffer, so scipy's FFT path spends almost all its work
            # computing samples the ``::factor`` slice then discards. A
            # strided sliding-window dot product costs ``n_out * L`` MACs
            # (multiply-accumulate; do a * b and add this to a running sum)
            # regardless of factor and wins once the FFT block is
            # kernel-dominated. Bit-comparable to the FFT path.
            kernel_1d = thiskernel.reshape(-1)
            windows = sliding_window_view(data, kernel_1d.size, axis=-1)[
                :, ::factor, :
            ]
            out = windows @ kernel_1d
        else:
            out = correlate(data, thiskernel, mode="valid")[..., ::factor]
    return out.reshape(outshape)

resample_simd(data0, outshape)

Resample via the sgnl_cpu_interp SIMD C extension.

Same kernels and output-length contract as resample_numpy (the extension generates the identical double-precision Lanczos-windowed sinc; downsample output is rescaled to this element's normalization convention via self._simd_down_scale).

Returns None when the input cannot be dispatched -- dtype outside the extension's native set, or a buffer shorter than the kernel (where resample_numpy defines the empty-output contract) -- in which case the caller falls back to resample_numpy.

Parameters:

Name	Type	Description	Default
data0	Array	Array, the data to be up/downsampled	required
outshape	tuple[int, ...]	tuple[int, ...], the shape of the output array	required

Returns:

Type	Description
Array | None	Array | None, the resampled data, or None to request fallback

Source code in src/sgnts/transforms/resampler.py

def resample_simd(self, data0: Array, outshape: tuple[int, ...]) -> Array | None:
    """Resample via the ``sgnl_cpu_interp`` SIMD C extension.

    Same kernels and output-length contract as ``resample_numpy`` (the
    extension generates the identical double-precision Lanczos-windowed
    sinc; downsample output is rescaled to this element's normalization
    convention via ``self._simd_down_scale``).

    Returns None when the input cannot be dispatched -- dtype outside the
    extension's native set, or a buffer shorter than the kernel (where
    ``resample_numpy`` defines the empty-output contract) -- in which
    case the caller falls back to ``resample_numpy``.

    Args:
        data0:
            Array, the data to be up/downsampled
        outshape:
            tuple[int, ...], the shape of the output array

    Returns:
        Array | None, the resampled data, or None to request fallback
    """
    data = data0.reshape(-1, data0.shape[-1])
    if data.dtype not in _SIMD_NATIVE_DTYPES:
        return None

    if self.outrate > self.inrate:
        if data.shape[-1] < 2 * self.half_length + 1:
            return None
        out = simd_upsample_transposed(
            data,
            factor=self.outrate // self.inrate,
            half_length=self.half_length,
        )
    else:
        factor = self.inrate // self.outrate
        if data.shape[-1] < self.kernel_length:
            return None
        out = simd_downsample_transposed(
            data,
            factor=factor,
            half_length=self.half_length // factor,
        )
        if self._simd_down_scale != 1.0:
            out = out * self._simd_down_scale
    return out.reshape(outshape)

resample_torch(data0, outshape)

Correlate the data with the kernel.

Parameters:

Name	Type	Description	Default
data0	Array	Array, the data to be up/downsampled	required
outshape	tuple[int, ...]	tuple[int, ...], the shape of the output array	required

Returns:

Type	Description
Array	Array, the resulting array of the up/downsamping

Source code in src/sgnts/transforms/resampler.py

def resample_torch(self, data0: Array, outshape: tuple[int, ...]) -> Array:
    """Correlate the data with the kernel.

    Args:
        data0:
            Array, the data to be up/downsampled
        outshape:
            tuple[int, ...], the shape of the output array

    Returns:
        Array, the resulting array of the up/downsamping
    """
    if not TORCH_AVAILABLE:
        raise ImportError(
            "PyTorch is not installed. Install it with 'pip install sgn-ts[torch]'"
        )

    if self.outrate > self.inrate:  # upsample
        if self.use_gstlal_cpu_upsample and GSTLAL_AVAILABLE:
            # Use gstlal (handles torch->numpy->torch conversion)
            data = data0.view(-1, data0.shape[-1])
            out = self.upsample_gstlal(data)
            return out.view(outshape)
        else:
            # Use PyTorch conv1d
            data = data0.view(-1, 1, data0.shape[-1])
            # The kernel is built to match the data's dtype/device.
            thiskernel = self._torch_kernel(data0)
            out = Fconv1d(data, thiskernel)
            out = out.mT.reshape(data.shape[0], -1)
            return out.view(outshape)
    else:  # downsample
        data = data0.view(-1, 1, data0.shape[-1])
        # The kernel is built to match the data's dtype/device.
        thiskernel = self._torch_kernel(data0)
        out = Fconv1d(data, thiskernel, stride=self.inrate // self.outrate)
        out = out.squeeze(1)

    return out.view(outshape)

upkernel(factor)

Compute the kernel for upsampling. Modified from gstlal_interpolator.c

This is a sinc windowed sinc function kernel The baseline kernel is defined as

\[\begin{align} g(k) &= \sin(\pi / f * (k-c)) / (\pi / f * (k-c)) * (1 - (k-c)^2 / c / c) & k != c \\ g(k) &= 1 & k = c \end{align}\]

Where f is the interpolation factor (must be power of 2, e.g., 2, 4, 8, ...) and c is defined as half the full kernel length.

You specify the half filter length at the original rate in samples, the kernel length is then given by:

kernel_length = half_length_at_original_rate * 2 * f + 1

Interpolation is then defined as a two step process. First the input data is zero filled to bring it up to the new sample rate, i.e., the input data, x, is transformed to x' such that:

x'[i] = x[i/f] if (i%f) == 0 = 0 if (i%f) > 0

y[i] = sum_{k=0}^{2c+1} x'[i-k] g[k]

Since more than half the terms in this series would be zero, the convolution is implemented by breaking up the kernel into f separate kernels each 1/f as large as the originalcalled z, i.e.,:

z[0][k/f] = g[kf] z[1][k/f] = g[kf+1] ... z[f-1][k/f] = g[k*f + f-1]

Now the convolution can be written as:

y[i] = sum_{k=0}^{2c/f+1} x[i/f] z[i%f][k]

which avoids multiplying zeros. Note also that by construction the sinc function has its zeros arranged such that z[0][:] had only one nonzero sample at its center. Therefore the actual convolution is:

y[i] = x[i/f] if i%f == 0 y[i] = sum_{k=0}^{2c/f+1} x[i/f] z[i%f][k] otherwise

Parameters:

Name	Type	Description	Default
factor	int	int, factor = outrate/inrate	required

Returns:

Type	Description
Array	Array, the upsampling kernel

Source code in src/sgnts/transforms/resampler.py

def upkernel(self, factor: int) -> Array:
    """Compute the kernel for upsampling. Modified from gstlal_interpolator.c

    This is a sinc windowed sinc function kernel
    The baseline kernel is defined as

    $$\\begin{align}
    g(k) &= \\sin(\\pi / f * (k-c)) /
            (\\pi / f * (k-c)) * (1 - (k-c)^2 / c / c)  & k != c \\\\
    g(k) &= 1 & k = c
    \\end{align}$$

    Where `f` is the interpolation factor (must be power of 2, e.g., 2, 4, 8, ...)
    and `c` is defined as half the full kernel length.

    You specify the half filter length at the original rate in samples,
    the kernel length is then given by:

        kernel_length = half_length_at_original_rate * 2 * f + 1

    Interpolation is then defined as a two step process.  First the
    input data is zero filled to bring it up to the new sample rate,
    i.e., the input data, x, is transformed to x' such that:

    x'[i] = x[i/f]	if (i%f) == 0
          = 0       if (i%f) > 0

    y[i] = sum_{k=0}^{2c+1} x'[i-k] g[k]

    Since more than half the terms in this series would be zero, the
    convolution is implemented by breaking up the kernel into f separate
    kernels each 1/f as large as the originalcalled z, i.e.,:

    z[0][k/f] = g[k*f]
    z[1][k/f] = g[k*f+1]
    ...
    z[f-1][k/f] = g[k*f + f-1]

    Now the convolution can be written as:

    y[i] = sum_{k=0}^{2c/f+1} x[i/f] z[i%f][k]

    which avoids multiplying zeros.  Note also that by construction the
    sinc function has its zeros arranged such that z[0][:] had only one
    nonzero sample at its center. Therefore the actual convolution is:

    y[i] = x[i/f]					if i%f == 0
    y[i] = sum_{k=0}^{2c/f+1} x[i/f] z[i%f][k]	otherwise


    Args:
        factor:
            int, factor = outrate/inrate

    Returns:
        Array, the upsampling kernel
    """
    kernel_length = int(2 * self.half_length * factor + 1)
    sub_kernel_length = int(2 * self.half_length + 1)

    # the domain should be the kernel_length divided by two
    c = kernel_length // 2
    x = np.arange(-c, c + 1)
    out = np.sinc(x / factor) * np.sinc(x / c)
    out = np.pad(out, (0, factor - 1))
    # FIXME: check if interleave same as no interleave
    vecs = out.reshape(-1, factor).T[:, ::-1]

    return vecs.reshape(int(factor), 1, sub_kernel_length)

upsample_gstlal(data)

Upsample using gstlal implementation.

Handles both numpy arrays and torch tensors.

Parameters:

Name	Type	Description	Default
data		Input data (numpy array or torch tensor), shape (-1, n_samples)	required

Returns:

Type	Description
	Upsampled data (same type as input), not reshaped

Source code in src/sgnts/transforms/resampler.py

def upsample_gstlal(self, data):
    """Upsample using gstlal implementation.

    Handles both numpy arrays and torch tensors.

    Args:
        data: Input data (numpy array or torch tensor), shape (-1, n_samples)

    Returns:
        Upsampled data (same type as input), not reshaped
    """
    # Check if input is torch tensor
    is_torch = TORCH_AVAILABLE and torch.is_tensor(data)
    if is_torch:
        # Convert torch -> numpy
        torch_device = data.device
        torch_dtype = data.dtype
        data_np = data.cpu().numpy()
    else:
        data_np = data

    # Call gstlal
    factor = self.outrate // self.inrate
    out_np = simd_upsample_transposed(
        data_np, factor=factor, half_length=self.half_length
    )

    # Convert back to torch if needed
    if is_torch:
        out = torch.from_numpy(out_np).to(torch_device).to(torch_dtype)
    else:
        out = out_np

    return out