pycbc icon indicating copy to clipboard operation
pycbc copied to clipboard

filter.overlap does not behave correctly when input series do not have the same length

Open Sbozzolo opened this issue 3 years ago • 0 comments

I would expect the function overlap to be symmetric in the input: overlap(h1, h2) = overlap(h2, h1). In its current implementation, this is not true. I am attaching a minimum working example to show that.

I believe that the reason for the asymmetry is in how the cutoff indices are selected (they are always extracted looking at delta_f and N of the second series). I suspect that the results that the function currently returns are wrong in both cases.

If the function is supposed to be used only for arrays with the same length, it would be better to check that and raise an error.

The overlaps as returned by the program below are 0.01255448903587639 -0.22526849934969062, where the numbers are overlap(1, 2) and overlap(2, 1).

from pycbc.filter import get_cutoff_indices, make_frequency_series, overlap
from pycbc.types import timeseries as pycbcts
from pycbc.waveform import get_td_waveform

fmin_gw = 50
fmax_gw = 100
delta_t = 1/4096

hp1, hc1 = get_td_waveform(approximant="IMRPhenomPv2",
                           mass1=10,
                           mass2=10,
                           delta_t=delta_t,
                           f_lower=40)

hp2, hc2 = get_td_waveform(approximant="IMRPhenomPv2",
                           mass1=10,
                           mass2=25,
                           delta_t=delta_t,
                           f_lower=40)

h1_pycbc = pycbcts.TimeSeries(hp1, delta_t=hp1.delta_t)
h2_pycbc = pycbcts.TimeSeries(hp2, delta_t=hp2.delta_t)

overlap12 = overlap(h1_pycbc, h2_pycbc, None, fmin_gw, fmax_gw)
overlap21 = overlap(h2_pycbc, h1_pycbc, None, fmin_gw, fmax_gw)

print(overlap12, overlap21)  

h1_freq = make_frequency_series(h1_pycbc)
h2_freq = make_frequency_series(h2_pycbc)

kmin1, kmax1 = get_cutoff_indices(fmin_gw, fmax_gw,
                                  h1_freq.delta_f, len(h1_freq))
kmin2, kmax2 = get_cutoff_indices(fmin_gw, fmax_gw,
                                  h2_freq.delta_f, len(h2_freq))

print(kmin1, kmax1)
print(kmin2, kmax2)
print(h1_freq.delta_f, h2_freq.delta_f)

Sbozzolo avatar Sep 16 '20 01:09 Sbozzolo