Convolve working with Interpolated Functions in Mathematica
I'm using Mathematica 7.
I have an interpolated function, here's an example:
pressures =
WeatherData["Chicago", "Pressure", {2010, 8}] //
DeleteCases[#, {_, _Missing}] & //
Map[{AbsoluteTime[#[[1]]], #[[2]]} &, #] & // Interpolation;
I'd like to compute it's derivative, which is straight forward:
dpressures = D[pressures[x], x]
Now, If you plot this funciton
Plot[3600*dpressures, {x, AbsoluteTime[{2010, 8, 2}], AbsoluteTime[{2010, 8, 30}]}]
(sorry, don't know how to post the image from within 开发者_开发问答Mathematica, and don't have time to figure it out.) You'll find that it's very noisy. So, I'd like to smooth it out. My first thought was to use Convolve, and integrate it against a Gaussian kernel, something like the following:
a = Convolve[PDF[NormalDistribution[0, 5], x], 3600*dpressures, x, y]
Returns
360 Sqrt[2/\[Pi]] Convolve[E^(-(x^2/50)), InterpolatingFunction[{{3.48961266 10^9, 3.49228746 10^9}},<>], ][x], x, y]
Which looks reasonable to me. Unfortunately, I believe I have made a mistake somewhere, because the result I get back does not seem to be evaluatable. That is:
a /. y -> AbsoluteTime[{2010, 8, 2}]
Returns
360 Sqrt[2/\[Pi]] Convolve[E^(-(x^2/50)), InterpolatingFunction[{{3.48961266 10^9, 3.49228746 10^9}},<>][x], x, 3489696000]]
Which just ain't what I was looking for I'm expecting a number between -1 and 1.
Convolve seeks a closed form for the convolution. You could try a numerical convolution, starting with something like
NConvolve[f_, g_, x_, y_?NumericQ] :=
NIntegrate[f (g /. x -> y - x), {x, -Infinity, Infinity}]
However, for this noisy non-smooth function the numerical integration will struggle. (It won't work with default settings, and would be slow even with carefully chosen settings.)
I suggest you operate directly on the underlying data instead of interpolating the noisy data.
The bounds of your time range:
In[89]:= {lower = Min[First[pressures]], upper = Max[First[pressures]]}
Out[89]= {3.48961*10^9, 3.49229*10^9}
Use your interpolation to get samples every hour*:
data = Table[pressures[x], {x, lower, upper, 3600}];
Now compare
ListLinePlot[Differences[data]]
with smoothed version over 5 hour window:
ListLinePlot[GaussianFilter[Differences[data], 5]]
- You may want to use InterpolationOrder -> 1 for noisy data.
精彩评论