{"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":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-04T21:46:01.376949Z","iopub.execute_input":"2022-10-04T21:46:01.377516Z","iopub.status.idle":"2022-10-04T21:46:04.093649Z","shell.execute_reply.started":"2022-10-04T21:46:01.3774Z","shell.execute_reply":"2022-10-04T21:46:04.089851Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Riroriro: Simulating gravitational waves and evaluating their detectability in Python\n\nAuthors: Wouter G. J. van Zeist, Héloïse F. Stevance, J. J. Eldridge\n\nhttps://doi.org/10.48550/arXiv.2103.06943\n\n\"Riroriro is a Python package to simulate the gravitational waveforms of binary mergers of black holes and/or neutron stars, and calculate several properties of these mergers and waveforms, specifically relating to their observability by gravitational wave detectors.\"\n\n\"The gravitational waveform simulation of Riroriro is based upon the methods of Buskirk and Babiuc-Hamilton (2019), a paper which describes a computational implementation of an earlier theoretical gravitational waveform model by Huerta et al. (2017), using post-Newtonian expansions and an approximation called the implicit rotating source to simplify the Einstein field equations and simulate gravitational waves. Riroriro's calculation of signal-to-noise ratios (SNR) of gravitational wave events is based on the methods of Barrett et al. (2018), with the simpler gravitational wave model Findchirp (Allen et al. (2012)) being used for comparison and calibration in these calculations.\"\n\nhttps://arxiv.org/abs/2103.06943","metadata":{}},{"cell_type":"markdown","source":"#Wouter van Zeist https://github.com/wvanzeist/riroriro\n\nCitation:\n@ARTICLE{2021JOSS....6.2968V,\n       author = {{van Zeist}, Wouter G.~J. and {Stevance}, H{\\'e}lo{\\\"i}se F. and {Eldridge}, J.~J.},\n        title = \"{Riroriro: Simulating gravitational waves and evaluating their detectability in Python}\",\n      journal = {The Journal of Open Source Software},\n     keywords = {Python, neutron stars, astronomy, gravitational waves, black holes, General Relativity and Quantum Cosmology, Astrophysics - High Energy Astrophysical Phenomena},\n         year = 2021,\n        month = mar,\n       volume = {6},\n       number = {59},\n          eid = {2968},\n        pages = {2968},\n          doi = {10.21105/joss.02968},\narchivePrefix = {arXiv},\n       eprint = {2103.06943},\n primaryClass = {gr-qc},\n       adsurl = {https://ui.adsabs.harvard.edu/abs/2021JOSS....6.2968V},\n      adsnote = {Provided by the SAO/NASA Astrophysics Data System}\n}","metadata":{}},{"cell_type":"markdown","source":"#All script by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nVote Geir's work. The original script.","metadata":{}},{"cell_type":"code","source":"!pip install riroriro","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-04T21:46:30.089094Z","iopub.execute_input":"2022-10-04T21:46:30.089581Z","iopub.status.idle":"2022-10-04T21:47:06.085625Z","shell.execute_reply.started":"2022-10-04T21:46:30.089542Z","shell.execute_reply":"2022-10-04T21:47:06.08433Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\"Riroriro is a set of Python modules containing functions to simulate the gravitational waveforms of mergers of black holes and/or neutron stars, and calculate several properties of these mergers and waveforms, specifically relating to their observability by gravitational wave detectors. Riroriro combines areas covered by previous gravitational wave models (such as gravitational wave simulation, SNR calculation, horizon distance calculation) into a single package with broader scope and versatility in Python, a programming language that is ubiquitous in astronomy. Aside from being a research tool, Riroriro is also designed to be easy to use and modify, and it can also be used as an educational tool for students learning about gravitational waves.\"\n\n\"The modules “inspiralfuns”, “mergerfirstfuns”, “matchingfuns”, “mergersecondfuns” and “gwexporter”, in that order, can be used to simulate the strain amplitude and frequency of a merger gravitational waveform. The module “snrcalculatorfuns” can compare such a simulated waveform to a detector noise spectrum to calculate a signal-to-noise ratio (SNR) for that signal for that detector. The module “horizondistfuns” calculates the horizon distance of a merger given its waveform, and the module “detectabilityfuns” evaluates the detectability of a merger given its SNR.\"\n\nMore information on the pip installation can be found here: https://pypi.org/project/riroriro/\n\nTutorials for Riroriro can be found here: https://github.com/wvanzeist/riroriro_tutorials\n\nFull documentation of each of the functions of Riroriro can be found here: https://wvanzeist.github.io/\n\nhttps://github.com/wvanzeist/riroriro","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nimport 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":"2022-10-04T21:47:09.21126Z","iopub.execute_input":"2022-10-04T21:47:09.211782Z","iopub.status.idle":"2022-10-04T21:47:11.665829Z","shell.execute_reply.started":"2022-10-04T21:47:09.211721Z","shell.execute_reply":"2022-10-04T21:47:11.664182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\"Riroriro is one of several Python packages associated with BPASS (Binary Population And Spectral Synthesis), a suite of programs that simulates the evolution of a population of binary and single-star systems from a wide range of initial conditions. Each of these associated packages are named after native animals of New Zealand. The riroriro (Gerygone igata, also known as the grey warbler) is a small bird that can be recognised by its distinctive melodious call but is rarely seen, similarly to how black hole binary mergers are detected by their gravitational wave signals rather than visually.\"\n\n\"The central website of BPASS, which also contains links to related programs, can be found here:\" https://bpass.auckland.ac.nz","metadata":{}},{"cell_type":"markdown","source":"#Generate gravitational waves\n\n\nRiroriro comes with a nice tutorial 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":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\n# 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":{"execution":{"iopub.status.busy":"2022-10-04T21:50:19.978431Z","iopub.execute_input":"2022-10-04T21:50:19.97909Z","iopub.status.idle":"2022-10-04T21:50:19.996374Z","shell.execute_reply.started":"2022-10-04T21:50:19.979029Z","shell.execute_reply":"2022-10-04T21:50:19.995123Z"},"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.\"\n\nhttps://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nSR = 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":{"execution":{"iopub.status.busy":"2022-10-04T21:51:34.899576Z","iopub.execute_input":"2022-10-04T21:51:34.900085Z","iopub.status.idle":"2022-10-04T21:51:34.908342Z","shell.execute_reply.started":"2022-10-04T21:51:34.900046Z","shell.execute_reply":"2022-10-04T21:51:34.907122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Helper function for plotting the last part of the GW signal (containing the chirp):","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\ndef 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":{"execution":{"iopub.status.busy":"2022-10-04T21:52:26.245213Z","iopub.execute_input":"2022-10-04T21:52:26.245715Z","iopub.status.idle":"2022-10-04T21:52:26.254066Z","shell.execute_reply.started":"2022-10-04T21:52:26.245669Z","shell.execute_reply":"2022-10-04T21:52:26.253084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Test the signal generation\n\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":"2022-10-04T21:53:00.554286Z","iopub.execute_input":"2022-10-04T21:53:00.554781Z","iopub.status.idle":"2022-10-04T21:54:38.040538Z","shell.execute_reply.started":"2022-10-04T21:53:00.554713Z","shell.execute_reply":"2022-10-04T21:54:38.03898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#How are gravitational waves detected?\n\n\"When a gravitational wave passes by Earth, it squeezes and stretches space. LIGO can detect this squeezing and stretching. Each LIGO observatory has two “arms” that are each more than 2 miles (4 kilometers) long. A passing gravitational wave causes the length of the arms to change slightly. The observatory uses lasers, mirrors, and extremely sensitive instruments to detect these tiny changes.\"\n\nhttps://spaceplace.nasa.gov/gravitational-waves/en/","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nfig = 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":{"execution":{"iopub.status.busy":"2022-10-04T21:54:42.995132Z","iopub.execute_input":"2022-10-04T21:54:42.99566Z","iopub.status.idle":"2022-10-04T21:54:43.954713Z","shell.execute_reply.started":"2022-10-04T21:54:42.995618Z","shell.execute_reply":"2022-10-04T21:54:43.953431Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#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":"2022-10-04T21:55:18.615291Z","iopub.execute_input":"2022-10-04T21:55:18.617713Z","iopub.status.idle":"2022-10-04T21:55:46.058381Z","shell.execute_reply.started":"2022-10-04T21:55:18.617606Z","shell.execute_reply":"2022-10-04T21:55:46.056686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Compare with the original:","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nfig = 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":{"execution":{"iopub.status.busy":"2022-10-04T21:55:54.086436Z","iopub.execute_input":"2022-10-04T21:55:54.086909Z","iopub.status.idle":"2022-10-04T21:55:55.075482Z","shell.execute_reply.started":"2022-10-04T21:55:54.086869Z","shell.execute_reply":"2022-10-04T21:55:55.07458Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#I didn't screw the code! The simulated wave data also has a frequency vector. Let's visualize that:","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nfig, 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":{"execution":{"iopub.status.busy":"2022-10-04T21:56:39.829711Z","iopub.execute_input":"2022-10-04T21:56:39.830963Z","iopub.status.idle":"2022-10-04T21:56:41.200672Z","shell.execute_reply.started":"2022-10-04T21:56:39.830915Z","shell.execute_reply":"2022-10-04T21:56:41.199472Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#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\n\n\"We can now do a constant-Q transform (here using librosa) to visualize the now familiar chirp part of the gravitational wave signal.\"","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nhop_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":{"execution":{"iopub.status.busy":"2022-10-04T21:58:32.002301Z","iopub.execute_input":"2022-10-04T21:58:32.002969Z","iopub.status.idle":"2022-10-04T21:58:33.411472Z","shell.execute_reply.started":"2022-10-04T21:58:32.002922Z","shell.execute_reply":"2022-10-04T21:58:33.409398Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Chirp frequency vs. system mass\n\n\"Now 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":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nhop_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":{"execution":{"iopub.status.busy":"2022-10-04T21:59:47.05643Z","iopub.execute_input":"2022-10-04T21:59:47.056996Z","iopub.status.idle":"2022-10-04T22:09:40.13415Z","shell.execute_reply.started":"2022-10-04T21:59:47.056956Z","shell.execute_reply":"2022-10-04T22:09:40.132633Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Chirp frequency vs. mass ratio\n\nNext mass ratio (q parameter).","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nhop_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":{"execution":{"iopub.status.busy":"2022-10-04T22:09:51.140175Z","iopub.execute_input":"2022-10-04T22:09:51.140884Z","iopub.status.idle":"2022-10-04T22:19:05.659238Z","shell.execute_reply.started":"2022-10-04T22:09:51.140821Z","shell.execute_reply":"2022-10-04T22:19:05.657568Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Amplitude vs distance\n\n\"Do gravitational waves follow the general Inverse-square law (double the distance and get 0.25 of the amplitude)?\"","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nhop_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":"2022-10-04T22:19:21.615166Z","iopub.execute_input":"2022-10-04T22:19:21.615583Z","iopub.status.idle":"2022-10-04T22:25:38.218708Z","shell.execute_reply.started":"2022-10-04T22:19:21.615547Z","shell.execute_reply":"2022-10-04T22:25:38.217448Z"},"_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.","metadata":{}},{"cell_type":"markdown","source":"#Signal position within the 2s time series\n\n\"In 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":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nfig = 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":{"execution":{"iopub.status.busy":"2022-10-04T22:27:59.742662Z","iopub.execute_input":"2022-10-04T22:27:59.743115Z","iopub.status.idle":"2022-10-04T22:28:01.230518Z","shell.execute_reply.started":"2022-10-04T22:27:59.743077Z","shell.execute_reply":"2022-10-04T22:28:01.22967Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Create noise signals from PSD\n\n\"So, 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\n\"Make sure the number of samples in your given PSD is M = N/2 + 1, where N is desired FFT size\nGive each spectral component a random phase, uniformly distributed between 0 and 360 degrees (or 0 and 2Pi radians)\"\n\n\"Multiply the PSD amplitudes with random numbers with variance 1 Perform inverse FFT to obtain the time series. But first we need to obtain the PSD for the three detectors. There are some unofficial PSDs in the riroriro tutorial project. So we fetch them below. Also we install a package for spectral resampling called SpectRes.\"\n\nhttps://github.com/wvanzeist/riroriro_tutorials\n\nhttps://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":"2022-10-04T22:30:31.897877Z","iopub.execute_input":"2022-10-04T22:30:31.899169Z","iopub.status.idle":"2022-10-04T22:31:04.676001Z","shell.execute_reply.started":"2022-10-04T22:30:31.899121Z","shell.execute_reply":"2022-10-04T22:31:04.674229Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#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":"2022-10-04T22:31:51.556441Z","iopub.execute_input":"2022-10-04T22:31:51.557071Z","iopub.status.idle":"2022-10-04T22:31:51.812566Z","shell.execute_reply.started":"2022-10-04T22:31:51.557018Z","shell.execute_reply":"2022-10-04T22:31:51.811055Z"},"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":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nfrom 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":"2022-10-04T22:32:41.641852Z","iopub.execute_input":"2022-10-04T22:32:41.642371Z","iopub.status.idle":"2022-10-04T22:32:41.658268Z","shell.execute_reply.started":"2022-10-04T22:32:41.642331Z","shell.execute_reply":"2022-10-04T22:32:41.65678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Generate a few signals:","metadata":{}},{"cell_type":"code","source":"#Code by Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\nfig = 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":"2022-10-04T22:33:09.3021Z","iopub.execute_input":"2022-10-04T22:33:09.302555Z","iopub.status.idle":"2022-10-04T22:33:12.635252Z","shell.execute_reply.started":"2022-10-04T22:33:09.302518Z","shell.execute_reply":"2022-10-04T22:33:12.633779Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Conclusion\n\n\"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.\"\n\nBy Geir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals","metadata":{}},{"cell_type":"markdown","source":"#Acknowledgements:\n\nWouter van Zeist https://github.com/wvanzeist/riroriro\n\nGeir Drange https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals\n\n","metadata":{}}]}