Fourier spectral smoothing method

When it comes to smoothing spectra, Savitzky–Golay (SG) is the method to immediately springs to mind. SG is implemented in most commercial chemometrics packages and works reasonably well in most circumstances. This unfortunately leads to it being implemented without much thinking, just as a fixed recipe. The aim of this post is to show an alternative to SG, the Fourier spectral smoothing method, that is a smoothing method based on the Fourier Transform. We’ll see one case where a Fourier smoothing method can give some advantage over SG.

This post is another instalment of our series on data pre-processing. Feel free to take a look at some of our other posts on the topic.

• Two scatter correction techniques for NIR spectroscopy
• NIR data correlograms with Seaborn
• Savitzky–Golay smoothing method

In the previous list, the post on the SG method is the only pre-requisite to understand this tutorial. We’ll build on some of the concepts (and code) of that post here, so feel free take a look at it before moving on.

In this tutorial we are going to briefly discuss the theory behind Fourier Transform (FT) smoothing. I’m going to ask you to bear with me on this one; I generally try not to get too much into the maths, but in this case I will make an exception. An easy discussion on the properties of the FT may be new for some of you, and it’s absolutely worth the time.

Secondly, we are going to understand how we can smooth data via FT, how we can differentiate data with FT, and obviously how to combine the two concepts to do derivative smoothing.

As in the SG post, the data used in this tutorial are taken from the paper Abundance and composition of kaolinite on Mars: Information from NIR spectra of rocks from acid-alteration environments, Riotinto, SE Spain by J. Cuadros and collaborators.

The dataset is freely available for download here.

Finally, and to avoid possible misunderstandings, note that this article is not about FT infrared spectroscopy, which is a totally different topic. Here we are going to talk about FT as a mathematical tools to perform some analysis on a spectrum.

Theory and main properties of the Fourier Transform

The Fourier Transform is arguably one of the mathematical topics that is available in countless shapes and form all over the internet. That should come at no surprise, given that FT use spans a great deal of technical applications.

Therefore, for the sake of this humble tutorial, I’m not going to repeat topics you can easily find elsewhere: my aim here is to present the main properties of the FT in the light of our specific application, that is smoothing spectral data.

In this spirit, I will bravely skip the definition of FT altogether. Let’s just use a symbolic representation of the FT that goes like this. Let’s consider a function f(x), which is our spectrum. We represent the FT of this function with the symbol F(q). That is we use the same letter in lower-case for the original function and upper-case for the Fourier-transformed function. Also, we associate the coordinate q to the corresponding coordinate in the original space, x. In symbols:

f(x) \rightarrow F(q).

It goes without saying then, that f(x) and F(q) live in different spaces (or domains), and that such spaces are related. This is the key idea behind using FT for spectral pre-processing. Smoothing and differentiation are operations that can be easily done on F(q) in a way that is uniquely related to f(x). Moreover, it turns out that these operations are more easily accomplished in the Fourier domain, that is more easily accomplished on F(q) rather than f(x).

Alright, here comes the second important theory bit of today’s tutorial: the Convolution Theorem. Smoothing is a form of convolution. As we (very briefly) discussed in the SG post, a rolling average of a spectrum is the convolution of that spectrum with a window function. If we write the convolution of two functions as f(x) \ast g(x), the convolution theorem states that

f(x) \ast g(x) \rightarrow F(q)G(q).

That is the FT of the convolution product of two spectra is the simple product of the Fourier transforms of these two spectra. We’ll see how to use a Gaussian function multiplier in Fourier domain to perform a smoothing of the original spectrum.

Now, for those of you who made this far, I promise this is going to be the last theory concept we are going to throw on the table. This one’s called Fourier derivative theorem. If we denote with f'(x) and f''(x) respectively the first and second derivative of our spectrum, the derivative theorem states that

f'(x) \rightarrow iqF(q),

f''(x) \rightarrow -q^{2}F(q).

In the previous equations i is the imaginary unit. The Fourier derivative theorem holds for any derivative order (not just first and second), however we’ll stick to those as they are the most relevant for us.

If you stare for a second at the equations above, you’ll quickly understand that a derivative operation becomes a simple multiplication in Fourier space. More specifically is a multiplication by a linear ramp for the first derivative case and a quadratic ramp for the second derivative. In any case, it can be implemented very easily with a few lines of code.

At this point it is clear how you can implement differentiation and smoothing in a single step. Just need to chain the derivative and the convolution theorem. Specifically, for the second derivative, we are looking at doing something like this in the Fourier space

-q^{2}F(q)G(q).

Fourier spectral smoothing method in Python

OK, I really, really, appreciate you made it through the theory part. It’s finally time to implement these concepts in Python. As usual, we begin with the imports and by loading the data

As we have noted in the previous post, the data is contained in an Excel spreadsheet. By loading the files without additional options, the first tab of the data file is returned. Again, for this tutorial we are going to use only the very first spectrum (labelled ‘C1’) of the data set.

We are going to implement the following workflow: 1) Define a smoothing filter in Fourier space; 2) Calculate the FT of the spectrum; 4) Multiply the two 5) Take the inverse FT of the product to obtain the smoothed spectrum.

The first thing to address is the smoothing filter. We are going to define a super-Gaussian (sometimes called generalised Gaussian):

G(q) = A \exp \left[ \left( - \frac{(x-x_{0})^{2}}{2 \sigma^{2}}\right)^{m}\right] .

The parameter A is a normalisation factor, which doesn’t have to be specified in the code. When m=1 the function above becomes a conventional Gaussian. For larger values of  m it tends to a flat-top function. Unlike a flat-top function however, the super-Gaussian is always continuous, without abrupt jumps (Fourier transforms don’t like sharp edges). For more information check the Wikipedia page on the Gaussian function.

Here’s the code to generate a super-Gaussian, using the appropriate Scipy function.

We have defined a simple Gaussian window with standard deviation of 40 pixels, then we have shifted this function by half the width of the spectrum, so that its maximum value will be at the first pixel. The shift is required to align the function with the fast Fourier Transform operation (FFT), which we are going to see next.

Note: to increase the amount of smoothing one has to decrease the value of \sigma .

The FT is numerically implemented through some version of the Fast Fourier Transform (FFT) algorithm.

Implementations of this and related algorithm is available in the Discrete Fourier Transform routine of Numpy.

The most important trick, any time you want to use the FFT on a generic function, is to make that function symmetric. Let me take a minute to explain this point, because it’s quite important,  but often overlooked.

The FFT operation implicitly assumes that the array we are going to transform is periodic. That is that it repeats itself an infinite number of times in both direction. Our initial spectrum looks like this:

Notice that we have deliberately plotted the spectrum as a function of the pixel number, not the wavelength. Now, if you imagine this function to be repeated (tiled) an infinite number of times in both directions, you can realise that there are going to be abrupt changes in the signal any time a new repetition begins. That’s because the initial and the final value of this array are not equal.

We can however concatenate a flipped version of the array itself to make it symmetric.

which produces the following new array

The initial and the final values of XX are equal, and there won’t be abrupt jumps when the function is tiled.

Perfect, we are now ready to perform the Fourier smoothing, keeping in mind that we have to redefine the win  array to match the new size of the XX  array:

We have done the FFT of the symmetric array XX in the first line. In the second line we have multiplied the result by the window function (Gaussian), inverse FFT, took the real_part of the result and finally selected only the first half of the resulting array. Taking the real part is required as the FFT generally produces complex arrays, while we know that we expect a real-valued array. Selecting only the first half of the result is required to restore the original array size, because we had doubled it when we produced a symmetric array.

OK, here’s a comparison of the Fourier smoothing with the original data and an equivalent SG smoothing over a limited wavelength interval.

Green and red curve track pretty close together. In this case Fourier smoothing and conventional SG smoothing appear to produce a similar result. In the next section we are going to see a case where FT may help.

Fourier derivative smoothing

Arguably, the best advantage of using a Fourier-based derivative smoothing is the fact that, if properly done, we avoid losing points at the extrema of the spectrum.

But let’s go step by step. We’ll work out the case for the second derivative here. To apply the formula above for the derivative smoothing, we need to define an array representing the Fourier coordinate q. Here’s the code

The array qq has the same dimension as XX and it is shifted to be symmetric about zero. To be consistent with the FFT definition in Python, we also need to multiply by a factor 2\pi.

What is left to do is to define a new window and apply the formula

Note, we have decreased the window width, as we need to increase the amount of smoothing when dealing with second derivatives.

Here’s the comparison of our Fourier derivative smoothing with a SG-based method. We selected the initial wavelength interval to highlight the fact that Fourier smoothing avoid loosing the information at the beginning of the array.

Take a look at how we can make use of the full spectrum, the contrast of the second derivative peaks is better and the curve appears generally smoother than the SG counterpart.

Wrapping up

The topic of this tutorial was spectral smoothing and derivatives using a Fourier-based method. We discussed some of the basic properties of the Fourier transform (FT) and how they can be employed to perform numerical differentiation and smoothing of spectra.

Fourier smoothing can be considered a good alternative to Savitzky–Golay (SG) smoothing, as it overcome some of the limitations inherently present in SG. The maths of the FT-based smoothing is arguably more involved than SG. Especially when dealing with derivatives however, Fourier methods have the advantage of not wasting any data point at the edges of the spectrum, as it happens when we do conventional numerical differentiation.

Hope you enjoyed the post. If you found it useful, please share it with a friend or colleagues and help us grow. Bye for now and thanks for reading!