c语言 - 有没有办法在另一个维度上计算二维 FFT 的 1D FFT，而无需使用英特尔 mkl 进行转置?

Question

我想使用 mkl 来计算存储为 1D 数组的 2D 数组的一维 FFT。 例如

for (int j=0; j<NJ; j++) //rows
{
for (int i=0; i<NI; i++) //columns
{
Pre_2D_array[i+j*NI].x=1.0;
Pre_2D_array[i+j*NI].y=2.0;
}
}

我想在行维度上计算Pre_2D_array的一维 FFT。我能想到的唯一方法是重塑阵列并像这样做 FFT，

for (int i=0; i<NI; i++) //columns
{
for (int j=0; j<NJ; j++) //rows
{
2D_array[j+i*NJ]=Pre_2D_array[i+j*NI];
}
}
DFTI_DESCRIPTOR_HANDLE desc_x = 0;
DftiCreateDescriptor(&desc_x, DFTI_PREC, DFTI_COMPLEX, 1, NJ);
DftiSetValue(desc_x, DFTI_NUMBER_OF_TRANSFORMS, NI);
DftiSetValue(desc_x, DFTI_INPUT_DISTANCE,  NJ);
DftiCommitDescriptor(desc_x);
DftiComputeForward(desc_x, 2D_array);

这能得到正确的答案。但是当数组很大时，对原始数组进行转置(重塑)会浪费太多时间。 有没有办法在不重塑阵列的情况下进行FFT？或者有什么快速的方法可以尽快重塑阵列？

CPUINFO 是：

processor   : 0
vendor_id   : GenuineIntel
cpu family  : 6
model       : 79
model name  : Intel(R) Xeon(R) CPU E5-2648L v4 @ 1.80GHz
stepping    : 1
microcode   : 0xb000022
cpu MHz     : 1795.882
cache size  : 35840 KB
physical id : 0
siblings    : 14
core id     : 0
cpu cores   : 14
apicid      : 0
initial apicid  : 0
fpu     : yes
fpu_exception   : yes
cpuid level : 20
wp      : yes
flags       : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc aperfmperf eagerfpu pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 fma cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch arat epb xsaveopt pln pts dtherm tpr_shadow vnmi flexpriority ept vpid fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm rdseed adx smap
bogomips    : 3591.76
clflush size    : 64
cache_alignment : 64
address sizes   : 46 bits physical, 48 bits virtual
power management:

Answer 1

FFTW 库在 fftw_plan_many_dft() 等函数中引入了istride或ostride参数，以避免转置数组。该页面上的最后一个示例是第二个维度上的 DFT。

同样，英特尔数学核心函数库引入了数据布局配置参数，如DFTI_INPUT_STRIDES和DFTI_OUTPUT_STRIDES或DFTI_NUMBER_OF_TRANSFORMS。

第二维度上的DFT可能看起来像(我没有测试过)：

DftiCreateDescriptor(&desc_x, DFTI_PREC, DFTI_COMPLEX, 1, NJ);
DftiSetValue(desc_x, DFTI_NUMBER_OF_TRANSFORMS, NI);
DftiSetValue(desc_x, DFTI_INPUT_STRIDES, &NI);
DftiSetValue(desc_x, DFTI_OUTPUT_STRIDES, &NI);
DftiSetValue(desc_x, DFTI_INPUT_DISTANCE,  1);
DftiSetValue(desc_x, DFTI_OUTPUT_DISTANCE,  1);
DftiCommitDescriptor(desc_x);

对于就地转换 (DFTI_PLACEMENT=DFTI_INPLACE)，将忽略DFTI_OUTPUT_STRIDES。

Answer 2

据我所知，英特尔 MKL 无法在数据元素之间跨步执行数据中的 FFT。

但是，FFTW确实如此。 根据 FFTW 文档的4.4.1 高级复杂 DFT：

fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany,
fftw_complex *in, const int *inembed,
int istride, int idist,
fftw_complex *out, const int *onembed,
int ostride, int odist,
int sign, unsigned flags);

此例程计划多个多维复杂 DFT，并且 扩展fftw_plan_dft例程(请参阅复杂 DFT)以计算 有多少个变换，每个变换都有等级等级和大小n。另外 转换数据不必是连续的，但可以布置在 记忆与任意的步伐。考虑到这些可能性，fftw_plan_many_dft添加新参数howmany、{i,o}nembed、{i,o}stride，{i,o}dist。FFTW 基本接口(请参阅复杂 DFT) 提供专门用于等级 1、2 和 3 的例程，但 高级界面仅处理常规排名情况。

howmany是要计算的转换的(非负)数。这 生成的计划计算howmany转换，其中 第 k 个变换位于位置in+k*idist(在 C 指针算术中)， 它的输出位于位置out+k*odist.在此获得的计划 方式通常比多次调用 FFTW 更快 单个转换。基本fftw_plan_dft接口对应 以howmany=1(在这种情况下，将忽略 dist 参数)。

每个howmany变换都有秩和大小n，如 基本界面。此外，高级界面允许输入 并将每个转换的输出数组作为 较大的等级数组，由inembed和onembed描述 参数，分别。{i,o}nembed必须是长度的数组rank，并且n元素上应小于或等于{i,o}nembed.传递nembed参数的NULL是等效的 传递n(即相同的物理和逻辑维度，如 基本界面。

步幅参数指示输入的j-th元素或 输出数组分别位于j*istride或j*ostride处。 (对于多维数组，j是普通的行主索引。 当与howmany循环中的k-th转换结合使用时，从 上面，这意味着 (j，k) - 第 1 个元素位于j*stride+k*dist处。 (基本fftw_plan_dft界面对应于步幅为 1。

对于就地变换，输入和输出步幅和距离 参数应相同;否则，计划人员可能会返回NULL.

该页面方便地提供了一个(有点令人困惑的)在二维数组的列上执行一维FFT的示例：

使用 10 行和 3 列转换 2d 数组的每一列：

int rank = 1; /* not 2: we are computing 1d transforms */
int n[] = {10}; /* 1d transforms of length 10 */
int howmany = 3;
int idist = odist = 1;
int istride = ostride = 3; /* distance between two elements in 
the same column */
int *inembed = n, *onembed = n;

有关更多示例，请参阅如何在转置的数据数组上使用fftw_plan_many_dft？

Answer 3

您不能沿数据的更高维度进行 1D FFT。您需要先进行转置，以便 FT 维度是数据在 RAM 中连续的维度。

但是，它并不像您想象的那么糟糕。在多核机器上，您可以轻松设置一些线程，其唯一工作是预先/后期排列FT数据。