Convolution of two probability distributions using FFT

[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:

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

 

[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

 

[jamie_m]

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!