"""Container classes for spectral data."""

# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

from copy import deepcopy
from functools import partial
from inspect import signature

import numpy as np

from .._fiff.meas_info import ContainsMixin, Info
from .._fiff.pick import _pick_data_channels, _picks_to_idx, pick_info
from ..channels.channels import UpdateChannelsMixin
from ..channels.layout import _merge_ch_data, find_layout
from ..defaults import (
    _BORDER_DEFAULT,
    _EXTRAPOLATE_DEFAULT,
    _INTERPOLATION_DEFAULT,
    _handle_default,
)
from ..html_templates import _get_html_template
from ..utils import (
    GetEpochsMixin,
    _build_data_frame,
    _check_method_kwargs,
    _check_pandas_index_arguments,
    _check_pandas_installed,
    _check_sphere,
    _time_mask,
    _validate_type,
    fill_doc,
    legacy,
    logger,
    object_diff,
    repr_html,
    verbose,
    warn,
)
from ..utils.check import (
    _check_fname,
    _check_option,
    _import_h5io_funcs,
    _is_numeric,
    check_fname,
)
from ..utils.misc import _pl
from ..utils.spectrum import _get_instance_type_string, _split_psd_kwargs
from ..viz.topo import _plot_timeseries, _plot_timeseries_unified, _plot_topo
from ..viz.topomap import _make_head_outlines, _prepare_topomap_plot, plot_psds_topomap
from ..viz.utils import (
    _format_units_psd,
    _get_plot_ch_type,
    _make_combine_callable,
    _plot_psd,
    _prepare_sensor_names,
    plt_show,
)
from .multitaper import _psd_from_mt, psd_array_multitaper
from .psd import _check_nfft, psd_array_welch


class SpectrumMixin:
    """Mixin providing spectral plotting methods to sensor-space containers."""

    @legacy(alt=".compute_psd().plot()")
    @verbose
    def plot_psd(
        self,
        fmin=0,
        fmax=np.inf,
        tmin=None,
        tmax=None,
        picks=None,
        proj=False,
        reject_by_annotation=True,
        *,
        method="auto",
        average=False,
        dB=True,
        estimate="power",
        xscale="linear",
        area_mode="std",
        area_alpha=0.33,
        color="black",
        line_alpha=None,
        spatial_colors=True,
        sphere=None,
        exclude="bads",
        ax=None,
        show=True,
        n_jobs=1,
        verbose=None,
        **method_kw,
    ):
        """%(plot_psd_doc)s.

        Parameters
        ----------
        %(fmin_fmax_psd)s
        %(tmin_tmax_psd)s
        %(picks_good_data_noref)s
        %(proj_psd)s
        %(reject_by_annotation_psd)s
        %(method_plot_psd_auto)s
        %(average_plot_psd)s
        %(dB_plot_psd)s
        %(estimate_plot_psd)s
        %(xscale_plot_psd)s
        %(area_mode_plot_psd)s
        %(area_alpha_plot_psd)s
        %(color_plot_psd)s
        %(line_alpha_plot_psd)s
        %(spatial_colors_psd)s
        %(sphere_topomap_auto)s

            .. versionadded:: 0.22.0
        exclude : list of str | 'bads'
            Channels names to exclude from being shown. If 'bads', the bad
            channels are excluded. Pass an empty list to plot all channels
            (including channels marked "bad", if any).

            .. versionadded:: 0.24.0
        %(ax_plot_psd)s
        %(show)s
        %(n_jobs)s
        %(verbose)s
        %(method_kw_psd)s

        Returns
        -------
        fig : instance of Figure
            Figure with frequency spectra of the data channels.

        Notes
        -----
        %(notes_plot_psd_meth)s
        """
        init_kw, plot_kw = _split_psd_kwargs(plot_fun=Spectrum.plot)
        return self.compute_psd(**init_kw).plot(**plot_kw)

    @legacy(alt=".compute_psd().plot_topo()")
    @verbose
    def plot_psd_topo(
        self,
        tmin=None,
        tmax=None,
        fmin=0,
        fmax=100,
        proj=False,
        *,
        method="auto",
        dB=True,
        layout=None,
        color="w",
        fig_facecolor="k",
        axis_facecolor="k",
        axes=None,
        block=False,
        show=True,
        n_jobs=None,
        verbose=None,
        **method_kw,
    ):
        """Plot power spectral density, separately for each channel.

        Parameters
        ----------
        %(tmin_tmax_psd)s
        %(fmin_fmax_psd_topo)s
        %(proj_psd)s
        %(method_plot_psd_auto)s
        %(dB_spectrum_plot_topo)s
        %(layout_spectrum_plot_topo)s
        %(color_spectrum_plot_topo)s
        %(fig_facecolor)s
        %(axis_facecolor)s
        %(axes_spectrum_plot_topo)s
        %(block)s
        %(show)s
        %(n_jobs)s
        %(verbose)s
        %(method_kw_psd)s Defaults to ``dict(n_fft=2048)``.

        Returns
        -------
        fig : instance of matplotlib.figure.Figure
            Figure distributing one image per channel across sensor topography.
        """
        init_kw, plot_kw = _split_psd_kwargs(plot_fun=Spectrum.plot_topo)
        return self.compute_psd(**init_kw).plot_topo(**plot_kw)

    @legacy(alt=".compute_psd().plot_topomap()")
    @verbose
    def plot_psd_topomap(
        self,
        bands=None,
        tmin=None,
        tmax=None,
        ch_type=None,
        *,
        proj=False,
        method="auto",
        normalize=False,
        agg_fun=None,
        dB=False,
        sensors=True,
        show_names=False,
        mask=None,
        mask_params=None,
        contours=0,
        outlines="head",
        sphere=None,
        image_interp=_INTERPOLATION_DEFAULT,
        extrapolate=_EXTRAPOLATE_DEFAULT,
        border=_BORDER_DEFAULT,
        res=64,
        size=1,
        cmap=None,
        vlim=(None, None),
        cnorm=None,
        colorbar=True,
        cbar_fmt="auto",
        units=None,
        axes=None,
        show=True,
        n_jobs=None,
        verbose=None,
        **method_kw,
    ):
        """Plot scalp topography of PSD for chosen frequency bands.

        Parameters
        ----------
        %(bands_psd_topo)s
        %(tmin_tmax_psd)s
        %(ch_type_topomap_psd)s
        %(proj_psd)s
        %(method_plot_psd_auto)s
        %(normalize_psd_topo)s
        %(agg_fun_psd_topo)s
        %(dB_plot_topomap)s
        %(sensors_topomap)s
        %(show_names_topomap)s
        %(mask_evoked_topomap)s
        %(mask_params_topomap)s
        %(contours_topomap)s
        %(outlines_topomap)s
        %(sphere_topomap_auto)s
        %(image_interp_topomap)s
        %(extrapolate_topomap)s
        %(border_topomap)s
        %(res_topomap)s
        %(size_topomap)s
        %(cmap_topomap)s
        %(vlim_plot_topomap_psd)s
        %(cnorm)s

            .. versionadded:: 1.2
        %(colorbar_topomap)s
        %(cbar_fmt_topomap_psd)s
        %(units_topomap)s
        %(axes_spectrum_plot_topomap)s
        %(show)s
        %(n_jobs)s
        %(verbose)s
        %(method_kw_psd)s

        Returns
        -------
        fig : instance of Figure
            Figure showing one scalp topography per frequency band.
        """
        init_kw, plot_kw = _split_psd_kwargs(plot_fun=Spectrum.plot_topomap)
        return self.compute_psd(**init_kw).plot_topomap(**plot_kw)

    def _set_legacy_nfft_default(self, tmin, tmax, method, method_kw):
        """Update method_kw with legacy n_fft default for plot_psd[_topo]().

        This method returns ``None`` and has a side effect of (maybe) updating
        the ``method_kw`` dict.
        """
        if method == "welch" and method_kw.get("n_fft") is None:
            tm = _time_mask(self.times, tmin, tmax, sfreq=self.info["sfreq"])
            method_kw["n_fft"] = min(np.sum(tm), 2048)


class BaseSpectrum(ContainsMixin, UpdateChannelsMixin):
    """Base class for Spectrum and EpochsSpectrum."""

    def __init__(
        self,
        inst,
        method,
        fmin,
        fmax,
        tmin,
        tmax,
        picks,
        exclude,
        proj,
        remove_dc,
        *,
        n_jobs,
        verbose=None,
        **method_kw,
    ):
        # arg checking
        self._sfreq = inst.info["sfreq"]
        if np.isfinite(fmax) and (fmax > self.sfreq / 2):
            raise ValueError(
                f"Requested fmax ({fmax} Hz) must not exceed ½ the sampling "
                f'frequency of the data ({0.5 * inst.info["sfreq"]} Hz).'
            )
        # method
        self._inst_type = type(inst)
        method = _validate_method(method, _get_instance_type_string(self))
        psd_funcs = dict(welch=psd_array_welch, multitaper=psd_array_multitaper)
        # triage method and kwargs. partial() doesn't check validity of kwargs,
        # so we do it manually to save compute time if any are invalid.
        psd_funcs = dict(welch=psd_array_welch, multitaper=psd_array_multitaper)
        _check_method_kwargs(psd_funcs[method], method_kw, msg=f'PSD method "{method}"')
        self._psd_func = partial(psd_funcs[method], remove_dc=remove_dc, **method_kw)

        # apply proj if desired
        if proj:
            inst = inst.copy().apply_proj()
        self.inst = inst

        # prep times and picks
        self._time_mask = _time_mask(inst.times, tmin, tmax, sfreq=self.sfreq)
        self._picks = _picks_to_idx(
            inst.info, picks, "data", exclude, with_ref_meg=False
        )

        # add the info object. bads and non-data channels were dropped by
        # _picks_to_idx() so we update the info accordingly:
        self.info = pick_info(inst.info, sel=self._picks, copy=True)

        # assign some attributes
        self.preload = True  # needed for __getitem__, never False
        self._method = method
        # self._dims may also get updated by child classes
        self._dims = (
            "channel",
            "freq",
        )
        if method_kw.get("average", "") in (None, False):
            self._dims += ("segment",)
        if self._returns_complex_tapers(**method_kw):
            self._dims = self._dims[:-1] + ("taper",) + self._dims[-1:]
        # record data type (for repr and html_repr)
        self._data_type = (
            "Fourier Coefficients"
            if method_kw.get("output") == "complex"
            else "Power Spectrum"
        )
        # set nave (child constructor overrides this for Evoked input)
        self._nave = None

    def __eq__(self, other):
        """Test equivalence of two Spectrum instances."""
        return object_diff(vars(self), vars(other)) == ""

    def __getstate__(self):
        """Prepare object for serialization."""
        inst_type_str = _get_instance_type_string(self)
        out = dict(
            method=self.method,
            data=self._data,
            sfreq=self.sfreq,
            dims=self._dims,
            freqs=self.freqs,
            inst_type_str=inst_type_str,
            data_type=self._data_type,
            info=self.info,
            nave=self.nave,
            weights=self.weights,
        )
        return out

    def __setstate__(self, state):
        """Unpack from serialized format."""
        from ..epochs import Epochs
        from ..evoked import Evoked
        from ..io import Raw

        self._method = state["method"]
        self._data = state["data"]
        self._freqs = state["freqs"]
        self._dims = state["dims"]
        self._sfreq = state["sfreq"]
        self.info = Info(**state["info"])
        self._data_type = state["data_type"]
        self._nave = state.get("nave")  # objs saved before #11282 won't have `nave`
        self._weights = state.get("weights")  # objs saved before #12747 won't have
        self.preload = True
        # instance type
        inst_types = dict(Raw=Raw, Epochs=Epochs, Evoked=Evoked, Array=np.ndarray)
        self._inst_type = inst_types[state["inst_type_str"]]

    def __repr__(self):
        """Build string representation of the Spectrum object."""
        inst_type_str = _get_instance_type_string(self)
        # shape & dimension names
        dims = " × ".join(
            [f"{dim[0]} {dim[1]}s" for dim in zip(self.shape, self._dims)]
        )
        freq_range = f"{self.freqs[0]:0.1f}-{self.freqs[-1]:0.1f} Hz"
        return (
            f"<{self._data_type} (from {inst_type_str}, "
            f"{self.method} method) | {dims}, {freq_range}>"
        )

    @repr_html
    def _repr_html_(self, caption=None):
        """Build HTML representation of the Spectrum object."""
        inst_type_str = _get_instance_type_string(self)
        units = [f"{ch_type}: {unit}" for ch_type, unit in self.units().items()]
        t = _get_html_template("repr", "spectrum.html.jinja")
        t = t.render(spectrum=self, inst_type=inst_type_str, units=units)
        return t

    def _check_values(self):
        """Check PSD results for correct shape and bad values."""
        assert len(self._dims) == self._data.ndim, (self._dims, self._data.ndim)
        assert self._data.shape == self._shape
        # TODO: should this be more fine-grained (report "chan X in epoch Y")?
        ch_dim = self._dims.index("channel")
        dims = list(range(self._data.ndim))
        dims.pop(ch_dim)
        # take min() across all but the channel axis
        # (if the abs becomes memory intensive we could iterate over channels)
        use_data = self._data
        if use_data.dtype.kind == "c":
            use_data = np.abs(use_data)
        bad_value = use_data.min(axis=tuple(dims)) == 0
        bad_value &= ~np.isin(self.ch_names, self.info["bads"])
        if bad_value.any():
            chs = np.array(self.ch_names)[bad_value].tolist()
            s = _pl(bad_value.sum())
            warn(f'Zero value in spectrum for channel{s} {", ".join(chs)}', UserWarning)

    def _returns_complex_tapers(self, **method_kw):
        return self.method == "multitaper" and method_kw.get("output") == "complex"

    def _compute_spectra(self, data, fmin, fmax, n_jobs, method_kw, verbose):
        # make the spectra
        result = self._psd_func(
            data, self.sfreq, fmin=fmin, fmax=fmax, n_jobs=n_jobs, verbose=verbose
        )
        # assign ._data (handling unaggregated multitaper output)
        if self._returns_complex_tapers(**method_kw):
            fourier_coefs, freqs, weights = result
            self._data = fourier_coefs
            self._weights = weights
        else:
            psds, freqs = result
            self._data = psds
            self._weights = None
        # assign properties (._data already assigned above)
        self._freqs = freqs
        # this is *expected* shape, it gets asserted later in _check_values()
        # (and then deleted afterwards)
        self._shape = (len(self.ch_names), len(self.freqs))
        # append n_welch_segments (use "" as .get() default since None considered valid)
        if method_kw.get("average", "") in (None, False):
            n_welch_segments = _compute_n_welch_segments(data.shape[-1], method_kw)
            self._shape += (n_welch_segments,)
        # insert n_tapers
        if self._returns_complex_tapers(**method_kw):
            self._shape = self._shape[:-1] + (self._weights.size,) + self._shape[-1:]
        # we don't need these anymore, and they make save/load harder
        del self._picks
        del self._psd_func
        del self._time_mask

    @property
    def _detrend_picks(self):
        """Provide compatibility with __iter__."""
        return list()

    @property
    def ch_names(self):
        return self.info["ch_names"]

    @property
    def data(self):
        return self._data

    @property
    def freqs(self):
        return self._freqs

    @property
    def method(self):
        return self._method

    @property
    def nave(self):
        return self._nave

    @property
    def weights(self):
        return self._weights

    @property
    def sfreq(self):
        return self._sfreq

    @property
    def shape(self):
        return self._data.shape

    def copy(self):
        """Return copy of the Spectrum instance.

        Returns
        -------
        spectrum : instance of Spectrum
            A copy of the object.
        """
        return deepcopy(self)

    @fill_doc
    def get_data(
        self, picks=None, exclude="bads", fmin=0, fmax=np.inf, return_freqs=False
    ):
        """Get spectrum data in NumPy array format.

        Parameters
        ----------
        %(picks_good_data_noref)s
        %(exclude_spectrum_get_data)s
        %(fmin_fmax_psd)s
        return_freqs : bool
            Whether to return the frequency bin values for the requested
            frequency range. Default is ``False``.

        Returns
        -------
        data : array
            The requested data in a NumPy array.
        freqs : array
            The frequency values for the requested range. Only returned if
            ``return_freqs`` is ``True``.
        """
        picks = _picks_to_idx(
            self.info, picks, "data_or_ica", exclude=exclude, with_ref_meg=False
        )
        fmin_idx = np.searchsorted(self.freqs, fmin)
        fmax_idx = np.searchsorted(self.freqs, fmax, side="right")
        freq_picks = np.arange(fmin_idx, fmax_idx)
        freq_axis = self._dims.index("freq")
        chan_axis = self._dims.index("channel")
        # normally there's a risk of np.take reducing array dimension if there
        # were only one channel or frequency selected, but `_picks_to_idx`
        # always returns an array of picks, and np.arange always returns an
        # array of freq bin indices, so we're safe; the result will always be
        # 2D.
        data = self._data.take(picks, chan_axis).take(freq_picks, freq_axis)
        if return_freqs:
            freqs = self._freqs[fmin_idx:fmax_idx]
            return (data, freqs)
        return data

    @fill_doc
    def plot(
        self,
        *,
        picks=None,
        average=False,
        dB=True,
        amplitude=False,
        xscale="linear",
        ci="sd",
        ci_alpha=0.3,
        color="black",
        alpha=None,
        spatial_colors=True,
        sphere=None,
        exclude=(),
        axes=None,
        show=True,
    ):
        """%(plot_psd_doc)s.

        Parameters
        ----------
        %(picks_all_data_noref)s

            .. versionchanged:: 1.5
                In version 1.5, the default behavior changed so that all
                :term:`data channels` (not just "good" data channels) are shown by
                default.
        average : bool
            Whether to average across channels before plotting. If ``True``, interactive
            plotting of scalp topography is disabled, and parameters ``ci`` and
            ``ci_alpha`` control the style of the confidence band around the mean.
            Default is ``False``.
        %(dB_spectrum_plot)s
        amplitude : bool
            Whether to plot an amplitude spectrum (``True``) or power spectrum
            (``False``).

                .. versionchanged:: 1.8
                    In version 1.8, the default changed to ``amplitude=False``.
        %(xscale_plot_psd)s
        ci : float | 'sd' | 'range' | None
            Type of confidence band drawn around the mean when ``average=True``. If
            ``'sd'`` the band spans ±1 standard deviation across channels. If
            ``'range'`` the band spans the range across channels at each frequency. If a
            :class:`float`, it indicates the (bootstrapped) confidence interval to
            display, and must satisfy ``0 < ci <= 100``. If ``None``, no band is drawn.
            Default is ``sd``.
        ci_alpha : float
            Opacity of the confidence band. Must satisfy ``0 <= ci_alpha <= 1``. Default
            is 0.3.
        %(color_plot_psd)s
        alpha : float | None
            Opacity of the spectrum line(s). If :class:`float`, must satisfy
            ``0 <= alpha <= 1``. If ``None``, opacity will be ``1`` when
            ``average=True`` and ``0.1`` when ``average=False``. Default is ``None``.
        %(spatial_colors_psd)s
        %(sphere_topomap_auto)s
        %(exclude_spectrum_plot)s

            .. versionchanged:: 1.5
                In version 1.5, the default behavior changed from ``exclude='bads'`` to
                ``exclude=()``.
        %(axes_spectrum_plot_topomap)s
        %(show)s

        Returns
        -------
        fig : instance of matplotlib.figure.Figure
            Figure with spectra plotted in separate subplots for each channel type.
        """
        # Must nest this _mpl_figure import because of the BACKEND global
        # stuff
        from ..viz._mpl_figure import _line_figure, _split_picks_by_type

        # arg checking
        ci = _check_ci(ci)
        _check_option("xscale", xscale, ("log", "linear"))
        sphere = _check_sphere(sphere, self.info)
        # defaults
        scalings = _handle_default("scalings", None)
        titles = _handle_default("titles", None)
        units = _handle_default("units", None)

        _validate_type(amplitude, bool, "amplitude")
        estimate = "amplitude" if amplitude else "power"

        logger.info(f"Plotting {estimate} spectral density ({dB=}).")

        # split picks by channel type
        picks = _picks_to_idx(
            self.info, picks, "data", exclude=exclude, with_ref_meg=False
        )
        (picks_list, units_list, scalings_list, titles_list) = _split_picks_by_type(
            self, picks, units, scalings, titles
        )
        # prepare data (e.g. aggregate across dims, convert complex to power)
        psd_list = [
            self._prepare_data_for_plot(
                self._data.take(_p, axis=self._dims.index("channel"))
            )
            for _p in picks_list
        ]
        # initialize figure
        fig, axes = _line_figure(self, axes, picks=picks)
        # don't add ylabels & titles if figure has unexpected number of axes
        make_label = len(axes) == len(fig.axes)
        # Plot Frequency [Hz] xlabel only on the last axis
        xlabels_list = [False] * (len(axes) - 1) + [True]
        # plot
        _plot_psd(
            self,
            fig,
            self.freqs,
            psd_list,
            picks_list,
            titles_list,
            units_list,
            scalings_list,
            axes,
            make_label,
            color,
            area_mode=ci,
            area_alpha=ci_alpha,
            dB=dB,
            estimate=estimate,
            average=average,
            spatial_colors=spatial_colors,
            xscale=xscale,
            line_alpha=alpha,
            sphere=sphere,
            xlabels_list=xlabels_list,
        )
        plt_show(show, fig)
        return fig

    @fill_doc
    def plot_topo(
        self,
        *,
        dB=True,
        layout=None,
        color="w",
        fig_facecolor="k",
        axis_facecolor="k",
        axes=None,
        block=False,
        show=True,
    ):
        """Plot power spectral density, separately for each channel.

        Parameters
        ----------
        %(dB_spectrum_plot_topo)s
        %(layout_spectrum_plot_topo)s
        %(color_spectrum_plot_topo)s
        %(fig_facecolor)s
        %(axis_facecolor)s
        %(axes_spectrum_plot_topo)s
        %(block)s
        %(show)s

        Returns
        -------
        fig : instance of matplotlib.figure.Figure
            Figure distributing one image per channel across sensor topography.
        """
        if layout is None:
            layout = find_layout(self.info)

        psds, freqs = self.get_data(return_freqs=True)
        # prepare data (e.g. aggregate across dims, convert complex to power)
        psds = self._prepare_data_for_plot(psds)
        if dB:
            psds = 10 * np.log10(psds)
            y_label = "dB"
        else:
            y_label = "Power"
        show_func = partial(
            _plot_timeseries_unified, data=[psds], color=color, times=[freqs]
        )
        click_func = partial(_plot_timeseries, data=[psds], color=color, times=[freqs])
        picks = _pick_data_channels(self.info)
        info = pick_info(self.info, picks)
        fig = _plot_topo(
            info,
            times=freqs,
            show_func=show_func,
            click_func=click_func,
            layout=layout,
            axis_facecolor=axis_facecolor,
            fig_facecolor=fig_facecolor,
            x_label="Frequency (Hz)",
            unified=True,
            y_label=y_label,
            axes=axes,
        )
        plt_show(show, block=block)
        return fig

    @fill_doc
    def plot_topomap(
        self,
        bands=None,
        ch_type=None,
        *,
        normalize=False,
        agg_fun=None,
        dB=False,
        sensors=True,
        show_names=False,
        mask=None,
        mask_params=None,
        contours=6,
        outlines="head",
        sphere=None,
        image_interp=_INTERPOLATION_DEFAULT,
        extrapolate=_EXTRAPOLATE_DEFAULT,
        border=_BORDER_DEFAULT,
        res=64,
        size=1,
        cmap=None,
        vlim=(None, None),
        cnorm=None,
        colorbar=True,
        cbar_fmt="auto",
        units=None,
        axes=None,
        show=True,
    ):
        """Plot scalp topography of PSD for chosen frequency bands.

        Parameters
        ----------
        %(bands_psd_topo)s
        %(ch_type_topomap_psd)s
        %(normalize_psd_topo)s
        %(agg_fun_psd_topo)s
        %(dB_plot_topomap)s
        %(sensors_topomap)s
        %(show_names_topomap)s
        %(mask_evoked_topomap)s
        %(mask_params_topomap)s
        %(contours_topomap)s
        %(outlines_topomap)s
        %(sphere_topomap_auto)s
        %(image_interp_topomap)s
        %(extrapolate_topomap)s
        %(border_topomap)s
        %(res_topomap)s
        %(size_topomap)s
        %(cmap_topomap)s
        %(vlim_plot_topomap_psd)s
        %(cnorm)s
        %(colorbar_topomap)s
        %(cbar_fmt_topomap_psd)s
        %(units_topomap)s
        %(axes_spectrum_plot_topomap)s
        %(show)s

        Returns
        -------
        fig : instance of Figure
            Figure showing one scalp topography per frequency band.
        """
        ch_type = _get_plot_ch_type(self, ch_type)
        if units is None:
            units = _handle_default("units", None)
        unit = units[ch_type] if hasattr(units, "keys") else units
        scalings = _handle_default("scalings", None)
        scaling = scalings[ch_type]

        (
            picks,
            pos,
            merge_channels,
            names,
            ch_type,
            sphere,
            clip_origin,
        ) = _prepare_topomap_plot(self, ch_type, sphere=sphere)
        outlines = _make_head_outlines(sphere, pos, outlines, clip_origin)

        psds, freqs = self.get_data(picks=picks, return_freqs=True)
        # prepare data (e.g. aggregate across dims, convert complex to power)
        psds = self._prepare_data_for_plot(psds)
        psds *= scaling**2

        if merge_channels:
            psds, names = _merge_ch_data(psds, ch_type, names, method="mean")

        names = _prepare_sensor_names(names, show_names)
        return plot_psds_topomap(
            psds=psds,
            freqs=freqs,
            pos=pos,
            bands=bands,
            ch_type=ch_type,
            normalize=normalize,
            agg_fun=agg_fun,
            dB=dB,
            sensors=sensors,
            names=names,
            mask=mask,
            mask_params=mask_params,
            contours=contours,
            outlines=outlines,
            sphere=sphere,
            image_interp=image_interp,
            extrapolate=extrapolate,
            border=border,
            res=res,
            size=size,
            cmap=cmap,
            vlim=vlim,
            cnorm=cnorm,
            colorbar=colorbar,
            cbar_fmt=cbar_fmt,
            unit=unit,
            axes=axes,
            show=show,
        )

    def _prepare_data_for_plot(self, data):
        # handle unaggregated Welch
        if "segment" in self._dims:
            logger.info("Aggregating Welch estimates (median) before plotting...")
            data = np.nanmedian(data, axis=self._dims.index("segment"))
        # handle unaggregated multitaper (also handles complex -> power)
        elif "taper" in self._dims:
            logger.info("Aggregating multitaper estimates before plotting...")
            data = _psd_from_mt(data, self.weights)

        # handle complex data (should only be Welch remaining)
        if np.iscomplexobj(data):
            data = (data * data.conj()).real  # Scaling may be slightly off

        # handle epochs
        if "epoch" in self._dims:
            # XXX TODO FIXME decide how to properly aggregate across repeated
            # measures (epochs) and non-repeated but correlated measures
            # (channels) when calculating stddev or a CI. For across-channel
            # aggregation, doi:10.1007/s10162-012-0321-8 used hotellings T**2
            # with a correction factor that estimated data rank using monte
            # carlo simulations; seems like we could use our own data rank
            # estimation methods to similar effect. Their exact approach used
            # complex spectra though, here we've already converted to power;
            # not sure if that makes an important difference? Anyway that
            # aggregation would need to happen in the _plot_psd function
            # though, not here... for now we just average like we always did.

            # only log message if averaging will actually have an effect
            if data.shape[0] > 1:
                logger.info("Averaging across epochs before plotting...")
            # epoch axis should always be the first axis
            data = data.mean(axis=0)

        return data

    @verbose
    def save(self, fname, *, overwrite=False, verbose=None):
        """Save spectrum data to disk (in HDF5 format).

        Parameters
        ----------
        fname : path-like
            Path of file to save to.
        %(overwrite)s
        %(verbose)s

        See Also
        --------
        mne.time_frequency.read_spectrum
        """
        _, write_hdf5 = _import_h5io_funcs()
        check_fname(fname, "spectrum", (".h5", ".hdf5"))
        fname = _check_fname(fname, overwrite=overwrite, verbose=verbose)
        out = self.__getstate__()
        write_hdf5(fname, out, overwrite=overwrite, title="mnepython")

    @verbose
    def to_data_frame(
        self, picks=None, index=None, copy=True, long_format=False, *, verbose=None
    ):
        """Export data in tabular structure as a pandas DataFrame.

        Channels are converted to columns in the DataFrame. By default,
        an additional column "freq" is added, unless ``index='freq'``
        (in which case frequency values form the DataFrame's index).

        Parameters
        ----------
        %(picks_all)s
        index : str | list of str | None
            Kind of index to use for the DataFrame. If ``None``, a sequential
            integer index (:class:`pandas.RangeIndex`) will be used. If a
            :class:`str`, a :class:`pandas.Index` will be used (see Notes). If
            a list of two or more string values, a :class:`pandas.MultiIndex`
            will be used. Defaults to ``None``.
        %(copy_df)s
        %(long_format_df_spe)s
        %(verbose)s

        Returns
        -------
        %(df_return)s

        Notes
        -----
        Valid values for ``index`` depend on whether the Spectrum was created
        from continuous data (:class:`~mne.io.Raw`, :class:`~mne.Evoked`) or
        discontinuous data (:class:`~mne.Epochs`). For continuous data, only
        ``None`` or ``'freq'`` is supported. For discontinuous data, additional
        valid values are ``'epoch'`` and ``'condition'``, or a :class:`list`
        comprising some of the valid string values (e.g.,
        ``['freq', 'epoch']``).
        """
        # check pandas once here, instead of in each private utils function
        pd = _check_pandas_installed()  # noqa
        # triage for Epoch-derived or unaggregated spectra
        from_epo = _get_instance_type_string(self) == "Epochs"
        unagg_welch = "segment" in self._dims
        unagg_mt = "taper" in self._dims
        # arg checking
        valid_index_args = ["freq"]
        if from_epo:
            valid_index_args += ["epoch", "condition"]
        index = _check_pandas_index_arguments(index, valid_index_args)
        # get data
        picks = _picks_to_idx(self.info, picks, "all", exclude=())
        data = self.get_data(picks)
        if copy:
            data = data.copy()
        # reshape
        if unagg_mt:
            data = np.moveaxis(data, self._dims.index("freq"), -2)
        if from_epo:
            n_epochs, n_picks, n_freqs = data.shape[:3]
        else:
            n_epochs, n_picks, n_freqs = (1,) + data.shape[:2]
        n_segs = data.shape[-1] if unagg_mt or unagg_welch else 1
        data = np.moveaxis(data, self._dims.index("channel"), -1)
        # at this point, should be ([epoch], freq, [segment/taper], channel)
        data = data.reshape(n_epochs * n_freqs * n_segs, n_picks)
        # prepare extra columns / multiindex
        mindex = list()
        default_index = list()
        if from_epo:
            rev_event_id = {v: k for k, v in self.event_id.items()}
            _conds = [rev_event_id[k] for k in self.events[:, 2]]
            conditions = np.repeat(_conds, n_freqs * n_segs)
            epoch_nums = np.repeat(self.selection, n_freqs * n_segs)
            mindex.extend([("condition", conditions), ("epoch", epoch_nums)])
            default_index.extend(["condition", "epoch"])
        freqs = np.tile(np.repeat(self.freqs, n_segs), n_epochs)
        mindex.append(("freq", freqs))
        default_index.append("freq")
        if unagg_mt or unagg_welch:
            name = "taper" if unagg_mt else "segment"
            seg_nums = np.tile(np.arange(n_segs), n_epochs * n_freqs)
            mindex.append((name, seg_nums))
            default_index.append(name)
        # build DataFrame
        df = _build_data_frame(
            self, data, picks, long_format, mindex, index, default_index=default_index
        )
        return df

    def units(self, latex=False):
        """Get the spectrum units for each channel type.

        Parameters
        ----------
        latex : bool
            Whether to format the unit strings as LaTeX. Default is ``False``.

        Returns
        -------
        units : dict
            Mapping from channel type to a string representation of the units
            for that channel type.
        """
        units = _handle_default("si_units", None)
        return {
            ch_type: _format_units_psd(units[ch_type], power=True, latex=latex)
            for ch_type in sorted(self.get_channel_types(unique=True))
        }


@fill_doc
class Spectrum(BaseSpectrum):
    """Data object for spectral representations of continuous data.

    .. warning:: The preferred means of creating Spectrum objects from
                 continuous or averaged data is via the instance methods
                 :meth:`mne.io.Raw.compute_psd` or
                 :meth:`mne.Evoked.compute_psd`. Direct class instantiation
                 is not supported.

    Parameters
    ----------
    inst : instance of Raw or Evoked
        The data from which to compute the frequency spectrum.
    %(method_psd_auto)s
        ``'auto'`` (default) uses Welch's method for continuous data
        and multitaper for :class:`~mne.Evoked` data.
    %(fmin_fmax_psd)s
    %(tmin_tmax_psd)s
    %(picks_good_data_noref)s
    %(exclude_psd)s
    %(proj_psd)s
    %(remove_dc)s
    %(reject_by_annotation_psd)s
    %(n_jobs)s
    %(verbose)s
    %(method_kw_psd)s

    Attributes
    ----------
    ch_names : list
        The channel names.
    freqs : array
        Frequencies at which the amplitude, power, or fourier coefficients
        have been computed.
    %(info_not_none)s
    method : ``'welch'``| ``'multitaper'``
        The method used to compute the spectrum.
    nave : int | None
        The number of trials averaged together when generating the spectrum. ``None``
        indicates no averaging is known to have occurred.
    weights : array | None
        The weights for each taper. Only present if spectra computed with
        ``method='multitaper'`` and ``output='complex'``.

        .. versionadded:: 1.8

    See Also
    --------
    EpochsSpectrum
    SpectrumArray
    mne.io.Raw.compute_psd
    mne.Epochs.compute_psd
    mne.Evoked.compute_psd

    References
    ----------
    .. footbibliography::
    """

    def __init__(
        self,
        inst,
        method,
        fmin,
        fmax,
        tmin,
        tmax,
        picks,
        exclude,
        proj,
        remove_dc,
        reject_by_annotation,
        *,
        n_jobs,
        verbose=None,
        **method_kw,
    ):
        from ..io import BaseRaw

        # triage reading from file
        if isinstance(inst, dict):
            self.__setstate__(inst)
            return
        # do the basic setup
        super().__init__(
            inst,
            method,
            fmin,
            fmax,
            tmin,
            tmax,
            picks,
            exclude,
            proj,
            remove_dc,
            n_jobs=n_jobs,
            verbose=verbose,
            **method_kw,
        )
        # get just the data we want
        if isinstance(self.inst, BaseRaw):
            start, stop = np.where(self._time_mask)[0][[0, -1]]
            rba = "NaN" if reject_by_annotation else None
            data = self.inst.get_data(
                self._picks, start, stop + 1, reject_by_annotation=rba
            )
            if np.any(np.isnan(data)) and method == "multitaper":
                raise NotImplementedError(
                    'Cannot use method="multitaper" when reject_by_annotation=True. '
                    'Please use method="welch" instead.'
                )

        else:  # Evoked
            data = self.inst.data[self._picks][:, self._time_mask]
        # set nave
        self._nave = getattr(inst, "nave", None)
        # compute the spectra
        self._compute_spectra(data, fmin, fmax, n_jobs, method_kw, verbose)
        # check for correct shape and bad values
        self._check_values()
        del self._shape  # calculated from self._data henceforth
        # save memory
        del self.inst

    def __getitem__(self, item):
        """Get Spectrum data.

        Parameters
        ----------
        item : int | slice | array-like
            Indexing is similar to a :class:`NumPy array<numpy.ndarray>`; see
            Notes.

        Returns
        -------
        %(getitem_spectrum_return)s

        Notes
        -----
        Integer-, list-, and slice-based indexing is possible:

        - ``spectrum[0]`` gives all frequency bins in the first channel
        - ``spectrum[:3]`` gives all frequency bins in the first 3 channels
        - ``spectrum[[0, 2], 5]`` gives the value in the sixth frequency bin of
          the first and third channels
        - ``spectrum[(4, 7)]`` is the same as ``spectrum[4, 7]``.

        .. note::

           Unlike :class:`~mne.io.Raw` objects (which returns a tuple of the
           requested data values and the corresponding times), accessing
           :class:`~mne.time_frequency.Spectrum` values via subscript does
           **not** return the corresponding frequency bin values. If you need
           them, use ``spectrum.freqs[freq_indices]`` or
           ``spectrum.get_data(..., return_freqs=True)``.
        """
        from ..io import BaseRaw

        self._parse_get_set_params = partial(BaseRaw._parse_get_set_params, self)
        return BaseRaw._getitem(self, item, return_times=False)


def _check_data_shape(data, info, freqs, dim_names, weights, is_epoched):
    if data.ndim != len(dim_names):
        raise ValueError(
            f"Expected data to have {len(dim_names)} dimensions, got {data.ndim}."
        )

    allowed_dims = ["epoch", "channel", "freq", "segment", "taper"]
    if not is_epoched:
        allowed_dims.remove("epoch")
    # TODO maybe we should be nice and allow plural versions of each dimname?
    for dim in dim_names:
        _check_option("dim_names", dim, allowed_dims)
    if "channel" not in dim_names or "freq" not in dim_names:
        raise ValueError("Both 'channel' and 'freq' must be present in `dim_names`.")

    if list(dim_names).index("channel") != int(is_epoched):
        raise ValueError(
            f"'channel' must be the {'second' if is_epoched else 'first'} dimension of "
            "the data."
        )
    want_n_chan = _pick_data_channels(info, exclude=()).size
    got_n_chan = data.shape[list(dim_names).index("channel")]
    if got_n_chan != want_n_chan:
        raise ValueError(
            f"The number of channels in `data` ({got_n_chan}) must match the number of "
            f"good + bad data channels in `info` ({want_n_chan})."
        )

    # given we limit max array size and ensure channel & freq dims present, only one of
    # taper or segment can be present
    if "taper" in dim_names:
        if dim_names[-2] != "taper":  # _psd_from_mt assumes this (called when plotting)
            raise ValueError(
                "'taper' must be the second to last dimension of the data."
            )
        # expect weights for each taper
        actual = None if weights is None else weights.size
        expected = data.shape[list(dim_names).index("taper")]
        if actual != expected:
            raise ValueError(
                f"Expected size of `weights` to be {expected} to match 'n_tapers' in "
                f"`data`, got {actual}."
            )
    elif "segment" in dim_names and dim_names[-1] != "segment":
        raise ValueError("'segment' must be the last dimension of the data.")

    # freq being in wrong position ruled out by above checks
    want_n_freq = freqs.size
    got_n_freq = data.shape[list(dim_names).index("freq")]
    if got_n_freq != want_n_freq:
        raise ValueError(
            f"The number of frequencies in `data` ({got_n_freq}) must match the number "
            f"of elements in `freqs` ({want_n_freq})."
        )


@fill_doc
class SpectrumArray(Spectrum):
    """Data object for precomputed spectral data (in NumPy array format).

    Parameters
    ----------
    data : ndarray, shape (n_channels, [n_tapers], n_freqs, [n_segments])
        The spectra for each channel.
    %(info_not_none)s
    %(freqs_tfr_array)s
    dim_names : tuple of str
        The name of the dimensions in the data, in the order they occur. Must contain
        ``'channel'`` and ``'freq'``;  if data are unaggregated estimates, also include
        either a ``'segment'`` (e.g., Welch-like algorithms) or ``'taper'`` (e.g.,
        multitaper algorithms) dimension. If including ``'taper'``, you should also pass
        a ``weights`` parameter.

        .. versionadded:: 1.8
    weights : ndarray | None
        Weights for the ``'taper'`` dimension, if present (see ``dim_names``).

        .. versionadded:: 1.8
    %(verbose)s

    See Also
    --------
    mne.create_info
    mne.EvokedArray
    mne.io.RawArray
    EpochsSpectrumArray

    Notes
    -----
    %(notes_spectrum_array)s

        .. versionadded:: 1.6
    """

    @verbose
    def __init__(
        self,
        data,
        info,
        freqs,
        dim_names=("channel", "freq"),
        weights=None,
        *,
        verbose=None,
    ):
        # (channel, [taper], freq, [segment])
        _check_option("data.ndim", data.ndim, (2, 3))  # only allow one extra dimension

        _check_data_shape(data, info, freqs, dim_names, weights, is_epoched=False)

        self.__setstate__(
            dict(
                method="unknown",
                data=data,
                sfreq=info["sfreq"],
                dims=dim_names,
                freqs=freqs,
                inst_type_str="Array",
                data_type=(
                    "Fourier Coefficients"
                    if np.iscomplexobj(data)
                    else "Power Spectrum"
                ),
                info=info,
                weights=weights,
            )
        )


@fill_doc
class EpochsSpectrum(BaseSpectrum, GetEpochsMixin):
    """Data object for spectral representations of epoched data.

    .. warning:: The preferred means of creating Spectrum objects from Epochs
                 is via the instance method :meth:`mne.Epochs.compute_psd`.
                 Direct class instantiation is not supported.

    Parameters
    ----------
    inst : instance of Epochs
        The data from which to compute the frequency spectrum.
    %(method_psd)s
    %(fmin_fmax_psd)s
    %(tmin_tmax_psd)s
    %(picks_good_data_noref)s
    %(exclude_psd)s
    %(proj_psd)s
    %(remove_dc)s
    %(n_jobs)s
    %(verbose)s
    %(method_kw_psd)s

    Attributes
    ----------
    ch_names : list
        The channel names.
    freqs : array
        Frequencies at which the amplitude, power, or fourier coefficients
        have been computed.
    %(info_not_none)s
    method : ``'welch'``| ``'multitaper'``
        The method used to compute the spectrum.
    weights : array | None
        The weights for each taper. Only present if spectra computed with
        ``method='multitaper'`` and ``output='complex'``.

        .. versionadded:: 1.8

    See Also
    --------
    EpochsSpectrumArray
    Spectrum
    mne.Epochs.compute_psd

    References
    ----------
    .. footbibliography::
    """

    def __init__(
        self,
        inst,
        method,
        fmin,
        fmax,
        tmin,
        tmax,
        picks,
        exclude,
        proj,
        remove_dc,
        *,
        n_jobs,
        verbose=None,
        **method_kw,
    ):
        # triage reading from file
        if isinstance(inst, dict):
            self.__setstate__(inst)
            return
        # do the basic setup
        super().__init__(
            inst,
            method,
            fmin,
            fmax,
            tmin,
            tmax,
            picks,
            exclude,
            proj,
            remove_dc,
            n_jobs=n_jobs,
            verbose=verbose,
            **method_kw,
        )
        # get just the data we want
        data = self.inst._get_data(picks=self._picks, on_empty="raise")[
            :, :, self._time_mask
        ]
        # compute the spectra
        self._compute_spectra(data, fmin, fmax, n_jobs, method_kw, verbose)
        self._dims = ("epoch",) + self._dims
        self._shape = (len(self.inst),) + self._shape
        # check for correct shape and bad values
        self._check_values()
        del self._shape
        # we need these for to_data_frame()
        self.event_id = self.inst.event_id.copy()
        self.events = self.inst.events.copy()
        self.selection = self.inst.selection.copy()
        # we need these for __getitem__()
        self.drop_log = deepcopy(self.inst.drop_log)
        self._metadata = self.inst.metadata
        # save memory
        del self.inst

    def __getitem__(self, item):
        """Subselect epochs from an EpochsSpectrum.

        Parameters
        ----------
        item : int | slice | array-like | str
            Access options are the same as for :class:`~mne.Epochs` objects,
            see the docstring of :meth:`mne.Epochs.__getitem__` for
            explanation.

        Returns
        -------
        %(getitem_epochspectrum_return)s
        """
        return super().__getitem__(item)

    def __getstate__(self):
        """Prepare object for serialization."""
        out = super().__getstate__()
        out.update(
            metadata=self._metadata,
            drop_log=self.drop_log,
            event_id=self.event_id,
            events=self.events,
            selection=self.selection,
        )
        return out

    def __setstate__(self, state):
        """Unpack from serialized format."""
        super().__setstate__(state)
        self._metadata = state["metadata"]
        self.drop_log = state["drop_log"]
        self.event_id = state["event_id"]
        self.events = state["events"]
        self.selection = state["selection"]

    def average(self, method="mean"):
        """Average the spectra across epochs.

        Parameters
        ----------
        method : 'mean' | 'median' | callable
            How to aggregate spectra across epochs. If callable, must take a
            :class:`NumPy array<numpy.ndarray>` of shape
            ``(n_epochs, n_channels, n_freqs)`` and return an array of shape
            ``(n_channels, n_freqs)``. Default is ``'mean'``.

        Returns
        -------
        spectrum : instance of Spectrum
            The aggregated spectrum object.
        """
        _validate_type(method, ("str", "callable"), "method")
        method = _make_combine_callable(
            method, axis=0, valid=("mean", "median"), keepdims=False
        )
        if not callable(method):
            raise ValueError(
                '"method" must be a valid string or callable, '
                f"got a {type(method).__name__} ({method})."
            )
        # averaging unaggregated spectral estimates are not supported
        if "segment" in self._dims:
            raise NotImplementedError(
                "Averaging individual Welch segments across epochs is not "
                "supported. Consider averaging the signals before computing "
                "the Welch spectrum estimates."
            )
        if "taper" in self._dims:
            raise NotImplementedError(
                "Averaging multitaper tapers across epochs is not supported. Consider "
                "averaging the signals before computing the complex spectrum."
            )
        # serialize the object and update data, dims, and data type
        state = super().__getstate__()
        state["nave"] = state["data"].shape[0]
        state["data"] = method(state["data"])
        state["dims"] = state["dims"][1:]
        state["data_type"] = f'Averaged {state["data_type"]}'
        defaults = dict(
            method=None,
            fmin=None,
            fmax=None,
            tmin=None,
            tmax=None,
            picks=None,
            exclude=(),
            proj=None,
            remove_dc=None,
            reject_by_annotation=None,
            n_jobs=None,
            verbose=None,
        )
        return Spectrum(state, **defaults)


@fill_doc
class EpochsSpectrumArray(EpochsSpectrum):
    """Data object for precomputed epoched spectral data (in NumPy array format).

    Parameters
    ----------
    data : ndarray, shape (n_epochs, n_channels, [n_tapers], n_freqs, [n_segments])
        The spectra for each channel in each epoch.
    %(info_not_none)s
    %(freqs_tfr_array)s
    %(events_epochs)s
    %(event_id)s
    dim_names : tuple of str
        The name of the dimensions in the data, in the order they occur. Must contain
        ``'channel'`` and ``'freq'``;  if data are unaggregated estimates, also include
        either a ``'segment'`` (e.g., Welch-like algorithms) or ``'taper'`` (e.g.,
        multitaper algorithms) dimension. If including ``'taper'``, you should also pass
        a ``weights`` parameter.

        .. versionadded:: 1.8
    weights : ndarray | None
        Weights for the ``'taper'`` dimension, if present (see ``dim_names``).

        .. versionadded:: 1.8
    %(verbose)s

    See Also
    --------
    mne.create_info
    mne.EpochsArray
    SpectrumArray

    Notes
    -----
    %(notes_spectrum_array)s

        .. versionadded:: 1.6
    """

    @verbose
    def __init__(
        self,
        data,
        info,
        freqs,
        events=None,
        event_id=None,
        dim_names=("epoch", "channel", "freq"),
        weights=None,
        *,
        verbose=None,
    ):
        # (epoch, channel, [taper], freq, [segment])
        _check_option("data.ndim", data.ndim, (3, 4))  # only allow one extra dimension

        if list(dim_names).index("epoch") != 0:
            raise ValueError("'epoch' must be the first dimension of `data`.")
        if events is not None and data.shape[0] != events.shape[0]:
            raise ValueError(
                f"The first dimension of `data` ({data.shape[0]}) must match the first "
                f"dimension of `events` ({events.shape[0]})."
            )

        _check_data_shape(data, info, freqs, dim_names, weights, is_epoched=True)

        self.__setstate__(
            dict(
                method="unknown",
                data=data,
                sfreq=info["sfreq"],
                dims=dim_names,
                freqs=freqs,
                inst_type_str="Array",
                data_type=(
                    "Fourier Coefficients"
                    if np.iscomplexobj(data)
                    else "Power Spectrum"
                ),
                info=info,
                events=events,
                event_id=event_id,
                metadata=None,
                selection=np.arange(data.shape[0]),
                drop_log=tuple(tuple() for _ in range(data.shape[0])),
                weights=weights,
            )
        )


def read_spectrum(fname):
    """Load a :class:`mne.time_frequency.Spectrum` object from disk.

    Parameters
    ----------
    fname : path-like
        Path to a spectrum file in HDF5 format, which should end with ``.h5`` or
        ``.hdf5``.

    Returns
    -------
    spectrum : instance of Spectrum
        The loaded Spectrum object.

    See Also
    --------
    mne.time_frequency.Spectrum.save
    """
    read_hdf5, _ = _import_h5io_funcs()
    _validate_type(fname, "path-like", "fname")
    fname = _check_fname(fname=fname, overwrite="read", must_exist=False)
    # read it in
    hdf5_dict = read_hdf5(fname, title="mnepython")
    defaults = dict(
        method=None,
        fmin=None,
        fmax=None,
        tmin=None,
        tmax=None,
        picks=None,
        exclude=(),
        proj=None,
        remove_dc=None,
        reject_by_annotation=None,
        n_jobs=None,
        verbose=None,
    )
    Klass = EpochsSpectrum if hdf5_dict["inst_type_str"] == "Epochs" else Spectrum
    return Klass(hdf5_dict, **defaults)


def _check_ci(ci):
    ci = "sd" if ci == "std" else ci  # be forgiving
    if _is_numeric(ci):
        if not (0 < ci <= 100):
            raise ValueError(f"ci must satisfy 0 < ci <= 100, got {ci}")
        ci /= 100.0
    else:
        _check_option("ci", ci, [None, "sd", "range"])
    return ci


def _compute_n_welch_segments(n_times, method_kw):
    # get default values from psd_array_welch
    _defaults = dict()
    for param in ("n_fft", "n_per_seg", "n_overlap"):
        _defaults[param] = signature(psd_array_welch).parameters[param].default
    # override defaults with user-specified values
    for key, val in _defaults.items():
        _defaults.update({key: method_kw.get(key, val)})
    # sanity check values / replace `None`s with real numbers
    n_fft, n_per_seg, n_overlap = _check_nfft(n_times, **_defaults)
    # compute expected number of segments
    step = n_per_seg - n_overlap
    return (n_times - n_overlap) // step


def _validate_method(method, instance_type):
    """Convert 'auto' to a real method name, and validate."""
    if method == "auto":
        method = "welch" if instance_type.startswith("Raw") else "multitaper"
    _check_option("method", method, ("welch", "multitaper"))
    return method