0

I am trying to do a 3D FFT with the FFTW library, but I have some difficulties with the inverse transformation.

At first I do the foreword transformation via:

fftwf_plan_dft_3d(_dimensions[0], _dimensions[1], _dimensions[2], (fftwf_complex*)_inputBuffer, (fftwf_complex*)_outputBuffer, FFTW_FORWARD, FFTW_ESTIMATE);

Although my data is real data I am using the complex to complex transformation, as want to replace it later by an opencl fft which only supports complex to complex transformations.

In 3D fourier space I do a very simple low pass filter:

for all x, y, z:

// global position of the current bin
int gid = (y * w + x) + (z * w * h);

// position of the symmetric bin
vec3 conPos(M - x - 1, N - y - 1, L - z - 1);

// global position of the symmetric element
int conGid = (conPos.y * w + conPos.x) + (conPos.z * w * h);

if (sqrt(x * x + y * y + z * z) > 500)
{
    complex[gid].real = 0.0f;
    complex[gid].imag = 0.0f;
    complex[conGid].real = 0.0f;
    complex[conGid].imag = 0.0f;
}

At last the inverse transformation:

fftwf_plan_dft_3d(_dimensions[0], _dimensions[1], _dimensions[2], (fftwf_complex*)_inputBuffer, (fftwf_complex*)_outputBuffer, FFTW_BACKWARD, FFTW_ESTIMATE);
// normalization ...

The Result is not as I would expect it. After the inverse transformation the imaginary parts are not all zero as they are supposed to be.

As far as I see it, after the forward transformation of real data only the half of the total buffer size is used and there are no conjugate complex values in the other half. (see: c2c with real data) If this is the case I have to compute them on my own before the backward transformation, but I could not found a hint in the fftw docs which half is computed and which is not.

I have written a very simple 2D-Test-Case to view this symmetry in the fourier space:

int w = 4;
int h = 4;
int size  = w * h;

cl_float rawImage[16] = ...; // loading image

fftwf_complex *complexImage = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * size);
fftwf_complex *freqBuffer = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * size);

for (int i = 0; i < size; i++)
{
    complexImage[i][0] = rawImage[i]; complexImage[i][1] = 0.0f;
}

fftwf_plan forward = fftwf_plan_dft_2d(w, h, complexImage, freqBuffer, FFTW_FORWARD, FFTW_ESTIMATE);

fftwf_execute(forward);

for (int y = 0; y < h; y++)
{
    for (int x = 0; x < w; x++)
    {
        int gid = y * w + x;
        qDebug() << gid << "real:" << freqBuffer[gid][0] << "imag:" << freqBuffer[gid][1];
    }
}

This gives me the following output:

gid
0    real 3060 imag 0 
1    real 510 imag 510 
2    real 0 imag 0 
3    real 510 imag -510 
4    real 510 imag 510 
5    real 0 imag -510 
6    real 0 imag 0 
7    real -510 imag 0 
8    real 0 imag 0 
9    real 0 imag 0 
10   real 0 imag 0 
11   real 0 imag 0 
12   real 510 imag -510 
13   real -510 imag 0 
14   real 0 imag 0 
15   real 0 imag 510

As far as I see it there are no symmetric values. Why?

It would be nice if someone could give me a hint.

Greetings

Wolf

DerHandwerk
  • 99
  • 3
  • 10

2 Answers2

1

That link is misleading. The DFT of a real-only signal does not (in general) result in half the output samples being zero. It simply imposes a (conjugate) symmetry on them.

So in your filter code, you are only manipulating half of the output values that you should be. Every time you manipulate output bin n, you also need to manipulate bin N-n (where N is the length of the DFT), in order to maintain the symmetry that will give you a real-only result when you apply the inverse DFT.

My advice would be to tackle a much simpler problem first - a 1D filter. Once you have that correct, then it should be easy to scale up to 3D.

Oliver Charlesworth
  • 267,707
  • 33
  • 569
  • 680
  • That sounds reasonable. Does that mean in case of 1D FFT I only need to got from n = 0 to N / 2 since in every step I filter complex[n] and complex[N-n]? – DerHandwerk Apr 29 '12 at 19:46
  • @DerHandwerk: Not sure I understand what you're asking. A length-N 1D DFT will only have N/2+1 independent output samples; the relationship is y[n] = y[N-n]* (where "*" denotes [complex conjugate](http://en.wikipedia.org/wiki/Complex_conjugate)). – Oliver Charlesworth Apr 29 '12 at 19:52
  • Ok as far as I see it for 1D FFT my low pass has to do the following: `if amplitude > value then y[n] = 0; y[N - n - 1] = 0; end` for n = 0 to N - 1 – DerHandwerk Apr 29 '12 at 20:32
1

If you want a strictly real result after an inverse FFT (minus the usual numerical noise from using finite size arithmetic), you have to make sure the the input data you feed a full IFFT is completely conjugate symmetric (last half of the vector is the mirrored complex conjugate of the first half). It does not appear that you have forced your data to be that way.

hotpaw2
  • 70,107
  • 14
  • 90
  • 153