Aluode/PerceptionLabPortable

1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import numpy as np6 7from ..utils import _pl, logger, verbose8 9 10@verbose11def peak_finder(x0, thresh=None, extrema=1, verbose=None):12    """Noise-tolerant fast peak-finding algorithm.13 14    Parameters15    ----------16    x0 : 1d array17        A real vector from the maxima will be found (required).18    thresh : float | None19        The amount above surrounding data for a peak to be20        identified. Larger values mean the algorithm is more selective in21        finding peaks. If ``None``, use the default of22        ``(max(x0) - min(x0)) / 4``.23    extrema : {-1, 1}24        1 if maxima are desired, -1 if minima are desired25        (default = maxima, 1).26    %(verbose)s27 28    Returns29    -------30    peak_loc : array31        The indices of the identified peaks in x0.32    peak_mag : array33        The magnitude of the identified peaks.34 35    Notes36    -----37    If repeated values are found the first is identified as the peak.38    Conversion from initial Matlab code from:39    Nathanael C. Yoder (ncyoder@purdue.edu)40 41    Examples42    --------43    >>> import numpy as np44    >>> from mne.preprocessing import peak_finder45    >>> t = np.arange(0, 3, 0.01)46    >>> x = np.sin(np.pi*t) - np.sin(0.5*np.pi*t)47    >>> peak_locs, peak_mags = peak_finder(x) # doctest: +SKIP48    >>> peak_locs # doctest: +SKIP49    array([36, 260]) # doctest: +SKIP50    >>> peak_mags # doctest: +SKIP51    array([0.36900026, 1.76007351]) # doctest: +SKIP52    """53    x0 = np.asanyarray(x0)54    s = x0.size55 56    if x0.ndim >= 2 or s == 0:57        raise ValueError("The input data must be a non empty 1D vector")58 59    if thresh is None:60        thresh = (np.max(x0) - np.min(x0)) / 461        logger.debug(f"Peak finder automatic threshold: {thresh:0.2g}")62 63    assert extrema in [-1, 1]64 65    if extrema == -1:66        x0 = extrema * x0  # Make it so we are finding maxima regardless67 68    dx0 = np.diff(x0)  # Find derivative69    # This is so we find the first of repeated values70    dx0[dx0 == 0] = -np.finfo(float).eps71    # Find where the derivative changes sign72    ind = np.where(dx0[:-1:] * dx0[1::] < 0)[0] + 173 74    # Include endpoints in potential peaks and valleys75    x = np.concatenate((x0[:1], x0[ind], x0[-1:]))76    ind = np.concatenate(([0], ind, [s - 1]))77    del x078 79    #  x only has the peaks, valleys, and endpoints80    length = x.size81    min_mag = np.min(x)82 83    if length > 2:  # Function with peaks and valleys84        # Set initial parameters for loop85        temp_mag = min_mag86        found_peak = False87        left_min = min_mag88 89        # Deal with first point a little differently since tacked it on90        # Calculate the sign of the derivative since we took the first point91        # on it does not necessarily alternate like the rest.92        signDx = np.sign(np.diff(x[:3]))93        if signDx[0] <= 0:  # The first point is larger or equal to the second94            ii = -195            if signDx[0] == signDx[1]:  # Want alternating signs96                x = np.concatenate((x[:1], x[2:]))97                ind = np.concatenate((ind[:1], ind[2:]))98                length -= 199 100        else:  # First point is smaller than the second101            ii = 0102            if signDx[0] == signDx[1]:  # Want alternating signs103                x = x[1:]104                ind = ind[1:]105                length -= 1106 107        # Preallocate max number of maxima108        maxPeaks = int(np.ceil(length / 2.0))109        peak_loc = np.zeros(maxPeaks, dtype=np.int64)110        peak_mag = np.zeros(maxPeaks)111        c_ind = 0112        # Loop through extrema which should be peaks and then valleys113        while ii < (length - 1):114            ii += 1  # This is a peak115            # Reset peak finding if we had a peak and the next peak is bigger116            # than the last or the left min was small enough to reset.117            if found_peak and (118                (x[ii] > peak_mag[-1]) or (left_min < peak_mag[-1] - thresh)119            ):120                temp_mag = min_mag121                found_peak = False122 123            # Make sure we don't iterate past the length of our vector124            if ii == length - 1:125                break  # We assign the last point differently out of the loop126 127            # Found new peak that was lager than temp mag and threshold larger128            # than the minimum to its left.129            if (x[ii] > temp_mag) and (x[ii] > left_min + thresh):130                temp_loc = ii131                temp_mag = x[ii]132 133            ii += 1  # Move onto the valley134            # Come down at least thresh from peak135            if not found_peak and (temp_mag > (thresh + x[ii])):136                found_peak = True  # We have found a peak137                left_min = x[ii]138                peak_loc[c_ind] = temp_loc  # Add peak to index139                peak_mag[c_ind] = temp_mag140                c_ind += 1141            elif x[ii] < left_min:  # New left minima142                left_min = x[ii]143 144        # Check end point145        if (x[-1] > temp_mag) and (x[-1] > (left_min + thresh)):146            peak_loc[c_ind] = length - 1147            peak_mag[c_ind] = x[-1]148            c_ind += 1149        elif not found_peak and temp_mag > min_mag:150            # Check if we still need to add the last point151            peak_loc[c_ind] = temp_loc152            peak_mag[c_ind] = temp_mag153            c_ind += 1154 155        # Create output156        peak_inds = ind[peak_loc[:c_ind]]157        peak_mags = peak_mag[:c_ind]158    else:  # This is a monotone function where an endpoint is the only peak159        x_ind = np.argmax(x)160        peak_mags = x[x_ind]161        if peak_mags > (min_mag + thresh):162            peak_inds = ind[x_ind]163        else:164            peak_mags = []165            peak_inds = []166 167    # Change sign of data if was finding minima168    if extrema < 0:169        peak_mags *= -1.0170 171    # ensure output type array172    if not isinstance(peak_inds, np.ndarray):173        peak_inds = np.atleast_1d(peak_inds).astype("int64")174 175    if not isinstance(peak_mags, np.ndarray):176        peak_mags = np.atleast_1d(peak_mags).astype("float64")177 178    # Plot if no output desired179    if len(peak_inds) == 0:180        logger.info("No significant peaks found")181    else:182        logger.info(f"Found {len(peak_inds)} significant peak{_pl(peak_inds)}")183 184    return peak_inds, peak_mags185