ナンクル力学系

id:salamannの卒論のために沢山FFTしてパワースペクトルのグラフ描くスクリプトをpythonで書いてみたやつ。id:salamannは卒論の謝辞に俺のこと書くべき！あとFFTのパワースペクトルについてちゃんと調べるべき！

細かい説明書こうかとも思ったけど、コメントいっぱい入れたので分かる気がする。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import scipy
import scipy.io
import scipy.fftpack
import pylab
import os

def fft_power(x, cols, sep):
    u"各列のパワースペクトルを計算する関数"
    n = len(x[0])# データの長さ
    pwr = []# 空リスト
    for c in cols:# 各列ごとにfftしてスペクトルに
        y = scipy.fftpack.fft(x[c]) * 2.0 / n
        pwr.append(# リストpwrの最後尾に追加
            scipy.absolute( y[0:n/2+2]**2 )# パワースペクトルの計算
            )
    return scipy.array(pwr)# リストからscipy.arrayに変換

def plot_all(pwr, cols, name, dir):
    u"各列ごとのプロットと、まとめてプロットをする"
    # 各列ごとのプロット
    for p,c in zip(pwr,cols):# zip([1,2,3],[a,b,c])は[(1,a),(2,b),(3,c)]を
        pylab.clf()# pyalbのキャンパスクリアー
        #fig=pylab.figure(figsize=(6,6))# 図の作成
        # x=1,2,3,... y=p[1],p[2],p[3],... でプロット
        pylab.plot(xrange(1,100),p[1:100],label="column: "+str(c))
        pylab.legend()# 凡例を出力
        pylab.savefig( os.path.join(dir,name+"_"+str(c)+".png"))
        #pylab.close(fig)
    pylab.clf()# pyalbのキャンパスクリアー
    # まとめてプロット
    #fig=pylab.figure(figsize=(6,6))# 図の作成
    for p,c in zip(pwr,cols):# zip([1,2,3],[a,b,c])は[(1,a),(2,b),(3,c)]を
        # x=1,2,3,... y=p[1],p[2],p[3],... でプロット
        pylab.plot(xrange(1,100),p[1:100],label="column: "+str(c))
    pylab.legend()# 凡例を出力
    pylab.savefig( os.path.join(dir,name+".png" ))
    #pylab.close(fig)

def main(data_name, cols, num_row, dir_suffix, sep):
    global data
    # data_nameから".csv"を取り除いて、後ろにdir_suffixをつける
    data_dir = data_name.replace(".csv","") + dir_suffix
    if not os.path.isdir(data_dir):
        # data_dirが存在してなければ作成
        os.mkdir(data_dir)
    else:
        print "Already exists: " + data_dir
    # ファイルからデータ読み込み
    print "read %s" % data_name
    data = scipy.io.read_array(data_name, separator=sep)
    stop_row = (len(data)/num_row) * num_row # stop_rowはnum_rowの倍数に（端数の行を処理しないように）
    for i in range(0, stop_row, num_row):
        print "row %d to %d." % (i, i+num_row)
        x = data[i:i+num_row].transpose()# i行からi+num_row行までスライスして転置
        print "FFT..."
        pwr = fft_power(x, cols, sep)
        print "plot..."
        plot_all(pwr, cols, data_name + "_%010d" % i, data_dir)
        # データの保存
        scipy.io.write_array(
            os.path.join(data_dir,"power_%010d.csv" % i), # data_dir以下のpower.csvのパス
            pwr.transpose(), # pwrを転置
            separator="," # カンマ区切り
            )
        # メモリ食いすぎかもしれないので開放(これでできたっけw)
        del pwr
        del x

### 以下は適当にいじってくれ
# 出力するディレクトリのsuffix
data_dir_suffix = "_data"
# FFT解析する列
data_cols=[1,2,3,4,5,6]
# FFT解析する行数の指定
data_row_num = 5000 # 10000だと1secだっけ？
# 図の作成と設定
pylab.figure(figsize=(6,6))

# 実行例 (ipythonに入って、run salamann_fft.pyとしたあと、)
# main("604_fft.csv", data_cols, data_row_num, data_dir_suffix, ",") #カンマ区切りの場合
# main("604_fft.csv", data_cols, data_row_num, data_dir_suffix, None) #スペース区切りの場合

実はこのスクリプトの最初のバージョンには問題があって、plotに時間かかるばかりかメモリを食いまくってpythonが強制終了してしまっていたたらしい。問題は、plot_allの中の今はコメントアウトしてる「fig=pylab.figure(figsize=(6,6))# 図の作成」の部分で毎回figureを作っていたから。これを無くしてpylab.figureを処理の前に実行すればよかったみたい。この関数一回実行すれば同じfigureを使いまわせるんでした＞＜（普段使わないので。。。という言い訳w）

しかし、

• 毎回ちゃんとpylab.closeしてるのになんでメモリ増えるの？
• figure作るってなんでこんなコストかかるの？

とかは良く分からない。この関数をループの中で呼ばなきゃ良いんだろうが。

—

参考：SciPy/FFTPACK : パワースペクトル

pythonの科学計算用ライブラリscipyのarrayの使い方が忘れやすいのでメモ．

scipy.arrayを使って行列から部分行列にアクセスする時には，

• array[ 行開始:行終了, 列開始:列終了 ]
  または
• array[ 行開始:行終了 ][ 列開始:列終了 ]
  ※ だたし，例えばarray[:][0:1]としても列ベクトルは取れない（array[:,0:1]なら大丈夫）ので，前者にそろえるのが良いと思います．

とします．省略可能なのは，

• 開始の0
• 終了が行か列の最後

ipythonでの実行例：

In [1]: import scipy

In [2]: a=scipy.array([[11,12,13],[21,22,23],[31,32,33]])

In [3]: a
Out[3]:
array([[11, 12, 13],
[21, 22, 23],
[31, 32, 33]])

In [4]: a[1:3,0:2]
Out[4]:
array([[21, 22],
[31, 32]])

In [5]: a[1:3,:2]
Out[5]:
array([[21, 22],
[31, 32]])

In [6]: a[1:,:2]
Out[6]:
array([[21, 22],
[31, 32]])