Approximately Fit Data Without FindFit

Question

I have some noisy data, and I'd like to analyze the derivative of the function the data approximately matches. Interpolation is normally great, but doesn't work very well at all here.

data = Table[{x, Sin[x] + RandomReal[{-.15, .15}]}, {x, 0, 2 \[Pi], .05}];

xi = Interpolation[data];
Plot[xi'[t], {t, 0, 6.25}];

As you can see, the derivative (as expected when you think of how interpolation works) is not at all Cos[x]. In this case, because it is so obvious that the curve is of the Sin[x] family, it would be trivial to fit it and take the derivative of the fit. In the data I actually collected, I do not know the function, and assuming some high-order polynomial, while it may fit well, is almost certainly an incorrect assumption.

Is there a way to plot the derivative cleanly without requiring other knowledge about the function? Is it possible to tell Interpolation that it doesn't need to hit every point exactly?

Answer 1

As explained in this question, you can do a non-parametric fit to your data using B-Splines, and differentiate this fit:

pts = Table[{x, Sin[2 Pi x] + RandomReal[{-.15, .15}]}, {x, 0, 
    1, .0125}];

kfun[n_, d_] := 
  Join[ConstantArray[0, d], Range[0, 1, 1/(n - d)], 
   ConstantArray[1, d]]; 
uparam[pts_] := N[Range[0, 1, 1/(Length[pts] - 1)]];
mbasis[pts_, n_, d_] := 
  With[{param = uparam[pts]}, 
   Table[BSplineBasis[{d, kfun[n, d]}, j - 1, param[[i]]], {i, 
     Length[param]}, {j, n}]];

Clear[ctrlpts];
ctrlpts[lambda_: 0] := 
  With[{mat = mbasis[pts, 25, 3], 
    reg = SparseArray[{{i_, i_} -> 
        2., {i_, j_} /; Abs[i - j] == 1 -> -1.}, {25, 25}, 0.]}, 
   LinearSolve[Transpose[mat].mat + 10^(lambda) Transpose[reg].reg, 
    Transpose[mat].(Last /@ pts)]];

Show[ListPlot[pts, AxesLabel -> {x, y}], 
 ListLinePlot[{First /@ pts, mbasis[pts, 25, 3].ctrlpts[0.25]} // 
   Transpose, PlotStyle -> Red]]

Note that the fit does not go through all the points as you requested. Here we consider an explicit penalty function. The idea here is that we find the best (spline) weights subject to a prior corresponding to a roughness penalty (which allows us to tune how smooth the spine function should be, which involves adding a tunable cost to unsmooth spline).

The fit is now controlled by the relative weight of the penalty (given as an argument to ctrlpts). We can differentiate it:

df[x_] = BSplineFunction[ctrlpts[1], SplineDegree -> 3]'[x];
Plot[{df[x], 2 Pi Cos[2 Pi x]}, {x, 0.05, 0.95}]

There are known methods (such as cross validation) to estimate automatically what the proper amount of smoothing should be, depending on what it is you want to estimate (the function, its derivative, its second derivative etc.).

Note that the behaviour of your basis function at the edge of the requested interval needs to be addressed depending on what a proper boundary should be. For instance, the Fourier filtering method presented by others is formally equivalent to this B-Sline fit, while assuming periodic boundary condition and a particular choice of Wiener filter.

Answer 2

You want to remove high-frequency noise while retaining the low-frequency signal. This is a job for a bandpass filter. A simple one is the MovingAverage, which you can apply like so:

xsi = Interpolation[MovingAverage[data, 20]]
Plot[{Derivative[1][xsi][t], Cos[t]}, {t, 0, 6.25},  PlotRange -> {-1.1, 1.1}]

Answer 3

As Xerxes rightfully says LowPassFilter would be a good one if you have v9. A poor man's filter would be the following:

With[{x = data\[Transpose][[1]], y = data\[Transpose][[2]], ld = Length@data},
 Table[
  ListPlot[
   Chop[
    InverseFourier[(Boole[Abs[# - Round[(ld + 1)/2]] > num] & /@ Range[ld]) Fourier[y]]
    ]
   ], {num, 10, 60, 10}
  ]
 ]

It works by performing a DFT, multiplying the highest frequency components with 0, and doing an inverse DFT.

Alternatively, ListConvolve could be used:

Transpose[{
  data\[Transpose][[1]], 
  ListConvolve[{1, 1, 2, 1, 1}/6, data\[Transpose][[2]], {3, 3}]}
] // ListPlot

You could play with various kernels to see how it suits your data.