Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
1  
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. –  Jagra Feb 11 '13 at 14:28
    
You might want to look at mathematica.stackexchange.com/a/10997/1089 –  chris Feb 12 '13 at 20:13
add comment

3 Answers

up vote 6 down vote accepted

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

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.

share|improve this answer
    
Very nice! Can you explain what you mean by "the relative weight of the penalty"? –  0xFE Feb 13 '13 at 0:53
1  
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). –  chris Feb 13 '13 at 8:33
add comment

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

share|improve this answer
5  
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"]}]]. –  Xerxes Feb 11 '13 at 6:33
    
Your comment works better than your answer - why not to put it in the answer with plots? ;) –  Vitaliy Kaurov Feb 11 '13 at 6:38
2  
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. –  Xerxes Feb 11 '13 at 7:04
    
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. –  0xFE Feb 12 '13 at 1:51
add comment

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

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:

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.

Mathematica graphics

share|improve this answer
    
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. –  0xFE Feb 12 '13 at 1:46
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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