1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Convolution of two probability distributions using FFT

  1. Nov 8, 2012 #1
    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 (Text):

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


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


        for (ulx=0; ulx < n; ulx++)
            p_star_q[ulx] /= n;


    Last edited: Nov 8, 2012
  2. jcsd
  3. Nov 9, 2012 #2

    Stephen Tashi

    User Avatar
    Science Advisor

    Perhaps it's "circular convolution"

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

    [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]
  4. Nov 9, 2012 #3
    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!
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook