/*----------------------------------------------------------------------------
複素指数関数を使った離散フーリエ変換と離散逆フーリエ変換
discrete Fourier transform and inverse discrete Fourier transform
using complex exponential function
complex_DFT.c
----------------------------------------------------------------------------*/
#include "complex_DFT.h"
int main( void )
{
int i;
int n;
int num;
double rr; // 複素数の実部, real part of a complex number
double ii; // 複素数の虚部, imaginary part of a complex number
double power_spec;
double temp;
double complex *x; // 入力データ, input data/
double complex *F; // 離散フーリエ変換, discrete Fourier transform
/*----------------------------------------------------------------------*/
if ( scanf( "%d", &num ) != 1 ) {
fputs( "cannot get the data number\n", stderr );
return -1;
}
if (( x = (complex *)malloc( (num) * sizeof(complex) ) ) == NULL ) {
fputs( "cannot allocate memory for x\n", stderr );
return -1;
}
if (( F = (complex *)malloc( (num) * sizeof(complex) )) == NULL ) {
fputs( "cannot allocate memory for F\n", stderr );
return -1;
}
/*----------------------------------------------------------------------*/
// 実データには虚部が存在しないため、実部のみ読み込み。
// As there is no imaginary part in real data, input the read value
ii = 0.0;
for ( i = 0; i < num; i++ ) {
if ( scanf( "%lf", &rr ) != 1 ) {
fputs( "cannot read the data\n", stderr );
return -1;
}
x[i] = rr + I * ii;
}
/*
虚数単位は数学では小文字の i を用いるが、C言語では大文字の I を使用する。
Although imaginary unit uses the lower case i in mathematics,
the uppercase letter I is used in C language.
*/
/*----------------------------------------------------------------------*/
// 離散フーリエ変換
// データxを入力し、フーリエ係数Fを出力
// discrete Fourier transform
// input the data x, output the Fourier coefficient F
complex_fourier( num, x, F );
printf("# discrete Fourier transform\n");
printf("# Frequency real imaginary Power Spectrum Amplitude Spectrum Phase Spectrum\n");
/*----------------------------------------------------------------
conj:
複素数の虚部の符号を反転させて複素共役を計算する。
Calculate the complex conjugate by reversing the sign of
the imaginary part of the complex number.
z = 3.0 + 2.0 * I
conj(z) = 3.0 - 2.0 * I
----------------------------------------------------------------*/
n = 0;
power_spec = F[n] * conj( F[n] );
printf("# %d", n );
// printf("# %lf", n * SAMPLING_FREQ/(double)(num) );
printf(" %lf", creal( F[n] ) ); /* フーリエ係数の実部, real part of Fourier coefficient */
printf(" %lf", cimag( F[n] ) ); /* フーリエ係数の虚部, 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 ( creal( F[n] ) == 0.0 ) {
printf(" %lf", 0.0 );
printf("\n");
} else {
temp = cimag( F[n] )/(double)creal( F[n] );
printf(" %lf", atan( temp ) );
// printf(" %.6g", atan( temp ) );
printf("\n");
}
for( n = 1; n < num; n++ ) {
power_spec = F[n] * conj( F[n] );
printf(" %d", n );
// printf(" %lf", n * SAMPLING_FREQ/(double)(num) );
printf(" %lf", creal( F[n] ) ); /* フーリエ係数の実部, real part of Fourier coefficient */
printf(" %lf", cimag( F[n] ) ); /* フーリエ係数の虚部, 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 ( creal( F[n] ) == 0.0 ) {
printf(" %lf", 0.0 );
printf("\n");
} else {
temp = cimag( F[n] )/(double)creal( F[n] );
printf(" %lf", atan( temp ) );
// printf(" %.6g", atan( temp ) );
printf("\n");
}
}
/*----------------------------------------------------------------------*/
// データxを初期化
// initialise data x
for( n = 0; n < num; n++ ) {
x[n] = 0.0;
}
// 逆離散フーリエ変換
// フーリエ係数Fを入力し、データxを出力
// inverse discrete Fourier transform
// input the Fourier coefficient F, and output the data x
complex_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, creal( x[n] ), cimag( x[n] ) );
}
/*----------------------------------------------------------------------*/
free( x );
free( F );
return 0;
}