/*---------------------------------------------------------------------------------------------------
ウィーナー・ヒンチンの定理（巡回的データを用いて）
～離散フーリエ変換を使って相互相関関数を計算～
Wiener-Khinchin theorem (using cyclic data)
~calculate the cross-correlation function using the discrete Fourier transform~
WKT_CC.c
---------------------------------------------------------------------------------------------------*/
#include "WKT_CC.h"
int main( void )
{
int i;
int n;
int num;
int num2;
double sum_Cxy;
double diffx;
double diffy;
double *Cxy; // データから直接求める相互相関関数, cross-correlation function obtained directly from the data
double *xx; // 相互相関関数を求めるためデータ, data for calculating the cross-correlation function
double *yy; // 相互相関関数を求めるためデータ, data for calculating the cross-correlation function
complex *power_spec; // データから直接求めるパワースペクトル, power spectrum obtained directly from data
// power_spec.real: パワーの実部, real part of power
// power_spec.imag: パワーの虚部, imaginary part of power
complex *x; // 読み込むデータ data to be read x(t) = x(t).real + i * x(t).imag
complex *y; // 読み込むデータ data to be read y(t) = y(t).real + i * y(t).imag
// x.real: データの実部, real part of data
// x.imag: データの虚部, imaginary part of data
complex *Fx; // 離散フーリエ変換 discrete Fourier transform F(w) = F(w).real+ i * F(w).imag
complex *Fy; // 離散フーリエ変換 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
/*----------------------------------------------------------------------*/
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 (( y = (complex *)malloc( (num) * sizeof(complex) ) ) == NULL ) {
fputs( "cannot allocate memory for y\n", stderr );
return -1;
}
if (( Fx = (complex *)malloc( (num) * sizeof(complex) )) == NULL ) {
fputs( "cannot allocate memory for Fx\n", stderr );
return -1;
}
if (( Fy = (complex *)malloc( (num) * sizeof(complex) )) == NULL ) {
fputs( "cannot allocate memory for Fy\n", stderr );
return -1;
}
if (( power_spec = (complex *)malloc( (num) * sizeof(complex))) == NULL ) {
fputs( "cannot allocate memory for power_spec\n", stderr );
return -1;
}
if (( Cxy = (double*)malloc( (num) * sizeof(double))) == NULL ) {
fputs( "cannot allocate memory for Cxy\n", stderr );
return -1;
}
num2 = num + num;
if (( xx = (double*)malloc( (num2) * sizeof(double))) == NULL ) {
fputs( "cannot allocate memory for xx\n", stderr );
return -1;
}
if (( yy = (double*)malloc( (num2) * sizeof(double))) == NULL ) {
fputs( "cannot allocate memory for yy\n", stderr );
return -1;
}
/*----------------------------------------------------------------------*/
// 実データには虚部が存在しないので、読み込んだ値はx[i].realとy[i].realに入力し、
// x[i].imagとy[i].imag にはゼロを入力
// As there is no imaginary part in real data, input the read value to x[i].real and y[i].real,
// and input zero to x[i].imag and y[i].imag
for ( i = 0; i < num; i++ ) {
if ( scanf( "%lf %lf", &x[i].real, &y[i].real ) != 2 ) {
fputs( "cannot read the data\n", stderr );
return -1;
}
x[i].imag = 0.0;
y[i].imag = 0.0;
}
// 巡回データの設定
// setting up cyclic data
/*
x.real = 1, 2, 3
xx = 1, 2, 3, 1, 2, 3
*/
for ( i = 0; i < num; i++ ){
xx[i] = x[i].real;
yy[i] = y[i].real;
xx[num+i] = xx[i];
yy[num+i] = yy[i];
}
/*----------------------------------------------------------------------*/
// 離散フーリエ変換
// データxを入力し、位相Fxを出力
// discrete Fourier transform
// input data x and output phase Fx
fourier( num, x, Fx );
// 離散フーリエ変換
// データyを入力し、位相Fyを出力
// discrete Fourier transform
// input data y and output phase Fy
fourier( num, y, Fy );
/*
for( n = 0; n < num; n++ ) {
printf("%lf %lf\n", Fx[n].real, Fx[n].imag );
}
printf("#----------------------------------------\n");
for( n = 0; n < num; n++ ) {
printf("%lf %lf\n", Fy[n].real, Fy[n].imag );
}
*/
// クロススペクトル：Fx[]の複素共役とFy[]の積
// Cross spectrum: product of complex conjugate of Fx[] and Fy[].
for( n = 0; n < num; n++ ) {
power_spec[n].real = Fx[n].real * Fy[n].real + Fx[n].imag * Fy[n].imag;
power_spec[n].imag = Fx[n].imag * Fy[n].real - Fx[n].real * Fy[n].imag;
Fx[n].real = power_spec[n].real;
Fx[n].imag = power_spec[n].imag;
}
/*----------------------------------------------------------------------
パワースペクトルを逆フーリエ変換し、相互相関関数を計算
inverse Fourier transform of the power spectrum
and calculation of the cross-correlation function
----------------------------------------------------------------------*/
/* フーリエ逆変換 */
/* 位相Fxを入力し、データxを出力 */
/* x[n].realに相互相関関数が入る */
/* x[n].imagの値はほぼゼロ */
// inverse Fourier transform
// input phase Fx, and output data x
// Fx is the power spectrum, and x is the cross-correlation function
inv_fourier( num, x, Fx );
/*----------------------------------------------------------------------*/
// 相互相関関数をデータから直接計算
// x[n].real: ウィーナー・ヒンチンの定理から求めた相互相関関数
// Cxy: データから直接求めた相互相関関数
// calculate cross-correlation function directly from the data
// x[n].real: cross-correlation function obtained by Wiener-Khinchin theorem
// Cxy: cross-correlation function obtained directly from the data
printf("#---------------------------------------------------------------------------\n");
printf("# ウィーナー・ヒンチンの定理\n");
printf("#---------------------------------------------------------------------------\n");
printf("# パワースペクトルを逆フーリエ変換し、相互相関関数を計算\n");
printf("# WKT CCF: パワースペクトルから求めた相互相関関数（ウィーナー・ヒンチンの定理）\n");
printf("# DATA CCF: データから直接求めた相互相関関数\n");
printf("# gap: WKT CCF - DATA CCF\n");
printf("#---------------------------------------------------------------------------\n");
printf("#---------------------------------------------------------------------------\n");
printf("# Wiener-Khinchin theorem (using cyclic data)\n");
printf("#---------------------------------------------------------------------------\n");
printf("# inverse Fourier transform of the power spectrum and calculation of the cross-correlation function\n");
printf("# WKT CCF: cross-correlation function obtained from the power spectrum (Wiener-Khinchin theorem)\n");
printf("# DATA CCF: cross-correlation function obtained directly from the data\n");
printf("# gap: WKT CCF - DATA CCF\n");
printf("#---------------------------------------------------------------------------\n");
printf("# lag WKT CCF DATA CCF gap\n");
for( n = 0; n < num; n++ ){
sum_Cxy = 0.0;
for( i = 0; i < num; i++ ){
sum_Cxy += ( xx[n+i] * yy[i] );
}
Cxy[n] = sum_Cxy; /* ### */
}
for( n = 0; n < num; n++ ){
diffx = x[n].real - Cxy[n];
printf(" %d %lf %lf %lf\n", n, x[n].real, Cxy[n], diffx );
}
printf("#---------------------------------------------------------------------------\n");
/*------------------------------------------------------------------------------------------
相互相関関数をフーリエ変換し、パワースペクトルを計算
Fourier transform the cross-correlation function and calculate the power spectrum
------------------------------------------------------------------------------------------*/
// ### で求めた相互相関関数を利用
// use the cross-correlation function obtained in ###
for( n = 0; n < num; n++ ){
x[n].real = Cxy[n];
x[n].imag = 0.0;
}
// 離散フーリエ変換
// データxを入力し、位相Fを出力
// x[]: 相互相関関数
// Fx[].real: パワースペクトル
// Fx[].imag: パワースペクトル
// discrete Fourier transform
// input data x and output phase F
// x is the cross-correlation function, and the F is the power spectrum
fourier( num, x, Fx );
printf("# 相互相関関数をフーリエ変換し、パワースペクトルを計算\n");
printf("# DATA PS: データから直接求めたパワースペクトル\n");
printf("# WKT PS: 相互相関関数から求めたパワースペクトル（ウィーナー・ヒンチンの定理）\n");
printf("# gap: DATA PS - WKT PS\n");
printf("#---------------------------------------------------------------------------\n");
printf("# Fourier transform the cross-correlation 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 cross-correlation function (Wiener-Khinchin theorem)\n");
printf("# gap: DATA PS - WKT PS\n");
printf("#---------------------------------------------------------------------------\n");
printf("# freq DATA PS (real) WKT PS (real) gap ((real) DATA PS (imag) WKT PS (imag) gap (imag)\n");
for( n = 0; n < num; n++ ) {
diffx = power_spec[n].real - Fx[n].real;
diffy = power_spec[n].imag - Fx[n].imag;
printf(" %d %lf %lf %lf %lf %lf %lf\n", n, power_spec[n].real, Fx[n].real, diffx, power_spec[n].imag, Fx[n].imag, diffy );
}
/*----------------------------------------------------------------------*/
free( x );
free( y );
free( xx );
free( yy );
free( Cxy );
free( Fx );
free( Fy );
free ( power_spec );
return 0;
}