/*--------------------------------------------------------------------------------
フーリエ・トランスフォーム・サロゲート法
Fourier transform (FT) surrogate method
FT_surrogate.c
*--------------------------------------------------------------------------------*/
#include "FT_surrogate.h"
int main( void )
{
FILE *fo;
char filename[20];
int i, k;
int num;
unsigned long int seed = MT_SEED;
double temp;
double pow_spec_surrogate;
double *pow_spec;
complex *x;
complex *F;
complex *xx; /* random shuffle surrogate data */
// x.real: 実部, real part
// x.imag: 虚部, imaginary part
/*--------------------------------------------------------------------------*/
// 乱数の初期設定
// initialisation of random numbers
init_genrand(seed);
/*--------------------------------------------------------------------------*/
if ( ( x = readdata( &num ) ) == NULL ) {
fputs( "readdata error\n", stderr );
return -1;
}
/*--------------------------------------------------------------------------*/
if (( xx = (complex *)malloc( num * sizeof(complex) )) == NULL ) {
fputs( "cannot allocate memory for F\n", stderr );
return -1;
}
if (( F = (complex *)malloc( num * sizeof(complex) )) == NULL ) {
fputs( "cannot allocate memory for F\n", stderr );
return -1;
}
if (( pow_spec = (double*)malloc((num) * sizeof(double))) == NULL ) {
fputs( "pow_spec cannot allocate memory\n", stderr );
return -1;
}
/*--------------------------------------------------------------------------*/
// 離散フーリエ変換
// データxを入力し、位相をFに出力
// discrete Fourier transform
// input the data x, output the phase F
fourier( num, x, F );
// 入力データのパワースペクトルの計算
// calculation of the power spectrum of the input data x
for( i = 0; i < num; i++ ){
pow_spec[i] = F[i].real * F[i].real + F[i].imag * F[i].imag;
// printf("# %lf %lf %.10lf\n", F[i].real, F[i].imag, pow_spec[i]);
}
/*==========================================================================*/
for( k = 0; k < SURROGATE_DATA_NUM; k++ ){
// ランダムな位相を得るため、RSサロゲートデータの作成
// genarate RS surrogate data to obtain random phases
/*------------------------------------------------------------------
オリジナルデータが入っているxを入力。
RSサロゲートデータをxxに格納。
xには何の操作も行わないので、入力した状態まま。
input the original data x
RS surrogate data is stored in xx.
As no operation is performed on x, it remains as entered.
------------------------------------------------------------------*/
if ( random_shuffle( num, x, xx ) != 0 ) {
fputs( "cannot generate random shuffle sorrogate data\n", stderr );
return -1;
}
/*--------------------------------------------------------------------------*/
// 離散フーリエ変換
// RSサロゲートデータxxを入力し、位相をFに出力
// discrete Fourier transform
// input the RS surrogate data xx, output the phase F
fourier( num, xx, F );
/*----------------------------------------------------------------------
RSサロゲートデータをフーリエ変換してランダムな位相を得た。
その位相のパワースペクトルが、オリジナルデータのパワースペクトル
と同じになるように変換する。
As the RS surrogate data was Fourier transformed, random phases
are obtained.
The power spectrum of that phases is transformed so that it is
the same as the power spectrum of the original data.
----------------------------------------------------------------------*/
for( i = 0; i < num; i++ ){
pow_spec_surrogate = F[i].real * F[i].real + F[i].imag * F[i].imag;
temp = sqrt( pow_spec[i]/(double)pow_spec_surrogate );
F[i].real = temp * F[i].real;
F[i].imag = temp * F[i].imag;
}
/*--------------------------------------------------------------------------*/
// 離散逆フーリエ変換
// 位相Fを入力し、データxxを出力
// データxxがFTサロゲートデータとなる。
// inverse discrete Fourier transform
// input the phase F, and output the data xx
// data xx becomes FT surrogate data.
inv_fourier( num, xx, F );
/*--------------------------------------------------------------------------*/
printf("# Now, generating the %dth FT surrogate data\n", k+1 );
if( k < 9 ){
sprintf(filename, "FT0%1d.dat", k+1);
}
if( k >= 9 ){
sprintf(filename, "FT%1d.dat", k+1);
}
if ((fo = fopen(filename,"w")) == NULL){
printf("file open error\n");
exit(-1);
}
fprintf(fo, "%d\n", num);
for( i = 0; i < num; i++ ){
// fprintf(fo, "%.12e\n", xx[i].real);
fprintf(fo, "%.10lf\n", xx[i].real);
// fprintf(fo, "%.10lf %.10lf\n", xx[i].real, xx[i].imag);
}
fclose(fo);
}
free( x );
free( xx );
free( F );
free( pow_spec );
return 0;
}
/*--------------------------------------------------------------------------*/
int random_shuffle( int num, complex *org_data, complex *RS )
{
int i;
int kk;
int *order;
unsigned long int j;
unsigned long int r;
if(( order = (int*)malloc((num) * sizeof(int))) == NULL ) {
fputs("cannot allocate memory for order\n", stderr );
return -1;
}
for( i = 0; i < num; i++ ) {
order[i] = i;
}
j = num;
for( i = 0; i < num; i++ ){
// 正の整数の一様分布, Uniform distribution of positive integers
r = genrand_int32();
r %= j;
kk = order[r];
order[r] = order[j-1];
j--;
RS[i] = org_data[kk];
// printf("%lf %lf\n", RS[i].real, RS[i].imag);
}
free( order );
return 0;
}