inverse fft of fft not returning expected data

Question

I'm trying to make sure FFTW does what I think it should do, but am having problems. I'm using OpenCV's cv::Mat. I made a test program that, given a Mat f, computes ifft(fft(f)) and compares the result to f. I would expect the difference between the two to be negligible, but there's a strange pattern in the data..

In this case, f is initialized to be an 8x8 array of floats with positive values less than 1.

Here's my test program code:

Mat f = .. //populate f
if (f.type() != CV_32FC1)
    DLOG << "Bad f type";
const int y = f.rows;
const int x = f.cols;
double* input = fftw_alloc_real(y * 2*(x/2 + 1));
// forward fft
fftw_plan plan = fftw_plan_dft_r2c_2d(x, y, input, (fftw_complex*)input, FFTW_MEASURE);
// inverse fft
fftw_plan iplan = fftw_plan_dft_c2r_2d(x, y, (fftw_complex*)input, input, FFTW_MEASURE);

// populate fftw data from f
for (int yi = 0; yi < y; ++yi)
{
    const float* yptr = f.ptr<float>(yi);
    for (int xi = 0; xi < x; ++xi)
        input[yi*x + xi] = (double)yptr[xi];
}

fftw_execute(plan); 
fftw_execute(iplan);

// put data into another cv::Mat for comparison
Mat check(y, x, f.type());
for (int yi = 0; yi < y; ++yi)
{
    float* yptr = check.ptr<float>(yi);
    for (int xi = 0; xi < x ; ++xi)
        yptr[xi] = (float)input[yi*x + xi];
}

DLOG << Util::summary(f, "f");
DLOG << f;
DLOG << Util::summary(check, "check");
DLOG << check;
Mat diff = f*x*y - check;
DLOG << Util::summary(diff, "diff");
DLOG << diff;

Where DLOG is my logger and Util::summary(cv::Mat m) just prints passed string and the dimensions, channels, min, and max of the mat.

Here's what the data looks like (output):

f: rows:8 cols:8 chans:1 min:0.00257996 max:0.4  
[0.050668437, 0.04509116, 0.033668514, 0.10986148, 0.12855141, 0.048241843, 0.12613985,.09731093;  
0.028602425, 0.0092236707, 0.037089188, 0.118964, 0.075040311, 0.40000001, 0.11959606, 0.071930833;  
0.0025799556, 0.051522054, 0.22233701, 0.052993439, 0.032000393, 0.12673819, 0.015244827, 0.044803992;  
0.13946071, 0.019708242, 0.0112687, 0.047459811, 0.019342113, 0.030085485, 0.018739942, 0.0098618753;  
0.041809395, 0.029681522, 0.026837418, 0.16038358, 0.29034778, 0.17247421, 0.1789207, 0.042179305;  
0.025630442, 0.017192598, 0.060540862, 0.1854037, 0.21287154, 0.04813192, 0.042614728, 0.034764063;  
0.0030835248, 0.018511582, 0.0071733585, 0.017076733, 0.064545207, 0.0026390438, 0.088922881, 0.045725599;
0.12798512, 0.23215951, 0.027465452, 0.03174505, 0.04352935, 0.025079668, 0.044403922, 0.035459157]

check: rows:8 cols:8 chans:1 min:-3.26489 max:25.6  
[3.24278, 2.8858342, 2.1547849, 7.0311346, 8.2272902, 3.0874779, 8.0729504, 6.2278996;  
0.30818239, 0, 2.373708, 7.6136961, 4.8025799, 25.6, 7.6541481, 4.6035733;  
0.16511716, 3.2974114, -3.2648909, 0, 2.0480251, 8.1112442, 0.97566891, 2.8674555;  
8.9254856, 1.2613275, 0.72119683, 3.0374279, -0.32588482, 0, 1.1993563, 0.63116002;  
2.6758013, 1.8996174, 1.7175947, 10.264549, 18.582258, 11.038349, 0.042666838, 0;  
1.6403483, 1.1003263, 3.8746152, 11.865837, 13.623778, 3.0804429, 2.7273426, 2.2249;  
0.44932228, 0, 0.45909494, 1.0929109, 4.1308932, 0.16889881, 5.6910644, 2.9264383;  
8.1910477, 14.858209, -0.071794562, 0, 2.7858784, 1.6050987, 2.841851, 2.2693861]  

diff: rows:8 cols:8 chans:1 min:-0.251977 max:17.4945  
[0, 0, 0, 0, 0, 0, 0, 0;  
1.5223728, 0.59031492, 0, 0, 0, 0, 0, 0;  
0, 0, 17.494459, 3.3915801, 0, 0, 0, 0;  
0, 0, 0, 0, 1.5637801, 1.9254711, 0, 0;  
0, 0, 0, 0, 0, 0, 11.408258, 2.6994755;  
0, 0, 0, 0, 0, 0, 0, 0;  
-0.2519767, 1.1847413, 0, 0, 0, 0, 0, 0;  
0, 0, 1.8295834, 2.0316832, 0, 0, 0, 0]

The difficult part for me is the nonzero entries in the diff matrix. I've accounted for the scaling FFTW does on the values and the padding needed to do an in-place fft on real data; what am I missing?

I find it surprising that the data could be off by a value of 17 (which is 66% of the max value), when there are so many zeros. Also, the data irregularities seem to form a diagonal pattern.

Answer 1

As you may have noticed when writting fftw_alloc_real(y * 2*(x/2 + 1)); fftw needs extra space in the x direction to store complex data. In your case, as x=8, it needs 2*(x/2+1)=10 reals.

http://www.fftw.org/doc/Real_002ddata-DFT-Array-Format.html#Real_002ddata-DFT-Array-Format

So...you should take care of this as you populate the input array or retreive values from it.

You way change

 input[yi*x + xi] = (double)yptr[xi];

for

 int xfft=2*(x/2 + 1);
 ...
 input[yi*xfft + xi] = (double)yptr[xi];

And

 yptr[xi] = (float)input[yi*x + xi];

for

 yptr[xi] = (float)input[yi*xfft + xi];

It should solve your problem since the non-nul points in your diff correspond to the extra padding.

Bye,