Fitting piecewise splines to noisy data

Question

I have a system that gives me a noisy data set similar to the one generated by this matlab/octave code. The y-axis represents the signal intensity and the x-axis represents spatial distance.

noise = randn(1,512);
points = gaussian(512,.01);
noisy = points + .02*(noise');
plot(noisy)

The actual signal will be asymmetric and not the result of a gaussian function(but close). I can't create the noise accurately by hand and the underlying signal algorithmically so image these two together.

The goal is to construct an interpolating curve, using bicubic splines, based on only one data set of this signal. The signal cannot be captured multiple times and averaged together to remove the noise. Since the signal is noisy, I can't simply pick points every N samples as the interpolation points so I need to do some signal conditioning first.

I will not be able to use an exact model of the underlying physical process but the sources of noise are known and I can make some assumptions regarding the maximum slope and SNR. The interpolation is used to correct for manufacturing limitations so the curve is intended to be smooth(exclusive of noise) but for cost reasons is not. Therefore we don't expect any large deviations from this example.

My initial approaches used an M-point moving average filter to smooth the signal, then used polyfit handle the ends of the signal. This worked but could potentially reduce the center peak, if M is large enough. I've also since learned smoothing signals before interpolation is generally discouraged.

What would be a better approach to achieve this?

Answer 1

It seems like you don't want to make any prior assumptions about the exact shape of the signal, other than assuming that it is ``smooth''. This is best attacked using a non-parametric function fitting/regression approach. Use either a) spline fitting or b) kernel smoothing.

Two reasons why this is a good idea:

1. These methods are optimal in the sense that no other non-parametric filtering technique (moving average, etc.) can give you a better worst case error decay rate. This answers your main question: there is no better method if you care about minimax error.

2. Spline fitting and kernel regression codes are very easily available.

A quick Google search gives me the following: In Octave interp1 has a 'cubic' option: http://www.gnu.org/software/octave/doc/interpreter/One_002ddimensional-Interpolation.html

Open source code for kernel smoothing: http://code.google.com/p/pmtk3/source/browse/trunk/matlabTools/stats/kernelSmoothingRegression.m?r=2645&spec=svn2776

Answer 2

Not knowing the underlying model, you might just try a few things and see what best passes your fitness metric. Try interpolating the entire data set with polynomials of various degrees well below N. Try low-pass filters, maybe applied in the frequency domain to an FFT of your data set, with various cut-off frequencies and transition widths to roll-off any spectrum capable of producing a slope over the maximum expected.

You can calculate the maximum slope of each trial, and subtract the trial interpolation from the data to get a hypothetical noise signal. Iterate the polynomial degree and/or the filter passband until you get a maximim slope and a noise level that meets your fitness metric.

You might even be able to automate the above with either a successive approximation or genetic/evolutionary algorithm.