code: def csd_tfr(
      epochs_tfr,
      tmin=None,
      tmax=None,
      picks=None,
      projs=None,
      verbose=None,
  ):
      """Compute covariance matrices across frequencies for TFR epochs.
  
      Parameters
      ----------
      epochs_tfr : EpochsTFR
          The time-frequency resolved epochs over which to compute the
          covariance.
      tmin : float | None
          Minimum time instant to consider, in seconds. If ``None`` start at
          first sample.
      tmax : float | None
          Maximum time instant to consider, in seconds. If ``None`` end at last
          sample.
      %(picks_good_data_noref)s
      projs : list of Projection | None
          List of projectors to store in the CSD object. Defaults to ``None``,
          which means the projectors defined in the EpochsTFR object will be
          copied.
      %(verbose)s
  
      Returns
      -------
      res : instance of CrossSpectralDensity
          Cross-spectral density restricted to selected channels.
      """
      _validate_type(epochs_tfr, EpochsTFR)
      epochs_tfr, projs = _prepare_csd(
          epochs_tfr, tmin, tmax, picks, projs
      )
      X = epochs_tfr.data
      times = epochs_tfr.times
      n_channels, n_freqs = (
          len(epochs_tfr.ch_names),
          epochs_tfr.freqs.size,
      )
      data = np.zeros(
          (n_channels * (n_channels + 1) // 2, n_freqs),
          dtype=np.complex128,
      )
  
      # Slice X to the requested time window
      tstart = (
          None if tmin is None else np.searchsorted(times, tmin - 1e-10)
      )
      tstop = (
          None if tmax is None else np.searchsorted(times, tmax + 1e-10)
      )
      X = X[:, :, :, tstart:tstop]
  
      for idx, epochs_data in enumerate(X):
          # This is equivalent to:
          # csds = np.vstack([np.mean(epochs_data[[i]] * epochs_data_conj[i:],
          #                           axis=2) for i in range(n_channels)])
          # There is a redundancy in the calculation here because we don't really
          # need the lower triangle of the matrix, but it should still be faster
          # than a loop (hopefully!).
          csds = np.einsum(
              "xft,yft->xyf", epochs_data, np.conj(epochs_data)
          )
          csds = csds[np.triu_indices(n_channels) + (slice(None),)]
          csds /= epochs_data.shape[-1]
  
          # Scaling by sampling frequency for compatibility with Matlab
          csds /= epochs_tfr.info["sfreq"]
          data += csds
  
      # scale to compute mean
      data /= len(epochs_tfr)
  
      # TO DO: EpochTFR should store n_fft to be consistent
      return CrossSpectralDensity(
          data=data,
          ch_names=epochs_tfr.ch_names,
          tmin=tmin,
          tmax=tmax,
          frequencies=epochs_tfr.freqs,
          n_fft=None,
          projs=projs,
      )
  
imports: import copy
  import numbers
  import numpy
  from scipy.fft import rfftfreq
  from .._fiff.pick import _picks_to_idx, pick_channels
  from ..parallel import parallel_func
  from ..time_frequency.multitaper import _compute_mt_params, _csd_from_mt, _mt_spectra, _psd_from_mt_adaptive
  from ..utils import ProgressBar, _check_fname, _import_h5io_funcs, _validate_type, copy_function_doc_to_method_doc, logger, verbose, warn
  from ..viz.misc import plot_csd
  from .tfr import EpochsTFR, _cwt_array, _get_nfft, morlet
label_desc: Use image node_id 18 for calling a global function with example usage: csd_tfr(epochs_tfr, tmin, tmax, picks, projs, verbose) and returns: CrossSpectralDensity