Skip to content

coherence

Coherence

Compute the frequency-dependent coherence between two light curves or GP models.

This class estimates the coherence spectrum, which quantifies the degree of linear correlation between two time series as a function of frequency. Coherence values range from 0 to 1, with values near 1 indicating a strong linear relationship at that frequency.

Inputs can be either LightCurve objects or trained GaussianProcess models from this package. If GP models are provided and posterior samples already exist, those are used. If no samples exist, 1000 GP realizations will be generated automatically on a 1000-point grid.

If both inputs are GP models, the coherence is computed for each sample pair and the mean and standard deviation across samples are returned. Otherwise, coherence is computed on the raw input light curves.

Poisson noise bias correction is supported and may be enabled to correct for uncorrelated noise.

Parameters:

Name Type Description Default
lc_or_model1 LightCurve or GaussianProcess

First input light curve or trained GP model.

 required
lc_or_model2 LightCurve or GaussianProcess

Second input light curve or trained GP model.

 required
fmin float or auto

Minimum frequency for the coherence spectrum. If 'auto', uses the lowest nonzero FFT frequency.

 'auto'
fmax float or auto

Maximum frequency. If 'auto', uses the Nyquist frequency.

 'auto'
num_bins int

Number of frequency bins.

 None
bin_type str

Type of frequency binning ('log' or 'linear').

 'log'
bin_edges array - like

Custom frequency bin edges.

 []
subtract_noise_bias bool

Whether to subtract Poisson noise bias from the coherence spectrum.

 True
bkg1 float

Background count rate for lightcurve 1 (used in noise bias correction).

 0
bkg2 float

Background count rate for lightcurve 2.

 0

Attributes:

Name Type Description
freqs array - like

Frequency bin centers.
freq_widths array - like

Widths of each frequency bin.
cohs array - like

Coherence values.
coh_errors array - like

Uncertainties in the coherence values.
Source code in stela_toolkit/coherence.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
class Coherence:
    """
    Compute the frequency-dependent coherence between two light curves or GP models.

    This class estimates the coherence spectrum, which quantifies the degree of linear correlation
    between two time series as a function of frequency. Coherence values range from 0 to 1,
    with values near 1 indicating a strong linear relationship at that frequency.

    Inputs can be either LightCurve objects or trained GaussianProcess models from this package.
    If GP models are provided and posterior samples already exist, those are used.
    If no samples exist, 1000 GP realizations will be generated automatically on a 1000-point grid.

    If both inputs are GP models, the coherence is computed for each sample pair and the
    mean and standard deviation across samples are returned. Otherwise, coherence is computed
    on the raw input light curves.

    Poisson noise bias correction is supported and may be enabled to correct for uncorrelated noise.

    Parameters
    ----------
    lc_or_model1 : LightCurve or GaussianProcess
        First input light curve or trained GP model.

    lc_or_model2 : LightCurve or GaussianProcess
        Second input light curve or trained GP model.

    fmin : float or 'auto', optional
        Minimum frequency for the coherence spectrum. If 'auto', uses the lowest nonzero FFT frequency.

    fmax : float or 'auto', optional
        Maximum frequency. If 'auto', uses the Nyquist frequency.

    num_bins : int, optional
        Number of frequency bins.

    bin_type : str, optional
        Type of frequency binning ('log' or 'linear').

    bin_edges : array-like, optional
        Custom frequency bin edges.

    subtract_noise_bias : bool, optional
        Whether to subtract Poisson noise bias from the coherence spectrum.

    bkg1 : float, optional
        Background count rate for lightcurve 1 (used in noise bias correction).

    bkg2 : float, optional
        Background count rate for lightcurve 2.

    Attributes
    ----------
    freqs : array-like
        Frequency bin centers.

    freq_widths : array-like
        Widths of each frequency bin.

    cohs : array-like
        Coherence values.

    coh_errors : array-like
        Uncertainties in the coherence values.
    """

    def __init__(self,
                 lc_or_model1,
                 lc_or_model2,
                 fmin='auto',
                 fmax='auto',
                 num_bins=None,
                 bin_type="log",
                 bin_edges=[],
                 subtract_noise_bias=True,
                 bkg1=0,
                 bkg2=0):

        input_data = _CheckInputs._check_lightcurve_or_model(lc_or_model1)
        if input_data['type'] == 'model':
            self.times1, self.rates1 = input_data['data']
        else:
            self.times1, self.rates1, _ = input_data['data']

        input_data = _CheckInputs._check_lightcurve_or_model(lc_or_model2)
        if input_data['type'] == 'model':
            self.times2, self.rates2 = input_data['data']
        else:
            self.times2, self.rates2, _ = input_data['data']
        _CheckInputs._check_input_bins(num_bins, bin_type, bin_edges)

        if not np.allclose(self.times1, self.times2):
            raise ValueError("The time arrays of the two light curves must be identical.")

        # Use absolute min and max frequencies if set to 'auto'
        self.dt = np.diff(self.times1)[0]
        self.fmin = np.fft.rfftfreq(len(self.rates1), d=self.dt)[1] if fmin == 'auto' else fmin
        self.fmax = np.fft.rfftfreq(len(self.rates1), d=self.dt)[-1] if fmax == 'auto' else fmax  # nyquist frequency

        self.num_bins = num_bins
        self.bin_type = bin_type
        self.bin_edges = bin_edges

        self.bkg1 = bkg1
        self.bkg2 = bkg2

        # Check if the input rates are for multiple realizations
        # this needs to be corrected for handling different shapes and dim val1 != dim val2
        # namely for multiple observations
        if len(self.rates1.shape) == 2 and len(self.rates2.shape) == 2:
            coherence_spectrum = self.compute_stacked_coherence()
        else:
            coherence_spectrum = self.compute_coherence(subtract_noise_bias=subtract_noise_bias)

        self.freqs, self.freq_widths, self.cohs, self.coh_errors = coherence_spectrum

    def compute_coherence(self, times1=None, rates1=None, times2=None, rates2=None, subtract_noise_bias=True):
        """
        Compute the coherence spectrum between two light curves.

        Parameters
        ----------
        times1, rates1 : array-like, optional
            Time and rate values for the first light curve. Defaults to object attributes.
        times2, rates2 : array-like, optional
            Time and rate values for the second light curve. Defaults to object attributes.
        subtract_noise_bias : bool, optional
            Whether to subtract the estimated noise bias.

        Returns
        -------
        freqs : array-like
            Frequency bin centers.
        freq_widths : array-like
            Frequency bin widths.
        coherence : array-like
            Coherence spectrum.
        None
            Reserved for compatibility (with the stacked method).
        """

        times1 = self.times1 if times1 is None else times1
        rates1 = self.rates1 if rates1 is None else rates1
        times2 = self.times2 if times2 is None else times2
        rates2 = self.rates2 if rates2 is None else rates2

        lc1 = LightCurve(times=times1, rates=rates1)
        lc2 = LightCurve(times=times2, rates=rates2)
        cross_spectrum = CrossSpectrum(
            lc1, lc2,
            fmin=self.fmin, fmax=self.fmax,
            num_bins=self.num_bins, bin_type=self.bin_type, bin_edges=self.bin_edges
        )

        power_spectrum1 = PowerSpectrum(
            lc1,
            fmin=self.fmin, fmax=self.fmax,
            num_bins=self.num_bins, bin_type=self.bin_type, bin_edges=self.bin_edges
        )
        power_spectrum2 = PowerSpectrum(
            lc2,
            fmin=self.fmin, fmax=self.fmax,
            num_bins=self.num_bins, bin_type=self.bin_type, bin_edges=self.bin_edges
        )

        ps1 = power_spectrum1.powers
        ps2 = power_spectrum2.powers
        cs = cross_spectrum.cs

        if subtract_noise_bias:
            bias = self.compute_bias(ps1, ps2)
        else:
            bias = 0

        cohs = (np.abs(cs) ** 2 - bias) / (ps1 * ps2)

        M = self.count_frequencies_in_bins()
        coh_errors = np.sqrt((1 - cohs) / (2 * cohs * M))

        return power_spectrum1.freqs, power_spectrum1.freq_widths, cohs, coh_errors

    def compute_stacked_coherence(self):
        """
        Compute the coherence from stacked realizations of the light curves.

        For multiple realizations (GP samples), this method computes the
        coherence for each pair of realizations and returns the mean and standard deviation.

        Returns
        -------
        freqs : array-like
            Frequency bin centers.
        freq_widths : array-like
            Frequency bin widths.
        coherence_mean : array-like
            Mean coherence spectrum across realizations.
        coherence_std : array-like
            Standard deviation of the coherence across realizations.
        """

        coherences = []
        for i in range(self.rates1.shape[0]):
            coherence_spectrum = self.compute_coherence(times1=self.times1, rates1=self.rates1[i],
                                                        times2=self.times2, rates2=self.rates2[i],
                                                        subtract_noise_bias=False
                                                    )
            freqs, freq_widths, coherence, _ = coherence_spectrum
            coherences.append(coherence)

        coherences = np.vstack(coherences)
        coherences_mean = np.mean(coherences, axis=0)
        coherences_std = np.std(coherences, axis=0)

        return freqs, freq_widths, coherences_mean, coherences_std

    def compute_bias(self, power_spectrum1, power_spectrum2):
        """
        Estimate the Poisson noise bias for the coherence calculation. 

        Parameters
        ----------
        power_spectrum1 : array-like
            Power spectrum of the first light curve.
        power_spectrum2 : array-like
            Power spectrum of the second light curve.

        Returns
        -------
        bias : array-like
            Estimated noise bias per frequency bin.
        """

        mean1 = np.mean(self.rates1)
        mean2 = np.mean(self.rates2)

        pnoise1 = 2 * (mean1 + self.bkg1) / mean1 ** 2
        pnoise2 = 2 * (mean2 + self.bkg2) / mean2 ** 2

        bias = (
            pnoise2 * (power_spectrum1 - pnoise1)
            + pnoise1 * (power_spectrum2 - pnoise2)
            + pnoise1 * pnoise2
        )
        num_freq = self.count_frequencies_in_bins()
        bias /= num_freq
        return bias

    def plot(self, freqs=None, freq_widths=None, cohs=None, coh_errors=None, **kwargs):
        """
        Plot the coherence spectrum.

        Parameters
        ----------
        **kwargs : dict
            Additional keyword arguments for plot customization (e.g., xlabel, xscale).
        """

        freqs = self.freqs if freqs is None else freqs
        freq_widths = self.freq_widths if freq_widths is None else freq_widths
        cohs = self.cohs if cohs is None else cohs
        coh_errors = self.coh_errors if coh_errors is None else coh_errors

        kwargs.setdefault('xlabel', 'Frequency')
        kwargs.setdefault('ylabel', 'Coherence')
        kwargs.setdefault('xscale', 'log')
        Plotter.plot(
            x=freqs, y=cohs, xerr=freq_widths, yerr=coh_errors, **kwargs
        )

    def count_frequencies_in_bins(self, fmin=None, fmax=None, num_bins=None, bin_type=None, bin_edges=[]):
        """
        Counts the number of frequencies in each frequency bin.
        Wrapper method to use FrequencyBinning.count_frequencies_in_bins with class attributes.
        """

        return FrequencyBinning.count_frequencies_in_bins(
            self, fmin=fmin, fmax=fmax, num_bins=num_bins, bin_type=bin_type, bin_edges=bin_edges
        )

compute_bias(power_spectrum1, power_spectrum2)

Estimate the Poisson noise bias for the coherence calculation.

Parameters:

Name Type Description Default
power_spectrum1 array - like

Power spectrum of the first light curve.

 required
power_spectrum2 array - like

Power spectrum of the second light curve.

 required

Returns:

Name Type Description
bias array - like

Estimated noise bias per frequency bin.
Source code in stela_toolkit/coherence.py 
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
def compute_bias(self, power_spectrum1, power_spectrum2):
    """
    Estimate the Poisson noise bias for the coherence calculation. 

    Parameters
    ----------
    power_spectrum1 : array-like
        Power spectrum of the first light curve.
    power_spectrum2 : array-like
        Power spectrum of the second light curve.

    Returns
    -------
    bias : array-like
        Estimated noise bias per frequency bin.
    """

    mean1 = np.mean(self.rates1)
    mean2 = np.mean(self.rates2)

    pnoise1 = 2 * (mean1 + self.bkg1) / mean1 ** 2
    pnoise2 = 2 * (mean2 + self.bkg2) / mean2 ** 2

    bias = (
        pnoise2 * (power_spectrum1 - pnoise1)
        + pnoise1 * (power_spectrum2 - pnoise2)
        + pnoise1 * pnoise2
    )
    num_freq = self.count_frequencies_in_bins()
    bias /= num_freq
    return bias

compute_coherence(times1=None, rates1=None, times2=None, rates2=None, subtract_noise_bias=True)

Compute the coherence spectrum between two light curves.

Parameters:

Name Type Description Default
times1 array - like

Time and rate values for the first light curve. Defaults to object attributes.

 None
rates1 array - like

Time and rate values for the first light curve. Defaults to object attributes.

 None
times2 array - like

Time and rate values for the second light curve. Defaults to object attributes.

 None
rates2 array - like

Time and rate values for the second light curve. Defaults to object attributes.

 None
subtract_noise_bias bool

Whether to subtract the estimated noise bias.

 True

Returns:

Name Type Description
freqs array - like

Frequency bin centers.
freq_widths array - like

Frequency bin widths.
coherence array - like

Coherence spectrum.
None

Reserved for compatibility (with the stacked method).
Source code in stela_toolkit/coherence.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
def compute_coherence(self, times1=None, rates1=None, times2=None, rates2=None, subtract_noise_bias=True):
    """
    Compute the coherence spectrum between two light curves.

    Parameters
    ----------
    times1, rates1 : array-like, optional
        Time and rate values for the first light curve. Defaults to object attributes.
    times2, rates2 : array-like, optional
        Time and rate values for the second light curve. Defaults to object attributes.
    subtract_noise_bias : bool, optional
        Whether to subtract the estimated noise bias.

    Returns
    -------
    freqs : array-like
        Frequency bin centers.
    freq_widths : array-like
        Frequency bin widths.
    coherence : array-like
        Coherence spectrum.
    None
        Reserved for compatibility (with the stacked method).
    """

    times1 = self.times1 if times1 is None else times1
    rates1 = self.rates1 if rates1 is None else rates1
    times2 = self.times2 if times2 is None else times2
    rates2 = self.rates2 if rates2 is None else rates2

    lc1 = LightCurve(times=times1, rates=rates1)
    lc2 = LightCurve(times=times2, rates=rates2)
    cross_spectrum = CrossSpectrum(
        lc1, lc2,
        fmin=self.fmin, fmax=self.fmax,
        num_bins=self.num_bins, bin_type=self.bin_type, bin_edges=self.bin_edges
    )

    power_spectrum1 = PowerSpectrum(
        lc1,
        fmin=self.fmin, fmax=self.fmax,
        num_bins=self.num_bins, bin_type=self.bin_type, bin_edges=self.bin_edges
    )
    power_spectrum2 = PowerSpectrum(
        lc2,
        fmin=self.fmin, fmax=self.fmax,
        num_bins=self.num_bins, bin_type=self.bin_type, bin_edges=self.bin_edges
    )

    ps1 = power_spectrum1.powers
    ps2 = power_spectrum2.powers
    cs = cross_spectrum.cs

    if subtract_noise_bias:
        bias = self.compute_bias(ps1, ps2)
    else:
        bias = 0

    cohs = (np.abs(cs) ** 2 - bias) / (ps1 * ps2)

    M = self.count_frequencies_in_bins()
    coh_errors = np.sqrt((1 - cohs) / (2 * cohs * M))

    return power_spectrum1.freqs, power_spectrum1.freq_widths, cohs, coh_errors

compute_stacked_coherence()

Compute the coherence from stacked realizations of the light curves.

For multiple realizations (GP samples), this method computes the coherence for each pair of realizations and returns the mean and standard deviation.

Returns:

Name Type Description
freqs array - like

Frequency bin centers.
freq_widths array - like

Frequency bin widths.
coherence_mean array - like

Mean coherence spectrum across realizations.
coherence_std array - like

Standard deviation of the coherence across realizations.
Source code in stela_toolkit/coherence.py 
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def compute_stacked_coherence(self):
    """
    Compute the coherence from stacked realizations of the light curves.

    For multiple realizations (GP samples), this method computes the
    coherence for each pair of realizations and returns the mean and standard deviation.

    Returns
    -------
    freqs : array-like
        Frequency bin centers.
    freq_widths : array-like
        Frequency bin widths.
    coherence_mean : array-like
        Mean coherence spectrum across realizations.
    coherence_std : array-like
        Standard deviation of the coherence across realizations.
    """

    coherences = []
    for i in range(self.rates1.shape[0]):
        coherence_spectrum = self.compute_coherence(times1=self.times1, rates1=self.rates1[i],
                                                    times2=self.times2, rates2=self.rates2[i],
                                                    subtract_noise_bias=False
                                                )
        freqs, freq_widths, coherence, _ = coherence_spectrum
        coherences.append(coherence)

    coherences = np.vstack(coherences)
    coherences_mean = np.mean(coherences, axis=0)
    coherences_std = np.std(coherences, axis=0)

    return freqs, freq_widths, coherences_mean, coherences_std

count_frequencies_in_bins(fmin=None, fmax=None, num_bins=None, bin_type=None, bin_edges=[])

Counts the number of frequencies in each frequency bin. Wrapper method to use FrequencyBinning.count_frequencies_in_bins with class attributes.

Source code in stela_toolkit/coherence.py 
278
279
280
281
282
283
284
285
286
def count_frequencies_in_bins(self, fmin=None, fmax=None, num_bins=None, bin_type=None, bin_edges=[]):
    """
    Counts the number of frequencies in each frequency bin.
    Wrapper method to use FrequencyBinning.count_frequencies_in_bins with class attributes.
    """

    return FrequencyBinning.count_frequencies_in_bins(
        self, fmin=fmin, fmax=fmax, num_bins=num_bins, bin_type=bin_type, bin_edges=bin_edges
    )

plot(freqs=None, freq_widths=None, cohs=None, coh_errors=None, **kwargs)

Plot the coherence spectrum.

Parameters:

Name Type Description Default
**kwargs dict

Additional keyword arguments for plot customization (e.g., xlabel, xscale).

 {}
Source code in stela_toolkit/coherence.py 
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
def plot(self, freqs=None, freq_widths=None, cohs=None, coh_errors=None, **kwargs):
    """
    Plot the coherence spectrum.

    Parameters
    ----------
    **kwargs : dict
        Additional keyword arguments for plot customization (e.g., xlabel, xscale).
    """

    freqs = self.freqs if freqs is None else freqs
    freq_widths = self.freq_widths if freq_widths is None else freq_widths
    cohs = self.cohs if cohs is None else cohs
    coh_errors = self.coh_errors if coh_errors is None else coh_errors

    kwargs.setdefault('xlabel', 'Frequency')
    kwargs.setdefault('ylabel', 'Coherence')
    kwargs.setdefault('xscale', 'log')
    Plotter.plot(
        x=freqs, y=cohs, xerr=freq_widths, yerr=coh_errors, **kwargs
    )