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