Re: Derivative of experimental data

*To*: mathgroup at smc.vnet.net*Subject*: [mg124722] Re: Derivative of experimental data*From*: "Oleksandr Rasputinov" <oleksandr_rasputinov at ymail.com>*Date*: Thu, 2 Feb 2012 04:56:58 -0500 (EST)*Delivered-to*: l-mathgroup@mail-archive0.wolfram.com*References*: <jgauad$eq8$1@smc.vnet.net>

On Wed, 01 Feb 2012 08:49:17 -0000, Gabriel Landi <gtlandi at gmail.com> wrote: > Dear MathGroup users, > > I have a question which is very important for my current research, and > which involves not only Mathematica, but computer science in general. > > I have experimental data which is not very noisy, a small example of > which may be downloaded here. > > Basically I need to compute the derivative of this data. Here are my > options so far: > > Option 1: Finite differencing. Its terrible since the noise enhances > dramatically. > Option 2: Fitting some arbitrary function. The problem is that the > general functional form of the data changes from experiment to > experiment, so it is not possible to find a function which fits > adequately in all cases. > Option 3: Savitzky-Golay filters (self-implemented in Mathematica, based > on the discussion in Numerical Recipes, 3rd Ed.). It doesn't seem to > make much of a difference; probably because my data is not really that > noisy. > Option 4: Smoothing Splines filter. I am currently using Mr. Ludsteck > package HPFilter. So far it is by far the best outcome. However, it is > not free of some wild oscillations that are clearly non-analytical and > which are giving me quite the headache. > > Any suggestions are more than welcome. > I really appreciate any help I can get. > > Best regards, > > Gabriel Landi Ordinarily I would suggest the Savitzky-Golay method for this problem, but since you mention it is not really necessary for your data, perhaps you can use a straightforward interpolation? len = 25; data = Transpose[{Range[0., len]/len, RandomReal[{-1, 1}, len + 1]}]; ifun = Interpolation[ data, InterpolationOrder -> 3, Method -> "Spline" ]; Show[{ Plot[{ifun[x], 1/len ifun'[x]}, {x, 0, 1}], ListPlot[data, PlotStyle -> Red] }] (* looks OK *) If this doesn't work, I think the Savitzky-Golay filter is most probably the way to go, but you might need to adjust the filter window width for best results. In case it would be helpful, here is my implementation of the filter. Since I generally only need it for smoothing (i.e. without differentiation) I haven't spent any time figuring out the correct procedure for dealing with the end-points when taking derivatives, although it may work well enough as-is. ClearAll[SavitzkyGolayHatMatrix]; SavitzkyGolayHatMatrix[ windowHalfWidth_Integer, polynomialOrder_Integer, derivativeOrder_Integer, opts___ ] /; 0 <= derivativeOrder <= polynomialOrder <= 2 windowHalfWidth := Block[{ wprec = WorkingPrecision /. {opts}~Join~ Options[SavitzkyGolayHatMatrix, WorkingPrecision], vandermondeMatrix }, vandermondeMatrix = Table[ If[i == j == 0, 1, i^j], {i, -windowHalfWidth, +windowHalfWidth}, {j, 0, polynomialOrder} ]~N~wprec; (-1)^derivativeOrder derivativeOrder! Drop[vandermondeMatrix, None, -derivativeOrder].Drop[ PseudoInverse[vandermondeMatrix], derivativeOrder] ]; Options[SavitzkyGolayHatMatrix] = {WorkingPrecision -> Infinity}; Attributes[SavitzkyGolayHatMatrix] = {NHoldAll}; ClearAll[SavitzkyGolayFilter]; SavitzkyGolayFilter[ data_?VectorQ, windowHalfWidth_Integer, polynomialOrder_Integer, derivativeOrder_Integer, opts___ ] /; 0 <= derivativeOrder <= polynomialOrder <= 2 windowHalfWidth < Length[data] := Block[{ wprec = WorkingPrecision /. {opts}~Join~ Options[SavitzkyGolayFilter, WorkingPrecision], dprec = Precision[data], ndat, hatMatrix }, ndat = If[wprec < dprec, N[data, wprec], data]; hatMatrix = SavitzkyGolayHatMatrix[ windowHalfWidth, polynomialOrder, derivativeOrder, WorkingPrecision -> Min[wprec, dprec] ]; Flatten[{ ListConvolve[#, Take[ndat, 2 windowHalfWidth + 1]] & /@ Reverse@Take[hatMatrix, -windowHalfWidth], ListConvolve[hatMatrix[[windowHalfWidth + 1]], ndat], ListConvolve[#, Take[ndat, -(2 windowHalfWidth + 1)]] & /@ Reverse@Take[hatMatrix, windowHalfWidth] }] ]; Options[SavitzkyGolayFilter] = {WorkingPrecision -> Infinity}; Attributes[SavitzkyGolayFilter] = {NHoldRest}; Example: data = Tan@TriangleWave@Range[0.20, 5.80, 0.01]; (* signal *) data += RandomReal[{-0.05, 0.05}, Length[data]]; (* noise *) ListPlot[ SavitzkyGolayFilter[data, 8 (* window half-width *), 3 (* polynomial order *), 1 (* derivative order *) ], Joined -> True, PlotRange -> All ] (* looks OK *) For wide windows with high polynomial orders, calculating the filter coefficients by this method can be subject to numerical instability. If so, you can set WorkingPrecision to whatever value is necessary (up to and including Infinity, although this uses rational arithmetic and so is slow).