/*---------------------------------------------------------------------------------------------------
ウィーナー・ヒンチンの定理（非巡回的データを用いて）
～離散フーリエ変換を使って自己相関関数を計算～
Wiener-Khinchin theorem (using non-cyclic data)
~calculate the auto-correlation function using the discrete Fourier transform~
WKT_AC.c
---------------------------------------------------------------------------------------------------*/
#include "WKT_AC.h"
int main( void )
{
int i;
int n;
int num;
int num2;
double sum_Cxx;
double diff;
double *Cxx; // データから直接求める自己相関関数, auto-correlation function obtained directly from the data
double *xx; // 自己相関関数を求めるためデータ, data for calculating the auto-correlation function
double *power_spec; // データから直接求めるパワースペクトル, power spectrum obtained directly from data
complex *x; // 読み込むデータ data to be read x(t) = x(t).real + i * x(t).imag
// x.real: データの実部, real part of data
// x.imag: データの虚部, imaginary part of data
complex *F; // 離散フーリエ変換, discrete Fourier transform F(w) = F(w).real + i * F(w).imag
// F.real: フーリエ変換の実部, real part of the Fourier transform
// F.imag: フーリエ変換の虚部, imaginary part of the Fourier transform
/*----------------------------------------------------------------------*/
// read the number of the input data
if ( scanf( "%d", &num ) != 1 ) {
fputs( "cannot get the data number\n", stderr );
return -1;
}
num2 = num + num;
if (( x = (complex *)malloc( (num2) * sizeof(complex) ) ) == NULL ) {
fputs( "cannot allocate memory for x\n", stderr );
return -1;
}
if (( F = (complex *)malloc( (num2) * sizeof(complex) )) == NULL ) {
fputs( "cannot allocate memory for F\n", stderr );
return -1;
}
if (( power_spec = (double*)malloc( (num2) * sizeof(double))) == NULL ) {
fputs( "cannot allocate memory for power_spec\n", stderr );
return -1;
}
if (( xx = (double*)malloc( (num2) * sizeof(double))) == NULL ) {
fputs( "cannot allocate memory for xx\n", stderr );
return -1;
}
if (( Cxx = (double*)malloc( (num2) * sizeof(double))) == NULL ) {
fputs( "cannot allocate memory for Cxx\n", stderr );
return -1;
}
/*----------------------------------------------------------------------*/
// 非巡回データの設定
// setting up non-cyclic data
for ( i = 0; i < num2; i++ ){
x[i].real = 0.0;
x[i].imag = 0.0;
}
for ( i = 0; i < num; i++ ) {
if ( scanf( "%lf", &x[i].real ) != 1 ) {
fputs( "cannot read the data\n", stderr );
return -1;
}
}
/*
input data = 1, 2, 3
x.real = xx = 1, 2, 3, 0, 0, 0
*/
for ( i = 0; i < num2; i++ ){
xx[i] = x[i].real;
}
/*----------------------------------------------------------------------*/
/*
入力されたデータ数は num である。しかし、fourier()の num2 を num とすると、
フーリエ変換の際にデータxが巡回データとして認識される。
そのため、強制的に非巡回データを入力するため num2 を使用する。
The number of input data is num. However, if num2 in fourier() is set to num,
the Fourier transform recognises data x as cyclic data.
Therefore, num2 is used to force the input of non-cyclic data.
*/
// 離散フーリエ変換
// データxを入力し、位相Fを出力
// discrete Fourier transform
// input data x and output phase F
fourier( num2, x, F );
// パワースペクトルの計算
// F[n].realにパワースペクトルを入力、F[n].imagにゼロを入力
// calculate power spectrum
// input the power spectrum to F[n].real, and input zero to F[n].imag
for( n = 0; n < num2; n++ ) {
power_spec[n] = F[n].real * F[n].real + F[n].imag * F[n].imag;
F[n].real = power_spec[n];
F[n].imag = 0.0;
}
/*----------------------------------------------------------------------
パワースペクトルを逆フーリエ変換し、自己相関関数を計算
inverse Fourier transform of the power spectrum
and calculation of the autocorrelation function
----------------------------------------------------------------------*/
// 逆フーリエ変換
// 位相Fを入力し、データxを出力
// 入力するF.realはパワースペクトル、出力されるx[n].realに自己相関関数が入る
// inverse Fourier transform
// input phase F, and output data x
// input data F is the power spectrum, and the output x is the auto-correlation function
inv_fourier( num2, x, F );
/*----------------------------------------------------------------------*/
// 自己相関関数をデータから直接計算
// x[n].real: ウィーナー・ヒンチンの定理から求めた自己相関関数
// Cxx: データから直接求めた自己相関関数
// calculate auto-correlation function directly from the data
// x[n].real: auto-correlation function obtained by Wiener-Khinchin theorem
// Cxx: auto-correlation function obtained directly from the data
printf("#---------------------------------------------------------------------------\n");
printf("# ウィーナー・ヒンチンの定理（非巡回的データを使用）\n");
printf("#---------------------------------------------------------------------------\n");
printf("# パワースペクトルを逆フーリエ変換し、自己相関関数を計算\n");
printf("# WKT ACF: パワースペクトルから求めた自己相関関数（ウィーナー・ヒンチンの定理）\n");
printf("# DATA ACF: データから直接求めた自己相関関数\n");
printf("# gap: WKT ACF - DATA ACF\n");
printf("#---------------------------------------------------------------------------\n");
printf("#---------------------------------------------------------------------------\n");
printf("# Wiener-Khinchin theorem (using non-cyclic data)\n");
printf("#---------------------------------------------------------------------------\n");
printf("# inverse Fourier transform of the power spectrum and calculation of the auto-correlation function\n");
printf("# WKT ACF: auto-correlation function obtained from the power spectrum (Wiener-Khinchin theorem)\n");
printf("# DATA ACF: auto-correlation function obtained directly from the data\n");
printf("# gap: WKT ACF - DATA ACF\n");
printf("#---------------------------------------------------------------------------\n");
printf("# lag WKT ACF DATA ACF gap\n");
for( n = 0; n < num; n++ ){
sum_Cxx = 0.0;
for( i = 0; i < num; i++ ){
sum_Cxx += ( xx[n+i] * xx[i] );
}
Cxx[n] = sum_Cxx; /* ### */
}
for( n = 0; n < num; n++ ){
diff = x[n].real - Cxx[n];
printf(" %d %lf %lf %g\n", n, x[n].real, Cxx[n], diff );
}
printf("#---------------------------------------------------------------------------\n");
/*------------------------------------------------------------------------------------------
自己相関関数をフーリエ変換し、パワースペクトルを計算
Fourier transform the autocorrelation function and calculate the power spectrum
------------------------------------------------------------------------------------------*/
// 自己相関関数をフーリエ変換するため、自己相関関数をデータの半分で対称化
// To Fourier transform the autocorrelation function, the autocorrelation function
// is symmetrized by half of the data
// ### で求めた自己相関関数を利用
// use the auto-correlation function obtained in ###
for( n = 1; n < num; n++ ){
Cxx[num2-n] = Cxx[n];
}
for( n = 0; n < num2; n++ ){
x[n].real = Cxx[n];
x[n].imag = 0.0;
}
// 離散フーリエ変換
// データxを入力し、位相Fを出力
// 入力するxは自己相関関数、出力されるFはパワースペクトルが入る
// discrete Fourier transform
// input data x and output phase F
// input data x is the auto-correlation function, and the output F is power spectrum
fourier( num2, x, F );
printf("# 自己相関関数をフーリエ変換し、パワースペクトルを計算\n");
printf("# DATA PS: データから直接求めたパワースペクトル\n");
printf("# WKT PS: 自己相関関数から求めたパワースペクトル（ウィーナー・ヒンチンの定理）\n");
printf("# gap: DATA PS - WKT PS\n");
printf("#---------------------------------------------------------------------------\n");
printf("# Fourier transform the autocorrelation function and calculate the power spectrum\n");
printf("# DATA PS: power spectrum obtained directly from the data\n");
printf("# WKT PS: power spectrum obtained from the auto-correlation function (Wiener-Khinchin theorem)\n");
printf("# gap: DATA PS - WKT PS\n");
printf("#---------------------------------------------------------------------------\n");
printf("# freq DATA PS WKT PS gap\n");
for( n = 0; n < num; n++ ) {
diff = power_spec[n] - F[n].real;
printf(" %d %lf %lf %g\n", n, power_spec[n], F[n].real, diff );
}
/*----------------------------------------------------------------------*/
free( x );
free( xx );
free( Cxx );
free( F );
free( power_spec );
return 0;
}