今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]);//位相
}
}
No.2ベストアンサー
- 回答日時:
正しくは、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のパソコンだったので、サインテーブルを作って積和演算の部分はマシン語で書いたり、いろいろやりました。
きれいなスペクトルが出ると感動しますよ。ご健闘を祈ります。
プログラム関する問題は解決しました、ありがとうございました。
しかし実験データから振幅の大きな周波数の周りにそれよりは小さな振幅がでているので窓関数処理をしたほうがよいみたいです。ご指摘ありがとうございました。
No.1
- 回答日時:
Σ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>を使ってみてもいいんじゃないでしょうか?
ご指摘ありがとうございます。
どうもnとkの単純な間違いだったみたいなので、お騒がせして申し訳ありませんでした。いまは正常にプログラムは走っていて、これから画像処理のほうに応用していく予定でいき、FFTにも手をだして行きたいと思っています。ありがとうございました。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- C言語・C++・C# 10個の実数に対する降順ソート結果を出力するプログラムを作りたいのですが、以下のプログラムをどう直せ 1 2022/07/09 22:16
- C言語・C++・C# c言語でユーザ関数を利用して複素数のべき乗と絶対値の数列を計算するプログラムが作りたいです。 3 2023/01/29 22:13
- C言語・C++・C# LU分解法のピボッティングについて(C言語/gcc-9) 3 2022/07/11 23:10
- 数学 単振り子とルンゲ・タック法 1 2022/07/15 00:05
- C言語・C++・C# LU分解法のピボット選択機能実装について(C言語・gcc-9) 1 2022/07/22 15:20
- 経済学 資本移動や価格変動のない次のような固定為替レート・モデルを考える。 C = 10 + 0.8 Y I 3 2022/06/21 20:50
- C言語・C++・C# Cのdoubleの浮動小数点表示について 3 2023/04/17 13:14
- C言語・C++・C# C言語初心者 構造体 課題について 1 2023/03/10 19:30
- C言語・C++・C# プログラミングの授業の課題です 1 2023/01/17 22:15
- Visual Basic(VBA) VBAプログラミング 2 2022/11/27 12:07
関連するカテゴリからQ&Aを探す
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
c言語で、繰り返し文の中で、0....
-
プログラムでの数字につく”f”の...
-
C言語のプログラムで#include<m...
-
C言語で直角三角形の斜辺を求め...
-
sin(x)の近似について
-
doubleの変数にintとintの割り...
-
Cで3乗根を求める方法
-
C言語初心者 構造体 課題について
-
至急です! マクロ定義で #defi...
-
float型とdouble型の変数の違い...
-
C言語を実行すると-infが出てき...
-
2分法で方程式の複数の解を自...
-
MATLABで画像のヒストグラムを...
-
浮動小数点数が表示されないん...
-
関数におけるif文とreturn文に...
-
2次方程式の解を求めるプログ...
-
double型とint型で三分の一乗の...
-
floating point not loadedとは?
-
C言語 関数プロトタイプ宣言の...
-
C言語で表記についの質問です
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
プログラムでの数字につく”f”の...
-
doubleの変数にintとintの割り...
-
C 開放してるのにエラー(doubl...
-
Cで3乗根を求める方法
-
float型とdouble型の変数の違い...
-
至急です! マクロ定義で #defi...
-
C言語の型による処理速度の違い
-
int とdoubleの比較
-
関数におけるif文とreturn文に...
-
C言語初心者 構造体 課題について
-
c言語のコンパイルエラー canno...
-
C言語 関数プロトタイプ宣言の...
-
C言語を実行すると-infが出てき...
-
float?数字の後にLがつくもの
-
数値を指数部と仮数部に分離したい
-
difftime()について
-
浮動小数点数が表示されないん...
-
たくさんの数の平均を求める方...
-
DWORDの警告
-
-1.#IND00と出てしまうのですが...
おすすめ情報