OKWAVEのAI「あい」が美容・健康の悩みに最適な回答をご提案!
-PR-
解決
済み

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

  • 暇なときにでも
  • 質問No.100563
  • 閲覧数896
  • ありがとう数4
  • 気になる数0
  • 回答数2
  • コメント数0

お礼率 40% (2/5)

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

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

    x[]はデータ
    Nはデータ数
    A[]は振幅
    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]);//位相

}
}
通報する
  • 回答数2
  • 気になる
    質問をブックマークします。
    マイページでまとめて確認できます。

質問者が選んだベストアンサー

  • 回答No.2
レベル11

ベストアンサー率 48% (146/301)

正しくは、X(k)=1/N{Σx(n)*e^(-j2πkn/N)}
(ただし、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のパソコンだったので、サインテーブルを作って積和演算の部分はマシン語で書いたり、いろいろやりました。
きれいなスペクトルが出ると感動しますよ。ご健闘を祈ります。
 
お礼コメント
tani92

お礼率 40% (2/5)

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

その他の回答 (全1件)

  • 回答No.1
レベル10

ベストアンサー率 40% (54/135)

Σ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_i ...続きを読む
Σ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>を使ってみてもいいんじゃないでしょうか?
お礼コメント
tani92

お礼率 40% (2/5)

ご指摘ありがとうございます。
どうもnとkの単純な間違いだったみたいなので、お騒がせして申し訳ありませんでした。いまは正常にプログラムは走っていて、これから画像処理のほうに応用していく予定でいき、FFTにも手をだして行きたいと思っています。ありがとうございました。
投稿日時 - 2001-07-09 01:18:32


このQ&Aのテーマ
このQ&Aで解決しましたか?
関連するQ&A
-PR-
-PR-
こんな書き方もあるよ!この情報は知ってる?あなたの知識を教えて!
このQ&Aにはまだコメントがありません。
あなたの思ったこと、知っていることをここにコメントしてみましょう。

その他の関連するQ&A、テーマをキーワードで探す

キーワードでQ&A、テーマを検索する
-PR-
-PR-
-PR-

特集


いま みんなが気になるQ&A

関連するQ&A

-PR-

ピックアップ

-PR-
ページ先頭へ