前回の解答

DFTの関数は以下のようになる．

void dft(double* in, Complex* out, int N) /*function of fft*/
{
　Complex tmp;
　int omega, k;

　for(omega=0; omega<N; omega++) {
　　for(k=0; k<N; k++) {
　　　// 積分の中身の計算
　　　tmp = Cexp( ToComplex(0.0 , 2.0 * pi * omega * k / N) );
　　　tmp = Cmul( ToComplex(in[k], 0.0), tmp);
　　　// Sigmaの計算
　　　out[omega].re += tmp.re;
　　　out[omega].im += tmp.im;
　　}
　}
}

なお，以下のような関数を作ると処理が楽になる．

複素数の中身を初期化（実部，虚部ともゼロを代入）

void Cclear(Complex* a, int n) /*複素数初期化*/
{
　int i;

　for(i=0;i<n;i++){
　　a[i].re = 0;
　　a[i].im = 0;
　}
}

複素数内のデータを表示する

void Cprint(double *a,int n) /*複素数表示*/
{
　int i;
　for(i=0;i<n;i++){
　　printf("%d %lf\n",i,a[i]);
　}
　printf("\n");
}

DFTの結果の入った配列を読み込む．但し，パワースペクトルなので2乗する．

void Sprint(Complex *a, int n) /*スペクトル表示*/
{
　int i;
　for(i=0;i<n/2;i++){
　　printf("%d %lf\n",i,Cabs(a[i])*Cabs(a[i]));
　}
　printf("\n");
}