Forums

Matlab and VC++ values differ when determining FFTusing FFTW

Started by chiraag August 27, 2009
Hi there,

I have recently started using FFTW to determine fft in C++, so I am really
not good at it. I am currently trying to just evaluate my fft values in
both matlab and C++. I have used FFTW to determine the fft in C++ and fft2
in matlab for a 4*4 matrix [int values]. My values were same in both Matlab
and C++. My real requirement is to use fft on complex matrices. But my
values no longer match in both. There are variations in the values when I
do this. I have used the same fft plan for both int and complex matrices.

I have pasted the portion of my code here. Please take a look.

#include <iostream>
#include "stdafx.h"

#include "cvaux.h"
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>

#include <math.h>
#include "complex.h"

#include "fftw++.h"

double  x,y,r;
Complex  UR[960][1280];
Complex  H_r[960][1280];
Complex B[960][1280];
Complex b[960][1280];

void main()

{int Nx=1280,Ny=960;
 double lamda=532;
 double pinhole_ccd_D=63;
 double pinhole_ccd_D_rec=pinhole_ccd_D*100;
IplImage *img=0,*img2=0;
img=cvLoadImage("f1.bmp",0);//Please use any image
int height=img->height,
width=img->width,
step=img->widthStep;
CvScalar H;

for (m=0;m<Nx;m++)
{ for(n=0;n<Ny;n++)
{x = (1 + m- Nx / 2);
y = (1 + n- Ny / 2);
r= pow((x * dx),2) + pow((y * dy), 2);
UR[n][m]=Complex((cos(k*sqrt(r+pinDrec))/sqrt(r+pinDrec)),(-sin(k*sqrt(r+pinDrec))/sqrt(r+pinDrec)));
H=cvGet2D(img,n,m);
B[n][m]= Complex((double) H.val[0],0);
H_r[n][m]=Complex((B[n][m].Re()*UR[n][m].Re()),(B[n][m].Re()*UR[n][m].Im()));
}}
   uchar *img_data;
   Complex *img2_data;
   fftw_complex   *data_in;
   fftw_complex   *fft;        
   fftw_plan      plan_f;
//initialize arrays for fftw operations 
data_in = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width *
height );
 fft     = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width *
height );

// create plans 
plan_f = fftw_plan_dft_1d( width * height, data_in, fft,  FFTW_FORWARD, 
FFTW_ESTIMATE );
// load img1's data to fftw input 
    for(int m = 0, k1 = 0 ; m < width ; m++ ) 
	{for( n = 0 ; n < height ; n++ ) 
		{
            data_in[k1][0] = H_r[n][m].Re();
            data_in[k1][1] = H_r[n][m].Im();
            k1++;
        }
    }
    
    // perform FFT
    fftw_execute( plan_f );

   for( m = 0,k1=0; m < width ; m++ ) 
{ for( n = 0 ; n < height ; n++ )
{ b[m][n]=Complex(fft[k1][0],fft[k1][1]);
k1++;}}

   	
Could some one please guide me in this regard??
Thank you.
				



On 27 Aug, 16:55, "chiraag" <sinara...@gmail.com> wrote:
> Hi there, > > I have recently started using FFTW to determine fft in C++, so I am really > not good at it. I am currently trying to just evaluate my fft values in > both matlab and C++. I have used FFTW to determine the fft in C++ and fft2 > in matlab for a 4*4 matrix [int values]. My values were same in both Matlab > and C++. My real requirement is to use fft on complex matrices. But my > values no longer match in both. There are variations in the values when I > do this. I have used the same fft plan for both int and complex matrices.
Different implementations of the same algorithm will always (except for the most trivial cases) produce different results. One would expect variations in the results, on the relative order of eps = 1e-6 for single precision floating point numbers, or eps = 1e-15 for double precision floating point numbers. 'On the order of eps' means that one might see relative variations up to 10*eps or maybe 50*eps. If the variations are larger, there might be implementation errors in one of the routines. Rune
chiraag <sinara009@gmail.com> wrote:

> But my values no longer match in both. There are variations > in the values when I do this. I have used the same fft plan > for both int and complex matrices.
Can you quantify for us how the two results differ? i.e. are they different by just round-off errors; are they the wrong sign or scaled differently; or is one or the other outputs flat-out wrong? Steve
>chiraag <sinara009@gmail.com> wrote: > >> But my values no longer match in both. There are variations >> in the values when I do this. I have used the same fft plan >> for both int and complex matrices. > >Can you quantify for us how the two results differ? i.e. >are they different by just round-off errors; are they >the wrong sign or scaled differently; or is one or the >other outputs flat-out wrong? > >Steve
On the assumption that they're more than just off by a slight amount, there are a few other things: 1) Loading images can be tricky (RGB vs BGR, strides, row-major vs column-major, etc). Have you verified that your input array is the same as in MATLAB? For an image of that size, your eyes will likely get tired of comparing... 2) Question 3.18 http://www.fftw.org/faq/section3.html#vbetalia might be appropriate. Are you trying to compile FFTW with MSVC? I have successfully *used* FFTW with MSVC, but I used a precompiled binary of FFTW for that very reason (I usually prefer to compile stuff from source, but I made an exception for FFTW). In the same vein, have you tried compiling your program on *nix instead? OpenCV is cross-platform... 3) Your code is very difficult to read. It lacks proper indentation, and it looks like certain variables (e.g., "step") are unused. If you clean it up some, you might get more people to review it. You also might have a better chance of understanding it if you have to look at it again in a year or two. 4) If all else fails, you might just try a canned set of maybe a half-dozen values through a 1d transform, and see if that agrees.
I have already tried the FFT plan[1d transform] for a 4*4 matrix and my
values do agree with each other in both Matlab and VC++.

The matrix used for the same was 
a[4][4]= {0,2,4,6,
          1,3,5,7,
          4,6,8,10,
          9,11,13,15};
The fft values that i obtained in both matlab and c++ is 
fft={104+0i,-16+16i,-16+0i,-16-16i,
-16+32i,0+0i,0+0i,0+0i,
-24+0i,0+0i,0+0i,0+0i,
0+0i,0+0i,0+0i,0+0i};

I am using the same fft plan for the complx matrix. I have already
verified that the image values. They are the same. Till the point before i
execute the fft plan, the values in both matlab and c++ are the same with
very small differences. Once I execute my plan, the changes occur. I have
take the default precision in FFTW.

I will show some of the values after fft implementation.

FFTW-C++                          
-1501149.47126+398668.5919i
MATLAB VALUE
1.50114859e+06+3.98671903e+05i

the firt value alone is right, there after the values seem to differ.

FFTW C++ VALUE
1106677.08963712-55975.9601i
MATLAB VALUE
1.106787211e+06-5.3278153e+04
FFTW C++ VALUE
-1115382.849286-1096920.917146i
MATLAB VALUE
-1.110085929120e+06-1.1025766182e+06i
FFTW C++ VALUE
-12266619.02617079+849319.819840i
MATLAB VALUE
-1.272957528e+06+8.4017428300i
FFTW C++ VALUE
-83882.284090035+1079329.51661278i
MATLAB VALUE
-9.421581101e+04+1.0785670834e+06i

These are just the first five values in the matrix. As you can see there
has already been a drift in values.I am using the same plan however with no
considerations of single or double precision. Just used the default. I
tried using long double precision by adding l to all the fftw_ calls. But
when i do this error props up. Could you please tell me where I could be
going wrong??

Thank you.

>
On 28 Aug, 05:41, "chiraag" <sinara...@gmail.com> wrote:

> I will show some of the values after fft implementation. > > FFTW-C++ =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 > -1501149.47126+398668.5919i > MATLAB VALUE > 1.50114859e+06+3.98671903e+05i > > the firt value alone is right, there after the values seem to differ.
The computed results seem to deviate somwhere around the 6th or 7th significant digit, which happens to be the numerical accuracy of the single-precision floating point format. The very first thing I would do, is to verify that both matlab and FFTW use the same floating-point format for their computations. After that, I would check the documentation for FFTW and verifyc that the libraries are used as intended. Some of the functions in FFTW require the input data to be formatted in ways that might not be obvious to a first-time user of the library. Rune
These are just the first five values in the matrix. As you can see there
>has already been a drift in values.I am using the same plan however with
no
>considerations of single or double precision. Just used the default. I >tried using long double precision by adding l to all the fftw_ calls.
But
>when i do this error props up. Could you please tell me where I could be >going wrong??
Another thing you might try, right after running the transform in C, is to run the inverse transform in C and see if you get your original data back (fabs(old-new) < thresh_relative*eps*fmax(fabs(old),fabs(new)) or something similar).
>1) Loading images can be tricky (RGB vs BGR, strides, row-major vs
column-major, etc). One more I forgot to mention: top-down vs bottom-up DIBs (BMP's). Annoyingly, bottom-up is the default (legacy OS/2 stuff, believe it or not), meaning the bottom row is the first one saved in the file...usually. Instead of hard-coding a matrix, load a 4x4 grayscale BMP using the same OpenCV logic. You might also post the MATLAB code, if it's short. At what point in the code did you check that the images were equal? It looks like you're concatenating each row of the image. Unless you're comparing raw numbers (e.g., fprintf a CSV file in C and use csvread or dlmread in MATLAB) right before the transform, I wouldn't necessarily conclude the problem is with the use of FFTW. That said, if the precision you're using is the same as MATLAB, the fact that MATLAB uses some version of FFTW should be of some help.