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

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"?

Answer 3

In line with the OPs request for a comparison of several different methods here's a comparison of six different ways to filter (or convolve) a data set x with a kernel h: convolution, correlation, the frequency domain method, a direct time-domain method such as might be programmed in C or Java, and a vectorized version such as would be common in Matlab or Octave. 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 Fourier Transforms that the product of the transforms is equal to the convolution of the time domain signals. The following calculate the Fourier 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 response 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}

And finally, here is a "vectorized" version such as might be programmed in Matlab

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}

How do all these methods stack up in terms of computation/performance?

h = RandomReal[{-1, 1}, 100];
x = RandomReal[{-1, 1}, 10^5];
n = Length[x] + Length[h] - 1;
xPad = PadRight[x, n]; 
{{"ListConvolve", First@Timing@ListConvolve[h, xPad, {1, 1}]},
 {"ListCorrelate", First@Timing@ListCorrelate[Reverse[h], xPad, {-1, -1}]},
 {"Fourier", First@Timing[
    ffth = Fourier[PadRight[h, n], FourierParameters -> {1, -1}];
    fftx = Fourier[PadRight[x, n], FourierParameters -> {1, -1}];
    InverseFourier[ffth fftx, FourierParameters -> {1, -1}]]},
 {"time domain", First@Timing[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]}]]},
 {"java", First@Timing[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]}];]},
 {"vectorized", First@Timing[
    Table[Inner[Times, Reverse[h], z[[i ;; i + Length[h] - 1]], Plus], {i, 1, Length[x]}]]}}

The output on my machine is:

{{"ListConvolve", 6.391073}, {"ListCorrelate", 6.495283}, {"Fourier", 0.212509}, 
 {"time domain", 2.737165}, {"java", 65.056947}, {"vectorized", 8.214995}}

Thus the ListConvolve, ListCorrelate and the vectorized versions are all about the same. The Do loop of the time domain method is significantly faster, and the Fourier method is the fastest by quite a margin. The direct loops of the "java" method are much slower.