diff --git a/doc/changes/dev/13521.bugfix.rst b/doc/changes/dev/13521.bugfix.rst
new file mode 100644
index 00000000000..bba18190331
--- /dev/null
+++ b/doc/changes/dev/13521.bugfix.rst
@@ -0,0 +1 @@
+Fix bug preventing :func:`mne.time_frequency.read_spectrum` from reading saved :class:`mne.time_frequency.Spectrum` objects created from :meth:`mne.time_frequency.EpochsSpectrum.average`, by `Thomas Binns`_.
\ No newline at end of file
diff --git a/mne/io/edf/edf.py b/mne/io/edf/edf.py
index b03a80d5763..b1ebe3ce622 100644
--- a/mne/io/edf/edf.py
+++ b/mne/io/edf/edf.py
@@ -2117,8 +2117,11 @@ def read_raw_bdf(
         >>> events[:, 2] &= (2**16 - 1)  # doctest:+SKIP
 
     The above operation can be carried out directly in :func:`mne.find_events`
-    using the ``mask`` and ``mask_type`` parameters (see
-    :func:`mne.find_events` for more details).
+    using the ``mask`` parameter as follows:
+
+        >>> events = mne.find_events(..., mask=2**16 - 1)  # doctest:+SKIP
+
+    See :func:`mne.find_events` for more details.
 
     It is also possible to retrieve system codes, but no particular effort has
     been made to decode these in MNE. In case it is necessary, for instance to
diff --git a/mne/time_frequency/spectrum.py b/mne/time_frequency/spectrum.py
index 0caa7c0ae20..0d0ce0c30c8 100644
--- a/mne/time_frequency/spectrum.py
+++ b/mne/time_frequency/spectrum.py
@@ -1752,7 +1752,7 @@ def read_spectrum(fname):
         n_jobs=None,
         verbose=None,
     )
-    Klass = EpochsSpectrum if hdf5_dict["inst_type_str"] == "Epochs" else Spectrum
+    Klass = EpochsSpectrum if "epoch" in hdf5_dict["dims"] else Spectrum
     return Klass(hdf5_dict, **defaults)
 
 
diff --git a/mne/time_frequency/tests/test_spectrum.py b/mne/time_frequency/tests/test_spectrum.py
index ddb36404f4e..b1ad677352d 100644
--- a/mne/time_frequency/tests/test_spectrum.py
+++ b/mne/time_frequency/tests/test_spectrum.py
@@ -190,6 +190,12 @@ def test_spectrum_io(inst, tmp_path, request, evoked):
     orig.save(fname)
     loaded = read_spectrum(fname)
     assert orig == loaded
+    # Test Spectrum from EpochsSpectrum.average() can be read (gh-13521)
+    if isinstance(inst, BaseEpochs):
+        origavg = orig.average()
+        origavg.save(fname, overwrite=True)
+        loadedavg = read_spectrum(fname)
+        assert origavg == loadedavg
 
 
 def test_spectrum_copy(raw_spectrum):
@@ -320,6 +326,7 @@ def test_epochs_spectrum_average(epochs_spectrum, method):
     avg_spect = epochs_spectrum.average(method=method)
     assert avg_spect.shape == epochs_spectrum.shape[1:]
     assert avg_spect._dims == ("channel", "freq")  # no 'epoch'
+    assert repr(avg_spect).startswith("<Averaged Power Spectrum (from Epochs")
 
 
 @pytest.mark.parametrize("inst", ("raw_spectrum", "epochs_spectrum", "evoked"))
diff --git a/mne/time_frequency/tests/test_tfr.py b/mne/time_frequency/tests/test_tfr.py
index 72c8a68ae4c..ed6ddd6da82 100644
--- a/mne/time_frequency/tests/test_tfr.py
+++ b/mne/time_frequency/tests/test_tfr.py
@@ -674,6 +674,11 @@ def test_tfr_io(inst, average_tfr, request, tmp_path):
     with tfr.info._unlock():
         tfr.info["meas_date"] = want
     assert tfr_loaded == tfr
+    # test AverageTFR from EpochsTFR.average() can be read (gh-13521)
+    tfravg = tfr.average()
+    tfravg.save(fname, overwrite=True)
+    tfravg_loaded = read_tfrs(fname)
+    assert tfravg == tfravg_loaded
     # test with taper dimension and weights
     n_tapers = 3  # anything >= 1 should do
     weights = np.ones((n_tapers, tfr.shape[2]))  # tapers x freqs
@@ -1217,6 +1222,10 @@ def test_averaging_epochsTFR():
     ):
         tapered.average()
 
+    # Test repr from original instance info is preserved
+    avgpower = power.average()
+    assert repr(avgpower).startswith("<Average Power from Epochs")
+
 
 def test_averaging_freqsandtimes_epochsTFR():
     """Test that EpochsTFR averaging freqs methods work."""
diff --git a/mne/time_frequency/tfr.py b/mne/time_frequency/tfr.py
index 3a69c93bc8d..f64680845c4 100644
--- a/mne/time_frequency/tfr.py
+++ b/mne/time_frequency/tfr.py
@@ -3313,7 +3313,6 @@ def average(self, method="mean", *, dim="epochs", copy=False):
         state["freqs"] = freqs
         state["times"] = times
         if dim == "epochs":
-            state["inst_type_str"] = "Evoked"
             state["nave"] = n_epochs
             state["comment"] = f"{method} of {n_epochs} EpochsTFR{_pl(n_epochs)}"
             out = AverageTFR(inst=state)