正規分布に従う乱数を分数階積分してフラクショナルブラウン運動を生成し，csv形式に出力するプログラム．
下記のC++11のコンパイラで動作を確認
- Embarcadero Technologies Inc. bcc32c version 3.3.1
- Microsoft Visual Studio Community 2015 Version 14.0.25431.01 Update 3
//>>bcc32c fBm_generator.cpp //>>fBm_generator.exe //または //>>fBm_generator.exe [Hurst指数] #include<iostream> #include<math.h> #include <random> #include <fstream> #include <stdlib.h> ////////設定項目//////// //Hurst指数 下記の値または，コマンドライン引数の一番目 #define HURST_EXPONENT (0.5) //出力する時系列データの点数 #define DATALEN (65536) //fBmの時系列データを出力するファイル名．ただし，出力形式はcsv #define fBm_FILENAME ("fBm.csv") //fBmの生成に使用する乱数の時系列データを出力するファイル名．ただし，出力形式はcsv #define WHITE_NOISE_FILENAME ("white.csv") ////////計算プログラム//////// //分数階積分する関数"FractionalIntegral"は下記のサイトから引用，一部改変した． //http://www.geocities.co.jp/Technopolis-Mars/1795/study/fractint.html #define MAXKAISUU s //積分演算をおこなう時に利用する過去の情報の最大数 #define m_chk_mean false //最終的に平均値を零とする処理をおこなうかどうか #define m_DT 1.0 //時系列のサンプリング時間間隔 //入力信号をそのまま微分する //*data・・・微分をおこなう対象データのポインタ //sekibun・・・積分階数 //s・・・データの大きさ void FractionalIntegral(double *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); } } //正規分布に従う乱数を発生させる関数"randnormal"は下記のサイトを引用，改変した． //https://cpprefjp.github.io/reference/random.html double randnormal(double mean, double stdDev) { // メルセンヌ・ツイスター法による擬似乱数生成器を、 // ハードウェア乱数をシードにして初期化 std::random_device seed_gen; std::mt19937 engine(seed_gen()); // 正規分布 // 平均mean、標準偏差stdDevで分布させる std::normal_distribution<> dist(mean, stdDev); return dist(engine); } int main(int argc, char *argv[]) { double signal[DATALEN]; //正規分布に従う乱数を発生 for (int i = 0; i < DATALEN; ++i) { signal[i] = randnormal(0, 1); } std::cout << "乱数の一番目の値：" << signal[0] << std::endl; //正規分布に従う乱数を保存 std::ofstream filew(WHITE_NOISE_FILENAME); for (int i = 0; i < DATALEN; ++i) { filew << signal[i] << "\n"; } //Hurst指数に0.5を加算し，乱数時系列信号を分数階積分する階数integrateOrderを決定 double integrateOrder = HURST_EXPONENT + 0.5; if (argc > 1) { char *strtod_error; double input = strtod(argv[1], &strtod_error); if (*strtod_error != '\0') { std::cout << "Hurst指数として設定できない文字列：" << strtod_error << std::endl; } integrateOrder = input + 0.5; } std::cout << "Hurst指数：" << integrateOrder - 0.5 << std::endl; //分数階積分を実行 FractionalIntegral(signal, integrateOrder, DATALEN); //フラクショナルブラウン運動信号を保存 std::ofstream file(fBm_FILENAME); for (int i = 0; i < DATALEN; ++i) { file << signal[i] << "\n"; } file.close(); filew.close(); return 0; }
