{"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":"# Fast Fourier Transform & Denoising\ncomparacion entre tecnicas para limpiar ruido\n- Usar promedios para suavizar\n- Usar la FFT (Transformada Rapida de Fourier)\n\n## La FFT\nTengamos $x, y[x)$ como un punto de data, donde se pueden expresar los datos como una combinacion de componentes sinusoidales de frecuencia $\\omega_k = \\frac{2\\pi}{L}k$,\n\nSe ha de tener en cuenta que los datos son discretos, por tanto, tengamos $x_n, y_n$ como un punto de data\n$$y[x) \\rightarrow y_n$$\n\ncon el cambio a discretos\n$$ y_n = \\frac{1}{N} \\sum_{k=0}^{N-1} c_k \\exp(i \\frac{2\\pi x_n}{L}k)  $$\n$$ y_n = \\frac{1}{N} \\sum_{k=0}^{N-1} c_k \\exp(i \\omega_k x_n)  $$\n\ndonde:\n$$ x_n = \\frac{n}{N}L$$\nde esa manera, $x$ va tomando las distintas posiciones equiespaciadas con el indice $n$, asi dejamos de depender del intervalo $L$ y solo de la cantidad de datos $N$ y el indice $n$\n\n$$\ny_n = \\frac{1}{N} \\sum_{k=0}^{N-1} c_k \\exp(i \\frac{2\\pi k}{N} n) \n$$","metadata":{"_uuid":"63108ed7f840eba7d84f1fb5e6774a66316c1b4c"}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom numpy.fft import *\nimport pyarrow.parquet as pq\nimport matplotlib.pyplot as plt\n\nsns.set_style(\"whitegrid\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-11-23T15:29:55.780483Z","iopub.execute_input":"2022-11-23T15:29:55.780780Z","iopub.status.idle":"2022-11-23T15:29:57.149173Z","shell.execute_reply.started":"2022-11-23T15:29:55.780741Z","shell.execute_reply":"2022-11-23T15:29:57.147959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1 - Loading Data","metadata":{"_uuid":"2791aa2e427c195b5434f37ed425bc8be5775767"}},{"cell_type":"markdown","source":"### Signals","metadata":{"_uuid":"f9db7c1298259905c3a6276d2b09a76625049835"}},{"cell_type":"code","source":"# la data se guarda como pandas\n# son 999 signals distintas\n# N = 800000 muestras\nsignals = pq.read_table('../input/train.parquet', columns=[str(i) for i in range(999)]).to_pandas()\nprint(signals.shape)\n# data vs signal type\nsignals\n","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.status.busy":"2022-11-23T15:38:05.418079Z","iopub.execute_input":"2022-11-23T15:38:05.418445Z","iopub.status.idle":"2022-11-23T15:38:12.148062Z","shell.execute_reply.started":"2022-11-23T15:38:05.418374Z","shell.execute_reply":"2022-11-23T15:38:12.147054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# la traspuesta y un reshape, \n# en filas deja el dato de la señal\n# hemos de separar en 3, cada 3 datos de señal pertenecen a una medicion de señal\n    # pues se refiere a los 3 cables en la powerline, cada uno indicado con una fase\n# en columnas deja las observaciones\nsignals = np.array(signals).T.reshape((999//3, 3, 800000))\nprint(signals.shape)\n\nsignals","metadata":{"_uuid":"1255c51404246d23465d50575049db937fd7c6f8","execution":{"iopub.status.busy":"2022-11-23T15:39:03.924291Z","iopub.execute_input":"2022-11-23T15:39:03.924957Z","iopub.status.idle":"2022-11-23T15:39:06.270464Z","shell.execute_reply.started":"2022-11-23T15:39:03.924884Z","shell.execute_reply":"2022-11-23T15:39:06.269283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(signals[0, 2, :])","metadata":{"execution":{"iopub.status.busy":"2022-11-23T15:41:39.065425Z","iopub.execute_input":"2022-11-23T15:41:39.065739Z","iopub.status.idle":"2022-11-23T15:41:39.071288Z","shell.execute_reply.started":"2022-11-23T15:41:39.065693Z","shell.execute_reply":"2022-11-23T15:41:39.070331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 10))\nplt.title('Como se ve la data')\nplt.ylabel('Medicion')\nplt.xlabel('')\nplt.plot(signals[0, 0, :], label='Fase o Cable 0')\nplt.plot(signals[0, 1, :], label='Fase 1')\nplt.plot(signals[0, 2, :], label='Fase 2')\n\nplt.legend()\nplt.show()","metadata":{"_uuid":"832d8fb176ea537856fd6eae5bfb3425874a1fbd","execution":{"iopub.status.busy":"2022-11-23T15:32:43.435100Z","iopub.execute_input":"2022-11-23T15:32:43.435845Z","iopub.status.idle":"2022-11-23T15:32:48.291646Z","shell.execute_reply.started":"2022-11-23T15:32:43.435763Z","shell.execute_reply":"2022-11-23T15:32:48.290572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Target","metadata":{"_uuid":"6fb88dc514013b35e5915a5e520aab7ef07babaa"}},{"cell_type":"code","source":"# en target podemos observar donde se encuentran los fallos en las powerlines\ntrain_df = pd.read_csv('../input/metadata_train.csv')\ntrain_df","metadata":{"_uuid":"2221c7c5dd351386cbaa8162e622e535f64c3a1b","execution":{"iopub.status.busy":"2022-11-23T15:40:05.842361Z","iopub.execute_input":"2022-11-23T15:40:05.842653Z","iopub.status.idle":"2022-11-23T15:40:05.874741Z","shell.execute_reply.started":"2022-11-23T15:40:05.842608Z","shell.execute_reply":"2022-11-23T15:40:05.873919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# tomamos cada 3 datos, de manera que tendremos las mediciones\ntarget = train_df['target'][::3]\n","metadata":{"_uuid":"013f3ab108163f9272295192892ed44e7ec71e11","execution":{"iopub.status.busy":"2022-11-23T15:44:47.231185Z","iopub.execute_input":"2022-11-23T15:44:47.231829Z","iopub.status.idle":"2022-11-23T15:44:47.236708Z","shell.execute_reply.started":"2022-11-23T15:44:47.231761Z","shell.execute_reply":"2022-11-23T15:44:47.236058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 10))\nsns.countplot(target)\nplt.show()","metadata":{"_uuid":"226e535bbaf27f2f663c09fedddd555689ff11a8","execution":{"iopub.status.busy":"2022-11-23T15:44:10.512305Z","iopub.execute_input":"2022-11-23T15:44:10.512602Z","iopub.status.idle":"2022-11-23T15:44:10.876018Z","shell.execute_reply.started":"2022-11-23T15:44:10.512557Z","shell.execute_reply":"2022-11-23T15:44:10.874995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2 - Smoothing by mean\nThe idea is to reduce the length and the noise of the signal by merging $k$ neighbour values into their average.","metadata":{"_uuid":"b4f641105de1800b26081e6b0fc0dd93d912e452"}},{"cell_type":"code","source":"def sample(signal, kernel_size):\n    sampled = np.zeros((signal.shape[0], signal.shape[1], signal.shape[2]//kernel_size))\n    for i in range(signal.shape[2]//kernel_size):\n        begin = kernel_size * i\n        end = min(kernel_size * (i + 1), signal.shape[2])\n        sampled[:, :, i] = np.mean(signal[:, :, begin:end], axis=2)\n    return sampled","metadata":{"_uuid":"1f48858b8e21bba91768295c1382c1717601a6b2","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sampled = sample(signals, 100)","metadata":{"_uuid":"ca3b12709e3603dcc89625f0ff4adc948a9795e5","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 10))\nplt.plot(sampled[0, 0, :], label='Phase 0')\nplt.plot(sampled[0, 1, :], label='Phase 1')\nplt.plot(sampled[0, 2, :], label='Phase 2')\nplt.legend()\nplt.show()","metadata":{"_uuid":"728175d3acac3e7779e486f6df5a54139711212b","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3 - Fast Fourier Transform denoising\n\n#### A little bit of maths ...\nThe Fourier Transform of an 1D signal $x$ of length $n$ is the following : \n\n> ### $\\mathscr{f}_j = \\sum_{k=0}^{n-1} x_k e^{\\frac{2\\pi i}{n} jk} , ~~\\forall j=0, ... , n-1$ \n\nThe idea is to represent the signal in the complex space, It is roughly a sum of sinusoïdal functions. And there is one coefficient per frequency present in the signal.\n\nThe frequency takes the following values : \n- $f = \\frac{1}{dn} [0, 1, \\ldots ,   \\frac{n}{2}-1,  -\\frac{n}{2}, \\ldots , -1] $  if $n$ is even\n- $f =\\frac{1}{dn}  [0, 1, \\ldots,  \\frac{n-1}{2}, -\\frac{n-1}{2}, \\ldots, -1] $   if $n$ is odd\n\n#### Denoising algorithm\nThe denoising steps are the following :\n- Apply the fft to the signal\n- Compute the frequencies associated with each coefficient\n- Keep only the coefficients which have a low enough frequency (in absolute)\n- Compute the inverse fft\n","metadata":{"_uuid":"8594f0c680ec7ec26b214be18d0e5e5b72a19ae4","trusted":true}},{"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":{"_uuid":"bf7d1a88706a1dfebfe70675806ed799d529536f","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Testing some thresholds","metadata":{"_uuid":"481f1062cadb64df8776fc2b92fe8202c8e34d8a"}},{"cell_type":"code","source":"filtered = filter_signal(signals[0, 0, :], threshold=1e3)","metadata":{"_uuid":"a22aad3bdbf3369a2a180e4692ba29a99d4be078","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 10))\nplt.plot(signals[0, 0, :], label='Raw')\nplt.plot(filtered, label='Filtered')\nplt.legend()\nplt.title(\"FFT Denoising with threshold = 1e3\", size=15)\nplt.show()","metadata":{"_uuid":"ac58d384afa698f801f4a4f49d4550b48835dae3","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filtered = filter_signal(signals[0, 0, :], threshold=1e5)","metadata":{"_uuid":"92e910fc06f227234c97d61a10409dcb8df13025","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 10))\nplt.plot(signals[0, 0, :], label='Raw')\nplt.plot(filtered, label='Filtered')\nplt.legend()\nplt.title(\"FFT Denoising with threshold = 1e5\", size=15)\nplt.show()","metadata":{"_uuid":"1139e9e41ca5914d1c5857c1bd615c4a8c593349","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filtered = filter_signal(signals[0, 0, :], threshold=1e7)","metadata":{"_uuid":"0ba9d1babc5176696201591b8192e22a51409422","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 10))\nplt.plot(signals[0, 0, :], label='Raw')\nplt.plot(filtered, label='Filtered')\nplt.legend()\nplt.title(\"FFT Denoising with threshold = 1e7\", size=15)\nplt.show()","metadata":{"_uuid":"6de0fa39dbbead16f9207681867b984433b6275d","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Other uses of the fft ...\nThe fft coefficients can be used as features to represent the signal. I'll try that in a future kernel.\n\nHowever, there are too many of them so some more processing has to be done before feeding them into a classifier. Denoising being a solution.","metadata":{"_uuid":"6e0fc79ce57e97e66cdc036207bb34ea0dc9de82"}},{"cell_type":"markdown","source":"Thats all for now,\n#### *Thanks for reading !*\nAny feedback is appreciated","metadata":{"_uuid":"e745528af3801510dce81ea056e1d2137142d578"}}]}