Closed moble closed 3 years ago
Converting to the frequency domain is very tricky. You have to take a few more steps:
Putting these together, I would do something more like this:
import numpy as np
from scipy import signal
import sxs
import matplotlib.pyplot as plt
waveform = sxs.load('SXS:BBH:0123/Lev/rhOverM', extrapolation_order=2)
# Ignore junk radiation (data before t=0)
t_start = 0.0
# Interpolate to uniform time steps. Here, I'm just choosing the smallest
# timestep in the input data. If you want more physical results, you'll
# need to scale to the appropriate total mass, and a corresponding physical
# timestep.
t_end = waveform.t[-1]
dt = np.min(np.diff(waveform.t))
t_uniform = np.arange(t_start, t_end, dt)
h = waveform.interpolate(t_uniform)
# Evaluate in a particular direction. I'm arbitrarily choosing a direction
# (and implicitly a polarization angle).
theta, phi = 0.0, 0.0
h = h.evaluate(theta, phi)
# Take the real part (like a single detector would measure)
h = h.real
# Multiply by a window function. I like Tukey windows because, if you make
# them narrow enough, they leave the inspiral and merger data unchanged.
# Adjust the 0.12 below as needed for a particular waveform.
window = signal.windows.tukey(h.shape[0], alpha=0.12, sym=True)
h = h * window
# Now we can convert to the frequency domain
h_tilde = np.fft.fft(h)
freq = np.fft.fftfreq(h.size, dt)
# And plot it
plt.loglog(freq, np.abs(h_tilde))
plt.xlabel('Frequency (units of $1/M$)')
plt.ylabel(r'$|\tilde{h}|$');
The result looks like this:
The actual data starts around a frequency of 6e-3. A physical binary system would extend to earlier times and frequencies, but our numerical waveforms are finite, and that happens to be where this system starts. At smaller frequencies, you basically have nonsense that you should ignore (or you could try harder to window the data to make it smaller). Then, the around 6e-3, we start to see frequencies that are in the actual data. At first, you see some oscillations due to remaining Gibbs effects, but then it becomes smoother and more physically realistic. (Note that lower frequencies have larger amplitude in the frequency domain, even though they have smaller amplitude in the time domain, because the signal spends longer at lower frequencies, so more power gets built up there.) You get a Fourier transform that slopes downward to around 1e-1, where it sort of flattens out a little, and then drops steeply. That little flat part around 1e-1 is the merger, and the drop is because of ringdown — the physical system just never really oscillates at any higher frequencies. Above about 2e-1, you just have numerical noise that you could maybe reduce with better windowing.
Hi Mr Boyle,
When I do the waveform.evaluate, it means that I do the sum of all modes of (h_lm * -2Y_lm)
in the some direction (theta, phi)? And If I want to evaluate the (2, 2) mode of waveform, so I just can slice the waveform to
mode (2, 2) and then evaluate the waveform?
And if I want to add the pattern function to the waveform, if it is right that I can do like this:
t_start = 0.0
t_end = waveform.t[-1]
dt = np.min(np.diff(waveform.t))
t_uniform = np.arange(t_start, t_end, dt)
h = waveform.interpolate(t_uniform)
h = f_plus * h.real+f_cross * h.imag
#for scaling the waveform with mass and distance
t=t_uniform*( G * mass / c ** 3)
h=h* (G * mass / R / c ** 2)
and then can I do the fft to transform the waveform to the frequency domain?
Thank you!
When I do the waveform.evaluate, it means that I do the sum of all modes of
(h_lm * -2Y_lm)
in the some direction (theta, phi)?
That's right.
And If I want to evaluate the (2, 2) mode of waveform, so I just can slice the waveform to mode (2, 2) and then evaluate the waveform?
You can't slice to a single mode, because the code assumes that for any l
in the (l, m)
modes, all possible m
values are present in the data. (You can try such a slice and get a result, which is a bug that's currently open, but the result won't work properly.) Instead, you should make a copy of the waveform, and just set the modes you don't want to 0. Also, you probably want both the (2, 2)
and (2, -2)
modes, so here's what I would do:
# Copy the entire waveform
waveform2 = waveform.copy()
# Set ALL data to 0
waveform2.data[...] = 0.0
# Copy back the modes you want to keep
waveform2.data[:, waveform2.index(2, -2)] = waveform.data[:, waveform.index(2, -2)]
waveform2.data[:, waveform2.index(2, 2)] = waveform.data[:, waveform.index(2, 2)]
pattern function
You mean polarization function? Basically, yes, but you have to be careful because — as is standard — we represent the waveform as the complex quantity h
, which is actually h_plus - 1j * h_cross
. Note the minus sign! That is, h_plus
is the same as h.real
, but h_cross
is -h.imag
.
Thank you very much for giving me the suggestions and answering my questions. But now I also have a very simple question when I try to calculate the SNR for aLIGO.
For calculating the SNR, I used the code as you teach me:
import matplotlib.pyplot as plt
import numpy as np
import sxs
from scipy import signal
extrapolation_order = 2
waveform = sxs.load("SXS:BBH:0123/Lev/rhOverM", extrapolation_order=extrapolation_order, download=False)
t_start = 0.0
t_end = waveform.t[-1]
dt = np.min(np.diff(waveform.t))
t_uniform = np.arange(t_start, t_end, dt)
h = waveform.interpolate(t_uniform)
#Evaluate the waveform
theta, phi = 0.0, 0.0
h = h.evaluate(theta, phi)
#For the polarization function of LIGO, and I just choose the angle (theta, phi, psi)
#as (np.pi, np.pi, np.pi)
f_plus = 0.5 * (1 + np.cos(np.pi) ** 2.0) * np.cos(2.0 * np.pi) * np.cos(2.0 * np.pi) \
- np.cos(np.pi) * np.sin(2.0 * np.pi) * np.sin(2.0 * np.pi)
f_cross = 0.5 * (1 + np.cos(np.pi) ** 2.0) * np.cos(2.0 * np.pi) * np.sin(2.0 * np.pi) \
+ np.cos(np.pi) * np.sin(2.0 * np.pi) * np.cos(2.0 * np.pi)
h_ligo = f_plus * h.real + f_cross * (- h.imag)
#Then I do the fft with Tukey window and scale the waveform with mass and distance
#like GW150914 event
c = 299792458 # m/s
G = 6.67259e-11 # m^(3)kg^(-1)s^(-2)
R = 410 * 1e6 * 3.08567758149137e16 # m
mass = 65 * 1.98847e30 #kg
h = h_ligo
window = signal.windows.tukey(h.shape[0], alpha=0.12, sym=True)
h = h * window
freqs = np.fft.fftfreq(h.size, dt)[0:int(len(h)/2)] / (G * mass / c ** 3)
h_f = np.fft.fft(h)[0:int(len(h)/2)] * (G * mass / R / c ** 2)
The result of the fft:
For the aLIGO noise curve I just use the aLIGO design curve :
asd = np.loadtxt('/mnt/hgfs/CK3/ligo.txt').transpose()
Then I interpolate the frequency domain waveform to the aLIGO frequency like this:
from scipy import interpolate
freqh = freqs
hfi= h_f
hf_interpolate = interpolate.interp1d(freqh, hfi, kind='cubic')
#Then I calculate the SNR^2 for the aLIGO used SNR^2=4 * df * np.sum(|h_f|^2/S_n)
habs = (hf_interpolate(asd[0])) * np.conjugate(hf_interpolate(asd[0]))
z1 = np.real(habs)/(asd[1] ** 2)
SNR2 = 4 * (asd[0][1]-asd[0][0]) * np.sum(z1)
And for the result of SNR square is about 2.525e12
which is impossible. So I think I must doing wrong or I miss some details of the waveform, but I can't figure it out.
I'm very sorry for writing such a large piece of codes for this simple question to bother you, and I'm so appreciate your help in fft and SXS waveform.
Thank you!
There are a few things to note here.
(asd[0][1]-asd[0][0])
isn't enough to do the integral correctly. You can check this for your data by plotting np.diff(asd[0])
. If it's basically constant, you're fine, but otherwise you'd probably be better off constructing another spline and using it to integrate.asd_interp = interp1d(asd[0], asd[1], kind='cubic', bounds_error=False, fill_value=(np.inf, np.inf))(freqs)
This also has the nice effect of allowing you to skip item 2 above, because freqs
is evenly spaced in frequency.
rfft
and rfftfreq
. sxs
module comes with a few useful constants that are kept up to date, have somewhat higher accuracy, and have more descriptive names — plus you reduce the chances of typos if you use them.Putting this all together, I would do something more like this:
R = 410_000_000 * sxs.parsec_in_meters
mass = 65 * sxs.m_sun_in_seconds
s_dt = dt * mass
s = h_ligo * mass * sxs.speed_of_light / R
window = signal.windows.tukey(s.shape[0], alpha=0.12, sym=True)
freqs = np.fft.rfftfreq(s.size, s_dt)
s_f = s_dt * np.fft.rfft(s * window)
asd_interp = interpolate.interp1d(
asd[0], asd[1], kind='cubic', bounds_error=False, fill_value=(np.inf, np.inf)
)(freqs)
SNR = np.sqrt(4 * (freqs[1] - freqs[0]) * np.sum(abs(s_f)**2 / asd_interp**2))
The result is a much more reasonable SNR of 86 with the latest aLIGO design curve. Noting that this is optimally oriented, and comparing to results from GWTC-2, this looks about right.
I check the result of interpolating the waveform onto the noise curve, the noise of waveform is really amplified, and the Gibbs effect also increase. And the np.diff(asd[0])
is actually not a constant, thank you for pointing out my mistake.
I really appreciate your help and suggestions, these really help me a lot.
Hi Mr. Boyle, sorry to bother you again. I get some other waveforms like NRSurdq4 and IMRPhenomD, I want use some functions of sxs
package like sxs.waveforms.memory.add_memory
and waveform.evaluate
, so how can I transform the waveform I get to the WaveformModes
of sxs
?
Thank you!
You would have to get the raw data (modes and time) as numpy arrays, and then construct a WaveformModes
object by filling out the parameters relevant to your data.
Note that IMRPhenom waveforms are already in the frequency domain, so you would have to inverse-fft them to get to the time domain, to do things like add_memory
.
Sorry for bothering, I reinstall my operation system and conda install -c conda-forge sxs, then I import sxs and the program report a error:
AttributeError: Failed in nopython mode pipeline (step: nopython rewrites) module ‘spherical’ has no attribute '_linear_matrix_offset’.
How can I solve this problem?
You must be using the wrong python executable, so you'll have to track this down on your own, because I don't have enough information to figure it out. Basically, whatever version of python you're actually using has an old version of the spherical
module; recent versions of that package don't even have _linear_matrix_offset
. Maybe you didn't un-install your old version of conda correctly, or maybe your PATH is pointing to the wrong python
executable, or any number of things that I can't guess from here.
Thank you very much!
Sorry for bothering you again and again with these simple questions, I want to hit the wall with my head for my stupid. And I'm too foolish to understand how to construct WaveformModes
object use IMRPhenom waveform, so I decide to write memory calculation program by myself with refer Eq. (17b) in K. Mitman et al..
If I don't work with a wrong understanding, the product of two SWSH of spin-weighted -2 is -2Ylm * \bar(-2Ylm) = 0Ylm
, so for calculating memory just of (2, 2) mode, I do some easy calculation like:
0.125 * ethbar**2 * D^-1 * 0Ylm * integrate(hlm.dot * \bar(hlm.dot))
with
D^-1 * 0Ylm = (0.125 * (l-1) * l * (l+1) * (l+2)) * 0Ylm
and
ehtbar**2 * 0Ylm = \sqrt((l-1) * (l + 1 + 1)) * \sqrt(l * (l + 1)) * -2Ylm
and for the (2, 2) mode I finally get
h_mem = \sqrt(4) * \sqrt(6) * -2Y_22 * \integrate(hdot * \bar(hdot))
For the python code, which is really hard for me, and I always make mistakes in programming. I write the code like this:
import matplotlib.pyplot as plt
import numpy as np
import sxs
extrapolation_order = 3
w = sxs.load("SXS:BBH:0305/Lev/rhOverM", extrapolation_order=extrapolation_order, download=False)
# I sliced waveform to (2, 2) mode
w22 = w[:, w.index(2, 2)]
# and I interpolate the waveform with w.t
from scipy import interpolate
w22_inter = interpolate(w22.t, w22, "Cubic")
# try to calculate memory of (2, 2) mode
dt = w.t[1] - w.t[0]
hdot = np.gradient(h, dt)
hmode = hdot * np.conjugate(hdot)
Dinverse = ((1/8) * (2-1) * 2 * (2+1) * (2+2)) ** (-1)
ethbar = -np.sqrt((2+0)*(2-0+1)) *(-np.sqrt((2-1)*(2+1+1)))
# doing the integrate of homde with time
hint = dt * np.cumsum(np.real(hmode))*ethbar * Dinverse/8
for the result compared with waveform.memory.add_memory, it looks like this:
I think I'm doing bad with interpolate the waveform, and I also read your code about interpolate, spherical and so on, these codes are really hard for me to understand with my python level.
Thank you for helping me again and again, and sorry for bothering you again and again.
the product of two SWSH of spin-weighted -2 is
-2Ylm * \bar(-2Ylm) = 0Ylm
, so for calculating memory just of (2, 2) mode
That's wrong. You can use the expression for spin-weighted spherical harmonics in terms of Wigner D matrices to evaluate the product, and get an expression involving Clebsch-Gordan coefficients (or equivalently Wigner 3-j). But if you're integrating the product, you can just use orthonormality.
Also, I haven't followed everything that you're doing in your code, but remember that the time step is not constant, so you can't integrate in time with cumsum
.
Hi Mr. Boyle,
If I'm not work with wrong calculation, for the product of SWSH -2Y22 * conjugate(-2Y22)
, I will get:
-2Y22 * 2Y2(-2) = sum(conjugate(0YL0) * (-1)**(2 + 2 + L) * sqrt(25 * (2*L+1)/4*pi) * two Wigner 3-j symbols),
for L in [0, 5],
with l1 = 2, l2 = 2, l3 = L, m1 = 2, m2 = -2, m3 = 0, s1 = -2, s2 = 2, s3 = 0.
so the result of this product is the sum of five sYlm
with some coefficients, and I also use your package spherical
to evaluate the Wigner 3-j symbol, only L = 5
the symbol is equal to zero. So why the result of calculating displacement memory is only 0Y20
left?
So why the result of calculating displacement memory is only
0Y20
left?
That's not the only thing left, though it is usually the largest. You just have to finish the calculation and plug in some numbers to see why.
[A user asked this question through a separate channel recently. It'll be easier to answer it here.]
I need to do the fft to transform time domain waveform to the frequency domain. But I tried many times and I can't get the correct waveform of frequency domain. For the python code to do fft, I used:
For the results of doing fft of (2, 2) mode, it have much higher oscillation in the high frequency domain with the loglog plot of absolute value of h(f), and I think the trend of waveform with increasing of frequency is also incorrect.