/*--------------------------------------------------------------------------------
離散フーリエ変換と逆離散フーリエ変換
discrete Fourier transform and inverse discrete Fourier transform
DFT.c
*--------------------------------------------------------------------------------*/
#include "DFT.h"
int main( void )
{
int n;
int num;
double power_spec; // パワースペクトル, power spectrum
complex *x; // data to be read x(t) = x(t).real + i * x(t).imag
// x.real: データの実部, real part of the data
// x.imag: データの虚部, imaginary part of the data
complex *F; // F(w) = F(w).real + i * F(w).imag
// F.real: 離散フーリエ変換の実部, real part of discrete Fourier transform
// F.imag: 離散フーリエ変換の虚部, imaginary part of discrete Fourier transform
// 実データには虚部が存在しないので、読み込んだ値はx[t].real に入力し、x[t].imag にはゼロを入力
// As there is no imaginary part in real data, input the read value to x[t].real,
// and input zero to x[t].imag
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;
}
/*----------------------------------------------------------------------*/
// 離散フーリエ変換
// データxを入力し、フーリエ係数Fを出力
// discrete Fourier transform
// input the data x, output the Fourier coefficient F
fourier( num, x, F );
printf("# discrete Fourier transform\n");
printf("# Frequency real imaginary Power Spectrum Amplitude spectrum Phase spectrum\n"); /* */
n = 0;
power_spec = F[n].real * F[n].real + F[n].imag * F[n].imag;
printf("# %d", n );
// printf("# %lf", n * SAMPLING_FREQ/(double)(num) );
printf(" %lf", F[n].real ); /* フーリエ係数の実部, real part of Fourier coefficient */
printf(" %lf", F[n].imag ); /* フーリエ係数の虚部, imaginary part of Fourier coefficient */
printf(" %lf", power_spec ); /* パワースペクトル, power spectrum */
printf(" %lf", 2.0 * sqrt(power_spec)/(double)num ); /* 振幅スペクトル, amplitude spectrum */
// 位相スペクトル
// phase spectrum
if ( F[n].real == 0.0 ) {
printf(" %lf", 0.0 );
printf("\n");
} else {
printf(" %lf", atan(F[n].imag/(double)F[n].real) );
// printf(" %.6g", atan(F[n].imag/(double)F[n].real) );
printf("\n");
}
for( n = 1; n < num; n++ ) {
power_spec = F[n].real * F[n].real + F[n].imag * F[n].imag;
printf(" %d", n );
// printf(" %lf", n * SAMPLING_FREQ/(double)(num) );
printf(" %lf", F[n].real );
printf(" %lf", F[n].imag );
printf(" %lf", power_spec );
printf(" %lf", 2.0 * sqrt(power_spec)/(double)num );
// 位相スペクトル
// phase spectrum
if ( F[n].real == 0.0 ) {
printf(" %lf", 0.0 );
printf("\n");
} else {
printf(" %lf", atan(F[n].imag/(double)F[n].real) );
// printf(" %.6g", atan(F[n].imag/(double)F[n].real) );
printf("\n");
}
}
/*----------------------------------------------------------------------*/
// データxを初期化
// initialise data x
for( n = 0; n < num; n++ ) {
x[n].real = 0.0;
x[n].imag = 0.0;
}
// 逆離散フーリエ変換
// フーリエ係数Fを入力し、データxを出力
// inverse discrete Fourier transform
// input the Fourier coefficient F, and output the data x
inv_fourier( num, x, F );
printf("\n");
printf("# inverse discrete Fourier transform\n");
printf("# t real imaginary\n");
for( n = 0; n < num; n++ ) {
printf(" %d %lf %lf\n", n, x[n].real, x[n].imag);
}
/*----------------------------------------------------------------------*/
free( x );
free( F );
return 0;
}