Repository: StingraySoftware/stingray

Issue:
Check from_lc_iterable and from_lightcurve for a single light curve

Discussion:
pep8speaks: Hello @matteobachetti! Thanks for updating 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/crossspectrum.py`](https://github.com/StingraySoftware/stingray/blob/60295b194d5022d555ae0c3edb48b4488df2fef1/stingray/crossspectrum.py):

> [Line 756:13](https://github.com/StingraySoftware/stingray/blob/60295b194d5022d555ae0c3edb48b4488df2fef1/stingray/crossspectrum.py#L756): [E122](https://duckduckgo.com/?q=pep8%20E122) continuation line missing indentation or outdented










codecov[bot]: # [Codecov](https://codecov.io/gh/StingraySoftware/stingray/pull/624?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) Report
> Merging [#624](https://codecov.io/gh/StingraySoftware/stingray/pull/624?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) (4cdddb0) into [main](https://codecov.io/gh/StingraySoftware/stingray/commit/e736acd21b007cee35b30bec4cf62b78ae8539e4?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) (e736acd) will **decrease** coverage by `3.52%`.
> The diff coverage is `95.12%`.

[![Impacted file tree graph](https://codecov.io/gh/StingraySoftware/stingray/pull/624/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/624?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware)

```diff
@@            Coverage Diff             @@
##             main     #624      +/-   ##
==========================================
- Coverage   96.67%   93.15%   -3.53%     
==========================================
  Files          40       40              
  Lines        7670     7713      +43     
==========================================
- Hits         7415     7185     -230     
- Misses        255      528     +273     
```


| [Impacted Files](https://codecov.io/gh/StingraySoftware/stingray/pull/624?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) | Coverage Δ | |
|---|---|---|
| [stingray/fourier.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvZm91cmllci5weQ==) | `98.36% <ø> (-1.37%)` | :arrow_down: |
| [stingray/powerspectrum.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvcG93ZXJzcGVjdHJ1bS5weQ==) | `97.26% <93.47%> (-2.11%)` | :arrow_down: |
| [stingray/crossspectrum.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvY3Jvc3NzcGVjdHJ1bS5weQ==) | `97.59% <96.10%> (-1.15%)` | :arrow_down: |
| [stingray/deadtime/model.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvZGVhZHRpbWUvbW9kZWwucHk=) | `64.64% <0.00%> (-33.34%)` | :arrow_down: |
| [stingray/pulse/\_\_init\_\_.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvcHVsc2UvX19pbml0X18ucHk=) | `76.92% <0.00%> (-23.08%)` | :arrow_down: |
| [stingray/utils.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvdXRpbHMucHk=) | `79.69% <0.00%> (-19.35%)` | :arrow_down: |
| [stingray/stats.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvc3RhdHMucHk=) | `81.48% <0.00%> (-17.99%)` | :arrow_down: |
| [stingray/pulse/pulsar.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvcHVsc2UvcHVsc2FyLnB5) | `91.81% <0.00%> (-6.77%)` | :arrow_down: |
| [stingray/pulse/search.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvcHVsc2Uvc2VhcmNoLnB5) | `89.93% <0.00%> (-6.72%)` | :arrow_down: |
| [stingray/filters.py](https://codecov.io/gh/StingraySoftware/stingray/pull/624/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvZmlsdGVycy5weQ==) | `94.44% <0.00%> (-4.63%)` | :arrow_down: |
| ... and [15 more](https://codecov.io/gh/StingraySoftware/stingray/pull/624/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/624?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/624?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware). Last update [e736acd...4cdddb0](https://codecov.io/gh/StingraySoftware/stingray/pull/624?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/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 23458c14..26b5c9f1 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -497,10 +497,9 @@ class Crossspectrum(object):
 
     Other Parameters
     ----------------
-    gti: 2-d float array
-        ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
-        This choice overrides the GTIs in the single light curves. Use with
-        care!
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
 
     lc1: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects
         For backwards compatibility only. Like ``data1``, but no
@@ -550,6 +549,11 @@ class Crossspectrum(object):
 
     nphots2: float
         The total number of photons in light curve 2
+
+    legacy: bool
+        Use the legacy machinery of AveragedCrossspectrum. This might be useful to compare
+        with old results, and is also needed to use light curve lists as an input, to
+        conserve the spectra of each segment, or to use the large_data option.
     """
 
     def __init__(
@@ -563,7 +567,8 @@ class Crossspectrum(object):
         power_type="real",
         dt=None,
         fullspec=False,
-        skip_checks=False
+        skip_checks=False,
+        legacy=False
     ):
         self._type = None
         # for backwards compatibility
@@ -586,21 +591,17 @@ class Crossspectrum(object):
                 fullspec=fullspec,
             )
 
-        norm = norm.lower()
         self.dt = dt
+        norm = norm.lower()
+        self.norm = norm
 
         if not good_input:
-            self.freq = None
-            self.power = None
-            self.power_err = None
-            self.df = None
-            self.dt = None
-            self.nphots1 = None
-            self.nphots2 = None
-            self.m = 1
-            self.n = None
-            self.norm = norm
-            return
+            return self._initialize_empty()
+
+        if not legacy and data1 is not None and data2 is not None:
+            return self._initialize_from_any_input(
+                data1, data2, dt=dt, norm=norm, power_type=power_type,
+                fullspec=fullspec, gti=gti)
 
         if not isinstance(data1, EventList):
             lc1 = data1
@@ -634,6 +635,7 @@ class Crossspectrum(object):
         gti=None,
         lc1=None,
         lc2=None,
+        segment_size=None,
         power_type="real",
         dt=None,
         fullspec=False,
@@ -712,11 +714,51 @@ class Crossspectrum(object):
 
         dt_is_invalid = (dt is None) or (dt <= np.finfo(float).resolution)
 
-        if (isinstance(data1, EventList) or isinstance(data2, EventList)) and dt_is_invalid:
-            raise ValueError(
-                "If using event lists, please specify the bin time to generate lightcurves."
+        if segment_size is None:
+            # checks to be run for non-averaged spectra
+            if gti is not None and len(gti) > 1:
+                raise TypeError("Non-averaged cross spectra need a single GTI")
+
+        if type(data1) != type(data2):
+            raise TypeError("Input data have to be of the same kind")
+
+        if isinstance(data1, EventList):
+            if dt_is_invalid:
+                raise ValueError(
+                    "If using event lists, please specify the bin time to generate lightcurves."
+                )
+        elif isinstance(data1, Lightcurve):
+            if (data1.err_dist.lower() != data2.err_dist.lower()):
+                simon(
+                    "Your lightcurves have different statistics."
+                    "The errors in the Crossspectrum will be incorrect."
             )
 
+            # If dt differs slightly, its propagated error must not be more than
+            # 1/100th of the bin
+            if not np.isclose(data1.dt, data2.dt, rtol=0.01 * data1.dt / data1.tseg):
+                raise StingrayError("Light curves do not have same time binning dt.")
+
+            if data1.tseg != data2.tseg:
+                simon(
+                    "Lightcurves do not have same tseg. This means that the data"
+                    "from the two channels are not completely in sync. This "
+                    "might or might not be an issue. Keep an eye on it."
+                )
+        elif isinstance(data1, (list, tuple)):
+            if not isinstance(data1[0], Lightcurve) or not isinstance(data2[0], Lightcurve):
+                raise TypeError("Inputs lists have to contain light curve objects")
+
+            if (data1[0].err_dist.lower() != data2[0].err_dist.lower()):
+                simon(
+                    "Your lightcurves have different statistics."
+                    "The errors in the Crossspectrum will be incorrect."
+                )
+        elif isinstance(data1, (Generator, Iterator)):
+            pass
+        else:
+            raise TypeError("Input data are invalid")
+
         return True
 
     @property
@@ -947,7 +989,7 @@ class Crossspectrum(object):
         bin_cs.nphots1 = self.nphots1
         bin_cs.power_err = binerr
 
-        if hasattr(self, "unnorm_power"):
+        if hasattr(self, "unnorm_power") and self.unnorm_power is not None:
             unnorm_power_err = None
             if hasattr(self, "unnorm_power_err") and self.unnorm_power_err is not None:
                 unnorm_power_err = self.unnorm_power_err
@@ -1148,7 +1190,7 @@ class Crossspectrum(object):
         new_spec.m = nsamples * self.m
         new_spec.dt = self.dt
 
-        if hasattr(self, "unnorm_power"):
+        if hasattr(self, "unnorm_power") and self.unnorm_power is not None:
             unnorm_power_err = None
             if hasattr(self, "unnorm_power_err") and self.unnorm_power_err is not None:
                 unnorm_power_err = self.unnorm_power_err
@@ -1427,7 +1469,8 @@ class Crossspectrum(object):
             power_type="all",
             silent=False,
             fullspec=False,
-            use_common_mean=True):
+            use_common_mean=True,
+            gti=None):
         """Calculate AveragedCrossspectrum from two event lists
 
         Parameters
@@ -1462,6 +1505,9 @@ class Crossspectrum(object):
         power_type : str, default 'all'
             If 'all', give complex powers. If 'abs', the absolute value; if 'real',
             the real part
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Good Time intervals. Defaults to the common GTIs from the two input
+            objects
         """
 
         return crossspectrum_from_events(
@@ -1473,7 +1519,8 @@ class Crossspectrum(object):
             power_type=power_type,
             silent=silent,
             fullspec=fullspec,
-            use_common_mean=use_common_mean)
+            use_common_mean=use_common_mean,
+            gti=gti)
 
     @staticmethod
     def from_lightcurve(
@@ -1484,7 +1531,8 @@ class Crossspectrum(object):
             power_type="all",
             silent=False,
             fullspec=False,
-            use_common_mean=True):
+            use_common_mean=True,
+            gti=None):
         """Calculate AveragedCrossspectrum from two light curves
 
         Parameters
@@ -1516,6 +1564,9 @@ class Crossspectrum(object):
         power_type : str, default 'all'
             If 'all', give complex powers. If 'abs', the absolute value; if 'real',
             the real part
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Good Time intervals. Defaults to the common GTIs from the two input
+            objects
         """
         return crossspectrum_from_lightcurve(
             lc1,
@@ -1525,7 +1576,8 @@ class Crossspectrum(object):
             power_type=power_type,
             silent=silent,
             fullspec=fullspec,
-            use_common_mean=use_common_mean)
+            use_common_mean=use_common_mean,
+            gti=gti)
 
     @staticmethod
     def from_lc_iterable(
@@ -1537,7 +1589,8 @@ class Crossspectrum(object):
             power_type="all",
             silent=False,
             fullspec=False,
-            use_common_mean=True):
+            use_common_mean=True,
+            gti=None):
         """Calculate AveragedCrossspectrum from two light curves
 
         Parameters
@@ -1572,6 +1625,9 @@ class Crossspectrum(object):
         power_type : str, default 'all'
             If 'all', give complex powers. If 'abs', the absolute value; if 'real',
             the real part
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Good Time intervals. Defaults to the common GTIs from the two input
+            objects
         """
 
         return crossspectrum_from_lc_iterable(
@@ -1583,8 +1639,84 @@ class Crossspectrum(object):
             power_type=power_type,
             silent=silent,
             fullspec=fullspec,
-            use_common_mean=use_common_mean)
+            use_common_mean=use_common_mean,
+            gti=gti)
+
+    def _initialize_from_any_input(
+            self, data1, data2, dt=None, segment_size=None, norm="frac",
+            power_type="all", silent=False, fullspec=False, gti=None,
+            use_common_mean=True):
+        """Initialize the class, trying to understand the input types.
+
+        The input arguments are the same as ``__init__()``. Based on the type
+        of ``data1``, this method will call the appropriate
+        ``crossspectrum_from_XXXX`` function, and initialize ``self`` with
+        the correct attributes.
+        """
+        if isinstance(data1, EventList):
+            spec = crossspectrum_from_events(
+                data1,
+                data2,
+                dt,
+                segment_size,
+                norm=norm,
+                power_type=power_type,
+                silent=silent,
+                fullspec=fullspec,
+                use_common_mean=use_common_mean,
+                gti=gti,
+            )
+        elif isinstance(data1, Lightcurve):
+            spec = crossspectrum_from_lightcurve(
+                data1,
+                data2,
+                segment_size,
+                norm=norm,
+                power_type=power_type,
+                silent=silent,
+                fullspec=fullspec,
+                use_common_mean=use_common_mean,
+                gti=gti,
+            )
+            spec.lc1 = data1
+            spec.lc2 = data2
+        elif isinstance(data1, (tuple, list)):
+            dt = data1[0].dt
+            # This is a list of light curves.
+            spec = crossspectrum_from_lc_iterable(
+                data1,
+                data2,
+                dt,
+                segment_size,
+                norm=norm,
+                power_type=power_type,
+                silent=silent,
+                fullspec=fullspec,
+                gti=gti,
+                use_common_mean=use_common_mean,
+            )
+        else:  # pragma: no cover
+            raise TypeError(f"Bad inputs to Crosssspectrum: {type(data1)}")
 
+        for key, val in spec.__dict__.items():
+            setattr(self, key, val)
+        return
+
+    def _initialize_empty(self):
+        """Set all attributes to None."""
+        self.freq = None
+        self.power = None
+        self.power_err = None
+        self.unnorm_power = None
+        self.unnorm_power_err = None
+        self.df = None
+        self.dt = None
+        self.nphots1 = None
+        self.nphots2 = None
+        self.m = 1
+        self.n = None
+        self.fullspec = None
+        return
 
 class AveragedCrossspectrum(Crossspectrum):
     """
@@ -1614,10 +1746,9 @@ class AveragedCrossspectrum(Crossspectrum):
 
     Other Parameters
     ----------------
-    gti: 2-d float array
-        ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
-        This choice overrides the GTIs in the single light curves. Use with
-        care!
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
 
     dt : float
         The time resolution of the light curve. Only needed when constructing
@@ -1667,6 +1798,10 @@ class AveragedCrossspectrum(Crossspectrum):
         with old results, and is also needed to use light curve lists as an input, to
         conserve the spectra of each segment, or to use the large_data option.
 
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
+
     Attributes
     ----------
     freq: numpy.ndarray
@@ -1697,10 +1832,9 @@ class AveragedCrossspectrum(Crossspectrum):
     nphots2: float
         The total number of photons in the second (reference) light curve
 
-    gti: 2-d float array
-        ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
-        They are calculated by taking the common GTI between the
-        two light curves
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
     """
 
     def __init__(
@@ -1744,23 +1878,15 @@ class AveragedCrossspectrum(Crossspectrum):
                 fullspec=fullspec,
                 segment_size=segment_size
             )
-        self.norm = norm.lower()
+        norm = norm.lower()
+        self.norm = norm
         self.dt = dt
         self.save_all = save_all
         self.segment_size = segment_size
         self.show_progress = not silent
 
         if not good_input:
-            self.freq = None
-            self.power = None
-            self.power_err = None
-            self.df = None
-            self.dt = None
-            self.nphots1 = None
-            self.nphots2 = None
-            self.m = 1
-            self.n = None
-            return
+            return self._initialize_empty()
 
         if isinstance(data1, Generator):
             warnings.warn(
@@ -1780,49 +1906,10 @@ class AveragedCrossspectrum(Crossspectrum):
             legacy = True
 
         if not legacy and data1 is not None and data2 is not None:
-            if isinstance(data1, EventList):
-                spec = crossspectrum_from_events(
-                    data1,
-                    data2,
-                    dt,
-                    segment_size,
-                    norm=norm,
-                    power_type=power_type,
-                    silent=silent,
-                    fullspec=fullspec,
-                    use_common_mean=use_common_mean,
-                )
-            elif isinstance(data1, Lightcurve):
-                spec = crossspectrum_from_lightcurve(
-                    data1,
-                    data2,
-                    segment_size,
-                    norm=norm,
-                    power_type=power_type,
-                    silent=silent,
-                    fullspec=fullspec,
-                    use_common_mean=use_common_mean,
-                )
-                spec.lc1 = data1
-                spec.lc2 = data2
-            else:
-                dt = data1[0].dt
-                # This is a list of light curves.
-                spec = crossspectrum_from_lc_iterable(
-                    data1,
-                    data2,
-                    dt,
-                    segment_size,
-                    norm=norm,
-                    power_type=power_type,
-                    silent=silent,
-                    fullspec=fullspec,
-                    use_common_mean=use_common_mean,
-                )
-
-            for key, val in spec.__dict__.items():
-                setattr(self, key, val)
-            return
+            return self._initialize_from_any_input(
+                data1, data2, dt=dt, segment_size=segment_size, gti=gti, norm=norm,
+                power_type=power_type, silent=silent, fullspec=fullspec,
+                use_common_mean=use_common_mean)
 
         log.info("Using legacy interface.")
 
@@ -1886,7 +1973,7 @@ class AveragedCrossspectrum(Crossspectrum):
 
         Crossspectrum.__init__(
             self, data1, data2, norm, gti=gti, power_type=power_type, dt=dt,
-            fullspec=fullspec, skip_checks=True
+            fullspec=fullspec, skip_checks=True, legacy=legacy
         )
 
         return
@@ -1913,7 +2000,7 @@ class AveragedCrossspectrum(Crossspectrum):
         ...
         ValueError: segment_size must be finite!
         """
-        good = Crossspectrum.initial_checks(self, data1, **kwargs)
+        good = Crossspectrum.initial_checks(self, data1, segment_size=segment_size, **kwargs)
         if not good:
             return False
         if isinstance(self, AveragedCrossspectrum) and segment_size is None and data1 is not None:
@@ -2000,13 +2087,6 @@ class AveragedCrossspectrum(Crossspectrum):
         assert isinstance(lc1, Lightcurve)
         assert isinstance(lc2, Lightcurve)
 
-        if lc1.tseg != lc2.tseg:
-            simon(
-                "Lightcurves do not have same tseg. This means that the data"
-                "from the two channels are not completely in sync. This "
-                "might or might not be an issue. Keep an eye on it."
-            )
-
         # If dt differs slightly, its propagated error must not be more than
         # 1/100th of the bin
         if not np.isclose(lc1.dt, lc2.dt, rtol=0.01 * lc1.dt / lc1.tseg):
@@ -2084,6 +2164,7 @@ class AveragedCrossspectrum(Crossspectrum):
                     norm=self.norm,
                     power_type=self.power_type,
                     fullspec=self.fullspec,
+                    legacy=True
                 )
 
             cs_all.append(cs_seg)
@@ -2299,8 +2380,9 @@ def crossspectrum_from_time_array(
     ----------------
     segment_size : float
         The length, in seconds, of the light curve segments that will be averaged
-    gti : [[gti0, gti1], ...]
-        Common Good Time intervals
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
     norm : str, default "frac"
         The normalization of the periodogram. "abs" is absolute rms, "frac" is
         fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
@@ -2354,6 +2436,7 @@ def crossspectrum_from_events(
     silent=False,
     fullspec=False,
     use_common_mean=True,
+    gti=None,
 ):
     """Calculate AveragedCrossspectrum from two event lists
 
@@ -2388,6 +2471,9 @@ def crossspectrum_from_events(
     power_type : str, default 'all'
         If 'all', give complex powers. If 'abs', the absolute value; if 'real',
         the real part
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
 
     Returns
     -------
@@ -2395,7 +2481,8 @@ def crossspectrum_from_events(
         The output cross spectrum.
     """
 
-    gti = cross_two_gtis(events1.gti, events2.gti)
+    if gti is None:
+        gti = cross_two_gtis(events1.gti, events2.gti)
 
     return crossspectrum_from_time_array(
         events1.time,
@@ -2420,6 +2507,7 @@ def crossspectrum_from_lightcurve(
     silent=False,
     fullspec=False,
     use_common_mean=True,
+    gti=None
 ):
     """Calculate AveragedCrossspectrum from two light curves
 
@@ -2451,6 +2539,9 @@ def crossspectrum_from_lightcurve(
     power_type : str, default 'all'
         If 'all', give complex powers. If 'abs', the absolute value; if 'real',
         the real part
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
 
     Returns
     -------
@@ -2460,10 +2551,13 @@ def crossspectrum_from_lightcurve(
     force_averaged = segment_size is not None
     # Suppress progress bar for single periodogram
     silent = silent or (segment_size is None)
-    gti = cross_two_gtis(lc1.gti, lc2.gti)
+    if gti is None:
+        gti = cross_two_gtis(lc1.gti, lc2.gti)
 
-    err1 = lc1._counts_err
-    err2 = lc2._counts_err
+    err1 = err2 = None
+    if lc1.err_dist == "gauss":
+        err1 = lc1._counts_err
+        err2 = lc2._counts_err
 
     results = avg_cs_from_events(
         lc1.time,
@@ -2496,6 +2590,7 @@ def crossspectrum_from_lc_iterable(
     silent=False,
     fullspec=False,
     use_common_mean=True,
+    gti=None,
 ):
     """Calculate AveragedCrossspectrum from two light curves
 
@@ -2530,6 +2625,9 @@ def crossspectrum_from_lc_iterable(
     power_type : str, default 'all'
         If 'all', give complex powers. If 'abs', the absolute value; if 'real',
         the real part
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. The GTIs of the input light curves are
+        interesected with these.
 
     Returns
     -------
@@ -2541,13 +2639,23 @@ def crossspectrum_from_lc_iterable(
     # Suppress progress bar for single periodogram
     silent = silent or (segment_size is None)
 
+    common_gti = gti
+
     def iterate_lc_counts(iter_lc):
         for lc in iter_lc:
             if hasattr(lc, "counts"):
                 n_bin = np.rint(segment_size / lc.dt).astype(int)
 
+                gti = lc.gti
+                if common_gti is not None:
+                    gti = cross_two_gtis(common_gti, lc.gti)
+
+                err = None
+                if lc.err_dist == "gauss":
+                    err = lc.counts_err
+
                 flux_iterable = get_flux_iterable_from_segments(
-                    lc.time, lc.gti, segment_size, n_bin, fluxes=lc.counts, errors=lc._counts_err
+                    lc.time, gti, segment_size, n_bin, fluxes=lc.counts, errors=err
                 )
                 for out in flux_iterable:
                     yield out
diff --git a/stingray/fourier.py b/stingray/fourier.py
index 5b1dda05..544b9fdd 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -1187,6 +1187,7 @@ def avg_cs_from_iterables_quick(
                          "mean1": common_mean1,
                          "mean2": common_mean2,
                          "power_type": "all",
+                         "fullspec": False,
                          "segment_size": dt * n_bin})
 
     return results
@@ -1454,6 +1455,7 @@ def avg_cs_from_iterables(
                          "mean1": common_mean1,
                          "mean2": common_mean2,
                          "power_type": power_type,
+                         "fullspec": fullspec,
                          "variance": common_variance,
                          "variance1": common_variance1,
                          "variance2": common_variance2})
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 237e7c94..18a72819 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -92,10 +92,14 @@ class Powerspectrum(Crossspectrum):
     nphots: float
         The total number of photons in the light curve
 
+    legacy: bool
+        Use the legacy machinery of AveragedPowerspectrum. This might be useful to compare
+        with old results, and is also needed to use light curve lists as an input, to
+        conserve the spectra of each segment, or to use the large_data option.
     """
 
     def __init__(self, data=None, norm="frac", gti=None,
-                 dt=None, lc=None, skip_checks=False):
+                 dt=None, lc=None, skip_checks=False, legacy=False):
 
         self._type = None
         if lc is not None:
@@ -116,23 +120,18 @@ class Powerspectrum(Crossspectrum):
                 dt=dt
             )
 
-        self.norm = norm.lower()
+        norm = norm.lower()
+        self.norm = norm
         self.dt = dt
 
         if not good_input:
-            self.freq = None
-            self.power = None
-            self.power_err = None
-            self.df = None
-            self.dt = None
-            self.nphots = None
-            self.nphots1 = None
-            self.m = 1
-            self.n = None
-            return
+            return self._initialize_empty()
+
+        if not legacy and data is not None:
+            return self._initialize_from_any_input(data, dt=dt, norm=norm)
 
         Crossspectrum.__init__(self, data1=data, data2=data, norm=norm, gti=gti,
-                               dt=dt, skip_checks=True)
+                               dt=dt, skip_checks=True, legacy=legacy)
         self.nphots = self.nphots1
         self.dt = dt
 
@@ -411,13 +410,14 @@ class Powerspectrum(Crossspectrum):
         dt : float
             The time resolution of the intermediate light curves
             (sets the Nyquist frequency)
+
         Other parameters
         ----------------
         segment_size : float
             The length, in seconds, of the light curve segments that will be averaged
             Only relevant (and required) for AveragedPowerspectrum
-        gti : [[gti0, gti1], ...]
-            Good Time intervals
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Good Time intervals.
         norm : str, default "frac"
             The normalization of the periodogram. "abs" is absolute rms, "frac" is
             fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
@@ -438,7 +438,7 @@ class Powerspectrum(Crossspectrum):
 
     @staticmethod
     def from_events(events, dt, segment_size=None, norm="frac",
-                    silent=False, use_common_mean=True):
+                    silent=False, use_common_mean=True, gti=None):
         """Calculate AveragedPowerspectrum from an event list
 
         Parameters
@@ -466,15 +466,18 @@ class Powerspectrum(Crossspectrum):
             per-segment basis.
         silent : bool, default False
             Silence the progress bars
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Additional, optional, Good Time intervals, that get interesected with the 
+            GTIs of the input object.
         """
 
         return powerspectrum_from_events(
             events, dt, segment_size=segment_size, norm=norm,
-            silent=silent, use_common_mean=use_common_mean)
+            silent=silent, use_common_mean=use_common_mean, gti=gti)
 
     @staticmethod
     def from_lightcurve(lc, segment_size=None, norm="frac",
-                        silent=False, use_common_mean=True):
+                        silent=False, use_common_mean=True, gti=None):
         """Calculate AveragedPowerspectrum from a light curve
 
         Parameters
@@ -502,15 +505,18 @@ class Powerspectrum(Crossspectrum):
             per-segment basis.
         silent : bool, default False
             Silence the progress bars
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Additional, optional, Good Time intervals, that get interesected with the 
+            GTIs of the input object.
         """
 
         return powerspectrum_from_lightcurve(
             lc, segment_size=segment_size, norm=norm,
-            silent=silent, use_common_mean=use_common_mean)
+            silent=silent, use_common_mean=use_common_mean, gti=gti)
 
     @staticmethod
     def from_lc_iterable(iter_lc, dt, segment_size=None, norm="frac",
-                         silent=False, use_common_mean=True):
+                         silent=False, use_common_mean=True, gti=None):
         """Calculate AveragedCrossspectrum from two light curves
 
         Parameters
@@ -540,12 +546,78 @@ class Powerspectrum(Crossspectrum):
             per-segment basis.
         silent : bool, default False
             Silence the progress bars
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Good Time intervals.
         """
 
         return powerspectrum_from_lc_iterable(
             iter_lc, dt, segment_size=segment_size, norm=norm,
-            silent=silent, use_common_mean=use_common_mean)
+            silent=silent, use_common_mean=use_common_mean, gti=gti)
+
+    def _initialize_from_any_input(
+            self, data, dt=None, segment_size=None, norm="frac",
+            silent=False, use_common_mean=True, gti=None):
+        """Initialize the class, trying to understand the input types.
+
+        The input arguments are the same as ``__init__()``. Based on the type
+        of ``data``, this method will call the appropriate
+        ``powerspectrum_from_XXXX`` function, and initialize ``self`` with
+        the correct attributes.
+        """
+        if isinstance(data, EventList):
+            spec = powerspectrum_from_events(
+                data,
+                dt,
+                segment_size,
+                norm=norm.lower(),
+                silent=silent,
+                use_common_mean=use_common_mean,
+                gti=gti,
+            )
+        elif isinstance(data, Lightcurve):
+            spec = powerspectrum_from_lightcurve(
+                data,
+                segment_size,
+                norm=norm,
+                silent=silent,
+                use_common_mean=use_common_mean,
+                gti=gti,
+            )
+            spec.lc1 = data
+        elif isinstance(data, (tuple, list)):
+            if not isinstance(data[0], Lightcurve): # pragma: no cover
+                raise TypeError(f"Bad inputs to Powerspectrum: {type(data[0])}")
+            dt = data[0].dt
+            # This is a list of light curves.
+            spec = powerspectrum_from_lc_iterable(
+                data,
+                dt,
+                segment_size,
+                norm=norm,
+                silent=silent,
+                use_common_mean=use_common_mean,
+                gti=gti,
+            )
+        else: # pragma: no cover
+            raise TypeError(f"Bad inputs to Powerspectrum: {type(data)}")
 
+        for key, val in spec.__dict__.items():
+            setattr(self, key, val)
+        return
+
+    def _initialize_empty(self):
+        """Set all attributes to None."""
+        self.freq = None
+        self.power = None
+        self.power_err = None
+        self.unnorm_power = None
+        self.unnorm_power_err = None
+        self.df = None
+        self.dt = None
+        self.nphots1 = None
+        self.m = 1
+        self.n = None
+        return
 
 class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
     """
@@ -625,6 +697,10 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
     nphots: float
         The total number of photons in the light curve
 
+    legacy: bool
+        Use the legacy machinery of AveragedPowerspectrum. This might be useful to compare
+        with old results, and is also needed to use light curve lists as an input, to
+        conserve the spectra of each segment, or to use the large_data option.
     """
 
     def __init__(self, data=None, segment_size=None, norm="frac", gti=None,
@@ -654,23 +730,14 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
             )
 
         norm = norm.lower()
+        self.norm = norm
         self.dt = dt
         self.save_all = save_all
         self.segment_size = segment_size
         self.show_progress = not silent
 
         if not good_input:
-            self.freq = None
-            self.power = None
-            self.power_err = None
-            self.df = None
-            self.dt = None
-            self.nphots = None
-            self.nphots1 = None
-            self.m = 1
-            self.n = None
-            self.norm = norm
-            return
+            return self._initialize_empty()
 
         if isinstance(data, Generator):
             warnings.warn(
@@ -689,39 +756,9 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
             legacy = True
 
         if not legacy and data is not None:
-            if isinstance(data, EventList):
-                spec = powerspectrum_from_events(
-                    data,
-                    dt,
-                    segment_size,
-                    norm=norm.lower(),
-                    silent=silent,
-                    use_common_mean=use_common_mean,
-                )
-            elif isinstance(data, Lightcurve):
-                spec = powerspectrum_from_lightcurve(
-                    data,
-                    segment_size,
-                    norm=norm,
-                    silent=silent,
-                    use_common_mean=use_common_mean,
-                )
-                spec.lc1 = data
-            else:
-                dt = data[0].dt
-                # This is a list of light curves.
-                spec = powerspectrum_from_lc_iterable(
-                    data,
-                    dt,
-                    segment_size,
-                    norm=norm,
-                    silent=silent,
-                    use_common_mean=use_common_mean,
-                )
-
-            for key, val in spec.__dict__.items():
-                setattr(self, key, val)
-            return
+            return self._initialize_from_any_input(
+                data, dt=dt, segment_size=segment_size, norm=norm,
+                silent=silent, use_common_mean=use_common_mean)
 
         if large_data and data is not None:
             if not HAS_ZARR:
@@ -760,7 +797,8 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
             good = lengths >= segment_size
             data.gti = data.gti[good]
 
-        Powerspectrum.__init__(self, data, norm, gti=gti, dt=dt, skip_checks=True)
+        Powerspectrum.__init__(
+            self, data, norm, gti=gti, dt=dt, skip_checks=True, legacy=legacy)
 
         return
 
@@ -1110,7 +1148,7 @@ def powerspectrum_from_time_array(times, dt, segment_size=None, gti=None, norm="
 
 
 def powerspectrum_from_events(events, dt, segment_size=None, norm="frac",
-                              silent=False, use_common_mean=True):
+                              silent=False, use_common_mean=True, gti=None):
     """Calculate AveragedPowerspectrum from an event list
 
     Parameters
@@ -1137,6 +1175,9 @@ def powerspectrum_from_events(events, dt, segment_size=None, norm="frac",
         per-segment basis.
     silent : bool, default False
         Silence the progress bars
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Additional, optional, Good Time intervals, that get interesected with the 
+        GTIs of the input object.
 
     Returns
     -------
@@ -1151,7 +1192,8 @@ def powerspectrum_from_events(events, dt, segment_size=None, norm="frac",
 
 
 def powerspectrum_from_lightcurve(lc, segment_size=None, norm="frac",
-                                  silent=False, use_common_mean=True):
+                                  silent=False, use_common_mean=True,
+                                  gti=None):
     """Calculate AveragedPowerspectrum from a light curve
 
     Parameters
@@ -1178,6 +1220,9 @@ def powerspectrum_from_lightcurve(lc, segment_size=None, norm="frac",
         per-segment basis.
     silent : bool, default False
         Silence the progress bars
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Additional, optional, Good Time intervals, that get interesected with the 
+        GTIs of the input object.
 
     Returns
     -------
@@ -1187,17 +1232,21 @@ def powerspectrum_from_lightcurve(lc, segment_size=None, norm="frac",
     force_averaged = segment_size is not None
     # Suppress progress bar for single periodogram
     silent = silent or (segment_size is None)
+    err = None
+    if lc.err_dist == "gauss":
+        err = lc.counts_err
+
     table = avg_pds_from_events(
         lc.time, lc.gti, segment_size, lc.dt,
         norm=norm, use_common_mean=use_common_mean,
         silent=silent,
-        fluxes=lc.counts)
+        fluxes=lc.counts, errors=err)
 
     return _create_powerspectrum_from_result_table(table, force_averaged=force_averaged)
 
 
 def powerspectrum_from_lc_iterable(iter_lc, dt, segment_size=None, norm="frac",
-                                   silent=False, use_common_mean=True):
+                                   silent=False, use_common_mean=True, gti=None):
     """Calculate AveragedCrossspectrum from two light curves
 
     Parameters
@@ -1227,6 +1276,8 @@ def powerspectrum_from_lc_iterable(iter_lc, dt, segment_size=None, norm="frac",
         per-segment basis.
     silent : bool, default False
         Silence the progress bars
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. 
 
     Returns
     -------
@@ -1237,13 +1288,22 @@ def powerspectrum_from_lc_iterable(iter_lc, dt, segment_size=None, norm="frac",
     # Suppress progress bar for single periodogram
     silent = silent or (segment_size is None)
 
+    common_gti = gti
     def iterate_lc_counts(iter_lc):
         for lc in iter_lc:
             if hasattr(lc, "counts"):
                 n_bin = np.rint(segment_size / lc.dt).astype(int)
 
+                gti = lc.gti
+                if common_gti is not None:
+                    gti = cross_two_gtis(common_gti, lc.gti)
+                err = None
+                if lc.err_dist == "gauss":
+                    err = lc.counts_err
+
                 flux_iterable = get_flux_iterable_from_segments(
-                    lc.time, lc.gti, segment_size, n_bin, fluxes=lc.counts, errors=lc._counts_err
+                    lc.time, gti, segment_size, n_bin, fluxes=lc.counts,
+                    errors=err
                 )
                 for out in flux_iterable:
                     yield out


Test patch:
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 1ef200ca..90c8f62a 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -179,19 +179,18 @@ class TestAveragedCrossspectrumEvents(object):
     @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):
-
-        n = 10
+        n = 30
         lc_all = []
         for i in range(n):
             time = np.arange(0.0, 10.0, 10. / 10000)
-            counts = np.random.poisson(1000, size=time.shape[0])
+            counts = np.random.poisson(1000, size=time.size)
             lc = Lightcurve(time, counts)
             lc_all.append(lc)
 
         ps = AveragedCrossspectrum(lc_all, lc_all, 1.0, norm="leahy", legacy=legacy,
                                    use_common_mean=use_common_mean)
 
-        assert ps.m == 100
+        assert ps.m == 300
         assert np.isclose(np.mean(ps.power), 2.0, atol=1e-2, rtol=1e-2)
         assert np.isclose(np.std(ps.power), 2.0 / np.sqrt(ps.m), atol=0.1,
                           rtol=0.1)
@@ -385,7 +384,7 @@ class TestCoherence(object):
         t = np.arange(1280)
         a = np.random.poisson(100, len(t))
         lc = Lightcurve(t, a)
-        lc2 = Lightcurve(t, copy.copy(a))
+        lc2 = Lightcurve(t, copy.deepcopy(a))
 
         with pytest.warns(UserWarning) as record:
             c = AveragedCrossspectrum(lc, lc2, 128, use_common_mean=True)
@@ -597,10 +596,16 @@ class TestCrossspectrum(object):
         with pytest.raises(ValueError):
             cs = Crossspectrum(norm='frabs')
 
+    def test_init_with_wrong_lc_instance(self):
+        lc1_ = {"a": 1, "b": 2}
+        lc2_ = {"a": 1, "b": 2}
+        with pytest.raises(TypeError):
+            cs = Crossspectrum(lc1_, lc2_)
+
     def test_init_with_wrong_lc1_instance(self):
         lc_ = {"a": 1, "b": 2}
         with pytest.raises(TypeError):
-            cs = Crossspectrum(lc_, self.lc2)
+            cs = Crossspectrum(self.lc1, lc_)
 
     def test_init_with_wrong_lc2_instance(self):
         lc_ = {"a": 1, "b": 2}
@@ -612,25 +617,34 @@ class TestCrossspectrum(object):
         time = np.linspace(0.0, 1.0001, 10001)
         lc_ = Lightcurve(time, counts)
         with pytest.raises(StingrayError):
-            cs = Crossspectrum(self.lc1, lc_)
+            cs = Crossspectrum(self.lc1, lc_, legacy=True)
 
     def test_make_crossspectrum_diff_lc_stat(self):
-        lc_ = copy.copy(self.lc1)
+        lc_ = copy.deepcopy(self.lc1)
         lc_.err_dist = 'gauss'
 
         with pytest.warns(UserWarning) as record:
-            cs = Crossspectrum(self.lc1, lc_)
+            cs = Crossspectrum(self.lc1, lc_, legacy=True)
+        assert np.any(["different statistics" in r.message.args[0]
+                       for r in record])
+
+    def test_make_crossspectrum_diff_lc_iter_stat(self):
+        lc_ = copy.deepcopy(self.lc1)
+        lc_.err_dist = 'gauss'
+
+        with pytest.warns(UserWarning) as record:
+            cs = AveragedCrossspectrum([self.lc1], [lc_], segment_size=1, legacy=True)
         assert np.any(["different statistics" in r.message.args[0]
                        for r in record])
 
     def test_make_crossspectrum_bad_lc_stat(self):
-        lc1 = copy.copy(self.lc1)
+        lc1 = copy.deepcopy(self.lc1)
         lc1.err_dist = 'gauss'
-        lc2 = copy.copy(self.lc1)
+        lc2 = copy.deepcopy(self.lc1)
         lc2.err_dist = 'gauss'
 
         with pytest.warns(UserWarning) as record:
-            cs = Crossspectrum(lc1, lc2)
+            cs = Crossspectrum(lc1, lc2, legacy=True)
         assert np.any(["is not poisson" in r.message.args[0]
                        for r in record])
 
@@ -662,18 +676,19 @@ class TestCrossspectrum(object):
         assert type(new_cs) == type(self.cs)
         new_cs.time_lag()
 
-    def test_norm_abs(self):
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_norm_abs(self, legacy):
         # Testing for a power spectrum of lc1
-        cs = Crossspectrum(self.lc1, self.lc1, norm='abs')
+        cs = Crossspectrum(self.lc1, self.lc1, norm='abs', legacy=legacy)
         assert len(cs.power) == 4999
         assert cs.norm == 'abs'
         abs_noise = 2. * self.rate1  # expected Poisson noise level
         assert np.isclose(np.mean(cs.power[1:]), abs_noise)
 
-    def test_norm_leahy(self):
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_norm_leahy(self, legacy):
         # with pytest.warns(UserWarning) as record:
-        cs = Crossspectrum(self.lc1, self.lc1,
-                           norm='leahy')
+        cs = Crossspectrum(self.lc1, self.lc1, norm='leahy', legacy=legacy)
         assert len(cs.power) == 4999
         assert cs.norm == 'leahy'
         leahy_noise = 2.0  # expected Poisson noise level
@@ -818,8 +833,9 @@ class TestCrossspectrum(object):
         assert isinstance(pval, np.ndarray)
         assert pval.shape[0] == 2
 
-    def test_fullspec(self):
-        csT = Crossspectrum(self.lc1, self.lc2, fullspec=True)
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_fullspec(self, legacy):
+        csT = Crossspectrum(self.lc1, self.lc2, fullspec=True, legacy=legacy)
         assert csT.fullspec == True
         assert self.cs.fullspec == False
         assert csT.n == self.cs.n
@@ -935,7 +951,7 @@ class TestAveragedCrossspectrum(object):
 
         assert test_lc1.dt != test_lc2.dt
 
-        with pytest.raises(ValueError):
+        with pytest.raises(StingrayError):
             assert AveragedCrossspectrum(test_lc1, test_lc2, segment_size=1)
 
     def test_different_tseg(self):
@@ -1014,7 +1030,8 @@ class TestAveragedCrossspectrum(object):
             cs = AveragedCrossspectrum(self.lc1, self.lc2, segment_size=np.inf)
 
     @pytest.mark.parametrize("legacy", [False, True])
-    def test_with_iterable_of_lightcurves(self, legacy):
+    @pytest.mark.parametrize("err_dist", ["poisson", "gauss"])
+    def test_with_iterable_of_lightcurves(self, legacy, err_dist):
         def iter_lc(lc, n):
             "Generator of n parts of lc."
             t0 = int(len(lc) / n)
@@ -1027,6 +1044,9 @@ class TestAveragedCrossspectrum(object):
                     break
                 else:
                     i, t = t, t + t0
+        lc1 = copy.deepcopy(self.lc1)
+        lc2 = copy.deepcopy(self.lc2)
+        lc1.err_dist = lc2.err_dist = err_dist
         with pytest.warns(UserWarning) as record:
             cs = AveragedCrossspectrum(
                 iter_lc(self.lc1, 1), iter_lc(self.lc2, 1),
diff --git a/stingray/tests/test_largememory.py b/stingray/tests/test_largememory.py
index 7ca4fd4c..70f7cf69 100644
--- a/stingray/tests/test_largememory.py
+++ b/stingray/tests/test_largememory.py
@@ -456,9 +456,9 @@ class TestChunkPS(object):
     def test_invalid_data_to_pds(self):
         with pytest.raises(ValueError) as excinfo:
             AveragedPowerspectrum(
-                "sdfasfsa", segment_size=2048, large_data=True, silent=True
+                [self.lc1, self.lc1], segment_size=2048, large_data=True, silent=True
             )
-        assert "Invalid input data type: str" in str(excinfo.value)
+        assert "Invalid input data type: list" in str(excinfo.value)
 
     @pytest.mark.skipif("not HAS_ZARR")
     def test_events_to_cpds_unimplemented(self):
@@ -478,13 +478,13 @@ class TestChunkPS(object):
     def test_invalid_data_to_cpds(self):
         with pytest.raises(ValueError) as excinfo:
             AveragedCrossspectrum(
-                "sdfasfsa",
-                "sdfasfsa",
+                [self.lc1, self.lc1],
+                [self.lc2, self.lc2],
                 segment_size=4096,
                 large_data=True,
                 silent=True,
             )
-        assert "Invalid input data type: str" in str(excinfo.value)
+        assert "Invalid input data type: list" in str(excinfo.value)
 
     @pytest.mark.skipif("not HAS_ZARR")
     def test_calc_pds(self):
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index de6f3bf1..c8c82248 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -75,7 +75,8 @@ class TestAveragedPowerspectrumEvents(object):
     def test_from_lc_iter_works(self):
         pds_ev = AveragedPowerspectrum.from_lc_iterable(
             self.events.to_lc_iter(self.dt, self.segment_size),
-            segment_size=self.segment_size, dt=self.dt, norm="leahy", silent=True)
+            segment_size=self.segment_size, dt=self.dt, norm="leahy",
+            silent=True, gti=self.events.gti)
         assert np.allclose(self.leahy_pds.power, pds_ev.power)
 
     @pytest.mark.parametrize("err_dist", ["poisson", "gauss"])
@@ -110,6 +111,9 @@ class TestAveragedPowerspectrumEvents(object):
     def test_from_lc_iter_with_err_works(self):
         def iter_lc_with_errs(iter_lc):
             for lc in iter_lc:
+                # In order for error bars to be considered,
+                # err_dist has to be gauss.
+                lc.err_dist = "gauss"
                 lc._counts_err = np.zeros_like(lc.counts) + lc.counts.mean()**0.5
                 yield lc
 
@@ -120,6 +124,22 @@ class TestAveragedPowerspectrumEvents(object):
         power2 = self.leahy_pds.power.real
         assert np.allclose(power1, power2, rtol=0.01)
 
+    def test_from_lc_iter_with_err_ignored_with_wrong_err_dist(self):
+        def iter_lc_with_errs(iter_lc):
+            for lc in iter_lc:
+                # Not supposed to have error bars
+                lc.err_dist = "poisson"
+                # use a completely wrong error bar, for fun
+                lc._counts_err = np.zeros_like(lc.counts) + 14.2345425252462
+                yield lc
+
+        lccs = AveragedPowerspectrum.from_lc_iterable(
+            iter_lc_with_errs(self.events.to_lc_iter(self.dt, self.segment_size)),
+            segment_size=self.segment_size, dt=self.dt, norm='leahy', silent=True)
+        power1 = lccs.power.real
+        power2 = self.leahy_pds.power.real
+        assert np.allclose(power1, power2, rtol=0.01)
+
     def test_from_lc_iter_counts_only_works(self):
         def iter_lc_counts_only(iter_lc):
             for lc in iter_lc:
@@ -139,7 +159,7 @@ class TestAveragedPowerspectrumEvents(object):
             gti = np.array([[hdul[2].data["START"][0], hdul[2].data["STOP"][0]]])
 
             _ = AveragedPowerspectrum.from_time_array(
-                times,  segment_size=128, dt=self.dt, gti=gti, norm='none',
+                times, segment_size=128, dt=self.dt, gti=gti, norm='none',
                 use_common_mean=False)
 
     @pytest.mark.parametrize("norm", ["frac", "abs", "none", "leahy"])
@@ -147,9 +167,10 @@ class TestAveragedPowerspectrumEvents(object):
         lc = self.events.to_lc(self.dt)
         lc._counts_err = np.sqrt(lc.counts.mean()) + np.zeros_like(lc.counts)
         pds = AveragedPowerspectrum.from_lightcurve(
-            lc, segment_size=self.segment_size, norm=norm, silent=True)
+            lc, segment_size=self.segment_size, norm=norm, silent=True,
+            gti=lc.gti)
         pds_ev = AveragedPowerspectrum.from_events(
-            self.events, segment_size=self.segment_size, dt=self.dt, norm=norm, silent=True)
+            self.events, segment_size=self.segment_size, dt=self.dt, norm=norm, silent=True, gti=self.events.gti)
         for attr in ["power", "freq", "m", "n", "nphots", "segment_size"]:
             assert np.allclose(getattr(pds, attr), getattr(pds_ev, attr))
 
@@ -157,10 +178,11 @@ class TestAveragedPowerspectrumEvents(object):
         with pytest.raises(ValueError):
             assert AveragedPowerspectrum(self.lc, dt=self.dt)
 
-    def test_init_with_nonsense_segment(self):
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_init_with_nonsense_segment(self, legacy):
         segment_size = "foo"
         with pytest.raises(TypeError):
-            assert AveragedPowerspectrum(self.lc, segment_size, dt=self.dt)
+            assert AveragedPowerspectrum(self.lc, segment_size, dt=self.dt, legacy=legacy)
 
     def test_init_with_none_segment(self):
         segment_size = None
@@ -276,12 +298,13 @@ class TestPowerspectrum(object):
         assert ps.m == 1
         assert ps.n is None
 
-    def test_make_periodogram_from_lightcurve(self):
-        ps = Powerspectrum(self.lc)
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_make_periodogram_from_lightcurve(self, legacy):
+        ps = Powerspectrum(self.lc, legacy=legacy)
         assert ps.freq is not None
         assert ps.power is not None
         assert ps.power_err is not None
-        assert ps.df == 1.0 / self.lc.tseg
+        assert np.isclose(ps.df, 1.0 / self.lc.tseg)
         assert ps.norm == "frac"
         assert ps.m == 1
         assert ps.n == self.lc.time.shape[0]
@@ -300,15 +323,23 @@ class TestPowerspectrum(object):
         with pytest.raises(TypeError):
             assert Powerspectrum(self.lc.counts)
 
-    def test_init_with_nonsense_data(self):
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_init_with_nonsense_list(self, legacy):
         nonsense_data = [None for i in range(100)]
         with pytest.raises(TypeError):
-            assert Powerspectrum(nonsense_data)
+            assert Powerspectrum(nonsense_data, legacy=legacy)
+
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_init_with_nonsense_data(self, legacy):
+        nonsense_data = 34
+        with pytest.raises(TypeError):
+            assert Powerspectrum(nonsense_data, legacy=legacy)
 
-    def test_init_with_nonsense_norm(self):
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_init_with_nonsense_norm(self, legacy):
         nonsense_norm = "bla"
         with pytest.raises(ValueError):
-            assert Powerspectrum(self.lc, norm=nonsense_norm)
+            assert Powerspectrum(self.lc, norm=nonsense_norm, legacy=legacy)
 
     def test_init_with_wrong_norm_type(self):
         nonsense_norm = 1.0
@@ -475,8 +506,9 @@ class TestPowerspectrum(object):
         """
         pass
 
-    def test_rebin_makes_right_attributes(self):
-        ps = Powerspectrum(self.lc, norm="Leahy")
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_rebin_makes_right_attributes(self, legacy):
+        ps = Powerspectrum(self.lc, norm="Leahy", legacy=legacy)
         # replace powers
         ps.power = np.ones_like(ps.power) * 2.0
 
@@ -486,7 +518,7 @@ class TestPowerspectrum(object):
         assert bin_ps.freq is not None
         assert bin_ps.power is not None
         assert bin_ps.power is not None
-        assert bin_ps.df == rebin_factor * 1.0 / self.lc.tseg
+        assert np.isclose(bin_ps.df, rebin_factor * 1.0 / self.lc.tseg)
         assert bin_ps.norm.lower() == "leahy"
         assert bin_ps.m == 2
         assert bin_ps.n == self.lc.time.shape[0]
@@ -660,15 +692,17 @@ class TestAveragedPowerspectrum(object):
         with pytest.raises(ValueError):
             assert AveragedPowerspectrum(self.lc)
 
-    def test_init_with_nonsense_segment(self):
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_init_with_nonsense_segment(self, legacy):
         segment_size = "foo"
         with pytest.raises(TypeError):
-            assert AveragedPowerspectrum(self.lc, segment_size)
+            assert AveragedPowerspectrum(self.lc, segment_size, legacy=legacy)
 
-    def test_init_with_none_segment(self):
+    @pytest.mark.parametrize("legacy", [True, False])
+    def test_init_with_none_segment(self, legacy):
         segment_size = None
         with pytest.raises(ValueError):
-            assert AveragedPowerspectrum(self.lc, segment_size)
+            assert AveragedPowerspectrum(self.lc, segment_size, legacy=legacy)
 
     def test_init_with_inf_segment(self):
         segment_size = np.inf
@@ -732,7 +766,8 @@ class TestAveragedPowerspectrum(object):
         with pytest.warns(UserWarning) as record:
             cs = AveragedPowerspectrum(
                 iter_lc(self.lc, 1),
-                segment_size=1, legacy=legacy)
+                segment_size=1, legacy=legacy,
+                gti=self.lc.gti)
         message = "The averaged Power spectrum from a generator "
 
         assert np.any([message in r.message.args[0]
@@ -744,7 +779,16 @@ class TestAveragedPowerspectrum(object):
         lc.gti = gti
         lc_split = lc.split_by_gti()
 
-        cs = AveragedPowerspectrum(lc_split, segment_size=0.05, norm="leahy")
+        cs = AveragedPowerspectrum(lc_split, segment_size=0.05, norm="leahy",
+                                   gti=lc.gti)
+        cs_lc = AveragedPowerspectrum(lc, segment_size=0.05, norm="leahy",
+                                      gti=lc.gti)
+        for attr in ("power", "unnorm_power", "power_err", "unnorm_power_err",
+                     "freq"):
+            assert np.allclose(getattr(cs, attr), getattr(cs_lc, attr))
+
+        for attr in ("m", "n", "norm"):
+            assert getattr(cs, attr) == getattr(cs_lc, attr)
 
     @pytest.mark.parametrize('df', [2, 3, 5, 1.5, 1, 85])
     def test_rebin(self, df):