Convolution of two probability distributions using FFT


by jamie_m
Tags: convolution, convolution theorem, fft
jamie_m
jamie_m is offline
#1
Nov8-12, 05:12 PM
P: 15
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);
}
Phys.Org News Partner Mathematics news on Phys.org
Math modeling handbook now available
Hyperbolic homogeneous polynomials, oh my!
Researchers help Boston Marathon organizers plan for 2014 race
Stephen Tashi
Stephen Tashi is online now
#2
Nov9-12, 12:07 AM
Sci Advisor
P: 3,177
Perhaps it's "circular convolution"

[itex] (p*q)_n = \sum_{m=0}^3 { p_m q_{n-m} } [/itex]

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

Where we regard [itex] q_{-1} = q_3 , q_{-2}= q_2 [/itex] etc

[itex] = (0.1)(0.4) + (0.2)(0.1) + (0.3)(0.2) + (0.4)(0.3) = 0.24 [/itex]
jamie_m
jamie_m is offline
#3
Nov9-12, 06:36 AM
P: 15
Quote Quote by Stephen Tashi View Post
Perhaps it's "circular convolution"

[itex] (p*q)_n = \sum_{m=0}^3 { p_m q_{n-m} } [/itex]
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!


Register to reply

Related Discussions
Convolution of distributions on integers Set Theory, Logic, Probability, Statistics 0
Convolution of densities and distributions General Math 8
Can we deduce the product of distributions via convolution theorem Calculus 0
Distributions: Convolution product Calculus & Beyond Homework 4
Computing distributions by using convolution. Set Theory, Logic, Probability, Statistics 2