From ef714029f67b4dd1427e0d3a10689f2e810c2757 Mon Sep 17 00:00:00 2001 From: pupperemeritus Date: Wed, 23 Aug 2023 20:20:56 +0530 Subject: [PATCH] Edited docs(core.rst,lscs,lsps), modified function, docstrings, formatting,,new tests --- docs/core.rst | 17 ++++++++++++++ stingray/fourier.py | 9 ++------ stingray/lombscargle.py | 32 +++++--------------------- stingray/tests/test_fourier.py | 42 ++++++++++++++++++++++++++++++++++ 4 files changed, 67 insertions(+), 33 deletions(-) diff --git a/docs/core.rst b/docs/core.rst index e9f3353f7..56ec7088f 100644 --- a/docs/core.rst +++ b/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 \ No newline at end of file diff --git a/stingray/fourier.py b/stingray/fourier.py index 2eca42daf..9866dfcd1 100644 --- a/stingray/fourier.py +++ b/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 diff --git a/stingray/lombscargle.py b/stingray/lombscargle.py index e285e190a..22c29f20e 100644 --- a/stingray/lombscargle.py +++ b/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 . diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py index 2f127adfe..2a85f5191 100644 --- a/stingray/tests/test_fourier.py +++ b/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)