code: def plot_source_spectrogram(
      stcs,
      freq_bins,
      tmin=None,
      tmax=None,
      source_index=None,
      colorbar=False,
      show=True,
  ):
      """Plot source power in time-freqency grid.
  
      Parameters
      ----------
      stcs : list of SourceEstimate
          Source power for consecutive time windows, one SourceEstimate object
          should be provided for each frequency bin.
      freq_bins : list of tuples of float
          Start and end points of frequency bins of interest.
      tmin : float
          Minimum time instant to show.
      tmax : float
          Maximum time instant to show.
      source_index : int | None
          Index of source for which the spectrogram will be plotted. If None,
          the source with the largest activation will be selected.
      colorbar : bool
          If true, a colorbar will be added to the plot.
      show : bool
          Show figure if True.
  
      Returns
      -------
      fig : instance of Figure
          The figure.
      """
      import matplotlib.pyplot as plt
  
      # Input checks
      if len(stcs) == 0:
          raise ValueError("cannot plot spectrogram if len(stcs) == 0")
  
      stc = stcs[0]
      if tmin is not None and tmin < stc.times[0]:
          raise ValueError(
              "tmin cannot be smaller than the first time point "
              "provided in stcs"
          )
      if tmax is not None and tmax > stc.times[-1] + stc.tstep:
          raise ValueError(
              "tmax cannot be larger than the sum of the last time "
              "point and the time step, which are provided in stcs"
          )
  
      # Preparing time-frequency cell boundaries for plotting
      if tmin is None:
          tmin = stc.times[0]
      if tmax is None:
          tmax = stc.times[-1] + stc.tstep
      time_bounds = np.arange(tmin, tmax + stc.tstep, stc.tstep)
      freq_bounds = sorted(set(np.ravel(freq_bins)))
      freq_ticks = copy.deepcopy(freq_bounds)
  
      # Reject time points that will not be plotted and gather results
      source_power = []
      for stc in stcs:
          stc = stc.copy()  # copy since crop modifies inplace
          stc.crop(tmin, tmax - stc.tstep)
          source_power.append(stc.data)
      source_power = np.array(source_power)
  
      # Finding the source with maximum source power
      if source_index is None:
          source_index = np.unravel_index(
              source_power.argmax(), source_power.shape
          )[1]
  
      # If there is a gap in the frequency bins record its locations so that it
      # can be covered with a gray horizontal bar
      gap_bounds = []
      for i in range(len(freq_bins) - 1):
          lower_bound = freq_bins[i][1]
          upper_bound = freq_bins[i + 1][0]
          if lower_bound != upper_bound:
              freq_bounds.remove(lower_bound)
              gap_bounds.append((lower_bound, upper_bound))
  
      # Preparing time-frequency grid for plotting
      time_grid, freq_grid = np.meshgrid(time_bounds, freq_bounds)
  
      # Plotting the results
      fig = plt.figure(figsize=(9, 6), layout="constrained")
      plt.pcolor(
          time_grid,
          freq_grid,
          source_power[:, source_index, :],
          cmap="Reds",
      )
      ax = plt.gca()
  
      ax.set(
          title="Source power",
          xlabel="Time (s)",
          ylabel="Frequency (Hz)",
      )
  
      time_tick_labels = [str(np.round(t, 2)) for t in time_bounds]
      n_skip = 1 + len(time_bounds) // 10
      for i in range(len(time_bounds)):
          if i % n_skip != 0:
              time_tick_labels[i] = ""
  
      ax.set_xticks(time_bounds)
      ax.set_xticklabels(time_tick_labels)
      plt.xlim(time_bounds[0], time_bounds[-1])
      plt.yscale("log")
      ax.set_yticks(freq_ticks)
      ax.set_yticklabels([np.round(freq, 2) for freq in freq_ticks])
      plt.ylim(freq_bounds[0], freq_bounds[-1])
  
      plt.grid(True, ls="-")
      if colorbar:
          plt.colorbar()
  
      # Covering frequency gaps with horizontal bars
      for lower_bound, upper_bound in gap_bounds:
          plt.barh(
              lower_bound,
              time_bounds[-1] - time_bounds[0],
              upper_bound - lower_bound,
              time_bounds[0],
              color="#666666",
          )
  
      plt_show(show)
      return fig
  
imports: import copy
  import io
  import os
  import os.path
  import warnings
  from collections import defaultdict
  from glob import glob
  from itertools import cycle
  from pathlib import Path
  import numpy
  from scipy.signal import filtfilt, freqz, group_delay, lfilter, sosfilt, sosfiltfilt
  from .._fiff.constants import FIFF
  from .._fiff.pick import _DATA_CH_TYPES_SPLIT, _picks_by_type, pick_channels, pick_info, pick_types
  from .._fiff.proj import make_projector
  from .._freesurfer import _check_mri, _mri_orientation, _read_mri_info, _reorient_image
  from ..defaults import DEFAULTS
  from ..filter import estimate_ringing_samples
  from ..fixes import _safe_svd
  from ..rank import compute_rank
  from ..surface import read_surface
  from ..transforms import _frame_to_str, apply_trans
  from ..utils import _check_option, _mask_to_onsets_offsets, _on_missing, _pl, fill_doc, get_subjects_dir, logger, verbose, warn
  from .utils import _figure_agg, _get_color_list, _prepare_trellis, _validate_type, plt_show
label_desc: Use image node_id 3 for calling a global function with example usage: plot_source_spectrogram(stcs, freq_bins, tmin, tmax, source_index, colorbar, show) and returns: fig