{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"![logo](https://joss.theoj.org/logo_large.jpg)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"# Intro\nSo, I came across a paper called [Riroriro: Simulating gravitational waves and evaluatingtheir detectability in Python](https://joss.theoj.org/papers/10.21105/joss.02968) which describes a Python package for generating and analyzing simulated gravitational waves. The data set we have in this competition is dominated by detector noise. By exploring noise-free simulated waveforms we can gain a much better understanding of the properties of the signals. The host probably used a very similar (or the same?) numerical library to generate the gravitational waves in the dataset.  \n\nFirst, install Riroriro:","metadata":{}},{"cell_type":"code","source":"!pip install riroriro","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-08-31T19:38:25.525001Z","iopub.execute_input":"2021-08-31T19:38:25.525410Z","iopub.status.idle":"2021-08-31T19:39:01.343539Z","shell.execute_reply.started":"2021-08-31T19:38:25.525325Z","shell.execute_reply":"2021-08-31T19:39:01.342414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The source code can be found on [GitHub](https://github.com/wvanzeist/riroriro).","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport riroriro.inspiralfuns as ins\nimport riroriro.mergerfirstfuns as me1\nimport riroriro.matchingfuns as mat\nimport riroriro.mergersecondfuns as me2\nimport librosa\nimport librosa.display\nimport math\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2021-08-31T19:51:44.277084Z","iopub.execute_input":"2021-08-31T19:51:44.277493Z","iopub.status.idle":"2021-08-31T19:51:46.891927Z","shell.execute_reply.started":"2021-08-31T19:51:44.277441Z","shell.execute_reply":"2021-08-31T19:51:46.891056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generate gravitational waves\nRiroriro comes with a [nice tutorial](https://github.com/wvanzeist/riroriro_tutorials/blob/main/example_GW.ipynb) for how to generate GW signals, and below we simply copy all that code (except for size reduction) into a function:","metadata":{}},{"cell_type":"code","source":"# Parameters:\n# logMC: system mass (0.0-2.0)\n# q: mass ratio (0.1-1.0)\n# D: distance (Mpc)\n# merger_type: 'BH'=binary black hole merger, 'NS'=binary neutron star merger\n# flow: low frequency (Hz) \ndef gen_gw(logMc=1.4, q=0.8, D=100.0, flow=10.0, merger_type='BH'):\n    M, eta = ins.get_M_and_eta(logMc=logMc,q=q)\n    start_x = ins.startx(M,flow)\n    end_x = ins.endx(eta,merger_type)\n    x, xtimes, dt = ins.PN_parameter_integration(start_x,end_x,M,eta)\n    realtimes = ins.inspiral_time_conversion(xtimes,M)\n    i_phase, omega, freq = ins.inspiral_phase_freq_integration(x,dt,M)\n    r, rdot = ins.radius_calculation(x,M,eta)\n    A1, A2 = ins.a1_a2_calculation(r,rdot,omega,D,M,eta)\n    i_Aorth, i_Adiag = ins.inspiral_strain_polarisations(A1,A2,i_phase)\n    i_amp = ins.inspiral_strain_amplitude(i_Aorth,i_Adiag)\n    i_time = realtimes\n    i_omega = omega\n    sfin, wqnm = me1.quasi_normal_modes(eta)\n    alpha, b, C, kappa = me1.gIRS_coefficients(eta,sfin)\n    fhat, m_omega = me1.merger_freq_calculation(wqnm,b,C,kappa)\n    fhatdot = me1.fhat_differentiation(fhat)\n    m_time = me1.merger_time_conversion(M)\n    min_switch_ind = mat.min_switch_ind_finder(i_time,i_omega,m_time,m_omega)\n    final_i_index = mat.final_i_index_finder(min_switch_ind,i_omega,m_omega)\n    time_offset = mat.time_offset_finder(min_switch_ind,final_i_index,i_time,m_time)\n    i_m_time, i_m_omega = mat.time_frequency_stitching(min_switch_ind,final_i_index,time_offset,i_time,i_omega,m_time,m_omega)\n    i_m_freq = mat.frequency_SI_units(i_m_omega,M)\n    m_phase = me2.merger_phase_calculation(min_switch_ind,final_i_index,i_phase,m_omega)\n    i_m_phase = me2.phase_stitching(final_i_index,i_phase,m_phase)\n    m_amp = me2.merger_strain_amplitude(min_switch_ind,final_i_index,alpha,i_amp,m_omega,fhat,fhatdot)\n    i_m_amp = me2.amplitude_stitching(final_i_index,i_amp,m_amp)\n    m_Aorth, m_Adiag = me2.merger_polarisations(final_i_index,m_amp,m_phase,i_Aorth)\n    i_m_Aorth, i_m_Adiag = me2.polarisation_stitching(final_i_index,i_Aorth,i_Adiag,m_Aorth,m_Adiag)\n    return np.array(i_m_time), np.array(i_m_Aorth), np.array(i_m_Adiag), np.array(i_m_freq)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-09T15:28:25.462113Z","iopub.execute_input":"2021-08-09T15:28:25.462503Z","iopub.status.idle":"2021-08-09T15:28:25.474255Z","shell.execute_reply.started":"2021-08-09T15:28:25.462468Z","shell.execute_reply":"2021-08-09T15:28:25.473236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The function returns two waves that represent orthogonal/diagonal waves. The output timescale that is returned is non-linear, so to convert these signals into uniform sampled signals as in the dataset, we need to resample. The function below will resample the gravitational wave signals to 2048Hz. It is crude though, based on nearest sample, but good enough for studying spectrums. Interpolation would be more proper.","metadata":{}},{"cell_type":"code","source":"SR = 2048 # target sample rate (Hz)\n# Parameters:\n# dt: time series\n# amp: amplitude signal\n# seg: output sequence length (seconds)\ndef resample(dt, amp, seg=2.0):\n    end = dt[-1]\n    start = end - seg\n    d = np.zeros(int(SR*seg))\n    for i in range((int(SR*seg))):\n        t = start + i/SR\n        d[i] = amp[np.where(dt == dt[np.abs(dt-t).argmin()])[0][0]]\n    return d","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-09T15:28:26.767638Z","iopub.execute_input":"2021-08-09T15:28:26.767959Z","iopub.status.idle":"2021-08-09T15:28:26.773842Z","shell.execute_reply.started":"2021-08-09T15:28:26.767931Z","shell.execute_reply":"2021-08-09T15:28:26.772718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Another helper function for plotting the last part of the GW signal (containing the chirp):","metadata":{}},{"cell_type":"code","source":"def plot_sig(dt, sig1, sig2=None, seg=2.0):\n    end = dt[-1]\n    start = end - seg\n    plt.figure(1)\n    plt.plot(dt, sig1)\n    peak = np.max(np.abs(sig1))\n    plt.axis([start,end,np.min(sig1)-peak/10,np.max(sig1)+peak/10])\n    if sig2 is not None:\n        plt.plot(dt, sig2)\n    plt.xlabel('Time (s)')\n    plt.ylabel('Strain amplitude')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-09T15:28:30.040386Z","iopub.execute_input":"2021-08-09T15:28:30.040715Z","iopub.status.idle":"2021-08-09T15:28:30.048487Z","shell.execute_reply.started":"2021-08-09T15:28:30.040684Z","shell.execute_reply":"2021-08-09T15:28:30.04769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test the signal generation\nNow, let's test the code by emulating the first GW detected (GW150914).","metadata":{}},{"cell_type":"code","source":"m_time, m_Aorth, m_Adiag, m_freq = gen_gw(logMc=1.4, q=0.2)","metadata":{"execution":{"iopub.status.busy":"2021-08-09T15:28:32.856242Z","iopub.execute_input":"2021-08-09T15:28:32.856559Z","iopub.status.idle":"2021-08-09T15:30:33.759485Z","shell.execute_reply.started":"2021-08-09T15:28:32.85653Z","shell.execute_reply":"2021-08-09T15:30:33.758486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plot the signals, first the last two seconds, then the last 0.1s (ringdown part).","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(16,8))\nplt.subplot(2, 1, 1)\nplot_sig(m_time, m_Aorth, m_Adiag, seg=2)\nplt.subplot(2, 1, 2)\nplot_sig(m_time, m_Aorth, m_Adiag, seg=.1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-09T15:39:51.637507Z","iopub.execute_input":"2021-08-09T15:39:51.637916Z","iopub.status.idle":"2021-08-09T15:39:52.590038Z","shell.execute_reply.started":"2021-08-09T15:39:51.637875Z","shell.execute_reply":"2021-08-09T15:39:52.586377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, resample the signal to 2048Hz (only the orthogonal part for simplicity):","metadata":{}},{"cell_type":"code","source":"d1 = resample(m_time, m_Aorth, 2.0)","metadata":{"execution":{"iopub.status.busy":"2021-08-09T15:39:57.08652Z","iopub.execute_input":"2021-08-09T15:39:57.086841Z","iopub.status.idle":"2021-08-09T15:40:22.651771Z","shell.execute_reply.started":"2021-08-09T15:39:57.086813Z","shell.execute_reply":"2021-08-09T15:40:22.650768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compare with the original:","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(16,16))\nplt.subplot(2, 2, 1)\nplot_sig(m_time, m_Aorth, seg=2)\nplt.title('Original')\nplt.subplot(2, 2, 2)\nplot_sig(m_time, m_Aorth, seg=.1)\nplt.title('Original (zoomed)')\nplt.subplot(2, 2, 3)\nplt.plot(d1)\nplt.title('Resampled to 2048Hz')\nplt.subplot(2, 2, 4)\nplt.plot(d1[-205:])\nplt.title('Resampled to 2048Hz (zoomed)');","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-09T15:43:01.262604Z","iopub.execute_input":"2021-08-09T15:43:01.262943Z","iopub.status.idle":"2021-08-09T15:43:02.082651Z","shell.execute_reply.started":"2021-08-09T15:43:01.262914Z","shell.execute_reply":"2021-08-09T15:43:02.081904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Yes, we are happy with that! The simulated wave data also has a frequency vector. Let's visualize that:","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(12,8))\nplt.plot(m_time, m_freq, label=\"Min: {}Hz, Max: {}Hz\".format(int(np.min(m_freq)), int(np.max(m_freq))))\npeak = np.max(np.abs(m_freq))\nplt.axis([m_time[-1] - 2.0 if m_time[-1] >= 2.0 else m_time[0], m_time[-1], 0 , np.max(m_freq)+peak/10])\nax.legend()\nplt.xlabel('Time (s)')\nplt.ylabel('Frequency (Hz)');","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-09T15:43:07.813436Z","iopub.execute_input":"2021-08-09T15:43:07.813756Z","iopub.status.idle":"2021-08-09T15:43:08.873996Z","shell.execute_reply.started":"2021-08-09T15:43:07.813725Z","shell.execute_reply":"2021-08-09T15:43:08.873312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So the gravitational wave signal starts around 10Hz and rises rapidly to 181Hz at the end. The start frequency is actually defined by the 'flow' parameter.","metadata":{}},{"cell_type":"markdown","source":"# Spectrum\nWe can now do a constant-Q transform (here using librosa) to visualize the now familiar chirp part of the gravitational wave signal.\n","metadata":{}},{"cell_type":"code","source":"hop_length = 64\nC = np.abs(librosa.cqt(d1/np.max(d1), sr=SR, hop_length=hop_length, fmin=8, filter_scale=0.8, bins_per_octave=12))\nfig, ax = plt.subplots(figsize=(10,10))\nimg = librosa.display.specshow(librosa.amplitude_to_db(C, ref=np.max),\n                               sr=SR*2, hop_length=hop_length, bins_per_octave=12, ax=ax)\nax.set_title('Constant-Q power spectrum');","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-09T15:43:23.258618Z","iopub.execute_input":"2021-08-09T15:43:23.259108Z","iopub.status.idle":"2021-08-09T15:43:24.269381Z","shell.execute_reply.started":"2021-08-09T15:43:23.259074Z","shell.execute_reply":"2021-08-09T15:43:24.268686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Chirp frequency vs. system mass\nNow let's see how the system mass and mass ratio parameters effect the frequency content of the GW signal. We keep one parameter constant (middle value) while varying the other. Starting with system mass (logMc parameter).","metadata":{}},{"cell_type":"code","source":"hop_length = 64\nmass = [0.1, 0.7, 2.0]\nfig = plt.figure(figsize=(20,15))\nq = [0.45]\nfor m in range(len(mass)):\n    m_time, m_Aorth, _, m_freq = gen_gw(logMc=mass[m], q=q[0])\n    rd = resample(m_time, m_Aorth, 2.0)\n    # time series\n    ax = plt.subplot(len(mass), 4, 1+m*4)\n    plt.plot(rd)\n    plt.title('Signal (logMc={}, q={})'.format(mass[m], q[0]))\n    # zoomed times series (chirp)\n    ax = plt.subplot(len(mass), 4, 2+m*4)\n    plt.plot(rd[-205:])\n    plt.title('Signal chirp (zoomed)')\n    # frequency content (last 2s)\n    ax = plt.subplot(len(mass), 4, 3+m*4)\n    plt.plot(m_time, m_freq, label=\"Min: {}Hz, Max: {}Hz\".format(int(np.min(m_freq)), int(np.max(m_freq))))\n    peak = np.max(np.abs(m_freq))\n    plt.axis([m_time[-1] - 2.0 if m_time[-1] >= 2.0 else m_time[0], m_time[-1], 0 , np.max(m_freq)+peak/10])\n    ax.legend()\n    plt.xlabel('Time (s)')\n    plt.ylabel('Frequency (Hz)');\n    plt.title('Frequency content (last 2s)')\n    # Q-Transform\n    ax = plt.subplot(len(mass), 4, 4+m*4)\n    C = np.abs(librosa.cqt(rd/np.max(rd), sr=SR, hop_length=hop_length, fmin=8, filter_scale=0.8, bins_per_octave=12))\n    img = librosa.display.specshow(librosa.amplitude_to_db(C, ref=np.max),\n                                   sr=SR*2, hop_length=hop_length, bins_per_octave=12, ax=ax)\n    ax.set_title('Constant-Q power spectrum');","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-22T11:28:14.525785Z","iopub.execute_input":"2021-07-22T11:28:14.526313Z","iopub.status.idle":"2021-07-22T11:28:18.923718Z","shell.execute_reply.started":"2021-07-22T11:28:14.52628Z","shell.execute_reply":"2021-07-22T11:28:18.92236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Chirp frequency vs. mass ratio\nNext mass ratio (q parameter).","metadata":{}},{"cell_type":"code","source":"hop_length = 64\nmass = [1.0]\nfig = plt.figure(figsize=(20,15))\nq = [0.1, 0.45, 1.0]\nfor m in range(len(q)):\n    m_time, m_Aorth, _, m_freq = gen_gw(logMc=mass[0], q=q[m])\n    rd = resample(m_time, m_Aorth, 2.0)\n    # time series\n    ax = plt.subplot(len(q), 4, 1+m*4)\n    plt.plot(rd)\n    plt.title('Signal (logMc={}, q={})'.format(mass[0], q[m]))\n    # zoomed times series (chirp)\n    ax = plt.subplot(len(q), 4, 2+m*4)\n    plt.plot(rd[-205:])\n    plt.title('Signal chirp (zoomed)')\n    # frequency content (last 2s)\n    ax = plt.subplot(len(q), 4, 3+m*4)\n    plt.plot(m_time, m_freq, label=\"Min: {}Hz, Max: {}Hz\".format(int(np.min(m_freq)), int(np.max(m_freq))))\n    peak = np.max(np.abs(m_freq))\n    plt.axis([m_time[-1] - 2.0 if m_time[-1] >= 2.0 else m_time[0], m_time[-1], 0 , np.max(m_freq)+peak/10])\n    ax.legend()\n    plt.xlabel('Time (s)')\n    plt.ylabel('Frequency (Hz)');\n    plt.title('Frequency content (last 2s)')\n    # Q-Transform\n    ax = plt.subplot(len(q), 4, 4+m*4)\n    C = np.abs(librosa.cqt(rd/np.max(rd), sr=SR, hop_length=hop_length, fmin=8, filter_scale=0.8, bins_per_octave=12))\n    img = librosa.display.specshow(librosa.amplitude_to_db(C, ref=np.max),\n                                   sr=SR*2, hop_length=hop_length, bins_per_octave=12, ax=ax)\n    ax.set_title('Constant-Q power spectrum');","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-22T11:47:36.148353Z","iopub.execute_input":"2021-07-22T11:47:36.148942Z","iopub.status.idle":"2021-07-22T12:05:15.171639Z","shell.execute_reply.started":"2021-07-22T11:47:36.148892Z","shell.execute_reply":"2021-07-22T12:05:15.170392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Amplitude vs distance\nDo gravitational waves follow the general Inverse-square law (double the distance and get 0.25 of the amplitude)?","metadata":{}},{"cell_type":"code","source":"hop_length = 64\n\nfig = plt.figure(figsize=(20,15))\ndist = [100., 200., 400.]\nfor m in range(len(dist)):\n    m_time, m_Aorth, _, m_freq = gen_gw(logMc=1.4, q=0.2, D=dist[m])\n    rd = resample(m_time, m_Aorth, 2.0)\n    # time series\n    ax = plt.subplot(len(dist), 3, 1+m*3)\n    plt.plot(rd)\n    plt.title('Signal (D={} Mpc)'.format(int(dist[m])))\n    # zoomed times series (chirp)\n    ax = plt.subplot(len(dist), 3, 2+m*3)\n    plt.plot(rd[-205:])\n    plt.title('Signal chirp (zoomed)')\n    # Q-Transform\n    ax = plt.subplot(len(dist), 3, 3+m*3)\n    if m == 0:\n        smax = np.max(rd)\n    C = np.abs(librosa.cqt(rd/smax, sr=SR, hop_length=hop_length, fmin=8, filter_scale=0.8, bins_per_octave=12))\n    if m == 0:\n        Cmax = np.max(C)\n    img = librosa.display.specshow(librosa.amplitude_to_db(C, ref=Cmax), # was np.max\n                                   sr=SR*2, hop_length=hop_length, bins_per_octave=12, ax=ax)\n    ax.set_title('Constant-Q power spectrum');","metadata":{"execution":{"iopub.status.busy":"2021-08-09T15:48:14.789854Z","iopub.execute_input":"2021-08-09T15:48:14.791551Z","iopub.status.idle":"2021-08-09T15:55:34.652424Z","shell.execute_reply.started":"2021-08-09T15:48:14.79146Z","shell.execute_reply":"2021-08-09T15:55:34.651063Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Surprise - they do not! If the distance is doubled, we get 0.5 of the amplitude! Read more about this [here](https://www.forbes.com/sites/startswithabang/2019/03/02/ask-ethan-why-dont-gravitational-waves-get-weaker-like-the-gravitational-force-does/?sh=260758c72f58).","metadata":{}},{"cell_type":"markdown","source":"# Add noise to the signal\nNext we want to see how the signal looks after adding detector noise. Half the training data is simulated detector noise without signal (target=0), so all we have to do is to add our generated signal to one of those numpy files... Let's pick the first file with target=0:","metadata":{}},{"cell_type":"code","source":"noise = np.load('../input/g2net-gravitational-wave-detection/train/0/0/0/00001f4945.npy')\nplt.plot(noise[1,:]);","metadata":{"execution":{"iopub.status.busy":"2021-07-22T18:24:43.115729Z","iopub.execute_input":"2021-07-22T18:24:43.116493Z","iopub.status.idle":"2021-07-22T18:24:43.2878Z","shell.execute_reply.started":"2021-07-22T18:24:43.11642Z","shell.execute_reply":"2021-07-22T18:24:43.286883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then we add the GW150914 signal we created above to the noise data with different signal to noise ratios, and check how well the constant-Q transform captures the chirp:","metadata":{}},{"cell_type":"code","source":"gain = [0.0, 1.0, 1/2, 1/8, 1/16, 1/32]\nfig = plt.figure(figsize=(15,25))\n\ndef add_noise(gain, sig):\n    nd = noise[1,:] + gain*sig\n    return nd\n\nfor i in range(len(gain)):\n    nd = add_noise(gain[i], d1)\n    # time signal\n    ax = plt.subplot(len(gain), 2, 1+i*2)\n    plt.plot(nd)\n    if i == 0:\n        plt.title('Noise only')\n    else:\n        plt.title('Noise + signal, gain={}'.format(gain[i]))\n    # constant-Q transform\n    ax = plt.subplot(len(gain), 2, 2+i*2)\n    C = np.abs(librosa.cqt(nd/np.max(nd), sr=SR, hop_length=hop_length, fmin=8, filter_scale=0.8, bins_per_octave=12))\n    img = librosa.display.specshow(librosa.amplitude_to_db(C, ref=np.max),\n                                   sr=SR*2, hop_length=hop_length, bins_per_octave=12, ax=ax)\n    if i == 0:\n        ax.set_title('Constant-Q power spectrum')\n    else:\n        ax.set_title('Constant-Q power spectrum (signal @{}dB)'.format(int(10*math.log10(gain[i]))));","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-22T18:24:55.125673Z","iopub.execute_input":"2021-07-22T18:24:55.126242Z","iopub.status.idle":"2021-07-22T18:24:56.752016Z","shell.execute_reply.started":"2021-07-22T18:24:55.126187Z","shell.execute_reply":"2021-07-22T18:24:56.75088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"At -12dB the chirp is just barely visible!","metadata":{}},{"cell_type":"markdown","source":"# Signal position within the 2s time series\nIn the plots above, the GW signal ends at the very end of the 2s window. In the dataset the signals typically end somewhere in the last 0.5s part. This can be observed in Q-transforms of strong signals within the dataset. Below a few Q-Transforms are visualized with the signal ending at different times within the 2s window.","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(12,20))\nd2 = np.concatenate((d1, np.zeros(4096))) # pad with 2s of zeros\npos = [0.5, 0.7, 0.8, 0.9, 0.99]\nhop_length = 64\n\nfor i in range(len(pos)):\n    # time series\n    ax = plt.subplot(len(pos), 2, 1+i*2)\n    start = 4096 - int(4096*pos[i])\n    plt.plot(d2[start:start+4096])\n    plt.title('Signal ending at {}s'.format(pos[i]*2))\n    # Q-transform\n    ax = plt.subplot(len(pos), 2, 2+i*2)\n    C = np.abs(librosa.cqt(d2[start:start+4096]/np.max(d2), sr=SR, hop_length=hop_length, fmin=8, filter_scale=0.8, bins_per_octave=12))\n    img = librosa.display.specshow(librosa.amplitude_to_db(C, ref=np.max),\n                                   sr=SR*2, hop_length=hop_length, bins_per_octave=12, ax=ax)\n    ax.set_title('Constant-Q power spectrum')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-23T11:33:08.516213Z","iopub.execute_input":"2021-07-23T11:33:08.51663Z","iopub.status.idle":"2021-07-23T11:33:10.86844Z","shell.execute_reply.started":"2021-07-23T11:33:08.516595Z","shell.execute_reply":"2021-07-23T11:33:10.867629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create noise signals from PSD\nSo, let's say we have a power spectrum density of the detector noise - how can we generate random noise signals from that? The power spectral density is a representation of statistical power distribution across frequency range of a specific signal. The procedure is quite simple:  \n  * Make sure the number of samples in your given PSD is M = N/2 + 1, where N is desired FFT size\n  * Give each spectral component a random phase, uniformly distributed between 0 and 360 degrees (or 0 and 2Pi radians)\n  * Multiply the PSD amplitudes with random numbers with variance 1\n  * Perform inverse FFT to obtain the time series  \n  \nBut first we need to obtain the PSD for the three detectors. There are some unofficial PSDs in the [riroriro tutorial project](https://github.com/wvanzeist/riroriro_tutorials). So we fetch them below. Also we install a package for spectral resampling called [SpectRes](https://github.com/ACCarnall/SpectRes).","metadata":{}},{"cell_type":"code","source":"!git clone https://github.com/wvanzeist/riroriro_tutorials.git\n!pip install spectres","metadata":{"execution":{"iopub.status.busy":"2021-08-31T19:49:27.420176Z","iopub.execute_input":"2021-08-31T19:49:27.420874Z","iopub.status.idle":"2021-08-31T19:50:12.265809Z","shell.execute_reply.started":"2021-08-31T19:49:27.420830Z","shell.execute_reply":"2021-08-31T19:50:12.264548Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a look at one of those PSDs (Livo Livingston):","metadata":{}},{"cell_type":"code","source":"psd_liv = np.genfromtxt('riroriro_tutorials/noise_spectra/o3_h1.txt') #LIGO Livingston\nplt.plot(psd_liv[:,0], np.log10(psd_liv[:,1]));","metadata":{"execution":{"iopub.status.busy":"2021-08-31T19:52:39.041428Z","iopub.execute_input":"2021-08-31T19:52:39.041975Z","iopub.status.idle":"2021-08-31T19:52:39.277008Z","shell.execute_reply.started":"2021-08-31T19:52:39.041943Z","shell.execute_reply":"2021-08-31T19:52:39.276211Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The frequency vector starts at 10Hz and ends at 5kHz. And frequencies are not uniform (are they ever in astrophysics?). So we need to resample to our desired frequency interval of 0-1024Hz. We create a function for all the steps mentioned above:","metadata":{}},{"cell_type":"code","source":"from spectres import spectres\n\n# convert polar coordinates to rectangular format\ndef P2R(A, phi):\n    return A * (np.cos(phi) + np.sin(phi)*1j)\n\n# input is a 2 dimensional power spectrum: frequencies and amplitudes\ndef rand_wave(power_spectrum):\n    regrid = np.arange(0., 1024.5, .5) # note: 2049 length to get 4096 samples from irfft\n    # resample spectrum\n    PDS, _ = spectres(regrid, power_spectrum[:,0], power_spectrum[:,1],\n                      spec_errs=np.zeros(len(power_spectrum[:,0])), fill=0., verbose=False)\n    # add random phase and amplitude\n    ph = np.random.uniform(0.,2*np.pi,len(PDS)) # random phase\n    PDS *= np.random.randn(len(PDS)) # random amplitude\n    Z=P2R(PDS, ph) # polar to rectangular format\n    Z *= 100. # scale factor (to get about the right signal amplitude)\n    return np.fft.irfft(Z) ","metadata":{"execution":{"iopub.status.busy":"2021-08-31T19:56:48.085372Z","iopub.execute_input":"2021-08-31T19:56:48.085733Z","iopub.status.idle":"2021-08-31T19:56:48.096842Z","shell.execute_reply.started":"2021-08-31T19:56:48.085705Z","shell.execute_reply":"2021-08-31T19:56:48.095771Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then let's generate a few signals:","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(20,20))\nfor m in range(16):\n    ax = plt.subplot(4, 4, 1+m)\n    plt.plot(rand_wave(psd_liv));","metadata":{"execution":{"iopub.status.busy":"2021-08-31T20:04:03.154822Z","iopub.execute_input":"2021-08-31T20:04:03.155195Z","iopub.status.idle":"2021-08-31T20:04:07.154534Z","shell.execute_reply.started":"2021-08-31T20:04:03.155167Z","shell.execute_reply":"2021-08-31T20:04:07.153430Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks familiar! But we do not know anything about how the host generated the noise in the data files, so this bit is maybe not very useful for improving scores...","metadata":{}},{"cell_type":"markdown","source":"# Summary","metadata":{}},{"cell_type":"markdown","source":"From the plots above, we can see that the bigger the system mass, the lower the chirp frequency (as expected). And for small system masses, the chirp frequencies can reach many kHz. The mass ratio has less effect on the chirp, but the closer the masses are the higher the frequency. We have also seen how to add noise from data files with target=0, or even by creating new noise data form a PSD.","metadata":{}}]}