كيف يمكنني استخدام fftw_plan_many_dft على مجموعة من البيانات المنقولة؟

Question

لدي مجموعة ثنائية الأبعاد من البيانات المخزنة بتنسيق عمود رئيسي (نمط فورتران)، وأرغب في الحصول على التحويل السريع لكل منها صف.أرغب في تجنب تبديل المصفوفة (فهي ليست مربعة).على سبيل المثال، مجموعتي

fftw_complex* data = new fftw_complex[21*256];

يحتوي على إدخالات [r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255].

استطيع ان استخدم fftw_plan_many_dft لوضع خطة لحل كل من 21 FFTs الموجودة في data المصفوفة إذا كانت كذلك صف-الرائد، على سبيل المثال. [r0_val0, r0_val1,..., r0_val255, r1_val0,...,r20_val255]:

int main() {
  int N = 256;
  int howmany = 21;
  fftw_complex* data = new fftw_complex[N*howmany];
  fftw_plan p;

  // this plan is OK
  p = fftw_plan_many_dft(1,&N,howmany,data,NULL,1,N,data,NULL,1,N,FFTW_FORWARD,FFTW_MEASURE);
  // do stuff...

  return 0;
}

حسب الوثائق (القسم 4.4.1 من دليل 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);

ويجب أن أكون قادرًا على استخدام stride و dist المعلمات لتعيين الفهرسة.مما أستطيع أن أفهمه من الوثائق، تتم فهرسة الإدخالات في المصفوفة المراد تحويلها على أنها in + j*istride + k*idist أين j=0..n-1 و k=0..howmany-1.(مصفوفاتي هي 1D وهناك howmany منهم).ومع ذلك، ينتج التعليمة البرمجية التالية في seg.عيب (يحرر: طول الخطوة هو خطأ, ، انظر التحديث أدناه):

int main() {
  int N = 256;
  int howmany = 21;
  fftw_complex* data = new fftw_complex[N*howmany];
  fftw_plan p;

  // this call results in a seg. fault.
  p = fftw_plan_many_dft(1,&N,howmany,data,NULL,N,1,data,NULL,N,1,FFTW_FORWARD,FFTW_MEASURE);

  return 0;
}

تحديث:

لقد ارتكبت خطأ في اختيار طول الخطوة.المكالمة الصحيحة هي (طول الخطوة الصحيح هو howmany, ، لا N):

int main() {
  int N = 256;
  int howmany = 21;
  fftw_complex* data = new fftw_complex[N*howmany];
  fftw_plan p;

  // OK
  p = fftw_plan_many_dft(1,&N,howmany,data,NULL,howmany,1,data,NULL,howmany,1,FFTW_FORWARD,FFTW_MEASURE);
  // do stuff  

  return 0;
}

Answer 1

تعمل الوظيفة كما هو موثق.لقد ارتكبت خطأً في طول الخطوة، الذي ينبغي أن يكون في الواقع howmany في هذه الحالة.لقد قمت بتحديث السؤال ليعكس هذا.

أجد أن وثائق FFTW يصعب فهمها إلى حد ما بدون أمثلة (قد أكون أميًا...)، لذلك أقوم بنشر مثال أكثر تفصيلاً أدناه، مقارنة الاستخدام المعتاد لـ fftw_plan_dft_1d مع fftw_plan_many_dft.للتلخيص، في حالة howmany المصفوفات مع الطول N التي يتم تخزينها في كتلة متجاورة من الذاكرة المشار إليها باسم in, ، عناصر المصفوفة j لكل تحويل k نكون

*(in + j*istride + k*idist)

القطعتين التاليتين من التعليمات البرمجية متكافئتان.في الأول، يتم التحويل من بعض المصفوفات ثنائية الأبعاد بشكل صريح، وفي الثاني يتم إجراء التحويل fftw_plan_many_dft يتم استخدام المكالمة لفعل كل شيء في مكانه.

نسخة صريحة

int N, howmany;
// ...
fftw_complex* data = (fftw_complex*) fftw_malloc(N*howmany*sizeof(fftw_complex));
// ... load data with howmany arrays of length N 
int istride, idist;
// ... if data is column-major, set istride=howmany, idist=1
//    if data is row-major, set istride=1, idist=N

fftw_complex* in = (fftw_complex*) fftw_malloc(N*sizeof(fftw_complex));
fftw_complex* out = (fftw_complex*) fftw_malloc(N*sizeof(fftw_complex));
fftw_plan p = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_MEASURE);
for (int k = 0; k < howmany; ++k) {
  for (int j = 0; j < N; ++j) {
    in[j] = data[j*istride + k*idist];
  }
  fftw_execute(p);
  // do something with out
}

خطة كثيرة

int N, howmany;
// ...
fftw_complex* data = (fftw_complex*) fftw_malloc(N*howmany*sizeof(fftw_complex));
// ... load data with howmany arrays of length N 
int istride, idist;
// ... if data is column-major, set istride=howmany, idist=1
//    if data is row-major, set istride=1, idist=N

fftw_plan p = fftw_plan_many_dft(1,&N,howmany,data,NULL,howmany,1,data,NULL,howmany,1,FFTW_FORWARD,FFTW_MEASURE);
fftw_execute(p);