Edited docs(core.rst,lscs,lsps), modified function, docstrings, forma…

…tting,,new tests

docs/core.rst

@@ -85,3 +85,20 @@ Multi-taper Periodogram
    :maxdepth: 2
 
    notebooks/Multitaper/multitaper_example.ipynb
+
+
+Lomb Scargle Crossspectrum
+--------------------------
+.. toctree::
+   :maxdepth: 2
+
+   notebooks/LombScargle/LombScargleCrossspectrum_tutorial.ipynb
+
+
+Lomb Scargle Powerspectrum
+--------------------------
+
+.. toctree::
+   :maxdepth: 2
+
+   notebooks/LombScargle/LombScarglePowerspectrum_tutorial.ipynb

stingray/fourier.py

@@ -2126,11 +2126,6 @@ def impose_symmetry_lsft(
     freqs_new : np.array
         The new frequencies
     """
-    zero_present = np.amy(freqs == 0)
-    if not zero_present:
-        ft_res_new = np.concatenate([np.conjugate(np.flip(ft_res)), 0, ft_res])
-        freqs_new = np.concatenate([-np.flip(freqs), 0, freqs])
-    else:
-        ft_res_new = np.concatenate([np.conjugate(np.flip(ft_res))[1:], ft_res])
-        freqs_new = np.concatenate([-np.flip(freqs)[1:], freqs])
+    ft_res_new = np.concatenate([np.conjugate(np.flip(ft_res)), [0.0], ft_res])
+    freqs_new = np.concatenate([np.flip(-freqs), [0.0], freqs])
     return ft_res_new, freqs_new

stingray/lombscargle.py

@@ -475,31 +475,11 @@ def _ls_cross(self, lc1, lc2, freq=None, fullspec=False, method="fast", oversamp
             lsft1 = lsft_slow(lc1.counts, lc1.time, freq)
             lsft2 = lsft_slow(lc2.counts, lc2.time, freq)
         elif method == "fast":
-            lsft1 = lsft_fast(
-                lc1.counts,
-                lc1.time,
-                freq,
-                oversampling=oversampling,
-            )
-            lsft2 = lsft_fast(
-                lc2.counts,
-                lc2.time,
-                freq,
-                oversampling=oversampling,
-            )
+            lsft1 = lsft_fast(lc1.counts, lc1.time, freq, oversampling=oversampling)
+            lsft2 = lsft_fast(lc2.counts, lc2.time, freq, oversampling=oversampling)
         if fullspec:
-            lsft1, _ = impose_symmetry_lsft(
-                lsft1,
-                np.sum((lc1.counts)),
-                lc1.n,
-                freq,
-            )
-            lsft2, freq = impose_symmetry_lsft(
-                lsft2,
-                np.sum(lc2.counts),
-                lc2.n,
-                freq,
-            )
+            lsft1, _ = impose_symmetry_lsft(lsft1, np.sum((lc1.counts)), lc1.n, freq)
+            lsft2, freq = impose_symmetry_lsft(lsft2, np.sum(lc2.counts), lc2.n, freq)
         cross = np.multiply(lsft1, np.conjugate(lsft2))
         return freq, cross
 
@@ -590,10 +570,10 @@ class LombScarglePowerspectrum(LombScargleCrossspectrum):
         The light curve data to be Fourier-transformed.
 
     norm: {``frac``, ``abs``, ``leahy``, ``none``}, string, optional, default ``none``
-        The normalization of the cross spectrum.
+        The normalization of the power spectrum.
 
     power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``all``
-        Parameter to choose among complete, real part and magnitude of the cross spectrum.
+        Parameter to choose among complete, real part and magnitude of the power spectrum.
 
     fullspec: boolean, optional, default ``False``
         If False, keep only the positive frequencies, or if True, keep all of them .

stingray/tests/test_fourier.py

@@ -541,3 +541,45 @@ def test_unnormalize_poisson_noise(self, norm, power_type):
         noise_notnorm = poisson_level("none", self.meanrate, self.nph)
 
         assert np.isclose(noise_notnorm, unnorm_noise)
+
+
+def test_lsft_slow_fast():
+    np.random.seed(0)
+    rand = np.random.default_rng(42)
+    n = 1000
+    t = np.sort(rand.random(n)) * np.sqrt(n)
+    y = np.sin(2 * np.pi * 3.0 * t)
+    sub = np.min(y)
+    y -= sub
+    freqs = np.fft.fftfreq(n, np.median(np.diff(t, 1)))
+    freqs = freqs[freqs >= 0]
+    lsftslow = lsft_slow(y, t, freqs, sign=1)
+    lsftfast = lsft_fast(y, t, freqs, sign=1, oversampling=10)
+    assert np.argmax(lsftslow) == np.argmax(lsftfast)
+    assert round(freqs[np.argmax(lsftslow)], 1) == round(freqs[np.argmax(lsftfast)], 1) == 3.0
+    assert np.all(
+        np.all((lsftslow * np.conjugate(lsftslow)).imag == 0)
+        & np.all((lsftfast * np.conjugate(lsftfast)).imag == 0)
+    )
+
+
+def test_impose_symmetry_lsft():
+    np.random.seed(0)
+    rand = np.random.default_rng(42)
+    n = 1000
+    t = np.sort(rand.random(n)) * np.sqrt(n)
+    y = np.sin(2 * np.pi * 3.0 * t)
+    sub = np.min(y)
+    y -= sub
+    freqs = np.fft.fftfreq(n, np.median(np.diff(t, 1)))
+    freqs = freqs[freqs >= 0]
+    lsftslow = lsft_slow(y, t, freqs, sign=1)
+    lsftfast = lsft_fast(y, t, freqs, sign=1, oversampling=5)
+    imp_sym_slow, freqs_new_slow = impose_symmetry_lsft(lsftslow, 0, n, freqs)
+    imp_sym_fast, freqs_new_fast = impose_symmetry_lsft(lsftfast, 0, n, freqs)
+    assert imp_sym_slow.shape == imp_sym_fast.shape == freqs_new_fast.shape == freqs_new_slow.shape
+    assert np.all((imp_sym_slow.real) == np.flip(imp_sym_slow.real))
+    assert np.all((imp_sym_slow.imag) == -np.flip(imp_sym_slow.imag))
+    assert np.all((imp_sym_fast.real) == np.flip(imp_sym_fast.real))
+    assert np.all((imp_sym_fast.imag) == (-np.flip(imp_sym_fast.imag)))
+    assert np.all(freqs_new_slow == freqs_new_fast)