def test_csd_morlet():
    """Test computing cross-spectral density using Morlet wavelets."""
    epochs = _generate_coherence_data()
    sfreq = epochs.info['sfreq']

    # Compute CSDs by a variety of methods
    freqs = [10, 15, 22]
    n_cycles = [20, 30, 44]
    times = [(None, None), (1, 9)]
    as_arrays = [False, True]
    parameters = product(times, as_arrays)
    for (tmin, tmax), as_array in parameters:
        if as_array:
            csd = csd_array_morlet(epochs.get_data(), sfreq, freqs,
                                   t0=epochs.tmin, n_cycles=n_cycles,
                                   tmin=tmin, tmax=tmax,
                                   ch_names=epochs.ch_names)
        else:
            csd = csd_morlet(epochs, frequencies=freqs, n_cycles=n_cycles,
                             tmin=tmin, tmax=tmax)
        if tmin is None and tmax is None:
            assert csd.tmin == 0 and csd.tmax == 9.98
        else:
            assert csd.tmin == tmin and csd.tmax == tmax
        _test_csd_matrix(csd)

    # CSD diagonals should contain PSD
    tfr = tfr_morlet(epochs, freqs, n_cycles, return_itc=False)
    power = np.mean(tfr.data, 2)
    csd = csd_morlet(epochs, frequencies=freqs, n_cycles=n_cycles)
    assert_allclose(csd._data[[0, 3, 5]] * sfreq, power)

    # Test using plain convolution instead of FFT
    csd = csd_morlet(epochs, frequencies=freqs, n_cycles=n_cycles,
                     use_fft=False)
    assert_allclose(csd._data[[0, 3, 5]] * sfreq, power)

    # Test baselining warning
    epochs_nobase = epochs.copy()
    epochs_nobase.baseline = None
    epochs_nobase.info['highpass'] = 0
    with pytest.warns(RuntimeWarning, match='baseline'):
        csd = csd_morlet(epochs_nobase, frequencies=[10], decim=20)

def test_get_inst_data():
    """Test _get_inst_data."""
    raw = read_raw_fif(fname_raw)
    raw.crop(tmax=1.)
    assert_equal(_get_inst_data(raw), raw._data)
    raw.pick_channels(raw.ch_names[:2])

    epochs = _segment_raw(raw, 0.5)
    assert_equal(_get_inst_data(epochs), epochs._data)

    evoked = epochs.average()
    assert_equal(_get_inst_data(evoked), evoked.data)

    evoked.crop(tmax=0.1)
    picks = list(range(2))
    freqs = np.array([50., 55.])
    n_cycles = 3
    tfr = tfr_morlet(evoked, freqs, n_cycles, return_itc=False, picks=picks)
    assert_equal(_get_inst_data(tfr), tfr.data)

    assert_raises(TypeError, _get_inst_data, 'foo')

def morlet_analysis(epochs, n_cycles=4):
    """

    Parameters
    ----------
    epochs : list of epochs

    Returns
    -------
    result : numpy array
        The result of the multitaper analysis.

    """
    frequencies = np.arange(6., 30., 2.)
    # n_cycles = frequencies / 2.

    power, plv = tfr_morlet(epochs, freqs=frequencies, n_cycles=n_cycles,
                            return_itc=True,
                            verbose=True)

    return power, plv

###############################################################################
# Time-frequency analysis: power and inter-trial coherence
# --------------------------------------------------------
#
# We now compute time-frequency representations (TFRs) from our Epochs.
# We'll look at power and inter-trial coherence (ITC).
#
# To this we'll use the function :func:`mne.time_frequency.tfr_morlet`
# but you can also use :func:`mne.time_frequency.tfr_multitaper`
# or :func:`mne.time_frequency.tfr_stockwell`.

# define frequencies of interest (log-spaced)
freqs = np.logspace(*np.log10([6, 35]), num=8)
n_cycles = freqs / 2.  # different number of cycle per frequency
power, itc = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles, use_fft=True,
                        return_itc=True, decim=3, n_jobs=1)

###############################################################################
# Inspect power
# -------------
#
# .. note::
#     The generated figures are interactive. In the topo you can click
#     on an image to visualize the data for one sensor.
#     You can also select a portion in the time-frequency plane to
#     obtain a topomap for a certain time-frequency region.
power.plot_topo(baseline=(-0.5, 0), mode='logratio', title='Average power')
power.plot([82], baseline=(-0.5, 0), mode='logratio', title=power.ch_names[82])

fig, axis = plt.subplots(1, 2, figsize=(7, 4))
power.plot_topomap(ch_type='grad', tmin=0.5, tmax=1.5, fmin=8, fmax=12,

def test_compute_epochs_csd():
    """Test computing cross-spectral density from epochs
    """
    epochs, epochs_sin = _get_data()
    # Check that wrong parameters are recognized
    assert_raises(ValueError, compute_epochs_csd, epochs, mode='notamode')
    assert_raises(ValueError, compute_epochs_csd, epochs, fmin=20, fmax=10)
    assert_raises(ValueError, compute_epochs_csd, epochs, fmin=20, fmax=20.1)
    assert_raises(ValueError, compute_epochs_csd, epochs, tmin=0.15, tmax=0.1)
    assert_raises(ValueError, compute_epochs_csd, epochs, tmin=0, tmax=10)
    assert_raises(ValueError, compute_epochs_csd, epochs, tmin=10, tmax=11)

    data_csd_mt = compute_epochs_csd(epochs, mode='multitaper', fmin=8,
                                     fmax=12, tmin=0.04, tmax=0.15)
    data_csd_fourier = compute_epochs_csd(epochs, mode='fourier', fmin=8,
                                          fmax=12, tmin=0.04, tmax=0.15)

    # Check shape of the CSD matrix
    n_chan = len(data_csd_mt.ch_names)
    assert_equal(data_csd_mt.data.shape, (n_chan, n_chan))
    assert_equal(data_csd_fourier.data.shape, (n_chan, n_chan))

    # Check if the CSD matrix is hermitian
    assert_array_equal(np.tril(data_csd_mt.data).T.conj(),
                       np.triu(data_csd_mt.data))
    assert_array_equal(np.tril(data_csd_fourier.data).T.conj(),
                       np.triu(data_csd_fourier.data))

    # Computing induced power for comparison
    epochs.crop(tmin=0.04, tmax=0.15)
    tfr = tfr_morlet(epochs, freqs=[10], n_cycles=0.6, return_itc=False)
    power = np.mean(tfr.data, 2)

    # Maximum PSD should occur for specific channel
    max_ch_power = power.argmax()
    max_ch_mt = data_csd_mt.data.diagonal().argmax()
    max_ch_fourier = data_csd_fourier.data.diagonal().argmax()
    assert_equal(max_ch_mt, max_ch_power)
    assert_equal(max_ch_fourier, max_ch_power)

    # Maximum CSD should occur for specific channel
    ch_csd_mt = [np.abs(data_csd_mt.data[max_ch_power][i])
                 if i != max_ch_power else 0 for i in range(n_chan)]
    max_ch_csd_mt = np.argmax(ch_csd_mt)
    ch_csd_fourier = [np.abs(data_csd_fourier.data[max_ch_power][i])
                      if i != max_ch_power else 0 for i in range(n_chan)]
    max_ch_csd_fourier = np.argmax(ch_csd_fourier)
    assert_equal(max_ch_csd_mt, max_ch_csd_fourier)

    # Check a list of CSD matrices is returned for multiple frequencies within
    # a given range when fsum=False
    csd_fsum = compute_epochs_csd(epochs, mode='fourier', fmin=8, fmax=20,
                                  fsum=True)
    csds = compute_epochs_csd(epochs, mode='fourier', fmin=8, fmax=20,
                              fsum=False)
    freqs = [csd.frequencies[0] for csd in csds]

    csd_sum = np.zeros_like(csd_fsum.data)
    for csd in csds:
        csd_sum += csd.data

    assert(len(csds) == 2)
    assert(len(csd_fsum.frequencies) == 2)
    assert_array_equal(csd_fsum.frequencies, freqs)
    assert_array_equal(csd_fsum.data, csd_sum)

# Factor to down-sample the temporal dimension of the TFR computed by
# tfr_morlet.
decim = 2
freqs = np.arange(7, 30, 3)  # define frequencies of interest
n_cycles = freqs / freqs[0]
zero_mean = False  # don't correct morlet wavelet to be of mean zero
# To have a true wavelet zero_mean should be True but here for illustration
# purposes it helps to spot the evoked response.

###############################################################################
# Create TFR representations for all conditions
# ---------------------------------------------
epochs_power = list()
for condition in [epochs[k] for k in event_id]:
    this_tfr = tfr_morlet(condition, freqs, n_cycles=n_cycles,
                          decim=decim, average=False, zero_mean=zero_mean,
                          return_itc=False)
    this_tfr.apply_baseline(mode='ratio', baseline=(None, 0))
    this_power = this_tfr.data[:, 0, :, :]  # we only have one channel.
    epochs_power.append(this_power)

###############################################################################
# Setup repeated measures ANOVA
# -----------------------------
#
# We will tell the ANOVA how to interpret the data matrix in terms of factors.
# This is done via the factor levels argument which is a list of the number
# factor levels for each factor.

n_conditions = len(epochs.event_id)
n_replications = epochs.events.shape[0] // n_conditions

# Take only one channel
ch_name = 'MEG 1332'
epochs.pick_channels([ch_name])

evoked = epochs.average()

# Factor to down-sample the temporal dimension of the TFR computed by
# tfr_morlet. Decimation occurs after frequency decomposition and can
# be used to reduce memory usage (and possibly computational time of downstream
# operations such as nonparametric statistics) if you don't need high
# spectrotemporal resolution.
decim = 5
freqs = np.arange(8, 40, 2)  # define frequencies of interest
sfreq = raw.info['sfreq']  # sampling in Hz
tfr_epochs = tfr_morlet(epochs, freqs, n_cycles=4., decim=decim,
                        average=False, return_itc=False, n_jobs=1)

# Baseline power
tfr_epochs.apply_baseline(mode='logratio', baseline=(-.100, 0))

# Crop in time to keep only what is between 0 and 400 ms
evoked.crop(0., 0.4)
tfr_epochs.crop(0., 0.4)

epochs_power = tfr_epochs.data[:, 0, :, :]  # take the 1 channel

###############################################################################
# Compute statistic
# -----------------
threshold = 2.5
T_obs, clusters, cluster_p_values, H0 = \

#mne.io.Raw.plot(raw=raw_memmaped, duration=2, start=20, n_channels=20, scalings={'eeg': 8000}, remove_dc=True)

id = 1
events_mne = np.c_[np.array(events), np.zeros(len(events), dtype=int), id * np.ones(len(events), dtype=int)]
baseline = (-2.5, -2.3)
event_id = dict(left_paw=id)
epochs = mne.Epochs(raw_memmaped, events_mne, event_id, -3, 3, proj=True, picks=None, baseline=baseline, preload=True, reject=None)
averaged = epochs.average()

power = pickle.load( open(os.path.join(path, "Analysis\\tfr_power.p"), "rb"))

n_cycles = 3
frequencies = np.arange(5, 60, 3)

from mne.time_frequency import tfr_morlet
power, phase_lock = tfr_morlet(epochs, freqs=frequencies, n_cycles=n_cycles, decim=3000, n_jobs=10)



import gui_tfr_viewer
gui_tfr_viewer.TFR_Viewer(power)

box = (0, 0.8, 0, 1.1)
w, h = [.09, .05]

pos = [[ut.normList([x, y], normalizeTo=0.8, vMin=1, vMax=8)[0], ut.normList([x, y], vMin=1, vMax=16)[1], w, h] for [n, s, (x,y)] in cp.sort_index(0, by='Numbers', ascending=True).values]
layout = mne.layouts.Layout(box, pos, cp.sort_index(0, by='Numbers', ascending=True).Strings, cp.sort_index(0, by='Numbers', ascending=True).Numbers, '128ch')

power.plot_topo(picks=None, tmin=-3, tmax=3, fmin=5, fmax=60, vmin=-3e10, vmax=3e10, layout=layout, layout_scale=None)

    ax.set_title('Sim: Using S transform, width = {:0.1f}'.format(width))
plt.tight_layout()

###############################################################################
# Morlet Wavelets
# ===============
#
# Finally, show the TFR using morlet wavelets, which are a sinusoidal wave
# with a gaussian envelope. We can control the balance between spectral and
# temporal resolution with the ``n_cycles`` parameter, which defines the
# number of cycles to include in the window.

fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharey=True)
all_n_cycles = [1, 3, freqs / 2.]
for n_cycles, ax in zip(all_n_cycles, axs):
    power = tfr_morlet(epochs, freqs=freqs,
                       n_cycles=n_cycles, return_itc=False)
    power.plot([0], baseline=(0., 0.1), mode='mean', vmin=vmin, vmax=vmax,
               axes=ax, show=False, colorbar=False)
    n_cycles = 'scaled by freqs' if not isinstance(n_cycles, int) else n_cycles
    ax.set_title('Sim: Using Morlet wavelet, n_cycles = %s' % n_cycles)
plt.tight_layout()

###############################################################################
# Calculating a TFR without averaging over epochs
# -----------------------------------------------
#
# It is also possible to calculate a TFR without averaging across trials.
# We can do this by using ``average=False``. In this case, an instance of
# :class:`mne.time_frequency.EpochsTFR` is returned.

n_cycles = freqs / 2.

import numpy as np
from my_settings import (epochs_folder, tf_folder)
from mne.time_frequency import tfr_morlet
import sys

subject = sys.argv[1]

freqs = np.arange(8, 13, 1)  # define frequencies of interest
n_cycles = 4.  # freqs / 2.  # different number of cycle per frequency

sides = ["left", "right"]
conditions = ["ctl", "ent"]

epochs = mne.read_epochs(
    epochs_folder + "%s_trial_start-epo.fif" % subject,
    preload=False)
for cond in conditions:
    for side in sides:
        power, itc = tfr_morlet(epochs[cond + "/" + side],
                                freqs=freqs,
                                n_cycles=n_cycles,
                                use_fft=True,
                                return_itc=False,
                                decim=2,
                                average=False,
                                n_jobs=1)
        power.save(tf_folder + "%s_%s_%s-4-tfr.h5" % (subject, cond, side),
                   overwrite=True)
        itc.save(tf_folder + "%s_%s_%s-4-tfr.h5" % (subject, cond, side),
                 overwrite=True)

                                picks=picks, baseline=(None, 0),
                                reject=reject, preload=True)
epochs_condition_2.pick_channels([ch_name])

###############################################################################
# Factor to downsample the temporal dimension of the TFR computed by
# tfr_morlet. Decimation occurs after frequency decomposition and can
# be used to reduce memory usage (and possibly comptuational time of downstream
# operations such as nonparametric statistics) if you don't need high
# spectrotemporal resolution.
decim = 2
freqs = np.arange(7, 30, 3)  # define frequencies of interest
n_cycles = 1.5

tfr_epochs_1 = tfr_morlet(epochs_condition_1, freqs,
                          n_cycles=n_cycles, decim=decim,
                          return_itc=False, average=False)

tfr_epochs_2 = tfr_morlet(epochs_condition_2, freqs,
                          n_cycles=n_cycles, decim=decim,
                          return_itc=False, average=False)

tfr_epochs_1.apply_baseline(mode='ratio', baseline=(None, 0))
tfr_epochs_2.apply_baseline(mode='ratio', baseline=(None, 0))

epochs_power_1 = tfr_epochs_1.data[:, 0, :, :]  # only 1 channel as 3D matrix
epochs_power_2 = tfr_epochs_2.data[:, 0, :, :]  # only 1 channel as 3D matrix

###############################################################################
# Compute statistic
# -----------------

def test_compute_csd():
    """Test computing cross-spectral density from ndarray. """

    epochs = _get_data(mode='real')

    tmin = 0.04
    tmax = 0.15
    tmp = np.where(np.logical_and(epochs.times >= tmin,
                                  epochs.times <= tmax))[0]

    picks_meeg = mne.pick_types(epochs[0].info, meg=True, eeg=True, eog=False,
                                ref_meg=False, exclude='bads')

    epochs_data = [e[picks_meeg][:, tmp].copy() for e in epochs]
    n_trials = len(epochs)
    n_series = len(picks_meeg)
    X = np.concatenate(epochs_data, axis=0)
    X = np.reshape(X, (n_trials, n_series, -1))
    X_list = epochs_data

    sfreq = epochs.info['sfreq']

    # Check data types and sizes are checked
    diff_types = [np.random.randn(3, 5), "error"]
    err_data = [np.random.randn(3, 5), np.random.randn(2, 4)]
    assert_raises(ValueError, csd_array, err_data, sfreq)
    assert_raises(ValueError, csd_array, diff_types, sfreq)
    assert_raises(ValueError, csd_array, np.random.randn(3), sfreq)

    # Check that wrong parameters are recognized
    assert_raises(ValueError, csd_array, X, sfreq, mode='notamode')
    assert_raises(ValueError, csd_array, X, sfreq, fmin=20, fmax=10)
    assert_raises(ValueError, csd_array, X, sfreq, fmin=20, fmax=20.1)

    data_csd_mt, freqs_mt = csd_array(X, sfreq, mode='multitaper',
                                      fmin=8, fmax=12)
    data_csd_fourier, freqs_fft = csd_array(X, sfreq, mode='fourier',
                                            fmin=8, fmax=12)

    # Test as list too
    data_csd_mt_list, freqs_mt_list = csd_array(X_list, sfreq,
                                                mode='multitaper',
                                                fmin=8, fmax=12)
    data_csd_fourier_list, freqs_fft_list = csd_array(X_list, sfreq,
                                                      mode='fourier',
                                                      fmin=8, fmax=12)

    assert_array_equal(data_csd_mt, data_csd_mt_list)
    assert_array_equal(data_csd_fourier, data_csd_fourier_list)
    assert_array_equal(freqs_mt, freqs_mt_list)
    assert_array_equal(freqs_fft, freqs_fft_list)

    # Check shape of the CSD matrix
    n_chan = len(epochs.ch_names)
    assert_equal(data_csd_mt.shape, (n_chan, n_chan))
    assert_equal(data_csd_fourier.shape, (n_chan, n_chan))

    # Check if the CSD matrix is hermitian
    assert_array_equal(np.tril(data_csd_mt).T.conj(),
                       np.triu(data_csd_mt))
    assert_array_equal(np.tril(data_csd_fourier).T.conj(),
                       np.triu(data_csd_fourier))

    # Computing induced power for comparison
    epochs.crop(tmin=0.04, tmax=0.15)
    tfr = tfr_morlet(epochs, freqs=[10], n_cycles=0.6, return_itc=False)
    power = np.mean(tfr.data, 2)

    # Maximum PSD should occur for specific channel
    max_ch_power = power.argmax()
    max_ch_mt = data_csd_mt.diagonal().argmax()
    max_ch_fourier = data_csd_fourier.diagonal().argmax()
    assert_equal(max_ch_mt, max_ch_power)
    assert_equal(max_ch_fourier, max_ch_power)

    # Maximum CSD should occur for specific channel
    ch_csd_mt = np.abs(data_csd_mt[max_ch_power])
    ch_csd_mt[max_ch_power] = 0.
    max_ch_csd_mt = np.argmax(ch_csd_mt)
    ch_csd_fourier = np.abs(data_csd_fourier[max_ch_power])
    ch_csd_fourier[max_ch_power] = 0.
    max_ch_csd_fourier = np.argmax(ch_csd_fourier)
    assert_equal(max_ch_csd_mt, max_ch_csd_fourier)

    # Check a list of CSD matrices is returned for multiple frequencies within
    # a given range when fsum=False
    csd_fsum, freqs_fsum = csd_array(X, sfreq, mode='fourier', fmin=8,
                                     fmax=20, fsum=True)
    csds, freqs = csd_array(X, sfreq, mode='fourier', fmin=8, fmax=20,
                            fsum=False)

    csd_sum = np.sum(csds, axis=2)

    assert_equal(csds.shape[2], 2)
    assert_equal(len(freqs), 2)
    assert_array_equal(freqs_fsum, freqs)
    assert_array_equal(csd_fsum, csd_sum)

def test_compute_csd():
    """Test computing cross-spectral density from ndarray."""
    epochs = _get_real_data()

    tmin = 0.04
    tmax = 0.15
    tmp = np.where(np.logical_and(epochs.times >= tmin,
                                  epochs.times <= tmax))[0]

    picks_meeg = mne.pick_types(epochs[0].info, meg=True, eeg=True, eog=False,
                                ref_meg=False, exclude='bads')

    epochs_data = [e[picks_meeg][:, tmp].copy() for e in epochs]
    n_trials = len(epochs)
    n_series = len(picks_meeg)
    X = np.concatenate(epochs_data, axis=0)
    X = np.reshape(X, (n_trials, n_series, -1))
    X_list = epochs_data

    sfreq = epochs.info['sfreq']

    # Check data types and sizes are checked
    diff_types = [np.random.randn(3, 5), "error"]
    err_data = [np.random.randn(3, 5), np.random.randn(2, 4)]
    with warnings.catch_warnings(record=True):  # deprecation
        raises(ValueError, csd_array, err_data, sfreq)
        raises(ValueError, csd_array, diff_types, sfreq)
        raises(ValueError, csd_array, np.random.randn(3), sfreq)

        # Check that wrong parameters are recognized
        raises(ValueError, csd_array, X, sfreq, mode='notamode')
        raises(ValueError, csd_array, X, sfreq, fmin=20, fmax=10)
        raises(ValueError, csd_array, X, sfreq, fmin=20, fmax=20.1)

    # Test deprecation warning
    with warnings.catch_warnings(record=True) as ws:
        warnings.simplefilter('always')
        csd_mt = csd_array(X, sfreq, mode='multitaper', fmin=8, fmax=12)
    assert len([w for w in ws
                if issubclass(w.category, DeprecationWarning)]) == 1

    with warnings.catch_warnings(record=True):  # deprecation
        csd_fourier = csd_array(X, sfreq, mode='fourier', fmin=8, fmax=12)

    # Test as list too
    with warnings.catch_warnings(record=True):  # deprecation
        csd_mt_list = csd_array(X_list, sfreq, mode='multitaper',
                                fmin=8, fmax=12)
        csd_fourier_list = csd_array(X_list, sfreq, mode='fourier', fmin=8,
                                     fmax=12)

    assert_array_equal(csd_mt._data, csd_mt_list._data)
    assert_array_equal(csd_fourier._data, csd_fourier_list._data)
    assert_array_equal(csd_mt.frequencies, csd_mt_list.frequencies)
    assert_array_equal(csd_fourier.frequencies, csd_fourier_list.frequencies)

    # Check shape of the CSD matrix
    n_chan = len(epochs.ch_names)
    csd_mt_data = csd_mt.get_data()
    csd_fourier_data = csd_fourier.get_data()
    assert csd_mt_data.shape == (n_chan, n_chan)
    assert csd_fourier_data.shape == (n_chan, n_chan)

    # Check if the CSD matrix is hermitian
    assert_array_equal(np.tril(csd_mt_data).T.conj(),
                       np.triu(csd_mt_data))
    assert_array_equal(np.tril(csd_fourier_data).T.conj(),
                       np.triu(csd_fourier_data))

    # Computing induced power for comparison
    epochs.crop(tmin=0.04, tmax=0.15)
    tfr = tfr_morlet(epochs, freqs=[10], n_cycles=0.6, return_itc=False)
    power = np.mean(tfr.data, 2)

    # Maximum PSD should occur for specific channel
    max_ch_power = power.argmax()
    max_ch_mt = csd_mt_data.diagonal().argmax()
    max_ch_fourier = csd_fourier_data.diagonal().argmax()
    assert max_ch_mt == max_ch_power
    assert max_ch_fourier == max_ch_power

    # Maximum CSD should occur for specific channel
    ch_csd_mt = np.abs(csd_mt_data[max_ch_power])
    ch_csd_mt[max_ch_power] = 0.
    max_ch_csd_mt = np.argmax(ch_csd_mt)
    ch_csd_fourier = np.abs(csd_fourier_data[max_ch_power])
    ch_csd_fourier[max_ch_power] = 0.
    max_ch_csd_fourier = np.argmax(ch_csd_fourier)
    assert max_ch_csd_mt == max_ch_csd_fourier

    # Check a list of CSD matrices is returned for multiple frequencies within
    # a given range when fsum=False
    with warnings.catch_warnings(record=True):  # deprecation
        csd_fsum = csd_array(X, sfreq, mode='fourier', fmin=8, fmax=20,
                             fsum=True)
        csds = csd_array(X, sfreq, mode='fourier', fmin=8, fmax=20,
                         fsum=False)

    assert csds._data.shape[1] == 2
    assert len(csds.frequencies) == 2
#.........这里部分代码省略.........

def test_csd_epochs():
    """Test computing cross-spectral density from epochs."""
    epochs = _get_real_data()

    # Check that wrong parameters are recognized
    with warnings.catch_warnings(record=True):  # deprecation
        raises(ValueError, csd_epochs, epochs, mode='notamode')
        raises(ValueError, csd_epochs, epochs, fmin=20, fmax=10)
        raises(ValueError, csd_epochs, epochs, fmin=20, fmax=20.1)
        raises(ValueError, csd_epochs, epochs, tmin=0.15, tmax=0.1)
        raises(ValueError, csd_epochs, epochs, tmin=0, tmax=10)
        raises(ValueError, csd_epochs, epochs, tmin=10, tmax=11)

    # Test deprecation warning
    with warnings.catch_warnings(record=True) as ws:
        warnings.simplefilter('always')
        csd_mt = csd_epochs(epochs, mode='multitaper', fmin=8, fmax=12,
                            tmin=0.04, tmax=0.15)
    assert len([w for w in ws
                if issubclass(w.category, DeprecationWarning)]) == 1

    with warnings.catch_warnings(record=True):  # deprecation
        csd_fourier = csd_epochs(epochs, mode='fourier', fmin=8, fmax=12,
                                 tmin=0.04, tmax=0.15)

    # Check shape of the CSD matrix
    n_chan = len(csd_mt.ch_names)
    csd_mt_data = csd_mt.get_data()
    csd_fourier_data = csd_fourier.get_data()
    assert csd_mt_data.shape == (n_chan, n_chan)
    assert csd_fourier_data.shape == (n_chan, n_chan)

    # Check if the CSD matrix is hermitian
    assert_array_equal(np.tril(csd_mt_data).T.conj(),
                       np.triu(csd_mt_data))
    assert_array_equal(np.tril(csd_fourier_data).T.conj(),
                       np.triu(csd_fourier_data))

    # Computing induced power for comparison
    epochs.crop(tmin=0.04, tmax=0.15)
    tfr = tfr_morlet(epochs, freqs=[10], n_cycles=0.6, return_itc=False)
    power = np.mean(tfr.data, 2)

    # Maximum PSD should occur for specific channel
    max_ch_power = power.argmax()
    max_ch_mt = csd_mt_data.diagonal().argmax()
    max_ch_fourier = csd_fourier_data.diagonal().argmax()
    assert max_ch_mt == max_ch_power
    assert max_ch_fourier == max_ch_power

    # Maximum CSD should occur for specific channel
    ch_csd_mt = np.abs(csd_mt_data[max_ch_power])
    ch_csd_mt[max_ch_power] = 0.
    max_ch_csd_mt = np.argmax(ch_csd_mt)
    ch_csd_fourier = np.abs(csd_fourier_data[max_ch_power])
    ch_csd_fourier[max_ch_power] = 0.
    max_ch_csd_fourier = np.argmax(ch_csd_fourier)
    assert max_ch_csd_mt == max_ch_csd_fourier

    # Check a list of CSD matrices is returned for multiple frequencies within
    # a given range when fsum=False
    with warnings.catch_warnings(record=True):  # deprecation
        csd_fsum = csd_epochs(epochs, mode='fourier', fmin=8, fmax=20,
                              fsum=True)
        csds = csd_epochs(epochs, mode='fourier', fmin=8, fmax=20, fsum=False)
    assert len(csd_fsum.frequencies) == 1
    assert len(csds.frequencies) == 2
    assert_array_equal(csd_fsum.frequencies[0], csds.frequencies)

    csd_sum = csds._data.sum(axis=1, keepdims=True)
    assert_array_equal(csd_fsum._data, csd_sum)