{"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":"# Introduction to Signal Processing for Absolute Beginners\n\nIf you have no idea what spectrogram or frequency mean, and Fourier transform sounds like a magic spell, then you're like me! Let's get on a journey together to find out how to process signals, using G2Net competition data as examples.\n\nI will work on this notebook over the next several weeks, so you may come back when it's ready or join me on the learning journey. If you follow along, I'll appreciate comments with feedback or questions!","metadata":{"execution":{"iopub.status.busy":"2021-07-03T07:19:58.600322Z","iopub.execute_input":"2021-07-03T07:19:58.601184Z","iopub.status.idle":"2021-07-03T07:19:58.656806Z","shell.execute_reply.started":"2021-07-03T07:19:58.60106Z","shell.execute_reply":"2021-07-03T07:19:58.654923Z"}}},{"cell_type":"markdown","source":"![Signal Processing System](https://upload.wikimedia.org/wikipedia/commons/4/46/Signal_processing_system.png)","metadata":{}},{"cell_type":"markdown","source":"## Plan of Attack\n\nHow do we learn about this from scratch? Here's my plan:\n\n1. Start by tinkering! I already did this with a very naive approach first, and a some more informed but still very blind experiments later.\n\n2. Learn the basics! Let's understand some of the theory to figure out what we missed in our naive approach from the previous step, and come up with a more informed approach. \n\n3. Apply the learning in practice! After getting some understanding, let's put it to use to get some practise and verify if it works!","metadata":{}},{"cell_type":"markdown","source":"## What is Digital Signal Processing? \n\nI started by reading the first chapter from *Foundations of Digital Signal Processing* by Patrick Gaydecki (via google books). Here some things I have learned: \n\n* Digital Signal Processing involves manipulation of signals that have their origins in the analog world, such as audio, video, radar, thermal, magnetic or ultrasonic sensor systems\n* Signals are originally analog and continuous in nature\n* A transducer is a device that converts some form of energy to which it is designed to respond into an electrical signal. This signal is later converted into digital form. \n\nI also learned about some of the *founding fathers* of signal processing: \n* Jean-Baptiste Fourier (1768 - 1830) showed how any waveform could be represented as a series of weighted sine and cosine components of ascending frequency, which set the foundation for modern Fourier analysis and the Fourier transform. This allows us to analyse and process a signal. \n* Laplace developed the Laplace transform (not sure yet what that is)\n* Claud Shannon developed the science of information theory\n* Cooley and Turkey developed the FFT (Fast Fourier transform) in 1965, which allowed to obtain spectra in an efficient way.\n\nMy conclusion is that I need to take another step back, to high school physics, to understand more about waves. ","metadata":{}},{"cell_type":"markdown","source":"## What About Waves? \n\nIt seems that waves are a core concept necessary to understand signal processing. I can highly recommend Khan Academy videos such as this [introduction to waves](https://www.khanacademy.org/science/high-school-physics/waves-and-sound/introduction-to-waves/v/introduction-to-waves). Here are some things I have learned. \n\n![basic wave concepts](https://pbs.twimg.com/media/E5lsNiQWQAQ_Fq6?format=jpg&name=4096x4096)\n\nA **wave** is a disturbance propagating through space, usually transferring energy. In the example pictured above, imagine someone jerking a string up, down, and back up again. This is an example of *transverse* wave, because the medium (string) is moving up and down. If we repeat that movement and create a periodic *pulse* then we will create a *periodic wave*. \n\nThere are a few properties of *periodic waves* that we should know about: \n* **Amplitude** is the height of a wave, the maximum displacement measured from the equilibrium position\n* **Period** measures the number of seconds per each cycle\n* **Frequency** is the reverse of period, it counts how many cycles can fit in a second. The unit of frequency is **Hertz**, for example 10 cycles per second = 10 Hertz\n* **Wavelength** measures how far a wave has travelled after 1 period (1 cycle). This can be for example the distance between 2 neighboring peaks. \n* **Velocity** tells us how quickly a wave is moving to the right. Velocity equals distance over time, so we can say V = wavelength / period = wavelength * frequency\n\n### 1. Sound\n\nThe concepts above apply to sound, although it's a different type of wave. Let's start by considering how sound is produced via a loudspeaker...\n\n![sound waves](https://pbs.twimg.com/media/E5nmOMoXoBcnguC?format=jpg&name=large)\n\nA loudspeaker moves a diaphragm back and forth (oscillates), which pushes the air molecules to follow a similar back and forth movement. The molecules don't just move forward (like a ball that is kicked), but they move forward and backward like a wave. That wave propagates through air and we can see the compressed region of molecules moving forward. In contrast to the string example above, the movement is in the direction of propagated pulse, so we call this a *longitudinal wave*. It can be mapped on a diagram in a very similar way to *transverse ways* and the properties described above apply here as as well. \n\n* **Amplitude** measures the displacement of air molecule as it oscillates. Higher amplitude means higher volume, and lower amplited means lower volume. \n* **Frequency** corresponds to the sound pitch. Higher frequency corresponds to higher notes. Humans hear frequencies between 20Hz and 20,000Hz. \n* **Wavelength** measures the distance between two compressed regions of air. \n\nWe use **decibel scale** to measure loudness of sounds. To calculate decibel level, we need to know the Intensity of sound, which corresponds to Power (measured in Watts) over Area (measured in square meters). We divide Intensity by the threshold of hearing (1e-12 W/m2), take a log of this, and multiply by 10.\n\n#### Ultrasound \n\nSound waves can be useful beyond communication and entertainment! For example, sound waves with frequency higher than human hearing (above 20,000Hz) are used for medical imaging. A device is emiting high-frequency waves and captures the waves reflecteced by our internal organs to produce an image of our muscles, tendons and many internal organs. \n\n\n### 2. Electromagnetic Waves\n\nWhen an electric field changes in some region, that creates a changing magnetic field. This in turn leads to a change in the electric field. This can result in a chain reaction that propagates like a wave. These waves are called **electromagnetic waves**.\n\nIn the image below, you can see how an electromagnetic wave contains both electric and magnetic field vectors, orthogonal to each other, and how we can apply the concepts of amplitude, wavelength and frequency to this type of waves. An interesting fact about electromagnetic waves is that they can propagate in vacuum and don't require a medium. \n\n![](https://cdn.kastatic.org/ka-perseus-images/9b999f75e599f3ed46c0ed16586410ec1e5ffc34.png)\n\nBecause they have a constant speed (speed of light), there is a direct relationship between wavelenght and frequency. If we map all the possible frequencies/wavelengths on a chart, we call this electromagnetic spectrum. The visible light is a small part of that spectrum. Lower frequencies correspond to infrared, radio waves, microwaves etc. Higher frequencies correspond to ultraviolet, X-rays, gamma rays etc. \n\n![](https://cdn.kastatic.org/ka-perseus-images/7370593cc71daa2ccaca091cec088fa5fec6ca16.png)\n\n### 3. Gravitational Waves\n\nFinally, time to cover the waves we're supposed to detect in the G2Net competition! According to Wikipedia, \"**gravitational waves** are disturbances in the curvature of spacetime, generated by accelerated masses, that propagate as waves outward from their source at the speed of light\". \n\n![](https://upload.wikimedia.org/wikipedia/commons/thumb/2/2f/The_Gravitational_wave_spectrum_Sources_and_Detectors.jpg/2560px-The_Gravitational_wave_spectrum_Sources_and_Detectors.jpg)\n\n### Wavelets\n\nWhat happens if an oscillation is only temporary and doesn't display the periodic pattern, but still looks like a wave? You might have seen these from a seismograph recording or a heart monitor. Such brief oscillation is a **wavelet** and we can also analyse and use them with our digital signal processing tools. Here is an example of a wavelet. \n\n![wavelet example](https://upload.wikimedia.org/wikipedia/commons/thumb/0/0a/MorletWaveletMathematica.svg/1920px-MorletWaveletMathematica.svg.png)\n","metadata":{}},{"cell_type":"markdown","source":"## Frequency Domain\n\nI've been struggling to understand the concept of frequency domain - if you're like me, don't worry if this takes time to understand! Let's start with a rainbow. Rainbow shows us how white light can be decomposed into colors. Each color is associated with a different frequency. \n\n![rainbow](https://upload.wikimedia.org/wikipedia/commons/5/5c/Double-alaskan-rainbow.jpg)\n\nThis idea that a **signal can be a mix of different frequencies** is why we need the frequency domain to better understand a signal. And this is where our magic spell - **the Fourier transform** - comes in handy: it allows us to decompose any signal into a set of sine and cosine functions.\n\nSide note: humans can also analyse sounds and images in terms of their sinusoidal components, our ears and eyes have developed to distinguish different frequencies. \n\n#### Frequency vs. time domain\n\nSince we are dealing with digital representations of signals, we can obtain measurements at discrete points in time. A **sampling rate** tells us the number of samples per second (1/s unit is Hertz so we will use it here). For example, a typical sampling rate for audio signal is 44.1 kHz. \n\nSide note: if you fancy some more theory, check out [Nyquist–Shannon sampling theorem](https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem).\n\nIf we take a single sample (point it time), it's not possible to approximate the sinusoidal frequency function. The more points we have, the better the approximation. That's why we're talking about moving from time domain to frequency domain (Fourier analysis) or from frequency domain to time domain (Fourier synthesis). The analysis allows us to find the contribution of different frequencies in a given signal. The synthesis helps us create signals with known frequency content. \n\n### The Spectrogram\n\nOk, so we just learned that we can move from time domain to frequency domain with the help of Fourier transform. But wouldn't it make sense to see both time and frequency at the same time? Frequency can change over time, right? \n\nHere comes the **spectrogram**! The way to create it is by performing FFT on segments of the signal with a moving window. This is called **short-time Fourier transform**. The picture below illustrates it pretty well. \n\n![stft](https://www.researchgate.net/publication/346243843/figure/fig1/AS:961807523000322@1606324191138/Short-time-Fourier-transform-STFT-overview.png)\n\nFinally, we're at a point where we can write some code!","metadata":{"execution":{"iopub.status.busy":"2021-07-08T05:26:23.767447Z","iopub.execute_input":"2021-07-08T05:26:23.767995Z","iopub.status.idle":"2021-07-08T05:26:23.777736Z","shell.execute_reply.started":"2021-07-08T05:26:23.767888Z","shell.execute_reply":"2021-07-08T05:26:23.77672Z"}}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.mlab as mlab\nfrom scipy.interpolate import interp1d\nfrom scipy.signal import butter, filtfilt, iirdesign, zpk2tf, freqz\n\ndef id2path(id, is_test=False):\n    a, b, c = id[0], id[1], id[2]\n    if is_test: return f'../input/g2net-gravitational-wave-detection/test/{a}/{b}/{c}/{id}.npy'\n    return f'../input/g2net-gravitational-wave-detection/train/{a}/{b}/{c}/{id}.npy'","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-12T15:34:59.175790Z","iopub.execute_input":"2021-07-12T15:34:59.176158Z","iopub.status.idle":"2021-07-12T15:34:59.182722Z","shell.execute_reply.started":"2021-07-12T15:34:59.176127Z","shell.execute_reply":"2021-07-12T15:34:59.181808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's first look at our raw signal data. ","metadata":{}},{"cell_type":"code","source":"_id = '0021f9dd71' # credit for finding example with strong signal: https://www.kaggle.com/mistag/data-preprocessing-with-gwpy/\nx = np.load(id2path(_id, is_test=True))\nplt.figure(figsize=(12,4))\nplt.plot(x[0])\nplt.title('Detector 1')\nplt.xlabel('Time')\nplt.ylabel('Amplitude')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-12T15:25:46.577979Z","iopub.execute_input":"2021-07-12T15:25:46.578475Z","iopub.status.idle":"2021-07-12T15:25:46.747839Z","shell.execute_reply.started":"2021-07-12T15:25:46.578355Z","shell.execute_reply":"2021-07-12T15:25:46.746630Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's create a spectrogram by following [this tutorial](https://www.gw-openscience.org/GW150914data/GW150914_tutorial.html). We will start with a regular version, then we will apply whitening. It looks like whitening might not be helping though ([discussion](https://www.kaggle.com/c/g2net-gravitational-wave-detection/discussion/252138)).","metadata":{}},{"cell_type":"code","source":"fs = 2048 # sampling rate:\nNFFT = int(fs/16) # pick a shorter FTT time interval, like 1/16 of a second:\nNOVL = int(NFFT*15/16) # and with a lot of overlap, to resolve short-time features:\n# and choose a window that minimizes \"spectral leakage\" (https://en.wikipedia.org/wiki/Spectral_leakage)\nwindow = np.blackman(NFFT)\nspec_cmap='viridis'\nprint(f'fs: {fs}, NFFT: {NFFT}, NOVL: {NOVL}')\n\nplt.figure(figsize=(12,8))\nspec_H1, freqs, bins, im = plt.specgram(x[0], NFFT=NFFT, Fs=fs, window=window, noverlap=NOVL, cmap=spec_cmap)\nplt.xlabel('time')\nplt.ylabel('Frequency (Hz)')\nplt.axis([0.03, 2-0.03, 30, 500]) # according to discussions, the frequencies that interest us should be in this range\nplt.colorbar()\nplt.title('Detector 1')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-12T15:48:50.784648Z","iopub.execute_input":"2021-07-12T15:48:50.784977Z","iopub.status.idle":"2021-07-12T15:48:51.165185Z","shell.execute_reply.started":"2021-07-12T15:48:50.784948Z","shell.execute_reply":"2021-07-12T15:48:51.164229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can now see a gravitational wave signal on the spectrogram!","metadata":{}},{"cell_type":"markdown","source":"### Sources\n\n1. *Foundations of Digital Signal Processing* by Patrick Gaydecki (via google books)\n2. Khan Academy videos about waves\n3. MOOCs on Coursera: \n    *  https://www.coursera.org/learn/audio-signal-processing\n    *  https://www.coursera.org/learn/dsp1\n4. Preprocessing tutorials: https://www.kaggle.com/c/g2net-gravitational-wave-detection/discussion/250244","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}