How to use yasa - 10 common examples

assert hypno.size == data.size, 'Hypno must have same size as data.'
        unique_hypno = np.unique(hypno)
        logger.info('Number of unique values in hypno = %i', unique_hypno.size)
        # Check include
        assert include is not None, 'include cannot be None if hypno is given'
        include = np.atleast_1d(np.asarray(include))
        assert include.size >= 1, '`include` must have at least one element.'
        assert hypno.dtype.kind == include.dtype.kind, ('hypno and include '
                                                        'must have same dtype')
        if not np.in1d(hypno, include).any():
            logger.error('None of the stages specified in `include` '
                         'are present in hypno. Returning None.')
            return None

    # Check data amplitude
    data_trimstd = trimbothstd(data, cut=0.10)
    data_ptp = np.ptp(data)
    logger.info('Number of samples in data = %i', data.size)
    logger.info('Sampling frequency = %.2f Hz', sf)
    logger.info('Data duration = %.2f seconds', data.size / sf)
    logger.info('Trimmed standard deviation of data = %.4f uV', data_trimstd)
    logger.info('Peak-to-peak amplitude of data = %.4f uV', data_ptp)
    if not(1 < data_trimstd < 1e3 or 1 < data_ptp < 1e6):
        logger.error('Wrong data amplitude. Unit must be uV. Returning None.')
        return None

    if 'rel_pow' not in thresh.keys():
        thresh['rel_pow'] = 0.20
    if 'corr' not in thresh.keys():
        thresh['corr'] = 0.65
    if 'rms' not in thresh.keys():
        thresh['rms'] = 1.5

# Hilbert power (to define the instantaneous frequency / power)
    n = data_sigma.size
    nfast = next_fast_len(n)
    analytic = signal.hilbert(data_sigma, N=nfast)[:n]
    inst_phase = np.angle(analytic)
    inst_pow = np.square(np.abs(analytic))
    # inst_freq = sf / 2pi * 1st-derivative of the phase of the analytic signal
    inst_freq = (sf / (2 * np.pi) * np.ediff1d(inst_phase))

    # Let's define the thresholds
    if hypno is None:
        thresh_rms = mrms.mean() + thresh['rms'] * trimbothstd(mrms, cut=0.10)
    else:
        thresh_rms = mrms[mask].mean() + thresh['rms'] * \
            trimbothstd(mrms[mask], cut=0.10)

    # Avoid too high threshold caused by Artefacts / Motion during Wake.
    thresh_rms = min(thresh_rms, 10)
    idx_rel_pow = (rel_pow >= thresh['rel_pow']).astype(int)
    idx_mcorr = (mcorr >= thresh['corr']).astype(int)
    idx_mrms = (mrms >= thresh_rms).astype(int)
    idx_sum = (idx_rel_pow + idx_mcorr + idx_mrms).astype(int)

    # Make sure that we do not detect spindles in REM or Wake if hypno != None
    if hypno is not None:
        idx_sum[~mask] = 0

    # For debugging
    logger.info('Moving RMS threshold = %.3f', thresh_rms)
    logger.info('Number of supra-theshold samples for relative power = %i',
                idx_rel_pow.sum())

logger.info('Number of unique values in hypno = %i', unique_hypno.size)
        # Check include
        assert include is not None, 'include cannot be None if hypno is given'
        include = np.atleast_1d(np.asarray(include))
        assert include.size >= 1, '`include` must have at least one element.'
        assert hypno.dtype.kind == include.dtype.kind, ('hypno and include '
                                                        'must have same dtype')
        if not np.in1d(hypno, include).any():
            logger.error('None of the stages specified in `include` '
                         'are present in hypno. Returning None.')
            return None

    # Check data amplitude
    # times = np.arange(data.size) / sf
    data = np.vstack((loc, roc))
    loc_trimstd = trimbothstd(loc, cut=0.10)
    roc_trimstd = trimbothstd(roc, cut=0.10)
    loc_ptp, roc_ptp = np.ptp(loc), np.ptp(roc)
    logger.info('Number of samples in data = %i', data.shape[1])
    logger.info('Original sampling frequency = %.2f Hz', sf)
    logger.info('Data duration = %.2f seconds', data.shape[1] / sf)
    logger.info('Trimmed standard deviation of LOC = %.4f uV', loc_trimstd)
    logger.info('Trimmed standard deviation of ROC = %.4f uV', roc_trimstd)
    logger.info('Peak-to-peak amplitude of LOC = %.4f uV', loc_ptp)
    logger.info('Peak-to-peak amplitude of ROC = %.4f uV', roc_ptp)
    if not(1 < loc_trimstd < 1e3 or 1 < loc_ptp < 1e6):
        logger.error('Wrong LOC amplitude. Unit must be uV. Returning None.')
        return None
    if not(1 < roc_trimstd < 1e3 or 1 < roc_ptp < 1e6):
        logger.error('Wrong ROC amplitude. Unit must be uV. Returning None.')
        return None

method='corr', interp=True)
    _, mrms = moving_transform(data_sigma, data, sf, window=.3, step=.1,
                               method='rms', interp=True)

    # Hilbert power (to define the instantaneous frequency / power)
    n = data_sigma.size
    nfast = next_fast_len(n)
    analytic = signal.hilbert(data_sigma, N=nfast)[:n]
    inst_phase = np.angle(analytic)
    inst_pow = np.square(np.abs(analytic))
    # inst_freq = sf / 2pi * 1st-derivative of the phase of the analytic signal
    inst_freq = (sf / (2 * np.pi) * np.ediff1d(inst_phase))

    # Let's define the thresholds
    if hypno is None:
        thresh_rms = mrms.mean() + thresh['rms'] * trimbothstd(mrms, cut=0.10)
    else:
        thresh_rms = mrms[mask].mean() + thresh['rms'] * \
            trimbothstd(mrms[mask], cut=0.10)

    # Avoid too high threshold caused by Artefacts / Motion during Wake.
    thresh_rms = min(thresh_rms, 10)
    idx_rel_pow = (rel_pow >= thresh['rel_pow']).astype(int)
    idx_mcorr = (mcorr >= thresh['corr']).astype(int)
    idx_mrms = (mrms >= thresh_rms).astype(int)
    idx_sum = (idx_rel_pow + idx_mcorr + idx_mrms).astype(int)

    # Make sure that we do not detect spindles in REM or Wake if hypno != None
    if hypno is not None:
        idx_sum[~mask] = 0

    # For debugging

assert hypno.size == data.size, 'Hypno must have same size as data.'
        unique_hypno = np.unique(hypno)
        logger.info('Number of unique values in hypno = %i', unique_hypno.size)
        # Check include
        assert include is not None, 'include cannot be None if hypno is given'
        include = np.atleast_1d(np.asarray(include))
        assert include.size >= 1, '`include` must have at least one element.'
        assert hypno.dtype.kind == include.dtype.kind, ('hypno and include '
                                                        'must have same dtype')
        if not np.in1d(hypno, include).any():
            logger.error('None of the stages specified in `include` '
                         'are present in hypno. Returning None.')
            return None

    # Check data amplitude
    data_trimstd = trimbothstd(data, cut=0.10)
    data_ptp = np.ptp(data)
    logger.info('Number of samples in data = %i', data.size)
    logger.info('Sampling frequency = %.2f Hz', sf)
    logger.info('Data duration = %.2f seconds', data.size / sf)
    logger.info('Trimmed standard deviation of data = %.4f uV', data_trimstd)
    logger.info('Peak-to-peak amplitude of data = %.4f uV', data_ptp)
    if not(1 < data_trimstd < 1e3 or 1 < data_ptp < 1e6):
        logger.error('Wrong data amplitude. Unit must be uV. Returning None.')
        return None

    # Check if we can downsample to 100 or 128 Hz
    if downsample is True and sf > 128:
        if sf % 100 == 0 or sf % 128 == 0:
            new_sf = 100 if sf % 100 == 0 else 128
            fac = int(sf / new_sf)
            sf = new_sf

# Check include
        assert include is not None, 'include cannot be None if hypno is given'
        include = np.atleast_1d(np.asarray(include))
        assert include.size >= 1, '`include` must have at least one element.'
        assert hypno.dtype.kind == include.dtype.kind, ('hypno and include '
                                                        'must have same dtype')
        if not np.in1d(hypno, include).any():
            logger.error('None of the stages specified in `include` '
                         'are present in hypno. Returning None.')
            return None

    # Check data amplitude
    # times = np.arange(data.size) / sf
    data = np.vstack((loc, roc))
    loc_trimstd = trimbothstd(loc, cut=0.10)
    roc_trimstd = trimbothstd(roc, cut=0.10)
    loc_ptp, roc_ptp = np.ptp(loc), np.ptp(roc)
    logger.info('Number of samples in data = %i', data.shape[1])
    logger.info('Original sampling frequency = %.2f Hz', sf)
    logger.info('Data duration = %.2f seconds', data.shape[1] / sf)
    logger.info('Trimmed standard deviation of LOC = %.4f uV', loc_trimstd)
    logger.info('Trimmed standard deviation of ROC = %.4f uV', roc_trimstd)
    logger.info('Peak-to-peak amplitude of LOC = %.4f uV', loc_ptp)
    logger.info('Peak-to-peak amplitude of ROC = %.4f uV', roc_ptp)
    if not(1 < loc_trimstd < 1e3 or 1 < loc_ptp < 1e6):
        logger.error('Wrong LOC amplitude. Unit must be uV. Returning None.')
        return None
    if not(1 < roc_trimstd < 1e3 or 1 < roc_ptp < 1e6):
        logger.error('Wrong ROC amplitude. Unit must be uV. Returning None.')
        return None

    # Check if we can downsample to 100 or 128 Hz

sp_rms = np.zeros(n_sp)
    sp_osc = np.zeros(n_sp)
    sp_sym = np.zeros(n_sp)
    sp_abs = np.zeros(n_sp)
    sp_rel = np.zeros(n_sp)
    sp_sta = np.zeros(n_sp)

    # Number of oscillations (= number of peaks separated by at least 60 ms)
    # --> 60 ms because 1000 ms / 16 Hz = 62.5 ms, in other words, at 16 Hz,
    # peaks are separated by 62.5 ms. At 11 Hz, peaks are separated by 90 ms.
    distance = 60 * sf / 1000

    for i in np.arange(len(sp))[good_dur]:
        # Important: detrend the signal to avoid wrong peak-to-peak amplitude
        sp_x = np.arange(data[sp[i]].size, dtype=np.float64)
        sp_det = _detrend(sp_x, data[sp[i]])
        # sp_det = signal.detrend(data[sp[i]], type='linear')
        sp_amp[i] = np.ptp(sp_det)  # Peak-to-peak amplitude
        sp_rms[i] = _rms(sp_det)  # Root mean square
        sp_rel[i] = np.median(rel_pow[sp[i]])  # Median relative power

        # Hilbert-based instantaneous properties
        sp_inst_freq = inst_freq[sp[i]]
        sp_inst_pow = inst_pow[sp[i]]
        sp_abs[i] = np.median(np.log10(sp_inst_pow[sp_inst_pow > 0]))
        sp_freq[i] = np.median(sp_inst_freq[sp_inst_freq > 0])

        # Number of oscillations
        peaks, peaks_params = signal.find_peaks(sp_det,
                                                distance=distance,
                                                prominence=(None, None))
        sp_osc[i] = len(peaks)

sp_abs = np.zeros(n_sp)
    sp_rel = np.zeros(n_sp)
    sp_sta = np.zeros(n_sp)

    # Number of oscillations (= number of peaks separated by at least 60 ms)
    # --> 60 ms because 1000 ms / 16 Hz = 62.5 ms, in other words, at 16 Hz,
    # peaks are separated by 62.5 ms. At 11 Hz, peaks are separated by 90 ms.
    distance = 60 * sf / 1000

    for i in np.arange(len(sp))[good_dur]:
        # Important: detrend the signal to avoid wrong peak-to-peak amplitude
        sp_x = np.arange(data[sp[i]].size, dtype=np.float64)
        sp_det = _detrend(sp_x, data[sp[i]])
        # sp_det = signal.detrend(data[sp[i]], type='linear')
        sp_amp[i] = np.ptp(sp_det)  # Peak-to-peak amplitude
        sp_rms[i] = _rms(sp_det)  # Root mean square
        sp_rel[i] = np.median(rel_pow[sp[i]])  # Median relative power

        # Hilbert-based instantaneous properties
        sp_inst_freq = inst_freq[sp[i]]
        sp_inst_pow = inst_pow[sp[i]]
        sp_abs[i] = np.median(np.log10(sp_inst_pow[sp_inst_pow > 0]))
        sp_freq[i] = np.median(sp_inst_freq[sp_inst_freq > 0])

        # Number of oscillations
        peaks, peaks_params = signal.find_peaks(sp_det,
                                                distance=distance,
                                                prominence=(None, None))
        sp_osc[i] = len(peaks)

        # For frequency and amplitude, we can also optionally use these
        # faster alternatives. If we use them, we do not need to compute the

# Now we compute the PTP amplitude and keep only the good peaks
    sw_ptp = np.abs(data_filt[idx_neg_peaks]) + data_filt[idx_pos_peaks]
    good_ptp = np.logical_and(sw_ptp > amp_ptp[0], sw_ptp < amp_ptp[1])

    # If good_ptp is all False
    if all(~good_ptp):
        logger.warning('No slow-wave with good amplitude. Returning None.')
        return None

    sw_ptp = sw_ptp[good_ptp]
    idx_neg_peaks = idx_neg_peaks[good_ptp]
    idx_pos_peaks = idx_pos_peaks[good_ptp]

    # Now we need to check the negative and positive phase duration
    # For that we need to compute the zero crossings of the filtered signal
    zero_crossings = _zerocrossings(data_filt)
    # Make sure that there is a zero-crossing after the last detected peak
    if zero_crossings[-1] < max(idx_pos_peaks[-1], idx_neg_peaks[-1]):
        # If not, append the index of the last peak
        zero_crossings = np.append(zero_crossings,
                                   max(idx_pos_peaks[-1], idx_neg_peaks[-1]))

    # Find distance to previous and following zc
    neg_sorted = np.searchsorted(zero_crossings, idx_neg_peaks)
    previous_neg_zc = zero_crossings[neg_sorted - 1] - idx_neg_peaks
    following_neg_zc = zero_crossings[neg_sorted] - idx_neg_peaks
    neg_phase_dur = (np.abs(previous_neg_zc) + following_neg_zc) / sf

    # Distance (in samples) between the positive peaks and the previous and
    # following zero-crossings
    pos_sorted = np.searchsorted(zero_crossings, idx_pos_peaks)
    previous_pos_zc = zero_crossings[pos_sorted - 1] - idx_pos_peaks

f, t, Sxx = stft_power(data, sf, window=2, step=.2, band=freq_broad,
                           interp=False, norm=True)
    idx_sigma = np.logical_and(f >= freq_sp[0], f <= freq_sp[1])
    rel_pow = Sxx[idx_sigma].sum(0)

    # Let's interpolate `rel_pow` to get one value per sample
    # Note that we could also have use the `interp=True` in the `stft_power`
    # function, however 2D interpolation is much slower than
    # 1D interpolation.
    func = interp1d(t, rel_pow, kind='cubic', bounds_error=False,
                    fill_value=0)
    t = np.arange(data.size) / sf
    rel_pow = func(t)

    # Now we apply moving RMS and correlation on the sigma-filtered signal
    _, mcorr = moving_transform(data_sigma, data, sf, window=.3, step=.1,
                                method='corr', interp=True)
    _, mrms = moving_transform(data_sigma, data, sf, window=.3, step=.1,
                               method='rms', interp=True)

    # Hilbert power (to define the instantaneous frequency / power)
    n = data_sigma.size
    nfast = next_fast_len(n)
    analytic = signal.hilbert(data_sigma, N=nfast)[:n]
    inst_phase = np.angle(analytic)
    inst_pow = np.square(np.abs(analytic))
    # inst_freq = sf / 2pi * 1st-derivative of the phase of the analytic signal
    inst_freq = (sf / (2 * np.pi) * np.ediff1d(inst_phase))

    # Let's define the thresholds
    if hypno is None:
        thresh_rms = mrms.mean() + thresh['rms'] * trimbothstd(mrms, cut=0.10)