nlab study

非整数階微積分演算の高速化について

畳み込み演算による表記

　関数 f( t ) の q 階微分（-q 階積分）の計算式は，

と表される．ここでΓ関数による係数を関数 g( m ) と置き換えると，

が得られる．この式は明らかに f( t ) と g( m ) との畳み込み演算であり，f( t ) の複素スペクトルを F( k ) ， g( m ) の複素スペクトルを G( k ) とすると，非整数階の微積分演算後の関数 x( t ) は次式で与えられる．

ここで，Ｆ^-1[ ] は周波数変数 k の逆フーリエ変換を表す．
　実際の計算では３度のフーリエ変換として高速フーリエ変換(FFT)アルゴリズムが利用される．FFTのアルゴリズムについては詳細は述べないが， FFTによる非巡回的畳み込みの計算手順は次のようになる．

　１．　 g( m ) （ m = 0,1,･･･, N -1)，f( t ) （ t = 0, h ,･･･,( N -1) h )を準備する．
　２．　畳み込みの値 x( t ) ( t = 0, h ,･･･,( L -1) h )を計算するものとする．
　　　　ここで， L は 2N - 1 ≦ L を満足するように選ぶ．ただし， L = 2 ^integer とする．
　３．　係数 g( m ) とデータ f( t ) の後ろにゼロを付け加え，データの個数が L 個となるようにする．
　４．　係数 g( m ) とデータ f( t ) をデータ数 L のFFT演算を施し，複素スペクトル G( K ) とデータ F( k ) を求める．
　５．　周波数 k と K の単位が同じ（ k = K ）であると扱い，積 G( k ) F( k ) を求める．
　６．　５．の結果に逆FFTを施し 1/h^q を乗算して畳み込みの値 x( t ) を得る．

高速化の効果

　今までフラクショナル（非整数階）演算の悩みは時間がかかることであったが，この畳み込み演算の発想とFFTの利用によって処理時間に歴然とした差が現れた．もっと前から考えていたのだが，こんなに速くなるとは思わなかった．これで，非整数階微積分を利用した最尤推定法によるフラクタル次元の推定手法の優位性が向上した．
この高速演算処理はもっと他の処理にも応用できそうである．

プログラムリスト

上記の文字に対応していない場合があります．コメントを参照してください．


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

//*org_data・・・微積分をおこなう対象データのポインタ
//*ret_data・・・微積分演算後のデータのポインタ
//bisekibun・・・微積分階数（正：微分，負：積分）
//s・・・データの大きさ
//微積分する(厳密な処理を省いたバージョン)
void CFractionalDlg::FractionalDI_FFT(double*org_data, double *ret_data, double bisekibun, int s)
{
	register int t,kaisuu;
	register double f;

	//bisekibunが整数かどうか判別する（整数の時には動作が異なる）
	kaisuu=(bisekibun<0.0)? MAXKAISUU: ((bisekibun==(double)((unsigned int)bisekibun))? (unsigned int)bisekibun+1: MAXKAISUU);

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

	//データの振幅をシフトするかどうか
	f=org_data[0];                                         //シフト量
	if(f!=0.0 && bisekibun>=0.0){
		for(t=s; t; --t,ret_data[t]=org_data[t]-f);
	}else{
		for(t=s; t; --t,ret_data[t]=org_data[t]);
	}

	for(int L=2; L < s+kaisuu-1; L<<=1); //FFTのデータ点数の決定
	
	//Γ関数のフーリエ変換をする
	CFft* fft=new CFft(L);
	fft->CalcFFTsignal2(org_data,s);               //先にorg_dataのFFTを計算
	fft->CalcConvolution(ganma,kaisuu,ret_data,s); //ganmaのFFTと畳み込みの演算
	delete fft;

	double p=pow(m_DT,-bisekibun); //h^q  
	for(t=s; t ; --t,ret_data[t]*=p);
	if(bisekibun>0.0) ret_data[0]=0; //直流成分は0にする(この値だけはうまく処理できないため)

	if(m_chk_mean){
		for(f=0.0,t=s;t;) f+=ret_data[--t];
		for(f/=s,t=s;t;) ret_data[--t]-=f;
	}

	delete [] ganma;
}

FFTの演算のためのクラス


#include <math.h>
#include <complex>
using namespace std; //complex形を利用するのに必要

class CFft  
{
public:
	void LessErrorSam(double &sam, double add, double &error);
	void SetPast0(double *x,int suu); //実数データの収納
	void fft(bool inverse=false);
	void CalcFFTsignal2(double *realdata, int suu_real);
	void CalcConvolution(double* realdata,int suu_real,double* ret_data,int suu_ret);

	//取り扱いデータ
	double* sintbl; //サイン関数のテーブル（予め計算する）
	int* bitrev;    //符号反転テーブル
	complex<double> uj; //虚数
	complex<double> zero; //零
	complex<double> arg; //最小の回転角
	complex<double>* signal; //取り扱う複素データ
	complex<double>* signal2; //取り扱う複素データ
	int N,N2,N4,N8; //データ点数

	CFft(int n);	//コンストラクタ
	virtual ~CFft();	//デコンストラクタ

private:
	void make_bitrev();
	void make_sintbl();
};

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
#define PAI 3.14159265358979323846
CFft::CFft(int n) // n : FFTのデータ点数
{
	//FFTではnは必ず2の累乗でなくてはならない（DFTではOK）
	N=n,N2=n>>1,N4=n>>2,N8=n>>3; //nをどんどん割っていく

	bitrev=new int [n]; //符号反転テーブル
	make_bitrev();

	sintbl = new double [N + N4]; //サイン関数のテーブル
	make_sintbl();

	uj  =complex<double>(0.0,1.0); //虚数ベクトル
	zero=complex<double>(0.0,0.0); //ゼロベクトル

	arg= -uj*PAI*2.0/(double)n;    //最小の変化角度

	signal=new complex<double> [n]; //メインデータ
	signal2=NULL;                   //サブデータ
}

CFft::~CFft()
{
	delete [] bitrev;
	delete [] sintbl;
	delete [] signal;
	if(signal2!=NULL) delete [] signal2;
}

//signal2にデータのスペクトルが格納する
void CFft::CalcFFTsignal2(double *realdata, int suu_real)
{
	SetPast0(realdata,suu_real); 			//suu*2<=Nである必要がある
	if(signal2==NULL) delete [] signal2;
	signal2=new complex<double>[N];		//データ領域の確保
	fft(false);					//フーリエ変換
	for(int m=N; m ; --m,signal2[m]=signal[m]);	//データの格納
}

//signal2に元データのスペクトルが格納されているとき
//realdataとの畳み込みを計算する
//suu_realはrealdataの数
//suu_retはret_dataの個数
void CFft::CalcConvolution(double *realdata, int suu_real, double *ret_data, int suu_ret)
{
	SetPast0(realdata,suu_real);			//suu*2<=Nである必要がある
	fft(false);					//フーリエ変換
	for(int m=N; m ; --m,signal[m]*=signal2[m]);	//元データのスペクトルを掛ける
	fft(true);					//逆フーリエ変換
	for(m=suu_ret; m ; --m,ret_data[m]=signal[m].real()); //データの格納
}

/**************************************************************************************
予め計算しておく三角関数／ビット反転テーブルの作成
**************************************************************************************/
#define FFT_USE_SIN_FUNCTION //誤差を抑えるため，サインテーブルの計算に関数をそのまま使う
void CFft::make_sintbl()
{
	int i;
	register double t;

	#ifdef FFT_USE_SIN_FUNCTION
		for(i=N8;i;){ //後々，アセンブラを使って高速演算に書き直す予定
			--i,t=PAI*i/(double)N2,
			sintbl[   i]=sin(t),
			sintbl[N4-i]=cos(t);
		}
	#else
		register double c,s,dc,ds,t;
		register double e_c,e_s,e_dc,e_ds;

		t=sin(PAI/N);
		dc=2*t*t;
		ds=sqrt(dc*(2-dc));
		t=2*dc;
		c=sintbl[N4]=1.0;
		s=sintbl[ 0]=0.0;
		e_c=e_s=e_dc=e_ds=0.0;
		for(i=1;i<N8;i++){ //誤差を抑えた加算を行う
			LessErrorSam( c, -dc,e_c ); //c-=dc
			LessErrorSam(dc, t*c,e_dc); //dc+=t*c
			LessErrorSam( s,  ds,e_s ); //s+=ds
			LessErrorSam(ds,-t*s,e_ds); //ds-=t*s
			sintbl[   i]=s;
			sintbl[N4-i]=c;
		}
	#endif

	if (N8!=0) sintbl[N8]=sqrt(0.5);
	for(i=N4;i;)	--i,sintbl[N2-i]=sintbl[i];
	for(i=N2+N4;i;)	--i,sintbl[i+N2]=-sintbl[i];
}

void CFft::LessErrorSam(double &sam, double add, double &error)
{
	register double f;
	error+=add,		//積み残しにノルムを加える
	f=sam,			//前回までの和
	sam+=error,		//和を更新
	f-=sam,			//実際に積まれた値の符号反転
	error+=f;		//積み残し
}

//ビット反転テーブルの作成
void CFft::make_bitrev()
{
	int i=0,j=0,k;
	for(;;){
		bitrev[i]=j;
		if(++i >= N) break;
		k=N2;
		while(k<=j){ j-=k,k>>=1; }
		j+=k;
	}
}

/*****************************************************************************************
 FFTメインルーチン
 inverse==true のときには逆フーリエ変換である
******************************************************************************************/
void CFft::fft(bool inverse)
{
	int i,j,k,h,d,k2;
	complex<double> cs,dxy;

	for(i=0;i<N;i++)
		if(i<(j=bitrev[i])) cs=signal[i] , signal[i]=signal[j] , signal[j]=cs;

	if(inverse){ //逆フーリエ変換
		for(d=N, k=1;k<N;k=k2){
			k2=k << 1, d >>= 1;
			for(h=j=0; j<k; j++,h+=d){
				cs=complex<double>(sintbl[h+N4], sintbl[h]);
				for(i=j;i<N;i+=k2)
					dxy=cs*signal[i+k], signal[i+k]=signal[i]-dxy, signal[i]+=dxy;
			}
		}
		for(i=N;i;) signal[--i]/=N;
	}else{
		for(d=N, k=1;k<N;k=k2){
			k2=k << 1, d >>= 1;
			for(h=j=0; j<k; j++,h+=d){
				cs=complex<double>(sintbl[h+N4],-sintbl[h]);
				for(i=j;i<N;i+=k2)
					dxy=cs*signal[i+k], signal[i+k]=signal[i]-dxy, signal[i]+=dxy;
			}
		}
	}
}

言語は一応C++です．（C言語も混ざっていますが）．
バグがありましたら連絡してください．お願いします．

⇒非整数階微分
⇒非整数階積分
⇒非整数階微積分の統一表現

戻る

2003/6
by Hiroyuki Koga