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 UnitBoxe
s. 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 UnitBox
es 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
InterpolationOrder
s
- why is the
InterpolatingFunction
shifted by delta/2
when I use InterpolationOrder->0
and not the UnitBox
approach
GaussianFilter
is even simpler since you don't have to convolve manually and you don't have to build the Gaussian-kernel. $\endgroup$