Register to reply

Convolution of two probability distributions using FFT

by jamie_m
Tags: convolution, convolution theorem, fft
Share this thread:
jamie_m
#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
Heat distributions help researchers to understand curved space
Professor quantifies how 'one thing leads to another'
Team announces construction of a formal computer-verified proof of the Kepler conjecture
Stephen Tashi
#2
Nov9-12, 12:07 AM
Sci Advisor
P: 3,300
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
#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