Convolve working with Interpolated Functions in Mathematica

Question

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.

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