--- title: "Smoothing Out Noisy Signals Using Discrete Mollifier Transforms – Part 2" url: "https://www.bekk.christmas/post/2022/14/smoothing-out-noisy-signals-using-discrete-mollifier-transforms-part-2" authors: ["Simon Foldvik"] published: "2022-12-14" updated: "2024-11-28T14:18:29Z" language: "en-US" type: "article" tags: ["Algoritmer"] keywords: ["smoothing", "noise", "reduction", "transform", "operator", "time series", "signal", "processing", "algorithm"] --- # Smoothing Out Noisy Signals Using Discrete Mollifier Transforms – Part 2 In the [previous article](https://www.bekk.christmas/post/2022/12/smoothing-out-noisy-signals-using-discrete-mollifier-transforms) we considered how noise reduction filters may be useful in practical applications. In this article we walk through the mathematical and technical aspects of implementing the noise reduction filters used there. We conclude with final remarks comparing our regularization operators to the broader class of moving averages. The reader is assumed to be conversant in elementary calculus. ## In short Discrete mollifier transforms smooth noisy sampled signals by approximating the continuous mollification integral with a trapezoidal-rule-based convolution. The article works through the math, shows how to implement the filter efficiently, explains how to handle edge artifacts with padding, and relates the method to moving averages through an asymptotic comparison. - **For:** Developers or technically minded readers interested in signal smoothing and numerical implementation - **Good to know beforehand:** Elementary calculus _AI-generated_ ### Mathematical definition In this section we build up to the mathematical definition of our discrete mollifier transforms. Considerations of implementation are found in the appropriate sections below. #### The continuous case Pulling a rabbit out of the hat, consider the _standard mollifier with parameter epsilon_: ![Definition of standard mollifier with parameter epsilon](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/3bb2264acc5e19110203752667d7b3c94b978158-1458x218.png) the constant A > 0 being chosen so that ![Standard mollifiers are normalized](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/90c0e7c3bfd60250676b65f23fabe0cab9890d24-646x166.png) To be particular, ![Determination of normalizing constant for the standard mollifiers](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/62cf8b1bff6c1871f8931d4eacc71c0602d04fd9-1576x200.png) The standard mollifiers are smooth (infinitely differentiable) functions vanishing outside the neighbourhood of the origin with radius epsilon. ![Plot of a sample of standard mollifiers](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/4b104b4a5f6f0d0f7383f07ea3096f269127c381-1009x646.png) _A plot of the standard mollifiers for various epsilons. Note how the vanishing supports and the normalization constraint on the integrals force an ever sharper peak near zero._ The insight is that [convolution](https://en.wikipedia.org/wiki/Convolution) against the standard mollifiers produces smooth functions, even if the original function undergoing the convolution is highly irregular. To be precise, introduce the integral operator sending a function to its convolution against the standard mollifiers: ![Integral operator mapping a function to its convolution with a mollifier](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/f325ec1335e7a7aea0a1d5a6f58a78dd81b967d5-2278x178.png) We call this the _mollification_ of u with parameter epsilon. The convolution is well-defined whenever u is locally integrable in the Lebesgue sense, in which case its mollification is in fact a smooth function. We have thus obtained the smoothing operator in the continuous case: ![The mollification operator takes locally integrable functions to smooth functions](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/2fd34750092635112d41142f18e54cf2b93bac2c-878x174.png) It is from the good convergence properties of the mollifications back to the original function, as epsilon tends to zero, that the procedure gains its attractiveness. We will not delve into further details on that matter here, preferring instead to turn our attention to the discretization of the mollification operators themselves. #### The discrete case To discretize the mollification operators, we approximate the integrals defining the mollification of a locally integrable function u by means of the [trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule). Other schemes for numerical quadrature may be applied, but the trapezoidal method reduces in the final analysis down to a particularly nice convolutional form suitable for efficient implementation on a computer, and obtains satisfactory error bounds. Construe now u as being sampled at n equidistant points, indexed from 0 to n-1, and take epsilon to be k≥1 times the mesh size. Abusing notation ever so slightly, the integral we wish to approximate at index j is thus ![Convolutional integral to approximate](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/d930659b3c03fef9f60974e9311ee1c991979680-1176x178.png) Sparing the reader the inevitable juggling of indicies, the discrete mollification of u at index j works out to be ![Discrete mollification of u at index i](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/1b250d7f664fd4e7b9345469ce419507932a8bdf-1068x228.png) at least for indicies j satisfying 0 ≤ j-k < j+k ≤ n-1, that is, for k ≤ j ≤ n-k-1. It is this form we next seek to implement. ### Implementation We consider the naive implementation of the discrete mollification in Python. It is a simple exercise to port this snippet into any sensible programming language. Readers concerned with efficiency may consult the following subsection and optimize as they see fit. Translating the discrete mollification operator above into Python directly produces: ```python import numpy as np _NORMALIZATION_CONSTANT = 2.2522836210436354 def _mollifier(t: float) -> float: """Standard mollifier corresponding to epsilon = 1.""" return _NORMALIZATION_CONSTANT*np.exp(1.0/(t*t - 1.0)) if abs(t) < 1.0 else 0.0 def mollify(u: np.ndarray, k: int) -> np.ndarray: """Returns the discrete mollification of the given signal with window size `k`.""" if not isinstance(k, int): raise TypeError(f"Window size must be an integer: {k=}") if k < 1: raise ValueError(f"Window size must be at least 1: {k=}") if len(u.shape) > 1: raise ValueError(f"Input signal must be one-dimensional: {u}") n = u.size Mu = np.copy(u) for j in range(k, n - k): Mu[j] = sum(_mollifier(i/k)*u[j - i] for i in range(1 - k, k))/k return Mu ``` The effectiveness of the discrete smoothing operator defined above may be witnessed in the following figure. ![Comparison of true, noisy, and reconstructed signals](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/de893fc10e75f7c2ecc87c7ddc3239accdfc254f-1410x534.png) _Comparing the true signal (left) to the sampled/noisy signal (middle) and the reconstructed signal (right). The reconstruction on the right was obtained by feeding the middle signal through the discrete mollifier implemented above with k=8._ #### Optimizations Considering the convolutional form of the discrete mollification operators, the computational efficiency may be improved by precomputing the weights for various k and performing the convolution in a more optimized fashion, say through scipy.signal.fftconvolve. We undertake no such endevaour here. #### Correcting the tails The attentive reader may have noticed that the first and last k sample points of the original signal are left unaltered by our implementation, producing noisy artifacts at the tails of our filtered signal. This is undesirable and easy to fix. Remembering that the discrete mollification is only defined for indicies j satisfying k ≤ j ≤ (n-1)-k, we may simply pad our signal at both ends so as to achieve these inequalities for the portion of the padded signal corresponding to the original signal. In code: ```python import numpy as np _NORMALIZATION_CONSTANT = 2.2522836210436354 def _mollifier(t: float) -> float: """Standard mollifier corresponding to epsilon = 1.""" return _NORMALIZATION_CONSTANT*np.exp(1.0/(t*t - 1.0)) if abs(t) < 1.0 else 0.0 def _mollify(u: np.ndarray, k: int) -> np.ndarray: """Returns the discrete mollification of the given signal with window size `k`.""" if not isinstance(k, int): raise TypeError(f"Window size must be an integer: {k=}") if k < 1: raise ValueError(f"Window size must be at least 1: {k=}") if len(u.shape) > 1: raise ValueError(f"Input signal must be one-dimensional: {u}") n = u.size Mu = np.copy(u) for j in range(k, n - k): Mu[j] = sum(_mollifier(i/k)*u[j - i] for i in range(1 - k, k))/k return Mu def mollify(u: np.ndarray, k: int) -> np.ndarray: """Mollification filter with corrected tails.""" n = u.size padded = np.zeros(n + 2*k) padded[:k] = u[0] padded[k:k + n] = u padded[k + n:] = u[-1] mollified = _mollify(padded, k) return mollified[k:n + k] ``` Here is the result: ![Comparison before and after correcting filter tails](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/56529208b46709f175159d94fea7b4795f100edd-1357x529.png) _Comparison of the results of the discrete mollifier before (middle) and after (right) correcting noisy tail artifacts. The raw signal is shown on the left._ ### Comparison with the class of moving averages We now conclude with a brief comparison of our discrete mollifiers with the broader class of moving averages. A moving average is a convolutional operator of the form ![General form of a moving average operator](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/d7123e0353a68208faf59faa4b890c3e35d58978-749x156.png) with the weights summing to unity: ![Moving average weights summing to unity](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/cb2ba25342984fdf6a7eaa0e2b80f33be5a7024d-471x156.png) Even though our discrete mollifiers are not strictly of the moving average type, they are so asymptotically. Indeed, one may show the estimate ![Asymptotic estimate on sum of weights for the discrete mollifiers](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/b52666e176a422a54a81f0572f5063068bd420e8-1015x228.png) for some constant M≥0 independent of k, proving ![Quadratic asymptotics for the sum of weights in the discrete mollifiers](https://cdn.sanity.io/images/ah2n1vfr/bekk-blogg-prod/e7845950cebbc4501f4d14d07d701e090e5c7390-1072x228.png) ### Summary In this article we have undertaken the construction, analysis, and implementation of a family of discrete smoothing operators. Furthermore, the convolutional form of these operators paves the way for efficient implementation in any programming language, making noise reduction filters broadly available. We also related the discrete mollifiers defined here to the broader class of moving averages, discovering their asymptotic membership to this class. ## Relevant resources recommended by the author - [Smoothing Out Noisy Signals Using Discrete Mollifier Transforms – Part 1](https://www.bekk.christmas/post/2022/12/smoothing-out-noisy-signals-using-discrete-mollifier-transforms): A motivational discussion of the noise reduction filters designed here, exemplified by their use on electrocardiographic data.