# azure notebooksデフォルトのSciPyをpipで更新
!pip install -U scipy
Requirement already up-to-date: scipy in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (1.1.0)
# 本発表では以下を仮定
import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt
# version
print("numpy version", np.version.full_version)
print("scipy version", sp.version.full_version)
print("matplotlib version", mpl.__version__)
numpy version 1.14.2 scipy version 1.1.0 matplotlib version 2.2.2
本発表では基本的にnp.dotの代わりに@演算子を使用
# 本発表では以下を仮定
import scipy.io as spio
import scipy.special as special
import scipy.linalg as la
import scipy.fftpack as fft
import scipy.optimize as opt
import scipy.signal as sg
import scipy.interpolate as interp
import scipy.stats as stats
a = np.ones((3, 3))
spio.savemat('file.mat', {'a': a}) # savematには辞書を渡す
data = spio.loadmat('file.mat')
data['a']
array([[1., 1., 1.], [1., 1., 1.], [1., 1., 1.]])
a = np.ones(3)
spio.savemat('file.mat', {'a': a})
print(a)
print(spio.loadmat('file.mat')['a']) # matlabには1次元配列がない
print(spio.loadmat('file.mat', squeeze_me=True)["a"])
[1. 1. 1.] [[1. 1. 1.]] [1. 1. 1.]
scipy.special.jn()
scipy.special.ellpj()
scipy.special.gamma()
、対数ガンマ関数 scipy.special.gammaln()
scipy.special.erf()
special?
cimport scipy.special.cython_special
でcythonコード上でcimport%load_ext cython
The cython extension is already loaded. To reload it, use: %reload_ext cython
%%cython -a
cimport scipy.special.cython_special as csc
cdef:
double x = 1
double complex z = 1 + 1j
double si, ci, rgam
double complex cgam
rgam = csc.gamma(x)
print(rgam)
cgam = csc.gamma(z)
print(cgam)
csc.sici(x, &si, &ci)
print(si, ci)
Generated by Cython 0.28.3
Yellow lines hint at Python interaction.
Click on a line that starts with a "+
" to see the C code that Cython generated for it.
01:
02: cimport scipy.special.cython_special as csc
03:
04: cdef:
+05: double x = 1
__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_x = 1.0;
+06: double complex z = 1 + 1j
__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_z = __Pyx_c_sum_double(__pyx_t_double_complex_from_parts(1, 0), __pyx_t_double_complex_from_parts(0, 1.0));
07: double si, ci, rgam
08: double complex cgam
09:
+10: rgam = csc.gamma(x)
__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_rgam = __pyx_fuse_1__pyx_f_5scipy_7special_14cython_special_gamma(__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_x, 0);
+11: print(rgam)
__pyx_t_1 = PyFloat_FromDouble(__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_rgam); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 11, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_2 = __Pyx_PyObject_CallOneArg(__pyx_builtin_print, __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 11, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_2); __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+12: cgam = csc.gamma(z)
__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_cgam = __pyx_fuse_0__pyx_f_5scipy_7special_14cython_special_gamma(__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_z, 0);
+13: print(cgam)
__pyx_t_2 = __pyx_PyComplex_FromComplex(__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_cgam); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 13, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_2); __pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_builtin_print, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 13, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+14: csc.sici(x, &si, &ci)
__pyx_fuse_1__pyx_f_5scipy_7special_14cython_special_sici(__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_x, (&__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_si), (&__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_ci));
+15: print(si, ci)
__pyx_t_1 = PyFloat_FromDouble(__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_si); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 15, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_2 = PyFloat_FromDouble(__pyx_v_46_cython_magic_684790e4af39b5c80e39e5974eacdfbe_ci); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 15, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_2); __pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 15, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_3); __Pyx_GIVEREF(__pyx_t_1); PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_1); __Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_t_2); __pyx_t_1 = 0; __pyx_t_2 = 0; __pyx_t_2 = __Pyx_PyObject_Call(__pyx_builtin_print, __pyx_t_3, NULL); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 15, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_2); __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
scipy.linalg
はLAPACKを用いたより完璧なラッパーscipy.linalg
はnumpy.linalg
に存在しない機能を持つscipy.linalg.eig()
は一般化固有値問題を解くために第二の行列を引数に取ることが可能)linalg.det
で計算可能# 行列式の計算
A = np.array([[1, 2],
[3, 4]])
print(la.det(A)) # 1*4-2*3=-2
A = np.array([[3, 2],
[6, 4]])
print(la.det(A) ) # 3*4-2*6=0
print(la.det(np.ones((3, 4)))) # 正方行列ではない
-2.0 6.661338147750939e-16
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-13-c0c6957027f4> in <module>() 8 print(la.det(A) ) # 3*4-2*6=0 9 ---> 10 print(la.det(np.ones((3, 4)))) # 正方行列ではない ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/scipy/linalg/basic.py in det(a, overwrite_a, check_finite) 1030 a1 = _asarray_validated(a, check_finite=check_finite) 1031 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]: -> 1032 raise ValueError('expected square matrix') 1033 overwrite_a = overwrite_a or _datacopied(a1, a) 1034 fdet, = get_flinalg_funcs(('det',), (a1,)) ValueError: expected square matrix
linalg.inv
を用いる# 逆行列の計算
A = np.array([[1, 2],
[3, 4]])
A_inv = la.inv(A)
print(A_inv)
print(A @ A_inv) # AA^-1 ~= I
print(np.allclose(A @ A_inv, np.eye(2)))
[[-2. 1. ] [ 1.5 -0.5]] [[1.0000000e+00 0.0000000e+00] [8.8817842e-16 1.0000000e+00]] True
# 特異行列(行列式が0)の場合は逆行列が求められないのでLinAlgError例外を投げる
A = np.array([[3, 2],
[6, 4]])
la.inv(A)
--------------------------------------------------------------------------- LinAlgError Traceback (most recent call last) <ipython-input-15-6b5557458fe3> in <module>() 2 A = np.array([[3, 2], 3 [6, 4]]) ----> 4 la.inv(A) ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/scipy/linalg/basic.py in inv(a, overwrite_a, check_finite) 973 inv_a, info = getri(lu, piv, lwork=lwork, overwrite_lu=1) 974 if info > 0: --> 975 raise LinAlgError("singular matrix") 976 if info < 0: 977 raise ValueError('illegal value in %d-th argument of internal ' LinAlgError: singular matrix
linalg.solve
を用いるlinalg.solve
を用いた方が良いimport numpy as np
from scipy import linalg
A = np.array([[1, 2], [3, 4]])
b = np.array([[5], [6]])
print("逆行列による解法")
print(linalg.inv(A) @ b) # 遅い
print((A @ linalg.inv(A) @ b) - b) # 確認
print("solveによる解法")
print(np.linalg.solve(A, b)) # 高速
print((A @ np.linalg.solve(A, b)) - b) # 確認
逆行列による解法 [[-4. ] [ 4.5]] [[0.00000000e+00] [1.77635684e-15]] solveによる解法 [[-4. ] [ 4.5]] [[0.] [0.]]
linalg.norm(a, ord)
で行列やベクトルの様々なノルムを計算可能A = np.array([[1,2],[3,4]])
print(A)
print(la.norm(A))
print(la.norm(A, 'fro')) # デフォルト
print(np.sqrt(np.trace(A.T @ A)))
print(la.norm(A,1)) # L1ノルム(各列の絶対値和のうち最大のもの)
print(la.norm(A,-1)) # L1ノルム(各列の絶対値和のうち最小のもの)
print(la.norm(A,np.inf)) # L∞ノルム(各行の絶対値和のうち最大のもの)
[[1 2] [3 4]] 5.477225575051661 5.477225575051661 5.477225575051661 6.0 4.0 7.0
linalg.lstsq
linalg.pinv
linalg.pinv2
c1, c2 = 5.0, 2.0
i = np.r_[1:11]
xi = 0.1*i # 説明変数
yi = c1*np.exp(-xi) + c2*xi
zi = yi + 0.05 * np.max(yi) * np.random.randn(len(yi)) # 観測値
A = np.c_[np.exp(-xi)[:, np.newaxis],
xi[:, np.newaxis]] # 計画行列
print(A)
[[0.90483742 0.1 ] [0.81873075 0.2 ] [0.74081822 0.3 ] [0.67032005 0.4 ] [0.60653066 0.5 ] [0.54881164 0.6 ] [0.4965853 0.7 ] [0.44932896 0.8 ] [0.40656966 0.9 ] [0.36787944 1. ]]
# 計画行列と観測値を与えてモデルパラメータを推定
c, resid, rank, sigma = linalg.lstsq(A, zi)
# パラメータ、残差の総和、Aのrank, Aの特異値
print(c, resid, rank, sigma)
xi2 = np.r_[0.1:1.0:100j]
yi2 = c[0]*np.exp(-xi2) + c[1]*xi2
yi_true = c1*np.exp(-xi2) + c2*xi2
plt.plot(xi,zi,'x',xi2,yi2, xi2, yi_true)
plt.axis([0,1.1,3.0,5.5])
plt.xlabel('$x_i$')
plt.title('Data fitting with linalg.lstsq')
plt.show()
[5.22765058 1.77693552] 0.17102621042273095 2 [2.58763467 1.02933937]
linalg.eig
linalg.eigvals
eigh
, eigvalsh
が存在scipy.sparse.linalg
のeigs
, eighs
A = np.array([[1, 2], [3, 4]])
(l1, l2), v = la.eig(A)
print(l1, l2) # 固有値
print(v) # 各列ベクトルが固有ベクトル
print(la.norm(v, axis=0)) # 固有ベクトルは正規化されている
# Ax = λxとなることを確認
# |1 2| = |v1 v2 ||l1 0| = |l1*v1 l2*v2|
# |3 4| | 0 l2|
print(A@v - v@np.diag((l1,l2)))
(-0.3722813232690143+0j) (5.372281323269014+0j) [[-0.82456484 -0.41597356] [ 0.56576746 -0.90937671]] [1. 1.] [[ 0.00000000e+00+0.j -4.44089210e-16+0.j] [ 5.55111512e-17+0.j 0.00000000e+00+0.j]]
# 正方行列に対する特異値分解
A = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1])
U, s, Vh = la.svd(A)
Sig = np.diag(s)
print(A)
print(U.shape, s.shape, Sig.shape, Vh.shape)
print(U @ Sig @ Vh)
# 非正方行列に対する特異値分解
A = np.array([[1,2,3],[4,5,6]])
M,N = A.shape
U,s,Vh = linalg.svd(A)
Sig = la.diagsvd(s, M, N)
print("")
print(A)
print(U.shape, s.shape, Sig.shape, Vh.shape)
print(U@Sig@Vh)
[[1 1 2] [3 4 5] [6 7 9]] (3, 3) (3,) (3, 3) (3, 3) [[1. 1. 2.] [3. 4. 5.] [6. 7. 9.]] [[1 2 3] [4 5 6]] (2, 2) (2,) (2, 3) (3, 3) [[1. 2. 3.] [4. 5. 6.]]
scipy.fftpack.fft
でFFTを実行scipy.fftpack.fftfreq
で周波数軸の値を取得import IPython.display
fs = 48000.0; nfft = 2048
xs = np.arange(1.0 * fs) # 1[s]
ts = xs / fs
hz = 2000.0
ys = np.sin(2.0 * np.pi * hz * ts) + 0.1 * np.random.randn(int(fs))
fig = plt.figure(1, figsize=(12,4))
ax = fig.add_subplot(1, 1, 1)
ax.plot(ts[:nfft], ys[:nfft])
ax.set_xlabel("time [s]"); ax.set_ylabel("amplitude")
IPython.display.Audio(data=ys, rate=fs)
# 窓関数を取得
win = sg.get_window("hanning", nfft)
# FFTを実行
Ys = fft.fft(win * ys[:nfft], nfft)[:nfft//2+1]
# 周波数軸を取得
f = fft.fftfreq(nfft, d=1.0/fs)[:nfft//2+1]
f[-1] = fs//2 # デフォルトだと-24kHz
fig = plt.figure(1, figsize=(12,4));
ax = fig.add_subplot(1, 2, 1)
ax.plot(ts[:nfft], win * ys[:nfft])
ax.set_xlabel("time [s]"); ax.set_ylabel("amplitude")
ax = fig.add_subplot(1, 2, 2)
ax.semilogx(f, np.abs(Ys))
ax.set_xlabel("frequency [Hz]"); ax.set_ylabel("amplitude");
scipy.interpolate.interp1d
で補間関数を生成# 観測値
xs = np.linspace(0, 1, 10)
ys = np.sin(2.0 * np.pi * xs) + 0.01 * np.random.randn(len(xs))
# 異なる2種類の補間関数の生成
linear_interp = interp.interp1d(xs, ys)
cubic_interp = interp.interp1d(xs, ys, kind='cubic')
# 補間値を算出
xs_new = np.linspace(0, 1, 50)
ys_lin = linear_interp(xs_new)
ys_cubic = cubic_interp(xs_new)
fig = plt.figure(figsize=(12,6)); ax = fig.add_subplot(111)
ax.plot(xs, ys, "x-", label="obserbed")
ax.plot(xs_new, ys_lin, label="linear")
ax.plot(xs_new, ys_cubic, label="cubic")
ax.grid(); ax.legend()
<matplotlib.legend.Legend at 0x1d1fb80080>
numpy.random.randn
でホワイトノイズを生成scipy.signal.firwin
でローパスフィルタを設計scipy.signal.lfilter
でローパスフィルタをホワイトノイズに対して適用# ホワイトノイズの生成
time_s = 5.0
ys = 0.1 * np.random.randn(int(time_s * fs))
IPython.display.Audio(ys, rate=fs)
nlpf = 96
cutoff_hz = 4000.0
lpf = sg.firwin(nlpf, cutoff=cutoff_hz, fs=fs)
ys_lpf = sg.lfilter(lpf, [1.0], ys)
IPython.display.Audio(ys_lpf, rate=fs)
fig = plt.figure(1, figsize=(12,4));
ax = fig.add_subplot(1, 2, 1)
ax.plot(f, 20.0*np.log10(np.abs(fft.fft(ys, nfft)[:nfft//2+1])))
ax.set_xlabel("time [s]"); ax.set_ylabel("power [dB]")
ax = fig.add_subplot(1, 2, 2)
ax.plot(f, 20.0*np.log10(np.abs(fft.fft(ys_lpf, nfft))[:nfft//2+1]))
ax.set_xlabel("frequency [Hz]"); ax.set_ylabel("power [dB]")
Text(0,0.5,'power [dB]')
scipy.optimize.curve_fit
で曲線フィッティングscipy.optimize.minimize
で最小値を見つけるscipy.optimize.root
で方程式の根を見つける# ノイズの乗った正弦波を観測値とする
x_data = np.linspace(-5, 5, num=50)
y_data = 2.9 * np.sin(1.5 * x_data) + np.random.normal(size=50)
def test_func(x, a, b):
return a * np.sin(b * x)
# curve_fitを用いてa, bを推定
params, params_covariance = opt.curve_fit(test_func, x_data, y_data, p0=[2, 2])
print(params, params_covariance)
[3.05684189 1.48600684] [[ 0.055698 -0.0005898 ] [-0.0005898 0.00061353]]
fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(111)
ax.plot(x_data, y_data, "x")
ax.plot(x_data, test_func(x_data, params[0], params[1]))
[<matplotlib.lines.Line2D at 0x1d1fced470>]
# 最小値を見つけたい対象
def f(x):
return x**2 + 10*np.sin(x)
# 大域解が-1.3, 局所解が3.8
x = np.arange(-10.0, 10.0, 0.1)
fig = plt.figure(); ax = fig.add_subplot(111)
ax.plot(x, f(x))
[<matplotlib.lines.Line2D at 0x1d1fd7af60>]
result = opt.minimize(f, x0=0.0) # x0で探索開始点を与える
print(result)
print(result.x)
fun: -7.945823375615215 hess_inv: array([[0.08589237]]) jac: array([-1.1920929e-06]) message: 'Optimization terminated successfully.' nfev: 18 nit: 5 njev: 6 status: 0 success: True x: array([-1.30644012]) [-1.30644012]
# nfev(評価回数)が12に減る
opt.minimize(f, x0=0, method="L-BFGS-B")
fun: array([-7.94582338]) hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64> jac: array([-1.50990331e-06]) message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' nfev: 12 nit: 5 status: 0 success: True x: array([-1.30644014])
# 3.0から始めると局所解に陥る
result = opt.minimize(f, x0=3.0, method="L-BFGS-B")
print(result)
print(result.x)
fun: array([8.31558558]) hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64> jac: array([-1.77635684e-07]) message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' nfev: 14 nit: 6 status: 0 success: True x: array([3.83746709]) [3.83746709]
# `scipy.optimize.basinhopping`を用いて大域解を探す
opt.basinhopping(f, 1.0)
fun: -7.945823375615284 lowest_optimization_result: fun: -7.945823375615284 hess_inv: array([[0.08581295]]) jac: array([1.1920929e-07]) message: 'Optimization terminated successfully.' nfev: 15 nit: 3 njev: 5 status: 0 success: True x: array([-1.30644001]) message: ['requested number of basinhopping iterations completed successfully'] minimization_failures: 0 nfev: 1515 nit: 100 njev: 505 x: array([-1.30644001])
# 探索範囲に制約を与えることも可能
res = opt.minimize(f, x0=1,
bounds=((0, 10), ))
res.x
array([0.])
# f=0となる点をoptimize.rootで探す
root = opt.root(f, x0=1)
print(root)
print(root.x)
fjac: array([[-1.]]) fun: array([0.]) message: 'The solution converged.' nfev: 10 qtf: array([1.33310463e-32]) r: array([-10.]) status: 1 success: True x: array([0.]) [0.]
# もうひとつの根を探す
root = opt.root(f, x0=-3.0)
print(root)
print(root.x)
fjac: array([[-1.]]) fun: array([-1.77635684e-15]) message: 'The solution converged.' nfev: 8 qtf: array([-2.83790769e-11]) r: array([12.84592704]) status: 1 success: True x: array([-2.47948183]) [-2.47948183]
# サンプリングした点のヒストグラム(頻度分布)と確率分布を同時に描画
samples = np.random.normal(size=1000) # sampling
bins = np.arange(-4, 5)
histogram = np.histogram(samples, bins=bins, normed=True)[0] # plot histogram
bins = 0.5*(bins[1:] + bins[:-1])
pdf = stats.norm.pdf(bins) # normは正規分布の意味
fig = plt.figure(); ax = fig.add_subplot(111)
ax.plot(bins, histogram, label="histgram")
ax.plot(bins, pdf, label="probability density")
ax.grid(); ax.legend()
<matplotlib.legend.Legend at 0x1d1fc61048>
# 観測値を確率分布にフィッティング
loc, std = stats.norm.fit(samples)
print(loc, std)
xs = np.linspace(-3.0, 3.0, 1000)
f = lambda x: (1.0/np.sqrt(2.0*np.pi*std))*np.exp(-(x - loc)**2/ (2.0 * std**2))
plt.plot(xs, f(xs))
0.014729738495055104 0.9911940562256367
[<matplotlib.lines.Line2D at 0x1d1fad3048>]
# 平均、中央値、四分位数
print(np.mean(samples))
print(np.median(samples))
print(stats.scoreatpercentile(samples, 50))
print(stats.scoreatpercentile(samples, 90))
0.014729738495055104 0.03258817640410924 0.03258817640410924 1.300282756931229
scipy.spatial.distance
のcdist
, pdist
のみ紹介from scipy.spatial import distance
coords = [(35.0456, -85.2672),
(35.1174, -89.9711),
(35.9728, -83.9422),
(36.1667, -86.7833)]
print(np.sqrt((coords[0][0] - coords[1][0])**2
+(coords[0][1] - coords[1][1])**2))
distance.cdist(coords, coords, 'euclidean')
4.704447943170383
array([[0. , 4.70444794, 1.6171966 , 1.88558331], [4.70444794, 0. , 6.0892811 , 3.35605413], [1.6171966 , 6.0892811 , 0. , 2.84770898], [1.88558331, 3.35605413, 2.84770898, 0. ]])
# pdistは対角を除いた上三角成分のみ返す
presult = distance.pdist(coords)
print(presult)
# squareformで距離行列に変換
print(distance.squareform(presult))
[4.70444794 1.6171966 1.88558331 6.0892811 3.35605413 2.84770898] [[0. 4.70444794 1.6171966 1.88558331] [4.70444794 0. 6.0892811 3.35605413] [1.6171966 6.0892811 0. 2.84770898] [1.88558331 3.35605413 2.84770898 0. ]]
import scipy.misc
import matplotlib.cm as cm
face = scipy.misc.face()
print(face.shape)
print(face.max())
print(face.dtype)
fig = plt.figure(1)
ax = fig.add_subplot(1, 1, 1)
ax.imshow(face)
plt.show()
(768, 1024, 3) 255 uint8
face = scipy.misc.face(gray=True)
fig = plt.figure(1)
ax = fig.add_subplot(1, 1, 1)
ax.imshow(face, cmap=cm.gray)
plt.show()
# グレースケールのデモ用画像
ascent = scipy.misc.ascent()
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.imshow(ascent, cmap=cm.gray)
plt.show()
from scipy.misc import derivative
# x0の点におけるn次の中心差分を返す
# n==1の場合はf(x+h)-f(x-h)/2h
# 内部的にcentral_diff_weightを使う
# derivative(func, x0, dx=1.0, n=1, args=(), order=3)
f = lambda x: x**3 + x**2
derivative(f, 1.0, dx=1e-6)
4.999999999921734
print((f(1.0 + 1e-6) - f(1.0 - 1e-6)) / (2.0 * 1.e-6))
4.999999999921734
librosa.decompose.decompose
と同等の処理!pip install librosa
Requirement already satisfied: librosa in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.6.0) Requirement already satisfied: numpy>=1.8.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.14.2) Requirement already satisfied: resampy>=0.2.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.2.0) Requirement already satisfied: scikit-learn!=0.19.0,>=0.14.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.19.1) Requirement already satisfied: audioread>=2.0.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (2.1.5) Requirement already satisfied: six>=1.3 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.11.0) Requirement already satisfied: scipy>=0.14.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (1.1.0) Requirement already satisfied: joblib>=0.7.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (0.11) Requirement already satisfied: decorator>=3.0.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from librosa) (4.3.0) Requirement already satisfied: numba>=0.32 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from resampy>=0.2.0->librosa) (0.38.0) Requirement already satisfied: llvmlite>=0.23.0dev0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from numba>=0.32->resampy>=0.2.0->librosa) (0.23.0)
import librosa
import librosa.display
y, sr = librosa.load(librosa.util.example_audio_file())
IPython.display.Audio(data=y, rate=sr)
np.pad(np.arange(5), 2, "reflect")
array([2, 1, 0, 1, 2, 3, 4, 3, 2])
# strideトリックを使った信号を短時間波形フレームの配列に変換する関数
from numpy.lib.stride_tricks import as_strided
def split_frames(xs, frame_width=1024, noverlap=2):
nsample = np.alen(xs)
nshift = (frame_width // noverlap)
nframe = int(nsample / nshift) - noverlap + 1
nbyte = xs.dtype.itemsize
# as_strided(xs, ys_shape, new_stride)
xs_split = as_strided(xs, (nframe, frame_width), (nbyte * nshift, nbyte))
return xs_split
# 信号の短時間波形への切り出し
import matplotlib.gridspec as gridspec
fig = plt.figure(1, figsize=(12,8))
gs = gridspec.GridSpec(2, 21)
ax = fig.add_subplot(gs[0, :-1])
xs = np.arange(len(y)) / sr
ax.plot(xs[:sr], y[:sr])
ymax = np.max(y[:sr])
ax.set_xlim([0.0, 1.0]); ax.set_ylim([-ymax, ymax]);
ax.set_xlabel("time [s]"); ax.set_ylabel("amplitude"); ax.grid()
frame_width = sr * 0.1
xs = np.arange(frame_width) / sr
for i, y_cut in enumerate(split_frames(y, int(frame_width), 1)[:10]):
ax = fig.add_subplot(gs[1, 2*i:2*(i+1)])
ax.plot(xs, y_cut)
ax.set_ylim([-ymax, ymax]); ax.grid()
# 切り出し間隔を1/2オーバーラップした場合
import matplotlib.gridspec as gridspec
fig = plt.figure(1, figsize=(12,8))
gs = gridspec.GridSpec(3, 21)
ax = fig.add_subplot(gs[0, :-1])
ax.plot(y[:sr])
ymax = np.max(y[:sr])
ax.set_ylim([-ymax, ymax]); ax.grid()
frame_width = sr * 0.1
half_width = int(frame_width // 2)
xs = np.arange(frame_width)
for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
if i % 2 == 0:
ax = fig.add_subplot(gs[1, i:i+2])
ax.plot(xs[:half_width], y_cut[:half_width])
ax.plot(xs[half_width:], y_cut[half_width:])
else:
ax = fig.add_subplot(gs[2, i:i+2])
ax.plot(xs[half_width:], y_cut[half_width:])
ax.plot(xs[:half_width], y_cut[:half_width])
ax.set_ylim([-ymax, ymax]); ax.grid()
frame_width = int(sr * 0.1)
win = sg.get_window("hanning", frame_width)
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.plot(win)
ax.grid()
# 更に各短時間波形に窓掛けを行った場合
import matplotlib.gridspec as gridspec
fig = plt.figure(1, figsize=(12,8))
gs = gridspec.GridSpec(3, 21)
ax = fig.add_subplot(gs[0, :-1])
ax.plot(y[:sr])
ymax = np.max(y[:sr])
ax.set_ylim([-ymax, ymax]); ax.grid()
frame_width = sr * 0.1
half_width = int(frame_width // 2)
win = sg.get_window("hanning", int(frame_width))
xs = np.arange(frame_width)
for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
y_win = win * y_cut
if i % 2 == 0:
ax = fig.add_subplot(gs[1, i:i+2])
ax.plot(xs[:half_width], y_win[:half_width])
ax.plot(xs[half_width:], y_win[half_width:])
else:
ax = fig.add_subplot(gs[2, i:i+2])
ax.plot(xs[half_width:], y_win[half_width:])
ax.plot(xs[:half_width], y_win[:half_width])
ax.set_ylim([-ymax, ymax]); ax.grid()
# 1/4オーバーラップして窓掛けした各波形をそれぞれFFT
nfft = 2048
win = sg.get_window("hanning", nfft)
y_center = np.pad(y, int(nfft//2), "reflect")
fft_list = []
phase_list = []
for y_cut in split_frames(y_center, nfft, 4):
y_end = min(len(y_cut), nfft)
y_win = np.zeros(nfft)
y_win[:y_end] = win[:y_end] * y_cut[:y_end]
Y = fft.fft(y_win, nfft)
fft_list.append(Y[:nfft//2+1])
phase_list.append(np.exp(1.j * np.angle(Y[:nfft//2+1])))
my_D = np.array(fft_list, dtype=np.complex64).T
my_phase = np.array(phase_list, dtype=np.complex64).T
fig = plt.figure(1, figsize=(12,4))
ax = fig.add_subplot(2,1,1)
librosa.display.specshow(librosa.amplitude_to_db(np.abs(my_D), ref=np.max), y_axis="log")
<matplotlib.axes._subplots.AxesSubplot at 0x1d216b7320>
# 得られたスペクトログラムnp.abs(my_D.T)に対してNMFを実行
from sklearn.decomposition import NMF
model = NMF(n_components=8)
W = model.fit_transform(np.abs(my_D.T))
H = model.components_
components = H.T
activations = W.T
fig = plt.figure(2, figsize=(12,4))
ax = fig.add_subplot(1,2,1)
librosa.display.specshow(activations, x_axis='time')
ax.set_xlabel('Time'); ax.set_ylabel('Component'); ax.set_title('Activations')
ax = fig.add_subplot(1,2,2)
librosa.display.specshow(librosa.amplitude_to_db(components, ref=np.max), y_axis='log')
ax.set_xlabel('Component'); ax.set_ylabel('Frequency'); ax.set_title('Components')
Text(0.5,1,'Components')
# 1.e-12以下となっている要素の数をコンポーネントごとに算出
for i in range(8):
print(sum(np.where(components[:, i] < 1.e-12, 1, 0)))
479 24 30 918 8 28 266 217
# Play back the reconstruction
# Reconstruct a spectrogram by the outer product of component k and its activation
D_k = components @ activations
# invert the stft after putting the phase back in
#y_k = librosa.istft(D_k * phase)
y_k = librosa.istft(D_k * my_phase)
# And playback
print('Full reconstruction')
IPython.display.Audio(data=y_k, rate=sr)
Full reconstruction
components_ = np.zeros(components.shape)
components_[:, 0] = components[:, 0]
components_[:, 2] = components[:, 2]
components_[:, 6] = components[:, 6]
D_k = components_ @ activations
# invert the stft after putting the phase back in
y_k = librosa.istft(D_k * my_phase)
# And playback
print('Partly reconstruction')
IPython.display.Audio(data=y_k, rate=sr)
Partly reconstruction