{"cells":[{"metadata":{"trusted":true,"_uuid":"5226a95ba6e6ea70ec33bfe1334afdbfdb3d62e3"},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nimport pyarrow.parquet as pq\nimport gc\nimport pywt\nfrom statsmodels.robust import mad\nimport scipy\nfrom scipy import signal\nfrom scipy.signal import butter\n\nimport warnings\n\n# Suppress pandas future warnings, I am using different library versions locally\n# that do not raise warnings.\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n\ndata_dir = '../input'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ca160dc97e710c6ee002b2ff8a4b0c59d7079f6d"},"cell_type":"code","source":"print(scipy.__version__)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5cfda79d0648204193833b9a5aa1d4e8f4c90666"},"cell_type":"code","source":"metadata_train = pd.read_csv(data_dir + '/metadata_train.csv')\nmetadata_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4524f1ddb1414f8072a7f3d07189957ad6518532"},"cell_type":"code","source":"subset_train = pq.read_pandas(data_dir + '/train.parquet', columns=[str(i) for i in range(3)]).to_pandas()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a1f7dbcced9b5f139de0ffea754dba7fac08f896"},"cell_type":"code","source":"subset_train.memory_usage(index=True).sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1ac28211bd02087a9239e28ff13e6275a0fc30de"},"cell_type":"code","source":"# 800,000 data points taken over 20 ms\n# Grid operates at 50hz, 0.02 * 50 = 1, so 800k samples in 20 milliseconds will capture one complete cycle\nn_samples = 800000\n\n# Sample duration is 20 miliseconds\nsample_duration = 0.02\n\n# Sample rate is the number of samples in one second\n# Sample rate will be 40mhz\nsample_rate = n_samples * (1 / sample_duration)\n\ndef maddest(d, axis=None):\n    \"\"\"\n    Mean Absolute Deviation\n    \"\"\"\n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef 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    # scipy version 1.2.0\n    #sos = butter(10, low_freq, btype='hp', fs=sample_fs, output='sos')\n    \n    # scipy version 1.1.0\n    sos = butter(10, Wn=[norm_low_cutoff], btype='highpass', output='sos')\n    filtered_sig = signal.sosfilt(sos, x)\n\n    return filtered_sig\n\ndef 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    # Calculte 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":{"trusted":true,"_uuid":"a137e0bb2e22894ce620388954e24312224eb057"},"cell_type":"code","source":"train_length = 3\nfor i in range(train_length):\n    signal_id = str(i)\n    meta_row = metadata_train[metadata_train['signal_id'] == i]\n    measurement = str(meta_row['id_measurement'].values[0])\n    signal_id = str(meta_row['signal_id'].values[0])\n    phase = str(meta_row['phase'].values[0])\n    \n    subset_train_row = subset_train[signal_id]\n    \n    # Apply high pass filter with low cutoff of 10kHz, this will remove the low frequency 50Hz sinusoidal motion in the signal\n    x_hp = high_pass_filter(subset_train_row, low_cutoff=10000, sample_rate=sample_rate)\n    \n    # Apply denoising\n    x_dn = denoise_signal(x_hp, wavelet='haar', level=1)\n    \n    slice_size = 10000\n    font_size = 16\n    \n    fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(30, 10))\n    \n    ax[0, 0].plot(subset_train_row, alpha=0.5)\n    ax[0, 0].set_title(f\"m: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n    ax[0, 0].legend(['Original'], fontsize=font_size)\n    \n    # Show smaller slice of the signal to get a better idea of the effect the high pass frequency filter is having on the signal\n    ax[1, 0].plot(subset_train_row[:slice_size], alpha=0.5)\n    ax[1, 0].set_title(f\"m: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n    ax[1, 0].legend([f\"Original n: {slice_size}\"], fontsize=font_size)\n    \n    ax[0, 1].plot(x_hp, 'r', alpha=0.5)\n    ax[0, 1].set_title(f\"m: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n    ax[0, 1].legend(['HP filter'], fontsize=font_size)\n    ax[1, 1].plot(x_hp[:slice_size], 'r', alpha=0.5)\n    ax[1, 1].set_title(f\"m: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n    ax[1, 1].legend([f\"HP filter n: {slice_size}\"], fontsize=font_size)\n    \n    ax[0, 2].plot(x_dn, 'g', alpha=0.5)\n    ax[0, 2].set_title(f\"m: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n    ax[0, 2].legend(['HP filter and denoising'], fontsize=font_size)\n    ax[1, 2].plot(x_dn[:slice_size], 'g', alpha=0.5)\n    ax[1, 2].set_title(f\"m: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n    ax[1, 2].legend([f\"HP filter and denoising n: {slice_size}\"], fontsize=font_size)\n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"aec5f31af2a1dc55cf20452bd62bfcbd85484253"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}