Skip to content
Snippets Groups Projects
Commit 6d1a615e authored by Wei Changfeng's avatar Wei Changfeng :cancer:
Browse files

Delete whiten_signal.ipynb

parent bc8a7d06
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
This is a notebook to whiten signal genreated from "h_m16_L0.18_l2m2_r300.dat".
%% Cell type:code id: tags:
``` python
import json
import numpy as np
import matplotlib.pyplot as plt
```
%% Cell type:code id: tags:
``` python
with open('0.json','r') as f1:
waveform = json.load(f1)
#extract time series at V1,L1,and H1 from waveform file.
data_V1 = waveform['data']['V1']
data_L1 = waveform['data']['L1']
data_H1 = waveform['data']['H1']
times_V1 = waveform['times']['V1']
times_L1 = waveform['times']['L1']
times_H1 = waveform['times']['H1']
sample_rate = 256
dt=1./sample_rate
```
%% Cell type:code id: tags:
``` python
# Plot the waveform of L1.
plt.title('Original signal')
plt.xlabel('Time(s)')
plt.ylabel('Strain')
plt.plot(times_L1,data_L1)
```
%% Output
[<matplotlib.lines.Line2D at 0x7fd08c6a46a0>]
%% Cell type:code id: tags:
``` python
# Load ASD file of L1.
filename = './bilby_gw_detector_noise_curves_aLIGO_ZERO_DET_high_P_asd.txt'
freq = []
asd = []
with open(filename, '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.append(freq_tmp)
asd.append(asd_tmp)
pass
# Transfer data format from list to array.
freq = np.array(freq)
asd = np.array(asd)
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.
Nt = len(data_L1)
freq_rfft = np.fft.rfftfreq(Nt, dt)
asd_interp = np.interp(freq_rfft, freq, asd)
plt.loglog(freq, asd,'-o',label='L1 ASD')
plt.loglog(freq_rfft, asd_interp, 'x',label='Interpolation')
plt.grid('on')
plt.ylabel('ASD (strain/rtHz)')
plt.xlabel('Freq (Hz)')
plt.legend(loc='upper center')
# The frequency range is set here, which comes from the calculation of the LIGO detection range.
# The limit of the seismic isolation systems designed by aLIGO is 10Hz, which is actually not good for low frequency noise below 20Hz.
# According to Nyquist–Shannon sampling theorem, for half of the sampling rate(128Hz) and above,the signal is meaningless.
f_min = 20.
f_max = 128.
plt.axis([f_min, f_max, 1e-24, 1e-19])
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
# This step is not necessary. Just to show a potential issue: interpolation of ASD has some points out of the curve in left side,
# I am not sure if this has a big impact on the calculation of the whitening process.
plt.loglog(freq, asd,'-o',label='L1 ASD')
plt.loglog(freq_rfft, asd_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 0x7fd05c141eb8>
%% Cell type:code id: tags:
``` python
whitened_signal_L1 = whiten_signal(data_L1,asd_interp)
```
%% Cell type:code id: tags:
``` python
plt.plot(times_L1,whitened_signal_L1)
plt.title('Whitened signal')
plt.xlabel('Time(s)')
plt.ylabel('Strain')
```
%% Output
Text(0, 0.5, 'Strain')
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