def detect_Lacourse2018(dat_orig, s_freq, time, opts):
    """Spindle detection based on Lacourse et al., 2018
    
    Parameters
    ----------
    dat_orig : ndarray (dtype='float')
        vector with the data for one channel
    s_freq : float
        sampling frequency
    time : ndarray (dtype='float')
        vector with the time points for each sample
    opts : instance of 'DetectSpindle'
    
    Returns
    -------
    list of dict
        list of detected spindles
    dict
        'det_value_lo' with detection value, 'det_value_hi' is nan,
        'sel_value' is nan (for consistency with other methods)
    float
        spindle density, per 30-s epoch

    References
    ----------
    Lacourse, K. et al. J. Neurosci. Meth. (2018).
    """
    # Downsample z-score parameters, tolerance
    step = opts.windowing['step']
    if step:
        ds_freq = int(1 / step) # downsampled sampling frequency
        opts.zscore['dur'] *= opts.windowing['step']
        opts.tolerance *= opts.windowing['step']
        if opts.zscore['step']:
            opts.zscore['step'] *= opts.windowing['step']
    else:
        ds_freq = s_freq

    
    # Absolute sigma power
    dat_sigma = transform_signal(dat_orig, s_freq, 'double_sosbutter', 
                                 opts.det_butter)
    dat_det = transform_signal(dat_sigma, s_freq, 'moving_ms', opts.moving_ms)
    dat_det[dat_det <= 0] = 0.000000001
    abs_sig_pow = log10(dat_det)
    # Option to adapt the absolute threshold, for low-amplitude recordings
    if opts.abs_pow_thresh < 0:
        opts.abs_pow_thresh = (mean(abs_sig_pow) - 
                               opts.abs_pow_thresh * std(abs_sig_pow))
    abs_sig_pow = transform_signal(abs_sig_pow, ds_freq, 'smooth', opts.smooth)
    
    # Relative sigma power
    dat_det = transform_signal(dat_orig, s_freq, 'moving_power_ratio', 
                               opts.moving_power_ratio)
    dat_det[dat_det <= 0] = 0.000000001
    dat_det = log10(dat_det)
    rel_sig_pow = transform_signal(dat_det, s_freq, 'moving_zscore', 
                                   opts.zscore)
    rel_sig_pow = transform_signal(rel_sig_pow, ds_freq, 'smooth', opts.smooth)
    
    # Sigma covariance
    dat_broad = transform_signal(dat_orig, s_freq, 'double_sosbutter', 
                                 opts.det_butter2)
    dat_covar = transform_signal(dat_sigma, s_freq, 'moving_covar', 
                                 opts.moving_covar, dat2=dat_broad)
    dat_det = dat_covar.copy()
    dat_det[dat_det < 0] = 0 # negative covariances are discarded
    dat_det = log10(dat_det + 1) # add 1 to avoid -inf
    sigma_covar = transform_signal(dat_det, s_freq, 'moving_zscore', 
                                   opts.zscore)
    sigma_covar = transform_signal(sigma_covar, ds_freq, 'smooth', opts.smooth)
    
    # Sigma correlation
    dat_sd_broad = transform_signal(dat_broad, s_freq, 'moving_sd', 
                                    opts.moving_sd)
    dat_sd_sigma = transform_signal(dat_sigma, s_freq, 'moving_sd', 
                                    opts.moving_sd)
    dat_sd_broad[dat_sd_broad == 0] = 0.000000001
    dat_sd_sigma[dat_sd_sigma == 0] = 0.000000001
    sigma_corr = dat_covar / (dat_sd_broad * dat_sd_sigma)
    sigma_corr = transform_signal(sigma_corr, ds_freq, 'smooth', opts.smooth)

    # Thresholding
    abs_and_cov = logical_and(abs_sig_pow >= opts.abs_pow_thresh,
                              sigma_covar >= opts.covar_thresh)
    concensus = logical_and.reduce((rel_sig_pow >= opts.rel_pow_thresh,
                                    sigma_corr >= opts.corr_thresh,
                                    abs_and_cov))                                    
    events = detect_events(concensus, 'custom') # at s_freq * 0.1
    
    if events is not None:
        events = _merge_close(dat_sigma, events, time, opts.tolerance)
        events = _select_period(events, abs_and_cov)
        
        if opts.windowing['step']:
            events = events * (s_freq * opts.windowing['step']) # upsample
            events = asarray(around(events), dtype=int)
        
        events = within_duration(events, time, opts.duration)
        events = _merge_close(dat_sigma, events, time, opts.min_interval)
        events = remove_straddlers(events, time, s_freq)

        power_peaks = peak_in_power(events, dat_orig, s_freq, opts.power_peaks)
        powers = power_in_band(events, dat_orig, s_freq, opts.frequency)
        sp_in_chan = make_spindles(events, power_peaks, powers, dat_sigma,
                                   dat_orig, time, s_freq)

    else:
        lg.info('No spindle found')
        sp_in_chan = []

    values = {'abs_pow_thresh': opts.abs_pow_thresh, 
              'rel_pow_thresh': opts.rel_pow_thresh, 
              'covar_thresh': opts.covar_thresh,
              'corr_thresh': opts.corr_thresh}

    density = len(sp_in_chan) * s_freq * 30 / len(dat_orig)

    return sp_in_chan, values, density
def detect_Lacourse2018 ( dat_orig , s_freq , time , opts ) : step = opts . windowing [ 'step' ] if step : ds_freq = int ( 1 / step ) opts . zscore [ 'dur' ] *= opts . windowing [ 'step' ] opts . tolerance *= opts . windowing [ 'step' ] if opts . zscore [ 'step' ] : opts . zscore [ 'step' ] *= opts . windowing [ 'step' ] else : ds_freq = s_freq dat_sigma = transform_signal ( dat_orig , s_freq , 'double_sosbutter' , opts . det_butter ) dat_det = transform_signal ( dat_sigma , s_freq , 'moving_ms' , opts . moving_ms ) dat_det [ dat_det <= 0 ] = 0.000000001 abs_sig_pow = log10 ( dat_det ) if opts . abs_pow_thresh < 0 : opts . abs_pow_thresh = ( mean ( abs_sig_pow ) - opts . abs_pow_thresh * std ( abs_sig_pow ) ) abs_sig_pow = transform_signal ( abs_sig_pow , ds_freq , 'smooth' , opts . smooth ) dat_det = transform_signal ( dat_orig , s_freq , 'moving_power_ratio' , opts . moving_power_ratio ) dat_det [ dat_det <= 0 ] = 0.000000001 dat_det = log10 ( dat_det ) rel_sig_pow = transform_signal ( dat_det , s_freq , 'moving_zscore' , opts . zscore ) rel_sig_pow = transform_signal ( rel_sig_pow , ds_freq , 'smooth' , opts . smooth ) dat_broad = transform_signal ( dat_orig , s_freq , 'double_sosbutter' , opts . det_butter2 ) dat_covar = transform_signal ( dat_sigma , s_freq , 'moving_covar' , opts . moving_covar , dat2 = dat_broad ) dat_det = dat_covar . copy ( ) dat_det [ dat_det < 0 ] = 0 dat_det = log10 ( dat_det + 1 ) sigma_covar = transform_signal ( dat_det , s_freq , 'moving_zscore' , opts . zscore ) sigma_covar = transform_signal ( sigma_covar , ds_freq , 'smooth' , opts . smooth ) dat_sd_broad = transform_signal ( dat_broad , s_freq , 'moving_sd' , opts . moving_sd ) dat_sd_sigma = transform_signal ( dat_sigma , s_freq , 'moving_sd' , opts . moving_sd ) dat_sd_broad [ dat_sd_broad == 0 ] = 0.000000001 dat_sd_sigma [ dat_sd_sigma == 0 ] = 0.000000001 sigma_corr = dat_covar / ( dat_sd_broad * dat_sd_sigma ) sigma_corr = transform_signal ( sigma_corr , ds_freq , 'smooth' , opts . smooth ) abs_and_cov = logical_and ( abs_sig_pow >= opts . abs_pow_thresh , sigma_covar >= opts . covar_thresh ) concensus = logical_and . reduce ( ( rel_sig_pow >= opts . rel_pow_thresh , sigma_corr >= opts . corr_thresh , abs_and_cov ) ) events = detect_events ( concensus , 'custom' ) if events is not None : events = _merge_close ( dat_sigma , events , time , opts . tolerance ) events = _select_period ( events , abs_and_cov ) if opts . windowing [ 'step' ] : events = events * ( s_freq * opts . windowing [ 'step' ] ) events = asarray ( around ( events ) , dtype = int ) events = within_duration ( events , time , opts . duration ) events = _merge_close ( dat_sigma , events , time , opts . min_interval ) events = remove_straddlers ( events , time , s_freq ) power_peaks = peak_in_power ( events , dat_orig , s_freq , opts . power_peaks ) powers = power_in_band ( events , dat_orig , s_freq , opts . frequency ) sp_in_chan = make_spindles ( events , power_peaks , powers , dat_sigma , dat_orig , time , s_freq ) else : lg . info ( 'No spindle found' ) sp_in_chan = [ ] values = { 'abs_pow_thresh' : opts . abs_pow_thresh , 'rel_pow_thresh' : opts . rel_pow_thresh , 'covar_thresh' : opts . covar_thresh , 'corr_thresh' : opts . corr_thresh } density = len ( sp_in_chan ) * s_freq * 30 / len ( dat_orig ) return sp_in_chan , values , density
Spindle detection based on Lacourse et al., 2018
    
    Parameters
    ----------
    dat_orig : ndarray (dtype='float')
        vector with the data for one channel
    s_freq : float
        sampling frequency
    time : ndarray (dtype='float')
        vector with the time points for each sample
    opts : instance of 'DetectSpindle'
    
    Returns
    -------
    list of dict
        list of detected spindles
    dict
        'det_value_lo' with detection value, 'det_value_hi' is nan,
        'sel_value' is nan (for consistency with other methods)
    float
        spindle density, per 30-s epoch

    References
    ----------
    Lacourse, K. et al. J. Neurosci. Meth. (2018).