Hello,
This subject has been recently discussed on the mailing list but I think that it would be nice to replay the discussion with the tools of discourse (particularly math. formulas).
The subject is about the approximation of the convolution product of two causal functions x\rightarrow f(x) and x\rightarrow g(x), denoted here by h=f * g and
at discrete values of x given by x_i=i\delta, for integers i\geq 0 and \delta>0, given samples of f and g at x_i or Scilab functions computing their values at arbitrary values of x. A first remark is that even if f and g do not have this property, h is necessarily continuous and h(x_0)=h(0)=0. Then in both cases we have for i>0,
and the way we approximate the above inner integral determinates the ability to use a built-in Scilab function like convol(), which computes the discrete convolution between two vectors using the FFT.
Sampled case
In the first (sampled) case we only have the sampled values of f and g, that can be defined as f_i=f(x_i) and g_i=g(x_i) and if the inner integral is approximated by the left rectangle formula we have h_0=0 and for i>0
One can easily show that h_i is not equal to \delta times the i-th term of the discrete convolution f*g (where the sum is done up to i), unless f_0=0. Equivalently, if the right rectangle formula is used then we have the approximation
which is different from the i-th term of the discrete convolution f*g unless g_0=0. Hence, we can use convol() provided the rights arguments are given. We can illustrate this by an example, where f(x)=\mathbb{1}_{[0,1]}(x) and g(x)=\mathbb{1}_{[0,1]}(x)(1-x). We have in that case h(x)=\mathbb{1}_{[0,1]}(x)(x-x^2/2)+\mathbb{1}_{[1,2]}(x){(x-2)^2}/{2} and the following Scilab script compares h(x) and the approximations.
function y=f(x)
y=0+(x>=0 & x<=1);
end
function y=g(x)
y=(x>=0 & x<=1).*(1-x);
end
function y=h(x)
y=(x>=0 & x<=2).*((x<=1).*(x-1/2*x.^2)+(x>=1).*(0.5*(x-2).^2))
end
clf
x=linspace(0,2,1000);
Ī“=0.05;
xd=0:Ī“:1;
xdd=0:Ī“:2;
fd=f(xd);
gd=g(xd);
h_left=convol([0 fd(2:$)],gd)*Ī“;
h_right=convol([0 gd(2:$)],fd)*Ī“;
plot(xdd,h_left,'o',xdd,h_right,'o',x,h(x),"r")
legend("left rect.","right rect.","exact")
In fact we are tempted to improve this by taking the mean of these two approximations, which amounts to apply the trapezoidal rule, but the above graph shows that it wonāt improve the approximation in [1,2]. In fact in this interval the integrand t\rightarrow f(x-t)g(t) is not continuously differentiable, hence the trapezoidal rule does not attain its full order of approximation. We will see how to improve this in the next episode !
