## Convolution of two probability distributions using FFT

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:
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);
}

 Recognitions: Science Advisor 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$

 Quote by Stephen Tashi Perhaps it's "circular convolution" $(p*q)_n = \sum_{m=0}^3 { p_m q_{n-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!

 Tags convolution, convolution theorem, fft