Below you can see the code for convolution of two continuous functions. There is a function called fx which I took as the square root of a Gaussian distribution. The convolution is calculated using 2 methods. In one of them I use the built-in function conv() and in the other I use the definition of the convolution. In mat_conv1 and cont_conv1 the function fx is convolved with itself using the built-in and explicit calculation methods respectively. In mat_conv2 and cont_conv2 the two functions mat_conv1 and cont_conv1 are convolved with themselves using their respective methods. You can see the code below.
clear all
clc
a = 0;
b = 1;
c1 = 0.1;
c2 = 0.2;
d1 = 0.3;
d2 = 0.4;
ymin = -2;
ymax = 2;
ynum = 400;
ylist = linspace(ymin,ymax,ynum);
conv_1 = mat_conv2(ylist,d1,d2,c1,c2,a,b);
conv_2 = cont_conv2(ylist,d1,d2,c1,c2,a,b);
hold on
plot(ylist,conv_1,'DisplayName','conv()');
plot(ylist,conv_2','DisplayName','continuous');
legend('-DynamicLegend');
function res = fx(x,a,b)
res = sqrt(normpdf(x,a,b));
end
function res = mat_conv1(y,c1,c2,a,b)
dy = y(2)-y(1);
gx_1 = fx(y/c1,a,b);
gx_2 = fx(y/c2,a,b);
res = (1/(c1*c2))*conv(gx_1,gx_2,'same')*dy;
end
function res = mat_conv2(y,d1,d2,c1,c2,a,b)
dy = y(2)-y(1);
hx_1 = mat_conv1(y/d1,c1,c2,a,b);
hx_2 = mat_conv1(y/d2,c1,c2,a,b);
res = (1/(d1*d2))*conv(hx_1,hx_2,'same')*dy;
end
function res = cont_conv1(y,c1,c2,a,b)
dy = y(2)-y(1);
xx = -2:0.01:2;
for i=1:length(y)
gx_1 = fx(xx/c1,a,b);
gx_2 = fx((y(i)-xx)/c2,a,b);
yy = gx_1.*gx_2;
yy = trapz(xx,yy);
yy = (1/(c1*c2))*yy;
res(i) = yy;
end
end
function res = cont_conv2(y,d1,d2,c1,c2,a,b)
dy = y(2)-y(1);
xx = -2:0.01:2;
for i=1:length(y)
hx_1 = cont_conv1(xx/d1,c1,c2,a,b);
hx_2 = cont_conv1((y(i)-xx)/d2,c1,c2,a,b);
yy = hx_1.*hx_2;
yy = trapz(xx,yy);
yy = (1/(d1*d2))*yy;
res(i) = yy;
end
end
As you can see from the figure below the results match pretty neatly.
Now I only want to consider the positive domain. In other words, I want to take the convolution in the interval [0,2]. So I change ymin to zero and inside the functions cont_conv1 and cont_conv2 I start the vector xx from zero. The code is:
clear all
clc
a = 0;
b = 1;
c1 = 0.1;
c2 = 0.2;
d1 = 0.3;
d2 = 0.4;
ymin = 0;
ymax = 2;
ynum = 400;
ylist = linspace(ymin,ymax,ynum);
conv_1 = mat_conv2(ylist,d1,d2,c1,c2,a,b);
conv_2 = cont_conv2(ylist,d1,d2,c1,c2,a,b);
hold on
plot(ylist,conv_1,'DisplayName','conv()');
plot(ylist,conv_2','DisplayName','continuous');
legend('-DynamicLegend');
function res = fx(x,a,b)
res = sqrt(normpdf(x,a,b));
end
function res = mat_conv1(y,c1,c2,a,b)
dy = y(2)-y(1);
gx_1 = fx(y/c1,a,b);
gx_2 = fx(y/c2,a,b);
res = (1/(c1*c2))*conv(gx_1,gx_2,'same')*dy;
end
function res = mat_conv2(y,d1,d2,c1,c2,a,b)
dy = y(2)-y(1);
hx_1 = mat_conv1(y/d1,c1,c2,a,b);
hx_2 = mat_conv1(y/d2,c1,c2,a,b);
res = (1/(d1*d2))*conv(hx_1,hx_2,'same')*dy;
end
function res = cont_conv1(y,c1,c2,a,b)
dy = y(2)-y(1);
xx = 0:0.01:2;
for i=1:length(y)
gx_1 = fx(xx/c1,a,b);
gx_2 = fx((y(i)-xx)/c2,a,b);
yy = gx_1.*gx_2;
yy = trapz(xx,yy);
yy = (1/(c1*c2))*yy;
res(i) = yy;
end
end
function res = cont_conv2(y,d1,d2,c1,c2,a,b)
dy = y(2)-y(1);
xx = 0:0.01:2;
for i=1:length(y)
hx_1 = cont_conv1(xx/d1,c1,c2,a,b);
hx_2 = cont_conv1((y(i)-xx)/d2,c1,c2,a,b);
yy = hx_1.*hx_2;
yy = trapz(xx,yy);
yy = (1/(d1*d2))*yy;
res(i) = yy;
end
end
The result is given below.
As you can see the results are completely different. What am I doing wrong? How can I get rid of this discrepancy and use conv() to convolve the two functions in the positive domain?