def frequencies(self):
        """Get the central frequencies for the frequency bands, given the
           method of estimating the spectrum """

        self.method['Fs'] = self.method.get('Fs',self.input.sampling_rate)
        NFFT = self.method.get('NFFT',64)
        Fs = self.method.get('Fs')
        freqs = tsu.get_freqs(Fs,NFFT)
        lb_idx,ub_idx = tsu.get_bounds(freqs,self.lb,self.ub)
        
        return freqs[lb_idx:ub_idx]

    def frequencies(self):
        """Get the central frequencies for the frequency bands, given the
           method of estimating the spectrum """

        # Get the sampling rate from the seed time-series:
        self.method["Fs"] = self.method.get("Fs", self.seed.sampling_rate)
        NFFT = self.method.get("NFFT", 64)
        Fs = self.method.get("Fs")
        freqs = tsu.get_freqs(Fs, NFFT)
        lb_idx, ub_idx = tsu.get_bounds(freqs, self.lb, self.ub)

        return freqs[lb_idx:ub_idx]

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_phase_delay(time_series, lb=0, ub=None, csd_method=None):
    """
    The temporal delay calculated from the coherency phase spectrum.

    Parameters
    ----------

    time_series: float array
       The time-series data for which the delay is calculated.

    lb, ub: float
       Frequency boundaries (in Hz), for the domain over which the delays are
       calculated. Defaults to 0-max(f)

    csd_method : dict, optional.
       See :func:`get_spectra`

    Returns
    -------
    f : float array
       The mid-frequencies for the frequency bands over which the calculation
       is done.
    p : float array
       Pairwise temporal delays between time-series (in seconds).

    """
    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_idx == 0:
        lb_idx = 1

    p = np.zeros((time_series.shape[0], time_series.shape[0],
                  f[lb_idx:ub_idx].shape[-1]))

    for i in range(time_series.shape[0]):
        for j in range(i, time_series.shape[0]):
            p[i][j] = _coherency_phase_delay(f[lb_idx:ub_idx],
                                             fxy[i][j][lb_idx:ub_idx])
            p[j][i] = _coherency_phase_delay(
                                f[lb_idx:ub_idx],
                                fxy[i][j][lb_idx:ub_idx].conjugate())

    return f[lb_idx:ub_idx], p

def cache_fft(time_series, ij, lb=0, ub=None,
                  method=None, prefer_speed_over_memory=False,
                  scale_by_freq=True):
    """compute and cache the windowed FFTs of the time_series, in such a way
    that computing the psd and csd of any combination of them can be done
    quickly.

    Parameters
    ----------

    time_series : float array
       An ndarray with time-series, where time is the last dimension

    ij: list of tuples
      Each tuple in this variable should contain a pair of
      indices of the form (i,j). The resulting cache will contain the fft of
      time-series in the rows indexed by the unique elements of the union of i
      and j

    lb,ub: float
       Define a frequency band of interest, for which the fft will be cached

    method: dict, optional
        See :func:`get_spectra` for details on how this is used. For this set
        of functions, 'this_method' has to be 'welch'


    Returns
    -------
    freqs, cache

        where: cache =
             {'FFT_slices':FFT_slices,'FFT_conj_slices':FFT_conj_slices,
             'norm_val':norm_val}

    Notes
    -----

    - For these functions, only the Welch windowed periodogram ('welch') is
      available.

    - Detrending the input is not an option here, in order to save
      time on an empty function call.

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

    this_method = method.get('this_method', 'welch')

    if this_method == 'welch':
        NFFT = method.get('NFFT', 64)
        Fs = method.get('Fs', 2 * np.pi)
        window = method.get('window', mlab.window_hanning)
        n_overlap = method.get('n_overlap', int(np.ceil(NFFT / 2.0)))
    else:
        e_s = "For cache_fft, spectral estimation method must be welch"
        raise ValueError(e_s)
    time_series = utils.zero_pad(time_series, NFFT)

    #The shape of the zero-padded version:
    n_channels, n_time_points = time_series.shape

    # get all the unique channels in time_series that we are interested in by
    # checking the ij tuples
    all_channels = set()
    for i, j in ij:
        all_channels.add(i)
        all_channels.add(j)

    # for real time_series, ignore the negative frequencies
    if np.iscomplexobj(time_series):
        n_freqs = NFFT
    else:
        n_freqs = NFFT // 2 + 1

    #Which frequencies
    freqs = utils.get_freqs(Fs, NFFT)

    #If there are bounds, limit the calculation to within that band,
    #potentially include the DC component:
    lb_idx, ub_idx = utils.get_bounds(freqs, lb, ub)

    n_freqs = ub_idx - lb_idx
    #Make the window:
    if mlab.cbook.iterable(window):
        assert(len(window) == NFFT)
        window_vals = window
    else:
        window_vals = window(np.ones(NFFT, time_series.dtype))

    #Each fft needs to be normalized by the square of the norm of the window
    #and, for consistency with newer versions of mlab.csd (which, in turn, are
    #consistent with Matlab), normalize also by the sampling rate:

    if scale_by_freq:
        #This is the normalization factor for one-sided estimation, taking into
        #account the sampling rate. This makes the PSD a density function, with
        #units of dB/Hz, so that integrating over frequencies gives you the RMS
        #(XXX this should be in the tests!).
#.........这里部分代码省略.........

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 plot_snr_diff(tseries1, tseries2, lb=0, ub=None, fig=None,
                  ts_names=['1', '2'],
                  bandwidth=None, adaptive=False, low_bias=True):

    """
    Show distributions of differences between two time-series in the
    amount of snr (freq band by freq band) and information. For example,
    for comparing two stimulus conditions

    Parameters
    ----------
    tseries1, tseries2 : nitime TimeSeries objects
       These are the time-series to compare, with each of them having the
       dims: (n_channels, n_reps, time), where n_channels1 = n_channels2

    lb,ub: float
       Lower and upper bounds on the frequency range over which to
       calculate the information rate (default to [0,Nyquist]).

    fig: matplotlib figure object
       If you want to do this on already existing figure. Otherwise, a new
       figure object will be generated.

    ts_names: list of str
       Labels for the two inputs, to be used in plotting (defaults to
       ['1','2'])

    bandwidth, adaptive, low_bias: See :func:`nta.SNRAnalyzer` for details


    Returns
    -------

    A tuple containing:

    fig: a matplotlib figure object
        This figure displays:
        1. The histogram of the information differences between the two
        time-series
        2. The frequency-dependent SNR for the two time-series

    info1, info2: float arrays
        The frequency-dependent information rates (in bits/sec)

    s_n_r1, s_n_r2: float arrays
         The frequncy-dependent signal-to-noise ratios

    """
    if fig is None:
        fig = plt.figure()
    ax_scatter = fig.add_subplot(1, 2, 1)
    ax_snr = fig.add_subplot(1, 2, 2)

    SNR1 = []
    s_n_r1 = []
    info1 = []
    SNR2 = []
    info2 = []
    s_n_r2 = []

    #If you only have one channel, make sure that everything still works by
    #adding an axis
    if len(tseries1.data.shape) < 3:
        this1 = tseries1.data[np.newaxis, :, :]
        this2 = tseries2.data[np.newaxis, :, :]
    else:
        this1 = tseries1.data
        this2 = tseries2.data

    for i in range(this1.shape[0]):
        SNR1.append(nta.SNRAnalyzer(ts.TimeSeries(this1[i],
                                    sampling_rate=tseries1.sampling_rate),
                                bandwidth=bandwidth,
                                adaptive=adaptive,
                                low_bias=low_bias))
        info1.append(SNR1[-1].mt_information)
        s_n_r1.append(SNR1[-1].mt_snr)

        SNR2.append(nta.SNRAnalyzer(ts.TimeSeries(this2[i],
                                    sampling_rate=tseries2.sampling_rate),
                                bandwidth=bandwidth,
                                adaptive=adaptive,
                                low_bias=low_bias))

        info2.append(SNR2[-1].mt_information)
        s_n_r2.append(SNR2[-1].mt_snr)

    freqs = SNR1[-1].mt_frequencies

    lb_idx, ub_idx = tsu.get_bounds(freqs, lb, ub)
    freqs = freqs[lb_idx:ub_idx]

    info1 = np.array(info1)
    info_sum1 = np.sum(info1[:, lb_idx:ub_idx], -1)
    info2 = np.array(info2)
    info_sum2 = np.sum(info2[:, lb_idx:ub_idx], -1)

    ax_scatter.scatter(info_sum1, info_sum2)
    ax_scatter.errorbar(np.mean(info_sum1), np.mean(info_sum2),
                 yerr=np.std(info_sum2),
#.........这里部分代码省略.........

def plot_snr(tseries, lb=0, ub=None, fig=None):
    """
    Show the coherence, snr and information of an SNRAnalyzer

    Parameters
    ----------
    tseries: nitime TimeSeries object
       Multi-trial data in response to one stimulus/protocol with the dims:
       (n_channels,n_repetitions,time)

    lb,ub: float
       Lower and upper bounds on the frequency range over which to
       calculate (default to [0,Nyquist]).

    Returns
    -------

    A tuple containing:

    fig: a matplotlib figure object
        This figure displays:
        1. Coherence
        2. SNR
        3. Information
    """

    if fig is None:
        fig = plt.figure()

    ax_spectra = fig.add_subplot(1, 2, 1)
    ax_snr_info = fig.add_subplot(1, 2, 2)

    A = []
    info = []
    s_n_r = []
    coh = []
    noise_spectra = []
    signal_spectra = []
    #If you only have one channel, make sure that everything still works by
    #adding an axis
    if len(tseries.data.shape) < 3:
        this = tseries.data[np.newaxis, :, :]
    else:
        this = tseries.data

    for i in range(this.shape[0]):
        A.append(nta.SNRAnalyzer(ts.TimeSeries(this[i],
                                    sampling_rate=tseries.sampling_rate)))
        info.append(A[-1].mt_information)
        s_n_r.append(A[-1].mt_snr)
        coh.append(A[-1].mt_coherence)
        noise_spectra.append(A[-1].mt_noise_psd)
        signal_spectra.append(A[-1].mt_signal_psd)

    freqs = A[-1].mt_frequencies

    lb_idx, ub_idx = tsu.get_bounds(freqs, lb, ub)
    freqs = freqs[lb_idx:ub_idx]

    coh_mean = np.mean(coh, 0)
    snr_mean = np.mean(s_n_r, 0)
    info_mean = np.mean(info, 0)
    n_spec_mean = np.mean(noise_spectra, 0)
    s_spec_mean = np.mean(signal_spectra, 0)

    ax_spectra.plot(freqs, np.log(s_spec_mean[lb_idx:ub_idx]), label='Signal')
    ax_spectra.plot(freqs, np.log(n_spec_mean[lb_idx:ub_idx]), label='Noise')
    ax_spectra.set_xlabel('Frequency (Hz)')
    ax_spectra.set_ylabel('Spectral power (dB)')

    ax_snr_info.plot(freqs, snr_mean[lb_idx:ub_idx], label='SNR')
    ax_snr_info.plot(np.nan, np.nan, 'r', label='Info')
    ax_snr_info.set_ylabel('SNR')
    ax_snr_info.set_xlabel('Frequency (Hz)')
    ax_info = ax_snr_info.twinx()
    ax_info.plot(freqs, np.cumsum(info_mean[lb_idx:ub_idx]), 'r')
    ax_info.set_ylabel('Cumulative information rate (bits/sec)')
    return fig

inp_sampling_rate = 1000.0
lb = 0  # Hz
ub = 1000  # Hz

_, signal_spectra, _ = tsa.multi_taper_psd(signal, Fs=inp_sampling_rate, BW=None, adaptive=True, low_bias=True)
_, noise_spectra, _ = tsa.multi_taper_psd(
    restored - signal, Fs=inp_sampling_rate, BW=None, adaptive=True, low_bias=True
)


freqs = np.linspace(0, inp_sampling_rate / 2, signal.shape[-1] / 2 + 1)

f = plt.figure()
ax = f.add_subplot(1, 2, 1)
ax_snr_info = f.add_subplot(1, 2, 2)
lb_idx, ub_idx = tsu.get_bounds(freqs, lb, ub)
freqs = freqs[lb_idx:ub_idx]
snr = signal_spectra / noise_spectra

ax.plot(freqs, np.log(signal_spectra[lb_idx:ub_idx]), label="Signal")
ax.plot(freqs, np.log(noise_spectra[lb_idx:ub_idx]), label="Noise")

ax_snr_info.plot(freqs, snr[lb_idx:ub_idx], label="SNR")
ax_snr_info.plot(np.nan, np.nan, "r", label="Info")
ax_snr_info.set_ylabel("SNR")
ax_snr_info.set_xlabel("Frequency (Hz)")

ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Spectral power (dB)")

plt.show()