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).