- #1
jamie_m
- 14
- 0
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:
(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);
}
Last edited: