{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"# Denoising algorithms\n\n## Background\nBasically any digital or analog process which produces a signal is vulnerable to noise. For example digital microphones get noise from random electron excitement (which occurs at any temperature above absolute 0), and analog film cameras get noise from differences in the sizes of the grains of the exposed film strip. When trying to quantify this variance there is additionally the [measurement error](https://en.wikipedia.org/wiki/Observational_error) of the instrument itself to content with; while the underlying .\n\nAn important preprocessing step for any signal processing task is denoising. There is a broad literature on denoising algorithms and techniques dating back centuries, and a number of fundamental techniques (like Fast Fourier Transforms) that have fallen out of from this.\n\nIn Python the principal library for performing signal processing tasks is the `signal` submodule in the `scipy` package. In this notebook I briefly go through the denoising algorithms available. This notebook is somewhat an elaboration on [xhlulu -- exploring signal processing with scipy](https://www.kaggle.com/xhlulu/exploring-signal-processing-with-scipy)."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy import signal","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"13cd7a6c62e183001053901750b9505763e4087c"},"cell_type":"code","source":"train = pd.read_parquet(\"../input/train.parquet\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"17aadfd43e848bf1d7688f76b6386fa710c554d2"},"cell_type":"markdown","source":"We will use the following demo data:"},{"metadata":{"trusted":true,"_uuid":"dc37ece19dfdad573878cbbf084393e94fa3b583"},"cell_type":"code","source":"train.iloc[:, 0:3].plot(figsize=(12, 8))\nplt.axis('off'); pass","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f1a4830b026558a3659f4970d7720e7871025029"},"cell_type":"code","source":"df = train.iloc[:, 0:3]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bd93e1cb1ea0d55c2039b5bf14e8df3d76d75dab"},"cell_type":"markdown","source":"## Rolling average (or other statistic)"},{"metadata":{"trusted":true,"_uuid":"ff9b16aeca603d4f64611e7b0fe60b4de1a985f3"},"cell_type":"code","source":"df.rolling(100).mean().plot(figsize=(12, 8))\nplt.axis('off'); pass","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f4c227f8c3542e75530602076b2faffd1922e86b"},"cell_type":"markdown","source":"The simplest algorithm for denoising time-series data is taking a summary statistic using a **rolling window**. A rolling window collects observations into groups of $n$ size. The groups are shifted one observation at a time, creating a \"window\" that passes over the dataset, hence the name. Each observation is part of $n-1$ groups, except for observations very near the beginning or the end, which appear in fewer (although you may choose to anneal them). Any summary statistic can be used within the rolling window to aggregate the data, though the average is the most popular default option.\n\nRolling windows have the following advantages:\n* They are computationally easy and simple to understand.\n* Varying $n$ can greatly increase or decrease the \"smoothing factor\".\n* They are flexible, you can use any algorithm you want as the aggregation function (mean, median, min, max, a quantile...).\n* They allow you to extend the algorithms you are already familiar with in record-oriented contexts to the time-series context without requiring any additional work.\n\nRolling windows have the following disadvantages:\n* They clip the beginnings and ends of the observations, thus reducing the total number of observations, which can represent a significant chunk of the dataset if the window is large and/or the data is small.\n* They cannot capture macro periodicity without removing micro periodicity.\n\nThe simplest rolling window technique is to apply an average on a neighborhood. However, you can apply any algorithm on the window, including a very complex one. `pandas` has very strong support for rolling window techniques built into the library (taken from `scipy`). You can define a more sophisticated window: e.g. an exponentially weighted window, which greatly reduces the contribution of points further off from the center of the neighborhood. For the most part more sophisticated windows mostly make the result more robust to outliers. For more details see [the docs](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.rolling.html).\n\nOne rolling window technique worth particular mention is the [median filter](https://en.wikipedia.org/wiki/Median_filter). When applied to an image a median filter has the effect of smoothing (color) noise in the image whilst retaining boundaries in the image. This is obviously useful in image processing because noisy images will have many spurious boundaries, bad, but smeared images will have weak boundaries, also bad."},{"metadata":{"trusted":true,"_uuid":"ccc94cb15093da1c0081cd0940432403d84cb6e9"},"cell_type":"code","source":"df.rolling(100, win_type='gaussian').mean(std=df.std().mean()).plot(figsize=(12, 8))\nplt.axis('off'); pass","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8c18b22fc2e17c188bda9d272b1a9a6a572ae42e"},"cell_type":"markdown","source":"## Convolution\n\nConvolution is a mathematical operation which generates a new function (or pointwise distribution) that is a function of two prior functions (or pointwise distributions). Convolution is geometrically a graph of the area under the curve of two functions which are moved towards and then away from one another (\"convolved\"). For example, if the two functions are both squares, their convolution will be a triangle function beginning where the moving squares just begin to touch, maxing out where they overlap completely, and going back to zero where they just touch the other way. In the continuous case it is defined mathematically as:\n\n$$(f * g)(t) = \\int_{-\\infty}^{\\infty}f(x)g(x-\\tau)d\\tau$$\n\n[The visualization](https://en.wikipedia.org/wiki/Convolution#/media/File:Convolution_of_box_signal_with_itself2.gif) from the Wikipedia article for it has a great intuitive explanation of how it works.\n\nAny operation that can be expressed using a rolling window can be expressed using a convolution. For example, we can re-express a rolling mean with a window size of 10 by asking for the convolution of the function with a square window function (or \"kernel\") of length 10, divided by 10 (since area / length = average height):"},{"metadata":{"trusted":true,"_uuid":"3ef189aae10a78f92a8173b0f88cc078452bd068"},"cell_type":"code","source":"# pd.Series(np.repeat([0., 1., 0.], 10)).plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eee7421d64ab40c59483768da026b0af30047e70"},"cell_type":"code","source":"def apply_convolution(sig, window):\n    conv = np.repeat([0., 1., 0.], window)\n    filtered = signal.convolve(sig, conv, mode='same') / window\n    return filtered\n\ndf.apply(lambda srs: apply_convolution(srs, 100)).plot(figsize=(12, 8))\nplt.axis('off'); pass","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"91da1766af28618eb005e6c15a907843dd6316e7"},"cell_type":"markdown","source":"If this equivalence doesn't jump out at you right away, a good reference on it is [this](https://matthew-brett.github.io/teaching/smoothing_intro.html) and [this](https://matthew-brett.github.io/teaching/smoothing_as_convolution.html). In general you shouldn't need to use convolution to perform smoothing; the rolling window methods are a better \"form factor\", but you may be performing convolution under the hood, as convolution is an efficient array-linearizable operation and generic rolling window functions are not."},{"metadata":{"_uuid":"b81ca7feae32f5136c32ae81671ad1677413f3b8"},"cell_type":"markdown","source":"## Filters\n\nA large number of algorithms have been developed and published over the years for denoising signals. These algorithms are generally known as **filters**, or sometimes \"digital filters\" (to distinguish from electrical analog filters). Denoising is one of the fundamental tasks of signal processing, and it seems that different fiends have developed slightly different approaches, dependent slightly on the particularities of the kinds of problems they tend to run into.\n\nFor a list of filters available in `scipy` see the corresponding entry in the [documentation](https://scipy.github.io/devdocs/signal.html).\n\nA lot of filters are based on rolling averages or convolutions under the hood because those are the core operations for operating on data. Good filters just make smart mathematical adaptations that, for some subset of tasks, outperforms simply taking the average of the window.\n\nThe next section describes a _relatively_ simple example of such a filter."},{"metadata":{"_uuid":"26dc54de35ea217ce6eb8af7f276a12c5ca79027"},"cell_type":"markdown","source":"## Savitzky–Golay filter\n\nB-splines or \"basis splines\" approximate data as a smooth sequence of piecewise exponential functions. The functions being _exponential_ means that they have terms of the form:\n\n$$c_1x^n + c_2x^{n-1} + \\ldots + c_n$$\n\nThe functions being _piecewise_ means that each one is only used for a particular subsequence of the data. For example $f(x)$ for $[0, 0.5]$ and $g(x)$ for $(0.5, 1]$. The points where functions intersect are known as _knot points_.\n\nThe last major requirement of a B-spline is that it must be smooth, e.g. continuous, e.g. it must have a derivative defined at every point. Geometrically this means that the area around the knots must be smooth, just like the spline is smooth on the body of the function. So the functions must join together in a way that preserves local topology.\n\nApproximating a dataset using B-splines means solving for the exponential function which best fits the data within a local area between knots. You can use any general-purpose optimization algorithm to do this; [stochastic gradient descent](https://www.kaggle.com/residentmario/gradient-descent-with-linear-regression/) for example (see the [\"Curve fitting\"](https://en.wikipedia.org/wiki/B-spline#Curve_fitting) section of the Wikipedia article for a bit more detail). More complex functions can theoretically fit data within a region better, but require more work to optimize, as they have more parameters to tune. In practice cubic curves are sufficient for almost all applications.\n\nHere is an example of a curve approximated using B-splines:\n\n![](https://i.imgur.com/gckO3IS.png)\n\nB-splines and their slightly more advanced cousin the Bezier curve play an important role in computer graphics.\n\nB-splines are more powerful and compact than rolling averages, as they generate a output function, not just a bunch of numbers."},{"metadata":{"trusted":true,"_uuid":"d520e791fc4156b9aa859285ab33d688d7d833ed"},"cell_type":"code","source":"# TODO: recipe based on https://scipy-cookbook.readthedocs.io/items/Interpolation.html\n# I had difficulty getting this working because I did not know how to interpret the result\n# of signal.cspline1d, which claims to generate coefficients, not values; but those \n# coefficients cannot be immediately mapped onto data values.\n\n# coef = df.apply(lambda srs: signal.cspline1d(srs[::50].values, 10))\n# df_approx = signal.cspline1d_eval(coef.iloc[:, 0].values.astype('float64'), df.index.values.astype('float64'))\n# coef.head(100).apply(lambda t: t[0]*t.name**3 + t[1]*t.name**2 + t[2], axis='columns').plot()\n# ^ goes to \\infty","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2166e840f1b5d7aa5ec6811311a38b54602e126e"},"cell_type":"markdown","source":"The problem with using B-spline approximation directly (besides the fact that it really isn't all that well documented; I don't quite understand the application surface) is that determining where the knots should be is difficult. A good durable procedure for doing B-spline approximation was discovered in the 1960s and applied to the field of analytical chemistry (where it is still very dominant). It's called a Savitzky-Golay filter; you can read about it in [its Wikipedia article](https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter).\n\nBasically all it does is roll a window along the data, build a curve in the left-right neighborhood of the central point (e.g. if the window size is 9, there are 4 points to the left and 4 to the right), and then set the value to that curve's approximation of the value at the middle of the neighborhood. A good animation showing this in action is located [here](https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter#/media/File:Lissage_sg3_anim.gif)."},{"metadata":{"trusted":true,"_uuid":"aeab8ae3949d67581193b8c1371475cb3e6dbd42"},"cell_type":"code","source":"df.apply(lambda srs: signal.savgol_filter(srs.values, 99, 3)).plot(figsize=(12, 8))\nplt.axis('off'); pass","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aec345ee78ffe1e00098e4c66c9808a7e827b62f"},"cell_type":"markdown","source":"This algorithm is that it has an analytical solution, which makes it quick to calculate. But it is starting to be a bit difficult to grasp mathematically."},{"metadata":{"_uuid":"93478388a0d37a00bec1d441fb39de7aef4a2e8c"},"cell_type":"markdown","source":"## Machine learning models\n\nComplex filters take the idea of reducing the noise in a dataset far out into the realm of statistical modeling. However almost all of the relevant techniques accessible from e.g. `scipy` and to a regular programmer are relatively ancient, having been well-developed by the 1970s or thereabout. They predate the rise and rise of machine learning over the last few years.\n\nIf you think about it, modeling data is not that much different from denoising it. No useful model captures all of the variance in a dataset, so it learns to find the explainable variance and omit the incidental variance. Building a machine learning model on a dataset is not that different from denoising it. Therefore we can use a machine learning algorithm as an alternative to a filter, e.g. instead of applying a filter to denoise a dataset we can apply a machine learning model. This is a lot like pretraining in neural networks! For example, here is a simple model curve using the K nearest neighbors algorithm:"},{"metadata":{"trusted":true,"_uuid":"98e9b5450a323310e6819e55f12ac7f92d1cbbfd"},"cell_type":"code","source":"from sklearn.neighbors import KNeighborsRegressor\n\nclf = KNeighborsRegressor(n_neighbors=100, weights='uniform')\nclf.fit(df.index.values[:, np.newaxis], \n        df.iloc[:, 0])\ny_pred = clf.predict(df.index.values[:, np.newaxis])\nax = pd.Series(df.iloc[:, 0]).plot(color='lightgray')\npd.Series(y_pred).plot(color='black', ax=ax, figsize=(12, 8))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a59aec53bfd05b38a96a900a9d28a93c05797e3b"},"cell_type":"markdown","source":"The `KNeighborsRegressor` algorithm is just a rolling window along the dataset that sets the value of the model at every point to be the (potentially smoothed) average of the surrounding points. If you think about it, that's not all that different from the Savitzky–Golay filter from earlier, which does almost the same thing; the only difference is that Savitzky-Golay builds a cubic spline and uses that instead.\n\nThe biggest disadvantage of using machine learning algorithm to do smoothing like this is that machine learning algorithms (particularly parametric ones like kNN regression) are generally much slower than filters. Filters were designed in an era where compute was many orders of magnitude more expensive. They tend to make statistical assumptions about the underlying data which, although never entirely valid, leads to simple computations and fast implementations. Typical machine learning algorithms can take much longer to do the same work, but don't make assumptions about the data. Additionally they tend to be more familiar to the practioner.\n\n\n## Addendum\n\nAs the notebook [VSB Power Line Faults EDA + Feature Engineering](https://www.kaggle.com/timothycwillard/vsb-power-line-faults-eda-feature-engineering) points out, likely the strongest signal in the VSB Power Line Fault Detection competition dataset is the level of noise in various segments of the dataset. Denoising algorithms are relevant to the dataset because you can back out the level of noise in the dataset by comparing the original dataset with the smoothed version. Try using that as a feature for this competition."}],"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}