# Convolution with interpolated function

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.

-
Oh, there is always a way to increase the calculation time :) Perhaps you meant decrease? – Ajasja Nov 12 '12 at 10:21
As a minimal improvement, why not calculate the integral with the delta once and for all ? – b.gatessucks Nov 12 '12 at 10:40
@J. M. By the way: Thanks a lot for this careful edit and thanks for the alternative link to my example data! – partial81 Nov 12 '12 at 12:02
Why do the convolution with interpolating functions? Since the data is discrete, why not convolve the data directly with (a sampled version of) the Gaussian? The simplest way to do this is to use ListConvolve. It will be very fast. – bill s Nov 12 '12 at 12:39
@bills GaussianFilter is even simpler since you don't have to convolve manually and you don't have to build the Gaussian-kernel. – halirutan Nov 12 '12 at 14:47

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[
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[

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
-
 Thanks for this post! This seems to work well and very fast. For me it is fascinating to see that for r=3.8 you get a result which is very similar to the output of my calculation. Do you know how your r=3.8 and my 0.008 are connected mathematically? – partial81 Nov 13 '12 at 9:08 @partial81 Yes, I can probably explain what r you have to use, at least theoretically. Can you tell why you unaccepted this answer without writing a comment for the reason? – halirutan Nov 23 '12 at 1:51 It would be fantastic to see a theoretical reason how r of the GaussianFilter is connected to the standard deviation of the normal distribution which is convolved with another function. That’s also why I unaccepted your answer – I thought I see a reason but this was not correct (unfortunately I forgot to mention this in a comment). So, if you could explain or provide information how these two values are connected, I would be very happy. But I have doubts that there is a connection because I found no one in the books I took a look into. – partial81 Nov 23 '12 at 9:52 @partial81 See my updated answer. – halirutan Nov 25 '12 at 4:46 WOW! What a detailed answer! I did not expect to get something back like this! But due to your explanation I think I understood why it works so well with the GaussianFilter. Thanks a lot for your work and support!! – partial81 Nov 25 '12 at 12:12

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.

-
 Thanks for this post! I like it because I know with your solution what I am doing. But I hope to understand the GaussianFilter of @halirutan post better soon because it is so much faster. – partial81 Nov 13 '12 at 9:11