From b596182fd0f40ad17e864bce64bf81b3c328088e Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Mon, 21 Feb 2022 14:40:30 +0100 Subject: [PATCH] Add formulas for bispectrum in fourier.py --- stingray/bispectrum.py | 98 ++++++++++++++++++++++++++++++++++++++++++ stingray/fourier.py | 59 +++++++++++++++++++++++++ 2 files changed, 157 insertions(+) diff --git a/stingray/bispectrum.py b/stingray/bispectrum.py index 693207ffd..93b822ab9 100644 --- a/stingray/bispectrum.py +++ b/stingray/bispectrum.py @@ -15,6 +15,7 @@ from stingray import lightcurve import stingray.utils as utils +from stingray.fourier import get_flux_iterable_from_segments, bispectrum_and_bicoherence_from_iterable __all__ = ["Bispectrum"] @@ -450,3 +451,100 @@ def plot_phase(self, axis=None, save=False, filename=None): else: plt.savefig(filename) return plt + + +def bispectrum_and_bicoherence_from_time_array(times, gti, segment_size, dt): + """Calculate the bispectrum and the bicoherence from an array of times. + + Parameters + ---------- + times : float `np.array` + Array of times in the reference band + gti : [[gti00, gti01], [gti10, gti11], ...] + common good time intervals + segment_size : float + length of segments + dt : float + Time resolution of the light curves used to produce periodograms + + Returns + ------- + freqs : `np.array` + The frequency in each row or column + bispectrum: 2-d `np.array` + The unnormalized bispectrum + bicoherence: 2-d `np.array` + The bicoherence, calculated with the normalization from + Kim & Powers (1979), IEEE Trans. Plasma Sci., PS-7, 120 + """ + + n_bin = np.rint(segment_size / dt).astype(int) + dt = segment_size / n_bin + + flux_iterable = get_flux_iterable_from_segments( + times, gti, segment_size, n_bin + ) + + return bispectrum_and_bicoherence_from_iterable(flux_iterable, dt) + + +def bispectrum_and_bicoherence_from_events(events, segment_size, dt): + """Calculate the bispectrum and the bicoherence from an input event list. + + Parameters + ---------- + events : `EventList` + The input event list + segment_size : float + length of segments + dt : float + Time resolution of the light curves used to produce periodograms + + Returns + ------- + freqs : `np.array` + The frequency in each row or column + bispectrum: 2-d `np.array` + The unnormalized bispectrum + bicoherence: 2-d `np.array` + The bicoherence, calculated with the normalization from + Kim & Powers (1979), IEEE Trans. Plasma Sci., PS-7, 120 + """ + + times = events.time + gti = events.gti + + return bispectrum_and_bicoherence_from_time_array(times, gti, segment_size, dt) + + +def bispectrum_and_bicoherence_from_lightcurve(lightcurve, segment_size): + """Calculate the bispectrum and the bicoherence from an array of times. + + Parameters + ---------- + lightcurve : `Lightcurve` + Input light curve + segment_size : float + length of segments + dt : float + Time resolution of the light curves used to produce periodograms + + Returns + ------- + freqs : `np.array` + The frequency in each row or column + bispectrum: 2-d `np.array` + The unnormalized bispectrum + bicoherence: 2-d `np.array` + The bicoherence, calculated with the normalization from + Kim & Powers (1979), IEEE Trans. Plasma Sci., PS-7, 120 + """ + dt = lightcurve.dt + n_bin = np.rint(segment_size / dt).astype(int) + + flux_iterable = get_flux_iterable_from_segments( + lightcurve.time, lightcurve.gti, segment_size, + n_bin, fluxes=lightcurve.counts + ) + + return bispectrum_and_bicoherence_from_iterable(flux_iterable, dt) diff --git a/stingray/fourier.py b/stingray/fourier.py index 544b9fdd0..308913d7a 100644 --- a/stingray/fourier.py +++ b/stingray/fourier.py @@ -889,6 +889,65 @@ def get_flux_iterable_from_segments(times, gti, segment_size, n_bin=None, fluxes yield cts +from stingray.fourier import positive_fft_bins, fft, get_flux_iterable_from_segments +from tqdm import tqdm as show_progress + + +def bispectrum_and_bicoherence_from_iterable(flux_iterable, dt): + """Calculate the bispectrum and the bicoherence. + + Parameters + ---------- + flux_iterable : `iterable` of `np.array`s or of tuples (`np.array`, `np.array`) + Iterable providing either equal-length series of count measurements, or + of tuples (fluxes, errors). They must all be of the same length. + dt : float + Time resolution of the light curves used to produce periodograms + + Returns + ------- + freqs : `np.array` + The frequency in each row or column + bispectrum: 2-d `np.array` + The unnormalized bispectrum + bicoherence: 2-d `np.array` + The bicoherence, calculated with the normalization from + Kim & Powers (1979), IEEE Trans. Plasma Sci., PS-7, 120 + """ + + B = B1 = B2 = None + for flux in show_progress(flux_iterable): + if flux is None: + continue + + if isinstance(flux, tuple): + flux, _ = flux + + if B is None: + n_bin = flux.size + fgt0 = positive_fft_bins(n_bin) + Nft = fgt0.stop - fgt0.start + B = np.zeros((Nft, Nft), dtype=complex) + B1 = np.zeros((Nft, Nft)) + B2 = np.zeros((Nft, Nft)) + freqs = np.fft.fftfreq(n_bin, dt)[fgt0] + + ft = fft(flux)[fgt0] + + M = 0 + for i in range(B.shape[0]): + b1 = ft * ft[i] + b2 = np.roll(ft, -i) + B[i, :] += b1 * b2.conj() + B1[i, :] += (b1 * b1.conj()).real + B2[i, :] += (b2 * b2.conj()).real + M += 1 + + BC = (B * B.conj()).real / (B1 * B2) + + return freqs, B / M, BC + + def avg_pds_from_iterable(flux_iterable, dt, norm="frac", use_common_mean=True, silent=False): """Calculate the average periodogram from an iterable of light curves