{"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)\nimport h5py\nfrom datetime import datetime\nimport matplotlib.pyplot as plt\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-05T19:11:46.011579Z","iopub.execute_input":"2022-10-05T19:11:46.012004Z","iopub.status.idle":"2022-10-05T19:11:47.225925Z","shell.execute_reply.started":"2022-10-05T19:11:46.011971Z","shell.execute_reply":"2022-10-05T19:11:47.224792Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#PyCBC open-source Gravitational- Wave Astronomy Toolkit.\n\n\"PyCBC Inference: A Python-based parameter estimation toolkit for compact binary coalescence signals\"\n\nAuthors: C. M. Biwer, Collin D. Capano, Soumi De, Miriam Cabero, Duncan A. Brown, Alexander H. Nitz, V. Raymond - https://doi.org/10.1088/1538-3873/aaef0b\n\n\"The authors introduced new modules in the open-source PyCBC gravitational- wave astronomy toolkit that implement Bayesian inference for compact-object binary mergers. We review the Bayesian inference methods implemented and describe the structure of the modules. We demonstrate that the PyCBC Inference modules produce unbiased estimates of the parameters of a simulated population of binary black hole mergers. We show that the posterior parameter distributions obtained used our new code agree well with the published estimates for binary black holes in the first LIGO-Virgo observing run.\"\n\n\nhttps://doi.org/10.1088/1538-3873/aaef0b\n\nhttps://arxiv.org/abs/1807.10312","metadata":{}},{"cell_type":"markdown","source":"![](https://igoligo.files.wordpress.com/2016/02/20151116-kagra-fig4.jpg)https://igoligo.wordpress.com/2016/02/01/getting-started-science-mode-ready/","metadata":{}},{"cell_type":"code","source":"import librosa\nimport librosa.display\n\nfrom scipy.special import logit, expit\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import Normalize\n%matplotlib inline\n\nimport cv2\nfrom scipy.interpolate import interp1d\nimport pywt\n\ndef sigmoid(var):\n    return 1/(1+np.exp(-var))","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:32:37.905097Z","iopub.execute_input":"2022-10-05T18:32:37.905563Z","iopub.status.idle":"2022-10-05T18:32:39.585707Z","shell.execute_reply.started":"2022-10-05T18:32:37.905527Z","shell.execute_reply":"2022-10-05T18:32:39.58437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Code by Quentin R https://www.kaggle.com/code/qrendu/plot-spectrograms\n\nf = h5py.File('/kaggle/input/g2net-detecting-continuous-gravitational-waves/train/'+filename, 'r')","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:30:05.73489Z","iopub.execute_input":"2022-10-05T18:30:05.735357Z","iopub.status.idle":"2022-10-05T18:30:05.751145Z","shell.execute_reply.started":"2022-10-05T18:30:05.73532Z","shell.execute_reply":"2022-10-05T18:30:05.749517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install pycbc\nimport pycbc\nfrom pycbc.waveform import td_approximants, fd_approximants, get_td_waveform\nfrom pycbc.detector import Detector","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:30:44.658769Z","iopub.execute_input":"2022-10-05T18:30:44.659209Z","iopub.status.idle":"2022-10-05T18:32:08.786302Z","shell.execute_reply.started":"2022-10-05T18:30:44.659147Z","shell.execute_reply":"2022-10-05T18:32:08.785046Z"},"_kg_hide-output":true,"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Waveforms\n\nWhat waveforms can I generate?\n\nhttps://pycbc.org/pycbc/latest/html/waveform.html","metadata":{}},{"cell_type":"code","source":"# List of td approximants that are available\nprint(td_approximants())","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:33:06.356295Z","iopub.execute_input":"2022-10-05T18:33:06.356768Z","iopub.status.idle":"2022-10-05T18:33:06.363718Z","shell.execute_reply.started":"2022-10-05T18:33:06.356727Z","shell.execute_reply":"2022-10-05T18:33:06.362357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Method to generate FD waveform using TD approximant\n\nFrequency-domain (FD) Waveform and Time-domain (TD).","metadata":{}},{"cell_type":"code","source":"# List of fd approximants that are currently available\nprint(fd_approximants())","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:33:24.868211Z","iopub.execute_input":"2022-10-05T18:33:24.868768Z","iopub.status.idle":"2022-10-05T18:33:24.877229Z","shell.execute_reply.started":"2022-10-05T18:33:24.868724Z","shell.execute_reply":"2022-10-05T18:33:24.875751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Plotting Time Domain Waveforms","metadata":{}},{"cell_type":"code","source":"#Code by PyCBC https://pycbc.org/pycbc/latest/html/waveform.html\n\nimport matplotlib.pyplot as pp\nfrom pycbc.waveform import get_td_waveform\n\n\nfor apx in ['SEOBNRv4', 'IMRPhenomD']:\n    hp, hc = get_td_waveform(approximant=apx,\n                                 mass1=10,\n                                 mass2=10,\n                                 spin1z=0.9,\n                                 delta_t=1.0/4096,\n                                 f_lower=40)\n\n    pp.plot(hp.sample_times, hp, label=apx)\n\npp.ylabel('Strain')\npp.xlabel('Time (s)')\npp.legend()\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:33:29.240759Z","iopub.execute_input":"2022-10-05T19:33:29.241177Z","iopub.status.idle":"2022-10-05T19:33:30.077757Z","shell.execute_reply.started":"2022-10-05T19:33:29.24113Z","shell.execute_reply":"2022-10-05T19:33:30.076597Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Generating one waveform in multiple detectors","metadata":{}},{"cell_type":"code","source":"#Code by PyCBC https://pycbc.org/pycbc/latest/html/waveform.html\n\nimport matplotlib.pyplot as pp\nfrom pycbc.waveform import get_td_waveform\nfrom pycbc.detector import Detector\n\n\napx = 'SEOBNRv4'\n# NOTE: Inclination runs from 0 to pi, with poles at 0 and pi\n#       coa_phase runs from 0 to 2 pi.\nhp, hc = get_td_waveform(approximant=apx,\n                         mass1=10,\n                         mass2=10,\n                         spin1z=0.9,\n                         spin2z=0.4,\n                         inclination=1.23,\n                         coa_phase=2.45,\n                         delta_t=1.0/4096,\n                         f_lower=40)\n\ndet_h1 = Detector('H1')\ndet_l1 = Detector('L1')\ndet_v1 = Detector('V1')\n\n# Choose a GPS end time, sky location, and polarization phase for the merger\n# NOTE: Right ascension and polarization phase runs from 0 to 2pi\n#       Declination runs from pi/2. to -pi/2 with the poles at pi/2. and -pi/2.\nend_time = 1192529720\ndeclination = 0.65\nright_ascension = 4.67\npolarization = 2.34\nhp.start_time += end_time\nhc.start_time += end_time\n\nsignal_h1 = det_h1.project_wave(hp, hc,  right_ascension, declination, polarization)\nsignal_l1 = det_l1.project_wave(hp, hc,  right_ascension, declination, polarization)\nsignal_v1 = det_v1.project_wave(hp, hc,  right_ascension, declination, polarization)\n\npp.plot(signal_h1.sample_times, signal_h1, label='H1')\npp.plot(signal_l1.sample_times, signal_l1, label='L1')\npp.plot(signal_v1.sample_times, signal_v1, label='V1')\n\npp.ylabel('Strain')\npp.xlabel('Time (s)')\npp.legend()\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:35:08.15729Z","iopub.execute_input":"2022-10-05T19:35:08.157765Z","iopub.status.idle":"2022-10-05T19:35:09.045312Z","shell.execute_reply.started":"2022-10-05T19:35:08.157728Z","shell.execute_reply":"2022-10-05T19:35:09.044123Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Plotting GW phase and amplitude of TD waveform","metadata":{}},{"cell_type":"code","source":"#Code by PyCBC https://pycbc.org/pycbc/latest/html/waveform.html\n\nimport matplotlib.pyplot as pp\nfrom pycbc import waveform\n\n\nfor apx in ['SEOBNRv4', 'TaylorT4', 'IMRPhenomB']:\n    hp, hc = waveform.get_td_waveform(approximant=apx,\n                                 mass1=10,\n                                 mass2=10,\n                                 delta_t=1.0/4096,\n                                 f_lower=40)\n\n    hp, hc = hp.trim_zeros(), hc.trim_zeros()\n    amp = waveform.utils.amplitude_from_polarizations(hp, hc)\n    phase = waveform.utils.phase_from_polarizations(hp, hc)\n\n    pp.plot(phase, amp, label=apx)\n\npp.ylabel('GW Strain Amplitude')\npp.xlabel('GW Phase (radians)')\npp.legend(loc='upper left')\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:37:29.15874Z","iopub.execute_input":"2022-10-05T19:37:29.159255Z","iopub.status.idle":"2022-10-05T19:37:30.32876Z","shell.execute_reply.started":"2022-10-05T19:37:29.159214Z","shell.execute_reply":"2022-10-05T19:37:30.327473Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Plotting frequency evolution of TD waveform","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as pp\nfrom pycbc import waveform\n\n\nfor phase_order in [2, 3, 4, 5, 6, 7]:\n    hp, hc = waveform.get_td_waveform(approximant='SpinTaylorT4',\n                                 mass1=10, mass2=10,\n                                 phase_order=phase_order,\n                                 delta_t=1.0/4096,\n                                 f_lower=100)\n\n    hp, hc = hp.trim_zeros(), hc.trim_zeros()\n    amp = waveform.utils.amplitude_from_polarizations(hp, hc)\n    f = waveform.utils.frequency_from_polarizations(hp, hc)\n\n    pp.plot(f.sample_times, f, label=\"PN Order = %s\" % phase_order)\n\npp.ylabel('Frequency (Hz)')\npp.xlabel('Time (s)')\npp.legend(loc='upper left')\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:38:51.928868Z","iopub.execute_input":"2022-10-05T19:38:51.929312Z","iopub.status.idle":"2022-10-05T19:38:52.13117Z","shell.execute_reply.started":"2022-10-05T19:38:51.929276Z","shell.execute_reply":"2022-10-05T19:38:52.129328Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#By Giba's\n\ngwlist = ['SEOBNRv2', 'SEOBNRv2_opt', 'SEOBNRv4', 'SEOBNRv4_opt', 'SEOBNRv2T', 'SEOBNRv4T', ]\ngwlist = ['SEOBNRv2', 'SEOBNRv4' ]","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:34:03.196444Z","iopub.execute_input":"2022-10-05T18:34:03.196927Z","iopub.status.idle":"2022-10-05T18:34:03.203898Z","shell.execute_reply.started":"2022-10-05T18:34:03.196888Z","shell.execute_reply":"2022-10-05T18:34:03.202064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#PyCBC and Giba's code\n\nPyCBC https://pycbc.org/pycbc/latest/html/credit.html\n\nGiba https://www.kaggle.com/code/titericz/simulated-gw ","metadata":{}},{"cell_type":"code","source":"#Code by Giba https://www.kaggle.com/code/titericz/simulated-gw \n\nnp.random.seed(1)\n\n#Define the Detectors\ndet_h1 = Detector('H1')\ndet_l1 = Detector('L1')\ndet_v1 = Detector('V1')\n\nfor n in range(100):\n    #Define the GW params\n    gwapprox = np.random.choice( gwlist )\n    print(gwapprox)\n    hp, hc = get_td_waveform(approximant=gwapprox,\n                             mass1=16 + np.random.randint(0,10),\n                             mass2=16 + np.random.randint(0,10),\n                             delta_t=1.0/4096,\n                             spin1z=0.5 + np.random.rand()*0.5,\n                             spin2z=0.25 + np.random.rand()*0.5,\n                             inclination= 2 * np.pi * np.random.rand(),\n                             coa_phase= 2 * np.pi * np.random.rand(),\n                             phase_order = np.random.randint(2,8),\n                             f_lower=np.random.randint(24,64),\n                             distance=int(np.random.randint(1,1000)),\n                            )\n\n    \n    # Choose a GPS end time, sky location, and polarization phase for the merger\n    # NOTE: Right ascension and polarization phase runs from 0 to 2pi\n    #       Declination runs from pi/2. to -pi/2 with the poles at pi/2. and -pi/2.\n    end_time = 1192529720 + np.random.randint(1192529720//100000)\n    declination = np.pi * np.random.rand() - np.pi/2\n    right_ascension = 2 * np.pi * np.random.rand()\n    polarization = 2 * np.pi * np.random.rand()\n    hp.start_time += end_time\n    hc.start_time += end_time\n\n    signal_h1 = det_h1.project_wave(hp, hc,  right_ascension, declination, polarization)\n    signal_l1 = det_l1.project_wave(hp, hc,  right_ascension, declination, polarization)\n    signal_v1 = det_v1.project_wave(hp, hc,  right_ascension, declination, polarization)    \n    minlen = np.min( [len(signal_h1), len(signal_l1), len(signal_v1)] )\n    data = np.stack( (signal_h1[:minlen], signal_l1[:minlen], signal_v1[:minlen]),  ) * 1e19\n    print(data.shape)\n    \n    if data.shape[1]>4096:\n        data = data[:,data.shape[1]-4096:]\n        for N in range(80):\n            data[:,N] *= 1./(N+1)\n    \n    if len(hp)<4096:\n        for N in range(80):\n            data[:,N] *= 1./(N+1)\n        data = np.pad(data, ((0,0),(4096-data.shape[1],0)) )\n    \n    \n    plt.plot(data[0])\n    plt.plot(data[1])\n    plt.plot(data[2])\n    plt.show()\n    \n    cwt, freqs = pywt.cwt(data, scales=np.arange(1, 95, 0.62), wavelet='cmor1.5-0.95', sampling_period=1/2048, method='fft')\n    cwt = cwt.transpose(0,2,1)\n    print(cwt.shape)\n    cwt = np.log1p( np.abs(cwt) )\n    print( cwt.min(), cwt.max())\n    cwt -= cwt.min()\n    cwt /= cwt.max()\n    plt.imshow(cv2.resize(cwt,(256, 256)) )\n    plt.title(gwapprox)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:34:44.872855Z","iopub.execute_input":"2022-10-05T18:34:44.873966Z","iopub.status.idle":"2022-10-05T18:36:16.467304Z","shell.execute_reply.started":"2022-10-05T18:34:44.87391Z","shell.execute_reply":"2022-10-05T18:36:16.466016Z"},"_kg_hide-input":true,"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Signal Processing with GW150914\n\nHere are some interesting examples of how to process LIGO data using GW150914 as an example.\n\n#Plotting the whitened strain\n\nhttps://pycbc.org/pycbc/latest/html/gw150914.html","metadata":{}},{"cell_type":"code","source":"#Code by PyCBC https://pycbc.org/pycbc/latest/html/gw150914.html\n\nimport matplotlib.pyplot as pp\nfrom pycbc.filter import highpass_fir, lowpass_fir\nfrom pycbc.psd import welch, interpolate\nfrom pycbc.catalog import Merger\n\n\nfor ifo in ['H1', 'L1']:\n    # Read data and remove low frequency content\n    h1 = Merger(\"GW150914\").strain(ifo)\n    h1 = highpass_fir(h1, 15, 8)\n\n    # Calculate the noise spectrum\n    psd = interpolate(welch(h1), 1.0 / h1.duration)\n\n    # whiten\n    white_strain = (h1.to_frequencyseries() / psd ** 0.5).to_timeseries()\n\n    # remove some of the high and low\n    smooth = highpass_fir(white_strain, 35, 8)\n    smooth = lowpass_fir(white_strain, 300, 8)\n\n    # time shift and flip L1\n    if ifo == 'L1':\n        smooth *= -1\n        smooth.roll(int(.007 / smooth.delta_t))\n\n    pp.plot(smooth.sample_times, smooth, label=ifo)\n\npp.legend()\npp.xlim(1126259462.21, 1126259462.45)\npp.ylim(-150, 150)\npp.ylabel('Smoothed-Whitened Strain')\npp.grid()\npp.xlabel('GPS Time (s)')\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T18:56:27.40244Z","iopub.execute_input":"2022-10-05T18:56:27.402874Z","iopub.status.idle":"2022-10-05T18:56:33.263128Z","shell.execute_reply.started":"2022-10-05T18:56:27.402837Z","shell.execute_reply":"2022-10-05T18:56:33.261921Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Calculate the signal-to-noise\n\nI have no clue how to change the hdf5 files to the gwf of this snippet below.\n\nTherefore, I didn't calculated the noise in this Competition Dataset.","metadata":{}},{"cell_type":"code","source":"#Code by PyCBC https://pycbc.org/pycbc/latest/html/gw150914.html\n\nimport matplotlib.pyplot as pp\nfrom urllib.request import urlretrieve\nfrom pycbc.frame import read_frame\nfrom pycbc.filter import highpass_fir, matched_filter\nfrom pycbc.waveform import get_fd_waveform\nfrom pycbc.psd import welch, interpolate\n\n\n# Read data and remove low frequency content\nfname = 'H-H1_LOSC_4_V2-1126259446-32.gwf'\nurl = \"https://www.gw-openscience.org/GW150914data/\" + fname\nurlretrieve(url, filename=fname)\nh1 = read_frame('H-H1_LOSC_4_V2-1126259446-32.gwf', 'H1:LOSC-STRAIN')\nh1 = highpass_fir(h1, 15, 8)\n\n# Calculate the noise spectrum\npsd = interpolate(welch(h1), 1.0 / h1.duration)\n\n# Generate a template to filter with\nhp, hc = get_fd_waveform(approximant=\"IMRPhenomD\", mass1=40, mass2=32,\n                         f_lower=20, delta_f=1.0/h1.duration)\nhp.resize(len(h1) // 2 + 1)\n\n# Calculate the complex (two-phase SNR)\nsnr = matched_filter(hp, h1, psd=psd, low_frequency_cutoff=20.0)\n\n# Remove regions corrupted by filter wraparound\nsnr = snr[len(snr) // 4: len(snr) * 3 // 4]\n\npp.plot(snr.sample_times, abs(snr))\npp.ylabel('signal-to-noise')\npp.xlabel('GPS Time (s)')\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:02:09.467555Z","iopub.execute_input":"2022-10-05T19:02:09.468041Z","iopub.status.idle":"2022-10-05T19:02:11.255443Z","shell.execute_reply.started":"2022-10-05T19:02:09.467998Z","shell.execute_reply":"2022-10-05T19:02:11.254242Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Generating an Analytic PSD from lalsimulation\n\n\"A certain number of PSDs are built into lalsimulation, which you’ll be able to access through PyCBC. Below we show how to see which ones are available, and demonstrate how to generate one.\"\n\nhttps://pycbc.org/pycbc/latest/html/psd.html\n\n#Power Spectral Density (PSD)\n\n\"The power spectral density (PSD) then refers to the spectral energy distribution that would be found per unit time, since the total energy of such a signal over all time would generally be infinite.\"\n\nhttps://en.wikipedia.org/wiki/Spectral_density","metadata":{}},{"cell_type":"markdown","source":"![](https://slideplayer.com/slide/4757604/15/images/28/Power+spectrum+density.jpg)https://slideplayer.com/slide/4757604/","metadata":{}},{"cell_type":"code","source":"#Code by https://pycbc.org/pycbc/latest/html/psd.html\n\nimport matplotlib.pyplot as pp\nimport pycbc.psd\n\n\n# List the available analytic psds\n#print(pycbc.psd.get_lalsim_psd_list())\n\ndelta_f = 1.0 / 4\nflen = int(1024 / delta_f)\nlow_frequency_cutoff = 30.0\n\n# One can either call the psd generator by name\np1 = pycbc.psd.aLIGOZeroDetHighPower(flen, delta_f, low_frequency_cutoff)\n\n# or by using the name as a string.\np2 = pycbc.psd.from_string('aLIGOZeroDetLowPower', flen, delta_f, low_frequency_cutoff)\n\npp.plot(p1.sample_frequencies, p1, label='HighPower')\npp.plot(p2.sample_frequencies, p2, label='LowPower')\npp.legend()\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:23:02.968388Z","iopub.execute_input":"2022-10-05T19:23:02.968884Z","iopub.status.idle":"2022-10-05T19:23:03.192043Z","shell.execute_reply.started":"2022-10-05T19:23:02.968844Z","shell.execute_reply":"2022-10-05T19:23:03.190872Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Estimating the PSD of a time series","metadata":{}},{"cell_type":"code","source":"#Code by PyCBC https://pycbc.org/pycbc/latest/html/waveform.html\n\nimport matplotlib.pyplot as pp\nimport pycbc.noise\nimport pycbc.psd\n\n\n# generate some colored gaussian noise\nflow = 30.0\ndelta_f = 1.0 / 16\nflen = int(2048 / delta_f) + 1\npsd = pycbc.psd.aLIGOZeroDetHighPower(flen, delta_f, flow)\n\n### Generate 128 seconds of noise at 4096 Hz\ndelta_t = 1.0 / 4096\ntsamples = int(128 / delta_t)\nts = pycbc.noise.noise_from_psd(tsamples, delta_t, psd, seed=127)\n\n# Estimate the PSD\n# We'll choose 4 seconds PSD samples that are overlapped 50 %\nseg_len = int(4 / delta_t)\nseg_stride = int(seg_len / 2)\nestimated_psd = pycbc.psd.welch(ts,\n                      seg_len=seg_len,\n                      seg_stride=seg_stride)\n\npp.loglog(estimated_psd.sample_frequencies, estimated_psd, label='estimate')\npp.loglog(psd.sample_frequencies, psd, linewidth=3, label='known psd')\npp.xlim(xmin=flow, xmax=2000)\npp.ylim(1e-48, 1e-45)\npp.legend()\npp.grid()\npp.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:23:38.34579Z","iopub.execute_input":"2022-10-05T19:23:38.346239Z","iopub.status.idle":"2022-10-05T19:23:39.483343Z","shell.execute_reply.started":"2022-10-05T19:23:38.3462Z","shell.execute_reply":"2022-10-05T19:23:39.481746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Acknowledgements:\n\nPyCBC https://pycbc.org/pycbc/latest/html/waveform.html\n\nPyCBC Live -  Rapid Detection of Gravitational Waves from Compact Binary Mergers https://granite.phys.s.u-tokyo.ac.jp/kita/Seminar/presentation_2018_0608.pdf\n\nGiba https://www.kaggle.com/code/titericz/simulated-gw\n\nQuentin R https://www.kaggle.com/code/qrendu/plot-spectrograms\n\n","metadata":{}}]}