/*--------------------------------------------------------------------------------
高速フーリエ変換と逆高速フーリエ変換
Fast Fourier transform and inverse Fast Fourier transform
FFT.c
*--------------------------------------------------------------------------------*/
#include "FFT.h"
int main(void)
{
int k;
int num;
double power_spec;
double sf = (double)SAMPLING_FREQ;
double temp;
complex *x;
complex *F;
// x.real: 実部, real part
// x.imag: 虚部, imaginary part
/*--------------------------------------------------------------------------*/
// 入力データ読み込み
// read the input data
if ( ( x = readdata( &num ) ) == NULL ) {
fputs( "readdata error\n", stderr );
return -1;
}
if (( F = (complex *)malloc( num * sizeof(complex) )) == NULL ) {
fputs( "cannot allocate memory for F\n", stderr );
return -1;
}
/*--------------------------------------------------------------------------*/
// 入力データ数が2のべき乗 (2^n) であるかどうかを判定
// check whether the number of input data is a power of 2 (2^n)
if (( (num)&(num-1) ) != 0 ) {
printf("#-------------------------------------------------------\n");
printf(" the input data number must be power of two, 2^n.\n");
printf("#-------------------------------------------------------\n");
return -1;
}
/*--------------------------------------------------------------------------*/
// 高速フーリエ変換
// データxを入力し、フーリエ係数Fを出力
// Fast Fourier transform
// input the data x, output the Fourier coefficient F
FFT( num, x, F );
printf("# Fast Fourier Transform\n");
printf("# k real imaginary power\n");
k = 0;
power_spec = F[k].real * F[k].real + F[k].imag * F[k].imag;
printf("# %d %lf %lf %lf\n", k, F[k].real, F[k].imag, power_spec );
// printf("# %lf %lf %lf %.6g\n", k * SAMPLING_FREQ/(double)(num), F[k].real, F[k].imag, power_spec );
for( k = 1; k < num; k++ ) {
power_spec = F[k].real * F[k].real + F[k].imag * F[k].imag;
temp = k * sf/(double)(num);
printf(" %d %lf %lf %lf\n", k, F[k].real, F[k].imag, power_spec );
// printf(" %lf %lf %lf %.6g\n", temp, F[k].real, F[k].imag, power_spec );
}
/*--------------------------------------------------------------------------*/
// 逆高速フーリエ変換
// フーリエ係数Fを入力し、データをxに出力
// inverse Fast Fourier transform
// input the Fourier coefficient x, output the data x
inv_FFT( num, x, F );
printf("\n");
printf("# inverse Fast Fourier Transform\n");
printf("# k real imaginary\n");
for( k = 0; k < num; k++ ) {
printf(" %d %lf %lf\n", k, x[k].real, x[k].imag );
}
/*--------------------------------------------------------------------------*/
free( x );
free( F );
return 0;
}