{"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":"# GWpy Tutorial\n\nGWpy is a collaboration-driven Python package providing tools for studying data from ground-based gravitational-wave detectors.\n\nGWpy provides a user-friendly, intuitive interface to the common time-domain and frequency-domain data produced by the LIGO and Virgo instruments and their analysis, with easy-to-follow tutorials at each step.\n\n- https://gwpy.github.io/docs/v0.1/timeseries/index.html\n- https://www.kaggle.com/mistag/data-preprocessing-with-gwpy","metadata":{"execution":{"iopub.status.busy":"2021-08-10T14:44:20.645292Z","iopub.execute_input":"2021-08-10T14:44:20.645816Z","iopub.status.idle":"2021-08-10T14:44:20.656387Z","shell.execute_reply.started":"2021-08-10T14:44:20.645735Z","shell.execute_reply":"2021-08-10T14:44:20.655028Z"}}},{"cell_type":"code","source":"!python -m pip install gwpy\n!pip install astropy==4.2.1","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:37:21.542261Z","iopub.execute_input":"2021-08-10T17:37:21.542570Z","iopub.status.idle":"2021-08-10T17:37:52.875130Z","shell.execute_reply.started":"2021-08-10T17:37:21.542542Z","shell.execute_reply":"2021-08-10T17:37:52.874219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\n\nfrom gwpy.timeseries import TimeSeries\nfrom gwpy.plot import Plot\nfrom scipy import signal\nfrom sklearn.preprocessing import MinMaxScaler\nfrom PIL import Image\nfrom matplotlib import pyplot as plt\nfrom sklearn.preprocessing import normalize\nfrom sklearn.preprocessing import StandardScaler","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:37:52.877798Z","iopub.execute_input":"2021-08-10T17:37:52.878086Z","iopub.status.idle":"2021-08-10T17:37:54.976740Z","shell.execute_reply.started":"2021-08-10T17:37:52.878054Z","shell.execute_reply":"2021-08-10T17:37:54.975832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\n\ndata = TimeSeries.fetch_open_data('L1', 1126259446, 1126259478)\ndata","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:37:54.978774Z","iopub.execute_input":"2021-08-10T17:37:54.979171Z","iopub.status.idle":"2021-08-10T17:38:18.004383Z","shell.execute_reply.started":"2021-08-10T17:37:54.979132Z","shell.execute_reply":"2021-08-10T17:38:18.003440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_file(fname):\n    data = np.load(fname)\n    data = rescale(data)\n    d1 = TimeSeries(data[0,:], sample_rate=2048)\n    d2 = TimeSeries(data[1,:], sample_rate=2048)\n    d3 = TimeSeries(data[2,:], sample_rate=2048)  \n    return d1, d2, d3\n\n\ndef rescale(data):\n    # Normalize to [-1,1]\n    data = (data - np.mean(data, axis=1).reshape(3,-1)) / (np.max(data, axis=1) - np.min(data, axis=1)).reshape(3,-1) * 2\n    return data\n\n    \ndef plot_time_data(d1, d2, d3, title=''):\n    plot = Plot(d1, d2, d3, separate=True, sharex=True, figsize=[12, 8])\n    ax = plot.gca()   \n    ax.set_title(title)\n    ax.set_xlim(0,2)\n    ax.set_xlabel('Time [s]')\n    plot.show()\n        \n        \ndef preprocess(d1, d2, d3, bandpass=False, normalize=True, lf=35, hf=350):\n    white_d1 = d1.whiten(window=(\"tukey\",0.2))\n    white_d2 = d2.whiten(window=(\"tukey\",0.2))\n    white_d3 = d3.whiten(window=(\"tukey\",0.2))\n        \n    if bandpass: # bandpass filter\n        bp_d1 = white_d1.bandpass(lf, hf) \n        bp_d2 = white_d2.bandpass(lf, hf)\n        bp_d3 = white_d3.bandpass(lf, hf)\n        return bp_d1, bp_d2, bp_d3\n    else: # only whiten\n        return white_d1, white_d2, white_d3\n    \nsample_signals    = read_file('/kaggle/input/g2net-gravitational-wave-detection/train/0/0/0/00000e74ad.npy')\nsample_no_signals = read_file('/kaggle/input/g2net-gravitational-wave-detection/train/0/0/0/00001f4945.npy')","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:18.005862Z","iopub.execute_input":"2021-08-10T17:38:18.006109Z","iopub.status.idle":"2021-08-10T17:38:18.044116Z","shell.execute_reply.started":"2021-08-10T17:38:18.006084Z","shell.execute_reply":"2021-08-10T17:38:18.043068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_time_data(*sample_signals, title='raw')","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:18.045263Z","iopub.execute_input":"2021-08-10T17:38:18.045660Z","iopub.status.idle":"2021-08-10T17:38:21.595362Z","shell.execute_reply.started":"2021-08-10T17:38:18.045619Z","shell.execute_reply":"2021-08-10T17:38:21.594710Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_time_data(*preprocess(*sample_signals, bandpass=False), title='whiten')","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:21.596260Z","iopub.execute_input":"2021-08-10T17:38:21.596627Z","iopub.status.idle":"2021-08-10T17:38:25.335298Z","shell.execute_reply.started":"2021-08-10T17:38:21.596580Z","shell.execute_reply":"2021-08-10T17:38:25.334643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_time_data(*preprocess(*sample_signals, bandpass=True),  title='bandpass')","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:25.336249Z","iopub.execute_input":"2021-08-10T17:38:25.336598Z","iopub.status.idle":"2021-08-10T17:38:28.832596Z","shell.execute_reply.started":"2021-08-10T17:38:25.336572Z","shell.execute_reply":"2021-08-10T17:38:28.831425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.array( sample_signals[0].to_value().shape ) ","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:28.834498Z","iopub.execute_input":"2021-08-10T17:38:28.834838Z","iopub.status.idle":"2021-08-10T17:38:28.841192Z","shell.execute_reply.started":"2021-08-10T17:38:28.834807Z","shell.execute_reply":"2021-08-10T17:38:28.840209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Q-Transform","metadata":{}},{"cell_type":"code","source":"r1, r2, r3 = read_file('../input/g2net-gravitational-wave-detection/train/0/0/0/000a5b6e5c.npy') # this signal has target=1\np1, p2, p3 = preprocess(r1, r2, r3)\nhq = p2.q_transform(qrange=(16,32), frange=(30,400), logf=True, whiten=False)\nfig4 = hq.plot(figsize=[12, 10])\nax = fig4.gca()\nfig4.colorbar(label=\"Normalised energy\")\nax.grid(False)\nax.set_yscale('log')\nax.set_xlabel('Time [s]');","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:28.842963Z","iopub.execute_input":"2021-08-10T17:38:28.843334Z","iopub.status.idle":"2021-08-10T17:38:30.896482Z","shell.execute_reply.started":"2021-08-10T17:38:28.843294Z","shell.execute_reply":"2021-08-10T17:38:30.895596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Q_RANGE = (16,32)\nF_RANGE = (30,400)\n\ndef create_rgb(fname):\n    r1, r2, r3 = read_file(fname)\n    p1, p2, p3 = preprocess(r1, r2, r3)\n    hq1 = p1.q_transform(qrange=Q_RANGE, frange=F_RANGE, logf=True, whiten=False)\n    hq2 = p2.q_transform(qrange=Q_RANGE, frange=F_RANGE, logf=True, whiten=False)\n    hq3 = p3.q_transform(qrange=Q_RANGE, frange=F_RANGE, logf=True, whiten=False)\n    img = np.zeros([hq1.shape[0], hq1.shape[1], 3], dtype=np.uint8)\n    scaler = MinMaxScaler()\n    img[:,:,0] = 255*scaler.fit_transform(hq1)\n    img[:,:,1] = 255*scaler.fit_transform(hq2)\n    img[:,:,2] = 255*scaler.fit_transform(hq3)\n    return Image.fromarray(img).rotate(90, expand=1).resize((760,760))","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:30.897824Z","iopub.execute_input":"2021-08-10T17:38:30.898094Z","iopub.status.idle":"2021-08-10T17:38:30.905384Z","shell.execute_reply.started":"2021-08-10T17:38:30.898067Z","shell.execute_reply":"2021-08-10T17:38:30.904481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"create_rgb('../input/g2net-gravitational-wave-detection/train/0/0/0/000a5b6e5c.npy')","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:38:30.906500Z","iopub.execute_input":"2021-08-10T17:38:30.906816Z","iopub.status.idle":"2021-08-10T17:38:31.471580Z","shell.execute_reply.started":"2021-08-10T17:38:30.906787Z","shell.execute_reply":"2021-08-10T17:38:31.470503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Signal Processing\n\nhttps://gwpy.github.io/docs/v0.1/signal/index.html\n\n## rayleigh_spectrum([fftlength, overlap])\t\n\nhttps://gwpy.github.io/docs/stable/examples/spectrogram/rayleigh.html","metadata":{}},{"cell_type":"code","source":"r1, r2, r3 = read_file('../input/g2net-gravitational-wave-detection/train/0/0/0/000a5b6e5c.npy') # this signal has target=1\np1, p2, p3 = preprocess(r1, r2, r3)\nhq = p2.rayleigh_spectrogram(2)\nhq\nfig4 = hq.plot(figsize=[12, 10])\nax = fig4.gca()\nfig4.colorbar(label=\"Rayleigh Spectrogram\")\nax.grid(False)\nax.set_yscale('log')\nax.set_xlabel('Time [s]');","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:39:52.662352Z","iopub.execute_input":"2021-08-10T17:39:52.662953Z","iopub.status.idle":"2021-08-10T17:39:54.277696Z","shell.execute_reply.started":"2021-08-10T17:39:52.662913Z","shell.execute_reply":"2021-08-10T17:39:54.277017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Source: https://gwpy.github.io/docs/stable/examples/spectrogram/rayleigh.html\nrayleigh = p2.rayleigh_spectrogram(stride=.5, fftlength=.2, overlap=.1)\nprint(rayleigh)\n\nplot = rayleigh.plot(norm='log', vmin=0.25, vmax=4)\nax = plot.gca()\nax.set_yscale('log')\nax.set_ylim(30, 1500)\nax.set_title('Sensitivity of LIGO-Livingston around GW151226')\nax.colorbar(cmap='coolwarm', label='Rayleigh statistic')\nplot.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:47:32.948216Z","iopub.execute_input":"2021-08-10T17:47:32.948556Z","iopub.status.idle":"2021-08-10T17:47:34.864798Z","shell.execute_reply.started":"2021-08-10T17:47:32.948527Z","shell.execute_reply":"2021-08-10T17:47:34.863768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Coherence","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:43:44.410714Z","iopub.execute_input":"2021-08-10T17:43:44.411029Z","iopub.status.idle":"2021-08-10T17:43:44.414320Z","shell.execute_reply.started":"2021-08-10T17:43:44.411002Z","shell.execute_reply":"2021-08-10T17:43:44.413687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"r1, r2, r3 = read_file('../input/g2net-gravitational-wave-detection/train/0/0/0/000a5b6e5c.npy') # this signal has target=1\np1, p2, p3 = preprocess(r1, r2, r3)\nhq = p1.coherence(p2) # + p2.coherence(p3, 2) + p3.coherence(p1, 2)\nhq\n# fig4 = hq.plot(figsize=[12, 100])\n# ax = fig4.gca()\n# fig4.colorbar(label=\"Rayleigh Spectrogram\")\n# ax.grid(False)\n# ax.set_yscale('log')\n# ax.set_xlabel('Time [s]');","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:44:08.134318Z","iopub.execute_input":"2021-08-10T17:44:08.134669Z","iopub.status.idle":"2021-08-10T17:44:08.173342Z","shell.execute_reply.started":"2021-08-10T17:44:08.134638Z","shell.execute_reply":"2021-08-10T17:44:08.172456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## FrequencySeries Spectral Density\n\nA measurement of the acceleration spectral density (ASD) is the usual way to specify random vibration. The root mean square acceleration (Grms) is the square root of the area under the ASD curve in the frequency domain. The Grms value is typically used to express the overall energy of a particular random vibration event and is a statistical value used in mechanical engineering for structural design and analysis purposes.\n\nWhile the term power spectral density (PSD) is commonly used to specify a random vibration event, ASD is more appropriate when acceleration is being measured and used in structural analysis and testing.\n\n- https://gwpy.github.io/docs/v0.1/timeseries/index.html#gwpy.timeseries.TimeSeries.psd\n- https://en.wikipedia.org/wiki/Random_vibration","metadata":{}},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\ngwdata = p2\n# gwdata = TimeSeries.fetch('H1:LDAS-STRAIN', 'September 16 2010 06:40', 'September 16 2010 06:50')\nspectrum = gwdata.asd(fftlength=.2, overlap=.1)\nplot = spectrum.plot()\nax = plot.gca()\n# ax.set_xlim(40, 4000)\nax.set_ylabel(r'Gravitational-wave strain ASD [strain$/\\sqrt{\\mathrm{Hz}}$]')\n# ax.set_ylim(1e-23, 1e-19)\nplot.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:54:20.085166Z","iopub.execute_input":"2021-08-10T17:54:20.085480Z","iopub.status.idle":"2021-08-10T17:54:20.624176Z","shell.execute_reply.started":"2021-08-10T17:54:20.085453Z","shell.execute_reply":"2021-08-10T17:54:20.623304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\ngwdata = p2\n# gwdata = TimeSeries.fetch('H1:LDAS-STRAIN', 'September 16 2010 06:40', 'September 16 2010 06:50')\nspectrum = gwdata.psd(fftlength=.2, overlap=.1)\nplot = spectrum.plot()\nax = plot.gca()\n# ax.set_xlim(40, 4000)\nax.set_ylabel(r'Gravitational-wave strain PSD [strain$/\\sqrt{\\mathrm{Hz}}$]')\n# ax.set_ylim(1e-23, 1e-19)\nplot.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-10T17:54:33.773375Z","iopub.execute_input":"2021-08-10T17:54:33.773736Z","iopub.status.idle":"2021-08-10T17:54:34.421306Z","shell.execute_reply.started":"2021-08-10T17:54:33.773696Z","shell.execute_reply":"2021-08-10T17:54:34.420389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}