{"cells":[{"metadata":{"_uuid":"16ba27ea34a6c244349961c4e24aecab6ffcae26"},"cell_type":"markdown","source":"## Signal decomposition with Fast Fourier Transforms\n\nIn a previous notebook, [\"denoising algorithms\"](https://www.kaggle.com/residentmario/denoising-algorithms/), I discussed algorithms that can be used to compute a smoothed/simplified version of a time series dataset. The notebook focused in large part on filters: algorithms which use some kind of probabilistic or mathematical properties of the data to determine how the data ought to be smoothed.\n\nThe **Fast Fourier Transform** or FFT is another kind of signal analysis technique. It can also be used to simplify a signal, but it differs from most filters in that it is a decomposition technique, e.g. it can be used to transform a signal into a sum of other, simpler signals.\n\nThe F in FFT is for the **Fourier series**. A Fourier series is an infinite sum of sinusodial functions with coefficients which, taken as a sum, roughly equals an input function. In other words it's a way of decomposing a function $f(x)$ into a bunch of components $a_1(t) + a_2(t) + \\ldots$ in some domain $t$ that is dependent on $x$, which when taken as a whole can reconstruct $f(x)$.\n\nFourier series have deep mathematical significance, for instance they are related to the eigenvectors of the input matrix. However they are mathematically complex. For the purposes of application the relevant thing to know is that by computing a Fourier trainsform and then inverting it, you may collect the \"dominant\" signal on a given frequency band. For larger frequency bands this will equate to smoothing out the function; for small frequency bands this will equate to singling out the noise in the function.\n\nFor example:"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"%%time\n\nfrom numpy.fft import rfft, irfft, rfftfreq\nfrom scipy import fftpack\nimport pandas as pd\n\ntrain_meta_df = pd.read_csv(\"../input/metadata_train.csv\").set_index('signal_id')\ntrain = pd.read_parquet(\"../input/train.parquet\")\ny = train_meta_df.target","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1911fa1801883f758ad13b87f4203b3fd0858e57"},"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ndef low_pass(s, threshold=1e4):\n    fourier = rfft(s)\n    frequencies = rfftfreq(s.size, d=2e-2 / s.size)\n    fourier[frequencies > threshold] = 0\n    return irfft(fourier)\n\nlf_signal_1 = low_pass(train.iloc[:, 0])\nplt.plot(train.iloc[:, 0], color='lightgray')\nplt.plot(lf_signal_1, color='black')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6857d14e74550222716447a12bffed042c56aa2d"},"cell_type":"markdown","source":"Notice that this method:\n1. Computes the FFT of the time-series.\n2. Samples it along the given frequencies.\n3. Thresholds it, removing all low-frequency Fourier series from the result.\n\nWe can also go the other way and keep only the low-frequency Fourier series:"},{"metadata":{"trusted":true,"_uuid":"d28b7df9c398d16158633424ce6fab9cc30be228"},"cell_type":"code","source":"def high_pass(s, threshold=1e7):\n    fourier = rfft(s)\n    frequencies = rfftfreq(s.size, d=2e-2/s.size)\n    fourier[frequencies < threshold] = 0\n    return irfft(fourier)\n\nhf_signal_1 = high_pass(train.iloc[:,0], threshold=1e4)\n\nplt.plot(hf_signal_1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"59f7adb61c45c21e321c69f27517154e02d1f214"},"cell_type":"markdown","source":"If we sum these two signals we closely approximate the original dataset, to the degree that the difference doesn't appear on the plot:"},{"metadata":{"trusted":true,"_uuid":"5c08ce22719ef6eb2a161df014512ed6e1317920"},"cell_type":"code","source":"plt.plot(train.iloc[: 0], color='black')\nplt.plot(lf_signal_1 + hf_signal_1, color='lightgray')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"00d4d45da2f0604eeaf1f815710c813760a8ba12"},"cell_type":"markdown","source":"FFT does not compute the true infinite Fourier series (because that is infinite). Instead it computes a subsequence of coefficients for Fourier series spaces at certain frequency intervals. This is why the low-band filter isolates noise, as it includes only sinusoidals with low periodicity, whilst the high-band filter isolates the smoothed function, as it includes only sinusoidals with high periodicity. Here is a plot of our example frequency sample space:"},{"metadata":{"trusted":true,"_uuid":"0792b3d56eec8d84bbd772cd2fb05d7b068dec23"},"cell_type":"code","source":"plt.plot(rfftfreq(train.iloc[:, 0].size, d=2e-2 / train.iloc[:, 0].size))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d1d8d8648fa642354e7dcf25bf4c7804d33ade2c"},"cell_type":"code","source":"import numpy as np\n\ndef decompose_into_n_signals(srs, n):\n    fourier = rfft(srs)\n    frequencies = rfftfreq(srs.size, d=2e-2/srs.size)\n    out = []\n    for vals in np.array_split(frequencies, n):\n        ft_threshed = fourier.copy()\n        ft_threshed[(vals.min() > frequencies)] = 0\n        ft_threshed[(vals.max() < frequencies)] = 0        \n        out.append(irfft(ft_threshed))\n    return out\n\ndef plot_n_signals(sigs):\n    fig, axarr = plt.subplots(len(sigs), figsize=(12, 12))\n    for i, sig in enumerate(sigs):\n        plt.sca(axarr[i])\n        plt.plot(sig)\n    plt.gcf().suptitle(f\"Decomposition of signal into {len(sigs)} frequency bands\", fontsize=24)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"228a55988e9f6f7494eb3bc5800041bb0cb257d7"},"cell_type":"markdown","source":"If we decompose a complex signal into frequency bands (evenly spaced bands in this example, but you can choose whatever bands you'd like) we can see which frequencies dominate the information in the signal."},{"metadata":{"trusted":true,"_uuid":"a5f2b2dec9d833e5260ad5ec3cce09bdbfc91419"},"cell_type":"code","source":"plot_n_signals(decompose_into_n_signals(train.iloc[:,0], 5))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"04c720b933f2cac2271b025b489bb3265fa807ee"},"cell_type":"markdown","source":"In this case we see that there is a lot of information at the highest frequency band (indicating an underlying macro signal), then low signal at the intermediate band, then a lot of signal at low bands. This indicates that there is a strong macro signal and strong low-frequency signal, e.g. random noise, but less intermediate-strength structure."}],"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}