{"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://gwpy.github.io/images/gwpy_1200.png)","metadata":{}},{"cell_type":"markdown","source":"# Processing of gravitational data\nIn this notebook we will apply some signal processing to the gravitational data. Preprocessing is definately required on these signals, and thankfully there is a Python package called [GWpy](https://gwpy.github.io/docs/latest/index.html) that has all functions that are needed. Theory will not be discussed here in detail, but there is plenty of info and code on the topic:   \n  * [GW tutorials](https://www.gw-openscience.org/tutorials/)\n  * [Gravitational Wave Open Science Center](https://www.gw-openscience.org/software/)\n  * [A guide to LIGO–Virgo detector noise and extraction of transient gravitational-wave signals](https://iopscience.iop.org/article/10.1088/1361-6382/ab685e)\n  \n  \nFirst, install GWpy:","metadata":{}},{"cell_type":"code","source":"!python -m pip install gwpy\n!pip install astropy==4.2.1","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-13T13:35:53.830952Z","iopub.execute_input":"2021-07-13T13:35:53.83134Z","iopub.status.idle":"2021-07-13T13:36:10.18336Z","shell.execute_reply.started":"2021-07-13T13:35:53.831262Z","shell.execute_reply":"2021-07-13T13:36:10.182105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And import a few libraries.","metadata":{}},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\nfrom gwpy.plot import Plot\nimport numpy as np\nfrom scipy import signal\nfrom sklearn.preprocessing import MinMaxScaler\nfrom PIL import Image\nfrom matplotlib import pyplot as plt","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-13T13:36:10.184926Z","iopub.execute_input":"2021-07-13T13:36:10.185212Z","iopub.status.idle":"2021-07-13T13:36:12.201471Z","shell.execute_reply.started":"2021-07-13T13:36:10.18518Z","shell.execute_reply":"2021-07-13T13:36:12.200628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Read & plot files\nLet's define helper function to read numpy data and convert into GWpy TimeSeries format and plot the data. All the files are 2s recordings at 2048Hz sample rate from the 3 detectors.","metadata":{}},{"cell_type":"code","source":"def read_file(fname):\n    data = np.load(fname)\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\ndef plot_time_data(d1, d2, d3):\n    plot = Plot(d1, d2, d3, separate=True, sharex=True, figsize=[12, 8])\n    ax = plot.gca()\n    ax.set_xlim(0,2)\n    ax.set_xlabel('Time [s]')\n    plot.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:36:12.204516Z","iopub.execute_input":"2021-07-13T13:36:12.204867Z","iopub.status.idle":"2021-07-13T13:36:12.210195Z","shell.execute_reply.started":"2021-07-13T13:36:12.204841Z","shell.execute_reply":"2021-07-13T13:36:12.209632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Visualize a file:","metadata":{}},{"cell_type":"code","source":"d1, d2, d3 = read_file('../input/g2net-gravitational-wave-detection/train/0/0/0/000a5b6e5c.npy')\nplot_time_data(d1, d2, d3)","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:36:12.211086Z","iopub.execute_input":"2021-07-13T13:36:12.211438Z","iopub.status.idle":"2021-07-13T13:36:15.632709Z","shell.execute_reply.started":"2021-07-13T13:36:12.211412Z","shell.execute_reply":"2021-07-13T13:36:15.631645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocess\nThen we will follow the general processing steps outlined in [this article](https://iopscience.iop.org/article/10.1088/1361-6382/ab685e):  \n* Apply a window function (Tukey - tapered cosine window) to suppress spectral leakage\n* Whiten the spectrum\n* Bandpass","metadata":{}},{"cell_type":"markdown","source":"## Apply window function","metadata":{}},{"cell_type":"markdown","source":"The Tukey window looks like this:","metadata":{}},{"cell_type":"code","source":"window = signal.tukey(4096)\nplt.plot(window);","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:37:06.688238Z","iopub.execute_input":"2021-07-13T13:37:06.688554Z","iopub.status.idle":"2021-07-13T13:37:06.899215Z","shell.execute_reply.started":"2021-07-13T13:37:06.688528Z","shell.execute_reply":"2021-07-13T13:37:06.898378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look at the signal after windowing:","metadata":{}},{"cell_type":"code","source":"d1, d2, d3 = d1*window, d2*window, d3*window\nplot_time_data(d1, d2, d3)","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:37:09.678506Z","iopub.execute_input":"2021-07-13T13:37:09.678882Z","iopub.status.idle":"2021-07-13T13:37:12.815142Z","shell.execute_reply.started":"2021-07-13T13:37:09.678851Z","shell.execute_reply":"2021-07-13T13:37:12.814258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Take a look at the spectrum - there is a lot of low frequency noise:","metadata":{}},{"cell_type":"code","source":"fig2 = d1.asd(fftlength=2).plot(figsize=[12, 6])\nplt.xlim(10,1024)\nplt.ylim(1e-25, 1e-20);","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:37:19.464287Z","iopub.execute_input":"2021-07-13T13:37:19.464748Z","iopub.status.idle":"2021-07-13T13:37:20.100315Z","shell.execute_reply.started":"2021-07-13T13:37:19.464699Z","shell.execute_reply":"2021-07-13T13:37:20.099554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can easily highpass the data (here with 15Hz highpass frequency):","metadata":{}},{"cell_type":"code","source":"fig2b = d1.highpass(15).asd(fftlength=2).plot(figsize=[12, 6])\nplt.xlim(10,1024)\nplt.ylim(1e-25, 1e-20);","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:37:26.327638Z","iopub.execute_input":"2021-07-13T13:37:26.327998Z","iopub.status.idle":"2021-07-13T13:37:27.027281Z","shell.execute_reply.started":"2021-07-13T13:37:26.327963Z","shell.execute_reply":"2021-07-13T13:37:27.026302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Spectral whitening and bandpass filtering\nThis is super simple with GWpy:","metadata":{}},{"cell_type":"code","source":"white_data = d1.whiten(window=(\"tukey\",0.2)) # whiten-function has a built-in window function\nbp_data = white_data.bandpass(35, 350) # frequency range 35-350Hz\nfig3 = bp_data.plot(figsize=[12, 6])\nplt.xlim(0, 2)\nax = plt.gca()\nax.set_title('Whitened and bandpassed')\nax.set_xlabel('Time [s]');","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:38:34.538274Z","iopub.execute_input":"2021-07-13T13:38:34.538593Z","iopub.status.idle":"2021-07-13T13:38:35.581127Z","shell.execute_reply.started":"2021-07-13T13:38:34.538547Z","shell.execute_reply":"2021-07-13T13:38:35.580301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we have a preprocessed data that is ready for further analysis. First, let's define a function that combines all the steps above and outputs preprocessed data:","metadata":{}},{"cell_type":"code","source":"def preprocess(d1, d2, d3, bandpass=False, 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    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","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:49:25.551469Z","iopub.execute_input":"2021-07-13T13:49:25.551808Z","iopub.status.idle":"2021-07-13T13:49:25.559007Z","shell.execute_reply.started":"2021-07-13T13:49:25.551779Z","shell.execute_reply":"2021-07-13T13:49:25.558385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Q-Transform\nThe Q-Transform is related to the Fourier transform, and very closely related to a wavelet transform. The spectrogram is a possible candidate as input for a CNN model.","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-07-13T13:51:41.129285Z","iopub.execute_input":"2021-07-13T13:51:41.129719Z","iopub.status.idle":"2021-07-13T13:51:42.939124Z","shell.execute_reply.started":"2021-07-13T13:51:41.129688Z","shell.execute_reply":"2021-07-13T13:51:42.938343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Combine three channels into one RGB image\nSince we have 3 detectors, we can combine the Q-Transforms as RGB channels into one color image. Let's make a function for that:","metadata":{}},{"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-07-13T13:52:37.701088Z","iopub.execute_input":"2021-07-13T13:52:37.701401Z","iopub.status.idle":"2021-07-13T13:52:37.708778Z","shell.execute_reply.started":"2021-07-13T13:52:37.701374Z","shell.execute_reply":"2021-07-13T13:52:37.708181Z"},"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-07-13T13:52:40.843904Z","iopub.execute_input":"2021-07-13T13:52:40.844417Z","iopub.status.idle":"2021-07-13T13:52:41.4135Z","shell.execute_reply.started":"2021-07-13T13:52:40.844388Z","shell.execute_reply":"2021-07-13T13:52:41.412724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here is a very obvious one (chirp) from the test set:","metadata":{}},{"cell_type":"code","source":"create_rgb('../input/g2net-gravitational-wave-detection/test/0/0/2/0021f9dd71.npy')","metadata":{"execution":{"iopub.status.busy":"2021-07-13T13:52:46.742045Z","iopub.execute_input":"2021-07-13T13:52:46.742381Z","iopub.status.idle":"2021-07-13T13:52:47.178584Z","shell.execute_reply.started":"2021-07-13T13:52:46.742353Z","shell.execute_reply":"2021-07-13T13:52:47.177532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Awesome! (Or crazy?) Next step is to find out if we can train an image classifier with these images...","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}