{"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":"This code is directly copied from https://www.kaggle.com/mistag/data-preprocessing-with-gwpy by Geir Drange (please upvote his notebook), I'm publishing it only to show how I used this code to create a preprocessed dataset.","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 are 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","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-02T09:22:58.065402Z","iopub.execute_input":"2021-07-02T09:22:58.065979Z","iopub.status.idle":"2021-07-02T09:23:04.717325Z","shell.execute_reply.started":"2021-07-02T09:22:58.065847Z","shell.execute_reply":"2021-07-02T09:23:04.71628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2021-07-02T09:23:04.719022Z","iopub.execute_input":"2021-07-02T09:23:04.719305Z","iopub.status.idle":"2021-07-02T09:23:06.103143Z","shell.execute_reply.started":"2021-07-02T09:23:04.719273Z","shell.execute_reply":"2021-07-02T09:23:06.102266Z"},"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.","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-02T09:23:06.104706Z","iopub.execute_input":"2021-07-02T09:23:06.104997Z","iopub.status.idle":"2021-07-02T09:23:06.111743Z","shell.execute_reply.started":"2021-07-02T09:23:06.10497Z","shell.execute_reply":"2021-07-02T09:23:06.111075Z"},"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/0002b64784.npy')\nplot_time_data(d1, d2, d3)","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:06.112935Z","iopub.execute_input":"2021-07-02T09:23:06.113293Z","iopub.status.idle":"2021-07-02T09:23:09.719434Z","shell.execute_reply.started":"2021-07-02T09:23:06.113266Z","shell.execute_reply":"2021-07-02T09:23:09.718764Z"},"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-02T09:23:09.720447Z","iopub.execute_input":"2021-07-02T09:23:09.720857Z","iopub.status.idle":"2021-07-02T09:23:10.094992Z","shell.execute_reply.started":"2021-07-02T09:23:09.720799Z","shell.execute_reply":"2021-07-02T09:23:10.094055Z"},"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-02T09:23:10.09645Z","iopub.execute_input":"2021-07-02T09:23:10.09682Z","iopub.status.idle":"2021-07-02T09:23:13.465452Z","shell.execute_reply.started":"2021-07-02T09:23:10.096788Z","shell.execute_reply":"2021-07-02T09:23:13.464509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Take a look at the spectrum:","metadata":{}},{"cell_type":"code","source":"fig2 = d1.asd(fftlength=2).plot(figsize=[12, 6])\nplt.xlim(10,1024)\nplt.ylim(1e-24, 1e-19);","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:13.466712Z","iopub.execute_input":"2021-07-02T09:23:13.467003Z","iopub.status.idle":"2021-07-02T09:23:14.21251Z","shell.execute_reply.started":"2021-07-02T09:23:13.466975Z","shell.execute_reply":"2021-07-02T09:23:14.211391Z"},"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()\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-02T09:23:14.214093Z","iopub.execute_input":"2021-07-02T09:23:14.214496Z","iopub.status.idle":"2021-07-02T09:23:15.475891Z","shell.execute_reply.started":"2021-07-02T09:23:14.214455Z","shell.execute_reply":"2021-07-02T09:23:15.474727Z"},"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, lf=35, hf=350):\n    window = signal.tukey(4096)\n    d1, d2, d3 = d1*window, d2*window, d3*window\n    white_d1 = d1.whiten()\n    white_d2 = d2.whiten()\n    white_d3 = d3.whiten()\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","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:15.478734Z","iopub.execute_input":"2021-07-02T09:23:15.479078Z","iopub.status.idle":"2021-07-02T09:23:15.485589Z","shell.execute_reply.started":"2021-07-02T09:23:15.479047Z","shell.execute_reply":"2021-07-02T09:23:15.484417Z"},"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/0002b64784.npy')\np1, p2, p3 = preprocess(r1, r2, r3)\nhq = p1.q_transform(outseg=(0, 2))\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-02T09:23:15.487791Z","iopub.execute_input":"2021-07-02T09:23:15.488231Z","iopub.status.idle":"2021-07-02T09:23:17.708517Z","shell.execute_reply.started":"2021-07-02T09:23:15.488186Z","shell.execute_reply":"2021-07-02T09:23:17.70781Z"},"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":"def create_rgb(fname):\n    r1, r2, r3 = read_file(fname)\n    p1, p2, p3 = preprocess(r1, r2, r3)\n    hq1 = p1.q_transform(outseg=(0, 2))\n    hq2 = p2.q_transform(outseg=(0, 2))\n    hq3 = p3.q_transform(outseg=(0, 2))\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)","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:17.70951Z","iopub.execute_input":"2021-07-02T09:23:17.709933Z","iopub.status.idle":"2021-07-02T09:23:17.71719Z","shell.execute_reply.started":"2021-07-02T09:23:17.709886Z","shell.execute_reply":"2021-07-02T09:23:17.71632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"create_rgb('../input/g2net-gravitational-wave-detection/train/4/1/0/410031196b.npy')","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:17.718284Z","iopub.execute_input":"2021-07-02T09:23:17.718713Z","iopub.status.idle":"2021-07-02T09:23:18.409303Z","shell.execute_reply.started":"2021-07-02T09:23:17.71866Z","shell.execute_reply":"2021-07-02T09:23:18.40839Z"},"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":"def id2path(img_id, is_test):\n    a, b, c = img_id[0], img_id[1], img_id[2]\n    if is_test: return f'../input/g2net-gravitational-wave-detection/test/{a}/{b}/{c}/{img_id}.npy'\n    return f'../input/g2net-gravitational-wave-detection/train/{a}/{b}/{c}/{img_id}.npy'","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:18.41053Z","iopub.execute_input":"2021-07-02T09:23:18.410843Z","iopub.status.idle":"2021-07-02T09:23:18.415357Z","shell.execute_reply.started":"2021-07-02T09:23:18.410796Z","shell.execute_reply":"2021-07-02T09:23:18.414433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\ndf = pd.read_csv('../input/g2net-gravitational-wave-detection/sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:18.41689Z","iopub.execute_input":"2021-07-02T09:23:18.417204Z","iopub.status.idle":"2021-07-02T09:23:18.571318Z","shell.execute_reply.started":"2021-07-02T09:23:18.417174Z","shell.execute_reply":"2021-07-02T09:23:18.570259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nos.makedirs('test_images', exist_ok=True)","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:18.57286Z","iopub.execute_input":"2021-07-02T09:23:18.573265Z","iopub.status.idle":"2021-07-02T09:23:18.578196Z","shell.execute_reply.started":"2021-07-02T09:23:18.573224Z","shell.execute_reply":"2021-07-02T09:23:18.577199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def save_test_img(_id):\n    fname = id2path(_id, True)\n    im = create_rgb(fname)\n    im = im.resize((300,300), Image.BILINEAR)\n    im.save(f'test_images/{_id}.png', format=\"png\")","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:18.579589Z","iopub.execute_input":"2021-07-02T09:23:18.579933Z","iopub.status.idle":"2021-07-02T09:23:18.593322Z","shell.execute_reply.started":"2021-07-02T09:23:18.579858Z","shell.execute_reply":"2021-07-02T09:23:18.592631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = df[:56500]\nlen(df)","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:18.594456Z","iopub.execute_input":"2021-07-02T09:23:18.594892Z","iopub.status.idle":"2021-07-02T09:23:18.609148Z","shell.execute_reply.started":"2021-07-02T09:23:18.594847Z","shell.execute_reply":"2021-07-02T09:23:18.608429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/yasufuminakama/g2net-spectrogram-generation-train\nimport joblib\nfrom tqdm.auto import tqdm\n\n_ = joblib.Parallel(n_jobs=8)(\n    joblib.delayed(save_test_img)(_id) for _id in tqdm(df['id'].values)\n)\n","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:18.610203Z","iopub.execute_input":"2021-07-02T09:23:18.610644Z","iopub.status.idle":"2021-07-02T09:23:38.378521Z","shell.execute_reply.started":"2021-07-02T09:23:18.610599Z","shell.execute_reply":"2021-07-02T09:23:38.377564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import shutil\n\nshutil.make_archive('test_images', 'zip', 'test_images')\nshutil.rmtree('test_images')","metadata":{"execution":{"iopub.status.busy":"2021-07-02T09:23:38.379897Z","iopub.execute_input":"2021-07-02T09:23:38.380167Z","iopub.status.idle":"2021-07-02T09:23:41.726585Z","shell.execute_reply.started":"2021-07-02T09:23:38.380134Z","shell.execute_reply":"2021-07-02T09:23:41.725859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}