離散フーリエ変換のプログラムについて

質問

今ＤＦＴのプログラムをＣ言語で書いているんですが、うまく動いてくれません。
ＤＦＴの式はＸ(k)=1/Ｎ{Σx(k)*e^(-j2πkn/Ｎ)}のシグマの中をn=0からＮ-1まで足し合わせればいいと思っているのですがちがいますか。
下にプログラムを書きますのでお願いします。

void DFTcore(double x[],int N,double A[],double fai[],double a_rl[],double a_im[]){

　　　　ｘ[]はデータ
　　　　Ｎはデータ数
　　　　Ａ[]は振幅
　　　　fai[]は位相差
　　　　a_rl[]は実数部、a_im[]は虚数部

double w,a;
int k,n,t;

for(k=0;k<N;k++){　　//初期化
　　　a_rl[k]=0.0;　　//実数部
　　 a_im[k]=0.0;　　//虚数部
　　}
　　　　時間間隔は1秒
for(k=0;k<N;k++){
　　　　　　for(n=0;n<N;n++){
　 w=((2*PI)/N)*k*n;
　　　　　　　 a_rl[k]=a_rl[k]+x[k]*cos(w);
　 a_im[k]=-a_im[k]+x[k]*sin(w);
　}
　　a_rl[k]=a_rl[k]/(double)N;実数部
　　a_im[k]=-a_im[k]/(double)N;　虚数

　　A[k]=sqrt(a_rl[k]*a_rl[k]+a_im[k]*a_im[k]); //振幅
　　　fai[k]=atan(a_im[k]/a_rl[k]); //位相

}
}

回答No.2

anisol · ベストアンサー

正しくは、Ｘ(k)=1/Ｎ{Σx(n)*e^(-j2πkn/Ｎ)}
(ただし、k=f*N/f_sample , fは調べたい周波数、f_sampleはサンプリング周波数)
です。
とりあえずは、
a_rl[k]=a_rl[k]+x[n]*cos(w);
a_im[k]=a_im[k]+x[n]*sin(w);
とすれば動作確認はできます。
あとはご自分のプログラムの中でのkの意味(f*N/f_sampleは整数でなくてもよい)、f<(f_sample/2)、に注意し、プログラムを書きなおしてみてください。またきれいなスペクトルを得るには窓関数処理をお忘れなく。

数学、数式が大の苦手の私ですが、DFTはおぼろげに理解できたので、以前N88BASICという言語で(今は知っている人が少ないらしい)、DFTのプログラムを書いたことがあります。なにせクロック周波数12MHzのパソコンだったので、サインテーブルを作って積和演算の部分はマシン語で書いたり、いろいろやりました。
きれいなスペクトルが出ると感動しますよ。ご健闘を祈ります。

お礼

プログラム関する問題は解決しました、ありがとうございました。
しかし実験データから振幅の大きな周波数の周りにそれよりは小さな振幅がでているので窓関数処理をしたほうがよいみたいです。ご指摘ありがとうございました。

回答No.1

motsuan

Σx（t）exp(-i*w*t)/N 　（tについての和）
をやりたいんだと思うのですが、その場合
一つのwに対してtについての和をとるわけですから
a_rl[k]=a_rl[k]+x[k]*cos(w);
a_im[k]=-a_im[k]+x[k]*sin(w);
で a_rl[]、a_im[]、x[] の [] の中が
すべて k になっているのはおかしいですよね。
さらに
a_im[k]=-a_im[k]+x[k]*sin(w);
だと1回の代入のたびに符号を反転するようになっているので
たぶんめちゃくちゃになるでしょう。
a_im[k]=a_im[k]+x[k]*sin(-w)
　　　　=a_im[k]-x[k]*sin(w)
ということがしたいのではないでしょうか？
C++をつかっているようなので
STLの<complex>を使ってみてもいいんじゃないでしょうか？

お礼

ご指摘ありがとうございます。
どうもｎとｋの単純な間違いだったみたいなので、お騒がせして申し訳ありませんでした。いまは正常にプログラムは走っていて、これから画像処理のほうに応用していく予定でいき、ＦＦＴにも手をだして行きたいと思っています。ありがとうございました。