code: def psd_array_multitaper(
      x,
      sfreq,
      fmin=0.0,
      fmax=np.inf,
      bandwidth=None,
      adaptive=False,
      low_bias=True,
      normalization="length",
      remove_dc=True,
      output="power",
      n_jobs=None,
      *,
      max_iter=150,
      verbose=None,
  ):
      r"""Compute power spectral density (PSD) using a multi-taper method.
  
      The power spectral density is computed with DPSS
      tapers :footcite:p:`Slepian1978`.
  
      Parameters
      ----------
      x : array, shape=(..., n_times)
          The data to compute PSD from.
      sfreq : float
          The sampling frequency.
      %(fmin_fmax_psd)s
      bandwidth : float
          Frequency bandwidth of the multi-taper window function in Hz. For a
          given frequency, frequencies at ``± bandwidth / 2`` are smoothed
          together. The default value is a bandwidth of
          ``8 * (sfreq / n_times)``.
      adaptive : bool
          Use adaptive weights to combine the tapered spectra into PSD
          (slow, use n_jobs >> 1 to speed up computation).
      low_bias : bool
          Only use tapers with more than 90%% spectral concentration within
          bandwidth.
      %(normalization)s
      %(remove_dc)s
      output : str
          The format of the returned ``psds`` array, ``'complex'`` or
          ``'power'``:
  
          * ``'power'`` : the power spectral density is returned.
          * ``'complex'`` : the complex fourier coefficients are returned per
            taper.
      %(n_jobs)s
      %(max_iter_multitaper)s
      %(verbose)s
  
      Returns
      -------
      psds : ndarray, shape (..., n_freqs) or (..., n_tapers, n_freqs)
          The power spectral densities. All dimensions up to the last (or the
          last two if ``output='complex'``) will be the same as input.
      freqs : array
          The frequency points in Hz of the PSD.
      weights : ndarray
          The weights used for averaging across tapers. Only returned if
          ``output='complex'``.
  
      See Also
      --------
      csd_multitaper
      mne.io.Raw.compute_psd
      mne.Epochs.compute_psd
      mne.Evoked.compute_psd
  
      Notes
      -----
      .. versionadded:: 0.14.0
  
      References
      ----------
      .. footbibliography::
      """
      _check_option("normalization", normalization, ["length", "full"])
  
      # Reshape data so its 2-D for parallelization
      ndim_in = x.ndim
      x = np.atleast_2d(x)
      n_times = x.shape[-1]
      dshape = x.shape[:-1]
      x = x.reshape(-1, n_times)
  
      dpss, eigvals, adaptive = _compute_mt_params(
          n_times, sfreq, bandwidth, low_bias, adaptive
      )
      n_tapers = len(dpss)
      weights = np.sqrt(eigvals)[np.newaxis, :, np.newaxis]
  
      # decide which frequencies to keep
      freqs = rfftfreq(n_times, 1.0 / sfreq)
      freq_mask = (freqs >= fmin) & (freqs <= fmax)
      freqs = freqs[freq_mask]
      n_freqs = len(freqs)
  
      if output == "complex":
          psd = np.zeros(
              (x.shape[0], n_tapers, n_freqs), dtype="complex"
          )
      else:
          psd = np.zeros((x.shape[0], n_freqs))
  
      # Let's go in up to 50 MB chunks of signals to save memory
      n_chunk = max(50000000 // (len(freq_mask) * len(eigvals) * 16), 1)
      offsets = np.concatenate(
          (np.arange(0, x.shape[0], n_chunk), [x.shape[0]])
      )
      for start, stop in zip(offsets[:-1], offsets[1:]):
          x_mt = _mt_spectra(
              x[start:stop], dpss, sfreq, remove_dc=remove_dc
          )[0]
          if output == "power":
              if not adaptive:
                  psd[start:stop] = _psd_from_mt(
                      x_mt[:, :, freq_mask], weights
                  )
              else:
                  (
                      parallel,
                      my_psd_from_mt_adaptive,
                      n_jobs,
                  ) = parallel_func(_psd_from_mt_adaptive, n_jobs)
                  n_splits = min(stop - start, n_jobs)
                  out = parallel(
                      my_psd_from_mt_adaptive(
                          x, eigvals, freq_mask, max_iter
                      )
                      for x in np.array_split(x_mt, n_splits)
                  )
                  psd[start:stop] = np.concatenate(out)
          else:
              psd[start:stop] = x_mt[:, :, freq_mask]
  
      if normalization == "full":
          psd /= sfreq
  
      # Combining/reshaping to original data shape
      last_dims = (
          (n_freqs,) if output == "power" else (n_tapers, n_freqs)
      )
      psd.shape = dshape + last_dims
      if ndim_in == 1:
          psd = psd[0]
  
      if output == "complex":
          return psd, freqs, weights
      else:
          return psd, freqs
  
file_path: mne-python/mne/time_frequency/multitaper.py
imports: import numpy
  from scipy.fft import rfft, rfftfreq
  from scipy.integrate import trapezoid
  from scipy.signal import get_window
  from scipy.signal.windows import dpss
  from ..parallel import parallel_func
  from ..utils import _check_option, logger, verbose, warn
label_desc: Use image node_id 7 for calling a global function with example usage: psd_array_multitaper(x, sfreq, fmin, fmax, bandwidth, adaptive, low_bias, normalization, remove_dc, output, n_jobs) and returns: psd,  freqs,  weights, psd,  freqs