# Signal processing¶

Signal processing

kerrgeodesic_gw.signal_processing.fourier(signal)

Compute the Fourier transform of a time signal.

The Fourier transform of the time signal $$h(t)$$ is defined by

$\tilde{h}(f) = \int_{-\infty}^{+\infty} h(t)\, \mathrm{e}^{-2\pi \mathrm{i} f t} \, \mathrm{d}t$

INPUT:

• signal – list of pairs $$(t, h(t))$$, where $$t$$ is the time and $$h(t)$$ is the signal at $$t$$. NB: the sampling in $$t$$ must be uniform

OUTPUT:

• a numerical approximation (FFT) of the Fourier transform, as a list of pairs $$(f,\tilde{h}(f))$$

EXAMPLES:

Fourier transform of the $$h_+$$ signal from a particle orbiting at the ISCO of a Schwarzschild black hole:

sage: from kerrgeodesic_gw import h_particle_signal, fourier
sage: a, r0 = 0., 6.
sage: theta, phi = pi/2, 0
sage: sig = h_particle_signal(a, r0, theta, phi, 0., 800., nb_points=200)
sage: sig[:5]  # tol 1.0e-13
[(0.000000000000000, 0.1536656546005028),
(4.02010050251256, 0.03882960191740132),
(8.04020100502512, -0.0791773803745108),
(12.0603015075377, -0.18368354016995483),
(16.0804020100502, -0.25796165580196795)]
sage: ft = fourier(sig)
sage: ft[:5]   # tol 1.0e-13
[(-0.12437500000000001, (1.0901773159466113+0j)),
(-0.12313125000000001, (1.0901679813096925+0.021913447210468399j)),
(-0.12188750000000001, (1.0901408521337845+0.043865034665617614j)),
(-0.12064375000000001, (1.0900987733011678+0.065894415355341171j)),
(-0.11940000000000001, (1.0900473072498438+0.088044567109304195j))]


Plot of the norm of the Fourier transform:

sage: ft_norm = [(f, abs(hf)) for (f, hf) in ft]
sage: line(ft_norm, axes_labels=[r'$f M$', r'$|\tilde{h}_+(f)| r/(\mu M)$'],
....:      gridlines=True, frame=True, axes=False)
Graphics object consisting of 1 graphics primitive


The first peak is at the orbital frequency of the particle:

sage: f0 = n(1/(2*pi*r0^(3/2))); f0
0.0108291222393566


while the highest peak is at twice this frequency ($$m=2$$ mode).

kerrgeodesic_gw.signal_processing.max_detectable_radius(a, mu, theta, psd, BH_time_scale, distance, t_obs_yr=1, snr_threshold=10, r_min=None, r_max=200, m_min=1, m_max=None, approximation=None)

Compute the maximum orbital radius $$r_{0,\rm max}$$ at which a particle of given mass is detectable.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• mu – mass of the orbiting object, in solar masses

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• psd – function with a single argument ($$f$$) representing the detector’s one-sided noise power spectral density $$S_n(f)$$ (see e.g. lisa_detector.power_spectral_density())

• BH_time_scale – value of $$M$$ in the same time unit as $$S_n(f)$$; if $$S_n(f)$$ is provided in $$\mathrm{Hz}^{-1}$$, then BH_time_scale must be $$M$$ expressed in seconds.

• distance – distance $$r$$ to the detector, in parsecs

• t_obs_yr – (default: 1) observation period, in years

• snr_threshold – (default: 10) signal-to-noise value above which a detection is claimed

• r_min – (default: None) lower bound of the search interval for $$r_{0,\rm max}$$ (in units of $$M$$); if None then the ISCO value is used

• r_max – (default: 200) upper bound of the search interval for $$r_{0,\rm max}$$ (in units of $$M$$)

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• m_max – (default: None) upper bound in the summation over the Fourier mode $$m$$; if None, m_max is set to 10 for $$r_0 \leq 20 M$$ and to 5 for $$r_0 > 20 M$$

• approximation – (default: None) string describing the computational method for the signal; allowed values are

OUTPUT:

• Boyer-Lindquist radius (in units of $$M$$) of the most remote orbit for which the signal-to-noise ratio is larger than snr_threshold during the observation time t_obs_yr

EXAMPLES:

Maximum orbital radius for LISA detection of a 1 solar-mass object orbiting Sgr A* observed by LISA, assuming a BH spin $$a=0.9 M$$ and a vanishing inclination angle:

sage: from kerrgeodesic_gw import (max_detectable_radius, lisa_detector,
....:                              astro_data)
sage: a = 0.9
sage: mu = 1
sage: theta = 0
sage: psd = lisa_detector.power_spectral_density_RCLfit
sage: BH_time_scale = astro_data.SgrA_mass_s
sage: distance = astro_data.SgrA_distance_pc
sage: max_detectable_radius(a, mu, theta, psd, BH_time_scale, distance)  # tol 1.0e-13
46.98292057022975


Lowering the SNR threshold to 5:

sage: max_detectable_radius(a, mu, theta, psd, BH_time_scale, distance,  # tol 1.0e-13
....:                       snr_threshold=5)
53.503386416644645


Lowering the data acquisition time to 1 day:

sage: max_detectable_radius(a, mu, theta, psd, BH_time_scale, distance,  # tol 1.0e-13
....:                       t_obs_yr=1./365.25)
27.158725683511463


Assuming an inclination angle of $$\pi/2$$:

sage: theta = pi/2
sage: max_detectable_radius(a, mu, theta, psd, BH_time_scale, distance)  # tol 1.0e-13
39.81826500968452


Using the 1.5-PN approximation (a has to be zero and m_max has to be at most 5):

sage: a = 0
sage: max_detectable_radius(a, mu, theta, psd, BH_time_scale,  # tol 1.0e-13
....:                       distance, approximation='1.5PN', m_max=5)
39.74273792957649


Read a signal from a data file.

INPUT:

• filename – string; name of the file

• dirname – (default: None) string; name of directory where the file is located

OUTPUT:

• signal as a list of pairs $$(t, h(t))$$

EXAMPLES:

sage: from kerrgeodesic_gw import save_signal, read_signal
sage: sig0 = [(RDF(i/10), RDF(sin(i/5))) for i in range(5)]
sage: from sage.misc.temporary_file import tmp_filename
sage: filename = tmp_filename(ext='.dat')
sage: save_signal(sig0, filename)
sage: sig   # tol 1.0e-13
[(0.0, 0.0),
(0.1, 0.19866933079506122),
(0.2, 0.3894183423086505),
(0.3, 0.5646424733950354),
(0.4, 0.7173560908995228)]


A test:

sage: sig == sig0
True

kerrgeodesic_gw.signal_processing.save_signal(sig, filename, dirname=None)

Write a signal in a data file.

INPUT:

• sig – signal as a list of pairs $$(t, h(t))$$

• filename – string; name of the file

• dirname – (default: None) string; name of directory where the file is located

EXAMPLES:

sage: from kerrgeodesic_gw import save_signal, read_signal
sage: sig = [(RDF(i/10), RDF(sin(i/5))) for i in range(5)]
sage: sig   # tol 1.0e-13
[(0.0, 0.0),
(0.1, 0.19866933079506122),
(0.2, 0.3894183423086505),
(0.3, 0.5646424733950354),
(0.4, 0.7173560908995228)]
sage: from sage.misc.temporary_file import tmp_filename
sage: filename = tmp_filename(ext='.dat')
sage: save_signal(sig, filename)


A test:

sage: sig1 = read_signal(filename)
sage: sig1 == sig
True


Evaluate the signal-to-noise ratio of a signal observed in a detector of a given power spectral density.

The signal-to-noise ratio $$\rho$$ of the time signal $$h(t)$$ is computed according to the formula

(1)$\rho^2 = 4 \int_{0}^{+\infty} \frac{|\tilde{h}(f)|^2}{S_n(f)} \, \mathrm{d}f$

where $$\tilde{h}(f)$$ is the Fourier transform of $$h(t)$$ (see fourier()) and $$S_n(f)$$ is the detector’s one-sided noise power spectral density (see e.g. lisa_detector.power_spectral_density()).

INPUT:

• signal – list of pairs $$(t, h(t))$$, where $$t$$ is the time and $$h(t)$$ is the signal at $$t$$. NB: the sampling in $$t$$ must be uniform

• time_scale – value of $$t$$ unit in terms of $$S_n(f)$$ unit; if $$S_n(f)$$ is provided in $$\mathrm{Hz}^{-1}$$, then time_scale must be the unit of $$t$$ in signal expressed in seconds.

• psd – function with a single argument ($$f$$) representing the detector’s one-sided noise power spectral density $$S_n(f)$$

• fmin – lower bound used instead of $$0$$ in the integral (1)

• fmax – upper bound used instead of $$+\infty$$ in the integral (1)

• scale – (default: 1) scale factor by which $$h(t)$$ must be multiplied to get the actual signal

• interpolation – (default: 'linear') string describing the type of interpolation used to evaluate $$|h(f)|^2$$ from the list resulting from the FFT of signal; allowed values are

• 'linear': linear interpolation between two data points

• 'spline': cubic spline interpolation

• quad_epsrel – (default: 1.e-6) relative error tolerance in the computation of the integral (1)

• quad_limit – (default: 500) upper bound on the number of subintervals used in the adaptive algorithm to compute the integral (this corresponds to the argument limit of SciPy’s function quad)

OUTPUT:

• the signal-to-noise ratio $$\rho$$ computed via the integral (1), with the boundaries $$0$$ and $$+\infty$$ replaced by respectively fmin and fmax.

EXAMPLES:

Let us evaluate the SNR of the gravitational signal generated by a 1-solar mass object orbiting at the ISCO of Sgr A* observed by LISA during 1 day. We need the following functions:

sage: from kerrgeodesic_gw import (h_particle_signal, signal_to_noise,
....:                              lisa_detector, astro_data)


We model Sgr A* as a Schwarzschild black hole and we consider that the signal is constituted by the mode $$h_+(t)$$ observed in the orbital plane:

sage: a, r0 = 0., 6.
sage: theta, phi = pi/2, 0
sage: tmax = 24*3600/astro_data.MSgrA_s  # 1 day in units of Sgr A* mass
sage: h = h_particle_signal(a, r0, theta, phi, 0., tmax, mode='+',
....:                       nb_points=4000)


The signal-to-noise ratio is then computed as:

sage: time_scale = astro_data.MSgrA_s  # Sgr A* mass in seconds
sage: psd = lisa_detector.power_spectral_density_RCLfit
sage: fmin, fmax = 1e-5, 5e-3
sage: mu_ov_r = astro_data.Msol_m / astro_data.dSgrA  # mu/r
sage: signal_to_noise(h, time_scale, psd, fmin, fmax,     # tol 1.0e-13
....:                 interpolation='spline', scale=mu_ov_r)
7583.104124621172


Signal-to-noise for a signal computed at the quadrupole approximation:

sage: h = h_particle_signal(a, r0, theta, phi, 0., tmax, mode='+',
sage: signal_to_noise(h, time_scale, psd, fmin, fmax,     # tol 1.0e-13
....:                 interpolation='spline', scale=mu_ov_r)
5383.264004811493

kerrgeodesic_gw.signal_processing.signal_to_noise_particle(a, r0, theta, psd, t_obs, BH_time_scale, m_min=1, m_max=None, scale=1, approximation=None)

Evaluate the signal-to-noise ratio of gravitational radiation emitted by a single orbiting particle observed in a detector of a given power spectral density.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• psd – function with a single argument ($$f$$) representing the detector’s one-sided noise power spectral density $$S_n(f)$$ (see e.g. lisa_detector.power_spectral_density())

• t_obs – observation period, in the same time unit as $$S_n(f)$$

• BH_time_scale – value of $$M$$ in the same time unit as $$S_n(f)$$; if $$S_n(f)$$ is provided in $$\mathrm{Hz}^{-1}$$, then BH_time_scale must be $$M$$ expressed in seconds.

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• m_max – (default: None) upper bound in the summation over the Fourier mode $$m$$; if None, m_max is set to 10 for $$r_0 \leq 20 M$$ and to 5 for $$r_0 > 20 M$$

• scale – (default: 1) scale factor by which $$h(t)$$ must be multiplied to get the actual signal; this should by $$\mu/r$$, where $$\mu$$ is the particle mass and $$r$$ the radial coordinate of the detector

• approximation – (default: None) string describing the computational method for the signal; allowed values are

OUTPUT:

• the signal-to-noise ratio $$\rho$$

EXAMPLES:

Let us evaluate the SNR of the gravitational signal generated by a 1-solar mass object orbiting at the ISCO of Sgr A* observed by LISA during 1 day:

sage: from kerrgeodesic_gw import (signal_to_noise_particle,
....:                              lisa_detector, astro_data)
sage: a, r0 = 0., 6.
sage: theta = pi/2
sage: t_obs = 24*3600  # 1 day in seconds
sage: BH_time_scale = astro_data.SgrA_mass_s  # Sgr A* mass in seconds
sage: psd = lisa_detector.power_spectral_density_RCLfit
sage: mu_ov_r = astro_data.Msol_m / astro_data.dSgrA  # mu/r
sage: signal_to_noise_particle(a, r0, theta, psd, t_obs,  # tol 1.0e-13
....:                          BH_time_scale, scale=mu_ov_r)
7565.450402023329


sage: signal_to_noise_particle(a, r0, theta, psd, t_obs,  # tol 1.0e-13
....:                          BH_time_scale, scale=mu_ov_r,
5230.216539094391


Using the 1.5-PN approximation (m_max has to be at most 5):

sage: signal_to_noise_particle(a, r0, theta, psd, t_obs,  # tol 1.0e-13
....:                          BH_time_scale, scale=mu_ov_r,
....:                          approximation='1.5PN', m_max=5)
7601.135299165022


For large values of $$r_0$$, the 1.5-PN approximation and the quadrupole one converge:

sage: r0 = 100
sage: signal_to_noise_particle(a, r0, theta, psd, t_obs,  # tol 1.0e-13
....:                          BH_time_scale, scale=mu_ov_r,
0.003053021617751869
sage: signal_to_noise_particle(a, r0, theta, psd, t_obs,  # tol 1.0e-13
....:                          BH_time_scale, scale=mu_ov_r,
....:                          approximation='1.5PN')
0.003144006483338098

sage: r0 = 1000
sage: signal_to_noise_particle(a, r0, theta, psd, t_obs,  # tol 1.0e-13
....:                          BH_time_scale, scale=mu_ov_r,