フラクショナル積分の離散定義

プログラムリスト

#define MAXKAISUU s //積分演算をおこなう時に利用する過去の情報の最大数
#define m_chk_mean true //最終的に平均値を零とする処理をおこなうかどうか
#define m_DT 1.0 //時系列のサンプリング時間間隔

//入力信号をそのまま微分する
//*data・・・微分をおこなう対象データのポインタ
//sekibun・・・積分階数
//s・・・データの大きさ
void CFractionalDlg::FractionalIntegral(double *ret_data, double sekibun, int s)
{
	register int m,t;

	//sekibunが整数かどうか判別する（整数の時には動作が異なる）
	int kaisuu=(sekibun==(int)sekibun)? (int)sekibun: MAXKAISUU;
	double* ganma=new double[kaisuu+1];

	//ガンマ関数の値をストックする    (n+1)Ｃ(m+1)
	//ついでに符号も保存する
	ganma[0] = 1.0;
	for(m = 1; m <= kaisuu; m++)
		ganma[m] = - ganma[m - 1] * (sekibun - m + 1) / (double)m;

	//非整数階積分をする
	double p = pow(m_DT, sekibun); //h^d
	// for(t = s; t; data[--t] -= data[0]);  //積分の場合はいらない

	register double f;
	int M;
	for(t = 0; t < s; t++){
		M = (t < kaisuu)? t: kaisuu;
		for(f = data[t] * p, m = 1; m <= M; f -= ganma[m] * data[t - m], m++);
		data[t] = f;
	}
	delete [] ganma;

	//データの平均を０に
	if(m_chk_mean){
		for(f = 0.0, t = s; t; f += data[--t]);
		for(f /= s, t = s; t; data[--t] -= f);
	}
}