FFTW 实到实变换跨步阵列
FFTW real-to-real transform strided array
我有一个按行优先顺序存储的矩阵。我正在尝试使用 FFTW 计算子矩阵的 DCT,但我胡说八道。在接下来的段落中,我将描述问题和我的解决方案,我希望您能帮助理解为什么它不起作用。
给定一些 i
和 l
,我想计算一个子矩阵的 DCT,该子矩阵由 k mod l == i
的所有行 k
组成。例如,假设 l = 3
和 i = 2
。在下面的矩阵中,我要转换的子矩阵被标记为红色 (2 mod 3 = 2, 5 mod 3 = 2, 8 mod 3 = 2)。
源数组和目标数组具有相同的布局,转换后的矩阵应存储在目标数组中的相同位置。
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);
}
更新
我在 i=l=1
时针对简单情况进行了测试。即使在那种情况下,我也会胡说八道。我测试了一个 3x4 矩阵,它恰好是 DCT 基向量之一:
A(i,j) = cos((i + 0.5)*2*pi/3) * cos((j + 0.5)*3*pi/4)
我希望得到一个结果,其中除一个元素外所有元素都(接近)零。但是我得到的结果矩阵看起来像这样:
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
看起来很奇怪。
更新 2
我还测试了一个简单的矩阵,其中 (0,0) 元素为 1,其余为零。在这种情况下,也是 i=l=1
(子矩阵是整个矩阵)。这是结果:
2 2 2 0
1.73205 1.73205 1.73205 0
1 1 1 0
kind
应该是一个大小为 rank
的数组。参见例如http://www.fftw.org/doc/Guru-Real_002dto_002dreal-Transforms.html#Guru-Real_002dto_002dreal-Transforms.
我有一个按行优先顺序存储的矩阵。我正在尝试使用 FFTW 计算子矩阵的 DCT,但我胡说八道。在接下来的段落中,我将描述问题和我的解决方案,我希望您能帮助理解为什么它不起作用。
给定一些 i
和 l
,我想计算一个子矩阵的 DCT,该子矩阵由 k mod l == i
的所有行 k
组成。例如,假设 l = 3
和 i = 2
。在下面的矩阵中,我要转换的子矩阵被标记为红色 (2 mod 3 = 2, 5 mod 3 = 2, 8 mod 3 = 2)。
源数组和目标数组具有相同的布局,转换后的矩阵应存储在目标数组中的相同位置。
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);
}
更新
我在 i=l=1
时针对简单情况进行了测试。即使在那种情况下,我也会胡说八道。我测试了一个 3x4 矩阵,它恰好是 DCT 基向量之一:
A(i,j) = cos((i + 0.5)*2*pi/3) * cos((j + 0.5)*3*pi/4)
我希望得到一个结果,其中除一个元素外所有元素都(接近)零。但是我得到的结果矩阵看起来像这样:
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
看起来很奇怪。
更新 2
我还测试了一个简单的矩阵,其中 (0,0) 元素为 1,其余为零。在这种情况下,也是 i=l=1
(子矩阵是整个矩阵)。这是结果:
2 2 2 0
1.73205 1.73205 1.73205 0
1 1 1 0
kind
应该是一个大小为 rank
的数组。参见例如http://www.fftw.org/doc/Guru-Real_002dto_002dreal-Transforms.html#Guru-Real_002dto_002dreal-Transforms.