2019-02-18

pythonでレスラー方程式の解軌道を描画（ルンゲクッタ法）

pythonでルンゲクッタ法を実装した．
レスラー方程式を与えて描画する．

プログラム内関数のdxdt, dydt, dzdtを書き換えれば別の微分方程式も解ける．
下記プログラムはローレンツ方程式を解くようにも書き換えられる．

import numpy as np


def dxdt(a, b, c, x, y, z):  # 引数は，「関数のパラメータを並べたもの」，「変数を並べたもの」の順番で並べる．他の連立微分方程式と同じ形式で引数を並べる．
    return -y - z


def dydt(a, b, c, x, y, z):
    return x + a * y


def dzdt(a, b, c, x, y, z):
    return b + x * z - c * z


def _handler(func, *args):
    return func(*args)


pdelist = [dxdt, dydt, dzdt]  # 連立微分方程式を追加する場合は，関数名をこのlistにも追加すること．


def runge_kutta(funcs, param, vold, delta):
    """
    ルンゲクッタ法で時間を含まない微分方程式を解く
    :param funcs: 関数のlist．上記のpdelistのように並べる．
    :param param: パラメータを並べたndarray(np.array)
    :param vold: 現在時刻における未知関数の値を並べたndarray．pdelistに含まれる関数の引数のx, y, zに相当する．
    :param delta: 計算の1ステップあたりに進める時間
    :return: 計算を1ステップ進めた後の時刻における未知関数の値を並べたndarray．pdelistに含まれる関数の引数のx, y, zに相当する．
    """
    k1 = []
    k2 = []
    k3 = []
    k4 = []
    for fnc in funcs:
        k1.append(_handler(fnc, *param, *vold))

    k11 = vold + delta * 0.5 * np.array(k1)
    for fnc in funcs:
        k2.append(_handler(fnc, *param, *k11))

    k22 = vold + delta * 0.5 * np.array(k2)
    for fnc in funcs:
        k3.append(_handler(fnc, *param, *k22))

    k33 = vold + delta * np.array(k3)
    for fnc in funcs:
        k4.append(_handler(fnc, *param, *k33))
    return vold + (np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) +
                   np.array(k4)) * delta / 6


def main():
    import pickle
    import matplotlib.pyplot as plt
    param = np.array([0.2, 0.2, 7.1])  # レスラー方程式のパラメータ
    init = np.array([4, 4, 1])  # 微分方程式の初期値
    delta = 0.1  # 1ステップあたりの経過時間
    result = []
    vnew = init
    for i in range(1, 10000):
        result.append(vnew)
        vnew = runge_kutta(pdelist, param, vnew, delta)
    result.append(vnew)
    roessler = np.array(result)  # 解の変数．roessler[時刻(sample), 変数の番号]
    with open("signal_roessler.pickle", mode="wb") as f:  # 解の変数を保存
        pickle.dump(roessler, f)
    plt.plot(roessler[:, 0], roessler[:, 1])  # 解の変数の0番目(レスラー方程式のx)を横軸に，解の変数の1番目(レスラー方程式のy)を縦軸に取り，プロットする．
    plt.xlabel("x = roessler[:, 0]")
    plt.ylabel("y = roessler[:, 1]")


if __name__ == "__main__":
    import matplotlib.pyplot as plt

    main()
    plt.show()

上記プログラムを実行すると以下のグラフが表示される．

f:id:kazmus:20190218220859p:plain — レスラー方程式のx-y軌道．パラメータ(a,b,c)=(0.2, 0.2, 7.1)初期値(x,y,z)=(4, 4, 1)．4次のルンゲクッタ法で1ステップの時間を0.1とし，1万ステップ計算した．

2019-02-17

pythonで自己MIC(相互情報量)関数を計算

プログラミング

自己相関関数では，ある信号sig(t)とその信号を時間kだけ遅延させた信号sig(t+k)との相関係数を計算するが，
ここでは相関係数の代わりに相互情報量(MIC)を求める関数automicを示す．
MIC自体はminepyで計算する．

下記プログラムを実行すると，グラフが2つ表示される．

"""
「自己」相互情報量(MIC)関数を計算する．
自己相関関数において，遅延と相関係数の関係を求めるように，
本関数では，遅延とMICの関係を求める．
"""
from minepy import MINE  # MICの計算ライブラリをインポート
from typing import List  # 型表示における"List"を有効化


def automic(sig: List[float]) -> List[float]:
    """
    sig[0]（遅延させていない信号）とsig[k]（k点遅延させた信号）との相互情報量を求める．
    :rtype: List[float]
    :param sig: 1次元の数値配列．時系列信号など．
    :return: 1次元の数値配列mics．mics[k]は遅延kにおけるMIC．
    """
    mine = MINE()
    siglen = len(sig)
    mics = []
    for k in range(0, siglen - 1):
        mine.compute_score(sig[k:-1], sig[0:-k - 1])
        mics.append(mine.mic())
    return mics


if __name__ == "__main__":
    import numpy as np
    from matplotlib import pyplot as plt

    fs = 100  # サンプリング周波数
    f = 4  # 周波数
    tmax = 1  # 終了時刻[s]
    t = np.arange(0, tmax, 1 / fs)
    sig = np.sin(2 * np.pi * f * t)
    plt.plot(t, sig, marker=".")
    plt.xlabel("t [sample]")
    plt.ylabel("signal [a.u.]")
    plt.figure()
    mics = automic(sig)
    plt.plot(mics)
    plt.xlabel("delay k [sample]")
    plt.ylabel("MIC")
    plt.show()

f:id:kazmus:20190217232923p:plain — 上記プログラムでの入力信号

f:id:kazmus:20190217233012p:plain — 計算された自己MIC（相互情報量）関数．入力された信号と遅延させた信号が同位相と逆位相の場合でMICが1に近づいているため，8回ピークがある．ただし，最後のピークは入力波形が短くなりすぎたためノイズが生じている．

2019-02-12

matlabで振幅スペクトルの表示(fft)

プログラミング

matlabで高速フーリエ変換し，振幅スペクトルを表示する．

function showfft(sig,fs)
%sig:信号．1次元配列
%fs:信号のサンプリング周波数
N=length(sig);   % 信号のサンプル数
plot(fs*(0:N-1)/N,abs(fft(sig)/N))   %振幅スペクトルの描画:縦軸が振幅，横軸が周波数
title('amplitude spectrum')
xlabel('frequency [Hz]')
ylabel('amplitude [a.u.]')
end

2019-01-30

haskellで逆引き辞書ソート

プログラミング

単語を羅列したファイルから『逆引き広辞苑』みたいにソートする．

import System.IO

main :: IO ()
main = do
  hin <- openFile "test_words.txt" ReadMode  -- 読みこむファイル：単語リスト
  hout <- openFile "reverse_idx.txt" WriteMode  -- 書き込むファイル：逆引きソートした単語リスト
  contents <- hGetContents hin
  let
    rev_idx = reverseEachWords  -- もう一度各単語を逆にする． [[String]]
              $unduplicate  -- 単語リストをソートして↑ [[String]]
              $reverseEachWords  -- 単語リストの各単語を逆にして↑ [[String]]
              $words  -- 単語ごとにリストに分解して↑ [[String]]
              contents  -- 単語リストの文字列contentsを↑ [String]
  mapM_ (hPutStrLn hout) rev_idx
  hClose hout
  hClose hin
  


-- ソート，重複なくす
unduplicate::Ord a => [a] -> [a]
unduplicate []=[]
unduplicate (xs:xss)=unduplicate smaller++[xs]++unduplicate larger
  where
    smaller = [a | a <- xss, a < xs]
    larger  = [b | b <- xss, b > xs]


-- 単語リストの各要素にreverseを適用する
reverseEachWords :: [[a]] -> [[a]] -- [a]=String：1単語
reverseEachWords xs= map reverse xs