difference between matlab dct2 and FFTW lib(FFTW_REDFT10)
Abe Lau
I'm trying to rewrite an algorithm in C originally tested in
matlab using dct2 function. However, it seems to me that I don't get
the same results(apart from one is normalized and one is not) using
the dct in FFTW library when compared to what I get from dct2 in
matlab. I am wondering if anyone coulds give me a hint why the two
DCT will give different results.
For example:
for a 3x3 test array
[0.2 0.3 1
0 12 5
0.3 0.3 1]
dct2 (in matlab) gives
6.7000 -2.6536 -4.1719
-0.0408 -0.0500 -0.0289
-7.2832 2.4537 6.5500
however, in FFTW library using plan as follows:
fftw_plan_r2r_2d(3,3,fftw_disp_x,fftw_disp_x_w, FFTW_REDFT10,
FFTW_REDFT10, FFTW_ESTIMATE);
The documentation states that "FFTW_REDFT10" is THE DCT
I get:
126 28.4056 7.2
-74.4782 -20.4 -12.8172
21 -1.38564 18.6
The inverse DCT using their respective function -- idct2 in matlab
and FFTW_REDFT01 in FFTW give back the original array (just
un-normalized in FFTW).
Thanks so much if anyone could give me a hint what's the underlying
reason for the difference.
Cheers,
Abe
ste...@alum.mit.edu
> I'm trying to rewrite an algorithm in C originally tested in
> matlab using dct2 function. However, it seems to me that I don't get
> the same results(apart from one is normalized and one is not) using
> the dct in FFTW library when compared to what I get from dct2 in
> matlab. I am wondering if anyone coulds give me a hint why the two
> DCT will give different results.
>
> For example:
>
> for a 3x3 test array
> [0.2 0.3 1
> 0 12 5
> 0.3 0.3 1]
> I get:
> 126 28.4056 7.2
> -74.4782 -20.4 -12.8172
> 21 -1.38564 18.6
You're calling FFTW incorrectly, it looks like. (Your call to
fftw_plan_r2r_2d looks okay, so your bug is apparently elsewhere.) For
those inputs, FFTW gives:
80.4 -22.51666 -35.4
-0.3464102 -0.3 -0.1732051
-61.8 14.72243 39.3
This is the correct result. It differs from the Matlab function by
simple constant factors, as you can see if you simply compare the
definitions at:
http://www.mathworks.com/access/helpdesk/help/toolbox/images/dct2.html
http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
In particular, a MxN dct2 in Matlab and the corresponding REDFT10
transform in FFTW differ by:
1) an overall factor of 4 * sqrt(M/2) * sqrt(N/2)
2) an additional factor of sqrt(2) in the first row and first column
(and thus an additional factor of 2 in the upper-left output).
The former is simply an overall normalization. The latter is because
Matlab uses the unitary definition of the DCT-II, whereas FFTW uses the
definition that makes it equivalent to a real-even FFT of data shifted
by half a sample. None of these choices have a significant effect on
the application of the transform.
Regards,
Steven G. Johnson
Abe Lau
Thanks for telling me I've got something wrong with my FFTW code.
Yes I'm sorry I made a mistakes in the input when I did the test
array in debugging. I got the same results (apart from
normalization) after correcting my mistakes.
I didn't notice the
"2) an additional factor of sqrt(2) in the first row and first column
(and thus an additional factor of 2 in the upper-left output)."
as you have pointed out. That's the reason why I thought I got the
wrong results when I'm processing my image which is of size 512x512.
have been trying to work that out for a month but didn't notice this
additional factor in their definition.
Thanks once again. I really appreciate that.
Regards,
Abe Lau
Mihai
could you please tell me if you have used the fftw in a .c code?
I'm having trouble compiling a C code in Windows, actually the
trouble is linking the libfftw3-3.lib libraries.
could you please have a look at my post:
<http://newsreader.mathworks.com/WebX/.ef3877c?50@@>
Maybe you have an idea regarding my problems.
Thank you
mihai
Abe Lau
just got your email. seems that you solved the problem linking all
those library, but come to the point where the FFTW dct2 output is
different than what you expected (even after renormalization as
Steven has mentioned).
I do my C++ code on GNU/Linux platform using FFTW library from
fftw.org(linked statically). not really sure how to do that in MS
Windows though.
<http://www.fftw.org/fftw3_doc/Tutorial.html#Tutorial>
this tutorial has some samples of fftw use. Let me know if you wanna
have a look at the code I use in my C++ project.
Regards,
Abe Lau
Mihai
for me everything is working fine, I have a Simulink block which does
my own fft transform, and the results are the same as the matlab fft.
With the ifft I have some problems, using the same plan with
FFTW_BACKWARD gives me some strange results.
My real problem is MS Windows. I need to get the block working in
Windows too.
I could have a look at your C code because I could always learn
something. Here is my mdlOutputs, which uses the fftw:
static void mdlOutputs(SimStruct *S, int_T tid)
{
real_T *w_vect_r, *w_vect_i;
int_T i, j;
InputRealPtrsType uPtrs = ssGetInputPortRealSignalPtrs(S,0);
int_T width = ssGetOutputPortWidth(S,0);
const boolean_T cplx = (ssGetInputPortComplexSignal(S, 0) ==
COMPLEX_YES);
real_T *y = ssGetOutputPortRealSignal(S,0);
/* here begins the typical fftw code */
fftw_plan p;
fftw_complex *in;
fftw_complex *out;
/*memory allocation*/
out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*width);
in=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*width);
/*for filling the input array*/
for(i=0;i<width;i++){
in[i][0]=uPtrs[i][0];
in[i][1]=uPtrs[i][1];
}
/*plan creation and exection*/
p=fftw_plan_dft_1d(width,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
fftw_execute(p);
/*filling out the output array, which is the output of our
S-Function */
for(i=0;i<width;i++)
{
*y++ =out[i][0];
if (cplx)
{
*y++ = out[i][1]; /* typical way of adressing complex
numbers in matlab C MEX files */
}
}
/*resources deallocation*/
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
}
Another stupid thing is that after upgrading my ubuntu linux from
5.04 to 6.06 and Matlab from 7.0.1 to 7.1.0 I can't compile this
functions anymore. I posted the problem here:
Mihai, "strange things about mex compiling" #, 14 Jun 2006 2:17 pm </WebX?14@@.ef39656>
So, if you could give a tip on the inverse transformation or just
show me your code it would be nice. It's so good when you get some
answers in the forum.
Thank you
Mihai
Aishwarya Selvaraj
> fftw_plan_r2r_2d looks okay, so your bug is apparently elsewhere.) For
> those inputs, FFTW gives:
>
> 80.4 -22.51666 -35.4
> -0.3464102 -0.3 -0.1732051
> -61.8 14.72243 39.3
>
> This is the correct result. It differs from the Matlab function by
> simple constant factors, as you can see if you simply compare the
> definitions at:
>
> In particular, a MxN dct2 in Matlab and the corresponding REDFT10
> transform in FFTW differ by:
>
> 1) an overall factor of 4 * sqrt(M/2) * sqrt(N/2)
>
> 2) an additional factor of sqrt(2) in the first row and first column
> (and thus an additional factor of 2 in the upper-left output).
>
> The former is simply an overall normalization. The latter is because
> Matlab uses the unitary definition of the DCT-II, whereas FFTW uses the
> definition that makes it equivalent to a real-even FFT of data shifted
> by half a sample. None of these choices have a significant effect on
> the application of the transform.
I noticed the difference in the DCT result when compared to matlab dct2() function.
In you r post you mentioned the reason for the difference in value .
When I tried to find dct of a matrix of size mxn where m==n ;
The logic was right .
But when m is different from n . the logic was wrong .
Could you please throw some knowledge here .
Aishwarya Selvaraj
> fftw_plan_r2r_2d looks okay, so your bug is apparently elsewhere.) For
> those inputs, FFTW gives:
>
> 80.4 -22.51666 -35.4
> -0.3464102 -0.3 -0.1732051
> -61.8 14.72243 39.3
>
> This is the correct result. It differs from the Matlab function by
> simple constant factors, as you can see if you simply compare the
> definitions at:
>
> In particular, a MxN dct2 in Matlab and the corresponding REDFT10
> transform in FFTW differ by:
>
> 1) an overall factor of 4 * sqrt(M/2) * sqrt(N/2)
>
> 2) an additional factor of sqrt(2) in the first row and first column
> (and thus an additional factor of 2 in the upper-left output).
>
> The former is simply an overall normalization. The latter is because
> Matlab uses the unitary definition of the DCT-II, whereas FFTW uses the
> definition that makes it equivalent to a real-even FFT of data shifted
> by half a sample. None of these choices have a significant effect on
> the application of the transform.