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

[鳥取砂丘]

縦軸を振幅（出力電圧）として横軸を時間軸から周波数軸に変化させたいので、フーリエ変換プログラムを作ってみたいのですが・・・

#include <stdio.h>
#include <math.h>

#define PI	3.14
#define MAX	64		

int main (void)
{
	double x[MAX],y[MAX],value;
	int i,j;
	double tmp_x[MAX],tmp_y[MAX];
	FILE *fp,*fp2,*fp3;
	
	
	
	/* ファイル読込み */
	if((fp = fopen("text01.txt","r")) == NULL){
		return(-1);
	}
	if((fp2=fopen("outdate.txt","w")) == NULL){  ←　出力電圧
		return(-1);
	}
	if((fp3=fopen("syuuhasu.txt","w")) == NULL){　　←　周波数
		return(-1);
	}
	
	for(i=0; i<MAX; i++){
		fscanf(fp,"%lf\n",&value);
		x = value;	/* x :実数成分 */
		y = 0;		/* y :虚数成分 */
	}


	/* データを一旦コピー */
	for (i=0; i<MAX; i++) {
		tmp_x = x;
		tmp_y = y;
	}
	
	for (i=0; i<MAX; i++) {
		x = 0;
		y = 0;
		for (j=0; j<MAX; j++) {
			x += tmp_x[j]*cos(2*PI*i*j/MAX) + tmp_y[j]*sin(2*PI*i*j/MAX);
			y += -tmp_x[j]*sin(2*PI*i*j/MAX) + tmp_y[j]*cos(2*PI*i*j/MAX);
		}
	/* DFT */

	for(i=0; i<MAX; i++){

	fprintf(fp2,"%5.5f\n",x[i]);　←　これが出力電圧なはず？
	}
	for(i=0; i<MAX; i++){
		
	fprintf(fp3,"%5.5f\n", sqrt(x[i]*x[i]+y[i]*y[i])*2/MAX);　←これで時間を周波数にできる？
	}

	}
	fclose(fp);
	fclose(fp2);
	fclose(fp3);

return(0);

}

　

式は合っているはずなのですが周波数がこれでOKなのかが疑問です。

気になる点があればドシドシアドバイスください。

[管理人]

計算式については見ていないのですが、二重ループの中で両方iが使われているからおかしくないですか？
実際iはMAXになってぬけているので、外側のループは一度しか通っていないのではないかと思います。

しかし見た感じ二重ループする必要がなさそうなので、、括弧の対応が間違ってるんでしょうか？

[鳥取砂丘]

この部分でしょうか？

for (i=0; i<MAX; i++) {
x = 0;
y = 0;
for (j=0; j<MAX; j++) {
x += tmp_x[j]*cos(2*PI*i*j/MAX) + tmp_y[j]*sin(2*PI*i*j/MAX);
y += -tmp_x[j]*sin(2*PI*i*j/MAX) + tmp_y[j]*cos(2*PI*i*j/MAX);
}

これを

for (i=0; i<MAX; i++) {
x = 0;
y = 0;
　　　　 }

for (j=0; j<MAX; j++) {
x += tmp_x[j]*cos(2*PI*i*j/MAX) + tmp_y[j]*sin(2*PI*i*j/MAX);
y += -tmp_x[j]*sin(2*PI*i*j/MAX) + tmp_y[j]*cos(2*PI*i*j/MAX);
}

このようにするということですか？

[管理人]

計算式を知らないのでなんともいえませんが、そこでiのカウンタのfor文を閉じないと構文的におかしくないでしょうか？
後、memset関数を使えばfor文を使わなくても初期化できますよ。

[鳥取砂丘]

あまりCを理解してないのでforしかわからないです。すいません・・・。

とりあえずこれでやってみようと思います。

[管理人]

memsetは特に難しい関数じゃないので覚えておくと便利ですよ。
memsetは

・第一引数の先頭アドレスから
・第二引数の値で
・第三引数の量だけ埋める

という関数です。
memsetの仕様を見てみると
http://www9.plala.or.jp/sgwr-t/lib/memset.html

「
void *memset(void *buf, int ch, size_t n);

void *buf　：　セット先のメモリブロック
int ch　　：　セットする文字
size_t n　：　セットバイト数
」

と書いてありますよね。
第一引数に初期化したい配列名を書きます。２番目に何の値で初期化したいかその値を、３番目にその量を書きますから、今回の場合では
xを0で初期化して、その量はdouble型のサイズMAX分なのでこうなります。

memset(x,0,sizeof(double)*MAX);

別にforで初期化してもいいですけど、今回のようにループがあちこちあると見にくいような場合にはオススメです。

[しっぽ]

>　fprintf(fp2,"%5.5f\n",x);　←　これが出力電圧なはず？
違います。電圧スペクトルの実部です。

>　fprintf(fp3,"%5.5f\n", sqrt(x*x+y*y)*2/MAX);　←これで時間を周波数にできる？
これが、電圧スペクトルの振幅なのかという意味の質問ならYES（*2が必要かどうかは知りません）。

蛇足
>　x += tmp_x[j]*cos(2*PI*i*j/MAX) + tmp_y[j]*sin(2*PI*i*j/MAX);
>　y += -tmp_x[j]*sin(2*PI*i*j/MAX) + tmp_y[j]*cos(2*PI*i*j/MAX);
x、yはMAXで割る必要があるのでは。

[しっぽ]

訂正

sqrt(x*x+y*y)*2/MAX
でMAXで割っているので、下記は撤回します。

>　x += tmp_x[j]*cos(2*PI*i*j/MAX) + tmp_y[j]*sin(2*PI*i*j/MAX);
>　y += -tmp_x[j]*sin(2*PI*i*j/MAX) + tmp_y[j]*cos(2*PI*i*j/MAX);
x、yはMAXで割る必要があるのでは。

[tk-xleader]

->管理人さん

double型変数の配列にmemset()は危険だと思います。全ビット0の状態が必ずしもdouble型の0になるとは限らないようです。

[鳥取砂丘]

ということは・・・

for(i=0; i<MAX; i++){
fscanf(fp,"%lf\n",&value);
x = value; /* x :実数成分 */
y = 0; /* y :虚数成分 */
}


/* データを一旦コピー */
for (i=0; i<MAX; i++) {
tmp_x = x;
tmp_y = y;
}

for (i=0; i<MAX; i++) {
x = 0;
y = 0;
for (j=0; j<MAX; j++) {
x/MAX += tmp_x[j]*cos(2*PI*i*j/MAX) + tmp_y[j]*sin(2*PI*i*j/MAX);←改善
y/MAX += -tmp_x[j]*sin(2*PI*i*j/MAX) + tmp_y[j]*cos(2*PI*i*j/MAX);←改善
　　　　　　　}
/* DFT */

for(i=0; i<MAX; i++){

fprintf(fp2,"%.5f\n",sqrt(x[i]*x[i]+y[i]*y[i])*2/MAX);　←　これが出力？
}
for(i=0; i<MAX; i++){

fprintf(fp3,"%.5f\n", (double)i/100);←サンプリング件数/サンプリング周波数＝100から、これ　　　　　　　　　　　　　　　　　　　　　　　　　が周波数？
}

といった感じにすればいいのでしょうか？

二重ループの問題を指摘されたりもしましたがその点も改善すべきでしょうか？（管理人さんの指摘より）

[しっぽ]

>　x/MAX += tmp_x[j]*cos(2*PI*i*j/MAX) + tmp_y[j]*sin(2*PI*i*j/MAX);←改善
>　y/MAX += -tmp_x[j]*sin(2*PI*i*j/MAX) + tmp_y[j]*cos(2*PI*i*j/MAX);←改善
ゴメンなさい。私が変なことを言ったため、この部分はペケです。
元のままでいいです（スペクトルを書き込む段階でMAXで除しているから）。

>　fprintf(fp2,"%.5f\n",sqrt(x*x+y*y)*2/MAX);　←　これが出力？
電圧スペクトルの大きさだけ（位相は無視）でよいなら、YES。

>　fprintf(fp3,"%.5f\n", (double)i/100);←サンプリング件数/サンプリング周波数＝100から、
>　これが周波数？
時間軸でのサンプリング周期＝1（サンプリング周波数＝1）かつサンプリング件数=100という意味なら、
YES。

>　二重ループの問題を指摘されたりもしましたがその点も改善すべきでしょうか？（管理人さんの指摘より）
変更の必要ないです。好みの領域の問題だと思います。

[鳥取砂丘]

ありがとうございました。

うまくできたような感じです。

みなさんありがとうございました！