Discrete Convolution

Question

Ask to compute the convolution of 2 lists, I managed to do so, with what I feel is rather heavy :

Let my 2 lists be :

a = {1,2,3,4} b = {1,1,1,1,1,1};

The below function adds 0s on each part of one list given the Length of the other

bpad[listA_, listB_] :=
 Flatten@{
          Table[0, {Length@listA - 1}],
          listB,
          Table[0, {Length@listA - 1}]
          }

bpad[a,b]

Out = {0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0}

The following ones gives me the lists part I need to do the dot product with :

convParts[listA_, listB_] :=
         Table[
               Range[i, i + Length@listA - 1], 
               {i, Range[Length@bpad[listA, listB] - Length@listA + 1]}]

convParts[a, b]

Out= {{1, 2, 3, 4}, {2, 3, 4, 5}, {3, 4, 5, 6}, {4, 5, 6, 7}, {5, 6, 7,8}, {6, 7, 8, 9}, {7, 8, 9, 10}, {8, 9, 10, 11}, {9, 10, 11, 12}}

Finally the convolution itself :

convolve[listA_, listB_] :=
         Total@(bpad[listA, listB][[#]]*Reverse@listA) & /@ 
         convParts[listA, listB]

convolve[a,b]

Out= {1, 3, 6, 10, 10, 10, 9, 7,4}

How could I improve my solution, at each function or overall level.

Doing this for a class, it is, once again my opportunity to advertise for Mathematica against Matlab. I like in those case to show several solution

Answer 1

You could use ListConvolve:

ListConvolve[a, b, {1, -1}, 0]

concerning the padding:

ArrayPad[b, 3, 0]

And you could use Partition for the second of your steps:

Partition[Range[Length[ArrayPad[b, 3, 0]]], 3, 1]

Answer 2

In line with the OPs request for a comparison of several different methods here's a comparison of four different ways to filter (or convolve) a data set x with a kernel h: convolution, correlation, the frequency domain method, and a direct time domain method. The only difference (other than numerical factors) is in the way edge conditions are handled with padding. First we set up the data:

h = {1, -1, 2, -2, 3, -3}; 
x = {1, 2, 3, 4, 5, 6, -5, -4, -3, -2, -1};  
n = Length[x] + Length[h] - 1;
xPad = PadRight[x, n];

In the convolution method, the kernel h is thought of as the impulse response of a linear time-invariant system and the x is thought of as the input to that system. The convolution yConv is then the output of the system.

yConv = ListConvolve[h, x, {1, 1}, 0];
yConvPad = ListConvolve[h, xPad, {1, 1}]
{1, 1, 3, 3, 6, 6, -6, 6, -18, 6, -30, 6, 5, 5, 3, 3}

In the correlation method, the kernel h is thought of as a marker or mask and x is thought of as the data that is to be examined. The correlation yCorr is then how much like x the kernel is at each place in the sequence.

yCorr = ListCorrelate[Reverse[h], x, {-1, -1}, 0];
yCorrPad = ListCorrelate[Reverse[h], xPad, {-1, -1}]
{1, 1, 3, 3, 6, 6, -6, 6, -18, 6, -30, 6, 5, 5, 3, 3}

The Fourier method exploits the fact from Foureir Transforms that the product of the transfoms is equal to the convolution of the time domain signals. The following calculate the Fourer transform of h (ffth) and the Fourier transform of x (fftx), after padding to the same length. The element-by-element product is then inverse transformed, giving yFourier. which is numerically the same as the above methods.

ffth = Fourier[PadRight[h, n], FourierParameters -> {1, -1}];
fftx =  Fourier[PadRight[x, n], FourierParameters -> {1, -1}];
yFourier = InverseFourier[ffth fftx, FourierParameters -> {1, -1}]
{1., 1., 3., 3., 6., 6., -6., 6., -18., 6., -30., 6., 5., 5., 3., 3.}

In the time-domain method, the output of the system with impulse reponse h is calculated once for each time k, as the input takes on all values in x.

z = PadLeft[x, n];
yTim = ConstantArray[0, Length[x]];
Do[
   yTim[[k]] = Total[Reverse[h] z[[k ;; k + Length[h] - 1]]];
     , {k, 1, Length[x]}]
yTim
{1, 1, 3, 3, 6, 6, -6, 6, -18, 6, -30}

Here's a time domain version that's like one might program it in Java or C. Normally one would truncate the initial string of zeros.

z = PadLeft[x, n];
yJav = ConstantArray[0, n + 1];
Do[ 
  Do[
     yJav[[k]] = yJav[[k]] + h[[j]] z[[k - j]];
            , {j, 1, Length[h]}];
       , {k, Length[h] + 1, Length[x] + Length[h]}];
yJav
{0, 0, 0, 0, 0, 0, 1, 1, 3, 3, 6, 6, -6, 6, -18, 6, -30}

Table[Inner[Times, Reverse[h], z[[i ;; i + Length[h] - 1]], Plus], {i, 1, 
  Length[x]}]
{1, 1, 3, 3, 6, 6, -6, 6, -18, 6, -30}

Answer 3

Use the following:

ListConvolve[a, b, {1, -1}, 0]

(* ==> {1, 3, 6, 10, 10, 10, 9, 7, 4} *)

This says: align the first element of b with the 1st element of a, align the last element of b with the -1st (i.e. last) element of a, pad with 0s if necessary.

Have you tried searching the docs for "convolution"?