3

FFTWを使用して高速な合計を計算しようとしていますが、問題が発生しました。 

int numFreq3 = numFreq*numFreq*numFreq;

FFTW_complex* dummy_sq_fft = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* dummy_sq = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* orig = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_plan dummyPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
                  orig, dummy_sq_fft,
                  FFTW_FORWARD, FFTW_MEASURE );

FFTW_plan dummyInvPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
                      dummy_sq_fft, dummy_sq,
                      FFTW_BACKWARD, FFTW_MEASURE );

for(int i= 0; i < numFreq3; i++) {
  orig[ i ][ 0 ] = sparseProfile02[ 0 ][ i ][ 0 ];
  //img. part == 0
  orig[ i ][ 1 ]  = sparseProfile02[ 0 ][ i ] [ 1 ];
}

FFTW_execute(dummyPlan);
FFTW_execute(dummyInvPlan);

int count = 0; 
for(int i=0; i<numFreq3; i++) {
  double factor = dummy_sq[ i ][ 0 ]/sparseProfile02[ 0 ][ i ][ 0 ];

  if(factor < 0) {
    count++;
  }
}

std::cout<<"Count "<<count<<"\n";

FFTW_free(dummy_sq_fft);
FFTW_free(dummy_sq);
FFTW_destroy_plan(dummyPlan);
FFTW_destroy_plan(dummyInvPlan);

（ここで、sparseProfile02[0]はタイプFFTW_complex*であり、正の実数データのみが含まれています。）

ダミー_sq=IFFT（FFT（sparseProfile02 [0]））があるので、dummy_sq = n ^ 3*sparseProfile02が必要です。しかし、これは一部の時間にのみ当てはまります。実際、sparseProfile02グリッドの対応する値がゼロの場合は常に、dummy_sqグリッドの値は負になります（ただし、その逆はありません）。なぜこれが起こっているのか誰かが知っていますか？

4

2 に答える 2

2

ネクロマンシーに手を出す危険を冒して、fftw ドキュメント (ここ) で、fftw が正規化されないことが明示的に述べられていることに注意する必要があります。したがって、以前に変換された信号の逆変換の結果は、「n」でスケーリングされた元の信号になります。または信号の長さ。

問題になる可能性があります。

于 2011-10-24T05:36:31.087 に答える
1

FFT (順方向と逆方向) には丸め誤差があり、これが問題になっていると思います。一般に、ゼロがプロセスを通じて正確にゼロのままであることを期待すべきではありません (ただし、単純なテスト ケースではゼロになる可能性があります)。あなたのテストループでは、

fabs(dummy_sq[i][0] - numFreq*numFreq*numFreq*sparseProfile02[0][i][0])

あなたのデータの大きさに比べて大きいですか？

実数値を持つサイズ 2 の 1D FFT のみを使用した非常に単純な (病理学的) 例として:

ifft(fft([1e20, 1.0])) != [2e20, 2.0]

1.0 は、倍精度 FFT を使用する 1e20 で失われます。

また、sparseProfile02 でゼロ サンプルで除算すると、factor の NaN が得られる可能性もあります。

于 2010-09-16T06:53:26.673 に答える