Есть ли способ вычислить 1D БПФ из 2D БПФ в другом измерении без транспонирования с использованием Intel MKL?

Question

Я хочу использовать mkl для вычисления 1D БПФ двумерного массива, хранящегося как одномерный массив.например,

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;
   }
}

Я хочу вычислить 1D БПФ для Pre_2D_array в измерении строки.Единственный способ, которым я могу придумать, - это изменить форму массива и выполнить БПФ, как это,

   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);

Хотя это может дать правильный ответ.Но выполнение транспонирования (изменения формы) исходного массива тратит слишком много времени, когда массив большой.Есть ли способ сделать БПФ без изменения формы массива?Или любой быстрый способ изменить массив как можно быстрее?

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

Насколько я могу судить, Intel MKL не предоставляет возможности выполнять БПФ в данных с шагом между элементами данных.

Однако FFTW делает. Согласно 4.4.1 Расширенные сложные DFT документации FFTW :

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);

Эта процедура планирует несколько многомерных сложных ДПФ, и это расширяет подпрограмму fftw_plan_dft (см. сложные DFT) для вычисления Howmany преобразований, каждый из которых имеет ранг ранг и размер 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. Передача NULL для параметра nembed эквивалентна на прохождение n (т.е. те же физические и логические измерения, что и в базовый интерфейс.)

Параметры шага указывают, что элемент j-th ввода или Выходные массивы расположены на j*istride или j*ostride соответственно. (Для многомерного массива j - это обычный основной индекс строки.) В сочетании с преобразованием k-th в цикле howmany, из выше, это означает, что (j, k) -й элемент находится в j*stride+k*dist. (Базовый интерфейс fftw_plan_dft соответствует шагу 1.)

Для преобразований на месте, входной и выходной шаг и расстояние параметры должны быть одинаковыми; в противном случае планировщик может вернуться NULL.

На странице удобно представлен (несколько запутанный) пример выполнения 1D БПФ над столбцами двумерного массива:

Преобразование каждого столбца двумерного массива с 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;

Также см. Как использовать fftw_plan_many_dft для транспонированного массива данных? для дополнительных примеров.

Answer 2

Библиотека FFTW ввела аргументы istride или ostride в таких функциях, как fftw_plan_many_dft () , чтобы избежать транспонирования массива. Последний пример на этой странице - это ДПФ во втором измерении.

Аналогично, математическая библиотека Intel Intel вводит параметры конфигурации макета данных , такие как DFTI_INPUT_STRIDES и DFTI_OUTPUT_STRIDES или DFTI_NUMBER_OF_TRANSFORMS.

ДПФ над вторым измерением может выглядеть (я не проверял):

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 3

Вы не можете выполнять одномерные БПФ для более высоких измерений ваших данных.Сначала необходимо выполнить транспонирование, чтобы измерение FT было тем, в котором данные непрерывно располагаются в оперативной памяти.

Однако, это не так плохо, как вы думаете.На многоядерном компьютере вы можете легко настроить некоторые потоки, единственной задачей которых является предварительная / последующая организация данных FT.