Convolution is mathematically well understood, both in the continuous and discrete domains. There is a simple formula to comupte it. However, I recently needed to implement this in code and found that it was not all too easy. In this post, I tried to come up with a formula that can be used for programming convolution. In later posts about this, I will also try to analyze the computational complexity and spruce it up with some real working convolution code.

Firstly, what is discrete convolution? There is a straight-forward formula :

\displaystyle  h[n] \equiv (f*g)[n] = \sum_{k=-\infty}^{k =\infty} f[k]g[n-k]

For generality we will assume that {f} and {g} are discrete functions (i.e. functions from the Integers to the Reals). For example, {f} might be such that it is defined from {fl} to {fh} and it is zero outside. This is the most general representation for the signal. Similarly, let {g} be defined on {[gl,gh]} and zero outside. When we are asked to convolve these two signals, it is apparent that {fl,gl,fh,gh} would play some role in the computation. However, MATLAB for example, does not care about these values. So, it is slightly confusing how convolution is carried out. Let us examine.

In the definition of convolution, what are the values of {n} that are interesting? Before thinking about that, it is better to understand what values of {k} are important for our considerations. Fix any {n} and note that

  • { fl \leq k \leq fh}, because outside of this range, {f[k]} is zero and hence does not add to the convolution
  • Similarly, {gl \leq n - k \leq gh} which implies that {n - gh \leq k \leq n - gl}

Clearly, if {[fl,fh] \cap [n-gh, n-gl] = \phi}, the convolution is zero and we do not have to care about such an {n}. Otherwise, one can immediately see that the interesting values of {k} are :

\displaystyle [fl,fh] \cap [n-gh, n-gl] = \max (fl, n- gh) \leq k \leq \min (fh, n - gl)

Let us rewrite the formula again (for those {n} such that {[fl,fh] \cap [n-gh, n-gl] \neq \phi}):

\displaystyle  h[n] \equiv (f*g)[n] = \sum_{k= \max (fl, n- gh)}^{\min (fh, n - gl)} f[k]g[n-k]

Now, what are the values for {n} for which {[fl,fh] \cap [n-gh, n-gl] \neq \phi}? Consider, {  n = fl + gl}, then

\displaystyle  \begin{array}{rcl}  [fl,fh] \cap [fl+gl-gh, gl+gl-gl]  &=& \max (fl, fl+gl- gh) \leq k \leq \min (fh, fl+gl - gl) \\ 		 		 &=& fl \leq k \leq fl \\ 				 &=& fl \end{array}

When {n < fl + gl}, {\min (fh, n - gl)  < fl} and hence the convolution is zero. The upper limit on {n} can similarly be found to be {n = fh + gh}. Therefore, the full definition of convolution useful for computer implementation is:

\displaystyle  h[n] \equiv (f*g)[n] = \sum_{k= \max (fl, n- gh)}^{\min  (fh, n - gl)} f[k]g[n-k] ~~\forall n \in [fl+gl,fh+gh]

Upon further thought, it can be easily seen that {fl,fh,gl,gh} are not important for the core convolution. If the signal {(f, [fl,fh])} is replaced by {(a, [0, fh-fl])} and {(g,  [gl,gh])} is replaced by {(b, [0,  gh-gl])} there would no difference in the computation of the main convolution function. Only the range of {n} for which the convolution is defined change. This is the primary reason why MATLAB does not care about these things.