フーリエ変換のC言語プログラムについて
正弦波（およびガウス性雑音)をフーリエ変換(離散)→逆フーリエ変換するというプログラムを組みました。正弦波をフーリエ変換すると実部は２回ピークがくるはずなのですが、すべて「0.000000」または「-0.000000」と表示されてしまいます。虚部は正常なようで実装の仕方もさほど違わないので、何が問題なのかわからずにいます。念のためコードはすべて載せますが、該当箇所は関数Fourierの
fp = fopen("reX.txt", "w"); //書き込む
あたりです。問題点を教えていただけないでしょうか。お願いします。
//gauss.txt, sin.txt:発生させたガウス性雑音＆正弦波
//reX, imX:フーリエ変換の実部と虚部
//re-1, im-1:逆フーリエ変換の実部と虚部
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define PI 3.14159265358979323846
#define N 256
//n:長さ, w:角周波数, p:位相(phase), a:振幅
void SinCurve(int n, double w, double p, double a)
{
FILE *fp;
double x;
int t;
fp = fopen("sin.txt", "w"); //書き込むので"w"
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(t = 0; t < n; t++)
{
x = a * sin( w*(double)t + p );
fprintf(fp, "%f\n", x);
}
}
fclose(fp);
}
//n:長さ, s:分散, m:平均
void Gauss(int n, double s, double m)
{
FILE *fp;
double x, x1, x2, y1;
int t;
srand((unsigned) time(NULL));
fp = fopen("gauss.txt", "w"); //書き込むので"w"
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(t = 0; t < n; t++)
{
x1 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0);
x2 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0);
y1 = pow(-2.0*log(x1), 0.5) * cos(2.0*PI*x2);
y1 = s * y1 + m;
fprintf(fp, "%f\n", y1);
}
}
fclose(fp);
}
//ファイル名sのデータをフーリエ変換し、実部と虚部をreX, imXに保存
void Fourier(int num, char *s)
{
FILE *fp;
int k, n;
double largeX, x[N+100], t;
fp = fopen(s, "r"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
// printf("%s\n", s);
for(k = 0; k < num; k++)
{
fscanf(fp, "%lf", &x[k]);
printf("x[%d]=%f\n", k, x[k]);
}
}
fp = fopen("reX.txt", "w"); //書き込む
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
largeX = 0.0;
t = 2.0*PI*(double)k / (double)N;
for(n = 0; n < num; n++)
{
largeX += x[n] * cos((double)n*t);
// printf("%f\n", largeX);
}
fprintf(fp, "%f\n", largeX);
printf("reX[%d]=%f\n", k, largeX);
}
}
fp = fopen("imX.txt", "w"); //書き込む
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
largeX = 0.0;
t = 2.0*PI*k / (double)N;
for(n = 0; n < num; n++)
{
largeX -= x[n] * sin(n*t);
}
fprintf(fp, "%f\n", largeX);
}
}
fclose(fp);
}
void InverseFourier(int num)
{
FILE *fp;
int k, n;
double a[N+100], b[N+100], x, t;
//a:reX, b:imX
fp = fopen("reX.txt", "r"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
fscanf(fp, "%lf", &a[k]);
// printf("a[%d]=%f\n", k, a[k]);
}
}
fp = fopen("imX.txt", "r"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
fscanf(fp, "%lf", &b[k]);
// printf("b[%d]=%f\n", k, b[k]);
}
}
fp = fopen("re-1.txt", "w"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(n = 0; n < num; n++)
{
x = 0.0;
t = 2.0*PI*(double)n / (double)N;
for(k = 0; k < num; k++)
{
x +=a[k] *cos(k*t) - b[k] *sin(k*t);
}
x /= (double)N;
fprintf(fp, "%f\n", x);
// printf("x[%d]=%f\n", n, x);
}
}
/*
fp = fopen("im-1.txt", "w"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(n = 0; n < num; n++)
{
x = 0.0;
for(k = 0; k < num; k++)
{
t = 2.0*PI*(double)k / (double)N;
x = x + a[k] *sin(n*t) + b[k] *cos(n*t);
}
x /= (double)N;
fprintf(fp, "%f\n", x);
}
}
*/
fclose(fp);
}
int main(void)
{
SinCurve(N, PI/8.0, 0.0, 1.0);
// Gauss(N, 1.0, 0.0);
Fourier(N, "sin.txt");
// Fourier(N, "gauss.txt");
InverseFourier(N);
return 0;
}
お礼
ありがとうございます。できました。％ｌｆとしたおかげかポインタの概念がよくわかってなかったのかわかりませんが、とりあえずできました。ありがとうございました。