3

I am testing the library just using a discrete cosine transform.

#include <fftw3.h>
void dump_vector(int n, double* vec) {
    for(int i = 0; i < n; i++)
        printf("%f ", vec[i]);
    printf("\n");
}
int main()
{
    double a[] = {0.5, 0.6, 0.7, 0.8};
    double b[] = {0, 0, 0, 0};
    printf("Original vector\n");
    dump_vector(4, a);
    fftw_plan plan = fftw_plan_r2r_1d(4, a, a, FFTW_REDFT10, FFTW_ESTIMATE);
    fftw_execute(plan);
    printf("DCT\n");
    dump_vector(4, a);
    fftw_plan plani = fftw_plan_r2r_1d(4, a, a, FFTW_REDFT10, FFTW_ESTIMATE);
    fftw_execute(plani);
    printf("IDCT\n");
    dump_vector(4, a);
    return 0;
}

I hope to obtain the same a, or maybe an approximation, but my ouput is the following:

Original vector
0.500000 0.600000 0.700000 0.800000 
DCT
5.200000 -0.630864 0.000000 -0.044834 
IDCT
9.048603 9.208347 8.182682 5.179908
Chris Nolet
  • 8,714
  • 7
  • 67
  • 92
FacundoGFlores
  • 7,858
  • 12
  • 64
  • 94

1 Answers1

2

See the documentation of fftw about real to real transforms

FFTW_REDFT10 computes an REDFT10 transform, i.e. a DCT-II (sometimes called “the” DCT). (Logical N=2*n, inverse is FFTW_REDFT01.)

Hence, the flag FFTW_REDFT01 must be used for the inverse transform instead of FFTW_REDFT10.

Moreover, FFTW does not rescale the output of the transform. Consequently, the ouput must be divised by either n or n*n, where n is the length of the vector. (I will test it in a couple of minutes...)

EDIT: the scaling factor is neither n nor n*n it's 2*n...

francis
  • 9,525
  • 2
  • 25
  • 41