強弱と高低

Python3をメインに

matplotlibでOverflowError: In draw_path: Exceeded cell block limitというエラーが出る

引っかかったところ

matploblibで比較的大きなデータをプロットしようとしたとき, plt.savefigをしようとして OverflowError: In draw_path: Exceeded cell block limit というエラーが出た.

解決策

結構簡単で,

import matplotlib as mpl
mpl.rcParams['agg.path.chunksize'] = 100000

この文を加えるだけ. 100000という数字は大きければとりあえずOKなようだ.

経験的モード分解(Empirical Mode Decomposition)とは

経験的モード分解とは

Huangらによって提案された時系列信号を非定常・非線形な時間-周波数な空間に変換する手法のことで, 大まかな理解として, もともとの信号 x(t) を固有モード関数IMF(Intrinsic Mode Function)と残滓に分解するということを頭に入れておけば良い.

分解するということ

信号を分解する手法は従来から存在しており, FFTやウェーブレット変換がある. これは大学の授業でも扱うくらい古典的で強い手法だ. では何がEMDと違うのか.

FFTやウェーブレット変換はsin関数等の基底を仮定して周波数領域へと変換する一方, EMDは所謂ヒューリスティックな変換を行う. 基底の仮定を置かずに信号を順々に分解していき閾値でもってループの終了とする.

アルゴリズム

これはすでに様々な論文が出回っており再掲するのも忍びない.が, 引用させてもらう.

1. 入力信号 x_0(t)の全ての極値を検出する
2. 極大点と極小点をそれぞれ補間し,上側包絡線 emax(t)と下側包絡線 emin(t)を得る
3. 上側包絡線と下側包絡線の局所平均m(t) = ( emin(t) + emax(t) ) / 2 を算出する
4. 入力と局所平均の差分 y(t) = x(t) - m(t)を入力とみなして 1~4 を繰り返す。
5. 局所平均の標準偏差が閾値以下になった時点で,y(t) を IMF とみなし,ループを終了する
6. 入力信号と IMF の残差 x1(t) = x0(t) - y(t)を新たな入力信号として 1~6 を繰り返す
7. 全ての IMF を抽出したら(x_n(t)の極値が1 つになったら)終了する

出典:http://old.acoust.ias.sci.waseda.ac.jp/publications/happyou/asj/asj-yatabe-2012sep.pdf

要するになにか

暴論ではあるが, EMDは入力信号を,時間軸を中心とした信号に分解することと言って良いかもしれない. 信号の平均が0になるような信号をIMFとしているため波形全体として上下のトレンドがあるような信号は存在しない.

例えば、下のようなトレンドのある信号は, sin波のような切片が存在しない平均面積0の波形に分解できるのだろう. http://touki-tousi.com/wp-content/uploads/2015/09/943578ab5a1148ce847c2eb36745c84d.jpg

何に使うのか

主成分が出てくることから, 支配的な振動を抽出できるという仮説が成り立ち,複雑にノイズが入り込んだ信号波形から,対象の長期信号を取得するような「変調」に使うことができる.そういった考え方の一つとしてノイズ除去がある.IMFの中にホワイトノイズのようなものがあったらそれを削除すれば元々の信号からノイズを削減できたことになる.

昨今の深層学習の成果では, EMDと同じように「信号の分解を何らかの仮定をおかずに行う」ことができてきている.例えば,時間方向での畳込みは,FFTのような操作と同等と解釈することができるということだ. こうした,データに語らせる方法が使えない場合において,EMDのようなヒューリスティックな手法は試してみる価値はあるのだろうと思う.

修正履歴

[2018-06-06] 深層学習への言及を追加した

numpy savetxtで複素数も扱うとき

savetxtを複素数でも扱う

普通savetxtは実数だけで扱っているので複素数のことなど気にも留めない。 そのまま複素数もいけるのかと思いきや、そうは問屋は卸さない。

問題

複素数を含む数をarrayに入れたまま保存すると、(1.00000000e + 1,0000000ej)なんて形になってしまう。この括弧が厄介でこれのせいでpandasでも呼び出せはするがアクセスができないなんてことになる。

解決

numpy savetxt complexとググれば答えは見つかりはする。

python - Writing and reading complex numbers using numpy.savetxt and numpy.loadtxt - Stack Overflow

ここでは

numpy.savetxt(save, array.reshape(1, array.shape[0]), newline = "\r\n", fmt = '%.4f%+.4fj '*4)

こうやるといい、と書いてある。しかしこのまま保存するとcsvにならず連続で4つの実数と虚数の組が並んでしまう。delimiterをホワイトスペースかカンマにするだけなのだが、お好きな方で。 カンマにするときは、

numpy.savetxt(savefilename, array, fmt = '%.4f%+.4fj, '*4, delimiter=',')

とすると綺麗にCSVになってくれる。 *4は、要するに列数だ。100列あるなら*100になるというわけ。また肝はfmt = '%.4f%+.4fj, '*4,でfmtの中にカンマを指定すれば思い通りに出力してくれる。

こんな感じでいじることができるので、任意の出力に対応することができる、らしい。複素数は今回のように出力するのが良いのではないだろうか。

MATLAB melbankmとは

melbankm

voiceboxに入っているフィルターかませる関数。これを元信号に乗算すればフィルターをかけたことになる。

詳細は以下のURL

Description of melbankm

使い方

melbankm(p,n,fs,fl,fh,w)をとり、

  • p : フィルターの数
  • n : FFTの一つの窓に含まれるデータ数
  • fs : サンプリング周波数
  • fl : 最小周波数
  • fh : 最大周波数
  • w : (任意)設定
melw = melbankm(p,n,fs,fl,fh,w)
(signal) * melw

これだけでいける。


利用者がいないのか、それとも英語だけでいけるからなのかわからないが英語くらいしかページがなかったし、あとはもう中国語になっていた。なので支援も兼ねて。

TeXworksのplatex有効化にいろいろ頑張っているようだけど・・・

TeXworksを日本語で使うにあたって

blog.livedoor.jp

のような所謂、platexのタイプセットを追加して・・・とある。もちろんUbuntuでの話。

だけどさ

例えばコマンド設定のところでlatex -interaction=nonstopmode %.tex ってある。これをさplatex -interaction=nonstopmode %.texにすれば動いてしまうわけですよ。 pをコマンドに追加すれば終了なわけで、/usr/local/binに須藤で追加できない環境下だったら先に進めないわけですよ。

こういう方法もある。目的はさっさとTeXworksを動かすことなので手段を間違えると人によっては詰んでしまう。

Pythonで並列処理入門

入門編でやれること

リストを順に処理するときに並列で処理したい!

data_list = ["a", "b", "c", "d"]

for i in data_list:
    #Do something

こういうのをコア数に応じて並列で処理してさくっと終えたい。ときに使える。

どうやるか

from multiprocessing import Pool

data_list = ["a", "b", "c", "d"]
def wrapper_main(nakami):
    #Do something

myprocess = Pool()
myprocess.map(wrapper_main, data_list)
myprocess.close()

実際に書く量はほとんど変わらない。もちろん、ひとつのリストの中身だけで完結することが大前提(並列処理って言ってるしね)。


少しずつ見ていきたい。

from multiprocessing import Pool

これで呼び出し

def wrapper_main(nakami):  
    #Do something

本来やりたい処理をここに入れる。変数nakamiはdata_listの中身のこと。ラッパーを書くときは、リストを引数に与えるのではなくリストの中身を与えた場合の処理を書く。

myprocess = Pool()
myprocess.map(wrapper_main, data_list)
myprocess.close()

これはほとんどおまじない。Poolの中に例えばPool(12)とすると12コアだけ使ってくれる。何も入れないとフルで使うことになる。 肝心なのはプロセスオブジェクト.map(関数名, リスト)という使い方。

これさえ守ればあっという間にマルチプロセスができる。

コアごとに違う処理を押し付けたい場合やwrapperの引数が二次元リストだったり複数あったりする場合は中級編として既存のブログ等を参考にすればできる。マルチスレッドとマルチプロセスの違いなども割愛。

pandasで条件にあった行を抽出してさらに列を指定する

やりたいこと

ある巨大行列のなかで、特定の列の値がXである行を抽出する。そしてさらに利用する列を指定する。

コード

import pandas as pd
df = pd.read_csv("test.csv", header=None)

df[df[2] == X].ix[:,5:90]

df[2] == Xが第2列の値がXである行を抽出してくれる。
更に加えて、.ix[;5:8]をすることで抽出した行の中から5列から7列目までの列を取り出してくれる。

とても便利ですねpandas