FFTW3でのインプレースFFT

Question

私はインプレースとアウトオブプレースのFFTをFFTW3でcで理解しようとしています。まず、一次元変換を使うことにしました。私はアウトオブプレイスを最初から複雑な変換に変換してから、その変換からの複雑な出力を使用してインプレース複合を実際の変換に行い、元の実際の入力データを元に戻したいと考えました。私は、しかし、私はちょうどなぜ思っていた？これを行うための私のサンプルコードは以下の通りです。

#include <stdlib.h> 
#include <stdio.h> 
#include <complex.h> 
#include <fftw3.h> 
#include <math.h> 

void one_dimensional_transform(void) 
{ 
    float *in; 
    fftwf_complex *out, *inplace; 
    int N = 8; 
    fftwf_plan plan; 
    int i; 
    FILE *fd; 

    // Do out of place transform. For this, in has N elements and out has N/2 + 1. 
    in = fftwf_alloc_real(N); 
    out = fftwf_alloc_complex(N/2 + 1); 
    inplace = fftwf_alloc_complex(2 * ((N/2) + 1)); 

    plan = fftwf_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE); 

    fd = fopen("input_data.txt", "w"); 

    for(i = 0; i < N; i++) 
    { 
     in[i] = sin(2.0 * M_PI * (2.0 * M_PI * i/N)); 
     fprintf(fd, "%f\n", in[i]); 
    } 
    fclose(fd); 

    fftwf_execute(plan); 
    fftwf_destroy_plan(plan); 
    fftwf_cleanup(); 

    // Do in-place transform 
    for(i = 0; i < N/2 + 1; i++) 
    { 
     inplace[i] = out[i]; 
    } 

    plan = fftwf_plan_dft_c2r_1d(N, (fftwf_complex *)inplace, (float *)inplace, FFTW_ESTIMATE); 

    // I think that I need to include the complex conjugates to the complex version of 
    // inplace (as this is what Mesinger does), and that's why we have the size of 
    // inplace being 2 * (N/2 + 1). 
    for(i = 1; i < N/2; i++) 
    { 
     inplace[N/2 + i] = conjf(inplace[N/2 - i]); 
    } 

    for(i = 0; i < 2 * (N/2 + 1); i++) 
    { 
     inplace[i] /= (float)N; 
     fftwf_execute(plan); 
     fftwf_destroy_plan(plan); 
     fftwf_cleanup(); 

     fd = fopen("inplace_data.txt", "w"); 

     for(i = 0; i < N; i++) 
     { 
      fprintf(fd, "%f\n", crealf(inplace[i])); 
     } 
     fclose(fd); 

     fftwf_free(in); 
     fftwf_free(out); 
     fftwf_free(inplace); 
} 

FFTW3のドキュメントから、実際のデータはサイズ2（N/2 + 1）でなければならないと言います。また、実際の複雑な変換は、エルミートの対称性のために複雑な配列の後半を切り捨てるので、これらの要素を明示的に入れなおすことにしましたが、必要かどうかはわかりません。私は、変換が正規化されていないと言いましたので、インプレース変換の出力をNで再調整しました。また、

for(i = 0; i < 2 * (N/2 + 1); i++) 
{ 
    inplace[i] /= (float)N; 
} 
fftwf_execute(plan); 
fftwf_destroy_plan(plan); 
fftwf_cleanup(); 

：

Answer 1

主な問題は、投稿のコードが実際に起因するあなたは、おそらく次のように実行することを意図し

for(i = 0; i < 2 * (N/2 + 1); i++) 
{ 
    inplace[i] /= (float)N; 
    fftwf_execute(plan); 
    fftwf_destroy_plan(plan); 
    fftwf_cleanup(); 
    ... 
} 

ループにfftwf_execute（！とfftw_destroy）を複数回実行されることです結果はfloatの配列に格納されているので、データを読み込んで保存する必要があります。

実際の部分を crealf(inplace[i])（2番目の値だけを取る）で抽出するのではなく、

のようになります。また、あなたのコメント

に関して最後に

、複雑に本物が原因エルミート対称性のために、複雑な配列の2番目の半分をカットオフを変換するので、私は明示的に戻って、これらの要素を置くことにしましたが、私それが必要かどうかわからない。

私はそれが必要ではないことを保証することができます。実際に変換する複合体は、エルミート対称性を持つという知識に基づいているため、スペクトルの後半は必要ありません。