code: def psd_array_welch(
      x,
      sfreq,
      fmin=0,
      fmax=np.inf,
      n_fft=256,
      n_overlap=0,
      n_per_seg=None,
      n_jobs=None,
      average="mean",
      window="hamming",
      remove_dc=True,
      *,
      output="power",
      verbose=None,
  ):
      """Compute power spectral density (PSD) using Welch's method.
  
      Welch's method is described in :footcite:t:`Welch1967`.
  
      Parameters
      ----------
      x : array, shape=(..., n_times)
          The data to compute PSD from.
      sfreq : float
          The sampling frequency.
      fmin : float
          The lower frequency of interest.
      fmax : float
          The upper frequency of interest.
      n_fft : int
          The length of FFT used, must be ``>= n_per_seg`` (default: 256).
          The segments will be zero-padded if ``n_fft > n_per_seg``.
      n_overlap : int
          The number of points of overlap between segments. Will be adjusted
          to be <= n_per_seg. The default value is 0.
      n_per_seg : int | None
          Length of each Welch segment (windowed with a Hamming window). Defaults
          to None, which sets n_per_seg equal to n_fft.
      %(n_jobs)s
      %(average_psd)s
  
          .. versionadded:: 0.19.0
      %(window_psd)s
  
          .. versionadded:: 0.22.0
      %(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
            window.
  
          .. versionadded:: 1.4.0
      %(verbose)s
  
      Returns
      -------
      psds : ndarray, shape (..., n_freqs) or (..., n_freqs, n_segments)
          The power spectral densities. If ``average='mean`` or
          ``average='median'``, the returned array will have the same shape
          as the input data plus an additional frequency dimension.
          If ``average=None``, the returned array will have the same shape as
          the input data plus two additional dimensions corresponding to
          frequencies and the unaggregated segments, respectively.
      freqs : ndarray, shape (n_freqs,)
          The frequencies.
  
      Notes
      -----
      .. versionadded:: 0.14.0
  
      References
      ----------
      .. footbibliography::
      """
      _check_option("average", average, (None, False, "mean", "median"))
      _check_option("output", output, ("power", "complex"))
      detrend = "constant" if remove_dc else False
      mode = "complex" if output == "complex" else "psd"
      n_fft = _ensure_int(n_fft, "n_fft")
      n_overlap = _ensure_int(n_overlap, "n_overlap")
      if n_per_seg is not None:
          n_per_seg = _ensure_int(n_per_seg, "n_per_seg")
      if average is False:
          average = None
  
      dshape = x.shape[:-1]
      n_times = x.shape[-1]
      x = x.reshape(-1, n_times)
  
      # Prep the PSD
      n_fft, n_per_seg, n_overlap = _check_nfft(
          n_times, n_fft, n_per_seg, n_overlap
      )
      win_size = n_fft / float(sfreq)
      logger.info("Effective window size : %0.3f (s)" % win_size)
      freqs = np.arange(n_fft // 2 + 1, dtype=float) * (sfreq / n_fft)
      freq_mask = (freqs >= fmin) & (freqs <= fmax)
      if not freq_mask.any():
          raise ValueError(
              f"No frequencies found between fmin={fmin} and fmax={fmax}"
          )
      freq_sl = slice(*(np.where(freq_mask)[0][[0, -1]] + [0, 1]))
      del freq_mask
      freqs = freqs[freq_sl]
  
      # Parallelize across first N-1 dimensions
      logger.debug(
          f"Spectogram using {n_fft}-point FFT on {n_per_seg} samples with "
          f"{n_overlap} overlap and {window} window"
      )
  
      parallel, my_spect_func, n_jobs = parallel_func(
          _spect_func, n_jobs=n_jobs
      )
      func = partial(
          spectrogram,
          detrend=detrend,
          noverlap=n_overlap,
          nperseg=n_per_seg,
          nfft=n_fft,
          fs=sfreq,
          window=window,
          mode=mode,
      )
      x_splits = [
          arr for arr in np.array_split(x, n_jobs) if arr.size != 0
      ]
      f_spect = parallel(
          my_spect_func(
              d,
              func=func,
              freq_sl=freq_sl,
              average=average,
              output=output,
          )
          for d in x_splits
      )
      psds = np.concatenate(f_spect, axis=0)
      shape = dshape + (len(freqs),)
      if average is None:
          shape = shape + (-1,)
      psds.shape = shape
      return psds, freqs
  
file_path: mne-python/mne/time_frequency/psd.py
imports: from functools import partial
  import numpy
  from scipy.signal import spectrogram
  from ..parallel import parallel_func
  from ..utils import _check_option, _ensure_int, logger, verbose
label_desc: Use image node_id 5 for calling a global function with example usage: psd_array_welch(x, sfreq, fmin, fmax, n_fft, n_overlap, n_per_seg, n_jobs, average, window, remove_dc) and returns: psds,  freqs