2D FFT Frequency Cut

ROOT
1

Hello,

I’ve successfully performed a FFT then iFFT on a 2D histogram as per:

Thanks John!

Now I would like to implement the ability to cut out selected frequencies. This was simple in the 1D case, where each element of the real and imaginary coefficient arrays are simply the coefficients of the corresponding frequency, when you use GetPointsComplex(double *re, double *im). Then you can make the desired cut by setting the corresponding array element to 0, and performing the backward transform on the modified array.

In the 2D case, however, I can’t figure out how the arrays correspond to frequency. Plotting the magnitude histogram, as outputed by the FFT function, makes sense, with a point (delta function) in 2D space corresponding to the frequency of the sine wave(s). But I can’t seem to recreate this histogram using the real and imaginary coefficient arrays. In other words, I don’t know how the arrays are filled so I can’t make a specific frequency cut.

I’ve attached my script which creates a 2D histogram using a sine function. It draws the image, the magnitude histogram (as output by ->FFT), and the real/imaginary/magnitude arrays as I have filled them (which may be incorrect). Finally, the unmodified coefficients are used to perform the iFFT and the resulting image is drawn (which is the same as the original image as expected).

Can anyone help me determine how the coefficient arrays are defined for the 2D FFT (how do they correspond to real frequencies) so that I can cut out frequencies from x and/or y?

Thank you,
Patrick

2DFFTCut.c (3.44 KB)

2

Hello,
in case of a multi-dimensional transform not all the coefficient are copied
due to the Hermitian symmetry ( re[i] = re[n-i] and im[i] = - im[n-i] ).

The convention used is explained in

fftw.org → Multi-Dimensional DFTs of Real Data

Therefore to get the right coefficient in the histogram, I would reccomend you to do as following, by using the method TVirtualFFT::GetPointComplex which retrieve the value point by point and deals correctly with the data organization.
If you need to use the arrays you should look in the implementation of that method
void TFFTRealComplex:GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput)

For example, you can then fill the histograms as follow:

  TH1 *re_coef = new TH2F("re_coef","re_coef",nBins[0],min[0],max[0],nBins[1],min[1],max[1]);
  TH1 *im_coef = new TH2F("im_coef","im_coef",nBins[0],min[0],max[0],nBins[1],min[1],max[1]);

  //Fill 2D histogram with real and imaginary coefficients 

   int ipoint[2]; 
  for (int i=0; i<nBins[0]; i++) {
    for (int j=0; j<nBins[1]; j++) {
       ipoint[0] = i; 
       ipoint[1] = j; 
       double a,b; 
       fft->GetPointComplex(ipoint, a, b);
       im_coef->SetBinContent(i+1,j+1, b);
       re_coef->SetBinContent(i+1,j+1, a);
    }
  }

by doing this the output looks reasonable to me

Best Regards

Lorenzo

3

Thanks for the advice, I got the histograms out correctly. However, now I’m having trouble reading them back in and performing the iFFT with them.

But besides that, I realized it is not a trivial task to cut frequencies from a 2D FFT. I will put this on hold (perhaps indefinitely) as the 1D projection suffices for now.