{"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":"While the use of Spectral Whitening did not lead to AUC improvments in my models I decided to document share this as experiments that fail are still valuable for the leassons that were learnt. \n\n## Spectral Whitening\n\nIn reviewing the examples given in [GWPy](https://gwpy.github.io/docs/stable/examples/spectrogram/spectrogram2.html) I noticed the use of a technique called \"spectral wightening\" that was central to detecting the GW signal in the noice. Based on this work I decided to try:\n1. Standardizing each signal, i.e., dividing by the max value of the signal.\n2. Apply spectral whitening\n3. Convert that to a spectrogram\n4. Train a CNN to recognise the \"chirp\"\n\nTo implement this I made the following choices:\n* For the CNN I used a pre-trained EfficientNet. To make this easier to code I used Keras\n* I first tried to use the spectrogram code in GWPy however I found that nnAudio ran about 100x faster, so like many other I switched to using nnAudio to create the spectrogram on the fly.\n* I tried to use the GWPy whitening function but this was too slow. As I could not find a whitening function in nnAudio or other PyTorch modules I wrote my own. First I had to understand what whitening was, then I had to code for it.\n\n\nI then ran this, using EffficientNetB3, and found that whitening the signal gave a 10% lower AUC.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom matplotlib import pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-14T06:09:40.955172Z","iopub.execute_input":"2021-07-14T06:09:40.955553Z","iopub.status.idle":"2021-07-14T06:09:40.959541Z","shell.execute_reply.started":"2021-07-14T06:09:40.955517Z","shell.execute_reply":"2021-07-14T06:09:40.958741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install gwpy> /dev/null","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-14T06:09:40.972379Z","iopub.execute_input":"2021-07-14T06:09:40.972903Z","iopub.status.idle":"2021-07-14T06:10:02.248294Z","shell.execute_reply.started":"2021-07-14T06:09:40.972858Z","shell.execute_reply":"2021-07-14T06:10:02.247132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# What a real GW event looks like\nLoosely following the example in [GWPy](https://gwpy.github.io/docs/stable/examples/spectrogram/spectrogram2.html) I read in a 2sec sample with known Gravity Wave event and create two spectrograms. \n1. The first, without whitening, the GW \"chirp\" can't be seen with the eye.\n2. The second, with whitening, clearly showsthe \"chip\" at the 0.4s mark starting at a low frequency (~30hx) and ramping up (~200hz) few hundrets of a second.\n \nLooking at how different these were I wondered if whitening could be used to improve my CNN model?","metadata":{"execution":{"iopub.status.busy":"2021-07-12T01:49:31.722147Z","iopub.execute_input":"2021-07-12T01:49:31.722468Z","iopub.status.idle":"2021-07-12T01:49:54.024843Z","shell.execute_reply.started":"2021-07-12T01:49:31.72244Z","shell.execute_reply":"2021-07-12T01:49:54.023636Z"}}},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\nlh = TimeSeries.fetch_open_data('H1', 1126259458+3.5, 1126259458+3.5+2,)\nspecgram = lh.spectrogram2(fftlength=1/16., overlap=15/256.) ** (1/2.)\nspecgram.plot(figsize=(8, 2));\nplt.show()\nspecgram = lh.whiten(2, 1).spectrogram2(fftlength=1/16., overlap=15/256.) ** (1/2.)\nspecgram.plot(figsize=(8, 2));\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-14T06:10:02.250487Z","iopub.execute_input":"2021-07-14T06:10:02.250966Z","iopub.status.idle":"2021-07-14T06:10:12.681731Z","shell.execute_reply.started":"2021-07-14T06:10:02.250919Z","shell.execute_reply":"2021-07-14T06:10:12.680708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Spectral Whitening\n\nMy first thought was just to use the GWPy tools however the spectrogram was 1 - 2 orders of magnitude slower than nnAudio. I switched to nnAudio and found I could generate the spectrograms of the fly, which was very convenient. \n\nHowever, I still had the problem that GWPy whitening  too slow so I set about building my own in torch to see how much faster this was. After much reading I implemented the following algorithm.","metadata":{}},{"cell_type":"code","source":"import torch\nfrom torch.fft import fft, rfft, ifft\n\ndef whiten(signal):\n    hann = torch.hann_window(len(signal), periodic=True, dtype=float)\n    spec = fft(torch.from_numpy(signal).float()* hann)\n    mag = torch.sqrt(torch.real(spec*torch.conj(spec))) \n\n    return torch.real(ifft(spec/mag)).numpy() * np.sqrt(len(signal)/2)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T06:10:12.683880Z","iopub.execute_input":"2021-07-14T06:10:12.684475Z","iopub.status.idle":"2021-07-14T06:10:13.942177Z","shell.execute_reply.started":"2021-07-14T06:10:12.684427Z","shell.execute_reply":"2021-07-14T06:10:13.941319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The charts of the whitened data below show the Torch version has similar but not quite the same output as GWPy. Comparing the GWPy [code](https://github.com/gwpy/gwpy/blob/v2.0.4/gwpy/timeseries/timeseries.py#L1669), the key difference is  in the last line my algorithm use multiplication in frequency-space and GWPy uses convolution with a hann window in time-space. I think its the use of the hann window that results a slighly different envelope. ","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(12, 4))\nplt.subplot(1, 3, 1)\nplt.plot(lh)\nplt.title(\"Raw Sample\")\nplt.subplot(1, 3, 2)\nplt.plot(lh.whiten(2,1))\nplt.title(\"GWPy Whitened\")\nplt.subplot(1, 3, 3)\nplt.title(\"nnAudio Whitened\")\nplt.plot(whiten(lh.value))\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-14T06:10:13.943603Z","iopub.execute_input":"2021-07-14T06:10:13.944076Z","iopub.status.idle":"2021-07-14T06:10:14.923636Z","shell.execute_reply.started":"2021-07-14T06:10:13.944023Z","shell.execute_reply":"2021-07-14T06:10:14.922401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What is most notable is the PYtorch is about 10x faster when using a CPU ","metadata":{}},{"cell_type":"code","source":"%timeit lh.whiten(2,0)\n%timeit whiten(lh.value)","metadata":{"execution":{"iopub.status.busy":"2021-07-14T06:10:14.925241Z","iopub.execute_input":"2021-07-14T06:10:14.925651Z","iopub.status.idle":"2021-07-14T06:10:29.715312Z","shell.execute_reply.started":"2021-07-14T06:10:14.925608Z","shell.execute_reply":"2021-07-14T06:10:29.714257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Comparing Spectrogram\nTo show that the new torch function works just as well, I create spectrograms of teh whitened sample using GWPy vs nnAudio/Torch","metadata":{}},{"cell_type":"code","source":"!pip install -q nnAudio","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-14T06:10:29.716698Z","iopub.execute_input":"2021-07-14T06:10:29.717068Z","iopub.status.idle":"2021-07-14T06:10:37.882321Z","shell.execute_reply.started":"2021-07-14T06:10:29.717031Z","shell.execute_reply":"2021-07-14T06:10:37.881172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from nnAudio.Spectrogram import CQT1992v2\ndef qgram(ts, \n           transform=CQT1992v2(sr=2048, fmin=20, fmax=2048/2, hop_length=64)): # tweak parameters to make a 64x64 image\n    image = transform(torch.from_numpy(ts).float()) # returns a tensor of spectrograms of shape = (num_samples, freq_bins,time_steps)\n    image = image.squeeze().numpy()\n\n    return image","metadata":{"execution":{"iopub.status.busy":"2021-07-14T06:10:37.885448Z","iopub.execute_input":"2021-07-14T06:10:37.885932Z","iopub.status.idle":"2021-07-14T06:10:37.925512Z","shell.execute_reply.started":"2021-07-14T06:10:37.885869Z","shell.execute_reply":"2021-07-14T06:10:37.924705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Again why not identical, we can see the tell tail signal of the 'chirp' in the spectrogram. However, when I used whitening in my modeling I found the AUC was substatially worse.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(12,7))\nplt.subplot(2,2,1)\nplt.imshow(lh.spectrogram2(fftlength=1/32., overlap=1/32-1/64))\nplt.title(\"GWPy Unwhitened\")\n\nplt.subplot(2,2,2)\nplt.imshow(lh.whiten(2,0).spectrogram2(fftlength=1/32., overlap=1/32-1/64.))\nplt.title(\"GWPy Whitened\")\n\nplt.subplot(2,2,3)\nplt.imshow(qgram(lh))\nplt.title(\"nnAudio Unwhitened\")\n\nplt.subplot(2,2,4)\nplt.imshow(qgram(whiten(lh)))\nplt.title(\"nnAudio Whitened\")\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-14T06:10:37.927447Z","iopub.execute_input":"2021-07-14T06:10:37.927928Z","iopub.status.idle":"2021-07-14T06:10:38.744514Z","shell.execute_reply.started":"2021-07-14T06:10:37.927884Z","shell.execute_reply":"2021-07-14T06:10:38.743741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Competition Data\nThe million dollar question is does this reveal the \"chirp\" signture for he competition data. \n\nIn the plots below we use the GWPy code, with the first column being the unwhitened spectrogram and the second whitened, but we don't get a dramatic result as we did with the real GW data.  Note the firs tthree are target =1 and the second three target =0","metadata":{}},{"cell_type":"code","source":"train_labels = pd.read_csv('../input/g2net-gravitational-wave-detection/training_labels.csv')\nsample_submission = pd.read_csv('../input/g2net-gravitational-wave-detection/sample_submission.csv')\n\ndef id2path(idx,is_train=True):\n    path = \"../input/g2net-gravitational-wave-detection\"\n    folder = 'train' if is_train else 'test'\n    return f'{path}/{folder}/{idx[0]}/{idx[1]}/{idx[2]}/{idx}.npy'","metadata":{"execution":{"iopub.status.busy":"2021-07-14T06:20:31.486848Z","iopub.execute_input":"2021-07-14T06:20:31.487557Z","iopub.status.idle":"2021-07-14T06:20:32.226531Z","shell.execute_reply.started":"2021-07-14T06:20:31.487516Z","shell.execute_reply":"2021-07-14T06:20:32.225510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idLst1 = train_labels[train_labels.target == 1].sample(3).values.tolist()\nidLst0 = train_labels[train_labels.target == 0].sample(3).values.tolist()\n\nfor _id, target in idLst1+idLst0:\n    wave = np.load(id2path(_id,is_train=True))[0]\n    ts  = TimeSeries(wave, sample_rate=2048)\n    plt.figure(figsize=(14, 5))\n    \n    plt.subplot(131)\n    plt.imshow(ts.spectrogram2(fftlength=1/32., overlap=1/32-1/64.))\n    plt.title(f\"Tgt={target} raw\")\n\n    plt.subplot(132)\n    plt.imshow(ts.whiten().spectrogram2(fftlength=1/32., overlap=1/32-1/64.))\n    plt.title(\"Whiten\")","metadata":{"execution":{"iopub.status.busy":"2021-07-14T07:32:23.458955Z","iopub.execute_input":"2021-07-14T07:32:23.459339Z","iopub.status.idle":"2021-07-14T07:32:26.439778Z","shell.execute_reply.started":"2021-07-14T07:32:23.459307Z","shell.execute_reply":"2021-07-14T07:32:26.438295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}