{"cells":[{"metadata":{},"cell_type":"markdown","source":"##                                                        Introduction\n### In this kernel, I will walk you through some methods used to denoise seismic signals (and other signals in general) and explain exactly how they work. We will look at <a href='#1'>Wavelet Denoising with High-Pass Filters</a> and <a href='#1'>Average Smoothing</a>. The first method is used to remove the artificial impulse and the latter one is used to remove general noise.\n\n<center><img src=\"https://i.imgur.com/hBPv3fh.png\" width=\"750px\"></center>\n\nThe main problem in the specific case of seismic signals is the fact that the signal we measure with a seismograph is not an accurate representation of the actual underground seismic signal we are trying to uncover. **In seismology, we (the people trying to measure the seismic signals) artificially generate signals called impulse signals. These impulse signals interact with the Earth's actual seismic signal (which is what we need) to produce the final signal which our seismograph picks up** (this same process takes place in the laboratory simulation of an earthquake). So, the real challenge is to uncover the actual signal from mixed seismogram (which is a combination of the Earth's impulse signal and the artificial impulse signal).\n\nThis actual underlying signal would be a better predictor of earthquake timing than the original raw signal, because it represents the actual seismic activity."},{"metadata":{},"cell_type":"markdown","source":"<center><img src=\"https://qph.fs.quoracdn.net/main-qimg-4033b6af35e7154e6b497adf0d93c2d9-c\" width=\"300px\"></center>"},{"metadata":{},"cell_type":"markdown","source":"The third signal from the right is what we have (with artificial impulse) and the last signal from the right is what we're looking for (without the artificial impulse)."},{"metadata":{},"cell_type":"markdown","source":"## Acknowledgements\nI would like to thank Jack for [his kernel](https://www.kaggle.com/jackvial/dwt-signal-denoising) and Theo for [his kernel](https://www.kaggle.com/theoviel/fast-fourier-transform-denoising) (both on denoising)."},{"metadata":{},"cell_type":"markdown","source":"### Import necessary libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport gc\nimport numpy as np\nfrom numpy.fft import *\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nimport pywt \nfrom statsmodels.robust import mad\nimport scipy\nfrom scipy import signal\nfrom scipy.signal import butter, deconvolve\nimport warnings\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Specify size of seismic signal segements and rate of sampling for high-pass filter"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"SIGNAL_LEN = 150000\nSAMPLE_RATE = 4000","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Read data"},{"metadata":{"trusted":true},"cell_type":"code","source":"seismic_signals = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Extract seismic data and targets and delete the original dataframe"},{"metadata":{"trusted":true},"cell_type":"code","source":"acoustic_data = seismic_signals.acoustic_data\ntime_to_failure = seismic_signals.time_to_failure\ndata_len = len(seismic_signals)\ndel seismic_signals\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Cut out segments of size 150k from the signal"},{"metadata":{"trusted":true},"cell_type":"code","source":"signals = []\ntargets = []\n\nfor i in range(data_len//SIGNAL_LEN):\n    min_lim = SIGNAL_LEN * i\n    max_lim = min([SIGNAL_LEN * (i + 1), data_len])\n    \n    signals.append(list(acoustic_data[min_lim : max_lim]))\n    targets.append(time_to_failure[max_lim])\n    \ndel acoustic_data\ndel time_to_failure\ngc.collect()\n    \nsignals = np.array(signals)\ntargets = np.array(targets)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The mean absolute deviation\nThis calculates the mean of the absolute values of the deviations of the individual numbers in the time series from the mean of the time series. It is a measure of entropy or disorder in the time series. The greater the MAD value, the more disorderly and unpredictable the time series is."},{"metadata":{"trusted":true},"cell_type":"code","source":"def maddest(d, axis=None):\n    \"\"\"\n    Mean Absolute Deviation\n    \"\"\"\n    \n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The butterworth high-pass filter with SOS filter\nA high-pass filter (HPF) is an electronic filter that passes signals with a frequency higher than a certain **cutoff frequency** and attenuates (reduces the amplitude of) signals with frequencies lower than the cutoff frequency. The amount of attenuation for each frequency depends on the filter design.\n\nAn SOS filter implements a usual [IIR filter](https://en.wikipedia.org/wiki/Infinite_impulse_response) on the [second-order sections](https://edoras.sdsu.edu/doc/matlab/toolbox/signal/basics27.html#973) of the signal.\n\nThese two filters are applied in succession in order to attenuate the unnecessary low frequency signals and prepare the signal for wavelet decomposition."},{"metadata":{"trusted":true},"cell_type":"code","source":"def high_pass_filter(x, low_cutoff=1000, SAMPLE_RATE=SAMPLE_RATE):\n    \"\"\"\n    From @randxie https://github.com/randxie/Kaggle-VSB-Baseline/blob/master/src/utils/util_signal.py\n    Modified to work with scipy version 1.1.0 which does not have the fs parameter\n    \"\"\"\n    \n    # nyquist frequency is half the sample rate https://en.wikipedia.org/wiki/Nyquist_frequency\n    nyquist = 0.5 * SAMPLE_RATE\n    norm_low_cutoff = low_cutoff / nyquist\n    \n    # Fault pattern usually exists in high frequency band. According to literature, the pattern is visible above 10^4 Hz.\n    sos = butter(10, Wn=[norm_low_cutoff], btype='highpass', output='sos')\n    filtered_sig = signal.sosfilt(sos, x)\n\n    return filtered_sig","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Wavelet Denoising\nThe filtered signals are now passed through [wavelet](https://en.wikipedia.org/wiki/Wavelet) decomposition (and the wavelet coefficients are obtained). It has to be understood that the raw signal we are working with is a convolution of the artificial and real impulse signals and this is why we need to \"wavelet decomposition\". The process is a sort of \"deconvolution\", which means it undoes the convolution process and uncovers the real Earth impulse from the mixed seismogram and the artificial impulse. We make use of the MAD value to understand the randomness in the signal and accordingly decide the minimum threshold for the wavelet coefficients in the time series. We filter out the low coefficients from the wavelet coefficients and reconstruct the real Earth signal from the remaining coefficients and that's it; we have successfully removed the impulse signal from the seismogram and obtained the real Earth signal."},{"metadata":{"trusted":true},"cell_type":"code","source":"def denoise_signal(x, wavelet='db4', level=1):\n    \"\"\"\n    1. Adapted from waveletSmooth function found here:\n    http://connor-johnson.com/2016/01/24/using-pywavelets-to-remove-high-frequency-noise/\n    2. Threshold equation and using hard mode in threshold as mentioned\n    in section '3.2 denoising based on optimized singular values' from paper by Tomas Vantuch:\n    http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    \"\"\"\n    \n    # Decompose to get the wavelet coefficients\n    coeff = pywt.wavedec(x, wavelet, mode=\"per\")\n    \n    # Calculate sigma for threshold as defined in http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    # As noted by @harshit92 MAD referred to in the paper is Mean Absolute Deviation not Median Absolute Deviation\n    sigma = (1/0.6745) * maddest(coeff[-level])\n\n    # Calculate the univeral threshold\n    uthresh = sigma * np.sqrt(2*np.log(len(x)))\n    coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])\n    \n    # Reconstruct the signal using the thresholded coefficients\n    return pywt.waverec(coeff, wavelet, mode='per')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualize the effect of high-pass filter and wavelet denoising on seismic signals"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=3, ncols=4, figsize=(30, 15))\n\nax[0, 0].plot(signals[0], 'crimson') \nax[0, 0].set_title('Original Signal', fontsize=16)\nax[0, 1].plot(high_pass_filter(signals[0], low_cutoff=10000, SAMPLE_RATE=4000000), 'mediumvioletred') \nax[0, 1].set_title('After High-Pass Filter', fontsize=16)\nax[0, 2].plot(denoise_signal(signals[0]), 'darkmagenta')\nax[0, 2].set_title('After Wavelet Denoising', fontsize=16)\nax[0, 3].plot(denoise_signal(high_pass_filter(signals[0], low_cutoff=10000, SAMPLE_RATE=4000000), wavelet='haar', level=1), 'indigo')\nax[0, 3].set_title('After High-Pass Filter and Wavelet Denoising', fontsize=16)\n\nax[1, 0].plot(signals[1], 'crimson') \nax[1, 0].set_title('Original Signal', fontsize=16)\nax[1, 1].plot(high_pass_filter(signals[1], low_cutoff=10000, SAMPLE_RATE=4000000), 'mediumvioletred') \nax[1, 1].set_title('After High-Pass Filter', fontsize=16)\nax[1, 2].plot(denoise_signal(signals[1]), 'darkmagenta')\nax[1, 2].set_title('After Wavelet Denoising', fontsize=16)\nax[1, 3].plot(denoise_signal(high_pass_filter(signals[1], low_cutoff=10000, SAMPLE_RATE=4000000), wavelet='haar', level=1), 'indigo')\nax[1, 3].set_title('After High-Pass Filter and Wavelet Denoising', fontsize=16)\n\nax[2, 0].plot(signals[2], 'crimson') \nax[2, 0].set_title('Original Signal', fontsize=16)\nax[2, 1].plot(high_pass_filter(signals[2], low_cutoff=10000, SAMPLE_RATE=4000000), 'mediumvioletred') \nax[2, 1].set_title('After High-Pass Filter', fontsize=16)\nax[2, 2].plot(denoise_signal(signals[2]), 'darkmagenta')\nax[2, 2].set_title('After Wavelet Denoising', fontsize=16)\nax[2, 3].plot(denoise_signal(high_pass_filter(signals[2], low_cutoff=10000, SAMPLE_RATE=4000000), wavelet='haar', level=1), 'indigo')\nax[2, 3].set_title('After High-Pass Filter and Wavelet Denoising', fontsize=16)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Clearly, we can see that the high-pass filter and wavelet denoising techniques are able to effectively denoise the signal **by removing the unnecesary artificial impulse and additional noise** from the seismic signals. The high-pass filter alone does not perform as well as the wavelet denoising alone, but **the high-pass filter followed by wavelet denoising outperforms both methods**. The wavelet method works well because it implements a sort of deconvolution process that helps remove the artificial impulse. **The high-pass filter can only filter noise and the wavelet denoising method can only remove the artificial signal, but together they can remove both**.  "},{"metadata":{},"cell_type":"markdown","source":"**As you can see, we have achieved the same shape as the deconvolved signal in the image earlier in the kernel.**"},{"metadata":{},"cell_type":"markdown","source":"## Average Smoothing\nIn this denoising method, we take a \"window\" with a fixed size (like 10). We first place the window at the beginning of the time series (first ten elements) and calculate the mean of that section. We now move the window across the time series in the forward direction by a particular \"stride\", calculate the mean of the new window and repeat the process, until we reach the end of the time series. All the means we calculated are concatenated into a new time series, which is the denoised signal. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def average_smoothing(signal, kernel_size, stride):\n    sample = []\n    start = 0\n    end = kernel_size\n    while end <= len(signal):\n        start = start + stride\n        end = end + stride\n        sample.append(np.mean(signal[start:end]))\n    return np.array(sample)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualize the effect of average smoothing on seismic signals (with differing kernel sizes)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=3, ncols=4, figsize=(30, 15))\n\nax[0, 0].plot(signals[0], 'mediumaquamarine') \nax[0, 0].set_title('Original Signal', fontsize=14)\nax[0, 1].plot(average_smoothing(signals[0], kernel_size=5, stride=5), 'mediumseagreen')\nax[0, 1].set_title('After Average Smoothing (kernel_size=5 and stride=5)', fontsize=14)\nax[0, 2].plot(average_smoothing(signals[0], kernel_size=7, stride=5), 'seagreen')\nax[0, 2].set_title('After Average Smoothing (kernel_size=7 and stride=5)', fontsize=14)\nax[0, 3].plot(average_smoothing(signals[0], kernel_size=9, stride=5), 'darkgreen')\nax[0, 3].set_title('After Average Smoothing (kernel_size=9 and stride=5)', fontsize=14)\n\nax[1, 0].plot(signals[1], 'mediumaquamarine') \nax[1, 0].set_title('Original Signal', fontsize=14)\nax[1, 1].plot(average_smoothing(signals[1], kernel_size=5, stride=5), 'mediumseagreen') \nax[1, 1].set_title('After Average Smoothing (kernel_size=5 and stride=5)', fontsize=14)\nax[1, 2].plot(average_smoothing(signals[1], kernel_size=7, stride=5), 'seagreen')\nax[1, 2].set_title('After Average Smoothing (kernel_size=7 and stride=5)', fontsize=14)\nax[1, 3].plot(average_smoothing(signals[1], kernel_size=9, stride=5), 'darkgreen')\nax[1, 3].set_title('After Average Smoothing (kernel_size=9 and stride=5)', fontsize=14)\n\nax[2, 0].plot(signals[2], 'mediumaquamarine') \nax[2, 0].set_title('Original Signal', fontsize=14)\nax[2, 1].plot(average_smoothing(signals[2], kernel_size=5, stride=5), 'mediumseagreen') \nax[2, 1].set_title('After Average Smoothing (kernel_size=5 and stride=5)', fontsize=14)\nax[2, 2].plot(average_smoothing(signals[2], kernel_size=7, stride=5), 'seagreen')\nax[2, 2].set_title('After Average Smoothing (kernel_size=7 and stride=5)', fontsize=14)\nax[2, 3].plot(average_smoothing(signals[2], kernel_size=9, stride=5), 'darkgreen')\nax[2, 3].set_title('After Average Smoothing (kernel_size=9 and stride=5)', fontsize=14)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the average smoothing method is **only able to remove some additional noise, but not the artificial impulse as it is only able to smooth the curve**."},{"metadata":{},"cell_type":"markdown","source":"**The high amplitude, unstable parts of the signal (with high peaks and valleys) are the actual signal we are looking for. The medium amplitude patches of the signal (between the high amplitude regions) represent the unnecessary noise and artificial impulse, and the wavelet method seems to do better at removing these patches. This denoising illustrated in the image at the beginning of the kernel.**"},{"metadata":{},"cell_type":"markdown","source":"**We can therefore conclude that the wavelet denoising method is overall more effective than the average smoothing method.**"},{"metadata":{},"cell_type":"markdown","source":"**I believe that denoising the signals can significantly boost the scores of all the models (NN or LightGBM).**"},{"metadata":{},"cell_type":"markdown","source":"### That's it ! Thanks for reading my kernel ! Hope you found it useful :)"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":4}