code: def test_csd_multitaper():
      """Test computing cross-spectral density using multitapers."""
      epochs = _generate_coherence_data()
      sfreq = epochs.info["sfreq"]
      _test_fourier_multitaper_parameters(
          epochs, csd_multitaper, csd_array_multitaper
      )
  
      # Compute CSDs using various parameters
      times = [(None, None), (1, 9)]
      as_arrays = [False, True]
      adaptives = [False, True]
      parameters = product(times, as_arrays, adaptives)
      for (tmin, tmax), as_array, adaptive in parameters:
          if as_array:
              csd = csd_array_multitaper(
                  epochs.get_data(copy=False),
                  sfreq,
                  epochs.tmin,
                  adaptive=adaptive,
                  fmin=9,
                  fmax=23,
                  tmin=tmin,
                  tmax=tmax,
                  ch_names=epochs.ch_names,
              )
          else:
              csd = csd_multitaper(
                  epochs,
                  adaptive=adaptive,
                  fmin=9,
                  fmax=23,
                  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
          csd = csd.mean([9.9, 14.9, 21.9], [10.1, 15.1, 22.1])
          _test_csd_matrix(csd)
  
      # Test equivalence with PSD
      spectrum = epochs.compute_psd(
          fmin=1e-3, normalization="full"
      )  # omit DC
      psd, psd_freqs = spectrum.get_data(return_freqs=True)
      csd = csd_multitaper(epochs)
      assert_allclose(psd_freqs, csd.frequencies)
      csd = np.array(
          [np.diag(csd.get_data(index=ii)) for ii in range(len(csd))]
      ).T
      assert_allclose(psd[0], csd)
  
      # For the next test, generate a simple sine wave with a known power
      times = np.arange(20 * sfreq) / sfreq  # 20 seconds of signal
      signal = np.sin(2 * np.pi * 10 * times)[
          None, None, :
      ]  # 10 Hz wave
      signal_power_per_sample = sum_squared(signal) / len(times)
  
      # Power per sample should not depend on time window length
      for tmax in [12, 18]:
          t_mask = times <= tmax
          n_samples = sum(t_mask)
          n_fft = len(times)
  
          # Power per sample should not depend on number of tapers
          for n_tapers in [1, 2, 5]:
              bandwidth = sfreq / float(n_samples) * (n_tapers + 1)
              csd_mt = (
                  csd_array_multitaper(
                      signal,
                      sfreq,
                      tmax=tmax,
                      bandwidth=bandwidth,
                      n_fft=n_fft,
                  )
                  .sum()
                  .get_data()
              )
              mt_power_per_sample = np.abs(csd_mt[0, 0]) * sfreq / n_fft
              assert (
                  abs(signal_power_per_sample - mt_power_per_sample)
                  < 0.001
              )
  
imports: import pickle
  from itertools import product
  from os import path
  import numpy
  import pytest
  from numpy.testing import assert_allclose, assert_array_equal
  from pytest import raises
  import mne
  from mne.channels import equalize_channels
  from mne.proj import Projection
  from mne.time_frequency import CrossSpectralDensity, csd_array_fourier, csd_array_morlet, csd_array_multitaper, csd_fourier, csd_morlet, csd_multitaper, csd_tfr, pick_channels_csd, read_csd, tfr_morlet
  from mne.time_frequency.csd import _sym_mat_to_vector, _vector_to_sym_mat
  from mne.utils import sum_squared
label_desc: Use image node_id 16 for calling a global function with example usage: test_csd_multitaper() without return types