Preprocessing

def to_differential(
    mats: List[np.ndarray],
    sr: Union[int, float],
    f: Dict[str, Union[int, float]],
) -> Tuple[List[np.ndarray], List[np.ndarray]]:
    """
    Exact Python replica of MATLAB `toDifferential`.

    Parameters
    ----------
    mats : list of ndarray
        `k` matrices, each shape (j, T) where j = number of monopolar
        EMG channels and T = samples.
    sr : float
        Sampling rate [Hz].
    f : dict
        Filter options – must contain keys
            'n'   : even Butterworth order (MATLAB style)
            'low' : lower cut-off frequency [Hz]
            'up'  : upper cut-off frequency [Hz]

    Returns
    -------
    dmat        : list of ndarray
        Filtered single-differential signals, each (j-1, T).
    dmatNoF     : list of ndarray
        Unfiltered single-differential signals, each (j-1, T).
    """

    # --------------- 1. sanity checks ----------
    if not mats:
        raise ValueError("`mats` must contain at least one matrix.")
    first_rows = mats[0].shape[0]
    if first_rows < 2:
        raise ValueError("Each matrix must have at least two rows (channels).")

    order = int(f["n"])
    fcl   = float(f["low"])
    fch   = float(f["up"])

    if order < 2 or order % 2:
        raise ValueError("Filter order 'n' must be an even integer ≥ 2.")

    dmat: List[np.ndarray]       = [None] * len(mats)
    dmatNoF: List[np.ndarray]    = [None] * len(mats)

    # --------------- 2. loop -----------------------
    for k, mat in enumerate(mats):
        mat = np.asarray(mat, dtype=np.float64)
        if mat.ndim != 2:
            raise ValueError("Each element of `mats` must be 2-D (channels × samples).")
        rows, T = mat.shape
        if rows != first_rows:
            raise ValueError("All matrices must have the same number of rows as mats[0].")

        diff_raw  = np.empty((rows - 1, T), dtype=np.float64)
        diff_filt = np.empty_like(diff_raw)

        for j in range(rows - 1):
            delta = mat[j + 1, :] - mat[j, :]
            diff_raw[j, :]  = delta
            diff_filt[j, :] = bandpass_filter(delta, order, fcl, fch, sr)

        dmat[k]       = diff_filt
        dmatNoF[k]    = diff_raw

    return dmat, dmatNoF

def emg_map_from_indices(emg_indices, rows, cols, order = 'column'):
    """
    EMG_MAP_FROM_INDICES builds an emg_map from a flat list of channel numbers

      emg_map = emg_map_from_indices(grid.emg_indices, grid.rows, grid.cols)

      INPUT
        emg_indices ... sequence of base-0 channel numbers, in acquisition
                         order, exactly rows*cols of them []
        rows        ... number of electrode rows of the grid []
        cols        ... number of electrode columns of the grid []
        order       ... char, 'column' (default) if the channels run down the
                         first column before moving to the next, 'row' if
                         they run along the first row []

      OUTPUT
        emg_map     ... 2-D array, cols-by-rows, base-0 channel numbers, ready
                         for map_to_columns []

      RAISES
        ValueError ... if rows or cols are not positive, if order is unknown,
                        or if len(emg_indices) != rows*cols

      *INFO* ... grids with an unwired position, e.g. HD08MM1305 with 64
                  electrodes in 5x13 = 65 positions, cannot be described this
                  way. Pass an explicit map with NaN at the empty position
                  instead - openhdemg's sort_rawemg and hdsemg-select's
                  electrode display grid both produce one
    """
    if rows < 1 or cols < 1:
        raise ValueError(f"rows and cols must be positive, got {rows} and {cols}.")
    if order not in ('column', 'row'):
        raise ValueError(f"order must be 'column' or 'row', got {order!r}.")

    idx = np.asarray(emg_indices, dtype=np.float64).ravel()
    if idx.size != rows * cols:
        raise ValueError(
            f"emg_indices holds {idx.size} channels but rows*cols is {rows * cols}. "
            "A grid with an unwired position needs an explicit map with NaN there."
        )

    if order == 'column':
        return idx.reshape(cols, rows)

    # Acquisition ran along the rows, so reshape row-wise and transpose
    return idx.reshape(rows, cols).T

def map_to_columns(emg_channels, emg_map):
    """
    MAP_TO_COLUMNS splits a channel matrix into one matrix per grid column

      columns = map_to_columns(emg.T, emg_map)

      INPUT
        emg_channels ... matrix, each row contains one channel's signal over
                          time (rows = channels, columns = samples). May hold
                          more channels than the map refers to, e.g. force or
                          path references [xV]
        emg_map      ... list of lists or 2-D array, outer index = grid
                          column, inner index = row, entries are base-0
                          channel numbers into emg_channels, NaN where no
                          electrode sits []

      OUTPUT
        columns      ... list of nCols matrices, each nRows-by-nSamples, in
                          map order. A NaN map entry becomes an all-NaN row
                          [xV]

      RAISES
        ValueError ... if the map is ragged or not 2-D, if an entry is not a
                        valid row index of emg_channels, or if a channel is
                        used more than once

      *INFO* ... EMGFile.data is samples-by-channels, so pass its transpose
      *INFO* ... excluding a channel means putting NaN in its map position;
                  global_amplitude then drops it, matching the MATLAB
                  '*INFO* NaN channels are ignored'
    """
    emg = np.asarray(emg_channels, dtype=np.float64)
    if emg.ndim != 2:
        raise ValueError(f"emg_channels must be 2-D (channels x samples), got {emg.ndim}-D.")

    grid = _as_map(emg_map)
    _check_indices(grid, emg.shape[0])

    # Build one matrix per column, NaN rows where the map has a gap
    columns = []
    for col in grid:
        block = np.full((col.size, emg.shape[1]), np.nan, dtype=np.float64)
        wired = ~np.isnan(col)
        block[wired] = emg[col[wired].astype(int)]
        columns.append(block)
    return columns