def fit_fooof_3d(fg, freqs, power_spectra, freq_range=None, n_jobs=1):
    """Fit FOOOF models across a 3d array of power spectra.

    Parameters
    ----------
    fg : FOOOFGroup
        Object to fit with, initialized with desired settings.
    freqs : 1d array
        Frequency values for the power spectra, in linear space.
    power_spectra : 3d array
        Power values, in linear space, with shape as: [n_conditions, n_power_spectra, n_freqs].
    freq_range : list of [float, float], optional
        Desired frequency range to fit. If not provided, fits the entire given range.
    n_jobs : int, optional, default: 1
        Number of jobs to run in parallel.
        1 is no parallelization. -1 uses all available cores.

    Returns
    -------
    fgs : list of FOOOFGroups
        Collected FOOOFGroups after fitting across power spectra, length of n_conditions.


    Examples
    --------
    Fit a 3d array of power spectra, assuming `freqs` and `spectra` are already defined:

    >>> from fooof import FOOOFGroup
    >>> fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_height=0.1)
    >>> fgs = fit_fooof_3d(fg, freqs, power_spectra, freq_range=[3, 30])  # doctest:+SKIP
    """
    fgs = []
    for cond_spectra in power_spectra:
        fg.fit(freqs, cond_spectra, freq_range, n_jobs)
        fgs.append(fg.copy())
    return fgs
