C++ FFTW forward backward DFTvalues get wrapped

Question

Hello StackOverflow community, i have a problem with the dft algorithm of the fftw library. All i want to do is to transform a certain pattern forward and backward to receive the input pattern again, of course there will be some sort of filtering in between the transformations later on.

So, what my program does atm is:

1. Create a test signal
2. Filter or "window" the test signal with a value of 1.0 or 0.5
3. Copy the test signal to a fftw_complex data type
4. Perform a forward and backward dft
5. Calculate the magnitude, which is called phase here
6. Copy and adjust data for display purposes, and finally display the images via OpenCV

My problem is that when is use no filtering my backward transformed image is wrapped somehow and i can't calculate the correct magnitude, which should be indentical to my input image / test signal. When i set the fitler/"window" to a value of 0.5 the backward transformation works fine, but my input image is just half as bright as it should be.

The following image illustrates my problem: (from top left to bottom right) 1. Input signal, 2. Real part of backward transformation, 3. From backward transformated data calculated magnitude, 4. Input signal multiplied with 0.5, 5. Real part of backward transformation, 6. From backward transformated data calculated magnitude. http://imageshack.com/a/img538/5426/nbL9YZ.png

Does anybody have an idea why the dft performs in that way?! It's kind of strange...

My code looks like this atm:

/***** parameters **************************************************************************/
int     imSize                                          = 256;
int     imN                                             = imSize * imSize;

char*   interferogram                                   = new char[imN];
double* spectrumReal                                    = new double[imN];
double* spectrumImaginary                               = new double[imN];
double* outputReal                                      = new double[imN];
double* outputImaginary                                 = new double[imN];
double* phase                                           = new double[imN];

char*   spectrumRealChar                                = new char[imN];
char*   spectrumImaginaryChar                           = new char[imN];
char*   outputRealChar                                  = new char[imN];
char*   outputImaginaryChar                             = new char[imN];
char*   phaseChar                                       = new char[imN];

Mat     interferogramMat                                = Mat(imSize, imSize, CV_8U, interferogram);
Mat     spectrumRealCharMat                             = Mat(imSize, imSize, CV_8U, spectrumRealChar);
Mat     spectrumImaginaryCharMat                        = Mat(imSize, imSize, CV_8U, spectrumImaginaryChar);
Mat     outputRealCharMat                               = Mat(imSize, imSize, CV_8U, outputRealChar);
Mat     outputImaginaryCharMat                          = Mat(imSize, imSize, CV_8U, outputImaginaryChar);
Mat     phaseCharMat                                    = Mat(imSize, imSize, CV_8U, phaseChar);


/***** compute interferogram ****************************************************************/
fill_n(interferogram, imN, 0);
double value = 0;
double window = 0;

for (int y = 0; y < imSize; y++)
{
    for (int x = 0; x < imSize; x++)
    {
        value = 127.5 + 127.5 * cos((2*PI) / 10000 * (pow(double(x - imSize/2), 2) + pow(double(y - imSize/2), 2)));

        window = 1;
        value *= window;

        interferogram[y * imSize + x] = (unsigned char)value;
    }
}


/***** create fftw arays and plans **********************************************************/
fftw_complex*       input;
fftw_complex*       spectrum;
fftw_complex*       output;
fftw_plan           p_fw;
fftw_plan           p_bw;

input               = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
spectrum            = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
output              = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
p_fw                = fftw_plan_dft_2d(imSize, imSize, input, spectrum, FFTW_FORWARD, FFTW_ESTIMATE);
p_bw                = fftw_plan_dft_2d(imSize, imSize, spectrum, output, FFTW_BACKWARD, FFTW_ESTIMATE);


/***** copy data ****************************************************************************/
for (int i = 0; i < imN; i++)
{
    input[i][0] = double(interferogram[i]) / 255.;
    input[i][1] = 0.;
    spectrum[i][0] = 0.;
    spectrum[i][1] = 0.;
    output[i][0] = 0.;
    output[i][1] = 0.;
}


/***** FPS algorithm ************************************************************************/
fftw_execute(p_fw);

fftw_execute(p_bw);

for (int i = 0; i < imN; i++)
{
    phase[i] = sqrt(pow(output[i][0], 2) + pow(output[i][1], 2));
}


/***** copy data ****************************************************************************/
for (int i = 0; i < imN; i++)
{
    spectrumReal[i] = spectrum[i][0];
    spectrumImaginary[i] = spectrum[i][1];

    outputReal[i] = output[i][0] / imN;
    outputImaginary[i] = output[i][1];
}

SaveCharImage(interferogram, imN, "01_interferogram_512px_8bit.raw");
SaveDoubleImage(spectrumReal, imN, "02_spectrum_real_512px_64bit.raw");
SaveDoubleImage(spectrumImaginary, imN, "03_spectrum_imaginary_512px_64bit.raw");
SaveDoubleImage(outputReal, imN, "03_output_real_512px_64bit.raw");

DoubleToCharArray(spectrumReal, spectrumRealChar, imSize);
DoubleToCharArray(spectrumImaginary, spectrumImaginaryChar, imSize);

DoubleToCharArray(outputReal, outputRealChar, imSize);
DoubleToCharArray(outputImaginary, outputImaginaryChar, imSize);

DoubleToCharArray(phase, phaseChar, imSize);


/***** show images **************************************************************************/

imshow("interferogram", interferogramMat);
imshow("spectrum real", spectrumRealCharMat);
imshow("spectrum imaginary", spectrumImaginaryCharMat);
imshow("out real", outputRealCharMat);
imshow("out imaginary", outputImaginaryCharMat);
imshow("phase", phaseCharMat);

int key = waitKey(0);

Answer 1

Here are some lines of your code :

char*   interferogram  = new char[imN];
...
double value = 0;
double window = 0;

for (int y = 0; y < imSize; y++)
{
   for (int x = 0; x < imSize; x++)
   {
      value = 127.5 + 127.5 * cos((2*PI) / 10000 * (pow(double(x - imSize/2), 2) + pow(double(y - imSize/2), 2)));

      window = 1;
      value *= window;

      interferogram[y * imSize + x] = (unsigned char)value;
   }
}

The problem is that a char is between -128 and 127, while unsigned char ranges from 0 to 255. In interferogram[y * imSize + x] = (unsigned char)value;, there is an implicit cast to char.

It does not affect the output if window=0.5, but it triggers a change if window=1 as value becomes higher than 127. This is exactly the problem that you noticed in your question !

It does not affect the first displayed image since CV_8U corresponds to unsigned char : interferogram is therefore cast back into a unsigned char*. Take a look at Can I turn unsigned char into char and vice versa? to know more about char to unsigned char cast.

The problem occurs at input[i][0] = double(interferogram[i]) / 255.; : if window=1, interferogram[i] may be negative and input[i][0] becomes negative.

Change all char to unsigned char and it should solve the problem.

You may also change

outputReal[i] = output[i][0] / imN;
outputImaginary[i] = output[i][1];

for

outputReal[i] = output[i][0];
outputImaginary[i] = output[i][1];

Calls to fftw seems to be fine.