Forums

fftw and dct/idct

Started by markus_h March 26, 2008
Hi,
I'm trying to convert a piece of Matlab code to C using the fftw library,
but I'm having some troubles with how to handle the Matlab function
"idct". I haven't found much documentation about this but I'm assuming
it's a type-II dct so I've been using the "FFTW_REDFT01" but I don't get
the same results as when I run the script in Matlab. Using google I have
found some vague hints that I might have to rescale the data (other than
with the size of the signal) with the square root of 2. I've been playing
around with some test data for a few hours, rescaling and using different
types of dct's, but I still can't make any sense of why I get different
results. 

The code: 
"
double temp[3];
	
fftw_plan idct =	fftw_plan_r2r_1d(3, temp, temp, FFTW_REDFT01,
FFTW_MEASURE);

temp[0]=2;
temp[1]=1;
temp[2]=0;

fftw_execute(idct);

for(int i = 0; i < 3; i++)
   temp[i] /=3;
"
results in the output 1.24402, 0.66667 and 0.089316 whereas in Matlab for
the array temp = [2; 1; 0] using idct(temp) I get 1.68161, 1.15470 and
0.44759

So, anyone have any idea why I get different results here or can shed some
light as to how the function idct works in Matlab compared to the fftw
equivalent?
Thanks,
Markus H


On Mar 26, 3:40 pm, "markus_h" <m_henningsso...@hotmail.com> wrote:
> I'm trying to convert a piece of Matlab code to C using the fftw library, > but I'm having some troubles with how to handle the Matlab function > "idct". I haven't found much documentation about this but I'm assuming > it's a type-II dct so I've been using the "FFTW_REDFT01" but I don't get > the same results as when I run the script in Matlab.
Hi Markus, there are two sources to check which immediately reveal the answer: the Matlab manual and the FFTW manual. In particular, if you look in the two manuals, they give the exact mathematical definitions of the transforms as computed by Matlab and FFTW: http://www.mathworks.com/access/helpdesk/help/toolbox/dspblks/ref/idct.html http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html If you compare the formulas, you see that they differ in two ways. For Matlab, most of the inputs are multiplied by sqrt(2/N) where N is the length of the transform, whereas in FFTW the same inputs are multipled by 2. In Matlab, the first input is multiplied by sqrt(1/N) wheras in FFTW the same input is unscaled. So, if you have an array x in Matlab and want to compute FFTW's REDFT01 transform, you would do: y = x; y(1) = x(1) / sqrt(2); idct(y) * sqrt(2*length(x)) Conversely, if you want to compute Matlab's idct in FFTW with its REDFT01 routine, you would multiply input[0] by sqrt(1.0/N) and the rest of the inputs by sqrt(0.5/N), where N is the length. Regards, Steven G. Johnson PS. Matlab's normalization is chosen to make the idct unitary. FFTW normalizes the transform, instead, to emphasize the connection to a DFT of real data of a certain symmetry.
>Hi Markus, there are two sources to check which immediately reveal the >answer: the Matlab manual and the FFTW manual. In particular, if you >look in the two manuals, they give the exact mathematical definitions >of the transforms as computed by Matlab and FFTW: > >http://www.mathworks.com/access/helpdesk/help/toolbox/dspblks/ref/idct.html >http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html > >If you compare the formulas, you see that they differ in two ways. >For Matlab, most of the inputs are multiplied by sqrt(2/N) where N is >the length of the transform, whereas in FFTW the same inputs are >multipled by 2. In Matlab, the first input is multiplied by sqrt(1/N) >wheras in FFTW the same input is unscaled. > >So, if you have an array x in Matlab and want to compute FFTW's >REDFT01 transform, you would do: > > y = x; y(1) = x(1) / sqrt(2); idct(y) * >sqrt(2*length(x)) > >Conversely, if you want to compute Matlab's idct in FFTW with its >REDFT01 routine, you would multiply input[0] by sqrt(1.0/N) and the >rest of the inputs by sqrt(0.5/N), where N is the length. > >Regards, >Steven G. Johnson > >PS. Matlab's normalization is chosen to make the idct unitary. FFTW >normalizes the transform, instead, to emphasize the connection to a >DFT of real data of a certain symmetry. >
Thanks Steven! That cleared up my idct problem. However now I've ran into similar problems with the idst. I've looked at the Matlab and fftw sites here http://www.mathworks.com/access/helpdesk/help/toolbox/pde/index.html?/access/helpdesk/help/toolbox/pde/ug/dst.html and here http://www.fftw.org/doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html#g_t1d-Real_002dodd-DFTs-_0028DSTs_0029 but I'm fairly new to these transforms and don't know how I would go about seeing the difference. And sorry to bother you again but google isn't being very good to me with this. Thanks again, Markus
On Mar 27, 6:27 am, "markus_h" <m_henningsso...@hotmail.com> wrote:
> Thanks Steven! That cleared up my idct problem. > However now I've ran into similar problems with the idst. I've looked at > the Matlab and fftw sites herehttp://www.mathworks.com/access/helpdesk/help/toolbox/pde/index.html?... > and herehttp://www.fftw.org/doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html#g_t1... > > but I'm fairly new to these transforms and don't know how I would go about > seeing the difference. And sorry to bother you again but google isn't being > very good to me with this.
By inspection of the formulas, Matlab's idst normalization differs from the corresponding transform in FFTW (RODFT00) by a factor of 1/(N +1); the dst normalization differs by a factor of 2. (Interestingly, Matlab didn't use the unitary normalization for the dst/idst, unlike for their dct/idct, and they defined it to be the DST- I rather than the DCT-II for the dct function.) Regards, Steven G. Johnson