# Approximate a convolution product **URL:** \ **Category:** Signal and image processing\ **Created:** [March 31, 2023, 5:16pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78 "2023-03-31T17:16:34Z")\ **Posts on this page:** 8\ **Page:** 1
**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\ **Post date:** [March 31, 2023, 5:16pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/1 "2023-03-31T17:16:34Z")
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 h(x)=\int\_0^x f(x-t)g(t)dt = \int\_0^x g(x-t)f(t)dt, 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, h(x\_i)=\sum\_{k=0}^{i-1}\int\_{x\_k}^{x\_{k+1}} f(x\_i-t)g(t)dt, 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 h\_i\approx \delta\sum\_{k=0}^{i-1}f\_{i-k}g\_k. 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 h\_i\approx \delta\sum\_{k=1}^{i}f\_{i-k}g\_k, 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. ```auto 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") ``` ![Figure n°0](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/b3302e20e642b0d3c31a975d06ebcd6570bf2a14.png) 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 ! ---
**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\ **Post date:** [April 3, 2023, 8:35am UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/2 "2023-04-03T08:35:11Z")
Hello, Like I said in episode 1 it is possible to improve the approximation we obtained in the sampled case, but in order to fully understand what could look as a “trick”, we need to consider the case when f(x) and g(x) can be evaluated at arbitrary values of x : ### Continuous case Since the integrand can be evaluated at x\_{k+1/2}=(x\_k+x\_{k+1})/2 we can use the midpoint rule, hence h(x\_i) can be approximated for i\>0 by \begin{align} h\_i&=\delta\sum\_{k=0}^{i-1} f(x\_i-x\_{k+1/2})g(x\_{k+1/2}), \\ &=\delta\sum\_{k=0}^{i-1} f(x\_{i-(k+1/2)})g(x\_{k+1/2}), \end{align} the i-th term of the discrete convolution of discrete sequences f and g given by f\_i=f(x\_{i+1/2}) and g\_i=g(x\_{i+1/2}) (the boundary value h\_0=0 is not computed by this formula). If we consider the same example as before, the following Scilab code computes the new approximation of x\rightarrow h(x). ```auto xdh=δ/2:δ:1; h_mid=[0 convol(f(xdh),g(xdh)) 0]*δ; plot(xdd,(h_left+h_right)/2,'o',xdd,h_mid,'o',x,h(x),'r') legend("trapezoidal rule","midpoint rule","exact") ``` ![Figure n°0](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/0bb8254b29462e61be083ec174d603f68b6bf706.png) We have added the approximation given by the trapezoidal rule. Compared to the latter, the midpoint rule does a perfect job. In the general case we could improve this by taking higher order quadrature formulas (see e.g. [https://en.wikipedia.org/wiki/Gaussian\_quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature)) if the regularity of the integrand diserves it, but that’s not the point here. ### Back to discrete case We cannot directly apply the midpoint rule, but we can approximate values of f and g at x\_{k+1/2} by linear interpolation: \begin{align}f(x\_{k+1/2})&\approx\frac{1}{2}(f(x\_{k})+f(x\_{k+1})),\\ g(x\_{k+1/2})&\approx\frac{1}{2}(g(x\_{k})+g(x\_{k+1})), \end{align} this gives the following results: ```auto fdh=(fd(1:$-1)+fd(2:$))/2; gdh=(gd(1:$-1)+gd(2:$))/2; plot(xdd,[0 convol(fdh,gdh)*δ 0],'o',x,h(x),"r") legend("midpoint-like rule","exact") ``` ![Figure n°1](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/aa958a7f69a7e33c1c8e6ef4f0201be6ec6c89f6.png) We have a perfect match (up to machine precision) again. In fact, this is due (as for the continous case) to the particular functions we have chosen, but this gives a good method in the general case. ---
**Author:** ![jbaudais](https://avatars.discourse-cdn.com/v4/letter/j/d78d45/32.png) [@jbaudais](https://scilab.discourse.group/u/jbaudais)\ **Post date:** [April 4, 2023, 3:36pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/3 "2023-04-04T15:36:15Z")
Hi, Take care when you use digital signal and analog signal. Both digital signal processing and analog signal processing give the same iff you correctly applied the Nyquist-Shannon sampling theorem. It is not enough to pick up samples anywhere to get digital signal processing: The rectangular or trapezoidal rules are approximations that fit more or less well depending one the processing used. But, Fourier transform, convolution, filter… are well defined in both analog and digital spaces thanks to Nyquist-Shannon (and others). There exists many articles, books… dedicated to that. I have in mind new approach on analog-digital equivalence presented on conference at the beginning of this millennium… If I find the reference,I will share them. ---
**Author:** ![jbaudais](https://avatars.discourse-cdn.com/v4/letter/j/d78d45/32.png) [@jbaudais](https://scilab.discourse.group/u/jbaudais)\ **Post date:** [April 4, 2023, 3:53pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/4 "2023-04-04T15:53:38Z")
Found: [EPFL | Biomedical Imaging Group | Unified Theory](http://bigwww.epfl.ch/publications/unser0508.html) but only as a starting point, I’m not able to replay the talk ☹ ---
**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\ **Post date:** [April 4, 2023, 4:56pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/5 "2023-04-04T16:56:33Z")
Hello jean-Yves, Thanks for your remark but that’s not (completely) the point here (maybe the topic should have been started elsewhere).The initial topic in the ML was about the correct computation of a convolution product between probability distributions in a statistics problem, hence, this is not signal processing. ---
**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\ **Post date:** [April 4, 2023, 5:26pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/6 "2023-04-04T17:26:10Z")
The talk seems VERY interesting (I just read the abstract), thanks ! ---
**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\ **Post date:** [April 4, 2023, 8:05pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/7 "2023-04-04T20:05:24Z")
Here are the links to the paper of M. User: > **[Cardinal Exponential Splines: Part I—Theory and Filtering Algorithms](https://infoscience.epfl.ch/record/63127)** > > Causal exponentials play a fundamental role in classical system theory. Starting from those elementary building blocks, we propose a complete and self-contained signal processing formulation of exponential splines defined on a uniform grid. We... > **[Cardinal Exponential Splines: Part II—Think Analog, Act Digital](https://infoscience.epfl.ch/record/63128)** > > By interpreting the Green-function reproduction property of exponential splines in signal processing terms, we uncover a fundamental relation that connects the impulse responses of allpole analog filters to their discrete counterparts. The link is... S. ---
**Author:** ![heinznabielek](https://avatars.discourse-cdn.com/v4/letter/h/e47c2d/32.png) [@heinznabielek](https://scilab.discourse.group/u/heinznabielek)\ **Post date:** [April 6, 2023, 8:08pm UTC](https://scilab.discourse.group/t/approximate-a-convolution-product/78/8 "2023-04-06T20:08:32Z")
EXAMPLE: - We have an object consisting of many (~10,000) small components that can fail in a statistical way during long-term operation at extreme conditions. - Our primary failure model #1 is described by Φ(t) going monotonically from zero to one at times from t=0 onwards. In the computer, this is realized in n timesteps. - ψ(t) is the time derivative of Φ. - Mechanism #2: There is a secondary failure mechanism that starts only when a component has failed by mechanism 1 and occurs after some delay. The delay function is g(t) and final expression for the evolution of the mechanism #2 failure is given by: - F(t) = Int\_{0}^{t} {ψ(τ) . g(t-τ)} dτ , a classical convolution integral where ψ=dΦ/dt is the failure rate in mechanism #1. To solve this problem, the Scilab code example (below) uses the routine “conv” successfully thanks to the support of Stéphane Mottelet and the Scilab community. ```auto m=2; TT=3e6; k=1E-7; // failure rate constant (s-1) in mechanism #1 n=500; // number of time steps in hourly intervals t=(1:n)'; // timesteps in hours ts= 3600*t; // timesteps in seconds PHI= 1 - exp(-((k*ts)^m)); // failure fraction in mechanism #1 plot("nl", t, PHI, 'r--'); xgrid function y=f(t) //failure rate by derivative mechanism #1 ts=3600*t;y=3600*m*k*((k*ts)^(m-1)) .* exp(-((k*ts)^m)); endfunction function y=g(t) // delay to failure mechanism #2 y= 1 - exp(-3600*t/TT); endfunction h=[0 conv(f((1/2:n)), g((1/2:n))) 0]; th=(0:n+n); plot(th(1:501), max(1e-8, h(1:501)), 'b-'); legend('primary failure fraction', 'secondary failure fraction', 4); xlabel('time (hours)'); ylabel('failure fraction'); ``` ![ConvolExampleHeinzNabielek](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/06191080e663add01a99d0bf208714eaa717ae9c.png)