FFTW3の謎の挙動
下記のようなプログラムを実行すると、out[1]からout[3]はちゃんと出るのだが、out[0]の出力だけ総0になってしまう。data[i]の中身を入れ替えてみてもout[0]だけ総0。なぜなんだ。
→解決しました。プラン作成時に入出力メモリが破壊されるのが原因。忘れてた…。ちゃんとconst使えってことなんですが、FFTW使ってると『ここは非constだけど非破壊』みたいなのをつい使ってしまい、だんだんいい加減に…
環境: ubuntu 10.04 1 LTS, g++ 4.4.3, fftw3 3.2.2-1
#include <complex> #include <iostream> #include <cmath> #include <fftw3.h> #include <valarray> using namespace std; template <typename T> T round2zero(T x){ return abs(x)<pow(2.0, -47) ? 0 : x ; } template <typename T> void print(const T &v){ for(unsigned i=0; i<v.size(); ++i){ cout<<round2zero(v[i])<<", "; } cout<<endl<<endl; } int main(int argc, char **argv) { const unsigned num=4, size=32; valarray<double> data[num]; valarray<complex<double> > out[num]; for(unsigned i=0; i<num; ++i){ data[i].resize(size); out[i].resize(size); // data[i][j]に波数(j+1)のsin波を代入 for(unsigned j=0; j<size; ++j) data[i][j]=sin(2*M_PI*(i+1)*j/static_cast<double>(size)); print(data[i]); } fftw_plan p[num]; for(unsigned i=0; i<num; ++i){ p[i] = fftw_plan_dft_r2c_1d(size, &data[i][0], reinterpret_cast<fftw_complex*>(&out[i][0]), FFTW_MEASURE); fftw_execute(p[i]); print(out[i]); } }