Created
January 30, 2013 15:42
-
-
Save non117/4674114 to your computer and use it in GitHub Desktop.
研究用ライブラリ
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# -*- 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 | |
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(result[-1]) | |
result.append(result[-1]) | |
result.append(result[-1]) | |
return result | |
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 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=[]): | |
pylab.clf() | |
for elem in lis: | |
pylab.plot(t, elem) | |
pylab.xlabel(label_t) | |
pylab.ylabel(label_y) | |
pylab.title(title) | |
pylab.grid = True | |
if legend: pylab.legend(legend, loc=0, frameon=False) | |
if axis:pylab.axis(axis) | |
pylab.savefig("%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(True) | |
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): | |
with open(filename, "wb") as f: | |
pickle.dump(obj, f) | |
def deserialize(filename): | |
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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment