Optimizing calculation in for loop for convolution

Hello everyone,

I'm trying to compute the convolution of a heaviside and a gaussian function. For that purpose i use the following code:

x=-100:0.01:100;
sigma=5;
tzero=10;
Hea=heaviside(x-tzero);

to=min(x)+3*sigma:0.01:max(x)-3*sigma;


for i=1:length(to)
    sliding_gaussian=exp(-((x-to(i))/sigma).^2);
    S=@(x) Hea.*sliding_gaussian;
    q=integral(S,-100,100,'ArrayValued',true);
    int(i)=abs(sum(q)*0.01);
end

But it's taking ages. I need to include that into a fit so the complete computation will take even longer.

I tried the matlab function "conv" but i guess i don't know how to use it, since it didn't return what i expected - I appreciate any example that would fit my purpose.

Thank you for your help.

Answer 1

[Bruno Luong]

Use error function ERF

x = -100:0.01:100;
sigma = 5;
tzero = 10;
to = min(x)+3*sigma:0.01:max(x)-3*sigma;
convfun = @(t) sqrt(pi)*sigma/2*(1+erf((t-tzero)/sigma));
int = convfun(to);

plot(to,int)
xlim([-10,30])

Answer 2

[Bjorn Gustavsson]

This should work:

Hea = heaviside(x-tzero); % your shifted Heaviside
sliding_gaussian=exp(-(x/sigma).^2); % Gaussian
sliding_gaussian = sliding_gaussian/sum(sliding_gaussian); % Normalize it
GHc = conv(Hea,sliding_gaussian,'same'); % use convolution
figure
plot(x,[Hea;GHc])
% The dip at the end is due to implicit zero-padding of Hea in conv
% If you have the image processing toolbox you can use imfilter 
% which does pretty much the same convolution-job as conv, but with
% more options of edge-handling:
B = imfilter(Hea,sliding_gaussian,'replicate');
hold on
plot(x,B)

HTH

[Au.]

Hi Bjorn,

Thank you for your answer and the explanation on why there is a dip using conv.

This feature was actually the problem since i needed to fit a lifetime.