/*--------------------------------------------------------------------------------
パーセバルの定理
Parseval's theorem
parseval.c
*--------------------------------------------------------------------------------*/
#include "parseval.h"
int main( void )
{
int n;
int num;
double sum;
double diff;
double power; // 時間軸上のパワー, power on the time axis
double power_spec; // 周波数軸上のパワー, power on frequency axis
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;
}
/*----------------------------------------------------------------------*/
// 時間軸上の時系列データのパワーの計算
// calculating the power of time series data on the time axis
sum = 0.0;
for( n = 0; n < num; n++ ) {
sum += ( x[n].real * x[n].real );
}
power = sum;
/*----------------------------------------------------------------------*/
// 離散フーリエ変換
// データxを入力し、位相Fを出力
// discrete Fourier transform
// input the data x, output the phase F
fourier( num, x, F );
// 周波数軸上のパワーの計算
// calculating the power on frequency axis
sum = 0.0;
for( n = 0; n < num; n++ ) {
sum += ( F[n].real * F[n].real + F[n].imag * F[n].imag );
}
power_spec = sum/(double)num;
/*----------------------------------------------------------------------*/
// 離散フーリエ変換を用いた場合は、パワースペクトルの平均値と
// 時間軸上のパワーは等しくなる。
// When the discrete Fourier transform is used, the mean value of the power spectrum
// obtained by the DFTand the power on the time axis are equal.
diff = power - power_spec;
printf("#-------------------------------------------------------\n");
printf("# 時間軸上のパワー：time power\n");
printf("# 周波数上のパワー：power spectrum\n");
printf("#-------------------------------------------------------\n");
printf("# time power power spectrum gap\n");
printf(" %lf %lf %g\n", power, power_spec, diff);
/*----------------------------------------------------------------------*/
free( x );
free( F );
return 0;
}