{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.layers import * # Keras is the most friendly Neural Network library, this Kernel use a lot of layers classes\nfrom keras.models import Model\nfrom tqdm import tqdm # Processing time measurement\nfrom sklearn.model_selection import train_test_split \nfrom keras import backend as K # The backend give us access to tensorflow operations and allow us to create the Attention class\nfrom keras import optimizers # Allow us to access the Adam class to modify some parameters\nfrom sklearn.model_selection import GridSearchCV, StratifiedKFold # Used to use Kfold to train our model\nfrom keras.callbacks import * # This object helps the model to train in a smarter way, avoiding overfitting\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\n\nimport gc\nimport pywt\nfrom statsmodels.robust import mad\nimport scipy\nfrom scipy import signal\nfrom scipy.signal import butter\n\nimport warnings","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Reading the data**\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import pyarrow.parquet as pq\nos.listdir('../input/vsb-power-line-fault-detection')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The metadata_train.csv includes signal_id, id_measurement, phase, and target (4 columns) for 8712 (rows) data.\nThe metadata_test.csv includes signal_id, id_measurement, and phase (3 columns) for 20.3k (rows) data.\nThe sample_submission.csv includes signal_id and target which add \"the predicted target value\" to the metadata_test.csv\n\nThe data test.parquet & The data train.parquet contain signal data: 1 column = 1 signal ~ 800k data (for 20 ms = 1 single complete grid cycle)\nThe data test.parquet: 20.3k column entries\nThe data train.parquet: 8712 column entries ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Reading the entire parquet file is a one liner. Parquet will handle the parallelization and recover the original int8 datatype.\n#train = pq.read_pandas('../input/vsb-power-line-fault-detection/train.parquet').to_pandas() \nThis code line need more memory of Kaggle than allocated.\nShould read a subset of that data first","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#subset_train = pq.read_pandas('../input/vsb-power-line-fault-detection/train.parquet', columns=[str(i) for i in range(1)]).to_pandas()\nsubset_train = pq.read_pandas('../input/vsb-power-line-fault-detection/train.parquet', columns=[str(i) for i in range(3)]).to_pandas()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot 1 signal at three phase","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(subset_train)\nplt.ylabel('signal')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"DWT Signal Denoising","execution_count":null},{"metadata":{"trusted":true},"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":{},"cell_type":"markdown","source":"Plotting the denoising signals","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"data_dir = '../input/vsb-power-line-fault-detection'\nmetadata_train = pd.read_csv(data_dir + '/metadata_train.csv')\nmetadata_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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 #plotting only 10k (in 80k) data\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":{"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":4,"nbformat_minor":4}