convolution.m

% compare different ways of convolution
% exp.decay (x) gauss
clear all;
h=0.05; %discretization step
x=-10:h:10;
y=zeros(1,length(x)); % non-convoluted signal
lambda=1.2;
for i=1:length(x)
    if x(i)>=0
        y(i)=lambda*exp(-x(i)*lambda); %int y dx = 1
    end
end
%============================================
% decay analytically convoluted with gauss,
% integral taken from Ryzhik, Gradshtein, p. 321 #3.322.2
sigma=0.25;
sigma2=sigma^2;
beta=sigma2/2;
gamma=lambda-x./sigma2;
analyt=lambda./sqrt(2*pi*sigma2).*sqrt(pi*beta).*exp(-x.^2./2./sigma2).*...
       exp(beta.*gamma.^2).*(1-erf(gamma.*sqrt(beta)));

plot(x,y,'r-',x,analyt,'b.');
hold on;

%============================================
% by-point convolution with internal matlab routine 'conv'
g=exp(-x.^2/(2*sigma2))./sqrt(2*pi*sigma2);
y1=conv(y,g);
x1=-20:h:20;
plot(x1,y1./length(x1)*40,'g+'); % normalization is L/Npoints
% to make this closer to analytical curve, decrease h
%============================================
% with fourier transform
f_y=fft(y);
f_g=fft(g);
f_conv=zeros(1,length(f_y));
f_conv=f_y.*f_g;
conv=abs(ifft(f_conv));
% NOTE: here first point is placed at x=h
% (not at x=0!)
ind=202; % x(202)==h
xnew=[x(ind:end),x(1:ind-1)];
plot(xnew,conv.*h,'r.');%,x(ind),conv(1)*h,'mx');
hold off;
legend('non-conv decay','analyt. conv. decay','by-point convolution','FFT conv decay');