卒論で使った研究用ライブラリ

researchlib.py

# -*- coding: utf-8 -*-
import cPickle as pickle

import pylab
#from mpl_toolkits.mplot3d.axes3d import Axes3D
from numpy import arange, sqrt, cos, sin, matrix, identity, pi, ones, convolve
from scipy import signal
from scipy.interpolate import UnivariateSpline

def squarediff(x, y):
    """ 二乗差分 """
    diff = x - y
    return diff * diff

def time_scale(x, freq):
    ''' 時間の長さを返す '''
    length = len(x)
    ts = arange(0, length/freq, 1/freq)
    return ts

def newton_diff(lis, h):
    ''' ニュートン法を使った微分 '''
    def fprime(h, f0, f1, f2, f3):
        ''' 微分の計算の一部 '''
        return (-11*f0 + 18*f1 - 9*f2 + 2*f3)/(6*h)
    
    result = []
    for i, x in enumerate(lis):
        if i+3 <= len(lis)-1:
            result.append(fprime(h, x, lis[i+1], lis[i+2], lis[i+3]))
    result.append((lis[-3] - lis[-2])/h)
    result.append((lis[-2] - lis[-1])/h)
    result.append(result[-1])
    return result

def simple_moving_average(x, n):
    ''' 単純移動平均, xが対象データ, nが窓の大きさ '''
    # n = 0.443*fs/fc
    w = ones(n)/float(n)
    return convolve(x, w, "valid")

def lpf(x, cutoff_freq, samp_rate):
    ''' バターワースフィルタを設計してLPFしてくれる '''
    norm_pass = cutoff_freq/(samp_rate/2)
    norm_stop = 1.5*norm_pass
    N, Wn = signal.buttord(wp=norm_pass, ws=norm_stop, gpass=2, gstop=30, analog=0)
    b, a = signal.butter(N, Wn, btype='low', analog=0, output='ba')
    y = signal.lfilter(b, a, x)
    return y


def spline(xs, ys):
    ''' スプライン補間 '''
    s = UnivariateSpline(xs, ys, s=1)
    return s(xs)

def distance(xd, yd, zd):
    ''' 
        ユークリッド距離 
        xd := x1 - x2
    '''
    return sqrt(xd * xd + yd * yd + zd * zd)

def identity4():
    ''' 4x4の単位行列 '''
    return matrix(identity(4))

def rotate_local(m, axis, deg):
    ''' axisを中心に, deg度回転 '''
    rad = deg*pi/180.
    rot = identity4()
    if axis == "x":
        a,b,c,d = 5,6,9,10
    elif axis == "y":
        a,b,c,d = 10,8,2,0
    elif axis == "z":
        a,b,c,d = 0,1,4,5
    rot.put(a, cos(rad))
    rot.put(b, -sin(rad))
    rot.put(c, sin(rad))
    rot.put(d, cos(rad))
    return m * rot

def rotate_world(m, axis, deg):
    ''' axisを中心に, deg度回転 '''
    rad = deg*pi/180.
    rot = identity4()
    if axis == "x":
        a,b,c,d = 5,6,9,10
    elif axis == "y":
        a,b,c,d = 10,8,2,0
    elif axis == "z":
        a,b,c,d = 0,1,4,5
    rot.put(a, cos(rad))
    rot.put(b, -sin(rad))
    rot.put(c, sin(rad))
    rot.put(d, cos(rad))
    return rot * m

def translate_world(m, x, y, z):
    ''' 並進 '''
    t = identity4()
    t.put(3, x)
    t.put(7, y)
    t.put(11, z)
    return t * m

def gen_graph(t, lis, label_t, label_y, title, legend=[], axis=[], filename=""):
    pylab.clf()
    for elem in lis:
        pylab.plot(t, elem)
    pylab.xlabel(label_t)
    pylab.ylabel(label_y)
    pylab.title(title)
    pylab.grid()
    if legend: pylab.legend(legend, loc=0, frameon=False)
    if axis:pylab.axis(axis)
    if filename:
        pylab.savefig(u"%s.png" % filename)
    else:
        pylab.savefig(u"%s.png" % title)

#def gen_3d_graph(t, x, y, label_t, label_x, label_y, title, legend=[]):
#    pylab.clf()
#    fig = pylab.figure()
#    ax = Axes3D(fig)
#    ax.set_title(title)
#    ax.set_xlabel(label_x)
#    ax.set_ylabel(label_y)
#    ax.set_zlabel(label_t)
#    ax.grid()
#    for i, v in enumerate(x):
#        ax.plot3D(x[i], y[i], t)
#    if legend: ax.legend(legend, loc=0, frameon=False)
#    pylab.savefig("%s.png" % title)

def serialize(obj, filename):
    '''pythonのオブジェクトをそのままdump'''
    with open(filename, "wb") as f:
        pickle.dump(obj, f)

def deserialize(filename):
    '''dumpしたのをメモリにロード'''
    with open(filename, "rb") as f:
        return pickle.load(f)

class Node:
    def __init__(self, x, y, cost, prev_node):
        self.coord = (x, y)
        self.cost = cost #その点までのコスト合計
        self.prev = prev_node #前の点

def dynamic_timewarp(a, b):
    m = len(a)
    n = len(b)
    cost_matrix = []
    # コストを求めてリストに詰め込む
    for x in range(m):
        array_ = []
        for y in range(n):
            diff = squarediff(a[x], b[y])
            array_.append(diff)
        cost_matrix.append(array_)
    
    # 初期化. 参照のコピーになるため, [[None]*n]*mとかしちゃ駄目
    path_matrix = []
    for x in range(m):
        temp = [None]*n
        path_matrix.append(temp)
    # 出発点の処理
    path_matrix[0][0] = Node(0, 0, cost_matrix[0][0], None)
    
    # 端っこの処理
    for x in range(1,m):
        prev_node = path_matrix[x-1][0]
        path_matrix[x][0] = Node(x, 0, prev_node.cost + cost_matrix[x][0], (x-1, 0))
    
    for y in range(1,n):
        prev_node = path_matrix[0][y-1]
        path_matrix[0][y] = Node(0, y, prev_node.cost + cost_matrix[0][y], (0, y-1))

    # 残りの経路を埋める処理
    for x in range(1,m):
        for y in range(1,n):
            min_node = min(path_matrix[x-1][y], path_matrix[x][y-1], path_matrix[x-1][y-1], key=lambda n:n.cost)
            path_matrix[x][y] = Node(x, y, min_node.cost + cost_matrix[x][y], min_node.coord)

    # 最小経路を取り出す
    node = path_matrix[m-1][n-1]
    min_path = []
    while(node.prev != None):
        min_path.append(node.coord)
        x, y = node.prev
        node = path_matrix[x][y]
    min_path.append((0, 0))

    # 結果
    result = [None]*m
    delay = [None]*m
    for i, j in min_path:
        result[i] = b[j]
        delay[i] = i - j
    
    return result, delay

try:
    from dtw import dtw as dynamic_timewarp
except ImportError:
    pass