I have a matrix stored in row-major order. I am trying to compute the DCT of a ub-matrix using FFTW, and I get nonsense. In the following paragraphs I will describe the problem and my solution, and I would like your help in understanding why it does not work.
Given some i and l, I would like to compute the DCT of a sub-matrix which consists of all rows k for which k mod l == i. For example, assume that l = 3 and i = 2. In the following matrix, the sub-matrix which I wish to transform is marked in red (2 mod 3 = 2, 5 mod 3 = 2, 8 mod 3 = 2).
The source and destination arrays have the same layout, and the transformed matrix should be stored in the same location in the destination array.
void transform(double* src, double* dest, size_t rows, size_t cols, size_t l, size_t i)
{
int rank = 2;
fftw_iodim64 dims[] = {
{ rows / l, l, l },
{ cols, rows, rows } };
fftw_r2r_kind kind = FFTW_REDFT10;
fftw_plan plan = fftw_plan_guru64_r2r(rank, dims, 0, NULL, src + l, dest + l, &kind, FFTW_ESTIMATE | FFTW_UNALIGNED | FFTW_PRESERVE_INPUT);
fftw_execute(plan);
fftw_destroy_plan(pla);
}
Update
I tested it for the simple case when i=l=1. Even in that case I get nonsense. I tested with a 3x4 matrix which is exactly one of the DCT basis vectors:
A(i,j) = cos((i + 0.5)*2*pi/3) * cos((j + 0.5)*3*pi/4)
I expected to get a result where all elements are (close to) zero except for one. But I get a resulting matrix that looks like this:
0 -2.22045e-016 1.33227e-015 2.22045e-016
2.22045e-016 -2.77556e-016 9.99201e-016 5.55112e-017
-8.88178e-016 -1.62359 7.83938 1.62359
Seems quite weird.
Update 2
I also tested with a simple matrix where the (0,0) element is 1 and the rest are zero. In this case, too, i=l=1 (the sub-matrix is the whole matrix). Here is the result:
2 2 2 0
1.73205 1.73205 1.73205 0
1 1 1 0
