有没有办法在另一个维度上计算 2D FFT 的 1D FFT 而无需使用 intel mkl 进行转置？

Question

我想使用 mkl 计算存储为一维数组的二维数组的一维 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。

同样，英特尔数学内核库引入了data layout configuration parameters，例如DFTI_INPUT_STRIDES and 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_OUTPUT_STRIDES 对于就地转换 (DFTI_PLACEMENT=DFTI_INPLACE) 会被忽略。

Answer 2

据我所知，英特尔 MKL 不提供对数据元素之间具有跨度的数据执行 FFT 的能力。

但是，FFTW 可以。每4.4.1 Advanced Complex DFTs of the FFTW documentation：

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 基本接口（参见复杂 DFTs) 提供了专门用于等级 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（即相同的物理和逻辑维度，如 基本界面。）

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

对于就地转换，输入和输出步幅和距离 参数应该相同；否则，规划者可能会返回 NULL.

该页面方便地提供了一个（有些令人困惑的）对二维数组的列执行一维 FFT 的示例：

将二维数组的每一列转换为 10 行 3 列：

   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;

有关更多示例，另请参阅How do I use fftw_plan_many_dft on a transposed array of data?。

Answer 3

您无法沿数据的更高维度进行一维 FFT。您需要先进行转置，以使 FT 维度成为数据在 RAM 中连续的维度。

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