Discrete convolution in Maxima

I'd like to compute probability density functions for the sum of independent variables, which is the convolution of the density functions of the variables. Toward that end, I'll make a discrete approximation of each density function and compute the convolution of those approximations. I'll come back to that in a future blog post; for now I'll just work on the convolution.

Discrete convolution can be implemented via the fast Fourier transform, since FT(conv(x, y)) = FT(x)*FT(y), where the multiplication is taken element by element. So just take FT(x), FT(y), multiply them together, and take the inverse transform. Voila, there's your convolution.

In the following implementation of conv, I've taken some details into account. Maxima's fft can handle only inputs whose length is a power of 2. The exact length of the convolution is length(x) + length(y) - 1, so we'll just chop off any extra elements. Maxima's fft has a scaling convention which needs to be adjusted in order to get the correct result. Also, I'll call realpart to chop off any imaginary parts (these can appear due to round-off error).

Discrete convolution via fast Fourier transform

load (fft) $

conv (a, b) :=

block

([n : length(a) + length(b) - 1],

inverse_fft (fft (zero_extend (a, next2(n))) * fft (zero_extend (b, next2(n)))),

realpart (next2(n) * %%),

rest (%%, n - next2(n))) $

zero_extend (a, n) := makelist (if i <= length(a) then a[i] else 0, i, 1, n) $

next2 (n) := block ([i : ilog2 (n)], if n = 2^i then n else 2^(i + 1)) $

ilog2 (n) := block ([i : 0], while n > 1 do (n : floor (n/2), i : i + 1), i) $

Let's try it out. I'll make a couple of random inputs.

foo (n) := makelist (random (10), i, 1, n) $

[x, y] : [foo (17), foo (11)];

⇒

$\left[\right[2,2,4,5,4,1,9,5,8,3,5,5,0,6,9,0,9],[4,7,6,9,9,3,9,6,6,6,3\left]\right]$

fpprintprec : 4 $

conv (x, y);

⇒

$[8.0,22.0,42.0,78.0,111.0,122.0,172.0,212.0,241.0,293.0,326.0,310.0,332.0,294.0,351.0,321.0,327.0,354.0,255.0,231.0,243.0,132.0,171.0,126.0,81.0,54.0,27.0]$

Well, that's plausible. I'll check it more carefully later, but first, is conv symmetric, at least?

lmax (abs (conv (x, y) - conv (y, x)));

⇒

$0.0$

Discrete convolution via summation

I'll also define the discrete convolution as a summation; this is slow compared to the fast Fourier transform. There are at least a couple of ways to do that; in naive_conv_1, I'll be careful to define the summation range, and in naive_conv_2, I'll let the summation index go where it will, but I'll protect the inputs against out-of-bounds indices. These two are equivalent -- which you think is more comprehensible is a matter of preference.

naive_conv_1 (a, b) :=

makelist

(sum

(b[j - i + 1] * a[i],

i,

max (1, j - length(b) + 1),

min (j, length(a))),

j,

1,

length(a) + length(b) - 1) $

naive_conv_2 (a, b) :=

block

([f : lambda ([aa, ii], if 1 <= ii and ii <= length (aa) then aa[ii] else 0)],

makelist

(sum

(f (b, j - i + 1) * f (a, i),

i,

1,

length(a)),

j,

1,

length(a) + length(b) - 1)) $

Let's try these functions.

naive_conv_1 (x, y);

⇒

$[8,22,42,78,111,122,172,212,241,293,326,310,332,294,351,321,327,354,255,231,243,132,171,126,81,54,27]$

That's comforting. If I handled it right, these two functions are equivalent.

is (naive_conv_1 (x, y) = naive_conv_2 (x, y));

⇒

$\mathrm{true}$

Comparing naive and fast convolution implementations

So far, so good. Is the naive implementation the same as the fast convolution function?

lmax (abs (naive_conv_1 (x, y) - conv (x, y)));

⇒

$5.684·{10}^{-14}$

Terrific. I'll let that be all for now.