背景
１次元、２次元の空間的なフーリエ空間で対数正規分布(log-normal)するデータを生成する python スクリプトを生成した。宇宙流体シミュレーションソフトウェア pluto (http://plutocode.ph.unito.it) への出力ファイルまで作成してくれる格好にした2Dの対数分布を生成する。
一次元の対数正規分布 log-normal するデータ生成方法
スクリプト
- http://www-x.phys.se.tmu.ac.jp/%7Esyamada/syamadatools/others/genlc_lognormal/GenRandomLCLogNormal_1D.py (python2系)
- https://colab.research.google.com/drive/1LniT5AF2tgtCUwXELtyqXYv2VaUX2ssl?usp=sharing
使い方
chmod +x ファイル で実行権限を与えてから、
./GenRandomLCLogNormal_1D.py
-h でヘルプを表示すると使い方が表示される。(現状python2系)
google Colab 上であれば、動作確認できます。
ドキュメント
単純にパワースペクトルから位相ランダムで生成したい場合はこちらを参照ください。
- https://qiita.com/yamadasuzaku/items/169d861c0a66c847a639
注意
非線形な領域、exp(x) = 1 + x が成り立つ範囲であれば、exp() しても、
パワースペクトルは変わらない。非線形性がでる範囲であれば、注意する
必要がある。exp(X) が取れる = X は無次元の量に規格化されているはず。その量が、
1に比して、大きいかどうかで、変わる。例えば、数10%の変動であれば、パワースペクトルに変化はないが、factor の変化や、桁で変化する場合は、パワースペクトルにも影響がでるので、ダイナミックレンジの問題とも言える。
２次元の対数正規分布 log-normal するデータ生成方法
スクリプト
使い方
chmod +x ファイル で実行権限を与えてから、
./GenRandomLCLogNormal_2D.py
-h でヘルプを表示すると使い方が表示される。例えば、
./GenRandomLCLogNormal_2D.py -k -2.5
とすると、べきを -2.5 に変更できます。
ドキュメント
補足事項
少し調整した箇所は、
cor_ifft_real = ifft_real/np.std(ifft_real) # just to suppress large values before exp()
print "..... divided by np.std(ifft_real)", np.std(ifft_real)
の箇所。exp()を取る前に、値をそこそこ小さく規格化してないと、exp をとった際に、値が激しく飛ぶことがあります。それで、強い理由はないが、np.std で、expを取る前の値の標準偏差で割ってます。grid 数は、512でベタ書きで hard cording した。optparse で変更できるようにする予定。nmax という変数。
結果の例
１次元の場合
4ks の時間で、時間分解能 0.2秒、カットオフ周波数 0.1 Hz で、べきが -1 から -2 に変化するパワースペクトルから生成した時系列データの例を示す。
まずは、入力のチェックをしている。
自分の入力したパワースペクトルを逆変換して、もとに戻ることを確認しているだけ。
黒のライトカーブが対数正規分布する一次元データになる。
２次元の場合
[syamada] $ ./GenRandomLCLogNormal_2D.py -p [~/work/software/sxs/sxspsp-20141229/others/genlc_lognormal]
---------------------------------------------
| STATUS : Start ./GenRandomLCLogNormal_2D.py
============== input parameters ====================
---- Index1 of input FFT = -1.0
---- Norm of input FFT = 1.0
---- Flag plotting = True
---- Output Filename = lognorm2d
---- mean = 0.0
---- sigma = 1.0
---- outfname = lognorm2d_index1-1p0
==========================================================
というパラメータで動作させた場合の例。
参考文献
イタレーションをしてないこと以外は、
- https://ams.confex.com/ams/11AR11CP/techprogram/paper_42772.htm
と同じ実装方法