{"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":"\nThere are many competitions on Kaggle that have data formats similar to this competition, which are high-frequency time-series data. Such as [LANL Earthquake Prediction](https://www.kaggle.com/competitions/LANL-Earthquake-Prediction), [VSB Power Line Fault Detection](https://www.kaggle.com/competitions/vsb-power-line-fault-detection)\n\nIn these competitions, some participants attempt to denoise the data and then use the processed features to improve their performance, and some of these solutions have indeed achieved better scores. Here, I will compile some denoising code. \n\nI think that in this competition, we can try to remove high-frequency noise data without processing the low-frequency data.\n\n","metadata":{}},{"cell_type":"markdown","source":"code reference:\n1. https://www.kaggle.com/code/tarunpaparaju/lanl-earthquake-prediction-signal-denoising\n2. https://www.kaggle.com/code/theoviel/fast-fourier-transform-denoising/notebook\n3. https://www.kaggle.com/code/ravishah1/parkinson-s-fog-accelerometer-data-eda","metadata":{}},{"cell_type":"code","source":"import os\n\nimport numpy as np\nfrom numpy.fft import *\n\nimport pandas as pd\nimport seaborn as sns\n\nimport pyarrow.parquet as pq\nimport matplotlib.pyplot as plt\n\njoin = os.path.join","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:11:24.705217Z","iopub.execute_input":"2023-05-25T09:11:24.706436Z","iopub.status.idle":"2023-05-25T09:11:24.714634Z","shell.execute_reply.started":"2023-05-25T09:11:24.706379Z","shell.execute_reply":"2023-05-25T09:11:24.713212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path_pj = '/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/'\n\npath_save = '/kaggle/working/'\n\npath_daori = path_pj \n\npath_dataori = path_daori\n\npath_train = join(path_daori, 'train')\npath_defog = join(path_train, 'defog')\n# path_notype = join(path_train, 'notype')\npath_tdcsfog = join(path_train, 'tdcsfog')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-05-25T09:11:24.717875Z","iopub.execute_input":"2023-05-25T09:11:24.718345Z","iopub.status.idle":"2023-05-25T09:11:24.730617Z","shell.execute_reply.started":"2023-05-25T09:11:24.718307Z","shell.execute_reply":"2023-05-25T09:11:24.729456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = pd.read_csv(join(path_tdcsfog, '003f117e14.csv'))","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:11:24.732169Z","iopub.execute_input":"2023-05-25T09:11:24.732668Z","iopub.status.idle":"2023-05-25T09:11:24.784877Z","shell.execute_reply.started":"2023-05-25T09:11:24.732628Z","shell.execute_reply":"2023-05-25T09:11:24.783511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## low-pass filter: simple FFT\n\nThis code implements a low-pass filter that takes an input signal, performs FFT on it, sets to zero the frequency components above a certain threshold, and then performs IFFT to obtain the output signal with high-frequency components removed. The function parameter `threshold` determines which high-frequency components are removed.","metadata":{}},{"cell_type":"code","source":"def filter_signal(signal, threshold=1e8):\n    fourier = rfft(signal)\n    frequencies = rfftfreq(signal.size, d=20e-3/signal.size)\n    fourier[frequencies > threshold] = 0\n    return irfft(fourier)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:11:24.786475Z","iopub.execute_input":"2023-05-25T09:11:24.787351Z","iopub.status.idle":"2023-05-25T09:11:24.794890Z","shell.execute_reply.started":"2023-05-25T09:11:24.787310Z","shell.execute_reply":"2023-05-25T09:11:24.794033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"signal = data['AccV']","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:11:24.797106Z","iopub.execute_input":"2023-05-25T09:11:24.798093Z","iopub.status.idle":"2023-05-25T09:11:24.814914Z","shell.execute_reply.started":"2023-05-25T09:11:24.798057Z","shell.execute_reply":"2023-05-25T09:11:24.813873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"threshold = 1e3\nfiltered = filter_signal(signal, threshold=threshold)\n\nplt.figure(figsize=(15, 10))\nplt.plot(signal, label='Raw')\nplt.plot(filtered, label='Filtered')\nplt.legend()\nplt.title(f\"FFT Denoising with threshold AccV = {threshold}\", size=15)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:21.837829Z","iopub.execute_input":"2023-05-25T09:12:21.838294Z","iopub.status.idle":"2023-05-25T09:12:22.316943Z","shell.execute_reply.started":"2023-05-25T09:12:21.838262Z","shell.execute_reply":"2023-05-25T09:12:22.314985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"threshold = 1e4\nfiltered = filter_signal(signal, threshold=threshold)\n\nplt.figure(figsize=(15, 10))\nplt.plot(signal, label='Raw')\nplt.plot(filtered, label='Filtered')\nplt.legend()\nplt.title(f\"FFT Denoising with threshold AccV = {threshold}\", size=15)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:34.146422Z","iopub.execute_input":"2023-05-25T09:12:34.146820Z","iopub.status.idle":"2023-05-25T09:12:34.612084Z","shell.execute_reply.started":"2023-05-25T09:12:34.146789Z","shell.execute_reply":"2023-05-25T09:12:34.610408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"threshold = 1e5\nfiltered = filter_signal(signal, threshold=threshold)\n\nplt.figure(figsize=(15, 10))\nplt.plot(signal, label='Raw')\nplt.plot(filtered, label='Filtered')\nplt.legend()\nplt.title(f\"FFT Denoising with threshold AccV = {threshold}\", size=15)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:37.036703Z","iopub.execute_input":"2023-05-25T09:12:37.037603Z","iopub.status.idle":"2023-05-25T09:12:37.424943Z","shell.execute_reply.started":"2023-05-25T09:12:37.037568Z","shell.execute_reply":"2023-05-25T09:12:37.423799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## low-pass&high-pass filter","metadata":{}},{"cell_type":"markdown","source":"The following code for plotting provides the original sequence, the sequence with low-frequency data removed, the sequence with high-frequency data removed, and the sequence with both low-frequency and high-frequency data removed. \n\nAs we can see from our data exploration, the time series data for this competition does not exhibit low-frequency trends or cyclicality. Therefore, I believe we only need to try to remove high-frequency noise data. Therefore, we should pay more attention the \"Original Signal\" and \"After Wavelet Denoising\" and params `low_cutoff` .","metadata":{}},{"cell_type":"code","source":"import pywt\n\nfrom scipy import signal\nfrom scipy.signal import butter, deconvolve","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:50.515959Z","iopub.execute_input":"2023-05-25T09:12:50.516485Z","iopub.status.idle":"2023-05-25T09:12:50.604873Z","shell.execute_reply.started":"2023-05-25T09:12:50.516447Z","shell.execute_reply":"2023-05-25T09:12:50.603490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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)\n\ndef high_pass_filter(x, low_cutoff=1000, SAMPLE_RATE=4000000):\n    \"\"\"\n    eliminate low freq data\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\n\n\ndef denoise_signal(x, wavelet='db4', level=3):\n    \"\"\"\n    \n    eliminate high freq data\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')","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:51.914301Z","iopub.execute_input":"2023-05-25T09:12:51.914702Z","iopub.status.idle":"2023-05-25T09:12:51.927855Z","shell.execute_reply.started":"2023-05-25T09:12:51.914674Z","shell.execute_reply":"2023-05-25T09:12:51.926686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def denoise_plot(signals, title='',start=None, stop=None, low_cutoff=10000, SAMPLE_RATE=150000, wavelet='haar', level=3):\n    fig, ax = plt.subplots(nrows=4, ncols=1, figsize=(30, 15))\n    \n    if title != '':\n        fig.suptitle(title, fontsize=20)\n    \n    ax[0].plot(signals, 'crimson') \n    ax[0].set_title('Original Signal', fontsize=16)\n\n    ax[1].plot(high_pass_filter(signals, low_cutoff=low_cutoff, SAMPLE_RATE=SAMPLE_RATE), 'mediumvioletred') \n    ax[1].set_title('After High-Pass Filter', fontsize=16)\n\n    ax[2].plot(denoise_signal(signals, wavelet=wavelet, level=level), 'darkmagenta')\n    ax[2].set_title('After Wavelet Denoising', fontsize=16)\n\n\n    ax[3].plot(denoise_signal(high_pass_filter(signals, low_cutoff=low_cutoff, SAMPLE_RATE=SAMPLE_RATE), wavelet=wavelet, level=level), 'indigo')\n    ax[3].set_title('After High-Pass Filter and Wavelet Denoising', fontsize=16)\n\n    for tk in range(4):\n\n        if start is not None and stop is not None:\n            if type(start) == int and type(stop) == int:\n                ax[tk].axvspan(start, stop, color='brown', alpha=0.2)\n            elif type(start) == list and type(stop) == list:\n                for i in range(len(start)):\n                    ax[tk].axvspan(start[i], stop[i], color='brown', alpha=0.2)\n            else:\n                raise NotImplementedError","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:54.033416Z","iopub.execute_input":"2023-05-25T09:12:54.035540Z","iopub.status.idle":"2023-05-25T09:12:54.054975Z","shell.execute_reply.started":"2023-05-25T09:12:54.035473Z","shell.execute_reply":"2023-05-25T09:12:54.053804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_subs_start_end(data):\n    start = None  \n    end = None  \n    \n    all_collect_start = []\n    all_collect_end = []\n    \n    for i in range(len(data) - 1):\n        if data[i+1] - data[i] != 1:  \n            if start is not None:  # \n                end = data[i]\n                print(f'Subsequence: {start}-{end}')\n                all_collect_start.append(start)\n                all_collect_end.append(end)\n                start = end = None  # \n        else:\n            if start is None:  # \n                start = data[i]\n\n    if start is not None and end is None:\n        end = data[-1]\n        print(f'Subsequence: {start}-{end}')\n        all_collect_start.append(start)\n        all_collect_end.append(end)\n    \n    return all_collect_start, all_collect_end","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:56.954590Z","iopub.execute_input":"2023-05-25T09:12:56.955166Z","iopub.status.idle":"2023-05-25T09:12:56.965744Z","shell.execute_reply.started":"2023-05-25T09:12:56.955099Z","shell.execute_reply":"2023-05-25T09:12:56.964108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = pd.read_csv(join(path_tdcsfog, '011322847a.csv'))","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:12:59.155475Z","iopub.execute_input":"2023-05-25T09:12:59.155889Z","iopub.status.idle":"2023-05-25T09:12:59.185772Z","shell.execute_reply.started":"2023-05-25T09:12:59.155859Z","shell.execute_reply":"2023-05-25T09:12:59.184486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_collect_start, all_collect_end = get_subs_start_end(data.loc[data['Turn'] == 1, 'Time'].values.tolist())","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:13:00.558992Z","iopub.execute_input":"2023-05-25T09:13:00.559419Z","iopub.status.idle":"2023-05-25T09:13:00.567981Z","shell.execute_reply.started":"2023-05-25T09:13:00.559387Z","shell.execute_reply":"2023-05-25T09:13:00.566743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"denoise_plot(data['AccV']\n             , title='AccV'\n             , start=all_collect_start\n             , stop=all_collect_end\n             , low_cutoff=5\n             , SAMPLE_RATE=128\n             , wavelet='Db3'\n             , level=2)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:13:22.883736Z","iopub.execute_input":"2023-05-25T09:13:22.884127Z","iopub.status.idle":"2023-05-25T09:13:24.607851Z","shell.execute_reply.started":"2023-05-25T09:13:22.884097Z","shell.execute_reply":"2023-05-25T09:13:24.606031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"denoise_plot(data['AccML']\n             , title='AccML'\n             , start=all_collect_start\n             , stop=all_collect_end\n             , low_cutoff=5\n             , SAMPLE_RATE=128\n             , wavelet='Db3'\n             , level=2)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:13:31.948253Z","iopub.execute_input":"2023-05-25T09:13:31.948656Z","iopub.status.idle":"2023-05-25T09:13:33.363317Z","shell.execute_reply.started":"2023-05-25T09:13:31.948626Z","shell.execute_reply":"2023-05-25T09:13:33.361942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"denoise_plot(data['AccAP']\n             , title='AccAP'\n             , start=all_collect_start\n             , stop=all_collect_end\n             , low_cutoff=5\n             , SAMPLE_RATE=128\n             , wavelet='Db3'\n             , level=2)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T09:13:33.855751Z","iopub.execute_input":"2023-05-25T09:13:33.856148Z","iopub.status.idle":"2023-05-25T09:13:35.191063Z","shell.execute_reply.started":"2023-05-25T09:13:33.856103Z","shell.execute_reply":"2023-05-25T09:13:35.189534Z"},"trusted":true},"execution_count":null,"outputs":[]}]}