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}];

Almost sin(x)

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

Not cos(x)

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?

  • 1
    $\begingroup$ This also looks like a question that a Laguerre filter could help. Based on Laguerre polynomials they need fewer inputs than a moving average and employ a feedback technique. I'll look around and see if I have an old Mma implmentation. $\endgroup$ – Jagra Feb 11 '13 at 14:28
  • $\begingroup$ You might want to look at mathematica.stackexchange.com/a/10997/1089 $\endgroup$ – chris Feb 12 '13 at 20:13


Since Version 12, Mathematica now incorporates a range of (underrated IMHO) regularisation methods to Fit and FindFit.

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}]];

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]]

Mathematica graphics

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}]

Mathematica graphics

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.

| improve this answer | |
  • $\begingroup$ Very nice! Can you explain what you mean by "the relative weight of the penalty"? $\endgroup$ – 0xFE Feb 13 '13 at 0:53
  • 1
    $\begingroup$ In the expression LinearSolve[Transpose[mat].mat +10^(lambda) Transpose[reg].reg,Transpose[mat].(Last /@ pts)], 10^lambda weights the penalty (or regularisation) relative to the goodness of fit term (Transpose[mat].mat). $\endgroup$ – chris Feb 13 '13 at 8:33

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}]

derivative with moving average

| improve this answer | |
  • 5
    $\begingroup$ You can knock down almost all the noise using xbi=Interpolation[Transpose[{data[[All,1]],LowpassFilter[data[[All,2]],8,Length[data],TukeyWindow,SampleRate->Floor[Pi/0.05],Padding->"Periodic"]}]]. $\endgroup$ – Xerxes Feb 11 '13 at 6:33
  • $\begingroup$ Your comment works better than your answer - why not to put it in the answer with plots? ;) $\endgroup$ – Vitaliy Kaurov Feb 11 '13 at 6:38
  • 2
    $\begingroup$ Well, the comment version is a lot more fiddly; look at all the magic numbers I had to put in it. Also, bandpass filters are only natively in v9, which I assume not everybody has yet. $\endgroup$ – Xerxes Feb 11 '13 at 7:04
  • $\begingroup$ I'm on v7 with access to v8, so no LowpassFilter. My original idea was actually something similar to MovingAverage. Unfortunately the Mathematica implementation of MovingAverage clobbers a value (maybe more than one?) each time. I ported the Matlab version of moving average and Nested it around 10 times with a window of 5. I wasn't happy with this solution, which is why I asked the question. $\endgroup$ – 0xFE Feb 12 '13 at 1:51

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},
    InverseFourier[(Boole[Abs[# - Round[(ld + 1)/2]] > num] & /@ Range[ld]) Fourier[y]]
   ], {num, 10, 60, 10}

Mathematica graphics

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

Alternatively, ListConvolve could be used:

  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.

Mathematica graphics

| improve this answer | |
  • $\begingroup$ Clever idea! Unfortunately when I tried plotting the derivative of the interpolation function of the filtered data, it was a total mess. It oscillated like crazy. $\endgroup$ – 0xFE Feb 12 '13 at 1:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.