Example of a multidimensional Fourier transform

I want to convolve the 3-dimensional gaussian with the 3-dimensional exponent (I posted a question in RooFit section, but it seems that RooFit can’t do that directly). I have Fourier transforms of the exponent and of the gaussian (they are easily calculated analytically). I want to compute the inverse Fourier transform of their product (which would correspond to the convolution of the original pdfs). How could that be done using FFT?

It seems that I should take TFFTComplexReal, which will transform my complex-valued function to the real function.

However, I can’t understand where I shall provide function values (those of the Fourier transform I want to invert - I repeat that I have analytical values for that).

In the tutorials the only example deals with 1-dimensional FFT (same for the documentation).

Thank you.

Hi,
You can look at how it is implemented in TF1Convolution::MakeFFTConv, see ROOT: hist/hist/src/TF1Convolution.cxx Source File .
One would need to extend the code for multi-dimensional functions

If you cannot understand how to extend this code, let me know and I could try to do it

Cheers

Lorenzo

1 Like

Hi Lorenzo,

thanks for the example. Now I understand what is going on much better!

However, I have a question. There is a dedicated algorithm in fftw, Complex Multi-Dimensional DFTs. Won’t it be more optimal than doing array-wise FFTs you are suggesting here? I mean the computational accuracy and speed (the handling of arrays in ROOT I understand now better).

Upd: now I understand that I can use the 3-dimensional algorithm right from your example.

Cheers,
Yaroslav

Hi,
ROOT uses already the multi-dimensional algorithms for the FFT , in case of the convolution from real to complex and then from complex to real.
See for example ROOT: math/fftw/src/TFFTRealComplex.cxx Source File

Lorenzo

1 Like

Thanks Lorenzo!

Now I understand how to make a 3d FFT.
It is important to use

void TFFTRealComplex::SetPoint(const Int_t *ipoint, ...)

(not the method with 1-dimensional int ipoint) to set function values (ipoint is a multidimensional vector of bins).

Another example of a convolution (in 1d) is in RooFFTConvPdf::fillCacheSlice.

However, I found that I don’t know how to update RooFit variables, so I decided not to use a 3-dimensional convolution, but use 3 one-dimensional convolutions (probably that is fine from the statistical point of view). The distribution of a projection is again a convolution, because each projection is a sum of projections (x = x_gs + x_exp).

Memory is probably still another concern for 3d convolutions: @Wouter_Verkerke wrote in 2009, that “10000 x 100 x 100 bins = 100 million bins” and takes about 7.8 Gb of memory.

Yaroslav

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.