Calculate significance by controlling for the false discovery rate.

    This function determines which of the p-values in `p` can be considered
    significant. Correction for multiple comparisons is performed by
    controlling the false discovery rate (FDR). The FDR is the maximum fraction
    of p-values that are wrongly considered significant [1]_.

    Parameters
    ----------
    p : array, shape (channels, channels, nfft)
        p-values.
    alpha : float
        Maximum false discovery rate.

    Returns
    -------
    s : array, dtype=bool, shape (channels, channels, nfft)
        Significance of each p-value.

    References
    ----------
    .. [1] Y. Benjamini, Y. Hochberg. Controlling the false discovery rate: a
           practical and powerful approach to multiple testing. J. Royal Stat.
           Soc. Series B 57(1): 289-300, 1995.

def significance_fdr(p, alpha):
    """Calculate significance by controlling for the false discovery rate.

    This function determines which of the p-values in `p` can be considered
    significant. Correction for multiple comparisons is performed by
    controlling the false discovery rate (FDR). The FDR is the maximum fraction
    of p-values that are wrongly considered significant [1]_.

    Parameters
    ----------
    p : array, shape (channels, channels, nfft)
        p-values.
    alpha : float
        Maximum false discovery rate.

    Returns
    -------
    s : array, dtype=bool, shape (channels, channels, nfft)
        Significance of each p-value.

    References
    ----------
    .. [1] Y. Benjamini, Y. Hochberg. Controlling the false discovery rate: a
           practical and powerful approach to multiple testing. J. Royal Stat.
           Soc. Series B 57(1): 289-300, 1995.
    """
    i = np.argsort(p, axis=None)
    m = i.size - np.sum(np.isnan(p))

    j = np.empty(p.shape, int)
    j.flat[i] = np.arange(1, i.size + 1)

    mask = p <= alpha * j / m

    if np.sum(mask) == 0:
        return mask

    # find largest k so that p_k <= alpha*k/m
    k = np.max(j[mask])

    # reject all H_i for i = 0...k
    s = j <= k

    return s