An Introduction to Smoothing

Smoothing is a process by which data points are averaged with their neighbours in a series, such as a time series, or image. This (usually) has the effect of blurring the sharp edges in the smoothed data. Smoothing is sometimes referred to as filtering, because smoothing has the effect of suppressing high frequency signal and enhancing low frequency signal. There are many different methods of smoothing, but here we discuss smoothing with a Gaussian kernel. We hope we will succeed in explaining this phrase in the explanation below.

This page is designed to be read in conjunction with the matlab commands that were used to make the figures. The code makes some of the ideas a little clearer, to those with some familiarity with matlab. The matlab commands are contained in the file smoothtalk.m

Some example data for smoothing

The numbers were generated with matlab, by creating 40 successive random numbers from a normal distribution. See the smoothtalk.m file for how to do this.

The Gaussian kernel

In the standard statistical way, we have defined the width of the Gaussian shape in terms of sigma. However, when the Gaussian is used for smoothing, it is usual to describe the width of the Gaussian with another related measure, the Full Width at Half Maximum (FWHM).

The FWHM is the width of the kernel, at half of the maximum of the height of the Gaussian. Thus, for the standard Gaussian above, the maximum height is ~0.4. The width of the kernel at 0.2 (on the Y axis) is the FWHM. As x = -1.175 and 1.175 when y = 0.2, the FWHM is in fact 2.35. The FWHM is related to sigma by the formula (matlab format):

FWHM = sigma * sqrt(8*log(2))

Smoothing with the kernel

With Gaussian smoothing, the function that is used is our Gaussian curve..

So, let us say that we are generating the new, smoothed value for the 14th value in our example data set. We are using a Gaussian with FWHM of 4 units on the x axis. To generate the Gaussian kernel average for this 14th data point, we first move the Gaussian shape to have its centre on 14 on the x axis. In order to make sure that we don't do an overall scaling of the values aftern smoothing, we then divide the values in the Gaussian curve by the total area under the curve, so that the values add up to one (see the smoothtalk.m file for how to do this in matlab):

We take the values of our resulting function (from the y axis), at each of the points in the data (the x axis). Thus we generate Gaussian function values for 13 12 11 etc, and 15 16 17 etc. This gives us a discrete Gaussian.

In fact the Gaussian values for 12 13 14 15 and 16 are:

    0.1174    0.1975    0.2349    0.1975    0.1174

and the data values for the same points are:

    1.0645    0.3893    0.3490   -0.6566   -0.1946

We then multiply the Gaussian values by the values of our data, and sum the results to get the new smoothed value for point 14. Thus, the new value for point 14 is ... + 0.1174*1.0645 + 0.1975*0.3893 + 0.2349*0.3490 + 0.1975*-0.6566 + 0.1174*-0.1946 + ...

We then store this new smoothed value for future use, and move on, to x = 15, and repeat the process, with the Gaussian kernel now centred over 15. If we do this for each point, we eventually get the smoothed version of our original data (see the smoothtalk.m file again for a simple but inefficient way of doing this):

Other kernels

This would have the effect of replacing each data point with a straight average of itself and the nieghbouring points.

Smoothing in 2D

and the same thing with discrete values for each pixel in the image.

We then proceed as before, multiplying the values of the kernel (as shown in the figure above) by the data in the image, to get the smoothed value for that point, and doing the same for every point on the image.

The procedure is the same for 3D data, except the kernel is rather more difficult to visualize, being something like a sphere with edges that fade out, as the cone fades out at the edges in the 2D case.

In fact, it turns out that we don't have to generate these 2D and 3D versions of the kernel for the computations, because we get the same result as we do by applying the full 2 or 3D kernel, if we simply apply a one dimensional smooth sequentially in the 2 or 3 dimenensions. Thus, for 2 dimensions, we could first smooth in the x direction, and then smooth the x-smoothed data, in the y direction.

Why smooth?

The next few images show the matched filter theorem in action. First we can generate a simulated signal in a one dimensional set of data, by creating a Gaussian with FWHM 8 pixels, centred over the 14th data point:

Next, we add some random noise to this signal:

We then smooth with a matching 8 pixel FWHM filter:

and recover our signal well from the noisy data.

Thus, we smooth with a filter that is of matched size to the activation we wish to detect. This is of particular relevance when comparing activation across subjects. Here, the anatomical variability between subjects will mean that the signal across subjects may be expected to be rather widely distributed over the cortical surface. In such a case it may be wiser to use a wide smoothing to detect this signal. In contrast, for a single subject experiment, where you want to detect (for example) a thalamic signal, which may be in the order of a few mm across, it would be wiser to use a very narrow smoothing, or even no smoothing.

Finding the signal for any smoothing level

Worsley KJ, Marret S, Neelin P, Evans AC (1996) Searching scale space for activation in PET images. Human Brain Mapping 4:74-90

Another promising method is to use wavelet transforms; see: Federico Turkheimer's wavelet introduction

Matthew Brett (FB) 19/8/99