#define TAU 1 // 埋め込みの遅れ時間
#define N 10000 // データの総点数
#define MAXJSPHERE (2*10*1000) // 相関積分の配列の要素数
// 2は2乗距離，10はディケード，1000は分割数を表す
#define chk_fast false // 高速化のためランダムサンプリングする
double m_data[N]; //１次元の時系列データ
//指定した次元までの相関積分を計算する
// Dimension：埋め込み次元
void CFractalDim::RunAnalysisCorrelationMS( int Dimension )
{
register int i, j, m, l;
register double f, g, p;
int ltaumax = (Dimension - 1) * TAU; //埋め込みに用いる過去のデータ点数
int n = N - ltaumax; //実際に埋め込むアトラクタの点数
data = &m_data[ltaumax]; //相関積分の計算に用いるデータの先頭番地
//m_data[0]〜m_data[ltaumax-1]は埋め込みに用いる
//データを正規化するための定数lengthを計算する
//データの最大最小値を求める
f = g = m_data[0];
for(i = 0; i < N; i++){
if(f < m_data[i]) f = m_data[i];
if(g > m_data[i]) g = m_data[i];
}
//以下の1024は適当に変更してよい
double length = f - g;
length = 1024.0 / length, length = length * length; //２乗距離に変換
//各次元における相関積分の値を格納する配列の定義
_int64** atrsuu = new _int64* [Dimension];
for(m = 0; m < Dimension; m++){
atrsuu[m] = new _int64[MAXJSPHERE + 1];
for(j = 0; j <= MAXJSPHERE; j++){
atrsuu[m][j] = 0; //初期化
}
}
//相関積分のインデックスでの最大半径の２乗を計算する
double lengthindex[MAXJSPHERE + 1];
//int z = (int)ceil((log10( x * length )) * 1000.0) の逆関数
for(j = 0; j <= MAXJSPHERE; j++){
lengthindex[ j ]=pow(2.0, j * 0.001) / length;
}
//ソーティングのためのクラス
CSort sort;
double* sortbox = new double [n * Dimension];
//高速化のため，ランダムサンプリングするかどうか？
if(chk_fast){
//ランダムサンプリングする場合
int k;
//１００分の１のデータを利用する
for(int loop=0; loop < n/100; loop++){
//中心点を決定する
i = (int)((n / (RAND_MAX + 1.0)) * (rand() + 0.5));
//各次元における，中心点からの２乗距離を計算する
for(j = 0; j < n; j++){
for(g = 0.0, m = l = 0; m < Dimension; l += TAU, m++){
f = data[i-l] - data[j-l], g += f * f, sortbox[m * n + j] = g;
}
}
//ヒストグラムの計算
for(m = 0; m < Dimension; m++){
//距離分布をクイックソートする(データはn個 [0, n-1]の範囲)
sort.QuickSortRiseD(&sortbox[m * n], 0, n-1);
//距離の短いものからヒストグラムのカウントをおこなう
for(k = j = 0; j <= MAXJSPHERE; j++){
for(f = lengthindex[j]; k < n && sortbox[m * n + k] <= f; atrsuu[m][j]++, k++);
}
//相関次元の場合は自分の点は計算しない
atrsuu[m][0]--; //自分の点を除外する
}
}
p = (double)(n - 1) * (double)(n / 100); //トータルの点の数
}else{
//ランダムサンプリングしない場合
int k;
double* len = new double[N];
double* plen = &len[ltaumax];
// n * nの相関積分マップの作成
// 条件 N=n+ltaumax , ltaumax=(Dimension-1)*tau;
for(i = 1; i < n; i++){
// 距離テーブルを対角線上に計算する
for(j = 0; j < N - i; f = m_data[i + j] - m_data[j], len[j] = f * f, j++);
//各次元の埋め込みによる距離関数を求める
for(j = 0; j < n-i; j++){
for(l = j, g = 0.0, m = 0; m < Dimension; l -= TAU, m++){
g += plen[l], sortbox[m * n + j]=g;
}
}
//距離関数を昇順ソートし，相関積分のインクリメントをする
for(m = 0; m < Dimension; m++){
sort.QuickSortRiseD(&sortbox[m * n], 0, n-i-1);
for(k = 0,j = 0;j <= MAXJSPHERE; j++){
for(f = lengthindex[j]; k < n-i && sortbox[m * n + k] <= f; atrsuu[m][j]++, k++);
}
}
}
delete [] len;
p = (double)(n - 1) * (double)n / 2.0;
}
delete [] sortbox;
//相関関数の積分
for(m = 0; m < Dimension; m++){
for(j = 0; j < MAXJSPHERE; j++){
atrsuu[m][j + 1] += atrsuu[m][j];
}
}
//ファイルへの出力（Ｅｘｃｅｌで開ける形式ってことで）
FILE* fp;
fp=fopen("corre.dat","w");
fprintf(fp, "↓log(ε)\t→Nm(ε)\tN=%d\tm=%d\tτ=%d\n", N, Dimension, TAU);
for(j = 0; j <= MAXJSPHERE; j+=10){
fprintf(fp,"%e", j / 1000.0 / 2.0);
for(m = 0; m < Dimension; m++){
if(atrsuu[m][j] != 0)
fprintf(fp,"\t%e",log10((double)atrsuu[m][j] / p) / log10(2.0));
else
fprintf(fp,"\t");
}
fprintf(fp,"\n");
}
fclose(fp);
//メモリの開放
for(m = 0; m < Dimension; m++){
delete [] atrsuu[m];
}
delete [] atrsuu;
}
言語は一応C++です．（C言語も混ざっていますが）．
相関積分の計算量
