Core

blipss.core.compare_cands

Cross-file comparison and clustering of FFA candidate periods.

Pipeline

1. Filter: For each input candidate file, retain only fundamental-flagged candidates whose S/N exceeds a pointing-specific threshold.
2. Merge: Concatenate filtered candidates across all files, tagged with a source file index.
3. Group: For each spectral channel, cluster candidate periods via Friends-of-Friends and keep the highest-S/N representative per cluster.
4. Encode: Build an N-file binary detection code for each surviving cluster, where the i-th digit is '1' if file i contributed a candidate to that cluster and '0' otherwise.

Results from all channels are merged and returned by group_candidates_by_channel.

filter_fundamental_candidates(channels, radiofreqs, phase_bins, boxcar_widths, periods, snrs, flags, snr_threshold)

Retain fundamental-flagged candidates whose S/N exceeds a detection threshold.

Parameters:

Name	Type	Description	Default
channels	NDArray[intp]	Spectral channel index of each candidate.	required
radiofreqs	NDArray[floating]	Radio frequency (MHz) of each candidate.	required
phase_bins	NDArray[uint]	Number of phase bins in the folded profile for each candidate.	required
boxcar_widths	NDArray[uint]	Best-fit boxcar widths in phase bins for each candidate.	required
periods	NDArray[floating]	Best-fit periods in seconds for each candidate.	required
snrs	NDArray[floating]	Peak signal-to-noise ratios for each candidate.	required
flags	NDArray[str_]	Harmonic classification label for each candidate.	required
snr_threshold	float	Minimum S/N for a candidate to be retained.	required

Returns:

Type	Description
NDArray[intp]	Tuple of (channels, radiofreqs, phase_bins, boxcar_widths, periods, snrs) for
NDArray[floating]	candidates flagged as fundamental with S/N at or above snr_threshold.

Source code in blipss/core/compare_cands.py

def filter_fundamental_candidates(
    channels: npt.NDArray[np.intp],
    radiofreqs: npt.NDArray[np.floating],
    phase_bins: npt.NDArray[np.uint],
    boxcar_widths: npt.NDArray[np.uint],
    periods: npt.NDArray[np.floating],
    snrs: npt.NDArray[np.floating],
    flags: npt.NDArray[np.str_],
    snr_threshold: float,
) -> tuple[
    npt.NDArray[np.intp],
    npt.NDArray[np.floating],
    npt.NDArray[np.uint],
    npt.NDArray[np.uint],
    npt.NDArray[np.floating],
    npt.NDArray[np.floating],
]:
    """
    Retain fundamental-flagged candidates whose S/N exceeds a detection threshold.

    Args:
        channels: Spectral channel index of each candidate.
        radiofreqs: Radio frequency (MHz) of each candidate.
        phase_bins: Number of phase bins in the folded profile for each candidate.
        boxcar_widths: Best-fit boxcar widths in phase bins for each candidate.
        periods: Best-fit periods in seconds for each candidate.
        snrs: Peak signal-to-noise ratios for each candidate.
        flags: Harmonic classification label for each candidate.
        snr_threshold: Minimum S/N for a candidate to be retained.

    Returns:
        Tuple of (channels, radiofreqs, phase_bins, boxcar_widths, periods, snrs) for
        candidates flagged as fundamental with S/N at or above ``snr_threshold``.
    """
    mask = (snrs >= snr_threshold) & (flags == FUNDAMENTAL_FLAG)
    return channels[mask], radiofreqs[mask], phase_bins[mask], boxcar_widths[mask], periods[mask], snrs[mask]

group_candidates_by_channel(file_index, channels, radiofreqs, phase_bins, boxcar_widths, periods, snrs, n_files, cluster_radius, n_jobs=1)

Cluster candidate periods within each spectral channel and assign a detection code to each.

Candidates are sorted by channel once up front, so each channel's rows are a contiguous slice rather than being re-selected with a fresh channels == ch scan of the full array per channel.

Parameters:

Name	Type	Description	Default
file_index	NDArray[intp]	Source file index of each merged candidate.	required
channels	NDArray[intp]	Spectral channel index of each merged candidate.	required
radiofreqs	NDArray[floating]	Radio frequency (MHz) of each merged candidate.	required
phase_bins	NDArray[uint]	Number of phase bins in the folded profile for each merged candidate.	required
boxcar_widths	NDArray[uint]	Best-fit boxcar widths in phase bins for each merged candidate.	required
periods	NDArray[floating]	Best-fit periods in seconds for each merged candidate.	required
snrs	NDArray[floating]	Peak signal-to-noise ratios for each merged candidate.	required
n_files	int	Total number of input files being compared.	required
cluster_radius	float	Friends-of-Friends clustering radius (s) applied to periods.	required
n_jobs	int	Number of worker processes for channel clustering. 1 (default) runs sequentially in-process; -1 uses all available CPU cores; any other positive value is used as-is.	1

Returns:

Type	Description
NDArray[intp]	Tuple of (channels, radiofreqs, phase_bins, boxcar_widths, periods, snrs, codes)
NDArray[floating]	merged across all channels, one entry per surviving period cluster.

Source code in blipss/core/compare_cands.py

def group_candidates_by_channel(
    file_index: npt.NDArray[np.intp],
    channels: npt.NDArray[np.intp],
    radiofreqs: npt.NDArray[np.floating],
    phase_bins: npt.NDArray[np.uint],
    boxcar_widths: npt.NDArray[np.uint],
    periods: npt.NDArray[np.floating],
    snrs: npt.NDArray[np.floating],
    n_files: int,
    cluster_radius: float,
    n_jobs: int = 1,
) -> tuple[
    npt.NDArray[np.intp],
    npt.NDArray[np.floating],
    npt.NDArray[np.uint],
    npt.NDArray[np.uint],
    npt.NDArray[np.floating],
    npt.NDArray[np.floating],
    npt.NDArray[np.str_],
]:
    """
    Cluster candidate periods within each spectral channel and assign a detection code to each.

    Candidates are sorted by channel once up front, so each channel's rows are a contiguous slice
    rather than being re-selected with a fresh ``channels == ch`` scan of the full array per channel.

    Args:
        file_index: Source file index of each merged candidate.
        channels: Spectral channel index of each merged candidate.
        radiofreqs: Radio frequency (MHz) of each merged candidate.
        phase_bins: Number of phase bins in the folded profile for each merged candidate.
        boxcar_widths: Best-fit boxcar widths in phase bins for each merged candidate.
        periods: Best-fit periods in seconds for each merged candidate.
        snrs: Peak signal-to-noise ratios for each merged candidate.
        n_files: Total number of input files being compared.
        cluster_radius: Friends-of-Friends clustering radius (s) applied to periods.
        n_jobs: Number of worker processes for channel clustering. 1 (default) runs sequentially
            in-process; -1 uses all available CPU cores; any other positive value is used as-is.

    Returns:
        Tuple of (channels, radiofreqs, phase_bins, boxcar_widths, periods, snrs, codes)
        merged across all channels, one entry per surviving period cluster.
    """
    order = np.argsort(channels, kind="stable")
    s_file_index = file_index[order]
    s_radiofreqs = radiofreqs[order]
    s_phase_bins = phase_bins[order]
    s_boxcar_widths = boxcar_widths[order]
    s_periods = periods[order]
    s_snrs = snrs[order]
    s_channels = channels[order]

    unique_channels, starts = np.unique(s_channels, return_index=True)
    bounds = np.append(starts, len(s_channels))

    n_jobs = (os.cpu_count() or 1) if n_jobs == -1 else max(n_jobs, 1)

    if n_jobs == 1:
        chunk_results = [
            _process_channel_range(
                unique_channels,
                bounds,
                s_file_index,
                s_radiofreqs,
                s_phase_bins,
                s_boxcar_widths,
                s_periods,
                s_snrs,
                n_files,
                cluster_radius,
                show_progress=True,
            )
        ]
    else:
        chunk_indices = [c for c in np.array_split(np.arange(len(unique_channels)), n_jobs) if len(c)]
        with ProcessPoolExecutor(max_workers=n_jobs) as executor:
            futures = [
                executor.submit(
                    _process_channel_range,
                    unique_channels[idx],
                    bounds[idx[0] : idx[-1] + 2] - bounds[idx[0]],
                    s_file_index[bounds[idx[0]] : bounds[idx[-1] + 1]],
                    s_radiofreqs[bounds[idx[0]] : bounds[idx[-1] + 1]],
                    s_phase_bins[bounds[idx[0]] : bounds[idx[-1] + 1]],
                    s_boxcar_widths[bounds[idx[0]] : bounds[idx[-1] + 1]],
                    s_periods[bounds[idx[0]] : bounds[idx[-1] + 1]],
                    s_snrs[bounds[idx[0]] : bounds[idx[-1] + 1]],
                    n_files,
                    cluster_radius,
                )
                for idx in chunk_indices
            ]
            chunk_results = [f.result() for f in tqdm(futures)]

    out_channels: list[int] = []
    out_radiofreqs: list[float] = []
    out_phase_bins: list[int] = []
    out_boxcar_widths: list[int] = []
    out_periods: list[float] = []
    out_snrs: list[float] = []
    out_codes: list[str] = []
    for (
        res_channels,
        res_radiofreqs,
        res_phase_bins,
        res_boxcar_widths,
        res_periods,
        res_snrs,
        res_codes,
    ) in chunk_results:
        out_channels.extend(res_channels)
        out_radiofreqs.extend(res_radiofreqs)
        out_phase_bins.extend(res_phase_bins)
        out_boxcar_widths.extend(res_boxcar_widths)
        out_periods.extend(res_periods)
        out_snrs.extend(res_snrs)
        out_codes.extend(res_codes)

    return (
        np.array(out_channels, dtype=np.intp),
        np.array(out_radiofreqs, dtype=np.float64),
        np.array(out_phase_bins, dtype=np.uint),
        np.array(out_boxcar_widths, dtype=np.uint),
        np.array(out_periods, dtype=np.float64),
        np.array(out_snrs, dtype=np.float64),
        np.array(out_codes, dtype=np.str_),
    )

blipss.core.compute_phase_resolved_ds

Atomic data-manipulation routines for the phase-resolved dynamic spectrum pipeline.

These functions implement the internal (n_channels, n_samples) array representation:

1. extract_waterfall_metadata: Pull start MJD, tsamp, and frequency axis from a Waterfall header.
2. align_band_orientation: Flip data and frequency axis so frequencies increase monotonically.
3. clip_channels: Select a contiguous channel sub-band by index.
4. fold_all_channels: Fold each channel time series in parallel with riptide and return the 2-D phase-resolved DS.

align_band_orientation(data, freqs_MHz, foff)

Flip data and frequency arrays so that channel index 0 corresponds to the lowest frequency.

When foff is negative the Waterfall stores channels in descending frequency order. This function normalises to ascending order so that downstream code can treat index 0 as the lowest frequency unconditionally.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	2-D data array of shape (n_channels, n_samples).	required
freqs_MHz	NDArray[floating]	1-D frequency array of shape (n_channels,).	required
foff	float	Channel bandwidth in MHz; negative indicates descending frequency order.	required

Returns:

Type	Description
tuple[NDArray[floating], NDArray[floating]]	Tuple of (data, freqs_MHz) flipped along the channel axis when foff < 0, otherwise unchanged.

Source code in blipss/core/compute_phase_resolved_ds.py

def align_band_orientation(
    data: npt.NDArray[np.floating],
    freqs_MHz: npt.NDArray[np.floating],
    foff: float,
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
    """
    Flip data and frequency arrays so that channel index 0 corresponds to the lowest frequency.

    When foff is negative the Waterfall stores channels in descending frequency order.
    This function normalises to ascending order so that downstream code can treat
    index 0 as the lowest frequency unconditionally.

    Args:
        data: 2-D data array of shape (n_channels, n_samples).
        freqs_MHz: 1-D frequency array of shape (n_channels,).
        foff: Channel bandwidth in MHz; negative indicates descending frequency order.

    Returns:
        Tuple of (data, freqs_MHz) flipped along the channel axis when foff < 0, otherwise unchanged.
    """
    if foff < 0:
        return np.flip(data, axis=0), np.flip(freqs_MHz)
    return data, freqs_MHz

clip_channels(data, freqs_MHz, start_ch, stop_ch)

Restrict data and frequency arrays to a contiguous range of channel indices.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	2-D data array of shape (n_channels, n_samples).	required
freqs_MHz	NDArray[floating]	1-D frequency array of shape (n_channels,).	required
start_ch	int	First channel index to include (inclusive).	required
stop_ch	int | None	Last channel index to exclude (exclusive); None retains all remaining channels.	required

Returns:

Type	Description
tuple[NDArray[floating], NDArray[floating]]	Tuple of (data[start_ch:stop_ch], freqs_MHz[start_ch:stop_ch]).

Raises:

Type	Description
ValueError	If start_ch is outside [0, n_channels) or stop_ch does not satisfy start_ch < stop_ch <= n_channels.

Source code in blipss/core/compute_phase_resolved_ds.py

def clip_channels(
    data: npt.NDArray[np.floating],
    freqs_MHz: npt.NDArray[np.floating],
    start_ch: int,
    stop_ch: int | None,
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
    """
    Restrict data and frequency arrays to a contiguous range of channel indices.

    Args:
        data: 2-D data array of shape (n_channels, n_samples).
        freqs_MHz: 1-D frequency array of shape (n_channels,).
        start_ch: First channel index to include (inclusive).
        stop_ch: Last channel index to exclude (exclusive); None retains all remaining channels.

    Returns:
        Tuple of (data[start_ch:stop_ch], freqs_MHz[start_ch:stop_ch]).

    Raises:
        ValueError: If start_ch is outside [0, n_channels) or stop_ch does not satisfy start_ch < stop_ch <= n_channels.
    """
    n_channels = len(data)
    if start_ch < 0 or start_ch >= n_channels:
        raise ValueError(f"start_ch {start_ch} is out of range [0, {n_channels})")
    if stop_ch is not None and (stop_ch <= start_ch or stop_ch > n_channels):
        raise ValueError(f"stop_ch {stop_ch} must satisfy {start_ch} < stop_ch <= {n_channels}")
    return data[start_ch:stop_ch], freqs_MHz[start_ch:stop_ch]

extract_waterfall_metadata(wat)

Extract the radio frequency axis, observation start MJD, and sampling interval from a Waterfall header.

Parameters:

Name	Type	Description	Default
wat	Waterfall	Blimpy Waterfall object containing the loaded filterbank data.	required

Returns:

Type	Description
NDArray[floating]	Tuple of
float	freqs_MHz: Radio frequency axis in MHz,
float	start_mjd: Observation start as MJD (UTC),
tuple[NDArray[floating], float, float]	tsamp: sampling interval in seconds.

Source code in blipss/core/compute_phase_resolved_ds.py

def extract_waterfall_metadata(
    wat: Waterfall,
) -> tuple[npt.NDArray[np.floating], float, float]:
    """
    Extract the radio frequency axis, observation start MJD, and sampling interval from a Waterfall header.

    Args:
        wat: Blimpy Waterfall object containing the loaded filterbank data.

    Returns:
        Tuple of
        - freqs_MHz: Radio frequency axis in MHz,
        - start_mjd: Observation start as MJD (UTC),
        - tsamp: sampling interval in seconds.
    """
    # Standard filterbank convention: fch1 is the centre frequency of the first channel; foff is the channel spacing.
    freqs_MHz: npt.NDArray[np.floating] = wat.header["fch1"] + np.arange(wat.header["nchans"]) * wat.header["foff"]
    start_mjd: float = wat.header["tstart"]
    tsamp: float = wat.header["tsamp"]
    return freqs_MHz, start_mjd, tsamp

fold_all_channels(data, tsamp, period, bins, do_deredden, rmed_width, n_workers=None)

Fold each spectral channel's time series in parallel to produce a phase-resolved dynamic spectrum.

Ensures a C-contiguous memory layout before dispatching channel rows to worker processes. Each channel is optionally detrended with a running-median filter, normalised to zero median and unit standard deviation, then folded using riptide.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	2-D data array of shape (n_channels, n_samples).	required
tsamp	float	Sampling interval in seconds.	required
period	float	Folding period in seconds.	required
bins	int	Number of phase bins in the folded profile.	required
do_deredden	bool	Apply running-median detrending before folding when True.	required
rmed_width	float	Running median window width in seconds; used only when do_deredden is True.	required
n_workers	int | None	Number of parallel worker processes; None uses all available CPUs.	None

Returns:

Type	Description
NDArray[floating]	2-D phase-resolved dynamic spectrum of shape (n_channels, bins).

Source code in blipss/core/compute_phase_resolved_ds.py

def fold_all_channels(
    data: npt.NDArray[np.floating],
    tsamp: float,
    period: float,
    bins: int,
    do_deredden: bool,
    rmed_width: float,
    n_workers: int | None = None,
) -> npt.NDArray[np.floating]:
    """
    Fold each spectral channel's time series in parallel to produce a phase-resolved dynamic spectrum.

    Ensures a C-contiguous memory layout before dispatching channel rows to worker processes.
    Each channel is optionally detrended with a running-median filter, normalised to
    zero median and unit standard deviation, then folded using riptide.

    Args:
        data: 2-D data array of shape (n_channels, n_samples).
        tsamp: Sampling interval in seconds.
        period: Folding period in seconds.
        bins: Number of phase bins in the folded profile.
        do_deredden: Apply running-median detrending before folding when True.
        rmed_width: Running median window width in seconds; used only when do_deredden is True.
        n_workers: Number of parallel worker processes; None uses all available CPUs.

    Returns:
        2-D phase-resolved dynamic spectrum of shape (n_channels, bins).
    """
    # np.flip returns negatively-strided views; contiguous layout avoids a per-worker copy on each row access.
    data = np.ascontiguousarray(data)
    n_channels = len(data)
    fold_fn = partial(
        _fold_single_channel,
        tsamp=tsamp,
        period=period,
        bins=bins,
        do_deredden=do_deredden,
        rmed_width=rmed_width,
    )
    # Batching amortizes IPC overhead; chunksize=1 would mean one round-trip per channel.
    chunksize = max(1, n_channels // ((n_workers or os.cpu_count() or 1) * 4))
    profiles = process_map(
        fold_fn,
        data,
        max_workers=n_workers,
        chunksize=chunksize,
        total=n_channels,
    )
    return np.array(profiles, dtype=np.float32)

blipss.core.harmonic_detection

Utilities for identifying and labeling harmonic and sub-harmonic relationships among a set of candidate periods.

Terminology

Given a fundamental period p0:

• The (N-1)th harmonic of p0 is a period P such that P ≈ p0 / N (N >= 2).
• The (N-1)th sub-harmonic of p0 is a period P such that P ≈ N * p0 (N >= 2).

Algorithm

label_harmonics uses a greedy approach: the unlabeled period with the highest S/N is selected as a fundamental ('F'), all of its harmonics ('H') and sub-harmonics ('S') present in the input are labeled and removed from consideration, then the process repeats on the remaining periods until every period has been assigned a label.

label_harmonics(periods, snrs, epsilon_harmonic=DEFAULT_EPSILON_HARMONIC, presorted=False)

label_harmonics(
    periods: npt.NDArray[np.floating],
    snrs: npt.NDArray[np.floating],
    epsilon_harmonic: float = ...,
    presorted: Literal[False] = ...,
) -> tuple[
    npt.NDArray[np.str_],
    npt.NDArray[np.floating],
    npt.NDArray[np.floating],
]

label_harmonics(
    periods: npt.NDArray[np.floating],
    snrs: npt.NDArray[np.floating],
    epsilon_harmonic: float = ...,
    *,
    presorted: Literal[True],
) -> npt.NDArray[np.str_]

Assign harmonic labels to an array of candidate periods.

Uses a greedy, highest-S/N-first algorithm (see module docstring). Each iteration selects the unlabeled period with the highest S/N as a fundamental, labels its harmonics and sub-harmonics, removes all of them from the active set, then repeats until no unlabeled periods remain.

A period P is the (N-1)th harmonic of p0 if |P - p0/N| <= epsilon_harmonic. A period P is the (N-1)th sub-harmonic of p0 if |P - N*p0| <= epsilon_harmonic.

Parameters:

Name	Type	Description	Default
periods	NDArray[floating]	Periods (s)	required
snrs	NDArray[floating]	S/N values associated with the above periods	required
epsilon_harmonic	float	Floating-point tolerance for harmonic period matching	DEFAULT_EPSILON_HARMONIC
presorted	bool	Set to True if periods are already sorted in descending S/N order. Skips the sort and returns only the flags array, which is useful when calling in a tight loop where the sort has already been done.	False

Returns:

Type	Description
tuple[NDArray[str_], NDArray[floating], NDArray[floating]] | NDArray[str_]	If presorted is False: tuple of (flags, sorted_periods, sorted_snrs)
tuple[NDArray[str_], NDArray[floating], NDArray[floating]] | NDArray[str_]	where periods and S/N arrays are in descending S/N order.
tuple[NDArray[str_], NDArray[floating], NDArray[floating]] | NDArray[str_]	If presorted is True: the flags array only (no copy of the sorted
tuple[NDArray[str_], NDArray[floating], NDArray[floating]] | NDArray[str_]	input arrays is returned).
tuple[NDArray[str_], NDArray[floating], NDArray[floating]] | NDArray[str_]	Flags are single characters: 'F' (fundamental), 'H' (harmonic), or
tuple[NDArray[str_], NDArray[floating], NDArray[floating]] | NDArray[str_]	'S' (sub-harmonic).

Source code in blipss/core/harmonic_detection.py

def label_harmonics(
    periods: npt.NDArray[np.floating],
    snrs: npt.NDArray[np.floating],
    epsilon_harmonic: float = DEFAULT_EPSILON_HARMONIC,
    presorted: bool = False,
) -> tuple[npt.NDArray[np.str_], npt.NDArray[np.floating], npt.NDArray[np.floating]] | npt.NDArray[np.str_]:
    """
    Assign harmonic labels to an array of candidate periods.

    Uses a greedy, highest-S/N-first algorithm (see module docstring). Each
    iteration selects the unlabeled period with the highest S/N as a
    fundamental, labels its harmonics and sub-harmonics, removes all of them
    from the active set, then repeats until no unlabeled periods remain.

    A period P is the (N-1)th harmonic of p0 if ``|P - p0/N| <= epsilon_harmonic``.
    A period P is the (N-1)th sub-harmonic of p0 if ``|P - N*p0| <= epsilon_harmonic``.

    Args:
        periods: Periods (s)
        snrs: S/N values associated with the above periods
        epsilon_harmonic: Floating-point tolerance for harmonic period matching
        presorted: Set to True if periods are already sorted in descending S/N
            order. Skips the sort and returns only the flags array, which is
            useful when calling in a tight loop where the sort has already been
            done.

    Returns:
        If presorted is False: tuple of (flags, sorted_periods, sorted_snrs)
        where periods and S/N arrays are in descending S/N order.

        If presorted is True: the flags array only (no copy of the sorted
        input arrays is returned).

        Flags are single characters: 'F' (fundamental), 'H' (harmonic), or
        'S' (sub-harmonic).
    """
    if not presorted:
        periods, snrs = _sort_by_snr(periods, snrs)

    temp_periods: list[npt.NDArray[np.floating]] = []
    temp_snrs: list[npt.NDArray[np.floating]] = []
    harm_flag: list[str] = []

    active = np.ones(len(periods), dtype=bool)

    while np.any(active):
        active_idx = np.where(active)[0]
        rem_p = periods[active_idx]

        p0 = float(rem_p[0])
        idx_subharm = _find_subharmonic_indices(rem_p, p0, epsilon_harmonic)
        idx_harm = _find_harmonic_indices(rem_p, p0, epsilon_harmonic)

        # idx_subharm and idx_harm are indices into rem_p (the active subset);
        # map them back to the original periods array via active_idx.
        local_idx = np.concatenate([idx_subharm, idx_harm])
        global_idx = active_idx[local_idx]

        temp_periods.append(periods[global_idx])
        temp_snrs.append(snrs[global_idx])

        # _find_subharmonic_indices guarantees idx_subharm[0] is p0 itself,
        # so the flag list always starts with exactly one fundamental flag.
        flags: list[str] = (
            [FUNDAMENTAL_FLAG] + [SUBHARMONIC_FLAG] * (len(idx_subharm) - 1) + [HARMONIC_FLAG] * len(idx_harm)
        )
        harm_flag.extend(flags)

        active[global_idx] = False

    if not temp_periods:
        empty: npt.NDArray[np.str_] = np.array([], dtype="<U1")
        if presorted:
            return empty
        return empty, np.array([]), np.array([])

    all_periods = np.concatenate(temp_periods)
    all_snrs = np.concatenate(temp_snrs)
    harm_flag_arr: npt.NDArray[np.str_] = np.array(harm_flag)

    sort_idx = np.argsort(-all_snrs)
    sorted_periods = all_periods[sort_idx]
    sorted_snrs = all_snrs[sort_idx]
    harm_flag_arr = harm_flag_arr[sort_idx]

    if presorted:
        return harm_flag_arr
    return harm_flag_arr, sorted_periods, sorted_snrs

blipss.core.inject_signal

Atomic data-manipulation routines for the real-data signal-injection pipeline.

These functions operate on the internal (n_channels, n_samples) array representation used throughout the inject_signal pipeline:

1. extract_data_array`: Load aWaterfall`` object into the internal layout.
2. compute_median_bandpass / compute_per_channel_std: Characterise the per-channel noise statistics of the real data; used to calibrate injected pulse amplitudes relative to the local bandpass.
3. pack_data_into_waterfall: Convert the modified array back to the sigproc (n_samples, n_ifs, n_channels) layout and store it in the Waterfall object prior to writing.

compute_median_bandpass(data)

Compute the per-channel median value across all time samples.

Used to estimate the baseline level of each spectral channel so that injected pulse amplitudes can be set relative to the local bandpass rather than an absolute scale.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	Array of shape (n_channels, n_samples).	required

Returns:

Type	Description
NDArray[floating]	1-D array of per-channel median values, shape (n_channels,).

Source code in blipss/core/inject_signal.py

def compute_median_bandpass(
    data: npt.NDArray[np.floating],
) -> npt.NDArray[np.floating]:
    """
    Compute the per-channel median value across all time samples.

    Used to estimate the baseline level of each spectral channel so that injected
    pulse amplitudes can be set relative to the local bandpass rather than an
    absolute scale.

    Args:
        data: Array of shape ``(n_channels, n_samples)``.

    Returns:
        1-D array of per-channel median values, shape ``(n_channels,)``.
    """
    return np.median(data, axis=1)

compute_per_channel_std(data)

Compute the per-channel standard deviation across all time samples.

Used alongside compute_median_bandpass to express the injected pulse amplitude in units of the local noise standard deviation (i.e., as an SNR).

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	Array of shape (n_channels, n_samples).	required

Returns:

Type	Description
NDArray[floating]	1-D array of per-channel standard deviations, shape (n_channels,).

Source code in blipss/core/inject_signal.py

def compute_per_channel_std(
    data: npt.NDArray[np.floating],
) -> npt.NDArray[np.floating]:
    """
    Compute the per-channel standard deviation across all time samples.

    Used alongside ``compute_median_bandpass`` to express the injected pulse
    amplitude in units of the local noise standard deviation (i.e., as an SNR).

    Args:
        data: Array of shape ``(n_channels, n_samples)``.

    Returns:
        1-D array of per-channel standard deviations, shape ``(n_channels,)``.
    """
    return np.std(data, axis=1)

extract_data_array(wat, if_channel=0)

Extract the 2-D data array, sample count, and sampling interval from a Waterfall object.

Converts from blimpy's on-disk (n_samples, n_ifs, n_channels) layout to the internal (n_channels, n_samples) representation used by the rest of the inject_signal pipeline.

Parameters:

Name	Type	Description	Default
wat	Waterfall	Blimpy Waterfall object containing the loaded filterbank data.	required
if_channel	int	Index of the polarisation (IF) channel to extract. Defaults to 0.	0

Returns:

Type	Description
NDArray[floating]	Tuple of:
int	data: 2-D array of shape (n_channels, n_samples)
float	n_samples: Number of time samples in the file
tuple[NDArray[floating], int, float]	tsamp: Sampling interval in seconds

Source code in blipss/core/inject_signal.py

def extract_data_array(
    wat: Waterfall,
    if_channel: int = 0,
) -> tuple[npt.NDArray[np.floating], int, float]:
    """
    Extract the 2-D data array, sample count, and sampling interval from a Waterfall object.

    Converts from blimpy's on-disk ``(n_samples, n_ifs, n_channels)`` layout to the
    internal ``(n_channels, n_samples)`` representation used by the rest of the
    inject_signal pipeline.

    Args:
        wat: Blimpy ``Waterfall`` object containing the loaded filterbank data.
        if_channel: Index of the polarisation (IF) channel to extract. Defaults to 0.

    Returns:
        Tuple of:
        - data: 2-D array of shape ``(n_channels, n_samples)``
        - n_samples: Number of time samples in the file
        - tsamp: Sampling interval in seconds
    """
    n_samples: int = wat.n_ints_in_file
    tsamp: float = wat.header["tsamp"]
    # wat.data has shape (n_samples, n_ifs, n_channels)
    data: npt.NDArray[np.floating] = wat.data[:, if_channel, :].T
    return data, n_samples, tsamp

pack_data_into_waterfall(data, wat, n_samples)

Reshape the modified data array back into sigproc layout and store it in the Waterfall object.

This is the inverse of extract_data_array: it converts from the internal (n_channels, n_samples) representation back to the (n_samples, n_ifs, n_channels) layout expected by blimpy's serialisers. n_ifs is read from the original Waterfall header so that multi-polarisation files are handled correctly.

Call this immediately before blipss.io.write_filterbank.write_waterfall.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	Array of shape (n_channels, n_samples) after signal injection.	required
wat	Waterfall	Blimpy Waterfall object whose .data attribute will be overwritten. The header is not modified.	required
n_samples	int	Number of time samples (as returned by extract_data_array).	required

Returns:

Type	Description
Waterfall	The same Waterfall object with its .data attribute updated.

Source code in blipss/core/inject_signal.py

def pack_data_into_waterfall(
    data: npt.NDArray[np.floating],
    wat: Waterfall,
    n_samples: int,
) -> Waterfall:
    """
    Reshape the modified data array back into sigproc layout and store it in the Waterfall object.

    This is the inverse of ``extract_data_array``: it converts from the internal
    ``(n_channels, n_samples)`` representation back to the ``(n_samples, n_ifs, n_channels)``
    layout expected by blimpy's serialisers. ``n_ifs`` is read from the original
    Waterfall header so that multi-polarisation files are handled correctly.

    Call this immediately before ``blipss.io.write_filterbank.write_waterfall``.

    Args:
        data: Array of shape ``(n_channels, n_samples)`` after signal injection.
        wat: Blimpy ``Waterfall`` object whose ``.data`` attribute will be overwritten.
            The header is not modified.
        n_samples: Number of time samples (as returned by ``extract_data_array``).

    Returns:
        The same ``Waterfall`` object with its ``.data`` attribute updated.
    """
    n_ifs: int = wat.header["nifs"]
    n_channels: int = wat.header["nchans"]
    wat.data = data.T.reshape((n_samples, n_ifs, n_channels))
    return wat

blipss.core.period_finding

Period search across spectral channels using the Fast Folding Algorithm (FFA).

Pipeline

For each spectral channel the following steps are applied in sequence:

1. FFA search: Dold the time series at trial periods and compute matched-filter S/N.
2. Threshold: Retain only periods whose peak S/N exceeds a detection threshold.
3. Cluster: Group nearby periods via Friends-of-Friends and keep the highest-S/N representative per cluster.
4. Label harmonics: Classify each surviving candidate as a fundamental, harmonic, or sub-harmonic.

Results from all channels are merged and returned by search_all_channels.

search_all_channels(data, start_channel, sampling_time_in_seconds, minimum_period_in_seconds, maximum_period_in_seconds, minimum_fold_periods, minimum_bins, maximum_bins, max_duty_cycle, do_deredden, running_median_width_in_seconds, snr_threshold, epsilon_fof, epsilon_harmonic, n_workers=None)

Search every spectral channel for periodic signals using the Fast Folding Algorithm.

Runs an FFA search on each channel of data, clusters and deduplicates period candidates, labels harmonics, and returns merged results across all channels.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	2D array of shape (n_channels, n_samples); each row is one spectral channel	required
start_channel	int	Global channel index of data[0]; offsets channel numbers in the output	required
sampling_time_in_seconds	float	Sampling time (s)	required
minimum_period_in_seconds	float	Minimum trial period (s) for the FFA search	required
maximum_period_in_seconds	float	Maximum trial period (s) for the FFA search	required
minimum_fold_periods	int	Minimum number of signal periods that must fit in the data	required
minimum_bins	int	Minimum number of phase bins across the full [0, 1] phase range; a folded profile may cover only a fraction of this range	required
maximum_bins	int	Maximum number of phase bins across the full [0, 1] phase range; a folded profile may cover only a fraction of this range	required
max_duty_cycle	float	Maximum duty cycle searched	required
do_deredden	bool	Whether to detrend each time series with a running median filter	required
running_median_width_in_seconds	float	Running median window width (s)	required
snr_threshold	float	Minimum matched-filtering S/N for a detection	required
epsilon_fof	float	Period tolerance for Friends-of-Friends clustering	required
epsilon_harmonic	float	Period tolerance for harmonic matching	required
n_workers	int | None	Number of parallel worker processes; 1 runs the serial tqdm loop, None or values >1 dispatch via process_map using all available CPUs (None) or the specified count	None

Returns:

Type	Description
NDArray[intp]	Tuple of (cand_channels, cand_periods, cand_snrs, cand_phase_bins,
NDArray[floating]	cand_boxcar_widths, cand_flags) across all channels.
NDArray[floating]	cand_channels: Spectral channel index of each candidate
NDArray[uint]	cand_periods: Trial periods (s) at which the peak S/N across all boxcar widths exceeded snr_threshold
NDArray[uint]	cand_snrs: Peak S/N for each candidate, maximised over all trial boxcar widths
NDArray[str_]	cand_phase_bins: Number of phase bins across the full [0, 1] phase range to be used to generate the folded profile for each candidate
tuple[NDArray[intp], NDArray[floating], NDArray[floating], NDArray[uint], NDArray[uint], NDArray[str_]]	cand_boxcar_widths: Width (in phase bins) of the boxcar matched filter that produced the peak S/N for each candidate; the implied duty cycle is cand_boxcar_widths / cand_phase_bins
tuple[NDArray[intp], NDArray[floating], NDArray[floating], NDArray[uint], NDArray[uint], NDArray[str_]]	cand_flags: Harmonic classification for each candidate ('F': fundamental, 'H': harmonic, 'S': sub-harmonic)

Source code in blipss/core/period_finding.py

def search_all_channels(
    data: npt.NDArray[np.floating],
    start_channel: int,
    sampling_time_in_seconds: float,
    minimum_period_in_seconds: float,
    maximum_period_in_seconds: float,
    minimum_fold_periods: int,
    minimum_bins: int,
    maximum_bins: int,
    max_duty_cycle: float,
    do_deredden: bool,
    running_median_width_in_seconds: float,
    snr_threshold: float,
    epsilon_fof: float,
    epsilon_harmonic: float,
    n_workers: int | None = None,
) -> tuple[
    npt.NDArray[np.intp],
    npt.NDArray[np.floating],
    npt.NDArray[np.floating],
    npt.NDArray[np.uint],
    npt.NDArray[np.uint],
    npt.NDArray[np.str_],
]:
    """
    Search every spectral channel for periodic signals using the Fast Folding Algorithm.

    Runs an FFA search on each channel of data, clusters and deduplicates period candidates,
    labels harmonics, and returns merged results across all channels.

    Args:
        data: 2D array of shape (n_channels, n_samples); each row is one spectral channel
        start_channel: Global channel index of data[0]; offsets channel numbers in the output
        sampling_time_in_seconds: Sampling time (s)
        minimum_period_in_seconds: Minimum trial period (s) for the FFA search
        maximum_period_in_seconds: Maximum trial period (s) for the FFA search
        minimum_fold_periods: Minimum number of signal periods that must fit in the data
        minimum_bins: Minimum number of phase bins across the full [0, 1] phase range;
            a folded profile may cover only a fraction of this range
        maximum_bins: Maximum number of phase bins across the full [0, 1] phase range;
            a folded profile may cover only a fraction of this range
        max_duty_cycle: Maximum duty cycle searched
        do_deredden: Whether to detrend each time series with a running median filter
        running_median_width_in_seconds: Running median window width (s)
        snr_threshold: Minimum matched-filtering S/N for a detection
        epsilon_fof: Period tolerance for Friends-of-Friends clustering
        epsilon_harmonic: Period tolerance for harmonic matching
        n_workers: Number of parallel worker processes; 1 runs the serial tqdm loop,
            None or values >1 dispatch via process_map using all available CPUs (None)
            or the specified count

    Returns:
        Tuple of (cand_channels, cand_periods, cand_snrs, cand_phase_bins,
        cand_boxcar_widths, cand_flags) across all channels.
        - cand_channels: Spectral channel index of each candidate
        - cand_periods: Trial periods (s) at which the peak S/N across all boxcar widths
            exceeded snr_threshold
        - cand_snrs: Peak S/N for each candidate, maximised over all trial boxcar widths
        - cand_phase_bins: Number of phase bins across the full [0, 1] phase range to be
            used to generate the folded profile for each candidate
        - cand_boxcar_widths: Width (in phase bins) of the boxcar matched filter that
            produced the peak S/N for each candidate; the implied duty cycle is
            cand_boxcar_widths / cand_phase_bins
        - cand_flags: Harmonic classification for each candidate ('F': fundamental,
            'H': harmonic, 'S': sub-harmonic)
    """
    # np.flip returns negatively-strided views; contiguous layout avoids a per-worker copy on each row access.
    data = np.ascontiguousarray(data)
    n_channels = len(data)

    cand_channels: list[int] = []
    cand_periods: list[float] = []
    cand_snrs: list[float] = []
    cand_phase_bins: list[int] = []
    cand_boxcar_widths: list[int] = []
    cand_flags: list[str] = []

    if n_workers == 1:
        for ch_idx in tqdm(range(n_channels)):
            result = _search_single_channel(
                data[ch_idx],
                sampling_time_in_seconds,
                minimum_period_in_seconds,
                maximum_period_in_seconds,
                minimum_fold_periods,
                minimum_bins,
                maximum_bins,
                max_duty_cycle,
                do_deredden,
                running_median_width_in_seconds,
                snr_threshold,
                epsilon_fof,
                epsilon_harmonic,
            )
            if result is None:
                continue
            periods, snrs, phase_bins, boxcar_widths, flags = result
            n = len(periods)
            cand_channels.extend([start_channel + ch_idx] * n)
            cand_periods.extend(periods.tolist())
            cand_snrs.extend(snrs.tolist())
            cand_phase_bins.extend(phase_bins.tolist())
            cand_boxcar_widths.extend(boxcar_widths.tolist())
            cand_flags.extend(flags.tolist())
    else:
        search_fn = partial(
            _search_single_channel,
            sampling_time_in_seconds=sampling_time_in_seconds,
            minimum_period_in_seconds=minimum_period_in_seconds,
            maximum_period_in_seconds=maximum_period_in_seconds,
            minimum_fold_periods=minimum_fold_periods,
            minimum_bins=minimum_bins,
            maximum_bins=maximum_bins,
            max_duty_cycle=max_duty_cycle,
            do_deredden=do_deredden,
            running_median_width_in_seconds=running_median_width_in_seconds,
            snr_threshold=snr_threshold,
            epsilon_fof=epsilon_fof,
            epsilon_harmonic=epsilon_harmonic,
        )
        # Batching amortizes IPC overhead; chunksize=1 would mean one round-trip per channel.
        chunksize = max(1, n_channels // ((n_workers or os.cpu_count() or 1) * 4))
        results = process_map(search_fn, data, max_workers=n_workers, chunksize=chunksize, total=n_channels)
        for ch_idx, result in enumerate(results):
            if result is None:
                continue
            periods, snrs, phase_bins, boxcar_widths, flags = result
            n = len(periods)
            cand_channels.extend([start_channel + ch_idx] * n)
            cand_periods.extend(periods.tolist())
            cand_snrs.extend(snrs.tolist())
            cand_phase_bins.extend(phase_bins.tolist())
            cand_boxcar_widths.extend(boxcar_widths.tolist())
            cand_flags.extend(flags.tolist())

    channels_arr, periods_arr, snrs_arr, phase_bins_arr, boxcar_widths_arr = _finalize_candidate_arrays(
        np.array(cand_channels),
        np.array(cand_periods),
        np.array(cand_snrs),
        np.array(cand_phase_bins),
        np.array(cand_boxcar_widths),
    )
    return channels_arr, periods_arr, snrs_arr, phase_bins_arr, boxcar_widths_arr, np.array(cand_flags)

blipss.core.plot_cands

Core logic for selecting and folding periodicity candidates for verification plots.

run_ffa_and_fold_channel(channel_data, tsamp, min_period, max_period, fpmin, bins_min, bins_max, ducy_max, do_deredden, rmed_width)

Run an FFA search on a single-channel time series, keeping both the detrended series and periodogram.

Parameters:

Name	Type	Description	Default
channel_data	NDArray[floating]	1D array of flux density samples for one spectral channel	required
tsamp	float	Sampling time (s)	required
min_period	float	Minimum trial period (s)	required
max_period	float	Maximum trial period (s)	required
fpmin	int	Minimum number of signal periods that must fit in the data duration	required
bins_min	int	Minimum number of phase bins across the full [0, 1] phase range	required
bins_max	int	Maximum number of phase bins across the full [0, 1] phase range	required
ducy_max	float	Maximum duty cycle searched	required
do_deredden	bool	Whether to detrend the time series with a running median filter	required
rmed_width	float	Running median window width (s)	required

Returns:

Type	Description
tuple[TimeSeries, Periodogram]	Tuple of (detrended_ts, periodogram) from the FFA search

Source code in blipss/core/plot_cands.py

def run_ffa_and_fold_channel(
    channel_data: npt.NDArray[np.floating],
    tsamp: float,
    min_period: float,
    max_period: float,
    fpmin: int,
    bins_min: int,
    bins_max: int,
    ducy_max: float,
    do_deredden: bool,
    rmed_width: float,
) -> tuple[TimeSeries, Periodogram]:
    """
    Run an FFA search on a single-channel time series, keeping both the detrended series and periodogram.

    Args:
        channel_data: 1D array of flux density samples for one spectral channel
        tsamp: Sampling time (s)
        min_period: Minimum trial period (s)
        max_period: Maximum trial period (s)
        fpmin: Minimum number of signal periods that must fit in the data duration
        bins_min: Minimum number of phase bins across the full [0, 1] phase range
        bins_max: Maximum number of phase bins across the full [0, 1] phase range
        ducy_max: Maximum duty cycle searched
        do_deredden: Whether to detrend the time series with a running median filter
        rmed_width: Running median window width (s)

    Returns:
        Tuple of (detrended_ts, periodogram) from the FFA search
    """
    raw_ts = TimeSeries.from_numpy_array(channel_data, tsamp=tsamp)
    detrended_ts, periodogram = ffa_search(
        raw_ts,
        period_min=min_period,
        period_max=max_period,
        fpmin=fpmin,
        bins_min=bins_min,
        bins_max=bins_max,
        ducy_max=ducy_max,  # riptide abbreviates "duty cycle" as "ducy"
        deredden=do_deredden,
        rmed_width=rmed_width,
        already_normalised=False,
    )
    return detrended_ts, periodogram

select_candidates_by_code(cand_channels, cand_periods, cand_bins, cand_codes, codes_plot)

Retain candidates whose binary detection code is one of the codes selected for plotting.

Parameters:

Name	Type	Description	Default
cand_channels	NDArray[intp]	Spectral channel index of each candidate	required
cand_periods	NDArray[floating]	Best-fit period (s) of each candidate	required
cand_bins	NDArray[uint]	Number of phase bins in the folded profile for each candidate	required
cand_codes	NDArray[str_]	Per-file binary detection code string for each candidate	required
codes_plot	Sequence[str]	Binary codes selected for plotting	required

Returns:

Type	Description
tuple[NDArray[intp], NDArray[floating], NDArray[uint], NDArray[str_]]	Tuple of (channels, periods, bins, codes) for candidates matching codes_plot

Source code in blipss/core/plot_cands.py

def select_candidates_by_code(
    cand_channels: npt.NDArray[np.intp],
    cand_periods: npt.NDArray[np.floating],
    cand_bins: npt.NDArray[np.uint],
    cand_codes: npt.NDArray[np.str_],
    codes_plot: Sequence[str],
) -> tuple[
    npt.NDArray[np.intp],
    npt.NDArray[np.floating],
    npt.NDArray[np.uint],
    npt.NDArray[np.str_],
]:
    """
    Retain candidates whose binary detection code is one of the codes selected for plotting.

    Args:
        cand_channels: Spectral channel index of each candidate
        cand_periods: Best-fit period (s) of each candidate
        cand_bins: Number of phase bins in the folded profile for each candidate
        cand_codes: Per-file binary detection code string for each candidate
        codes_plot: Binary codes selected for plotting

    Returns:
        Tuple of (channels, periods, bins, codes) for candidates matching codes_plot
    """
    mask = np.isin(cand_codes, codes_plot)
    return cand_channels[mask], cand_periods[mask], cand_bins[mask], cand_codes[mask]

blipss.core.simulate_data

Atomic data-generation and signal-injection routines for synthetic filterbank simulation

generate_white_noise_background(n_channels, n_samples, rng=None)

Generate a 2-D array of Gaussian white noise with shape (n_channels, n_samples).

Parameters:

Name	Type	Description	Default
n_channels	int	Number of spectral channels.	required
n_samples	int	Number of time samples.	required
rng	Generator | None	Random number generator. Defaults to a fresh unseeded Generator.	None

Returns:

Type	Description
NDArray[float64]	Array of shape (n_channels, n_samples) drawn from N(0, 1).

Source code in blipss/core/simulate_data.py

def generate_white_noise_background(
    n_channels: int,
    n_samples: int,
    rng: np.random.Generator | None = None,
) -> npt.NDArray[np.float64]:
    """
    Generate a 2-D array of Gaussian white noise with shape (n_channels, n_samples).

    Args:
        n_channels: Number of spectral channels.
        n_samples: Number of time samples.
        rng: Random number generator. Defaults to a fresh unseeded Generator.

    Returns:
        Array of shape (n_channels, n_samples) drawn from N(0, 1).
    """
    if rng is None:
        rng = np.random.default_rng()
    return rng.standard_normal((n_channels, n_samples))

inject_periodic_signal(data, sample_times, channel, period, duty_cycle, pulse_snr, initial_phase)

Add a boxcar pulse train in-place to a single channel of the data array.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	Array of shape (n_channels, n_samples) to modify in place.	required
sample_times	NDArray[floating]	1-D array of sample timestamps (s).	required
channel	int	Index of the spectral channel to inject the signal into.	required
period	float	Pulse repetition period (s).	required
duty_cycle	float	Fraction of the period during which the pulse is on; in (0, 1].	required
pulse_snr	float	Peak signal-to-noise ratio added to on-pulse samples.	required
initial_phase	float	Phase offset of the pulse centre (fraction of a period); in [0, 1).	required

Source code in blipss/core/simulate_data.py

def inject_periodic_signal(
    data: npt.NDArray[np.floating],
    sample_times: npt.NDArray[np.floating],
    channel: int,
    period: float,
    duty_cycle: float,
    pulse_snr: float,
    initial_phase: float,
) -> None:
    """
    Add a boxcar pulse train in-place to a single channel of the data array.

    Args:
        data: Array of shape (n_channels, n_samples) to modify in place.
        sample_times: 1-D array of sample timestamps (s).
        channel: Index of the spectral channel to inject the signal into.
        period: Pulse repetition period (s).
        duty_cycle: Fraction of the period during which the pulse is on; in (0, 1].
        pulse_snr: Peak signal-to-noise ratio added to on-pulse samples.
        initial_phase: Phase offset of the pulse centre (fraction of a period); in [0, 1).
    """
    pulse_phase = sample_times / period % 1.0
    on_pulse_mask = (pulse_phase >= initial_phase - 0.5 * duty_cycle) & (pulse_phase < initial_phase + 0.5 * duty_cycle)
    data[channel, on_pulse_mask] += pulse_snr

reshape_for_sigproc(data)

Transpose and add an IF axis to match the sigproc (n_samples, n_ifs, n_channels) layout.

Parameters:

Name	Type	Description	Default
data	NDArray[floating]	Array of shape (n_channels, n_samples).	required

Returns:

Type	Description
NDArray[floating]	Array of shape (n_samples, 1, n_channels).

Source code in blipss/core/simulate_data.py

def reshape_for_sigproc(data: npt.NDArray[np.floating]) -> npt.NDArray[np.floating]:
    """
    Transpose and add an IF axis to match the sigproc (n_samples, n_ifs, n_channels) layout.

    Args:
        data: Array of shape (n_channels, n_samples).

    Returns:
        Array of shape (n_samples, 1, n_channels).
    """
    return np.expand_dims(data.T, axis=1)