アルゴリズム

　自己アフィン的な時系列データを f ( t ) とするとき，そのスケーリング特性は，
(1)
と表される．ここで H はHurst指数である．台となる空間の次元を d とすると，そのフラクタル次元 D は次式で定義される．
(2)
　ここでは，時系列データ f ( t ) のパワースペクトル密度 S ( w ) からフラクタル次元 D _s を求める手法について述べる．
　f ( t )のフーリエ変換を F ( w ) とすると， S ( w ) は，
(3)
と表現される．ここで，< > はアンサンブル平均を表し，R ( t ) は自己相関関数で，そのスケーリング特性は，
(4)
となる．従って，式(3), (4)より，
(5)
が導かれ，式(2)より次式が得られる．
(6)
しかしながら，実際のパワースペクトル S ( w ) は図１のように，スペクトル成分が広範囲に振動するため， スケーリング指数 b を求めることは困難である． そこで我々は，周波数領域でのモーメント指数を導入し，その臨界特性を利用して，ベキ特性を抽出する．
　モーメント指数を a とすると，モーメント I _a は，
(7)
と定義される．ここで，式(7)の性質を考慮すると，
(8)
と表現できる．ここで， U = a - b + 1 である．さらに，log I _a の a による２階微分は，
(9)
となる．この式は，W → ∞ の極限では，U = 0 において発散するため，更に３階微分，
(10)
を計算し，これが0となる臨界モーメント指数 a _c を求める．ここで， I _a の n 階微分 I _a ^{( n)} は次式で与えられる．
(11)
このとき，U = a _c - b + 1 = 0より，臨界モーメント指数法によるフラクタル次元 D _M は次式で求められる．
(12)
　臨界モーメント指数法による a _c ，および，D _M の推定結果の一例を図２に示す．また， 母音データ{ /a/, /i/, /u/, /e/, /o/}に対する推定結果を表１に示す．

図１：母音/a/のパワースペクトル密度の一例


図２：母音/a/に対するモーメント指数の臨界特性の一例


表１：母音データ{ /a/, /i/, /u/, /e/, /o/}に対するフラクタル次元の推定結果

プログラムリスト

//臨界指数法処理ルーチン
#define max_arpha 20 //αの最大値・最小値 [α,-α]
#define step_arpha 400 //αの変化の刻み
double CFractalDim::CEMArphaC(double *x, double *y, int index)
{
	int j, loop;
	double* div = new double[step_arpha + 1];					//各ループの３階微分の値
	double* arpha = new double[step_arpha + 1];				//各ループのαの値
	double i_arpha, i_arpha1, i_arpha2, i_arpha3;
	double X, Y, f, g;

	for(loop = 0; loop <= step_arpha; loop++)					//各ループのαの値を求める
		arpha[loop] = max_arpha * (2.0 * loop / step_arpha - 1.0);

	for(loop = 0; loop <= step_arpha; loop++){
		i_arpha = i_arpha1 = i_arpha2 = i_arpha3 = 0.0;
		for(j = index - 2; j >= 0; j--){
			X = x[j] / x[0], Y = x[j + 1] / x[0],
			f = 0.5 * (Y - X) * y[j    ]* pow(X, arpha[loop]),
			g = 0.5 * (Y - X) * y[j + 1]* pow(Y, arpha[loop]),
			X = log(X), Y = log(Y), i_arpha  += (f + g),
			f *= X, g *= Y,         i_arpha1 += (f + g),
			f *= X, g *= Y,         i_arpha2 += (f + g),
			f *= X, g *= Y,         i_arpha3 += (f + g);
		}
		i_arpha1 /= i_arpha,
		i_arpha2 /= i_arpha,
		i_arpha3 /= i_arpha,
		div[loop] = i_arpha3 - 3.0 * i_arpha2 * i_arpha1 + 2.0 * i_arpha1 * i_arpha1 * i_arpha1;
	}

	double max = 0.0;
	double pos = 0, DIV, xp, arphac;
	bool check = true;

	for(loop = 0; loop < step_arpha - 1; loop++){
		if(div[loop] > max) max = div[loop], pos = loop;
	}

	//　0交差点を探索する
	for(loop = pos; loop < step_arpha - 1; loop++){
		if(div[loop] > 0.0 && div[loop + 1] >= 0.0){
			DIV = (div[loop + 1] - div[loop]) / (2.0 * max_arpha / step_arpha);
			pos = div[loop], xp = arpha[loop], f = arphac = -pos / DIV + xp;
			check = false;
			break;
		}
	}
	delete [] arpha;
	delete [] div;

	if(check){
		AfxMessageBox("自動次元推定ができませんでした");
		return 0.0;
	}

	return arphac;  //beta=a+1
}

参考文献