def xcorr(self):
        """The cross-correlation between every pairwise combination time-series
        in the object. Uses np.correlation('full').

        Returns
        -------

        TimeSeries : the time-dependent cross-correlation, with zero-lag
        at time=0

        """
        tseries_length = self.input.data.shape[0]
        t_points = self.input.data.shape[-1]
        xcorr = np.zeros((tseries_length,
                          tseries_length,
                          t_points * 2 - 1))
        data = self.input.data
        for i in range(tseries_length):
            data_i = data[i]
            for j in range(i, tseries_length):
                xcorr[i, j] = np.correlate(data_i,
                                          data[j],
                                          mode='full')

        idx = tril_indices(tseries_length, -1)
        xcorr[idx[0], idx[1], ...] = xcorr[idx[1], idx[0], ...]

        return ts.TimeSeries(xcorr,
                             sampling_interval=self.input.sampling_interval,
                             t0=-self.input.sampling_interval * t_points)

    def xcorr_norm(self):
        """The cross-correlation between every pairwise combination time-series
        in the object, where the zero lag correlation is normalized to be equal
        to the correlation coefficient between the time-series

        Returns
        -------

        TimeSeries : A TimeSeries object
            the time-dependent cross-correlation, with zero-lag at time=0

        """

        tseries_length = self.input.data.shape[0]
        t_points = self.input.data.shape[-1]
        xcorr = np.zeros((tseries_length,
                          tseries_length,
                          t_points * 2 - 1))
        data = self.input.data
        for i in range(tseries_length):
            data_i = data[i]
            for j in range(i, tseries_length):
                xcorr[i, j] = np.correlate(data_i,
                                          data[j],
                                          mode='full')
                xcorr[i, j] /= (xcorr[i, j, t_points])
                xcorr[i, j] *= self.corrcoef[i, j]

        idx = tril_indices(tseries_length, -1)
        xcorr[idx[0], idx[1], ...] = xcorr[idx[1], idx[0], ...]

        return ts.TimeSeries(xcorr,
                             sampling_interval=self.input.sampling_interval,
                             t0=-self.input.sampling_interval * t_points)

    def coherence_partial(self):
        """The partial coherence between data[i] and data[j], given data[k], as
        a function of frequency band"""

        tseries_length = self.input.data.shape[0]
        spectrum_length = self.spectrum.shape[-1]

        p_coherence = np.zeros((tseries_length,
                                tseries_length,
                                tseries_length,
                                spectrum_length))

        for i in range(tseries_length):
            for j in range(tseries_length):
                for k in range(tseries_length):
                    if j == k or i == k:
                        pass
                    else:
                        p_coherence[i][j][k] = tsa.coherence_partial_spec(
                            self.spectrum[i][j],
                            self.spectrum[i][i],
                            self.spectrum[j][j],
                            self.spectrum[i][k],
                            self.spectrum[j][k],
                            self.spectrum[k][k])

        idx = tril_indices(tseries_length, -1)
        p_coherence[idx[0], idx[1], ...] =\
                            p_coherence[idx[1], idx[0], ...].conj()

        return p_coherence

def coherency(time_series, csd_method=None):
    r"""
    Compute the coherency between the spectra of n-tuple of time series.
    Input to this function is in the time domain

    Parameters
    ----------

    time_series : n*t float array
       an array of n different time series of length t each

    csd_method : dict, optional.
       See :func:`get_spectra` documentation for details

    Returns
    -------

    f : float array
        The central frequencies for the frequency bands for which the spectra
        are estimated

    c : float array
        This is a symmetric matrix with the coherencys of the signals. The
        coherency of signal i and signal j is in f[i][j]. Note that f[i][j] =
        f[j][i].conj()

    Notes
    -----

    This is an implementation of equation (1) of Sun (2005):

    .. math::

        R_{xy} (\lambda) = \frac{f_{xy}(\lambda)}
        {\sqrt{f_{xx} (\lambda) \cdot f_{yy}(\lambda)}}

    F.T. Sun and L.M. Miller and M. D'Esposito (2005). Measuring temporal
    dynamics of functional networks using phase spectrum of fMRI
    data. Neuroimage, 28: 227-37.

    """
    if csd_method is None:
        csd_method = {'this_method': 'welch'}  # The default

    f, fxy = get_spectra(time_series, csd_method)

    #A container for the coherencys, with the size and shape of the expected
    #output:
    c = np.zeros((time_series.shape[0],
                  time_series.shape[0],
                  f.shape[0]), dtype=complex)  # Make sure it's complex

    for i in range(time_series.shape[0]):
        for j in range(i, time_series.shape[0]):
            c[i][j] = coherency_spec(fxy[i][j], fxy[i][i], fxy[j][j])

    idx = tril_indices(time_series.shape[0], -1)
    c[idx[0], idx[1], ...] = c[idx[1], idx[0], ...].conj()  # Make it symmetric

    return f, c

    def coherence(self):
        """
        The coherence between the different channels in the input TimeSeries
        object
        """

        #XXX Calculate this from the standard output, instead of recalculating
        #the coherence:

        tseries_length = self.input.data.shape[0]
        spectrum_length = self.spectrum.shape[-1]
        coherence = np.zeros((tseries_length,
                              tseries_length,
                              spectrum_length))

        for i in range(tseries_length):
            for j in range(i, tseries_length):
                coherence[i][j] = tsa.coherence_spec(self.spectrum[i][j],
                                                     self.spectrum[i][i],
                                                     self.spectrum[j][j])

        idx = tril_indices(tseries_length, -1)
        coherence[idx[0], idx[1], ...] = coherence[idx[1], idx[0], ...].conj()

        return coherence

def coherence(time_series, csd_method=None):
    r"""Compute the coherence between the spectra of an n-tuple of time_series.

    Parameters of this function are in the time domain.

    Parameters
    ----------
    time_series: float array
       an array of different time series with time as the last dimension

    csd_method: dict, optional
       See :func:`algorithms.spectral.get_spectra` documentation for details

    Returns
    -------
    f : float array
        The central frequencies for the frequency bands for which the spectra
        are estimated

    c : float array
        This is a symmetric matrix with the coherencys of the signals. The
        coherency of signal i and signal j is in f[i][j].

    Notes
    -----

    This is an implementation of equation (2) of Sun (2005):

    .. math::

        Coh_{xy}(\lambda) = |{R_{xy}(\lambda)}|^2 =
        \frac{|{f_{xy}(\lambda)}|^2}{f_{xx}(\lambda) \cdot f_{yy}(\lambda)}

    F.T. Sun and L.M. Miller and M. D'Esposito (2005). Measuring temporal
    dynamics of functional networks using phase spectrum of fMRI data.
    Neuroimage, 28: 227-37.

    """
    if csd_method is None:
        csd_method = {'this_method': 'welch'}  # The default

    f, fxy = get_spectra(time_series, csd_method)

    # A container for the coherences, with the size and shape of the expected
    # output:
    c = np.zeros((time_series.shape[0],
                  time_series.shape[0],
                  f.shape[0]))

    for i in xrange(time_series.shape[0]):
        for j in xrange(i, time_series.shape[0]):
            c[i][j] = coherence_spec(fxy[i][j], fxy[i][i], fxy[j][j])

    idx = tril_indices(time_series.shape[0], -1)
    c[idx[0], idx[1], ...] = c[idx[1], idx[0], ...].conj()  # Make it symmetric

    return f, c

def coherence_bavg(time_series, lb=0, ub=None, csd_method=None):
    r"""
    Compute the band-averaged coherence between the spectra of two time series.

    Input to this function is in the time domain.

    Parameters
    ----------
    time_series : float array
       An array of time series, time as the last dimension.

    lb, ub: float, optional
       The upper and lower bound on the frequency band to be used in averaging
       defaults to 1,max(f)

    csd_method: dict, optional.
       See :func:`get_spectra` documentation for details

    Returns
    -------
    c : float
       This is an upper-diagonal array, where c[i][j] is the band-averaged
       coherency between time_series[i] and time_series[j]
    """

    if csd_method is None:
        csd_method = {'this_method': 'welch'}  # The default

    f, fxy = get_spectra(time_series, csd_method)

    lb_idx, ub_idx = utils.get_bounds(f, lb, ub)

    if lb == 0:
        lb_idx = 1  # The lowest frequency band should be f0

    c = np.zeros((time_series.shape[0],
                  time_series.shape[0]))

    for i in range(time_series.shape[0]):
        for j in range(i, time_series.shape[0]):
            c[i][j] = _coherence_bavg(fxy[i][j][lb_idx:ub_idx],
                                      fxy[i][i][lb_idx:ub_idx],
                                      fxy[j][j][lb_idx:ub_idx])

    idx = tril_indices(time_series.shape[0], -1)
    c[idx[0], idx[1], ...] = c[idx[1], idx[0], ...].conj()  # Make it symmetric

    return c

    def coherency(self):
        """The standard output for this kind of analyzer is the coherency """
        data = self.input.data
        tseries_length = data.shape[0]
        spectrum_length = self.spectrum.shape[-1]

        coherency = np.zeros((tseries_length, tseries_length, spectrum_length), dtype=complex)

        for i in xrange(tseries_length):
            for j in xrange(i, tseries_length):
                coherency[i][j] = tsa.coherency_spec(self.spectrum[i][j], self.spectrum[i][i], self.spectrum[j][j])

        idx = tril_indices(tseries_length, -1)
        coherency[idx[0], idx[1], ...] = coherency[idx[1], idx[0], ...].conj()

        return coherency

def coherence_partial(time_series, r, csd_method=None):
    r"""
    Compute the band-specific partial coherence between the spectra of
    two time series.

    The partial coherence is the part of the coherence between x and
    y, which cannot be attributed to a common cause, r.

    Input to this function is in the time domain.

    Parameters
    ----------

    time_series: float array
       An array of time-series, with time as the last dimension.

    r: float array
        This array represents the temporal sequence of the common cause to be
        partialed out, sampled at the same rate as time_series

    csd_method: dict, optional
       See :func:`get_spectra` documentation for details


    Returns
    -------
    f: array,
        The mid-frequencies of the frequency bands in the spectral
        decomposition

    c: float array
       The frequency dependent partial coherence between time_series i and
       time_series j in c[i][j] and in c[j][i], with r partialed out


    Notes
    -----

    This is an implementation of equation (2) of Sun (2004):

    .. math::

        Coh_{xy|r} = \frac{|{R_{xy}(\lambda) - R_{xr}(\lambda)
        R_{ry}(\lambda)}|^2}{(1-|{R_{xr}}|^2)(1-|{R_{ry}}|^2)}

    F.T. Sun and L.M. Miller and M. D'Esposito (2004). Measuring interregional
    functional connectivity using coherence and partial coherence analyses of
    fMRI data Neuroimage, 21: 647-58.
    """

    if csd_method is None:
        csd_method = {'this_method': 'welch'}  # The default

    f, fxy = get_spectra(time_series, csd_method)

    # Initialize c according to the size of f:
    c = np.zeros((time_series.shape[0],
                  time_series.shape[0],
                  f.shape[0]), dtype=complex)

    for i in range(time_series.shape[0]):
        for j in range(i, time_series.shape[0]):
            f, fxx, frr, frx = get_spectra_bi(time_series[i], r, csd_method)
            f, fyy, frr, fry = get_spectra_bi(time_series[j], r, csd_method)
            c[i, j] = coherence_partial_spec(fxy[i][j], fxy[i][i],
                                                  fxy[j][j], frx, fry, frr)

    idx = tril_indices(time_series.shape[0], -1)
    c[idx[0], idx[1], ...] = c[idx[1], idx[0], ...].conj()  # Make it symmetric

    return f, c

def coherency_bavg(time_series, lb=0, ub=None, csd_method=None):
    r"""
    Compute the band-averaged coherency between the spectra of two time series.

    Input to this function is in the time domain.

    Parameters
    ----------
    time_series: n*t float array
       an array of n different time series of length t each

    lb, ub: float, optional
       the upper and lower bound on the frequency band to be used in averaging
       defaults to 1,max(f)

    csd_method: dict, optional.
       See :func:`get_spectra` documentation for details

    Returns
    -------
    c: float array
        This is an upper-diagonal array, where c[i][j] is the band-averaged
        coherency between time_series[i] and time_series[j]

    Notes
    -----

    This is an implementation of equation (A4) of Sun(2005):

    .. math::

        \bar{Coh_{xy}} (\bar{\lambda}) =
        \frac{\left|{\sum_\lambda{\hat{f_{xy}}}}\right|^2}
        {\sum_\lambda{\hat{f_{xx}}}\cdot sum_\lambda{\hat{f_{yy}}}}

    F.T. Sun and L.M. Miller and M. D'Esposito (2005). Measuring
    temporal dynamics of functional networks using phase spectrum of fMRI
    data. Neuroimage, 28: 227-37.
    """
    if csd_method is None:
        csd_method = {'this_method': 'welch'}  # The default

    f, fxy = get_spectra(time_series, csd_method)

    lb_idx, ub_idx = utils.get_bounds(f, lb, ub)

    if lb == 0:
        lb_idx = 1  # The lowest frequency band should be f0

    c = np.zeros((time_series.shape[0],
                  time_series.shape[0]), dtype=complex)

    for i in range(time_series.shape[0]):
        for j in range(i, time_series.shape[0]):
            c[i][j] = _coherency_bavg(fxy[i][j][lb_idx:ub_idx],
                                      fxy[i][i][lb_idx:ub_idx],
                                      fxy[j][j][lb_idx:ub_idx])

    idx = tril_indices(time_series.shape[0], -1)
    c[idx[0], idx[1], ...] = c[idx[1], idx[0], ...].conj()  # Make it symmetric

    return c

def coherence_regularized(time_series, epsilon, alpha, csd_method=None):
    r"""
    Same as coherence, except regularized in order to overcome numerical
    imprecisions

    Parameters
    ----------

    time_series: n-d float array
       The time series data for which the regularized coherence is calculated

    epsilon: float
       Small regularization parameter. Should be much smaller than any
       meaningful value of coherence you might encounter

    alpha: float
       large regularization parameter. Should be much larger than any
       meaningful value of coherence you might encounter (preferably much
       larger than 1).

    csd_method: dict, optional.
       See :func:`get_spectra` documentation for details

    Returns
    -------
    f: float array
       The central frequencies for the frequency bands for which the spectra
       are estimated

    c: n-d array
       This is a symmetric matrix with the coherencys of the signals. The
       coherency of signal i and signal j is in f[i][j].

    Returns
    -------
    frequencies, coherence

    Notes
    -----
    The regularization scheme is as follows:

    .. math::

        C_{x,y} = \frac{(\alpha f_{xx} + \epsilon)^2}
        {\alpha^{2}((f_{xx}+\epsilon)(f_{yy}+\epsilon))}

    """
    if csd_method is None:
        csd_method = {'this_method': 'welch'}  # The default

    f, fxy = get_spectra(time_series, csd_method)

    #A container for the coherences, with the size and shape of the expected
    #output:
    c = np.zeros((time_series.shape[0],
                  time_series.shape[0],
                  f.shape[0]), complex)

    for i in range(time_series.shape[0]):
        for j in range(i, time_series.shape[0]):
            c[i][j] = _coherence_reqularized(fxy[i][j], fxy[i][i],
                                             fxy[j][j], epsilon, alpha)

    idx = tril_indices(time_series.shape[0], -1)
    c[idx[0], idx[1], ...] = c[idx[1], idx[0], ...].conj()  # Make it symmetric

    return f, c

#.........这里部分代码省略.........
       to use. This parameters represents trade-off between frequency
       resolution (lower main lobe BW for the taper) and variance reduction
       (higher BW and number of averaged estimates).

    adaptive : {True, False}
       Use adaptive weighting to combine spectra
    low_bias : {True, False}
       Rather than use 2NW tapers, only use the tapers that have better than
       90% spectral concentration within the bandwidth (still using
       a maximum of 2NW tapers)
    sides : str (optional)   [ 'default' | 'onesided' | 'twosided' ]
         This determines which sides of the spectrum to return.  For
         complex-valued inputs, the default is two-sided, for real-valued
         inputs, default is one-sided Indicates whether to return a one-sided
         or two-sided

    Returns
    -------
    (freqs, csd_est) : ndarrays
        The estimatated CSD and the frequency points vector.
        The CSD{i,j}(f) are returned in a square "matrix" of vectors
        holding Sij(f). For an input array of (M,N), the output is (M,M,N)
    """
    # have last axis be time series for now
    N = s.shape[-1]
    rest_of = s.shape[:-1]
    M = int(np.product(rest_of))

    s = s.reshape(M, N)
    # de-mean this sucker
    s = utils.remove_bias(s, axis=-1)

    #Get the number of tapers from the sampling rate and the bandwidth:
    if BW is not None:
        NW = BW / (2 * Fs) * N
    else:
        NW = 4

    Kmax = int(2 * NW)

    dpss, eigvals = dpss_windows(N, NW, Kmax)
    if low_bias:
        keepers = (eigvals > 0.9)
        dpss = dpss[keepers]
        eigvals = eigvals[keepers]
        Kmax = len(dpss)

    # if the time series is a complex vector, a one sided PSD is invalid:
    if (sides == 'default' and np.iscomplexobj(s)) or sides == 'twosided':
        sides = 'twosided'
    elif sides in ('default', 'onesided'):
        sides = 'onesided'

    sig_sl = [slice(None)] * len(s.shape)
    sig_sl.insert(len(s.shape) - 1, np.newaxis)

    # tapered.shape is (M, Kmax, N)
    tapered = s[sig_sl] * dpss

    # compute the y_{i,k}(f)
    tapered_spectra = fftpack.fft(tapered)

    # compute the cross-spectral density functions
    last_freq = N / 2 + 1 if sides == 'onesided' else N

    if adaptive:
        w = np.empty(tapered_spectra.shape[:-1] + (last_freq,))
        nu = np.empty((M, last_freq))
        for i in xrange(M):
            w[i], nu[i] = utils.adaptive_weights(
                tapered_spectra[i], eigvals, sides=sides
                )
    else:
        weights = np.sqrt(eigvals).reshape(Kmax, 1)

    csdfs = np.empty((M, M, last_freq), 'D')
    for i in xrange(M):
        if adaptive:
            wi = w[i]
        else:
            wi = weights
        for j in xrange(i + 1):
            if adaptive:
                wj = w[j]
            else:
                wj = weights
            ti = tapered_spectra[i]
            tj = tapered_spectra[j]
            csdfs[i, j] = mtm_cross_spectrum(ti, tj, (wi, wj), sides=sides)

    upper_idc = triu_indices(M, k=1)
    lower_idc = tril_indices(M, k=-1)
    csdfs[upper_idc] = csdfs[lower_idc].conj()

    if sides == 'onesided':
        freqs = np.linspace(0, Fs / 2, N / 2 + 1)
    else:
        freqs = np.linspace(0, Fs, N, endpoint=False)

    return freqs, csdfs

#.........这里部分代码省略.........
    The number of points N, or a precomputed FFT Sk may be provided. By
    default, the CSD function returned is normalized so that the integral of
    the PSD is equal to the mean squared amplitude (mean energy) of s (see
    Notes).

    Parameters
    ---------

    s : ndarray
        Signals for which to estimate the CSD, time dimension in the last axis

    Fs: float (optional)
       The sampling rate. Defaults to 2*pi

    Sk : ndarray (optional)
        Precomputed FFT of rows of s

    NFFT : int (optional)
        Indicates an N-point FFT where N != s.shape[-1]

    sides : str (optional)   [ 'default' | 'onesided' | 'twosided' ]
        This determines which sides of the spectrum to return.
        For complex-valued inputs, the default is two-sided, for real-valued
        inputs, default is one-sided Indicates whether to return a one-sided
        or two-sided

    normalize : boolean (optional)
        Normalizes the PSD

    Returns
    -------

    freqs, csd_est : ndarrays
        The estimatated CSD and the frequency points vector.
        The CSD{i,j}(f) are returned in a square "matrix" of vectors
        holding Sij(f). For an input array that is reshaped to (M,N),
        the output is (M,M,N)

    Notes
    -----
    setting dw = 2*PI/N, then the integral from -PI, PI (or 0,PI) of PSD/(2PI)
    will be nearly equal to sxy(0), where sxx is the crosscovariance function
    of s1(n), s2(n). By definition, sxy(0) = E{s1(n)s2*(n)} ~
    (s1*s2.conj()).mean()
    """
    s_shape = s.shape
    s.shape = (np.prod(s_shape[:-1]), s_shape[-1])
    # defining an Sk_loc is a little opaque, but it avoids having to
    # reset the shape of any user-given Sk later on
    if Sk is not None:
        Sk_shape = Sk.shape
        N = Sk.shape[-1]
        Sk_loc = Sk.reshape(np.prod(Sk_shape[:-1]), N)
    else:
        if NFFT is not None:
            N = NFFT
        else:
            N = s.shape[-1]
        Sk_loc = fftpack.fft(s, n=N)
    # reset s.shape
    s.shape = s_shape

    M = Sk_loc.shape[0]
    norm = float(s.shape[-1])

    # if the time series is a complex vector, a one sided PSD is invalid:
    if (sides == 'default' and np.iscomplexobj(s)) or sides == 'twosided':
        sides = 'twosided'
    elif sides in ('default', 'onesided'):
        sides = 'onesided'

    if sides == 'onesided':
        # putative Nyquist freq
        Fn = N / 2 + 1
        # last duplicate freq
        Fl = (N + 1) / 2
        csd_mat = np.empty((M, M, Fn), 'D')
        freqs = np.linspace(0, Fs / 2, Fn)
        for i in xrange(M):
            for j in xrange(i + 1):
                csd_mat[i, j, 0] = Sk_loc[i, 0] * Sk_loc[j, 0].conj()
                csd_mat[i, j, 1:Fl] = 2 * (Sk_loc[i, 1:Fl] *
                                           Sk_loc[j, 1:Fl].conj())
                if Fn > Fl:
                    csd_mat[i, j, Fn - 1] = (Sk_loc[i, Fn - 1] *
                                             Sk_loc[j, Fn - 1].conj())

    else:
        csd_mat = np.empty((M, M, N), 'D')
        freqs = np.linspace(0, Fs / 2, N, endpoint=False)
        for i in xrange(M):
            for j in xrange(i + 1):
                csd_mat[i, j] = Sk_loc[i] * Sk_loc[j].conj()
    if normalize:
        csd_mat /= norm

    upper_idc = triu_indices(M, k=1)
    lower_idc = tril_indices(M, k=-1)
    csd_mat[upper_idc] = csd_mat[lower_idc].conj()
    return freqs, csd_mat