{"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":"import warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n\nfrom IPython.display import HTML","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:31.527014Z","iopub.execute_input":"2021-09-01T19:34:31.527408Z","iopub.status.idle":"2021-09-01T19:34:31.537614Z","shell.execute_reply.started":"2021-09-01T19:34:31.527328Z","shell.execute_reply":"2021-09-01T19:34:31.536931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Signal, where are you?\n\nHey Klagger! This is going to be a difficult competition. Why? You just can't create a spectrogram and use a machine learning model to do the classification. You need to do a good preprocessing! And this will be difficult. We are looking for the needle in the haystack as our data contains not only the signal + some noise but also a lot of signals that belong to other sources like instruments used during the experiments. Even if this data is simulated we can expect that the signal is hidden!  \n","metadata":{}},{"cell_type":"code","source":"HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/FlDtXIBrAYE\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:31.539028Z","iopub.execute_input":"2021-09-01T19:34:31.539464Z","iopub.status.idle":"2021-09-01T19:34:31.554824Z","shell.execute_reply.started":"2021-09-01T19:34:31.539418Z","shell.execute_reply":"2021-09-01T19:34:31.553746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Have fun!","metadata":{}},{"cell_type":"markdown","source":"# Sources\n\nWith this notebook I try to collect information for myself to get started with this competition and I hope that you find it useful. My work is built upon the work of fantastic people and here you can find the sources:\n\n* https://www.youtube.com/channel/UC4oFlSYpDywInX0lxpiBPwA\n* https://arxiv.org/pdf/1908.11170.pdf\n* https://www.gw-openscience.org/LVT151012data/LOSC_Event_tutorial_LVT151012.html","metadata":{}},{"cell_type":"markdown","source":"# Table of contents\n\nTo find novel approaches to extract the signal of gravitational waves we need to understand the data itself and what kind of methods are usually used. Consequently my notebook is mainly about preprocessing and data understanding. ;-) \n\n* [Preparation](#preparation)\n* [Understanding the data](#understanding)\n    * How does the data look like?\n    * What kind of signals can be found?\n* [What's our goal for data preprocessing?](#preprocessing_goals)\n* [What kind of preprocessing methods are usually used?](#preprocessing_steps)\n* [Why are we using transformations like Fourier or Q-transform?](#transformations)\n    * [Waves](#waves)\n    * [Inference of waves](#inference)\n    * [Fourier series and the Fourier transform](#fourier_series)\n    * [Extracting signal components of your desire](#signal_extraction)\n    * [Sliding on the PSD - signal, where are you?](#comp_psd)\n* [What happens during each preprocessing step?](#whathappens)\n    * [Apply a window function](#windowing)\n    * [Whitening](#whitening)\n    * [Bandpass filtering](#bandpass_filtering)\n* [Searching for a signal](#searching_signal)\n* [Understanding the sources of gravitational waves](#gw_sources)\n\nLet's go!","metadata":{}},{"cell_type":"markdown","source":"# Preparation <a class=\"anchor\" id=\"preparation\"></a>\n\nOk, we need to load the packages...","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set()\nplt.rcParams[\"axes.grid\"] = False\n\nimport matplotlib.mlab as mlab\nfrom scipy import signal\nfrom scipy.interpolate import interp1d\nfrom scipy.signal import butter, filtfilt, iirdesign, zpk2tf, freqz\n\nfrom IPython.display import HTML","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:31.559901Z","iopub.execute_input":"2021-09-01T19:34:31.560313Z","iopub.status.idle":"2021-09-01T19:34:32.483472Z","shell.execute_reply.started":"2021-09-01T19:34:31.560277Z","shell.execute_reply":"2021-09-01T19:34:32.482305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"... and I would like to select an example to try out the preprocessing methods that are usually used. Let's load the data together with the path of the npy files:","metadata":{}},{"cell_type":"code","source":"train_labels = pd.read_csv(\"../input/g2net-trainlabelpaths/training_labels_with_paths.csv\")\ntrain_labels.target.value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:32.485098Z","iopub.execute_input":"2021-09-01T19:34:32.485396Z","iopub.status.idle":"2021-09-01T19:34:34.391380Z","shell.execute_reply.started":"2021-09-01T19:34:32.485367Z","shell.execute_reply":"2021-09-01T19:34:34.390416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As I'm trying to find a signal, let's use a hot example:","metadata":{}},{"cell_type":"code","source":"example = 1\n\nexample_strain = np.load(train_labels[train_labels.target==1].iloc[example].filepath)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:34.393573Z","iopub.execute_input":"2021-09-01T19:34:34.393899Z","iopub.status.idle":"2021-09-01T19:34:34.472515Z","shell.execute_reply.started":"2021-09-01T19:34:34.393867Z","shell.execute_reply":"2021-09-01T19:34:34.471541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Understanding the data <a class=\"anchor\" id=\"understanding\"></a>\n\n## How does the data look like?","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(20,5))\n\nplt.plot(example_strain[0,:], c=\"firebrick\", label=\"detector 1\")\nplt.plot(example_strain[1,:], c=\"mediumseagreen\", label=\"detector 2\")\nplt.plot(example_strain[2,:], c=\"slateblue\", label=\"detector 3\")\nplt.title(\"Example\");\nplt.legend();","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:34.474246Z","iopub.execute_input":"2021-09-01T19:34:34.474592Z","iopub.status.idle":"2021-09-01T19:34:34.787952Z","shell.execute_reply.started":"2021-09-01T19:34:34.474560Z","shell.execute_reply":"2021-09-01T19:34:34.786549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Insights\n\nOk, the three signals originating from different detectors all look a bit different. How was this data generated? We are given...\n\n* detector noise from three real detectors (LIGO Hanford, LIGO Livingston, and Virgo). As far as I understand this noise is not simulated. \n* A simulated gravitational wave signal hidden in this noisy data in the case of hot targets. You can't see it with your eyes! ","metadata":{}},{"cell_type":"markdown","source":"## What kind of signals can be found?\n\nWe can see a lot of waves in the data above? What's their origin? The data we observe was collected using a large-scale Michelson interferometer. Watch this great video here and you see how streching of the arms of this interferometer causes the waves we found:\n\nhttps://www.ligo.caltech.edu/video/IFO-response \n\nThe interferometer is sensitive towards gravitational waves but unfortunately also for terrestral forces and displacements. This may also include vibrations of the instruments themselves etc.. This kind of forces cause streching of the interferometer arms and this leads to constructive interferance and the waves we can see above. \n\nThis video is also great to understand how streching of the arms leads to constructive inference:","metadata":{}},{"cell_type":"code","source":"HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/tQ_teIUb3tE\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:34.789968Z","iopub.execute_input":"2021-09-01T19:34:34.790556Z","iopub.status.idle":"2021-09-01T19:34:34.798277Z","shell.execute_reply.started":"2021-09-01T19:34:34.790469Z","shell.execute_reply":"2021-09-01T19:34:34.796736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see that the waves can cancel each other out or merge to add up their amplitudes. As there are always spacial forces acting on the arms we end up with data showing the waves above. But be careful, they do not obviously show our signal!","metadata":{}},{"cell_type":"markdown","source":"# What's our goal for data preprocessing? <a class=\"anchor\" id=\"preprocessing_goals\"></a>\n\nWe are looking for something that is hidden in noise and may be difficult to extract. Take a look and listen to well extracted (real - not simulated) signals:","metadata":{}},{"cell_type":"code","source":"HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/QyDcTbR-kEA\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:34.800296Z","iopub.execute_input":"2021-09-01T19:34:34.800778Z","iopub.status.idle":"2021-09-01T19:34:34.809295Z","shell.execute_reply.started":"2021-09-01T19:34:34.800731Z","shell.execute_reply":"2021-09-01T19:34:34.808536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* The data we can see above was measured by the LIGO Hanford and LIGO Livingston detectors and both show the same event roughly starting at 0.9 s.\n* You can see the signal as a bright stripe in fourier space (time-frequency domain) and you can hear it like a chirp. \n* The amplitudes of our waves do not look like our data above as it was already preprocessed. ;-)\n\nEven though this event can be seen nicely in the spectogram we can't expect that we will always see our simulated waves as above. ","metadata":{}},{"cell_type":"markdown","source":"# What kind of preprocessing methods are usually used? <a class=\"anchor\" id=\"preprocessing_steps\"></a>\n\nTake a look at these two notebooks:\n\n* https://colab.research.google.com/github/losc-tutorial/quickview/blob/master/index.ipynb (with a clearly visible signal in the spectogram)\n* https://www.gw-openscience.org/LVT151012data/LOSC_Event_tutorial_LVT151012.html (without a clearly visible but present signal in the spectogram!)\n\nBasically we can see that the following preprocessing steps are applied:\n\n* **Whitening using the Amplitude Spectral Density** to eliminate the noise that can be seen in the ASD in the range of low and high frequencies as well as close to spectral lines. I hope that we (me included) will understand this better during the following chapters. \n* Using a **bandpass filter** to supress further low frequency signals.\n\nBefore we try to understand what these steps are doing, I need to refresh my knowledge! ;-)","metadata":{}},{"cell_type":"markdown","source":"# Why are we using transformations like Fourier or Q-transform? <a class=\"anchor\" id=\"transformations\"></a>","metadata":{}},{"cell_type":"markdown","source":"## Waves <a class=\"anchor\" id=\"waves\"></a>\n\nWhat is meant by low frequency noise? :-) Perhaps you didn't know so far, so let's talk a bit about oscillations! For example by looking at a sine wave:","metadata":{}},{"cell_type":"code","source":"x = np.arange(0,2*np.pi, step=0.1)\ny = np.sin(x)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:34.810666Z","iopub.execute_input":"2021-09-01T19:34:34.811140Z","iopub.status.idle":"2021-09-01T19:34:34.818864Z","shell.execute_reply.started":"2021-09-01T19:34:34.811096Z","shell.execute_reply":"2021-09-01T19:34:34.817892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(20,5))\nplt.plot(y, '-o', c=\"tomato\")\nplt.ylabel(\"Elongation y\")\nplt.xlabel(\"Time t\")\nplt.title(\"A sine wave\");","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:34.821448Z","iopub.execute_input":"2021-09-01T19:34:34.821986Z","iopub.status.idle":"2021-09-01T19:34:35.046489Z","shell.execute_reply.started":"2021-09-01T19:34:34.821942Z","shell.execute_reply":"2021-09-01T19:34:35.045753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* In the range from 0 to $2\\pi$ the sinus wave show us one oscillation like a pendulum going from the middle to the left to the middle to the right and back. \n* The time that has gone during this oscillation is called the period $T$. \n* And the inverse of the period is called the frequency of the wave: $f=\\frac{1}{T}$\n* Further we can define the length of the wave by $\\lambda = c \\cdot T = \\frac{c}{f}$. This is called the wavelength and c is the velocity of light.\n* There is one more spacial feature about waves - the elongation y which shows us how far the pendulum (for example) is away from the middle. Furthermore we can define the amplitude $\\hat{y}$ as the maximum of the elongation.\n\n**What is meant now by a low frequency wave?**\n\nA low frequency means a high wavelength. Consequently when we say that we want to eliminate low frequency parts of our data, this means that we want to get rid of the big oscillations with high wavelengths we can see very clearly in our data! ;-) ","metadata":{}},{"cell_type":"markdown","source":"# Inference of waves <a class=\"anchor\" id=\"inference\"></a>\n\nWe can add waves like the sine wave above to create new waves that are superpositions of the other ones. This is also called inference (so it's not the inference you do with a ML model). Let's try out a few examples!\n\nWe need a \"time\":","metadata":{}},{"cell_type":"code","source":"dt = 10/5000\nprint(dt)\nt = np.linspace(0, 10, 5000)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:35.047985Z","iopub.execute_input":"2021-09-01T19:34:35.048397Z","iopub.status.idle":"2021-09-01T19:34:35.054378Z","shell.execute_reply.started":"2021-09-01T19:34:35.048355Z","shell.execute_reply":"2021-09-01T19:34:35.052923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And frequencies:","metadata":{}},{"cell_type":"code","source":"f1 = 10\nf2 = 1000","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:35.055896Z","iopub.execute_input":"2021-09-01T19:34:35.056231Z","iopub.status.idle":"2021-09-01T19:34:35.066092Z","shell.execute_reply.started":"2021-09-01T19:34:35.056200Z","shell.execute_reply":"2021-09-01T19:34:35.064790Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And the sine to compute the elongations y for our waves:","metadata":{}},{"cell_type":"code","source":"y1 = 0.3 * np.sin(2*np.pi*f1*t)\ny2 = np.sin(2*np.pi*f2*t)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:35.067905Z","iopub.execute_input":"2021-09-01T19:34:35.068235Z","iopub.status.idle":"2021-09-01T19:34:35.077539Z","shell.execute_reply.started":"2021-09-01T19:34:35.068207Z","shell.execute_reply":"2021-09-01T19:34:35.076327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Don't wonder, the stuff inside the sine wave is a bit more complicated than what we used above. Instead of $y=sin(t)$ it holds $y=r \\cdot sin(\\omega  t)$ with $\\omega=2 \\pi f$ being [the angular frequency](https://en.wikipedia.org/wiki/Angular_frequency) and $r$ being the radius that we can set to adjust the amplitude $\\hat{y}$ (max elongation). You can see that it's just a constant times frequency. So there is nothing new to learn. A higher frequency means more ups and downs of our wave within a fixed time period compared to a low frequency. ","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(2,1,figsize=(20,8))\nax[0].plot(t, y1, label=\"high freqency\")\nax[0].plot(t, y2, label=\"low frequency\")\nax[1].plot(t, y1+y2)\nax[0].legend();\n\nfor n in range(2):\n    ax[n].set_xlabel(\"time t\")\n    ax[n].set_ylabel(\"elongation y\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:35.078893Z","iopub.execute_input":"2021-09-01T19:34:35.079240Z","iopub.status.idle":"2021-09-01T19:34:35.483940Z","shell.execute_reply.started":"2021-09-01T19:34:35.079169Z","shell.execute_reply":"2021-09-01T19:34:35.482650Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ok, with these values we obtain one wave with high frequency and low radius and another wave with low frequency and higher radius. Adding up both we can see the same kind of \"problem\" like above - a \"big wave\" with a small one \"hidden\". ","metadata":{}},{"cell_type":"markdown","source":"## Fourier series and the Fourier transform <a class=\"anchor\" id=\"fourier_series\"></a>\n\nTo understand what's going on during preprocessing, I like to go back to fourier series. With a fourier series we can try to approximate a periodic function f(x) with sine and cosine waves. This function can be our data composed of the signals given by terrestrical forces and perhaps also the signal of a simulated gravitational wave. But it can also be the function we have built above using two simple sine waves:\n\n$$ f(x) = \\frac{a_{0}}{2} + \\sum_{k=1}^{\\infty} \\left[ a_{k} \\cos \\left(\\frac{k \\pi}{l} x\\right) + b_{k} \\sin \\left(\\frac{k \\pi}{l} x\\right) \\right] + R_{n}(x)$$","metadata":{}},{"cell_type":"markdown","source":"The part without the residual $R_{n}(x)$ is called the fourier series. We have already seen that we can try to built a similar kind of \"a hidden singal travelling on a big wave\" by only adding up two sine waves. So I think it's a nice idea to say that we can try to built periodic, unknown functions by using a infinite sum of sine and cosine waves.","metadata":{}},{"cell_type":"markdown","source":"Hah! Found something great!! I love this following video! It really helps to get an idea about fourier series and fourier transform. I hope you like it too. :D","metadata":{}},{"cell_type":"code","source":"HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/spUNpyF58BY\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:35.485341Z","iopub.execute_input":"2021-09-01T19:34:35.485792Z","iopub.status.idle":"2021-09-01T19:34:35.492859Z","shell.execute_reply.started":"2021-09-01T19:34:35.485752Z","shell.execute_reply":"2021-09-01T19:34:35.491627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Even if you don't understand it fully, this video should give an impression of the basic idea why we are doing the fourier transform here. And that's all to get started. ;-)","metadata":{}},{"cell_type":"markdown","source":"## Extracting signal components of your desire <a class=\"anchor\" id=\"signal_extraction\"></a>\n\nI don't know how you feel but I'm very curious to try out the idea of the Fourier transform to find coefficients for the simple waves we built above. We should be able to easily filter out both types of waves using the Power Spectral Density. If you need another video to get started, this one here is great:","metadata":{}},{"cell_type":"code","source":"HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/s2K1JfNR7Sc\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:35.494551Z","iopub.execute_input":"2021-09-01T19:34:35.495010Z","iopub.status.idle":"2021-09-01T19:34:35.505655Z","shell.execute_reply.started":"2021-09-01T19:34:35.494952Z","shell.execute_reply":"2021-09-01T19:34:35.504618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ok, let's try it:\n\n* Perform the Fast Fourier Transform with $y_{1}$, $y_{2}$ and $y_{1}+y_{2}$\n* Compute the Power Spectral Density to get the components\n* Show that $y_{1}+y_{2}$ is like PSDs of $y_{1}$ and $y_{2}$ added up.","metadata":{}},{"cell_type":"code","source":"# our simple sine waves and the iference of both:\nto_take = [y1, y2, y1+y2, y1+y2+0.5*np.random.randn(len(t))]\ncolors = [\"darkslateblue\", \"tomato\", \"darkorchid\", \"crimson\"]\n\nfig, ax = plt.subplots(2,4,figsize=(20,8))\nfor n in range(4):\n    \n    # perform the fast fourier transform\n    fhat = np.fft.fft(to_take[n], len(t))\n    \n    # compute the power spectral density\n    PSD = fhat * np.conj(fhat) / len(t)\n    \n    freq = 1/(dt*len(t)) * np.arange(len(t))\n    L = np.arange(1, np.floor(len(t)/2), dtype=\"int\")\n\n    # plot it to show the components of our waves:\n    ax[1,n].plot(freq[L], PSD[L], c=colors[n])\n    #ax[1,n].set_ylim(0,1500)\n    ax[1,n].set_yscale(\"log\")\n    ax[1,n].set_xlabel(\"Frequency Hz\")\n    ax[1,n].set_ylabel(\"\")\n    \n    ax[0,n].set_xlabel(\"time t\")\n    ax[0,n].set_ylabel(\"elongation y\")\n    \n    ax[0,n].plot(t, to_take[n], c=colors[n])\n\nax[0,0].set_title(\"$y_{1}=0.3 * \\sin(2*\\pi*30*t)$\");\nax[0,1].set_title(\"$y_{2}=1 * \\sin(2*\\pi*1000*t)$\");\nax[0,2].set_title(\"$y_{1}+y_{2}$\");","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:35.507107Z","iopub.execute_input":"2021-09-01T19:34:35.507453Z","iopub.status.idle":"2021-09-01T19:34:37.441585Z","shell.execute_reply.started":"2021-09-01T19:34:35.507401Z","shell.execute_reply":"2021-09-01T19:34:37.440781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Insights\n\nCool, isn't it? :-) We can clearly see that the wave that was built by adding up the two sine waves has both components of them in its PSD. This is awesome as we can now use such kind of PSD components to filter out signals of our desire! To get closer to the competition I have added a noisy version of $y1+y2$ that shows a lot of fluctuations and high frequency components in the PSD.\n\nLet's try it out using the competition data!","metadata":{}},{"cell_type":"markdown","source":"## Sliding on the PSD - signal, where are you? <a class=\"anchor\" id=\"comp_psd\"></a>","metadata":{}},{"cell_type":"markdown","source":"Ok, we are given:","metadata":{}},{"cell_type":"code","source":"sample_rate = 2048 #Hz (1/seconds)\ntime_span = 2 #seconds\nsamples_total = time_span * sample_rate\ndt = 1/(samples_total) #4096 points in total\ndt\n\nchannel = 1 #pick which detector you like to use","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:37.442605Z","iopub.execute_input":"2021-09-01T19:34:37.442890Z","iopub.status.idle":"2021-09-01T19:34:37.447384Z","shell.execute_reply.started":"2021-09-01T19:34:37.442861Z","shell.execute_reply":"2021-09-01T19:34:37.446470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can use the same idea as above to compute the Power Spectral Density:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(20,5))\n\nfhat = np.fft.fft(example_strain[channel,:], samples_total)\nPSD = fhat * np.conj(fhat) / samples_total\nfreq = 1/(dt*samples_total) * np.arange(samples_total)\n\n\nL = np.arange(1, np.floor(samples_total/2), dtype=\"int\")\nplt.plot(freq[L],PSD[L], '.-')\nplt.grid(\"on\")\nplt.xlabel(\"Frequency Hz\");\nplt.title(\"Power spectral density\");\nplt.yscale(\"log\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:37.448884Z","iopub.execute_input":"2021-09-01T19:34:37.449301Z","iopub.status.idle":"2021-09-01T19:34:37.772102Z","shell.execute_reply.started":"2021-09-01T19:34:37.449257Z","shell.execute_reply":"2021-09-01T19:34:37.770937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Insights\n\n* We can see that we have more fluctuations in the PSD for high frequencies. :-( That's bad as there should be our hidden signal as well.\n* The low frequencies belong to our \"big waves\" that are caused by instrumental vibrations or terrestrical forces etc. Here we won't find our signal.\n\nLet's pick a few component examples (or regions in the PSD) to show which kind of wave belongs to it. This way we can compare with the original data and get some feeling about its structure:","metadata":{}},{"cell_type":"code","source":"show_side_effects = True","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:37.773729Z","iopub.execute_input":"2021-09-01T19:34:37.774119Z","iopub.status.idle":"2021-09-01T19:34:37.779055Z","shell.execute_reply.started":"2021-09-01T19:34:37.774078Z","shell.execute_reply":"2021-09-01T19:34:37.777877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(6,1,figsize=(20,15))\n\nax[0].plot(example_strain[channel])\nax[1].plot(np.fft.ifft((PSD>1e-38)*fhat))\nax[2].plot(np.fft.ifft(((PSD>1e-40) & (PSD <= 1e-38))*fhat))\nax[3].plot(np.fft.ifft(((PSD>1e-42) & (PSD <= 1e-40))*fhat))\nax[4].plot(np.fft.ifft(((PSD>0.5e-42) & (PSD <= 1e-42))*fhat))\nax[5].plot(np.fft.ifft((PSD<=0.5e-42)*fhat))\n\nif not show_side_effects:\n    for n in range(3,6):\n        ax[n].set_xlim(20,2000)","metadata":{"_kg_hide-input":true,"_kg_hide-output":false,"execution":{"iopub.status.busy":"2021-09-01T19:34:37.780633Z","iopub.execute_input":"2021-09-01T19:34:37.781026Z","iopub.status.idle":"2021-09-01T19:34:38.717741Z","shell.execute_reply.started":"2021-09-01T19:34:37.780992Z","shell.execute_reply":"2021-09-01T19:34:38.716658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Insights\n\n* We can clearly filter out different components of our original data. Nonetheless it's still hard to see something like a signal caused by a gravitational wave.\n* Can you see that we have this interesting kind of side effects in the beginning and the end of each wave after doing the inverse Fourier transform? ","metadata":{}},{"cell_type":"markdown","source":"# Apply a window function <a class=\"anchor\" id=\"windowing\"></a>\n\nWe need to apply a window function to suppress an effect called spectral leakage. [Here you can find a nice explanation about this kind of effect!](https://dspillustrations.com/pages/posts/misc/spectral-leakage-zero-padding-and-frequency-resolution.html#:~:text=Spectral%20leakage%20occurs%20when%20a,frequencies%20after%20the%20DFT%20operation.&text=For%20a%20single%2Dtone%20signal,even%20when%20spectral%20leakage%20occurs.) I will try to add an example when I have more time. ","metadata":{}},{"cell_type":"code","source":"hp_window = 1\nhp_tukey_alpha = 0.125\nfband = [35.0, 200.0]","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:38.719055Z","iopub.execute_input":"2021-09-01T19:34:38.719364Z","iopub.status.idle":"2021-09-01T19:34:38.724248Z","shell.execute_reply.started":"2021-09-01T19:34:38.719327Z","shell.execute_reply":"2021-09-01T19:34:38.722922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"blackman_window = signal.blackman(int(samples_total*hp_window)) #signal.tukey(strain, alpha=1./8)\ntukey_window = signal.tukey(samples_total*hp_window, hp_tukey_alpha)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:38.726171Z","iopub.execute_input":"2021-09-01T19:34:38.726629Z","iopub.status.idle":"2021-09-01T19:34:38.737587Z","shell.execute_reply.started":"2021-09-01T19:34:38.726575Z","shell.execute_reply":"2021-09-01T19:34:38.736438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\n\nax[0].plot(blackman_window, c=\"black\")\nax[0].set_title(\"Blackman window\")\n\nax[1].plot(tukey_window, c=\"purple\")\nax[1].set_title(\"Tukey window\")\n\nfor n in range(2):\n    ax[n].set_ylabel(\"Amplitude\")\n    ax[n].set_xlabel(\"Sample\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:38.739548Z","iopub.execute_input":"2021-09-01T19:34:38.739977Z","iopub.status.idle":"2021-09-01T19:34:39.083011Z","shell.execute_reply.started":"2021-09-01T19:34:38.739934Z","shell.execute_reply":"2021-09-01T19:34:39.081682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a look how the window functions change our data:","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(3,1,figsize=(20,15))\n\nax[0].plot(example_strain[channel])\nax[0].set_title(\"Original data\")\n\nax[1].plot(example_strain[channel]*blackman_window)\nax[1].set_title(\"With blackman window applied\")\n\nax[2].plot(example_strain[channel]*tukey_window)\nax[2].set_title(\"With tukey window applied\");","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:39.087825Z","iopub.execute_input":"2021-09-01T19:34:39.088146Z","iopub.status.idle":"2021-09-01T19:34:39.869328Z","shell.execute_reply.started":"2021-09-01T19:34:39.088116Z","shell.execute_reply":"2021-09-01T19:34:39.868091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Hmm... what if our signal can be found in the beginning our end of the data? Then windowing would be very bad... wouldn't it?!","metadata":{}},{"cell_type":"code","source":"windowed_strain = example_strain[channel]*tukey_window","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:39.871475Z","iopub.execute_input":"2021-09-01T19:34:39.871888Z","iopub.status.idle":"2021-09-01T19:34:39.876855Z","shell.execute_reply.started":"2021-09-01T19:34:39.871850Z","shell.execute_reply":"2021-09-01T19:34:39.875627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fhat = np.fft.fft(windowed_strain, samples_total)\nPSD = fhat * np.conj(fhat) / samples_total\nfreq = 1/(dt*samples_total) * np.arange(samples_total)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:39.878587Z","iopub.execute_input":"2021-09-01T19:34:39.878985Z","iopub.status.idle":"2021-09-01T19:34:39.888996Z","shell.execute_reply.started":"2021-09-01T19:34:39.878942Z","shell.execute_reply":"2021-09-01T19:34:39.888022Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"L = np.arange(1, np.floor(samples_total/2), dtype=\"int\")\n\nplt.figure(figsize=(20,5))\nplt.plot(freq[L],PSD[L], '.-')\nplt.grid(\"on\")\nplt.xlabel(\"Frequency Hz\");\nplt.title(\"Power spectral density\");\nplt.yscale(\"log\")","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:45:32.561349Z","iopub.execute_input":"2021-09-01T19:45:32.561956Z","iopub.status.idle":"2021-09-01T19:45:32.914020Z","shell.execute_reply.started":"2021-09-01T19:45:32.561917Z","shell.execute_reply":"2021-09-01T19:45:32.913257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(6,1,figsize=(20,15))\n\nax[0].plot(example_strain[channel])\nax[1].plot(np.fft.ifft((PSD>1e-38)*fhat))\nax[2].plot(np.fft.ifft(((PSD>1e-40) & (PSD <= 1e-38))*fhat))\nax[3].plot(np.fft.ifft(((PSD>1e-42) & (PSD <= 1e-40))*fhat))\nax[4].plot(np.fft.ifft(((PSD>0.5e-42) & (PSD <= 1e-42))*fhat))\nax[5].plot(np.fft.ifft((PSD<=0.5e-42)*fhat))\n\nif not show_side_effects:\n    for n in range(3,6):\n        ax[n].set_xlim(20,2000)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:40.230317Z","iopub.execute_input":"2021-09-01T19:34:40.230743Z","iopub.status.idle":"2021-09-01T19:34:41.149381Z","shell.execute_reply.started":"2021-09-01T19:34:40.230702Z","shell.execute_reply":"2021-09-01T19:34:41.148204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It seems that the side effects have gone away. Hmm.","metadata":{}},{"cell_type":"markdown","source":"# Whitening <a class=\"anchor\" id=\"whitening\"></a>\n\n\n\n* The data is dominated by low frequency noise (the large oscillations we can see! ;-))\n* Divide fourier coefficients by an estimate of the amplitude spectral density of the noise\n* This way we perform a down-weighting of the frequencies where the noise is loud\n* To return to the time-domain do inverse Fourier transform\n\nWhy are we doing it this way? Using spline interpolation... seems like something that can be easily replaced with alternatives.","metadata":{}},{"cell_type":"code","source":"def whiten(strain, samples_total, dt):\n    # TODO: normalization \n    \n    fhat = np.fft.fft(strain, samples_total)\n    PSD = fhat * np.conj(fhat) / samples_total\n    freq = 1/(dt*samples_total) * np.arange(samples_total)\n    \n    # scipy interp1d interpolation\n    interp_psd = interp1d(freq, PSD, \"nearest\")\n    \n    w_fhat = fhat/np.sqrt(interp_psd(freq))\n    w_strain = np.fft.ifft(w_fhat)\n    return w_strain, interp_psd(freq)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:41.151022Z","iopub.execute_input":"2021-09-01T19:34:41.151365Z","iopub.status.idle":"2021-09-01T19:34:41.158564Z","shell.execute_reply.started":"2021-09-01T19:34:41.151331Z","shell.execute_reply":"2021-09-01T19:34:41.157314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"w_strain, ip = whiten(windowed_strain, samples_total, dt)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:41.160187Z","iopub.execute_input":"2021-09-01T19:34:41.160587Z","iopub.status.idle":"2021-09-01T19:34:41.175871Z","shell.execute_reply.started":"2021-09-01T19:34:41.160553Z","shell.execute_reply":"2021-09-01T19:34:41.174760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(2,1,figsize=(20,10))\nax[0].plot(np.log(ip[0:1024]), '-o')\nax[0].set_title(\"Interpolated PSD\")\nax[0].set_xlabel(\"Hz\")\nax[1].plot(w_strain, '-.')","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:41.177578Z","iopub.execute_input":"2021-09-01T19:34:41.177909Z","iopub.status.idle":"2021-09-01T19:34:41.617308Z","shell.execute_reply.started":"2021-09-01T19:34:41.177876Z","shell.execute_reply":"2021-09-01T19:34:41.616329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Still looks very noisy.","metadata":{"execution":{"iopub.status.busy":"2021-09-01T09:34:12.28034Z","iopub.execute_input":"2021-09-01T09:34:12.280674Z","iopub.status.idle":"2021-09-01T09:34:12.286137Z","shell.execute_reply.started":"2021-09-01T09:34:12.280643Z","shell.execute_reply":"2021-09-01T09:34:12.284995Z"}}},{"cell_type":"markdown","source":"# Apply a bandpass filter <a class=\"anchor\" id=\"bandpass_filtering\"></a>\n\n* https://en.wikipedia.org/wiki/Butterworth_filter","metadata":{}},{"cell_type":"code","source":"def bandpass(strain, fband, fs):\n    \"\"\"Bandpasses strain data using a butterworth filter.\n    \n    Args:\n        strain (ndarray): strain data to bandpass\n        fband (ndarray): low and high-pass filter values to use\n        fs (float): sample rate of data\n    \n    Returns:\n        ndarray: array of bandpassed strain data\n    \"\"\"\n    bb, ab = butter(4, [fband[0]*2./fs, fband[1]*2./fs], btype='band')\n    normalization = np.sqrt((fband[1]-fband[0])/(fs/2))\n    strain_bp = filtfilt(bb, ab, strain) / normalization\n    return strain_bp","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:41.618578Z","iopub.execute_input":"2021-09-01T19:34:41.618871Z","iopub.status.idle":"2021-09-01T19:34:41.625845Z","shell.execute_reply.started":"2021-09-01T19:34:41.618842Z","shell.execute_reply":"2021-09-01T19:34:41.624549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bandpassed_strain = bandpass(w_strain, fband, samples_total)","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:41.627257Z","iopub.execute_input":"2021-09-01T19:34:41.627610Z","iopub.status.idle":"2021-09-01T19:34:41.643812Z","shell.execute_reply.started":"2021-09-01T19:34:41.627578Z","shell.execute_reply":"2021-09-01T19:34:41.642468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(20,5))\nplt.plot(bandpassed_strain, '-')","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:34:41.645814Z","iopub.execute_input":"2021-09-01T19:34:41.646223Z","iopub.status.idle":"2021-09-01T19:34:41.863049Z","shell.execute_reply.started":"2021-09-01T19:34:41.646190Z","shell.execute_reply":"2021-09-01T19:34:41.861912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ok, here is still something not ok...","metadata":{}},{"cell_type":"markdown","source":"# Understanding the sources of gravitational waves\n\n\n\n","metadata":{}},{"cell_type":"code","source":"HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/p__5MuhvWK0\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-01T19:34:41.864889Z","iopub.execute_input":"2021-09-01T19:34:41.865308Z","iopub.status.idle":"2021-09-01T19:34:41.872530Z","shell.execute_reply.started":"2021-09-01T19:34:41.865268Z","shell.execute_reply":"2021-09-01T19:34:41.871462Z"},"trusted":true},"execution_count":null,"outputs":[]}]}