# Convolution of two probability distributions using FFT

1. Nov 8, 2012

### jamie_m

I've been trying to code an algorithm to compute the convolution of two probability distributions. using the FFT. This relies on the "convolution theorem":

(p*q)[z] = FFT^{-1}(FFT(p) \cdot FFT(q))

However, when I test it using the distributions

p={0.1, 0.2, 0.3, 0.4}
q={0.4, 0.3, 0.2, 0.1}

the result that's output - {0.24, 0.22, 0.24, 0.3} - doesn't correspond to any definition or generalised definition of convolution that I can find.

(I'm fairly sure that one of the 0.24s, for instance, results from p[0]q[2] + p[1]q[1] + p[2]q[0] + p[3]q[3] being computed and divided by n)

Can someone well-versed in FFTs (and/or convolutions) shed any light on what might be going on? I enclose the source code, which uses C++ and the FFTW package, below:

Code (Text):

void convolution_of_two_probability_distributions_with_fft(long double* p, long double* q, long double* p_star_q, unsigned long int n)
{
fftwl_plan fpcp;
fftwl_plan fpcq;
fftwl_plan f_reverse;
unsigned long int ulx;

fftwl_complex* transformed_p = (fftwl_complex*)fftwl_malloc(sizeof(fftwl_complex) * n);
fftwl_complex* transformed_q = (fftwl_complex*)fftwl_malloc(sizeof(fftwl_complex) * n);
fftwl_complex* transformed_p_dot_transformed_q = (fftwl_complex*)fftwl_malloc(sizeof(fftwl_complex) * n);

fpcp = fftwl_plan_dft_r2c_1d(n, p, transformed_p, FFTW_ESTIMATE);
fpcq = fftwl_plan_dft_r2c_1d(n, q, transformed_q, FFTW_ESTIMATE);
f_reverse = fftwl_plan_dft_c2r_1d(n, transformed_p_dot_transformed_q, p_star_q, FFTW_ESTIMATE);

fftwl_execute(fpcp);
fftwl_execute(fpcq);

for (unsigned long int ulx=0; ulx < n; ulx++)
{
//transformed_p_dot_transformed_q[ulx] = transformed_p[ulx] * transformed_q[ulx];
//There's no overloaded * for fftwl_complex.
transformed_p_dot_transformed_q[ulx][0] = (transformed_p[ulx][0] * transformed_q[ulx][0]) - (transformed_p[ulx][1] * transformed_q[ulx][1]);
transformed_p_dot_transformed_q[ulx][1] = (transformed_p[ulx][0] * transformed_q[ulx][1]) + (transformed_p[ulx][1] * transformed_q[ulx][0]);
}

fftwl_execute(f_reverse);

for (ulx=0; ulx < n; ulx++)
{
p_star_q[ulx] /= n;
}

fftwl_destroy_plan(fpcp);
fftwl_destroy_plan(fpcq);
fftwl_destroy_plan(f_reverse);

fftwl_free(transformed_p);
fftwl_free(transformed_q);
fftwl_free(transformed_p_dot_transformed_q);
}

Last edited: Nov 8, 2012
2. Nov 9, 2012

### Stephen Tashi

Perhaps it's "circular convolution"

$(p*q)_n = \sum_{m=0}^3 { p_m q_{n-m} }$

e.g.
$(p*q)_0 = \sum_{m = 0}^3 {p_m} q_{0-m}$
$= p_0 q_0 + p_1 q_{-1} + p_2 q_{-2} + p_3 q_{-3}$

Where we regard $q_{-1} = q_3 , q_{-2}= q_2$ etc

$= (0.1)(0.4) + (0.2)(0.1) + (0.3)(0.2) + (0.4)(0.3) = 0.24$

3. Nov 9, 2012

### jamie_m

I think you're right - this is consistent with the existing results and with another pair of distributions I've tested it with. Thanks for that!