Plain python implementation of the `scipy.signal.find_peaks` function.

find_peaks.py

def local_maxima_1D(x: list[float]) -> tuple[list[int], list[int], list[int]]:
    """
    Find local maxima in a 1D array.
    This function finds all local maxima in a 1D array and returns the indices
    for their edges and midpoints (rounded down for even plateau sizes).

    Parameters
    ----------
    x : list[float]
        The array to search for local maxima.

    Returns
    -------
    midpoints : list[int]
        Indices of midpoints of local maxima in `x`.
    left_edges : list[int]
        Indices of edges to the left of local maxima in `x`.
    right_edges : list[int]
        Indices of edges to the right of local maxima in `x`.

    Notes
    - A maxima is defined as one or more samples of equal value that are
      surrounded on both sides by at least one smaller sample.
    """
    midpoints = [0] * (len(x) // 2)
    left_edges = [0] * (len(x) // 2)
    right_edges = [0] * (len(x) // 2)
    m = 0  # Pointer to the end of valid area in allocated arrays
    i = 1  # Pointer to current sample, first one can't be maxima
    i_max = len(x) - 1  # Last sample can't be maxima
    while i < i_max:
        # Test if previous sample is smaller
        if x[i - 1] < x[i]:
            i_ahead = i + 1  # Index to look ahead of current sample

            # Find next sample that is unequal to x[i]
            while i_ahead < i_max and x[i_ahead] == x[i]:
                i_ahead += 1

            # Maxima is found if next unequal sample is smaller than x[i]
            if x[i_ahead] < x[i]:
                left_edges[m] = i
                right_edges[m] = i_ahead - 1
                midpoints[m] = (left_edges[m] + right_edges[m]) // 2
                m += 1
                # Skip samples that can't be maximum
                i = i_ahead
        i += 1

    # Keep only valid part of array memory.
    midpoints = midpoints[:m]
    left_edges = left_edges[:m]
    right_edges = right_edges[:m]

    return midpoints, left_edges, right_edges


def peak_prominences(
    x: list[float], peaks: list[int], wlen: int | None = None
) -> tuple[list[float], list[int], list[int]]:
    """
    Calculate the prominence of each peak in a signal.

    Parameters
    x: list[float] A signal with peaks.
    peaks: list[int] Indices of peaks in `x`.
    wlen: int A window length in samples (see `peak_prominences`) which is rounded up
        to the nearest odd integer. If smaller than 2 the entire signal `x` is
        used.

    Returns
    prominences : list[float]
        The calculated prominences for each peak in `peaks`.
    left_bases, right_bases : list[int]
        The peaks' bases as indices in `x` to the left and right of each peak.

    Raises
    ValueError
        If a value in `peaks` is an invalid index for `x`.

    Warns
    PeakPropertyWarning
        If a prominence of 0 was calculated for any peak.
    """
    if wlen is None:
        wlen = -1
    elif wlen <= 1:
        raise ValueError("wlen must be larger than 1")
    len_peaks = len(peaks)
    prominences = [0.0] * len_peaks
    left_bases = [0] * len_peaks
    right_bases = [0] * len_peaks

    for peak_nr in range(len_peaks):
        peak = peaks[peak_nr]
        i_min = 0
        i_max = len(x) - 1
        if not i_min <= peak <= i_max:
            raise ValueError(f"peak {peak} is not a valid index for `x`")

        if 2 <= wlen:
            # Adjust window around the evaluated peak (within bounds);
            # if wlen is even the resulting window length is implicitly
            # rounded to next odd integer
            i_min = max(peak - wlen // 2, i_min)
            i_max = min(peak + wlen // 2, i_max)

        # Find the left base in interval [i_min, peak]
        i = left_bases[peak_nr] = peak
        left_min = x[peak]
        while i_min <= i and x[i] <= x[peak]:
            if x[i] < left_min:
                left_min = x[i]
                left_bases[peak_nr] = i
            i -= 1

        # Find the right base in interval [peak, i_max]
        i = right_bases[peak_nr] = peak
        right_min = x[peak]
        while i <= i_max and x[i] <= x[peak]:
            if x[i] < right_min:
                right_min = x[i]
                right_bases[peak_nr] = i
            i += 1

        prominences[peak_nr] = x[peak] - max(left_min, right_min)
        if prominences[peak_nr] == 0:
            print(f"Prominence of 0 calculated for peak at index {peak}.")
    return prominences, left_bases, right_bases


# Create vanilla Python find_peaks function
def find_peaks(x: list[float], prominence: float = 0.001) -> list[int]:
    """Find peaks in a 1D array, given a prominence value."""
    peaks, _, _ = local_maxima_1D(x)
    prominences, _, _ = peak_prominences(x, peaks)
    # keep the peaks that have a prominence larger than the threshold
    peaks = [p for p, prom in zip(peaks, prominences) if prom > prominence]
    return peaks