Repository: StingraySoftware/stingray

Issue:
Roundtrip to Astropy tables, also for AveragedXYSpectrum and VarEnergySpectrum

Discussion:
pep8speaks: Hello @matteobachetti! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`stingray/base.py`](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/base.py):

> [Line 71:47](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/base.py#L71): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 309:30](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/base.py#L309): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ':'
> [Line 309:32](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/base.py#L309): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 310:5](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/base.py#L310): [E129](https://duckduckgo.com/?q=pep8%20E129) visually indented line with same indent as next logical line
> [Line 310:30](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/base.py#L310): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ':'
> [Line 310:32](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/base.py#L310): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals



* In the file [`stingray/tests/test_base.py`](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/tests/test_base.py):

> [Line 28:5](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/tests/test_base.py#L28): [E301](https://duckduckgo.com/?q=pep8%20E301) expected 1 blank line, found 0

* In the file [`stingray/tests/test_crossspectrum.py`](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/tests/test_crossspectrum.py):

> [Line 1327:33](https://github.com/StingraySoftware/stingray/blob/8b1d2a79045ab27fd3b3b64c63e15db3b0683590/stingray/tests/test_crossspectrum.py#L1327): [E225](https://duckduckgo.com/?q=pep8%20E225) missing whitespace around operator








codecov[bot]: # [Codecov](https://codecov.io/gh/StingraySoftware/stingray/pull/630?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) Report
> Merging [#630](https://codecov.io/gh/StingraySoftware/stingray/pull/630?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) (d3a8268) into [main](https://codecov.io/gh/StingraySoftware/stingray/commit/16ae8dc9c5630a17ba7dd2a0ecbf189db38be4ae?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) (16ae8dc) will **decrease** coverage by `6.67%`.
> The diff coverage is `82.97%`.

[![Impacted file tree graph](https://codecov.io/gh/StingraySoftware/stingray/pull/630/graphs/tree.svg?width=650&height=150&src=pr&token=FjWeFfhU9F&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware)](https://codecov.io/gh/StingraySoftware/stingray/pull/630?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware)

```diff
@@            Coverage Diff             @@
##             main     #630      +/-   ##
==========================================
- Coverage   97.02%   90.34%   -6.68%     
==========================================
  Files          40       41       +1     
  Lines        7664     7533     -131     
==========================================
- Hits         7436     6806     -630     
- Misses        228      727     +499     
```


| [Impacted Files](https://codecov.io/gh/StingraySoftware/stingray/pull/630?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) | Coverage Δ | |
|---|---|---|
| [stingray/io.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvaW8ucHk=) | `98.78% <ø> (+3.06%)` | :arrow_up: |
| [stingray/base.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvYmFzZS5weQ==) | `78.44% <78.44%> (ø)` | |
| [stingray/simulator/simulator.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvc2ltdWxhdG9yL3NpbXVsYXRvci5weQ==) | `99.43% <93.33%> (-0.57%)` | :arrow_down: |
| [stingray/crossspectrum.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvY3Jvc3NzcGVjdHJ1bS5weQ==) | `96.00% <100.00%> (-2.38%)` | :arrow_down: |
| [stingray/events.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvZXZlbnRzLnB5) | `100.00% <100.00%> (ø)` | |
| [stingray/lightcurve.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvbGlnaHRjdXJ2ZS5weQ==) | `96.92% <100.00%> (-0.93%)` | :arrow_down: |
| [stingray/powerspectrum.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvcG93ZXJzcGVjdHJ1bS5weQ==) | `95.49% <100.00%> (-3.90%)` | :arrow_down: |
| [stingray/simulator/transfer.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvc2ltdWxhdG9yL3RyYW5zZmVyLnB5) | `98.94% <100.00%> (+0.05%)` | :arrow_up: |
| [stingray/utils.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvdXRpbHMucHk=) | `95.71% <100.00%> (-3.32%)` | :arrow_down: |
| [stingray/varenergyspectrum.py](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvdmFyZW5lcmd5c3BlY3RydW0ucHk=) | `99.61% <100.00%> (+<0.01%)` | :arrow_up: |
| ... and [16 more](https://codecov.io/gh/StingraySoftware/stingray/pull/630/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/StingraySoftware/stingray/pull/630?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/StingraySoftware/stingray/pull/630?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware). Last update [16ae8dc...d3a8268](https://codecov.io/gh/StingraySoftware/stingray/pull/630?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware).


Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 3d79e075..a5c0a5f4 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -4,17 +4,23 @@ Changelog
 Not Yet Released
 ----------------
 
+- A new infrastructure for converting ``EventList`` and ``LightCurve`` objects into Astropy ``TimeSeries``
+- A new infrastructure for converting most ``Stingray`` classes into Astropy ``Table`` objects
+- Save and load of most ``Stingray`` classes to/from many different file formats (``pickle``, ``ECSV``, ``HDF5``, ``FITS``, and all formats compatible with Astropy Table)
+- Accept input ``EventList`` in ``DynamicalPowerSpectrum``
+- Implement Gardner-Done normalization (1 for perfect correlation, -1 for perfect anticorrelation) for ``Auto/Crosscorrelation``
+- Make modulation upper limit in ``(Averaged)Powerspectrum`` work with any normalization
 - Dropped support to Python <3.8
 - Multi-taper periodogram, including a Lomb-Scargle implementation for non-uniformly sampled data.
 - Upper limits on pulsations in periodograms and Z searches
-- New stingray.fourier module containing the basic timing products, usable on `numpy` arrays
-- Lots of performance improvements in the `AveragedCrossspectrum` and `AveragedPowerspectrum` classes
-- New methods in `Crossspectrum` and `Powerspectrum` to load data from specific inputs: `from_events`, `from_lightcurve`, `from_time_array`, `from_lc_list`
-    - `from_time_array` was also tested using memory-mapped event lists as inputs: useful in very large datasets
-- New and improved spectral timing methods: `ComplexCovarianceSpectrum`, `CovarianceSpectrum`, `LagSpectrum`, `RmsSpectrum`
+- New stingray.fourier module containing the basic timing products, usable on ``numpy`` arrays
+- Lots of performance improvements in the ``AveragedCrossspectrum`` and ``AveragedPowerspectrum`` classes
+- New methods in ``Crossspectrum`` and ``Powerspectrum`` to load data from specific inputs: ``from_events``, ``from_lightcurve``, ``from_time_array``, ``from_lc_list``
+    - ``from_time_array`` was also tested using memory-mapped event lists as inputs: useful in very large datasets
+- New and improved spectral timing methods: ``ComplexCovarianceSpectrum``, ``CovarianceSpectrum``, ``LagSpectrum``, ``RmsSpectrum``
 - Improved error bars on cross-spectral and spectral timing methods
 - Some deprecated features are now removed
-- PSDLogLikelihood now also works with a log-rebinned PDS
+- ``PSDLogLikelihood`` now also works with a log-rebinned PDS
 - PDS/CrossSp initially store the unnormalized power, and convert it on the fly when requested, to any normalization
 - Lots of bug fixes
 
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 816e7544..02ba12df 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -83,7 +83,9 @@ Code Reviews are super-useful: another contributor can review the code, which me
 - **Imports:**
     - Absolute imports are to be used in general. The exception to this is relative imports of the form `from . import modulename`, this convention makes it clearer what code is from the current sub-module as opposed to from another. It is best to use when referring to files within the same sub-module.
     - The import `numpy as np`, `import scipy as sp`, `import matplotlib as mpl`, and `import matplotlib.pyplot as plt` naming conventions should be used wherever relevant. `from packagename import *` should never be used, except as a tool to flatten the namespace of a module.
-- **Variable access in Classes:** Classes should either use direct variable access, or Python’s property mechanism for setting object instance variables. `get_value/set_value` style methods should be used only when getting and setting the values requires a computationally-expensive operation.
+- **Variable access in Classes:**
+    - Classes should either use direct variable access, or Python’s property mechanism for setting object instance variables. `get_value/set_value` style methods should be used only when getting and setting the values requires a computationally-expensive operation.
+    - Attribute names should be descriptive if possible, use names of desserts otherwise (e.g. for dummy test classes)
 - **super() function:** Classes should use the built-in `super()` function when making calls to methods in their super-class(es) unless there are specific reasons not to. `super()` should be used consistently in all sub-classes since it does not work otherwise.
 - **Multiple Inheritance:** Multiple inheritance should be avoided in general without good reason.
 - **__init__.py:** The `__init__.py` files for modules should not contain any significant implementation code.  `__init__.py` can contain docstrings and code for organizing the module layout, however if a module is small enough that it fits in one file, it should simply be a single file, rather than a directory with an  `__init__.py` file.
diff --git a/setup.cfg b/setup.cfg
index f7b4273c..25496e31 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -83,6 +83,7 @@ filterwarnings =
     ignore:table path was not set via the path= argument; using default path:UserWarning
     ignore:Large Datasets may not be processed efficiently:UserWarning
     ignore:.*is a deprecated alias for:DeprecationWarning
+    ignore:.*HIERARCH card will be created.*:
 
 ;addopts = --disable-warnings
 
diff --git a/stingray/base.py b/stingray/base.py
new file mode 100644
index 00000000..7bf8d934
--- /dev/null
+++ b/stingray/base.py
@@ -0,0 +1,628 @@
+"""Base classes"""
+from collections.abc import Iterable
+import pickle
+import warnings
+import copy
+
+import numpy as np
+from astropy.table import Table
+from astropy.time import Time, TimeDelta
+from astropy.units import Quantity
+
+
+class StingrayObject(object):
+    """This base class defines some general-purpose utilities.
+
+    The main purpose is to have a consistent mechanism for:
+
+    + round-tripping to and from Astropy Tables and other dataframes
+
+    + round-tripping to files in different formats
+
+    The idea is that any object inheriting :class:`StingrayObject` should,
+    just by defining an attribute called ``main_array_attr``, be able to perform
+    the operations above, with no additional effort.
+
+    ``main_array_attr`` is, e.g. ``time`` for :class:`EventList` and
+    :class:`Lightcurve`, ``freq`` for :class:`Crossspectrum`, ``energy`` for
+    :class:`VarEnergySpectrum`, and so on. It is the array with wich all other
+    attributes are compared: if they are of the same shape, they get saved as
+    columns of the table/dataframe, otherwise as metadata.
+    """
+
+    def __init__(cls, *args, **kwargs):
+        if not hasattr(cls, "main_array_attr"):
+            raise RuntimeError(
+                "A StingrayObject needs to have the main_array_attr attribute specified"
+            )
+
+    def array_attrs(self):
+        """List the names of the array attributes of the Stingray Object.
+
+        By array attributes, we mean the ones with the same size and shape as
+        ``main_array_attr`` (e.g. ``time`` in ``EventList``)
+        """
+
+        main_attr = getattr(self, getattr(self, "main_array_attr"))
+        if main_attr is None:
+            return []
+
+        return [
+            attr
+            for attr in dir(self)
+            if (
+                isinstance(getattr(self, attr), Iterable)
+                and np.shape(getattr(self, attr)) == np.shape(main_attr)
+            )
+        ]
+
+    def meta_attrs(self):
+        """List the names of the meta attributes of the Stingray Object.
+
+        By array attributes, we mean the ones with a different size and shape
+        than ``main_array_attr`` (e.g. ``time`` in ``EventList``)
+        """
+        array_attrs = self.array_attrs()
+        return [
+            attr
+            for attr in dir(self)
+            if (
+                attr not in array_attrs
+                and not attr.startswith("_")
+                # Use new assignment expression (PEP 572). I'm testing that
+                # self.attribute is not callable, and assigning its value to
+                # the variable attr_value for further checks
+                and not callable(attr_value := getattr(self, attr))
+                # a way to avoid EventLists, Lightcurves, etc.
+                and not hasattr(attr_value, "meta_attrs")
+            )
+        ]
+
+    def get_meta_dict(self):
+        """Give a dictionary with all non-None meta attrs of the object."""
+        meta_attrs = self.meta_attrs()
+        meta_dict = {}
+        for key in meta_attrs:
+            val = getattr(self, key)
+            if val is not None:
+                meta_dict[key] = val
+        return meta_dict
+
+    def to_astropy_table(self):
+        """Create an Astropy Table from a ``StingrayObject``
+
+        Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
+        ``EventList``) are converted into columns, while meta attributes
+        (``mjdref``, ``gti``, etc.) are saved into the ``meta`` dictionary.
+        """
+        data = {}
+        array_attrs = self.array_attrs()
+
+        for attr in array_attrs:
+            data[attr] = np.asarray(getattr(self, attr))
+
+        ts = Table(data)
+
+        ts.meta.update(self.get_meta_dict())
+
+        return ts
+
+    @classmethod
+    def from_astropy_table(cls, ts):
+        """Create a Stingray Object object from data in an Astropy Table.
+
+        The table MUST contain at least a column named like the
+        ``main_array_attr``.
+        The rest of columns will form the array attributes of the
+        new object, while the attributes in ds.attrs will
+        form the new meta attributes of the object.
+
+        It is strongly advisable to define such attributes and columns
+        using the standard attributes of the wanted StingrayObject (e.g.
+        ``time``, ``pi``, etc. for ``EventList``)
+        """
+        cls = cls()
+
+        if len(ts) == 0:
+            # return an empty object
+            return cls
+
+        array_attrs = ts.colnames
+
+        # Set the main attribute first
+        mainarray = np.array(ts[cls.main_array_attr])
+        setattr(cls, cls.main_array_attr, mainarray)
+
+        for attr in array_attrs:
+            if attr == cls.main_array_attr:
+                continue
+            setattr(cls, attr.lower(), np.array(ts[attr]))
+
+        for key, val in ts.meta.items():
+            setattr(cls, key.lower(), val)
+
+        return cls
+
+    def to_xarray(self):
+        """Create an ``xarray`` Dataset from a `StingrayObject`.
+
+        Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
+        ``EventList``) are converted into columns, while meta attributes
+        (``mjdref``, ``gti``, etc.) are saved into the ``ds.attrs`` dictionary.
+        """
+        from xarray import Dataset
+
+        data = {}
+        array_attrs = self.array_attrs()
+
+        for attr in array_attrs:
+            data[attr] = np.asarray(getattr(self, attr))
+
+        ts = Dataset(data)
+
+        ts.attrs.update(self.get_meta_dict())
+
+        return ts
+
+    @classmethod
+    def from_xarray(cls, ts):
+        """Create a `StingrayObject` from data in an xarray Dataset.
+
+        The dataset MUST contain at least a column named like the
+        ``main_array_attr``.
+        The rest of columns will form the array attributes of the
+        new object, while the attributes in ds.attrs will
+        form the new meta attributes of the object.
+
+        It is strongly advisable to define such attributes and columns
+        using the standard attributes of the wanted StingrayObject (e.g.
+        ``time``, ``pi``, etc. for ``EventList``)
+        """
+        cls = cls()
+
+        if len(ts[cls.main_array_attr]) == 0:
+            # return an empty object
+            return cls
+
+        array_attrs = ts.coords
+
+        # Set the main attribute first
+        mainarray = np.array(ts[cls.main_array_attr])
+        setattr(cls, cls.main_array_attr, mainarray)
+
+        for attr in array_attrs:
+            if attr == cls.main_array_attr:
+                continue
+            setattr(cls, attr, np.array(ts[attr]))
+
+        for key, val in ts.attrs.items():
+            if key not in array_attrs:
+                setattr(cls, key, val)
+
+        return cls
+
+    def to_pandas(self):
+        """Create a pandas ``DataFrame`` from a :class:`StingrayObject`.
+
+        Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
+        ``EventList``) are converted into columns, while meta attributes
+        (``mjdref``, ``gti``, etc.) are saved into the ``ds.attrs`` dictionary.
+        """
+        from pandas import DataFrame
+
+        data = {}
+        array_attrs = self.array_attrs()
+
+        for attr in array_attrs:
+            data[attr] = np.asarray(getattr(self, attr))
+
+        ts = DataFrame(data)
+
+        ts.attrs.update(self.get_meta_dict())
+
+        return ts
+
+    @classmethod
+    def from_pandas(cls, ts):
+        """Create an `StingrayObject` object from data in a pandas DataFrame.
+
+        The dataframe MUST contain at least a column named like the
+        ``main_array_attr``.
+        The rest of columns will form the array attributes of the
+        new object, while the attributes in ds.attrs will
+        form the new meta attributes of the object.
+
+        It is strongly advisable to define such attributes and columns
+        using the standard attributes of the wanted StingrayObject (e.g.
+        ``time``, ``pi``, etc. for ``EventList``)
+
+        """
+        cls = cls()
+
+        if len(ts) == 0:
+            # return an empty object
+            return cls
+
+        array_attrs = ts.columns
+
+        # Set the main attribute first
+        mainarray = np.array(ts[cls.main_array_attr])
+        setattr(cls, cls.main_array_attr, mainarray)
+
+        for attr in array_attrs:
+            if attr == cls.main_array_attr:
+                continue
+            setattr(cls, attr, np.array(ts[attr]))
+
+        for key, val in ts.attrs.items():
+            if key not in array_attrs:
+                setattr(cls, key, val)
+
+        return cls
+
+    @classmethod
+    def read(cls, filename, fmt=None, format_=None):
+        r"""Generic reader for :class`StingrayObject`
+
+        Currently supported formats are
+
+        * pickle (not recommended for long-term storage)
+        * any other formats compatible with the writers in
+          :class:`astropy.table.Table` (ascii.ecsv, hdf5, etc.)
+
+        Files that need the :class:`astropy.table.Table` interface MUST contain
+        at least a column named like the ``main_array_attr``.
+        The default ascii format is enhanced CSV (ECSV). Data formats
+        supporting the serialization of metadata (such as ECSV and HDF5) can
+        contain all attributes such as ``mission``, ``gti``, etc with
+        no significant loss of information. Other file formats might lose part
+        of the metadata, so must be used with care.
+
+        ..note::
+
+            Complex values can be dealt with out-of-the-box in some formats
+            like HDF5 or FITS, not in others (e.g. all ASCII formats).
+            With these formats, and in any case when fmt is ``None``, complex
+            values should be stored as two columns of real numbers, whose names
+            are of the format <variablename>.real and <variablename>.imag
+
+        Parameters
+        ----------
+        filename: str
+            Path and file name for the file to be read.
+
+        fmt: str
+            Available options are 'pickle', 'hea', and any `Table`-supported
+            format such as 'hdf5', 'ascii.ecsv', etc.
+
+        Returns
+        -------
+        obj: :class:`StingrayObject` object
+            The object reconstructed from file
+        """
+        if fmt is None and format_ is not None:
+            warnings.warn(
+                "The format_ keyword for read and write is deprecated. Use fmt instead",
+                DeprecationWarning,
+            )
+            fmt = format_
+
+        if fmt is None:
+            pass
+        elif fmt.lower() == "pickle":
+            with open(filename, "rb") as fobj:
+                return pickle.load(fobj)
+        elif fmt.lower() == "ascii":
+            fmt = "ascii.ecsv"
+
+        ts = Table.read(filename, format=fmt)
+
+        # For specific formats, and in any case when the format is not
+        # specified, make sure that complex values are treated correctly.
+        if fmt is None or "ascii" in fmt:
+            for col in ts.colnames:
+                if not (
+                    (is_real := col.endswith(".real"))
+                    or (is_imag := col.endswith(".imag"))
+                ):
+                    continue
+
+                new_value = ts[col]
+
+                if is_imag:
+                    new_value = new_value * 1j
+
+                # Make sure it's complex, even if we find the real part first
+                new_value = new_value + 0.0j
+
+                col_strip = col.replace(".real", "").replace(".imag", "")
+
+                if col_strip not in ts.colnames:
+                    # If the column without ".real" or ".imag" doesn't exist,
+                    # define it, and make sure it's complex-valued
+                    ts[col_strip] = new_value
+                else:
+                    # If it does exist, sum the new value to it.
+                    ts[col_strip] += new_value
+
+                ts.remove_column(col)
+
+        return cls.from_astropy_table(ts)
+
+    def write(self, filename, fmt=None, format_=None):
+        """Generic writer for :class`StingrayObject`
+
+        Currently supported formats are
+
+        * pickle (not recommended for long-term storage)
+        * any other formats compatible with the writers in
+          :class:`astropy.table.Table` (ascii.ecsv, hdf5, etc.)
+
+        ..note::
+
+            Complex values can be dealt with out-of-the-box in some formats
+            like HDF5 or FITS, not in others (e.g. all ASCII formats).
+            With these formats, and in any case when fmt is ``None``, complex
+            values will be stored as two columns of real numbers, whose names
+            are of the format <variablename>.real and <variablename>.imag
+
+        Parameters
+        ----------
+        filename: str
+            Name and path of the file to save the object list to.
+
+        fmt: str
+            The file format to store the data in.
+            Available options are ``pickle``, ``hdf5``, ``ascii``, ``fits``
+        """
+        if fmt is None and format_ is not None:
+            warnings.warn(
+                "The format_ keyword for read and write is deprecated. Use fmt instead",
+                DeprecationWarning,
+            )
+            fmt = format_
+        if fmt is None:
+            pass
+        elif fmt.lower() == "pickle":
+            with open(filename, "wb") as fobj:
+                pickle.dump(self, fobj)
+            return
+        elif fmt.lower() == "ascii":
+            fmt = "ascii.ecsv"
+
+        ts = self.to_astropy_table()
+        if fmt is None or "ascii" in fmt:
+            for col in ts.colnames:
+                if np.iscomplex(ts[col][0]):
+                    ts[f"{col}.real"] = ts[col].real
+                    ts[f"{col}.imag"] = ts[col].imag
+                    ts.remove_column(col)
+
+        try:
+            ts.write(filename, format=fmt, overwrite=True, serialize_meta=True)
+        except TypeError:
+            ts.write(filename, format=fmt, overwrite=True)
+
+
+class StingrayTimeseries(StingrayObject):
+    def to_astropy_timeseries(self):
+        """Save the ``StingrayTimeseries`` to an ``Astropy`` timeseries.
+
+        Array attributes (time, pi, energy, etc.) are converted
+        into columns, while meta attributes (mjdref, gti, etc.)
+        are saved into the ``meta`` dictionary.
+
+        Returns
+        -------
+        ts : `astropy.timeseries.TimeSeries`
+            A ``TimeSeries`` object with the array attributes as columns,
+            and the meta attributes in the `meta` dictionary
+        """
+        from astropy.timeseries import TimeSeries
+        from astropy.time import TimeDelta
+        from astropy import units as u
+
+        data = {}
+        array_attrs = self.array_attrs()
+
+        for attr in array_attrs:
+            if attr == "time":
+                continue
+            data[attr] = np.asarray(getattr(self, attr))
+
+        if data == {}:
+            data = None
+
+        if self.time is not None and np.size(self.time) > 0:
+            times = TimeDelta(self.time * u.s)
+            ts = TimeSeries(data=data, time=times)
+        else:
+            ts = TimeSeries()
+
+        ts.meta.update(self.get_meta_dict())
+
+        return ts
+
+    @classmethod
+    def from_astropy_timeseries(cls, ts):
+        """Create a `StingrayTimeseries` from data in an Astropy TimeSeries
+
+        The timeseries has to define at least a column called time,
+        the rest of columns will form the array attributes of the
+        new event list, while the attributes in table.meta will
+        form the new meta attributes of the event list.
+
+        It is strongly advisable to define such attributes and columns
+        using the standard attributes of EventList: time, pi, energy, gti etc.
+
+        Parameters
+        ----------
+        ts : `astropy.timeseries.TimeSeries`
+            A ``TimeSeries`` object with the array attributes as columns,
+            and the meta attributes in the `meta` dictionary
+
+        Returns
+        -------
+        ts : `StingrayTimeseries`
+            Timeseries object
+        """
+
+        time = ts["time"]
+        mjdref = None
+        if "mjdref" in ts.meta:
+            mjdref = ts.meta["mjdref"]
+
+        time, mjdref = interpret_times(time, mjdref)
+        cls.time = np.asarray(time)
+
+        array_attrs = ts.colnames
+        for key, val in ts.meta.items():
+            setattr(cls, key, val)
+
+        for attr in array_attrs:
+            if attr == "time":
+                continue
+            setattr(cls, attr, ts[attr])
+
+        return cls
+
+    def change_mjdref(self, new_mjdref):
+        """Change the MJD reference time (MJDREF) of the time series
+
+        The times of the time series will be shifted in order to be referred to
+        this new MJDREF
+
+        Parameters
+        ----------
+        new_mjdref : float
+            New MJDREF
+
+        Returns
+        -------
+        new_lc : :class:`StingrayTimeseries` object
+            The new time series, shifted by MJDREF
+        """
+        time_shift = (self.mjdref - new_mjdref) * 86400
+
+        ts = self.shift(time_shift)
+        ts.mjdref = new_mjdref
+        return ts
+
+    def shift(self, time_shift):
+        """Shift the time and the GTIs by the same amount
+
+        Parameters
+        ----------
+        time_shift: float
+            The time interval by which the light curve will be shifted (in
+            the same units as the time array in :class:`Lightcurve`
+
+        Returns
+        -------
+        ts : ``StingrayTimeseries`` object
+            The new time series shifted by ``time_shift``
+
+        """
+        ts = copy.deepcopy(self)
+        ts.time = np.asarray(ts.time) + time_shift
+        if hasattr(ts, "gti"):
+            ts.gti = np.asarray(ts.gti) + time_shift
+
+        return ts
+
+
+def interpret_times(time, mjdref=0):
+    """Understand the format of input times, and return seconds from MJDREF
+
+    Parameters
+    ----------
+    time : class:`astropy.Time`, class:`time.Time`, class:`astropy.TimeDelta`, class:`astropy.Quantity`, class:`np.array`
+        Input times.
+
+    Other Parameters
+    ----------------
+    mjdref : float
+        Input MJD reference of the times. Optional.
+
+    Returns
+    -------
+    time_s : class:`np.array`
+        Times, in seconds from MJDREF
+    mjdref : float
+        MJDREF. If the input time is a `time.Time` object and the input mjdref
+        is 0, it will be defined as the MJD of the input time.
+
+    Examples
+    --------
+    >>> import astropy.units as u
+    >>> time = Time(57483, format='mjd')
+    >>> newt, mjdref = interpret_times(time)
+    >>> newt == 0
+    True
+    >>> mjdref == 57483
+    True
+    >>> time = Time([57483], format='mjd')
+    >>> newt, mjdref = interpret_times(time)
+    >>> np.allclose(newt, 0)
+    True
+    >>> mjdref == 57483
+    True
+    >>> time = TimeDelta([3, 4, 5] * u.s)
+    >>> newt, mjdref = interpret_times(time)
+    >>> np.allclose(newt, [3, 4, 5])
+    True
+    >>> time = np.array([3, 4, 5])
+    >>> newt, mjdref = interpret_times(time, mjdref=45000)
+    >>> np.allclose(newt, [3, 4, 5])
+    True
+    >>> mjdref == 45000
+    True
+    >>> time = np.array([3, 4, 5] * u.s)
+    >>> newt, mjdref = interpret_times(time, mjdref=45000)
+    >>> np.allclose(newt, [3, 4, 5])
+    True
+    >>> mjdref == 45000
+    True
+    >>> newt, mjdref = interpret_times(1, mjdref=45000)
+    >>> newt == 1
+    True
+    >>> newt, mjdref = interpret_times(list, mjdref=45000)
+    Traceback (most recent call last):
+    ...
+    ValueError: Unknown time format: ...
+    >>> newt, mjdref = interpret_times("guadfkljfd", mjdref=45000)
+    Traceback (most recent call last):
+    ...
+    ValueError: Unknown time format: ...
+    """
+    if isinstance(time, TimeDelta):
+        out_times = time.to("s").value
+        return out_times, mjdref
+
+    if isinstance(time, Time):
+        mjds = time.mjd
+        if mjdref == 0:
+            if np.all(mjds > 10000):
+                if isinstance(mjds, Iterable):
+                    mjdref = mjds[0]
+                else:
+                    mjdref = mjds
+
+        out_times = (mjds - mjdref) * 86400
+        return out_times, mjdref
+
+    if isinstance(time, Quantity):
+        out_times = time.to("s").value
+        return out_times, mjdref
+
+    if isinstance(time, (tuple, list, np.ndarray)):
+        return time, mjdref
+
+    if not isinstance(time, Iterable):
+        try:
+            float(time)
+            return time, mjdref
+        except (ValueError, TypeError):
+            pass
+
+    raise ValueError(f"Unknown time format: {type(time)}")
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 26b5c9f1..9ed00d7e 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -14,6 +14,7 @@ from stingray.gti import bin_intervals_from_gtis, check_gtis, cross_two_gtis
 from stingray.largememory import createChunkedSpectra, saveData, HAS_ZARR
 from stingray.utils import genDataPath, rebin_data, rebin_data_log, simon
 
+from .base import StingrayObject
 from .events import EventList
 from .lightcurve import Lightcurve
 from .utils import show_progress
@@ -470,7 +471,10 @@ def cospectra_pvalue(power, nspec):
     return pval
 
 
-class Crossspectrum(object):
+class Crossspectrum(StingrayObject):
+    main_array_attr = "freq"
+    type = "crossspectrum"
+
     """
     Make a cross spectrum from a (binned) light curve.
     You can also make an empty :class:`Crossspectrum` object to populate with your
@@ -761,20 +765,6 @@ class Crossspectrum(object):
 
         return True
 
-    @property
-    def type(self):
-        if self._type is not None:
-            pass
-        elif self.__class__.__name__ in ["Powerspectrum", "AveragedPowerspectrum", "DynamicalPowerspectrum"]:
-            self._type = "powerspectrum"
-        elif self.__class__.__name__ in ["Crossspectrum", "AveragedCrossspectrum"]:
-            self._type = "crossspectrum"
-        return self._type
-
-    @type.setter
-    def type(self, value):
-        self._type = value
-
     def _make_auxil_pds(self, lc1, lc2):
         """
         Helper method to create the power spectrum of both light curves
@@ -983,10 +973,6 @@ class Crossspectrum(object):
         bin_cs.freq = binfreq
         bin_cs.power = bincs
         bin_cs.df = df
-        bin_cs.dt = self.dt
-        bin_cs.n = self.n
-        bin_cs.norm = self.norm
-        bin_cs.nphots1 = self.nphots1
         bin_cs.power_err = binerr
 
         if hasattr(self, "unnorm_power") and self.unnorm_power is not None:
@@ -1013,14 +999,6 @@ class Crossspectrum(object):
         if hasattr(self, "pds2"):
             bin_cs.pds2 = self.pds2.rebin(df=df, f=f, method=method)
 
-        try:
-            bin_cs.nphots2 = self.nphots2
-        except AttributeError:
-            if self.type == "powerspectrum":
-                pass
-            else:
-                raise AttributeError("Spectrum has no attribute named nphots2.")
-
         bin_cs.m = np.rint(step_size * self.m)
 
         return bin_cs
@@ -1719,6 +1697,7 @@ class Crossspectrum(object):
         return
 
 class AveragedCrossspectrum(Crossspectrum):
+    type = "crossspectrum"
     """
     Make an averaged cross spectrum from a light curve by segmenting two
     light curves, Fourier-transforming each segment and then averaging the
diff --git a/stingray/events.py b/stingray/events.py
index a8ef01b1..731a2510 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -13,6 +13,7 @@ import numpy as np
 import numpy.random as ra
 from astropy.table import Table
 
+from .base import StingrayObject, StingrayTimeseries
 from .filters import get_deadtime_mask
 from .gti import append_gtis, check_separate, cross_gtis, generate_indices_of_boundaries
 from .io import load_events_and_gtis
@@ -22,7 +23,7 @@ from .utils import assign_value_if_none, simon, interpret_times
 __all__ = ['EventList']
 
 
-class EventList(object):
+class EventList(StingrayTimeseries):
     """
     Basic class for event list data. Event lists generally correspond to individual events (e.g. photons)
     recorded by the detector, and their associated properties. For X-ray data where this type commonly occurs,
@@ -127,12 +128,13 @@ class EventList(object):
         The full header of the original FITS file, if relevant
 
     """
-
+    main_array_attr = "time"
     def __init__(self, time=None, energy=None, ncounts=None, mjdref=0, dt=0,
                  notes="", gti=None, pi=None, high_precision=False,
                  mission=None, instr=None, header=None, detector_id=None,
                  ephem=None, timeref=None, timesys=None,
                  **other_kw):
+        StingrayObject.__init__(self)
 
         self.energy = None if energy is None else np.asarray(energy)
         self.notes = notes
@@ -166,67 +168,6 @@ class EventList(object):
             if self.time.size != self.energy.size:
                 raise ValueError('Lengths of time and energy must be equal.')
 
-    def array_attrs(self):
-        """List the names of the array attributes of the event list.
-
-        By array attributes, we mean the ones with the same size and shape as
-        ``self.time``
-
-        Examples
-        --------
-        For an empty eventlist, this will be an empty list.
-        >>> ev = EventList()
-        >>> attrs = ev.array_attrs()
-        >>> len(attrs)
-        0
-
-        Otherwise, it will contain all attrs with the same shape as ``time``.
-        >>> ev = EventList(time=np.arange(5), pi=np.zeros(5), gti=[[45, 4]])
-        >>> attrs = ev.array_attrs()
-        >>> len(attrs)
-        2
-        >>> "pi" in attrs and "time" in attrs
-        True
-        >>> hasattr(ev, "gti") #  The gti attribute is set, however...
-        True
-        >>> "gti" in attrs #  The shape of self.gti is not equal to time! Not an array attr
-        False
-        """
-        if self.time is None:
-            return []
-
-        return [
-            attr for attr in dir(self)
-            if (
-                isinstance(getattr(self, attr), Iterable)
-                and np.shape(getattr(self, attr)) == self.time.shape
-            )
-        ]
-
-    def meta_attrs(self):
-        """List the names of the meta attributes of the event list.
-
-        By meta attributes, we mean the ones with a different size and shape
-        than ``self.time``
-
-        Examples
-        --------
-        >>> ev = EventList(time=np.arange(5), pi=np.zeros(5), gti=[[45, 4]])
-        >>> attrs = ev.meta_attrs()
-        >>> "pi" in attrs or "time" in attrs
-        False
-        >>> "gti" in attrs #  The shape of self.gti is not equal to time! Not an array attr
-        True
-        """
-        array_attrs = self.array_attrs()
-        return [
-            attr for attr in dir(self)
-            if (
-                attr not in array_attrs and not attr.startswith("_") and
-                not callable(getattr(self, attr))
-            )
-        ]
-
     def to_lc(self, dt, tstart=None, tseg=None):
         """
         Convert event list to a :class:`stingray.Lightcurve` object.
@@ -556,8 +497,8 @@ class EventList(object):
 
         return ev_new
 
-    @staticmethod
-    def read(filename, format_="pickle", **kwargs):
+    @classmethod
+    def read(cls, filename, fmt=None, format_=None, **kwargs):
         r"""Read a :class:`EventList` object from file.
 
         Currently supported formats are
@@ -581,20 +522,28 @@ class EventList(object):
         filename: str
             Path and file name for the file to be read.
 
-        format\_: str
+        fmt: str
             Available options are 'pickle', 'hea', and any `Table`-supported
             format such as 'hdf5', 'ascii.ecsv', etc.
 
+        Other parameters
+        ----------------
+        kwargs : dict
+            Any further keyword arguments to be passed to `load_events_and_gtis`
+            for reading in event lists in OGIP/HEASOFT format
+
         Returns
         -------
         ev: :class:`EventList` object
             The :class:`EventList` object reconstructed from file
         """
-        if format_.lower() == 'pickle':
-            with open(filename, 'rb') as fobj:
-                return pickle.load(fobj)
+        if fmt is None and format_ is not None:
+            warnings.warn(
+                "The format_ keyword for read and write is deprecated. "
+                "Use fmt instead", DeprecationWarning)
+            fmt = format_
 
-        if format_.lower() in ('hea', 'ogip'):
+        if fmt.lower() in ('hea', 'ogip'):
             evtdata = load_events_and_gtis(filename, **kwargs)
 
             evt = EventList(time=evtdata.ev_list,
@@ -615,45 +564,7 @@ class EventList(object):
                         setattr(evt, key.lower(), evtdata.additional_data[key])
             return evt
 
-        if format_.lower() == 'ascii':
-            format_ = 'ascii.ecsv'
-
-        ts = Table.read(filename, format=format_)
-        return EventList.from_astropy_table(ts)
-
-    def write(self, filename, format_='pickle'):
-        """
-        Write an :class:`EventList` object to file.
-
-        Possible file formats are
-
-        * pickle (not recommended for long-term storage)
-        * any other formats compatible with the writers in
-          :class:`astropy.table.Table` (ascii.ecsv, hdf5, etc.)
-
-        Parameters
-        ----------
-        filename: str
-            Name and path of the file to save the event list to..
-
-        format_: str
-            The file format to store the data in.
-            Available options are ``pickle``, ``hdf5``, ``ascii``, ``fits``
-        """
-        if format_ == 'pickle':
-            with open(filename, "wb") as fobj:
-                pickle.dump(self, fobj)
-            return
-
-        if format_ == 'ascii':
-            format_ = 'ascii.ecsv'
-
-        ts = self.to_astropy_table()
-        try:
-            ts.write(filename, format=format_, overwrite=True,
-                     serialize_meta=True)
-        except TypeError:
-            ts.write(filename, format=format_, overwrite=True)
+        return super().read(filename=filename, fmt=fmt)
 
     def filter_energy_range(self, energy_range, inplace=False, use_pi=False):
         """Filter the event list from a given energy range.
@@ -799,245 +710,3 @@ class EventList(object):
 
         return new_ev
 
-    def change_mjdref(self, new_mjdref):
-        """Change the MJD reference time (MJDREF) of the light curve.
-
-        Times will be now referred to this new MJDREF
-
-        Parameters
-        ----------
-        new_mjdref : float
-            New MJDREF
-
-        Returns
-        -------
-        new_lc : :class:`EventList` object
-            The new LC shifted by MJDREF
-        """
-        time_shift = (self.mjdref - new_mjdref) * 86400
-
-        new_ev = self.shift(time_shift)
-        new_ev.mjdref = new_mjdref
-        return new_ev
-
-    def shift(self, time_shift):
-        """
-        Shift the events and the GTIs in time.
-
-        Parameters
-        ----------
-        time_shift: float
-            The time interval by which the light curve will be shifted (in
-            the same units as the time array in :class:`Lightcurve`
-
-        Returns
-        -------
-        new_ev : lightcurve.Lightcurve object
-            The new event list shifted by ``time_shift``
-
-        """
-        new_ev = copy.deepcopy(self)
-        new_ev.time = new_ev.time + time_shift
-        new_ev.gti = new_ev.gti + time_shift
-
-        return new_ev
-
-    def get_meta_dict(self):
-        """Give a dictionary with all non-None meta attributes."""
-        meta_attrs = self.meta_attrs()
-        meta_dict = {}
-        for key in meta_attrs:
-            val = getattr(self, key)
-            if val is not None:
-                meta_dict[key] = val
-        return meta_dict
-
-    def to_astropy_timeseries(self):
-        """Save the event list to an Astropy timeseries.
-
-        Array attributes (time, pi, energy, etc.) are converted
-        into columns, while meta attributes (mjdref, gti, etc.)
-        are saved into the ``meta`` dictionary.
-        """
-        from astropy.timeseries import TimeSeries
-        from astropy.time import TimeDelta
-        from astropy import units as u
-        data = {}
-        array_attrs = self.array_attrs()
-
-        for attr in array_attrs:
-            if attr == "time":
-                continue
-            data[attr] = np.asarray(getattr(self, attr))
-
-        if data == {}:
-            data = None
-
-        if self.time is not None and self.time.size > 0:
-            times = TimeDelta(self.time * u.s)
-            ts = TimeSeries(data=data, time=times)
-        else:
-            ts = TimeSeries()
-
-        ts.meta.update(self.get_meta_dict())
-
-        return ts
-
-    @staticmethod
-    def from_astropy_timeseries(ts):
-        """Create an `EventList` object from data in an Astropy TimeSeries
-
-        The timeseries has to define at least a column called time,
-        the rest of columns will form the array attributes of the
-        new event list, while the attributes in table.meta will
-        form the new meta attributes of the event list.
-
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of EventList: time, pi, energy, gti etc.
-
-        """
-        kwargs = dict([(key.lower(), val) for (key, val) in ts.meta.items()])
-        ev = EventList(time=ts.time, **kwargs)
-        array_attrs = ts.colnames
-
-        for attr in array_attrs:
-            if attr == "time":
-                continue
-            setattr(ev, attr, ts[attr])
-
-        return ev
-
-    def to_astropy_table(self):
-        """Save the event list to an Astropy Table.
-
-        Array attributes (time, pi, energy, etc.) are converted
-        into columns, while meta attributes (mjdref, gti, etc.)
-        are saved into the ``meta`` dictionary.
-        """
-        data = {}
-        array_attrs = self.array_attrs()
-
-        for attr in array_attrs:
-            data[attr] = np.asarray(getattr(self, attr))
-
-        ts = Table(data)
-
-        ts.meta.update(self.get_meta_dict())
-
-        return ts
-
-    @staticmethod
-    def from_astropy_table(ts):
-        """Create an `EventList` object from data in an Astropy Table.
-
-        The table has to define at least a column called time,
-        the rest of columns will form the array attributes of the
-        new event list, while the attributes in table.meta will
-        form the new meta attributes of the event list.
-
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of EventList: time, pi, energy, gti etc.
-
-        """
-        kwargs = dict([(key.lower(), val) for (key, val) in ts.meta.items()])
-        ev = EventList(time=ts["time"], **kwargs)
-        array_attrs = ts.colnames
-
-        for attr in array_attrs:
-            if attr == "time":
-                continue
-            setattr(ev, attr, ts[attr])
-
-        return ev
-
-    def to_xarray(self):
-        """Save the event list to an xarray Dataset.
-
-        Array attributes (time, pi, energy, etc.) are converted
-        into columns, while meta attributes (mjdref, gti, etc.)
-        are saved into the ``ds.attrs`` dictionary.
-        """
-        from xarray import Dataset
-        data = {}
-        array_attrs = self.array_attrs()
-
-        for attr in array_attrs:
-            data[attr] = np.asarray(getattr(self, attr))
-
-        ts = Dataset(data)
-
-        ts.attrs.update(self.get_meta_dict())
-
-        return ts
-
-    @staticmethod
-    def from_xarray(ts):
-        """Create an `EventList` object from data in an xarray Dataset.
-
-        The dataset has to define at least a column called time,
-        the rest of columns will form the array attributes of the
-        new event list, while the attributes in ds.attrs will
-        form the new meta attributes of the event list.
-
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of EventList: time, pi, energy, gti etc.
-
-        """
-        array_attrs = ts.coords
-
-        kwargs = dict([(key.lower(), val)
-                      for (key, val) in ts.attrs.items() if key not in array_attrs])
-        ev = EventList(time=ts["time"].values, **kwargs)
-
-        for attr in array_attrs:
-            if attr == "time":
-                continue
-            setattr(ev, attr, ts[attr].values)
-
-        return ev
-
-    def to_pandas(self):
-        """Save the event list to a pandas DataFrame.
-
-        Array attributes (time, pi, energy, etc.) are converted
-        into columns, while meta attributes (mjdref, gti, etc.)
-        are saved into the ``ds.attrs`` dictionary.
-        """
-        from pandas import DataFrame
-        data = {}
-        array_attrs = self.array_attrs()
-
-        for attr in array_attrs:
-            data[attr] = np.asarray(getattr(self, attr))
-
-        ts = DataFrame(data)
-
-        ts.attrs.update(self.get_meta_dict())
-
-        return ts
-
-    @staticmethod
-    def from_pandas(ts):
-        """Create an `EventList` object from data in a pandas DataFrame.
-
-        The dataframe has to define at least a column called time,
-        the rest of columns will form the array attributes of the
-        new event list, while the attributes in ds.attrs will
-        form the new meta attributes of the event list.
-
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of EventList: time, pi, energy, gti etc.
-
-        """
-        array_attrs = ts.columns
-
-        kwargs = dict([(key.lower(), val)
-                      for (key, val) in ts.attrs.items() if key not in array_attrs])
-        ev = EventList(time=ts["time"].to_numpy(), **kwargs)
-
-        for attr in array_attrs:
-            if attr == "time":
-                continue
-            setattr(ev, attr, ts[attr].to_numpy())
-
-        return ev
diff --git a/stingray/io.py b/stingray/io.py
index 5332900d..a08fc2bb 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -888,545 +888,6 @@ def split_numbers(number, shift=0):
     return np.double(number_I) / 10**shift, np.double(number_F) / 10**shift
 
 
-def _save_pickle_object(object, filename):
-    """
-    Save a class object in pickle format.
-
-    Parameters
-    ----------
-    object: class instance
-        A class object whose attributes are saved in a
-        dictionary format
-
-    filename: str
-        Name of the file in which object is saved
-    """
-
-    with open(filename, "wb") as f:
-        pickle.dump(object, f)
-
-
-def _retrieve_pickle_object(filename):
-    """
-    Retrieves a pickled class object.
-
-    Parameters
-    ----------
-    filename: str
-        Name of the file in which object is saved
-
-    Returns
-    -------
-    data: class object
-    """
-
-    with open(filename, "rb") as f:
-        return pickle.load(f)
-
-
-def _save_hdf5_object(object, filename):
-    """
-    Save a class object in hdf5 format.
-
-    Parameters
-    ----------
-    object: class instance
-        A class object whose attributes are saved in a
-        dictionary format
-
-    filename: str
-        Name of the file in which object is saved
-    """
-
-    items = vars(object)
-    attrs = [name for name in items if items[name] is not None]
-
-    with h5py.File(filename, 'w') as hf:
-        for attr in attrs:
-            data = items[attr]
-
-            # If data is a single number, store as an attribute.
-            if _isattribute(data):
-                if isinstance(data, np.longdouble):
-                    data_I, data_F = split_numbers(data)
-                    names = [attr + '_I', attr + '_F']
-                    hf.attrs[names[0]] = data_I
-                    hf.attrs[names[1]] = data_F
-                else:
-                    hf.attrs[attr] = data
-
-            # If data is an array or list, create a dataset.
-            else:
-                try:
-                    if isinstance(data[0], np.longdouble):
-                        data_I, data_F = split_numbers(data)
-                        names = [attr + '_I', attr + '_F']
-                        hf.create_dataset(names[0], data=data_I)
-                        hf.create_dataset(names[1], data=data_F)
-                    else:
-                        hf.create_dataset(attr, data=data)
-                except IndexError:
-                    # To account for numpy arrays of type 'None' (0-d)
-                    pass
-
-
-def _retrieve_hdf5_object(filename):
-    """
-    Retrieves an hdf5 format class object.
-
-    Parameters
-    ----------
-    filename: str
-        The name of file with which object was saved
-
-    Returns
-    -------
-    data: dictionary
-        Loads the data from an hdf5 object file and returns
-        in dictionary format.
-    """
-
-    with h5py.File(filename, 'r') as hf:
-        dset_keys = hf.keys()
-        attr_keys = hf.attrs.keys()
-        data = {}
-
-        dset_copy = list(dset_keys)[:]
-        for key in dset_keys:
-
-            # Make sure key hasn't been removed
-            if key in dset_copy:
-                # Longdouble case
-                if key[-2:] in ['_I', '_F']:
-                    m_key = key[:-2]
-                    # Add integer and float parts
-                    data[m_key] = np.longdouble(hf[m_key + '_I'][()])
-                    data[m_key] += np.longdouble(hf[m_key + '_F'][()])
-                    # Remove integer and float parts from attributes
-                    dset_copy.remove(m_key + '_I')
-                    dset_copy.remove(m_key + '_F')
-                else:
-                    data[key] = hf[key][()]
-
-        attr_copy = list(attr_keys)[:]
-        for key in attr_keys:
-
-            # Make sure key hasn't been removed
-            if key in attr_copy:
-                # Longdouble case
-                if key[-2:] in ['_I', '_F']:
-                    m_key = key[:-2]
-                    # Add integer and float parts
-                    data[m_key] = np.longdouble(hf.attrs[m_key + '_I'])
-                    data[m_key] += np.longdouble(hf.attrs[m_key + '_F'])
-                    # Remove integer and float parts from attributes
-                    attr_copy.remove(m_key + '_I')
-                    attr_copy.remove(m_key + '_F')
-                else:
-                    data[key] = hf.attrs[key]
-
-    return data
-
-
-def _save_ascii_object(object, filename, fmt="%.18e", **kwargs):
-    """
-    Save an array to a text file.
-
-    Parameters
-    ----------
-    object : numpy.ndarray
-        An array with the data to be saved
-
-    filename : str
-        The file name to save to
-
-    fmt : str or sequence of strs, optional
-        Use for formatting of columns. See `numpy.savetxt` documentation
-        for details.
-
-    Other Parameters
-    ----------------
-    kwargs : any keyword argument taken by `numpy.savetxt`
-
-    """
-
-    try:
-        np.savetxt(filename, object, fmt=fmt, **kwargs)
-    except TypeError:
-        raise Exception("Formatting of columns not recognized! Use 'fmt' "
-                        "to format columns including strings or mixed types!")
-
-    pass
-
-
-def _retrieve_ascii_object(filename, **kwargs):
-    """
-    Helper function to retrieve ascii objects from file.
-    Uses astropy.Table for reading and storing the data.
-
-    Parameters
-    ----------
-    filename : str
-        The name of the file with the data to be retrieved.
-
-    Other Parameters
-    -----------------------------
-    usecols : {int | iterable}
-        The indices of the columns in the file to be returned.
-        By default, all columns will be returned
-
-    skiprows : int
-        The number of rows at the beginning to skip
-        By default, no rows will be skipped.
-
-    names : iterable
-        A list of column names to be attached to the columns.
-        By default, no column names are added, unless they are specified
-        in the file header and can be read by astropy.Table.read
-        automatically.
-
-    Returns
-    -------
-    data : astropy.Table object
-        An astropy.Table object with the data from the file
-    """
-    if not isinstance(filename, str):
-        raise TypeError("filename must be string!")
-
-    if 'usecols' in list(kwargs.keys()):
-        if np.size(kwargs['usecols']) != 2:
-            raise ValueError("Need to define two columns")
-        usecols = kwargs["usecols"]
-    else:
-        usecols = None
-
-    if 'skiprows' in list(kwargs.keys()):
-        assert isinstance(kwargs["skiprows"], int)
-        skiprows = kwargs["skiprows"]
-    else:
-        skiprows = 0
-
-    if "names" in list(kwargs.keys()):
-        names = kwargs["names"]
-    else:
-        names = None
-
-    data = Table.read(filename, data_start=skiprows,
-                      names=names, format="ascii")
-
-    if usecols is None:
-        return data
-    else:
-        colnames = np.array(data.colnames)
-        cols = colnames[usecols]
-
-        return data[cols]
-
-
-def _save_fits_object(object, filename, **kwargs):
-    """
-    Save a class object in fits format.
-
-    Parameters
-    ----------
-    object: class instance
-        A class object whose attributes would be saved in a dictionary format.
-
-    filename: str
-        The file name to save to
-
-    Additional Keyword Parameters
-    -----------------------------
-    tnames: str iterable
-        The names of HDU tables. For instance, in case of eventlist,
-        tnames could be ['EVENTS', 'GTI']
-
-    colsassign: dictionary iterable
-        This indicates the correct tables to which to assign columns
-        to. If this is None or if a column is not provided, it/they will
-        be assigned to the first table.
-
-        For example, [{'gti':'GTI'}] indicates that gti values should be
-        stored in GTI table.
-    """
-
-    tables = []
-
-    if 'colsassign' in list(kwargs.keys()):
-        colsassign = kwargs['colsassign']
-        iscolsassigned = True
-    else:
-        iscolsassigned = False
-
-    if 'tnames' in list(kwargs.keys()):
-        tables = kwargs['tnames']
-    else:
-        tables = ['MAIN']
-
-    items = vars(object)
-    attrs = [name for name in items if items[name] is not None]
-
-    cols = []
-    hdrs = []
-
-    for t in tables:
-        cols.append([])
-        hdrs.append(fits.Header())
-
-    for attr in attrs:
-        data = items[attr]
-
-        # Get the index of table to which column belongs
-        if iscolsassigned and attr in colsassign.keys():
-            index = tables.index(colsassign[attr])
-        else:
-            index = 0
-
-        # If data is a single number, store as metadata
-        if _isattribute(data):
-            if isinstance(data, np.longdouble):
-                # Longdouble case. Split and save integer and float parts
-                data_I, data_F = split_numbers(data)
-                names = [attr + '_I', attr + '_F']
-                hdrs[index][names[0]] = data_I
-                hdrs[index][names[1]] = data_F
-            else:
-                # Normal case. Save as it is
-                hdrs[index][attr] = data
-
-        # If data is an array or list, insert as table column
-        else:
-            try:
-                if isinstance(data[0], np.longdouble):
-                    # Longdouble case. Split and save integer and float parts
-                    data_I, data_F = split_numbers(data)
-                    names = [attr + '_I', attr + '_F']
-                    cols[index].append(
-                        fits.Column(name=names[0],
-                                    format='D',
-                                    array=data_I))
-                    cols[index].append(
-                        fits.Column(name=names[1],
-                                    format='D',
-                                    array=data_F))
-                else:
-                    # Normal case. Save as it is
-                    cols[index].append(
-                        fits.Column(name=attr,
-                                    format=_lookup_format(data[0]),
-                                    array=data))
-            except IndexError:
-                # To account for numpy arrays of type 'None' (0-d)
-                pass
-
-    tbhdu = fits.HDUList()
-
-    # Create binary tables
-    for i in range(0, len(tables)):
-        if len(cols[i]) > 0:
-            tbhdu.append(fits.BinTableHDU.from_columns(cols[i],
-                                                       header=hdrs[i],
-                                                       name=tables[i]))
-
-    tbhdu.writeto(filename)
-
-
-def _retrieve_fits_object(filename, **kwargs):
-    """
-    Retrieves a fits format class object.
-
-    Parameters
-    ----------
-    filename: str
-        The name of file with which object was saved
-
-    Other Parameters
-    ----------------
-    cols: str iterable
-        The names of columns to extract from fits tables.
-
-    Returns
-    -------
-    data: dictionary
-        Loads the data from a fits object file and returns
-        in dictionary format.
-    """
-
-    data = {}
-
-    if 'cols' in list(kwargs.keys()):
-        cols = [col.upper() for col in kwargs['cols']]
-    else:
-        cols = []
-
-    with fits.open(filename, memmap=False, ignore_missing_end=True) as hdulist:
-        fits_cols = []
-
-        # Get columns from all tables
-        for i in range(1, len(hdulist)):
-            fits_cols.append([h.upper() for h in hdulist[i].data.names])
-
-        for c in cols:
-            for i in range(0, len(fits_cols)):
-                # .upper() is used because `fits` stores values in upper case
-                hdr_keys = [h.upper() for h in hdulist[i + 1].header.keys()]
-
-                # Longdouble case. Check for columns
-                if c + '_I' in fits_cols[i] or c + '_F' in fits_cols[i]:
-                    if c not in data.keys():
-                        data[c] = np.longdouble(hdulist[i + 1].data[c + '_I'])
-                        data[c] += np.longdouble(hdulist[i + 1].data[c + '_F'])
-
-                # Longdouble case. Check for header keys
-                if c + '_I' in hdr_keys or c + '_F' in hdr_keys:
-                    if c not in data.keys():
-                        data[c] = \
-                            np.longdouble(hdulist[i + 1].header[c + '_I'])
-                        data[c] += \
-                            np.longdouble(hdulist[i + 1].header[c + '_F'])
-
-                # Normal case. Check for columns
-                elif c in fits_cols[i]:
-                    data[c] = hdulist[i + 1].data[c]
-
-                # Normal case. Check for header keys
-                elif c in hdr_keys:
-                    data[c] = hdulist[i + 1].header[c]
-        hdulist.close()
-    return data
-
-
-def _lookup_format(var):
-    """
-    Looks up relevant format in fits.
-
-    Parameters
-    ----------
-    var : object
-        An object to look up in the table
-
-    Returns
-    -------
-    lookup : str
-        The str describing the type of ``var``
-    """
-
-    lookup = {"<type 'int'>": "J", "<type 'float'>": "E",
-              "<type 'numpy.int64'>": "K", "<type 'numpy.float64'>": "D",
-              "<type 'numpy.float128'>": "D", "<type 'str'>": "30A",
-              "<type 'bool'": "L"}
-
-    form = type(var)
-
-    try:
-        return lookup[str(form)]
-    except KeyError:
-        # If an entry is not contained in lookup dictionary
-        return "D"
-
-
-def _isattribute(data):
-    """
-    Check if data is a single number or an array.
-
-    Parameters
-    ----------
-    data : object
-        The object to be checked.
-
-    Returns:
-        bool
-        True if the data is a single number, False if it is an iterable.
-    """
-
-    if isinstance(data, Iterable) and not isinstance(data, (str, bytes)):
-        return False
-    else:
-        return True
-
-
-def write(input_, filename, format_='pickle', **kwargs):
-    """
-    Pickle a class instance. For parameters depending on
-    ``format_``, see individual function definitions.
-
-    Parameters
-    ----------
-    object: a class instance
-        The object to be stored
-
-    filename: str
-        The name of the file to be created
-
-    format_: str
-        The format in which to store file. Formats supported
-        are ``pickle``, ``hdf5``, ``ascii`` or ``fits``
-    """
-
-    if format_ == 'pickle':
-        _save_pickle_object(input_, filename)
-
-    elif format_ == 'hdf5':
-        if _H5PY_INSTALLED:
-            _save_hdf5_object(input_, filename)
-        else:
-            utils.simon('h5py not installed, using pickle instead'
-                        'to save object.')
-            _save_pickle_object(input_, filename.split('.')[0] +
-                                '.pickle')
-
-    elif format_ == 'ascii':
-        _save_ascii_object(input_, filename, **kwargs)
-
-    elif format_ == 'fits':
-        _save_fits_object(input_, filename, **kwargs)
-
-    else:
-        utils.simon('Format not understood.')
-
-
-def read(filename, format_='pickle', **kwargs):
-    """
-    Return a saved class instance.
-
-    Parameters
-    ----------
-    filename: str
-        The name of the file to be retrieved.
-
-    format_: str
-        The format used to store file. Supported formats are
-        pickle, hdf5, ascii or fits.
-
-    Returns
-    -------
-    data : {``object`` | ``astropy.table`` | ``dict``}
-
-        * If ``format_`` is ``pickle``, an object is returned.
-        * If ``format_`` is ``ascii``, `astropy.table` object is returned.
-        * If ``format_`` is ``hdf5`` or 'fits``, a dictionary object is returned.
-    """
-
-    if format_ == 'pickle':
-        return _retrieve_pickle_object(filename)
-
-    elif format_ == 'hdf5':
-        if _H5PY_INSTALLED:
-            return _retrieve_hdf5_object(filename)
-        else:
-            utils.simon('h5py not installed, cannot read an'
-                        'hdf5 object.')
-
-    elif format_ == 'ascii':
-        return _retrieve_ascii_object(filename, **kwargs)
-
-    elif format_ == 'fits':
-        return _retrieve_fits_object(filename, **kwargs)
-
-    else:
-        utils.simon('Format not understood.')
-
-
 def savefig(filename, **kwargs):
     """
     Save a figure plotted by ``matplotlib``.
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 1fb5c5ab..e13803ee 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -15,6 +15,7 @@ from astropy.table import Table
 from astropy.time import TimeDelta, Time
 from astropy import units as u
 
+from stingray.base import StingrayTimeseries
 import stingray.utils as utils
 from stingray.exceptions import StingrayError
 from stingray.gti import (bin_intervals_from_gtis, check_gtis, create_gti_mask,
@@ -28,7 +29,7 @@ __all__ = ["Lightcurve"]
 valid_statistics = ["poisson", "gauss", None]
 
 
-class Lightcurve(object):
+class Lightcurve(StingrayTimeseries):
     """
     Make a light curve object from an array of time stamps and an
     array of counts.
@@ -164,12 +165,15 @@ class Lightcurve(object):
         The full header of the original FITS file, if relevant
 
     """
+    main_array_attr = "time"
 
     def __init__(self, time, counts, err=None, input_counts=True,
                  gti=None, err_dist='poisson', mjdref=0, dt=None,
                  skip_checks=False, low_memory=False, mission=None,
                  instr=None, header=None, **other_kw):
 
+        StingrayTimeseries.__init__(self)
+
         if other_kw != {}:
             warnings.warn(f"Unrecognized keywords: {list(other_kw.keys())}")
 
@@ -474,54 +478,6 @@ class Lightcurve(object):
                   "This could cause problems with Fourier transforms. "
                   "Please make the input time evenly sampled.")
 
-    def change_mjdref(self, new_mjdref):
-        """Change the MJD reference time (MJDREF) of the light curve.
-
-        Times will be now referred to this new MJDREF
-
-        Parameters
-        ----------
-        new_mjdref : float
-            New MJDREF
-
-        Returns
-        -------
-        new_lc : lightcurve.Lightcurve object
-            The new LC shifted by MJDREF
-        """
-        time_shift = -(new_mjdref - self.mjdref) * 86400
-
-        new_lc = self.shift(time_shift)
-        new_lc.mjdref = new_mjdref
-        return new_lc
-
-    def shift(self, time_shift):
-        """
-        Shift the light curve and the GTIs in time.
-
-        Parameters
-        ----------
-        time_shift: float
-            The time interval by which the light curve will be shifted (in
-            the same units as the time array in :class:`Lightcurve`
-
-        Returns
-        -------
-        new_lc : lightcurve.Lightcurve object
-            The new LC shifted by ``time_shift``
-
-        """
-        new_lc = Lightcurve(self.time + time_shift,
-                            self.counts,
-                            err=self._counts_err,
-                            gti=self.gti + time_shift,
-                            mjdref=self.mjdref,
-                            dt=self.dt,
-                            err_dist=self.err_dist,
-                            skip_checks=True)
-
-        return new_lc
-
     def _operation_with_other_lc(self, other, operation):
         """
         Helper method to codify an operation of one light curve with another (e.g. add, subtract, ...).
@@ -1639,39 +1595,8 @@ class Lightcurve(object):
         else:
             plt.show(block=False)
 
-    def write(self, filename, format_='pickle', **kwargs):
-        """
-        Write a :class:`Lightcurve` object to file. Currently supported formats are
-
-        * pickle (not recommended for long-term storage)
-        * HDF5
-        * ASCII
-
-        Parameters
-        ----------
-        filename: str
-            Path and file name for the output file.
-
-        format\_: str
-            Available options are 'pickle', 'hdf5', 'ascii'
-        """
-        if format_ == 'pickle':
-            with open(filename, "wb") as fobj:
-                pickle.dump(self, fobj)
-            return
-
-        if format_ == 'ascii':
-            format_ = 'ascii.ecsv'
-
-        ts = self.to_astropy_table()
-        try:
-            ts.write(filename, format=format_, overwrite=True,
-                     serialize_meta=True)
-        except TypeError:
-            ts.write(filename, format=format_, overwrite=True)
-
-    @staticmethod
-    def read(filename, format_='pickle', err_dist='gauss',
+    @classmethod
+    def read(cls, filename, fmt=None, format_=None, err_dist='gauss',
              skip_checks=False):
         """
         Read a :class:`Lightcurve` object from file.
@@ -1696,7 +1621,7 @@ class Lightcurve(object):
         filename: str
             Path and file name for the file to be read.
 
-        format\_: str
+        fmt: str
             Available options are 'pickle', 'hea', and any `Table`-supported
             format such as 'hdf5', 'ascii.ecsv', etc.
 
@@ -1715,21 +1640,18 @@ class Lightcurve(object):
         --------
         lc : :class:`Lightcurve` object
         """
-        if format_ == 'pickle':
-            with open(filename, 'rb') as fobj:
-                return pickle.load(fobj)
+        if fmt is None and format_ is not None:
+            warnings.warn(
+                "The format_ keyword for read and write is deprecated. "
+                "Use fmt instead", DeprecationWarning)
+            fmt = format_
 
-        if format_ == 'hea':
+        if fmt.lower() in ('hea', 'ogip'):
             data = lcurve_from_fits(filename)
             data.update({'err_dist': err_dist, 'skip_checks': skip_checks})
             return Lightcurve(**data)
 
-        if format_ == 'ascii':
-            format_ = 'ascii.ecsv'
-
-        ts = Table.read(filename, format=format_)
-        return Lightcurve.from_astropy_table(
-            ts, err_dist=err_dist, skip_checks=skip_checks)
+        return super().read(filename=filename, fmt=fmt)
 
     def split_by_gti(self, gti=None, min_points=2):
         """
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 057e8afd..e6aecf83 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -32,6 +32,7 @@ __all__ = ["Powerspectrum", "AveragedPowerspectrum", "DynamicalPowerspectrum"]
 
 
 class Powerspectrum(Crossspectrum):
+    type = "powerspectrum"
     """
     Make a :class:`Powerspectrum` (also called periodogram) from a (binned) light curve.
     Periodograms can be normalized by either Leahy normalization, fractional rms
@@ -619,7 +620,9 @@ class Powerspectrum(Crossspectrum):
         self.n = None
         return
 
+
 class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
+    type = "powerspectrum"
     """
     Make an averaged periodogram from a light curve by segmenting the light
     curve, Fourier-transforming each segment and then averaging the
@@ -877,6 +880,7 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
 
 
 class DynamicalPowerspectrum(AveragedPowerspectrum):
+    type = "powerspectrum"
     """
     Create a dynamical power spectrum, also often called a *spectrogram*.
 
diff --git a/stingray/simulator/simulator.py b/stingray/simulator/simulator.py
index 589a3713..5c776fb4 100644
--- a/stingray/simulator/simulator.py
+++ b/stingray/simulator/simulator.py
@@ -1,11 +1,10 @@
-
+import pickle
 from os import error
 import numpy as np
 import numbers
 import warnings
 from scipy import signal
 import astropy.modeling.models
-from stingray.io import write, read
 from stingray import utils
 from stingray import Lightcurve
 from stingray import AveragedPowerspectrum
@@ -591,41 +590,43 @@ class Simulator(object):
         return AveragedPowerspectrum(lc, seg_size).power
 
     @staticmethod
-    def read(filename, format_='pickle'):
+    def read(filename, fmt='pickle', format_=None):
         """
-        Imports Simulator object.
-
-        Parameters
-        ----------
-        filename : str
-            Name of the Simulator object to be read.
+        Reads transfer function from a 'pickle' file.
 
-        format_ : str
-            Available option is 'pickle.'
+        Parameter
+        ---------
+        fmt : str
+            the format of the file to be retrieved - accepts 'pickle'.
 
         Returns
         -------
-        object : `Simulator` object
+        data : class instance
+            `TransferFunction` object
         """
-        if format_ == 'pickle':
-            data = read(filename, 'pickle')
-            return data
+        if format_ is not None:
+            fmt = format_
+
+        if fmt == 'pickle':
+            with open(filename, "rb") as fobj:
+                return pickle.load(fobj)
+
         else:
-            raise KeyError("Format not supported.")
+            raise KeyError("Format not understood.")
 
-    def write(self, filename, format_='pickle'):
+    def write(self, filename, fmt="pickle", format_=None):
         """
-        Exports Simulator object.
+        Writes a transfer function to 'pickle' file.
 
         Parameters
         ----------
-        filename : str
-            Name of the Simulator object to be created.
-
-        format_ : str
-            Available options are 'pickle' and 'hdf5'.
-        """
-        if format_ == 'pickle':
-            write(self, filename)
+        fmt : str
+            the format of the file to be saved - accepts 'pickle'
+        """
+        if format_ is not None:
+            fmt = format_
+        if fmt == 'pickle':
+            with open(filename, "wb") as fobj:
+                pickle.dump(self, fobj)
         else:
-            raise KeyError("Format not supported.")
+            raise KeyError("Format not understood.")
diff --git a/stingray/simulator/transfer.py b/stingray/simulator/transfer.py
index 20b2a04b..5dd57725 100644
--- a/stingray/simulator/transfer.py
+++ b/stingray/simulator/transfer.py
@@ -1,8 +1,6 @@
-
+import pickle
 import numpy as np
 
-from stingray.io import read, write
-
 """
 Implementation of time and energy averaged responses from 2-d
 transfer functions.
@@ -181,13 +179,13 @@ class TransferFunction(object):
             plt.close()
 
     @staticmethod
-    def read(filename, format_='pickle'):
+    def read(filename, fmt='pickle', format_=None):
         """
         Reads transfer function from a 'pickle' file.
 
         Parameter
         ---------
-        format\_ : str
+        fmt : str
             the format of the file to be retrieved - accepts 'pickle'.
 
         Returns
@@ -195,28 +193,30 @@ class TransferFunction(object):
         data : class instance
             `TransferFunction` object
         """
+        if format_ is not None:
+            fmt = format_
 
-        object = read(filename, format_)
-
-        if format_ == 'pickle':
-            return object
+        if fmt == 'pickle':
+            with open(filename, "rb") as fobj:
+                return pickle.load(fobj)
 
         else:
             raise KeyError("Format not understood.")
 
-    def write(self, filename, format_='pickle'):
+    def write(self, filename, fmt="pickle", format_=None):
         """
         Writes a transfer function to 'pickle' file.
 
         Parameters
         ----------
-        format\_ : str
+        fmt : str
             the format of the file to be saved - accepts 'pickle'
         """
-
-        if format_ == 'pickle':
-            write(self, filename, format_)
-
+        if format_ is not None:
+            fmt = format_
+        if fmt == 'pickle':
+            with open(filename, "wb") as fobj:
+                pickle.dump(self, fobj)
         else:
             raise KeyError("Format not understood.")
 
diff --git a/stingray/utils.py b/stingray/utils.py
index 78951221..2b6123fd 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -10,11 +10,9 @@ from collections.abc import Iterable
 
 import numpy as np
 import scipy
-from astropy.time import Time, TimeDelta
-import astropy.units as u
-from astropy.units import Quantity
 from numpy import histogram as histogram_np
 from numpy import histogram2d as histogram2d_np
+from .base import interpret_times
 
 
 # If numba is installed, import jit. Otherwise, define an empty decorator with
@@ -219,7 +217,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
     xbin = np.arange(xedges[0], xedges[-1]+dx_new, dx_new)
 
     output = np.zeros(xbin.shape[0] - 1, dtype=type(y[0]))
-    outputerr = np.zeros(xbin.shape[0] - 1, dtype=type(y[0]))
+    outputerr = np.zeros(xbin.shape[0] - 1, dtype=type(yerr[0]))
     step_size = np.zeros(xbin.shape[0] - 1)
 
     all_x = np.searchsorted(xedges, xbin)
@@ -1062,84 +1060,6 @@ def genDataPath(dir_path):
         raise IOError(("Directory does not exist."))
 
 
-def interpret_times(time, mjdref=0):
-    """Get time interval in seconds from an astropy Time object
-
-    Examples
-    --------
-    >>> time = Time(57483, format='mjd')
-    >>> newt, mjdref = interpret_times(time)
-    >>> newt == 0
-    True
-    >>> mjdref == 57483
-    True
-    >>> time = Time([57483], format='mjd')
-    >>> newt, mjdref = interpret_times(time)
-    >>> np.allclose(newt, 0)
-    True
-    >>> mjdref == 57483
-    True
-    >>> time = TimeDelta([3, 4, 5] * u.s)
-    >>> newt, mjdref = interpret_times(time)
-    >>> np.allclose(newt, [3, 4, 5])
-    True
-    >>> time = np.array([3, 4, 5])
-    >>> newt, mjdref = interpret_times(time, mjdref=45000)
-    >>> np.allclose(newt, [3, 4, 5])
-    True
-    >>> mjdref == 45000
-    True
-    >>> time = np.array([3, 4, 5] * u.s)
-    >>> newt, mjdref = interpret_times(time, mjdref=45000)
-    >>> np.allclose(newt, [3, 4, 5])
-    True
-    >>> mjdref == 45000
-    True
-    >>> newt, mjdref = interpret_times(1, mjdref=45000)
-    >>> newt == 1
-    True
-    >>> newt, mjdref = interpret_times(list, mjdref=45000)
-    Traceback (most recent call last):
-    ...
-    ValueError: Unknown time format: ...
-    >>> newt, mjdref = interpret_times("guadfkljfd", mjdref=45000)
-    Traceback (most recent call last):
-    ...
-    ValueError: Unknown time format: ...
-    """
-    if isinstance(time, TimeDelta):
-        out_times = time.to('s').value
-        return out_times, mjdref
-
-    if isinstance(time, Time):
-        mjds = time.mjd
-        if mjdref == 0:
-            if np.all(mjds > 10000):
-                if isinstance(mjds, Iterable):
-                    mjdref = mjds[0]
-                else:
-                    mjdref = mjds
-
-        out_times = (mjds - mjdref) * 86400
-        return out_times, mjdref
-
-    if isinstance(time, Quantity):
-        out_times = time.to('s').value
-        return out_times, mjdref
-
-    if isinstance(time, (tuple, list, np.ndarray)):
-        return time, mjdref
-
-    if not isinstance(time, Iterable):
-        try:
-            float(time)
-            return time, mjdref
-        except (ValueError, TypeError):
-            pass
-
-    raise ValueError(f"Unknown time format: {type(time)}")
-
-
 def check_iterables_close(iter0, iter1, **kwargs):
     """Check that the values produced by iterables are equal.
 
diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py
index d67e8d8f..e5060a56 100644
--- a/stingray/varenergyspectrum.py
+++ b/stingray/varenergyspectrum.py
@@ -1,5 +1,6 @@
 import numpy as np
 import warnings
+from stingray.base import StingrayObject
 from stingray.gti import check_separate, cross_two_gtis
 
 from stingray.lightcurve import Lightcurve
@@ -123,7 +124,8 @@ def _decode_energy_specification(energy_spec):
     return energies
 
 
-class VarEnergySpectrum(metaclass=ABCMeta):
+class VarEnergySpectrum(StingrayObject, metaclass=ABCMeta):
+    main_array_attr = "energy"
     """
     Base class for variability-energy spectrum.
 
@@ -412,6 +414,18 @@ class VarEnergySpectrum(metaclass=ABCMeta):
     def _spectrum_function(self):
         pass
 
+    def from_astropy_table(self, *args, **kwargs):
+        raise NotImplementedError(
+            "from_XXXX methods are not implemented for VarEnergySpectrum")
+
+    def from_xarray(self, *args, **kwargs):
+        raise NotImplementedError(
+            "from_XXXX methods are not implemented for VarEnergySpectrum")
+
+    def from_pandas(self, *args, **kwargs):
+        raise NotImplementedError(
+            "from_XXXX methods are not implemented for VarEnergySpectrum")
+
 
 class RmsSpectrum(VarEnergySpectrum):
     """Calculate the rms-Energy spectrum.
diff --git a/tox.ini b/tox.ini
index a860e2f3..7b7030ce 100644
--- a/tox.ini
+++ b/tox.ini
@@ -57,6 +57,7 @@ deps =
 
     astropy3: astropy>=3.0.0,<4.0.0
     astropy4: astropy>=4.0.0,<5.0.0
+    astropy4: pyyaml
     astropy5: astropy>=5.0.0,<6.0.0
     astropylts: astropy~=5.0.0
 


Test patch:
diff --git a/stingray/simulator/tests/test_simulator.py b/stingray/simulator/tests/test_simulator.py
index 7a5e565f..67699acb 100644
--- a/stingray/simulator/tests/test_simulator.py
+++ b/stingray/simulator/tests/test_simulator.py
@@ -615,8 +615,9 @@ class TestSimulator(object):
     def test_io_with_unsupported_format(self):
         sim = Simulator(N=self.N, dt=self.dt, rms=self.rms, mean=self.mean)
         with pytest.raises(KeyError):
+            # Also use the deprecated format_, just because
             sim.write('sim.hdf5', format_='hdf5')
+        sim.write('sim.pickle', fmt='pickle')
         with pytest.raises(KeyError):
-            sim.write('sim.pickle', format_='pickle')
             sim.read('sim.pickle', format_='hdf5')
         os.remove('sim.pickle')
diff --git a/stingray/simulator/tests/test_transfer.py b/stingray/simulator/tests/test_transfer.py
index 19533d2d..97bd55d4 100644
--- a/stingray/simulator/tests/test_transfer.py
+++ b/stingray/simulator/tests/test_transfer.py
@@ -34,7 +34,7 @@ class TestSimulator(object):
         self.transfer.time_response(e0=3.3, e1=4.7)
 
     def test_time_response_with_incorrect_ranges(self):
-        """Test that incorrect energy ranges raises a 
+        """Test that incorrect energy ranges raises a
         Value Error.
         """
         with pytest.raises(ValueError):
@@ -72,8 +72,8 @@ class TestSimulator(object):
         os.remove('response.png')
 
     def test_io_with_pickle(self):
-        self.transfer.write('transfer.pickle', format_='pickle')
-        tr = self.transfer.read('transfer.pickle', format_='pickle')
+        self.transfer.write('transfer.pickle', fmt='pickle')
+        tr = self.transfer.read('transfer.pickle', fmt='pickle')
         assert (tr.data == self.transfer.data).all()
         os.remove('transfer.pickle')
 
diff --git a/stingray/tests/test_base.py b/stingray/tests/test_base.py
new file mode 100644
index 00000000..97425aa1
--- /dev/null
+++ b/stingray/tests/test_base.py
@@ -0,0 +1,228 @@
+import os
+import copy
+import pytest
+import numpy as np
+from stingray.base import StingrayObject, StingrayTimeseries
+
+_HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = True
+
+try:
+    import xarray
+    from xarray import Dataset
+except ImportError:
+    _HAS_XARRAY = False
+
+try:
+    import pandas
+    from pandas import DataFrame
+except ImportError:
+    _HAS_PANDAS = False
+
+try:
+    import h5py
+except ImportError:
+    _HAS_H5PY = False
+
+
+class DummyStingrayObj(StingrayObject):
+    main_array_attr = "guefus"
+
+    def __init__(self, dummy=None):
+        self.guefus = dummy
+        # StingrayObject.__init__(self)
+
+
+class DummyStingrayTs(StingrayTimeseries):
+    main_array_attr = "time"
+
+    def __init__(self, dummy=None):
+        self.guefus = dummy
+        # StingrayObject.__init__(self)
+
+
+# A StingrayObject with no main_array_attr
+class BadStingrayObj(StingrayObject):
+    def __init__(self, dummy=None):
+        self.guefus = dummy
+        StingrayObject.__init__(self)
+
+
+def _check_equal(so, new_so):
+    for attr in ["time", "guefus", "pardulas", "panesapa"]:
+        if not hasattr(so, attr):
+            assert not hasattr(new_so, attr)
+            continue
+        so_attr = at if (at := getattr(so, attr)) is not None else []
+        new_so_attr = at if (at := getattr(new_so, attr)) is not None else []
+
+        assert np.allclose(so_attr, new_so_attr)
+
+    for attr in ["mjdref", "pirichitus", "parafritus"]:
+        if not hasattr(so, attr):
+            assert not hasattr(new_so, attr)
+            continue
+        so_attr = getattr(so, attr)
+        new_so_attr = getattr(new_so, attr)
+        assert so_attr == new_so_attr
+
+
+class TestStingrayObject:
+    @classmethod
+    def setup_class(cls):
+        cls.arr = [4, 5, 2]
+        sting_obj = DummyStingrayObj(cls.arr)
+        sting_obj.pardulas = [3.0 + 1.0j, 2.0j, 1.0 + 0.0j]
+        sting_obj.pirichitus = 4
+        sting_obj.parafritus = "bonus!"
+        sting_obj.panesapa = [[41, 25], [98, 3]]
+        cls.sting_obj = sting_obj
+
+    def test_preliminary(self):
+        assert np.allclose(self.sting_obj.guefus, self.arr)
+
+    def test_instantiate_without_main_array_attr(self):
+        with pytest.raises(RuntimeError):
+            BadStingrayObj(self.arr)
+
+    def test_side_effects(self):
+        so = copy.deepcopy(self.sting_obj)
+        assert np.allclose(so.guefus, [4, 5, 2])
+        so.guefus = np.random.randint(0, 4, 3)
+        so.panesapa = np.random.randint(5, 9, (6, 2))
+        new_so = DummyStingrayObj([4, 5, 6])
+        assert not hasattr(new_so, "panesapa")
+
+    def test_astropy_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        assert np.allclose(so.guefus, [4, 5, 2])
+        so.guefus = np.random.randint(0, 4, 3)
+        # Set an attribute to a DummyStingrayObj. It will *not* be saved
+        so.stingattr = DummyStingrayObj([3, 4, 5])
+        ts = so.to_astropy_table()
+        new_so = DummyStingrayObj.from_astropy_table(ts)
+        _check_equal(so, new_so)
+        assert not hasattr(new_so, "stingattr")
+
+    @pytest.mark.skipif("not _HAS_XARRAY")
+    def test_xarray_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        assert np.allclose(so.guefus, [4, 5, 2])
+        so.guefus = np.random.randint(0, 4, 3)
+        ts = so.to_xarray()
+        new_so = DummyStingrayObj.from_xarray(ts)
+
+        _check_equal(so, new_so)
+
+    @pytest.mark.skipif("not _HAS_PANDAS")
+    def test_pandas_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        assert np.allclose(so.guefus, [4, 5, 2])
+        so.guefus = np.random.randint(0, 4, 3)
+        ts = so.to_pandas()
+        new_so = DummyStingrayObj.from_pandas(ts)
+
+        _check_equal(so, new_so)
+
+    def test_astropy_roundtrip_empty(self):
+        # Set an attribute to a DummyStingrayObj. It will *not* be saved
+        so = DummyStingrayObj([])
+        ts = so.to_astropy_table()
+        new_so = DummyStingrayObj.from_astropy_table(ts)
+        assert new_so.guefus is None or len(new_so.guefus) == 0
+
+    @pytest.mark.skipif("not _HAS_XARRAY")
+    def test_xarray_roundtrip_empty(self):
+        so = DummyStingrayObj([])
+        ts = so.to_xarray()
+        new_so = DummyStingrayObj.from_xarray(ts)
+
+        assert new_so.guefus is None or len(new_so.guefus) == 0
+
+    @pytest.mark.skipif("not _HAS_PANDAS")
+    def test_pandas_roundtrip_empty(self):
+        so = DummyStingrayObj([])
+        ts = so.to_pandas()
+        new_so = DummyStingrayObj.from_pandas(ts)
+
+        assert new_so.guefus is None or len(new_so.guefus) == 0
+
+    @pytest.mark.skipif("not _HAS_H5PY")
+    def test_hdf_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        so.write("dummy.hdf5")
+        new_so = so.read("dummy.hdf5")
+        os.unlink("dummy.hdf5")
+
+        _check_equal(so, new_so)
+
+    @pytest.mark.skipif("not _HAS_H5PY")
+    def test_hdf_roundtrip_give_old__format(self):
+        so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        with pytest.warns(DeprecationWarning):
+            so.write("dummy.hdf5", format_="hdf5")
+        with pytest.warns(DeprecationWarning):
+            new_so = DummyStingrayObj.read("dummy.hdf5", format_="hdf5")
+        os.unlink("dummy.hdf5")
+
+        _check_equal(so, new_so)
+
+    def test_file_roundtrip_fits(self):
+        so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        so.panesapa = np.random.randint(5, 9, (6, 2))
+        so.write("dummy.fits")
+        new_so = DummyStingrayObj.read("dummy.fits")
+        os.unlink("dummy.fits")
+        # panesapa is invalid for FITS header and got lost
+        assert not hasattr(new_so, "panesapa")
+        new_so.panesapa = so.panesapa
+        _check_equal(so, new_so)
+
+    @pytest.mark.parametrize("fmt", ["pickle", "ascii", "ascii.ecsv"])
+    def test_file_roundtrip(self, fmt):
+        so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        so.panesapa = np.random.randint(5, 9, (6, 2))
+        so.write(f"dummy.{fmt}", fmt=fmt)
+        new_so = DummyStingrayObj.read(f"dummy.{fmt}", fmt=fmt)
+        os.unlink(f"dummy.{fmt}")
+
+        _check_equal(so, new_so)
+
+
+class TestStingrayTimeseries:
+    @classmethod
+    def setup_class(cls):
+        cls.arr = [4, 5, 2]
+        sting_obj = DummyStingrayTs(cls.arr)
+        sting_obj.time = np.asarray([0, 1, 2])
+        sting_obj.mjdref = 59777.000
+        sting_obj.parafritus = "bonus!"
+        sting_obj.panesapa = np.asarray([[41, 25], [98, 3]])
+        sting_obj.gti = np.asarray([[-0.5, 2.5]])
+        cls.sting_obj = sting_obj
+
+    def test_astropy_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        # Set an attribute to a DummyStingrayObj. It will *not* be saved
+        ts = so.to_astropy_table()
+        new_so = DummyStingrayTs.from_astropy_table(ts)
+        _check_equal(so, new_so)
+
+    def test_astropy_ts_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        ts = so.to_astropy_timeseries()
+        new_so = DummyStingrayTs.from_astropy_timeseries(ts)
+        _check_equal(so, new_so)
+
+    def test_shift_time(self):
+        new_so = self.sting_obj.shift(1)
+        assert np.allclose(new_so.time - 1, self.sting_obj.time)
+        assert np.allclose(new_so.gti - 1, self.sting_obj.gti)
+
+    def test_change_mjdref(self):
+        new_so = self.sting_obj.change_mjdref(59776.5)
+        assert new_so.mjdref == 59776.5
+        assert np.allclose(new_so.time - 43200, self.sting_obj.time)
+        assert np.allclose(new_so.gti - 43200, self.sting_obj.gti)
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 90c8f62a..4c86d510 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -18,6 +18,25 @@ from stingray.fourier import poisson_level, raw_coherence
 from stingray.events import EventList
 import copy
 
+_HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = True
+
+try:
+    import xarray
+    from xarray import Dataset
+except ImportError:
+    _HAS_XARRAY = False
+
+try:
+    import pandas
+    from pandas import DataFrame
+except ImportError:
+    _HAS_PANDAS = False
+
+try:
+    import h5py
+except ImportError:
+    _HAS_H5PY = False
+
 np.random.seed(20160528)
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
@@ -1141,6 +1160,25 @@ class TestAveragedCrossspectrum(object):
             assert hasattr(new_cs.pds1, attr) and getattr(new_cs.pds1, attr).size == N
             assert hasattr(new_cs.pds2, attr) and getattr(new_cs.pds2, attr).size == N
 
+        for attr in cs1.meta_attrs():
+            if attr not in ["df", "gti", "m"]:
+                assert getattr(cs1, attr) == getattr(new_cs, attr)
+
+    @pytest.mark.parametrize("norm", ["leahy", "abs", "frac", "none"])
+    def test_rebin_factor_log_rebins_all_attrs(self, norm):
+        cs1 = AveragedCrossspectrum.from_lightcurve(self.lc1, self.lc2, segment_size=1,
+                                                    norm=norm)
+        new_cs = cs1.rebin_log(0.03)
+        N = new_cs.freq.size
+        for attr in ["power", "power_err", "unnorm_power", "unnorm_power_err"]:
+            assert hasattr(new_cs, attr) and getattr(new_cs, attr).size == N
+            assert hasattr(new_cs.pds1, attr) and getattr(new_cs.pds1, attr).size == N
+            assert hasattr(new_cs.pds2, attr) and getattr(new_cs.pds2, attr).size == N
+
+        for attr in cs1.meta_attrs():
+            if attr not in ["df", "gti", "m"]:
+                assert getattr(cs1, attr) == getattr(new_cs, attr)
+
     def test_rebin(self):
         with warnings.catch_warnings(record=True) as w:
             new_cs = self.cs.rebin(df=1.5)
@@ -1280,3 +1318,58 @@ class TestTimelagFunction(object):
 
         assert np.max(lag) <= np.pi
         assert np.min(lag) >= -np.pi
+
+
+class TestRoundTrip():
+    @classmethod
+    def setup_class(cls):
+        cls.cs = Crossspectrum()
+        cls.cs.freq = np.arange(10)
+        cls.cs.power = np.random.uniform(0, 10, 10) + 3j
+        cls.cs.m = 1
+        cls.cs.nphots1 = 34
+        cls.cs.nphots2 = 25
+
+    def _check_equal(self, so, new_so):
+        for attr in ["freq", "power"]:
+            assert np.allclose(getattr(so, attr), getattr(new_so, attr))
+
+        for attr in ["m", "nphots1", "nphots2"]:
+            assert getattr(so, attr) == getattr(new_so, attr)
+
+    def test_astropy_roundtrip(self):
+        so = self.cs
+        ts = so.to_astropy_table()
+        new_so = so.from_astropy_table(ts)
+        self._check_equal(so, new_so)
+
+    @pytest.mark.skipif('not _HAS_XARRAY')
+    def test_xarray_roundtrip(self):
+        so = self.cs
+        ts = so.to_xarray()
+        new_so = so.from_xarray(ts)
+
+        self._check_equal(so, new_so)
+
+    @pytest.mark.skipif('not _HAS_PANDAS')
+    def test_pandas_roundtrip(self):
+        so = self.cs
+        ts = so.to_pandas()
+        new_so = so.from_pandas(ts)
+
+        self._check_equal(so, new_so)
+
+    @pytest.mark.parametrize("fmt", ["pickle", "ascii", "ascii.ecsv", "fits", "hdf5"])
+    def test_file_roundtrip(self, fmt):
+        so = self.cs
+        fname = f"dummy.{fmt}"
+        if not _HAS_H5PY and fmt == "hdf5":
+            with pytest.raises(Exception) as excinfo:
+                so.write(fname, fmt=fmt)
+                assert h5py in str(excinfo.value)
+            return True
+        so.write(fname, fmt=fmt)
+        new_so = so.read(fname, fmt=fmt)
+        os.unlink(fname)
+
+        self._check_equal(so, new_so)
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index 9f2d8ba6..50eee0f6 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -44,6 +44,7 @@ class TestEvents(object):
 
     @classmethod
     def setup_class(self):
+        np.random.seed(57239875)
         self.time = [0.5, 1.5, 2.5, 3.5]
         self.counts = [3000, 2000, 2200, 3600]
         self.counts_flat = [3000, 3000, 3000, 3000]
@@ -113,13 +114,14 @@ class TestEvents(object):
         lc_sim = ev.to_lc(dt=lc.dt, tstart=lc.tstart, tseg=lc.tseg)
         assert np.all((lc - lc_sim).counts < 3 * np.sqrt(lc.counts))
 
-    def test_simulate_times_with_spline(self):
+    @pytest.mark.parametrize("use_spline", [True, False])
+    def test_simulate_times(self, use_spline):
         """Simulate photon arrival times, with use_spline option
         enabled.
         """
         lc = Lightcurve(self.time, self.counts_flat, gti=self.gti)
         ev = EventList()
-        ev.simulate_times(lc, use_spline=True)
+        ev.simulate_times(lc, use_spline=use_spline)
         lc_sim = ev.to_lc(dt=lc.dt, tstart=lc.tstart, tseg=lc.tseg)
         assert np.all((lc - lc_sim).counts < 3 * np.sqrt(lc.counts))
 
@@ -317,35 +319,44 @@ class TestEvents(object):
                 np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
         assert np.allclose(ev_new.gti, np.array([[5, 6]]))
 
+    @pytest.mark.skipif('not (_HAS_YAML)')
+    def test_io_warns(self):
+        ev = EventList(self.time)
+        with pytest.warns(DeprecationWarning):
+            ev.write('ascii_ev.ecsv', format_='pickle')
+
+        with pytest.warns(DeprecationWarning):
+            ev = ev.read('ascii_ev.ecsv', format_='pickle')
+
     @pytest.mark.skipif('not (_HAS_YAML)')
     def test_io_with_ascii(self):
         ev = EventList(self.time)
-        ev.write('ascii_ev.ecsv', format_='ascii')
-        ev = ev.read('ascii_ev.ecsv', format_='ascii')
+        ev.write('ascii_ev.ecsv', fmt='ascii')
+        ev = ev.read('ascii_ev.ecsv', fmt='ascii')
         print(ev.time, self.time)
         assert np.allclose(ev.time, self.time)
         os.remove('ascii_ev.ecsv')
 
     def test_io_with_pickle(self):
         ev = EventList(self.time, mjdref=54000)
-        ev.write('ev.pickle', format_='pickle')
-        ev = ev.read('ev.pickle', format_='pickle')
+        ev.write('ev.pickle', fmt='pickle')
+        ev = ev.read('ev.pickle', fmt='pickle')
         assert np.allclose(ev.time, self.time)
         os.remove('ev.pickle')
 
     @pytest.mark.skipif("not _H5PY_INSTALLED")
     def test_io_with_hdf5(self):
         ev = EventList(time=self.time, mjdref=54000)
-        ev.write('ev.hdf5', format_='hdf5')
+        ev.write('ev.hdf5', fmt='hdf5')
 
-        ev = ev.read('ev.hdf5', format_='hdf5')
+        ev = ev.read('ev.hdf5', fmt='hdf5')
         assert np.allclose(ev.time, self.time)
         os.remove('ev.hdf5')
 
     def test_io_with_fits(self):
         ev = EventList(time=self.time, mjdref=54000)
-        ev.write('ev.fits', format_='fits')
-        ev = ev.read('ev.fits', format_='fits')
+        ev.write('ev.fits', fmt='fits')
+        ev = ev.read('ev.fits', fmt='fits')
         assert np.allclose(ev.time, self.time)
         os.remove('ev.fits')
 
@@ -355,7 +366,7 @@ class TestEvents(object):
         """
         fname = os.path.join(datadir, 'monol_testA.evt')
         ev = EventList()
-        ev = ev.read(fname, format_='hea')
+        ev = ev.read(fname, fmt='hea')
         assert np.isclose(ev.mjdref, 55197.00076601852)
 
     def test_fits_with_additional(self):
@@ -364,7 +375,7 @@ class TestEvents(object):
         """
         fname = os.path.join(datadir, 'xmm_test.fits')
         ev = EventList()
-        ev = ev.read(fname, format_='hea', additional_columns=['PRIOR'])
+        ev = ev.read(fname, fmt='hea', additional_columns=['PRIOR'])
         assert hasattr(ev, 'prior')
 
     def test_timeseries_empty_evts(self):
diff --git a/stingray/tests/test_io.py b/stingray/tests/test_io.py
index 37640ba2..fcca3b7a 100644
--- a/stingray/tests/test_io.py
+++ b/stingray/tests/test_io.py
@@ -5,11 +5,11 @@ import matplotlib.pyplot as plt
 
 from astropy.tests.helper import pytest
 
-from ..io import read, write, split_numbers
+from ..io import split_numbers
 from ..io import ref_mjd
 from ..io import high_precision_keyword_read
 from ..io import load_events_and_gtis, read_mission_info
-from ..io import read_header_key, _retrieve_ascii_object
+from ..io import read_header_key
 
 import warnings
 
@@ -163,132 +163,6 @@ class TmpIOReadWrite(object):
 
 class TestFileFormats(object):
 
-    def test_pickle_read_write(self):
-        test_object = TmpIOReadWrite()
-        write(test_object, filename='test.pickle', format_='pickle')
-        assert read('test.pickle', 'pickle') is not None
-        os.remove('test.pickle')
-
-    def test_pickle_attributes(self):
-        """Test if pickle maintains class object attributes."""
-        test_object = TmpIOReadWrite()
-        write(test_object, filename='test.pickle', format_='pickle')
-        rec_object = read('test.pickle', 'pickle')
-        assert rec_object.number == test_object.number
-        assert rec_object.str == test_object.str
-        assert rec_object.list == test_object.list
-        assert (rec_object.array == test_object.array).all()
-        assert rec_object.long_number == test_object.long_number
-        assert (rec_object.long_array == test_object.long_array).all()
-
-        os.remove('test.pickle')
-
-    def test_pickle_functions(self):
-        """Test if pickle maintains class methods."""
-        test_object = TmpIOReadWrite()
-        write(test_object,'test.pickle', 'pickle')
-        assert read('test.pickle', 'pickle').test_operation() == test_object.number * 10
-        os.remove('test.pickle')
-
-    @skip_condition
-    def test_hdf5_write(self):
-        test_object = TmpIOReadWrite()
-        write(test_object, 'test.hdf5', 'hdf5')
-        os.remove('test.hdf5')
-
-    @skip_condition
-    def test_hdf5_read(self):
-        test_object = TmpIOReadWrite()
-        write(test_object, 'test.hdf5', 'hdf5')
-        read('test.hdf5','hdf5')
-        os.remove('test.hdf5')
-
-    @skip_condition
-    def test_hdf5_data_recovery(self):
-        test_object = TmpIOReadWrite()
-        write(test_object, 'test.hdf5', 'hdf5')
-        rec_object = read('test.hdf5','hdf5')
-
-        assert rec_object['number'] == test_object.number
-        assert rec_object['str'] == test_object.str
-        assert (rec_object['list'] == test_object.list).all()
-        assert (rec_object['array'] == np.array(test_object.array)).all()
-        assert rec_object['long_number'] == np.double(test_object.long_number)
-        assert (rec_object['long_array'] == np.double(np.array(test_object.long_array))).all()
-
-        os.remove('test.hdf5')
-
-    def test_save_ascii(self):
-        time = [1, 2, 3, 4]
-        counts = [2, 3, 41, 4]
-
-        write(np.array([time, counts]).T, "ascii_test.txt",
-              "ascii")
-
-        os.remove("ascii_test.txt")
-
-    def test_save_ascii_with_mixed_types(self):
-        time = ["bla", 1, 2, 3]
-        counts = [2,3,41,4]
-        with pytest.raises(Exception):
-             write(np.array([time, counts]).T,
-                   "ascii_test.txt", "ascii")
-
-    def test_save_ascii_with_format(self):
-        time = ["bla", 1, 2, 3]
-        counts = [2,3,41,4]
-        write(np.array([time, counts]).T,
-              filename="ascii_test.txt", format_="ascii",
-              fmt=["%s", "%s"])
-
-    def test_retrieve_bad(self):
-        with pytest.raises(TypeError):
-            _retrieve_ascii_object(1)
-
-    def test_read_ascii(self):
-        time = [1,2,3,4,5]
-        counts = [5,7,8,2,3]
-        np.savetxt("ascii_test.txt", np.array([time, counts]).T)
-        read("ascii_test.txt", "ascii")
-        os.remove("ascii_test.txt")
-
-    def test_fits_write(self):
-        test_object = TmpIOReadWrite()
-        with warnings.catch_warnings(record=True) as w:
-            write(test_object, 'test.fits', 'fits')
-        os.remove('test.fits')
-
-    def test_fits_read(self):
-        test_object = TmpIOReadWrite()
-        with warnings.catch_warnings(record=True) as w:
-            write(test_object, 'test.fits', 'fits')
-            read('test.fits','fits',cols=['array','number','long_number'])
-        os.remove('test.fits')
-
-    def test_fits_with_multiple_tables(self):
-        test_object = TmpIOReadWrite()
-        with warnings.catch_warnings(record=True) as w:
-            write(test_object, 'test.fits', 'fits', tnames=['EVENTS', 'GTI'],
-                colsassign={'number':'GTI', 'array':'GTI'})
-        os.remove('test.fits')
-
-    def test_fits_data_recovery(self):
-        test_object = TmpIOReadWrite()
-        with warnings.catch_warnings(record=True) as w:
-            write(test_object, 'test.fits', 'fits')
-            rec_object = read('test.fits', 'fits', cols = ['number', 'str', 'list',
-                'array','long_array','long_number'])
-
-        assert rec_object['NUMBER'] == test_object.number
-        assert rec_object['STR'] == test_object.str
-        assert (rec_object['LIST'] == test_object.list).all()
-        assert (rec_object['ARRAY'] == np.array(test_object.array)).all()
-        assert rec_object['LONG_NUMBER'] == np.double(test_object.long_number)
-        assert (rec_object['LONG_ARRAY'] == np.double(np.array(test_object.long_array))).all()
-
-        del rec_object
-        os.remove('test.fits')
-
     def test_savefig_without_plot(self):
         from ..io import savefig
         plt.close("all")
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index cf7f723b..71b4e71c 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -1021,7 +1021,7 @@ class TestLightcurve(object):
         fname = 'lcurveA.fits'
         with pytest.warns(UserWarning):
             lc = Lightcurve.read(os.path.join(datadir, fname),
-                                 format_='hea', skip_checks=True)
+                                 fmt='hea', skip_checks=True)
         ctrate = 1
         assert np.isclose(lc.countrate[0], ctrate)
 
@@ -1029,25 +1029,32 @@ class TestLightcurve(object):
         fname = 'lcurve_new.fits'
         with pytest.warns(UserWarning):
             lc = Lightcurve.read(os.path.join(datadir, fname),
-                                 format_='hea', skip_checks=True)
+                                 fmt='hea', skip_checks=True)
         ctrate = 0.91
 
         assert np.isclose(lc.countrate[0], ctrate)
         assert np.isclose(lc.mjdref, 55197.00076601852)
 
+    def test_io_warns(self):
+        lc = Lightcurve(self.times, self.counts)
+        with pytest.warns(DeprecationWarning):
+            lc.write('lc.pickle', format_='pickle')
+        with pytest.warns(DeprecationWarning):
+            lc.read('lc.pickle', format_='pickle')
+
     @pytest.mark.skipif('not _HAS_YAML')
     def test_io_with_ascii(self):
         lc = Lightcurve(self.times, self.counts)
-        lc.write('ascii_lc.ecsv', format_='ascii')
-        lc = lc.read('ascii_lc.ecsv', format_='ascii')
+        lc.write('ascii_lc.ecsv', fmt='ascii')
+        lc = lc.read('ascii_lc.ecsv', fmt='ascii')
         assert np.allclose(lc.time, self.times)
         assert np.allclose(lc.counts, self.counts)
         os.remove('ascii_lc.ecsv')
 
     def test_io_with_pickle(self):
         lc = Lightcurve(self.times, self.counts)
-        lc.write('lc.pickle', format_='pickle')
-        lc.read('lc.pickle', format_='pickle')
+        lc.write('lc.pickle', fmt='pickle')
+        lc.read('lc.pickle', fmt='pickle')
         assert np.allclose(lc.time, self.times)
         assert np.allclose(lc.counts, self.counts)
         assert np.allclose(lc.gti, self.gti)
@@ -1056,9 +1063,9 @@ class TestLightcurve(object):
     @pytest.mark.skipif('not _H5PY_INSTALLED')
     def test_io_with_hdf5(self):
         lc = Lightcurve(self.times, self.counts)
-        lc.write('lc.hdf5', format_='hdf5')
+        lc.write('lc.hdf5', fmt='hdf5')
 
-        data = lc.read('lc.hdf5', format_='hdf5')
+        data = lc.read('lc.hdf5', fmt='hdf5')
         assert np.allclose(data.time, self.times)
         assert np.allclose(data.counts, self.counts)
         assert np.allclose(data.gti, self.gti)
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index 3a383b83..2a8754e6 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -10,6 +10,25 @@ from stingray.events import EventList
 from stingray import Powerspectrum, AveragedPowerspectrum, \
     DynamicalPowerspectrum
 
+_HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = True
+
+try:
+    import xarray
+    from xarray import Dataset
+except ImportError:
+    _HAS_XARRAY = False
+
+try:
+    import pandas
+    from pandas import DataFrame
+except ImportError:
+    _HAS_PANDAS = False
+
+try:
+    import h5py
+except ImportError:
+    _HAS_H5PY = False
+
 np.random.seed(20150907)
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
@@ -69,7 +88,7 @@ class TestAveragedPowerspectrumEvents(object):
     def test_type_change(self):
         pds = copy.deepcopy(self.leahy_pds)
         assert pds.type == "powerspectrum"
-        pds._type = "astdfawerfsaf"
+        pds.type = "astdfawerfsaf"
         assert pds.type == "astdfawerfsaf"
 
     def test_from_events_works_ps(self):
@@ -246,14 +265,6 @@ class TestAveragedPowerspectrumEvents(object):
                                     norm="Leahy", dt=self.dt)
         bin_aps = aps.rebin_log(df)
 
-    def test_rebin_with_invalid_type_attribute(self):
-        new_df = 2
-        aps = AveragedPowerspectrum(self.lc, segment_size=1,
-                                    norm='leahy', dt=self.dt)
-        aps.type = 'invalid_type'
-        with pytest.raises(AttributeError):
-            assert aps.rebin(df=new_df)
-
     @pytest.mark.parametrize("use_common_mean", [True, False])
     @pytest.mark.parametrize("legacy", [True, False])
     def test_leahy_correct_for_multiple(self, legacy, use_common_mean):
@@ -837,14 +848,6 @@ class TestAveragedPowerspectrum(object):
                                     norm="Leahy")
         bin_aps = aps.rebin_log(df)
 
-    def test_rebin_with_invalid_type_attribute(self):
-        new_df = 2
-        aps = AveragedPowerspectrum(self.lc, segment_size=1,
-                                    norm='leahy')
-        aps.type = 'invalid_type'
-        with pytest.raises(AttributeError):
-            assert aps.rebin(df=new_df)
-
     @pytest.mark.parametrize("legacy", [True, False])
     def test_list_with_nonsense_component(self, legacy):
         n_lcs = 10
@@ -1080,3 +1083,61 @@ class TestDynamicalPowerspectrum(object):
         assert np.allclose(new_dps.freq, rebin_freq)
         assert np.allclose(new_dps.dyn_ps, rebin_dps, atol=0.00001)
         assert np.isclose(new_dps.df, df_new)
+
+
+class TestRoundTrip():
+    @classmethod
+    def setup_class(cls):
+        cls.cs = AveragedPowerspectrum()
+        cls.cs.freq = np.arange(10)
+        cls.cs.power = np.random.uniform(0, 10, 10)
+        cls.cs.m = 2
+        cls.cs.nphots1 = 34
+
+    def _check_equal(self, so, new_so):
+        for attr in ["freq", "power"]:
+            assert np.allclose(getattr(so, attr), getattr(new_so, attr))
+
+        for attr in ["m", "nphots1"]:
+            assert getattr(so, attr) == getattr(new_so, attr)
+
+    def test_astropy_roundtrip(self):
+        so = self.cs
+        ts = so.to_astropy_table()
+        new_so = so.from_astropy_table(ts)
+        self._check_equal(so, new_so)
+
+    @pytest.mark.skipif('not _HAS_XARRAY')
+    def test_xarray_roundtrip(self):
+        so = self.cs
+        ts = so.to_xarray()
+        new_so = so.from_xarray(ts)
+
+        self._check_equal(so, new_so)
+
+    @pytest.mark.skipif('not _HAS_PANDAS')
+    def test_pandas_roundtrip(self):
+        so = self.cs
+        ts = so.to_pandas()
+        new_so = so.from_pandas(ts)
+
+        self._check_equal(so, new_so)
+
+    @pytest.mark.skipif('not _HAS_H5PY')
+    def test_hdf_roundtrip(self):
+        so = self.cs
+        so.write("dummy.hdf5")
+        new_so = so.read("dummy.hdf5")
+        os.unlink("dummy.hdf5")
+
+        self._check_equal(so, new_so)
+
+    @pytest.mark.parametrize("fmt", ["pickle", "ascii", "ascii.ecsv", "fits"])
+    def test_file_roundtrip(self, fmt):
+        so = self.cs
+        fname = f"dummy.{fmt}"
+        so.write(fname, fmt=fmt)
+        new_so = so.read(fname, fmt=fmt)
+        # os.unlink(fname)
+
+        self._check_equal(so, new_so)
diff --git a/stingray/tests/test_varenergyspectrum.py b/stingray/tests/test_varenergyspectrum.py
index c30c4550..a0a78f21 100644
--- a/stingray/tests/test_varenergyspectrum.py
+++ b/stingray/tests/test_varenergyspectrum.py
@@ -10,6 +10,26 @@ from stingray.varenergyspectrum import ExcessVarianceSpectrum
 from stingray.lightcurve import Lightcurve
 
 from astropy.tests.helper import pytest
+from astropy.table import Table
+
+_HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = True
+
+try:
+    import xarray
+    from xarray import Dataset
+except ImportError:
+    _HAS_XARRAY = False
+
+try:
+    import pandas
+    from pandas import DataFrame
+except ImportError:
+    _HAS_PANDAS = False
+
+try:
+    import h5py
+except ImportError:
+    _HAS_H5PY = False
 
 np.random.seed(20150907)
 curdir = os.path.abspath(os.path.dirname(__file__))
@@ -398,3 +418,79 @@ class TestLagEnergySpectrum(object):
 
         assert np.all(np.isnan(lag.spectrum))
         assert np.all(np.isnan(lag.spectrum_error))
+
+
+class TestRoundTrip():
+    @classmethod
+    def setup_class(cls):
+        tstart = 0.0
+        tend = 100.0
+        nphot = 1000
+        alltimes = np.random.uniform(tstart, tend, nphot)
+        alltimes.sort()
+        cls.events = EventList(
+            alltimes, energy=np.random.uniform(0.3, 12, nphot), gti=[[tstart, tend]]
+        )
+        cls.vespec = DummyVarEnergy(
+            cls.events, [0.0, 10000], (0.5, 5, 10, "lin"), [0.3, 10], bin_time=0.1
+        )
+        cls.vespec.spectrum = np.zeros_like(cls.vespec.energy)
+        cls.vespec.spectrum_error = np.zeros_like(cls.vespec.energy)
+
+    def _check_equal(self, so, table):
+        for attr in ["energy", "spectrum", "spectrum_error"]:
+            assert np.allclose(getattr(so, attr), table[attr])
+
+        if hasattr(table, "meta"):
+            for attr in ["freq_interval"]:
+                assert getattr(so, attr) == table.meta[attr]
+        if hasattr(table, "attrs"):
+            for attr in ["freq_interval"]:
+                assert getattr(so, attr) == table.attrs[attr]
+
+    def test_astropy_export(self):
+        so = self.vespec
+        ts = so.to_astropy_table()
+        self._check_equal(so, ts)
+        with pytest.raises(NotImplementedError):
+            so.from_astropy_table(ts)
+
+    @pytest.mark.skipif('not _HAS_XARRAY')
+    def test_xarray_export(self):
+        so = self.vespec
+        ts = so.to_xarray()
+        self._check_equal(so, ts)
+        with pytest.raises(NotImplementedError):
+            so.from_xarray(ts)
+
+    @pytest.mark.skipif('not _HAS_PANDAS')
+    def test_pandas_export(self):
+        so = self.vespec
+        ts = so.to_pandas()
+        self._check_equal(so, ts)
+        with pytest.raises(NotImplementedError):
+            so.from_pandas(ts)
+
+    @pytest.mark.skipif('not _HAS_H5PY')
+    def test_hdf_export(self):
+        so = self.vespec
+        so.write("dummy.hdf5")
+        new_so = Table.read("dummy.hdf5")
+        os.unlink("dummy.hdf5")
+        self._check_equal(so, new_so)
+
+    @pytest.mark.parametrize("fmt", ["ascii.ecsv", "fits"])
+    def test_file_export(self, fmt):
+        so = self.vespec
+        so.write("dummy", fmt=fmt)
+        new_so = Table.read("dummy", format=fmt)
+        os.unlink("dummy")
+        self._check_equal(so, new_so)
+
+    @pytest.mark.parametrize("fmt", ["pickle"])
+    def test_file_export_pickle(self, fmt):
+        so = self.vespec
+        so.write("dummy", fmt=fmt)
+        new_so = so.read("dummy", fmt=fmt)
+        os.unlink("dummy")
+        self._check_equal(so, new_so.to_astropy_table())