Skip to content
Snippets Groups Projects
Commit c4074c8b authored by BNU computer's avatar BNU computer
Browse files

Upload New File

parent 71832034
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
This is a notebook to generate Korean BBH waveforms for VItamin box to test.
We produce random injection parameters first, and Minke will give waveforms. Then we whiten the waveforms and add Gaussian noise.
%% Cell type:code id: tags:
``` python
import minke.sources
import lalsimulation
import lal
import numpy as np
import minke.distribution
import matplotlib.pyplot as plt
```
%% Cell type:code id: tags:
``` python
nr_path = '../h_psi4/{}'
waveforms_merge = {
1: "h_m1_L0.9_l2m2_r300.dat",
2: "h_m2_L0.87_l2m2_r300.dat",
4: "h_m4_L0.5_l2m2_r300.dat",
8: "h_m8_L0.35_l2m2_r280.dat",
16: "h_m16_L0.18_l2m2_r300.dat"
}
nr_waveform = waveforms_merge[1]
duration = 1
sampling_frequency = 256
epoch = "0.25"
geocent_time = 0.25
num_samples = 4
locations = minke.distribution.uniform_sky(num_samples)
distances = np.random.uniform(1000, 3000, num_samples)
total_mass = 115
```
%% Cell type:code id: tags:
``` python
locations
```
%% Output
(array([3.01617093, 6.04781136, 0.75714575, 4.06214126]),
array([-0.45504638, -0.5027633 , 1.13408584, 0.83943672]),
array([4.16975228, 5.06352657, 0.0317775 , 3.93614481]))
%% Cell type:code id: tags:
``` python
distances
```
%% Output
array([1525.40619263, 1155.95688485, 2108.11537646, 2559.92018662])
%% Cell type:code id: tags:
``` python
def generate_for_detector(source, ifos, sampling_frequency, epoch, distance, total_mass, ra, dec, psi):
"""
Generate an injection for a given waveform.
Parameters
----------
source : ``minke.Source`` object
The source which should generate the waveform.
ifos : list
A list of interferometer initialisms, e.g. ['L1', 'H1']
sampling_frequency : int
The sampling_frequency in hertz for the generated signal.
epoch : str
The epoch (start time) for the injection.
Note that this should be provided as a string to prevent overflows.
distance : float
The distance for the injection, in megaparsecs.
total_mass : float
The total mass for the injected signal.
ra : float
The right ascension of the source, in radians.
dec : float
The declination of the source, in radians.
psi : float
The polarisation angle of the source, in radians.
Notes
-----
This is very much a draft implementation, and there are a number of
things which should be tidied-up before this is moved into Minke,
including bundling total mass with the source.
"""
nr_waveform = source.datafile
data = {}
data['data'] = {}
data['times'] = {}
data['meta'] = {}
data['epoch'] = epoch
data['meta']['ra'] = ra
data['meta']['dec'] = dec
data['meta']['psi'] = psi
data['meta']['distance'] = distance
data['meta']['total_mass'] = total_mass
data['meta']['sampling_frequency'] = sampling_frequency
data['meta']['waveform'] = nr_waveform
for ifo in ifos:
det = lalsimulation.DetectorPrefixToLALDetector(ifo)
hp, hx = source._generate(half=True, epoch=epoch, rate=sampling_frequency)[:2]
h_tot = lalsimulation.SimDetectorStrainREAL8TimeSeries(hp, hx, ra, dec, psi, det)
data['data'][ifo] = h_tot.data.data.tolist()
data['times'][ifo] = np.linspace(0, len(h_tot.data.data)*h_tot.deltaT, len(h_tot.data.data)).tolist()
return data
```
%% Cell type:code id: tags:
``` python
waveforms = []
ifos = ['V1', 'L1', 'H1']
for distance, location in zip(distances, np.vstack(locations).T):
hyper_object = minke.sources.Hyperbolic(
datafile = nr_path.format(nr_waveform),
total_mass = total_mass,
distance = distance,
)
hyper_object.datafile = nr_waveform
data = generate_for_detector(hyper_object, ifos, sampling_frequency, epoch, distance, total_mass, *location)
waveforms.append(data)
```
%% Cell type:code id: tags:
``` python
# We need to pad the beginning of the waveform with zeros
def pad_waveform(ifos,num_samples,waveforms):
for i in range(num_samples):
for ifo in ifos:
for j in range(256-len(waveforms[i]['data'][ifo])):
waveforms[i]['data'][ifo] = waveforms[i]['data'][ifo] + [0]
return waveforms
```
%% Cell type:code id: tags:
``` python
waveforms = pad_waveform(ifos,num_samples,waveforms)
```
%% Cell type:code id: tags:
``` python
times = np.linspace(0.25,1.25,256)
plt.plot(times,waveforms[0]['data']['L1'])
plt.title('L1')
plt.xlabel('times')
plt.ylabel('waveform')
```
%% Output
Text(0, 0.5, 'waveform')
%% Cell type:code id: tags:
``` python
# Load ASD file of Ligo and Virgo.
asd_file_Ligo = './aLIGO_ZERO_DET_high_P_asd.txt'
asd_file_Virgo = './AdV_asd.txt'
freq_Ligo = []
asd_Ligo = []
freq_Virgo = []
asd_Virgo = []
with open(asd_file_Ligo, 'r') as f1:
while True:
lines = f1.readline()
if not lines:
break
pass
freq_tmp, asd_tmp = [float(i) for i in lines.split()]
freq_Ligo.append(freq_tmp)
asd_Ligo.append(asd_tmp)
pass
# Transfer data format from list to array.
freq_Ligo = np.array(freq_Ligo)
asd_Ligo = np.array(asd_Ligo)
pass
with open(asd_file_Virgo, 'r') as f2:
while True:
lines = f2.readline()
if not lines:
break
pass
freq_tmp, asd_tmp = [float(i) for i in lines.split()]
freq_Virgo.append(freq_tmp)
asd_Virgo.append(asd_tmp)
pass
# Transfer data format from list to array.
freq_Virgo = np.array(freq_Virgo)
asd_Virgo= np.array(asd_Virgo)
pass
```
%% Cell type:code id: tags:
``` python
def whiten_signal(signal, asd):
signal_rfft = np.fft.rfft(signal)
whitened_signal_rfft = signal_rfft / asd
whitened_signal = np.fft.irfft(whitened_signal_rfft, n=Nt)
return whitened_signal
```
%% Cell type:code id: tags:
``` python
# Interpolate asd data linearly.
# Vitamin box requires sampling frequency = 256Hz, duration = 1 second.
Nt = 256
dt = duration/sampling_frequency
freq_Ligo_rfft = np.fft.rfftfreq(Nt, dt)
asd_Ligo_interp = np.interp(freq_Ligo_rfft, freq_Ligo, asd_Ligo)
plt.loglog(freq_Ligo, asd_Ligo,'-o',label='Ligo ASD')
plt.loglog(freq_Ligo_rfft, asd_Ligo_interp, 'x',label='Interpolation')
freq_Virgo_rfft = np.fft.rfftfreq(Nt, dt)
asd_Virgo_interp = np.interp(freq_Virgo_rfft, freq_Virgo, asd_Virgo)
plt.loglog(freq_Virgo, asd_Virgo,'-o',label='Virgo ASD')
plt.loglog(freq_Virgo_rfft, asd_Virgo_interp, 'x',label='Interpolation')
plt.grid('on')
plt.ylabel('ASD (strain/rtHz)')
plt.xlabel('Freq (Hz)')
plt.legend(loc='upper center')
```
%% Output
<matplotlib.legend.Legend at 0x7fefbea6e850>
%% Cell type:code id: tags:
``` python
asd_interp = {}
asd_interp['H1'] = asd_Ligo_interp
asd_interp['L1'] = asd_Ligo_interp
asd_interp['V1'] = asd_Virgo_interp
for i in range(num_samples):
for ifo in ifos:
waveforms[i]['data'][ifo] = whiten_signal(waveforms[i]['data'][ifo],asd_interp[ifo])
waveforms[i]['data'][ifo] += 1.0*np.random.randn(256)
```
%% Cell type:code id: tags:
``` python
times = np.linspace(geocent_time, geocent_time+1.0, 256)
plt.plot(times,waveforms[0]['data']['L1'])
plt.title('L1')
plt.xlabel('times')
plt.ylabel('waveform')
```
%% Output
Text(0, 0.5, 'waveform')
%% Cell type:code id: tags:
``` python
y_data_noisy_array=np.zeros((num_samples,3,256))
for i in range(num_samples):
y_data_noisy_array[i][0,:] = waveforms[i]['data']['H1']
y_data_noisy_array[i][1,:] = waveforms[i]['data']['L1']
y_data_noisy_array[i][2,:] = waveforms[i]['data']['V1']
```
%% Cell type:code id: tags:
``` python
import h5py
x_data_list = []
for i in range(num_samples):
with h5py.File(f'data_{i}.h5py','w') as f:
x_data_list = [57.5,57.5,distances[i],geocent_time,locations[2][i],locations[0][i],locations[1][i]]
x_data=np.array(x_data_list)
d1 = f.create_dataset('x_data',data = x_data)
d2 = f.create_dataset('y_data_noisy',data = y_data_noisy_array[i])
```
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment