Convolution with interpolated function

Question

I have some problems to calculate in reasonable speed the convolution of an interpolated function with a Gauss function.

I have here (ExampleData.txt, alternate Pastebin link) data which I interpolate linearly:

data = Import["ExampleData.txt", "Table"];
interPolFunc = Interpolation[data, InterpolationOrder -> 1]

First, I tried to work with following function (I have to add a Dirac delta function to my interpolated function):

firstTry[y_] := Integrate[
   PDF[NormalDistribution[0, 0.008],x-y]*(interPolFunc[x] + 4 DiracDelta[1-x]), {x, 0, 1.1}];

If I calculate some values with

firstTryTable = ParallelTable[{u, firstTry[u]}, {u, 0, 1.1, 0.005}];

I only get back the contribution of the Dirac delta function.

If I use this function

 secondTry[y_] :=
      Integrate[PDF[NormalDistribution[0, 0.008], x - y]*4*DiracDelta[1 - x],{x, 0, 1.1}]  +
      Integrate[PDF[NormalDistribution[0, 0.008], x - y]*interPolFunc[x],{x, 0, 1.1}];

and calculate values with

 secondTryTable = ParallelTable[{u, N[secondTry[u]]}, {u, 0, 1.1, 0.005}];

then the result seems to be reasonable.

So, my first question (not so important) is:
Why do I get back two different results?

The second is:
Is there any way to decrease the calculation time? It is really slow. I would be happy about some help.

Answer 1

I hope I see the essence here. You are interested in

the convolution of an interpolated function with a Gauss function

Your underlying data has regular spacings in x-direction and the convolution with a Gaussian is extremely fast implemented in GaussianFilter for discrete data. Why are you making it so complicated when the only thing you have to do is to filter your y-values?

data = Import[
   "https://github.com/downloads/stackmma/Attachments/ExampleData_14436.txt", "Table"];
ApplyGaussianFilter[data_, r_]:=Transpose[{#1, GaussianFilter[#2, r]}] & @@ Transpose[data];

ListLinePlot[ApplyGaussianFilter[data, 3]]

With r=3 it looks like the output of your calculation. If you want to adjust the setting for the GaussianFilter, check out its help-page.

And the timing is

First@AbsoluteTiming[ApplyGaussianFilter[data, 3]]

(* 0.001495 *)

Analytic vs. discrete approach

Let's say we have a short discrete signal. We pad it with enough zeroes to prevent boundary artefacts. Furthermore, for simplicity I will give no x-values for the signal. Instead, I assume that the first value is at x=1 and we have unit spacing between the values.

s = ArrayPad[RandomInteger[{0, 100}, 10], 25]
(*
{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 79,74,4,4,57,62,96,23,41,15,0,0,0,0,0,0,0,0,0,0,0,
 0,0,0,0,0,0,0,0,0,0,0,0,0,0}
*)

Analytic convolution

In your approach, you used an InterpolatingFunction which I won't do here. I will create a function from s by modulating shifted UnitBoxes. This is the same as when you use InterpolationOrder->0.

Assume you want to create a function wich is k in the interval between [5,6]. You could shift a UnitBox to 5.5 and multiply it by k to get the correct height. Now you could sum many of those modulated and shifted boxes and you get a different value for each interval. That's the trick we use to convert our signal s into an analytic function

analyticS=Sum[s[[i]]UnitBox[x-i+1/2],{i,Length[s]}]
(*15*UnitBox[-(69/2) + x] + 41*UnitBox[-(67/2) + x] + 
  23*UnitBox[-(65/2) + x] + 96*UnitBox[-(63/2) + x] + 
  62*UnitBox[-(61/2) + x] + 57*UnitBox[-(59/2) + x] + 
  4*UnitBox[-(57/2) + x] + 4*UnitBox[-(55/2) + x] + 
  74*UnitBox[-(53/2) + x] + 79*UnitBox[-(51/2) + x]
*)

You could now look at the plot of this and a e.g. ListInterpolation of s

ip = ListInterpolation[s, InterpolationOrder -> 0];
Plot[{analyticS, ip[x]}, {x, 1, Length[s]}, 
 PlotStyle -> {Automatic, Directive[Red, Dashed]}]

The good thing with analyticS is, that we can calculate the convolution analytically

analyticConv=With[{kernel=PDF[NormalDistribution[0,2],x]},
  Convolve[kernel,analyticS,x,t]
]
Plot[analyticConv, {t, 1, Length[s]}]
(*
  (1/(2*Sqrt[2*Pi]))*(15*Sqrt[2*Pi]*(-Erf[(-35 + t)/(2*Sqrt[2])] + 
  Erf[(-34 + t)/(2*Sqrt[2])]) + 41*Sqrt[2*Pi]*.....
*)

To compare this result with the discrete GaussianFilter we will discretizize this result. Note, that I sample it in the middle of each interval at 0.5, 1.5, ...

discrResult = Table[N@analyticConv, {t, 1/2, Length[s] - 1/2}]

Discrete convolution

What is the connection of this result with calling GaussianFilter directly on s is the question in the room. Several things have to be considered:

1.) When you call GaussianFilter[r] then a convolution matrix is created of radius r and the convolution matrix contains a discretized version of PDF[NormalDistribution[0, r/2], x]. Please note, that although the Gaussian decreases almost to zero inside the radius r of the convolution matrix, it does not vanish completely and the discretization will introduce further errors.

2.) If you want to see the discretization, which is used by GaussianFilter you can apply it to a unit impulse (which is a 1 padded by zeroes). Don't get confused by all the options I'm giving. It's only to show you the connection to a normal distribution:

GaussianFilter[ArrayPad[{1}, 2], 2, Method -> "Gaussian", "Standardization" -> False]
(* {0.053991, 0.241971, 0.398942, 0.241971, 0.053991} *)

You get the same result by discretizing a NormalDistribution with a sigma of r/2

Table[N@PDF[NormalDistribution[0, 1], x], {x, -2, 2}]
(* {0.053991, 0.241971, 0.398942, 0.241971, 0.053991} *)

This shows you 2 things: First, the GaussianFilter assumes that the distance between two values in your signal s is 1. If you, like in your case have another distance, you have to adjust sigma to get expected results. Second, you see, that the values of the filter are still at 0.053 at the boundary elements. To make the GaussianFilter more exact, you could use a larger radius r for a given sigma. You could use for instance GaussianFilter[s, {5, 1}] to get a filter with radius 5 and sigma of 1.

Let's apply such a GaussianFilter to s and compare the results. I'm using a large radius, to minimize truncation errors.

convDiscrete = GaussianFilter[s, {20, 2}, 
   Method -> "Gaussian", "Standardization" -> False]

Comparing the results

ListLinePlot[{convDiscrete, discrResult}, PlotRange -> All, 
 PlotStyle -> {Automatic, Directive[Red, Dashed]}]

This looks quite promising. Let's calculate the optimal sigma to see whether it can be improved

target[sigma_?NumericQ] := #.# &[(discrResult - GaussianFilter[s, {20, sigma},
  Method -> "Gaussian", "Standardization" -> False])];
NMinimize[{target[sigma], sigma > 0}, sigma]

(* {3.96549*10^-6, {sigma -> 2.02084}} *)

Visually, this makes no difference but it seems adding 2/100 makes the result more perfect.

Conclusion

What does that mean in your specific example? Let's first calculate the analytic convolution with sigma=0.008. delta is the space between each value and please note, that we now have to make the UnitBoxes smaller in width because of the small spacing delta.

data = Import[
   "https://github.com/downloads/stackmma/Attachments/ExampleData_14436.txt", "Table"];

delta = Subtract @@ data[[{2, 1}, 1]];
analyticData = Sum[elm[[2]]*UnitBox[(x - elm[[1]])/delta], {elm, data}];
convData = With[{kernel = PDF[NormalDistribution[0, 0.008], x]},
  Convolve[kernel, analyticData, x, t]
];

Let's take first a look whether the analyticData function represents your data correctly. The first entry which is not zero in your data is the 11th element. Lets make a Plot and draw a dashed line where data[[11]] is and check whether the column has a width of delta

Plot[{analyticData}, {x, 0, 0.1}, PlotRange -> All, 
 Epilog :> {Red, Dashed, Line[{#*{1, 0}, #*{1, 1.5}} &[data[[11]]]],
   Blue, Thickness[0.01], Dashing[0], 
   Line[{{# - delta/2, 0.8}, {# + delta/2, 0.8}} &[data[[11, 1]]]]}]

Now you have to calculate the appropriate sigma for the GaussianKernel. Remember, that we have to take care of the spacing which is not 1. You can calculate the parameter for the GaussianFilter using the sigma you used in the analytic convolution

$$ \sigma_g = \sigma/\Delta $$

this gives

0.008/delta

(* 1.68 *)

The rest of the approach is the same as in my first code-block at the very beginning of this answers. I just added the options to ensure that a real Gaussian is used:

Show[
 ListLinePlot@
    Transpose[{#1, 
      GaussianFilter[#2, {20, 1.68}, Method -> "Gaussian", 
       "Standardization" -> False]}] & @@ Transpose[data],
 Plot[{convData}, {t, 0, 1.1}, PlotStyle -> Directive[Red, Dashed]]
 ]

I hope this addition clears the situation so that you can now investigate on your own in

• what's the influence of different InterpolationOrders
• why is the InterpolatingFunction shifted by delta/2 when I use InterpolationOrder->0 and not the UnitBox approach

Answer 2

The ideas mentioned in comments and the prior response seem like good ways to go about this. As for the brute force direct method, for a reliable result you can precompute one part symbolically and handle the rest numerically.

ii[y_] = Integrate[
   PDF[NormalDistribution[0, 8/1000], x - y]*4 *DiracDelta[1 - x], {x,
     0, 11/10}];

firstTry[y_?NumberQ] := 
 ii[y] + NIntegrate[
   PDF[NormalDistribution[0, 0.008], x - y]*interPolFunc[x], {x, 0, 
    1.1}]
Timing[
 firstTryTable = Table[{u, firstTry[u]}, {u, 0, 1.1, 0.005}];]

(* Out[130]= {121.530000, Null} *)

ListLinePlot[firstTryTable]

will then give a plot similar to that shown by halirutan in the previous response.