{"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 is still a **WIP**\n","metadata":{}},{"cell_type":"markdown","source":"# Before we start","metadata":{}},{"cell_type":"markdown","source":"In this competition, I have discovered some new signal processing transformations (CQT mainly)\nand have explored again ones I konw already (CWT, STFT, and so on). \n\nSo I have decided to take this opportunity as a refresher and share what I have learned.\n\nNotice that most of the implementations will be using PyTorch.\n\nLet's go!","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-20T18:28:40.758473Z","iopub.execute_input":"2021-08-20T18:28:40.758864Z","iopub.status.idle":"2021-08-20T18:28:40.770499Z","shell.execute_reply.started":"2021-08-20T18:28:40.758785Z","shell.execute_reply":"2021-08-20T18:28:40.769353Z"}}},{"cell_type":"markdown","source":"# Fourier Transform","metadata":{}},{"cell_type":"markdown","source":"Let's start with the king of transformations: the [Fourier Transform](https://en.wikipedia.org/wiki/Fourier_transform). \n\nLet's suppose we have a temporal signal $x(t)$, i.e. a sequence of values that change over time. \n\n\nFor that, let's use one of the competition's signal. Notice that we get \nthree signals (one from each observatory) for the price of one file. Neat!","metadata":{"execution":{"iopub.status.busy":"2021-08-20T18:31:28.447111Z","iopub.execute_input":"2021-08-20T18:31:28.447577Z","iopub.status.idle":"2021-08-20T18:31:28.459484Z","shell.execute_reply.started":"2021-08-20T18:31:28.447533Z","shell.execute_reply":"2021-08-20T18:31:28.457532Z"}}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pylab as plt\n\n\n# Any file will do here:\nx = np.load(\"../input/g2net-gravitational-wave-detection/train/0/0/0/00000e74ad.npy\")\n\n\n\nfig, axes = plt.subplots(x.shape[0], 1, figsize=(12, 8))\nfor i, ax in enumerate(axes):\n    ax.plot(x[i, :])\nfig.suptitle(f\"Time series\")\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see, some patterns seem to be hidden within the signal. So the next question is: how can we extract these efficiently?","metadata":{}},{"cell_type":"markdown","source":"The first \"clever\" thing that has been done is to try to decompose the temporal signal\ninto a sum of functions that varies with frequencies. \n\nTo find these functions, we will use the Fourier transform. \n\nLet's see how it is computed and implemented in practice now.","metadata":{}},{"cell_type":"markdown","source":"# Fourier transform: some theory","metadata":{"execution":{"iopub.status.busy":"2021-08-21T12:44:49.901387Z","iopub.execute_input":"2021-08-21T12:44:49.901732Z","iopub.status.idle":"2021-08-21T12:44:49.905631Z","shell.execute_reply.started":"2021-08-21T12:44:49.9017Z","shell.execute_reply":"2021-08-21T12:44:49.904465Z"}}},{"cell_type":"markdown","source":"Here is how it is defined for a signal $x(t)$: \n    \n    \n$\\hat{f}(\\xi) = \\int_{-\\infty}^{\\infty} x(t)\\ e^{-2\\pi i t \\xi}\\,dt$\n\n\nthe resulting function is a complex valued one in the frequency domain.","metadata":{}},{"cell_type":"markdown","source":"For the code part, we will use the PyTorch [FFT](https://pytorch.org/docs/stable/fft.html) module.","metadata":{}},{"cell_type":"code","source":"import torch\nfrom torch.fft import fft\n\n# The FFT results are complex numbers.\nfft_x = fft(torch.from_numpy(x[0]))\n\n# Plot the module. \n \nfig = plt.figure(figsize=(12, 8))\nplt.plot(np.abs(fft_x).reshape(-1))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# From Fourier transform to DFT and FFT","metadata":{}},{"cell_type":"markdown","source":"Next step, we need to move from Fourier transform to discrete Fourier transform (or DFT for short).","metadata":{}},{"cell_type":"markdown","source":"# What about STFT?","metadata":{}},{"cell_type":"markdown","source":"STFT stands for short-term Fourier transform. It is computed by taking moving\nwindows.\n\n\n\nHere is the forumula:  $ and here is how to compute it in more details: \n\n1. Take \n2. Compute the DFT\n    \n    \nFor the implementation, we will be using `librosa.stft`","metadata":{}},{"cell_type":"markdown","source":"# Spectrogram: time vs frequency","metadata":{}},{"cell_type":"markdown","source":"From the STFT, we get a spectrogram: $ \\mathrm{spectrogram}(t,\\omega)=\\left|\\mathrm{STFT}(t,\\omega)\\right|^2 $","metadata":{}},{"cell_type":"markdown","source":"For that we will use some code samples from\nthis script: https://pytorch.org/tutorials/beginner/audio_preprocessing_tutorial.html","metadata":{}},{"cell_type":"markdown","source":"Alright, now that we know how to get frequencies, it is time to plot them vs timestamps.","metadata":{}},{"cell_type":"markdown","source":"Spectrogram: https://en.wikipedia.org/wiki/Spectrogram","metadata":{}},{"cell_type":"markdown","source":"We can use nnaudio, librosa, or even PyTorch. Here is the code: \n\n\nhttps://pytorch.org/docs/stable/generated/torch.stft.html\n\nhttps://librosa.org/doc/latest/generated/librosa.stft.html\n\nhttps://kinwaicheuk.github.io/nnAudio/v0.1.5/_autosummary/nnAudio.Spectrogram.STFT.html","metadata":{}},{"cell_type":"code","source":"import librosa\nfrom librosa.display import specshow\nimport numpy as np\n\n\nstft_x = librosa.stft(x[0, :])\nspectrogram_x = np.abs(stft_x) ** 2\n\n\nfig, ax = plt.subplots(figsize=(15, 3))\n\n\n\n\nspectrogram_x = librosa.amplitude_to_db(spectrogram_x, ref=np.max)\n\nim = specshow(spectrogram_x, y_axis='log', x_axis='time', ax=ax)\n\nfig.colorbar(im, format='%+2.0f dB')\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not too  much to see. So we need better representations.","metadata":{}},{"cell_type":"markdown","source":"# CQT","metadata":{}},{"cell_type":"markdown","source":"Here is how it is defined and computed in few steps:\n\n1. We start by computing the short-term Fourier transform (STFT in short): this is a Fourier transform with a moving window over the signal. As a window function, we often use the Hann or Hamming windows. Finally, notice that we will use the discrete version of the STFT. The equation should be:\n\n$X[k,m] = \\sum_{n=0}^{N-1} W[n-m] x[n] e^{-j 2 \\pi k n/N}.$\n\n2. We introduce the quality factor Q using the following formula: $Q = \\frac{f_k}{\\delta f_k}.$\n\n3. Now, the window length for the k-th bin is no longer fixed (N) by will vary depending on the filter (we use the previous definition to get a simplified expression):\n\n$N[k] = \\frac{f_\\text{s}}{\\delta f_k} = \\frac{f_\\text{s}}{f_k} Q.$\n\n4. The digital frequency $\\frac{2 \\pi k}{N}$ is also changed and the window function now depends on $k$ as well (via the window length $N[k]$)\n\n5. Finaly, replacing the different elements in the original formula, we get:\n\n$X[k] = \\frac{1}{N[k]} \\sum_{n=0}^{N[k]-1} W[k,n] x[n] e^{\\frac{-j2 \\pi Qn}{N[k]}}. $","metadata":{}},{"cell_type":"markdown","source":"Now, let's see how it can be coded. ","metadata":{}},{"cell_type":"code","source":"!pip install nnAudio","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch\nfrom nnAudio.Spectrogram import CQT1992v2\nimport matplotlib.pylab as plt\nimport numpy as np\n\n# sr is the sampling rate, it is 2048 Hz\n# fmax is half the sampling rate\n\ncqt_transform = CQT1992v2(sr=2048, fmin=20, fmax=1024, hop_length=64)\n\ndef run_cqt_transform(x: np.array) -> torch.Tensor:\n    # We stack the passed x since there are 3\n    # time series per file.\n    x = np.hstack(x)\n    # Normalize (is there a better way?)\n    x = x / np.max(x)\n    x = torch.from_numpy(x).float()\n    return cqt_transform(x)\n\n\nimg = run_cqt_transform(x)[0]\n\nfig, ax = plt.subplots(figsize=(12, 8)) \nax.imshow(img)\nax.set_title(f\"CQT\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks much better: here at least we can see something!\n    \nCan we do better? Let's try another transformation and check.","metadata":{}},{"cell_type":"markdown","source":"# CWT","metadata":{}},{"cell_type":"markdown","source":"The next transformation won't strictly speaking based on Fourier transform but\non similar ideas. It uses [wavelets](https://en.wikipedia.org/wiki/Wavelet) as the basis. \n\nThis is the CWT, short for continuous wavelet transform.","metadata":{"execution":{"iopub.status.busy":"2021-08-23T06:21:11.198897Z","iopub.execute_input":"2021-08-23T06:21:11.199306Z","iopub.status.idle":"2021-08-23T06:21:11.211781Z","shell.execute_reply.started":"2021-08-23T06:21:11.199194Z","shell.execute_reply":"2021-08-23T06:21:11.210382Z"}}},{"cell_type":"markdown","source":"Using the following PyTorch library: https://github.com/tomrunia/PyTorchWavelets\n=> hard to install!!!","metadata":{}},{"cell_type":"code","source":"!git clone https://github.com/ar4/PyTorchWavelets.git > /dev/null\n%cd PyTorchWavelets\n!pip install -r requirements.txt > /dev/null\n!python setup.py install > /dev/null","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here is how to implement it: ","metadata":{}},{"cell_type":"code","source":"# There is a notebook that explains how to do it.\nfrom wavelets_pytorch.transform import WaveletTransform, WaveletTransformTorch# PyTorch version\n\n\ndt = 0.1         # sampling frequency\ndj = 0.125       # scale distribution parameter\n\n# Batch of signals to process\n# batch = [batch_size x signal_length]\n\n# Initialize wavelet filter banks (scipy and torch implementation)\nwa_scipy = WaveletTransform(dt, dj)\n# Doesn't work yet?\nwa_torch = WaveletTransformTorch(dt, dj, cuda=True)\n\n# Performing wavelet transform (and compute scalogram)\ncwt_torch = wa_torch.cwt(torch.tensor(x[0], device=\"cuda:0\").to(torch.float32))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cwt_torch","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Scalogram: same thing but for wavelets","metadata":{"execution":{"iopub.status.busy":"2021-08-21T12:16:03.164482Z","iopub.execute_input":"2021-08-21T12:16:03.164864Z","iopub.status.idle":"2021-08-21T12:16:03.169177Z","shell.execute_reply.started":"2021-08-21T12:16:03.16482Z","shell.execute_reply":"2021-08-21T12:16:03.167943Z"}}},{"cell_type":"markdown","source":"# Mel and MFCC","metadata":{}},{"cell_type":"markdown","source":"If you want even more details, check the following [notebook](https://www.kaggle.com/yassinealouini/what-is-mfcc) I made few months ago.\n\n\nIt explains in more details how Mel spectrogram and MFCCs are computed and how to use them. \n\nNotice that, in short, they are made from STFT as well and the big difference is that\nit is made from a spectrum of a spectrum (a [cepstrum](https://en.wikipedia.org/wiki/Cepstrum)).","metadata":{}},{"cell_type":"markdown","source":"# Whitening","metadata":{}},{"cell_type":"markdown","source":"One additional useful method is whitening of the signal. To do so.","metadata":{}},{"cell_type":"code","source":"# Use the following: from scipy.cluster.vq import whiten","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# What about PSD?","metadata":{}},{"cell_type":"markdown","source":"## Definition","metadata":{}},{"cell_type":"markdown","source":"[PSD](https://en.wikipedia.org/wiki/Spectral_density) is short for power spectral density.","metadata":{}},{"cell_type":"markdown","source":"## PSD implementation","metadata":{}},{"cell_type":"markdown","source":"Here are scipy methods: https://scipy-lectures.org/intro/scipy/auto_examples/plot_spectrogram.html#generate-a-chirp-signal\n\n\nhttps://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html#scipy.signal.welch\n\nLet's see what this gives to the original signal.\n\n\n","metadata":{}},{"cell_type":"code","source":"from scipy.signal import welch","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fs = 10e5\n\nf, Pxx_den = welch(x[0, :], fs, nperseg=1024)\n\nplt.semilogy(f, Pxx_den)\n\n\nplt.xlabel('frequency (Hz)')\n\nplt.ylabel('PSD ($V^{2}/Hz$)')\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# How to do this with PyTorch?","metadata":{}},{"cell_type":"markdown","source":"Here is how to implement PSD using PyTorch.","metadata":{"execution":{"iopub.status.busy":"2021-10-15T06:10:35.500764Z","iopub.execute_input":"2021-10-15T06:10:35.501259Z","iopub.status.idle":"2021-10-15T06:10:35.508631Z","shell.execute_reply.started":"2021-10-15T06:10:35.501215Z","shell.execute_reply":"2021-10-15T06:10:35.506934Z"}}}]}