def detect_lines(s, tapers, p=None, **taper_kws):
    from nitime.algorithms import tapered_spectra, dpss_windows
    import scipy.stats.distributions as dists
    import scipy.ndimage as ndimage
    N = s.shape[-1]
    # Some boiler-plate --
    # 1) set up tapers
    # 2) perform FFT on all windowed series
    if not isinstance(tapers, np.ndarray):
        # then tapers is (NW, K)
        args = (N,) + tuple(tapers)
        dpss, eigvals = dpss_windows(*args)
        if taper_kws.pop('low_bias', False):
            keepers = (eigvals > 0.9)
            dpss = dpss[keepers]
        tapers = dpss
    # spectra is (n_arr, K, nfft)
    spectra = tapered_spectra(s, tapers, **taper_kws)
    nfft = spectra.shape[-1]
    spectra = spectra[..., :nfft//2 + 1]

    # Set up some data for the following calculations --
    #   get the DC component of the taper spectra
    K = tapers.shape[0]
    U0 = tapers.sum(axis=1)
    U_sq = np.sum(U0**2)
    #  first order linear regression for mu to explain spectra
    mu = np.sum( U0[:,None] * spectra, axis=-2 ) / U_sq

    # numerator of F-stat -- strength of regression
    numr = 0.5 * np.abs(mu)**2 * U_sq
    numr[...,0] = 1; # don't care about DC
    # denominator -- strength of residual
    spectra = np.rollaxis(spectra, -2, 0)
    U0.shape = (K,) + (1,) * (spectra.ndim-1)
    denomr = spectra - U0*mu
    denomr = np.sum(np.abs(denomr)**2, axis=0) / (2*K-2)
    denomr[...,0] = 1;
    f_stat = numr / denomr

    # look for lines in each F-spectrum
    if not p:
        # the number of simultaneous tests are nfft/2, so this puts
        # the expected value for false detection somewhere less than 1
        p = 1.0/nfft
    #thresh = dists.f.isf(p, 2, 2*K-2)
    thresh = dists.f.isf(p, 2, K-1)
    f_stat = np.atleast_2d(f_stat)
    mu = np.atleast_2d(mu)
    lines = ()
    for fs, m in zip(f_stat, mu):
        detected = np.where(fs > thresh)[0]
        # do a quick pass through the detected lines to reject multiple
        # hits within the 2NW resolution of the MT analysis -- approximate
        # 2NW by K
        ddiff = np.diff(detected)
        flagged_groups, last_group = ndimage.label( (ddiff < K) )
        for g in range(1,last_group+1):
            idx = np.where(flagged_groups==g)[0]
            idx = np.r_[idx, idx[-1]+1]
            # keep the super-threshold point with largest amplitude
            mx = np.argmax(np.abs(m[ detected[idx] ]))
            i_sv = detected[idx[mx]]
            detected[idx] = -1
            detected[idx[mx]] = i_sv
        detected = detected[detected>0]
        if len(detected):
            lines = lines + ( (detected/float(nfft), m[detected]), )
        else:
            lines = lines + ( (), )
    if len(lines) == 1:
        lines = lines[0]
    return lines


#-----------------------------------------------------------------------------
# Eigensystem utils
#-----------------------------------------------------------------------------

# If we can get it, we want the cythonized version
Detect the presence of line spectra in s using the F-test
dataset/ETHPy150Open nipy/nitime/nitime/utils.py/detect_lines