FrontPage/Python/Scipy+Numpy の編集

** numpy [#pbad2eee]

[[→ndarray:rest-term>http://rest-term.com/archives/2999/]]

** ベクトライズによる高速化 [#vdd41bc9]
http://postd.cc/python_meets_julia_micro_performance/
いったん配列を取得してから、それをループ処理し、要素1ずつにhex()とint()関数を適用するのは一見無駄なことのようにも見えます。ありがたいことに、Numpyはループというより、むしろ配列に対し関数を呼び出す方法を提供します。すなわちnumpy.vectorize()関数です。この関数は、一度に1つのオブジェクトを操作する関数を入力としてとります。これは配列に作用する新しい関数を返します。

vhex = np.vectorize(hex)
vint = np.vectorize(int)
def parse_int_numpy():
    n = np.random.randint(0,2**32-1,1000)
    s = vhex(n)
    m = vint(s,16)
    np.all(m == n)
    return s

** (重要) 配列操作の注意! [#w3deae23]
 a = [0,1,2,3,4,5]
 In[8]: a
 Out[8]: [0, 1, 2, 3, 4]
 In[9]: a[0:3]
 Out[9]: [0, 1, 2]
 
 window = 100
 start = 0
 end = 99

のとき

b = x[start : end+1] 
 もしくは
 b = x[start : start+window]

In[10]: window=100
 In[11]: start=0
 In[12]: end = 99
 In[13]: x = range(200)
 In[14]: b = x[start:end+1]
 In[15]: len(b)
 Out[15]: 100
 In[16]: b = x[start:start+window]
 In[17]: len(b) 
 Out[17]: 100

** インストール [#ue151d9b]
アップデートの際にはまったので備忘録を残します。numpyをアップデートしようとした際にMingw32Compilerのなんちゃらでエラーが発生。
 pip install numpy -U
アンインストール後、再度インストールしようとしてもエラー。
 pip install numpy
そこで, [[こちらから>http://www.lfd.uci.edu/~gohlke/pythonlibs/#numpy]]numpy‑1.9.3+mkl‑cp27‑none‑win_amd64.whlをダウンロード。適当なフォルダにおいて、'shift+右クリック'>'コマンドラインでひらく'をして,
 pip install 'numpy-1.9.3*****.whl'
.whlからインストールすると一発で上手くいった。良かった。。。
** ログの保存 [#l110da97]
 def debuglog(name='', data=[]):
    with open('log.csv', 'a') as f:
        writer = csv.writer(f, lineterminator='\n')
        writer.writerow(name)
        writer.writerow(data)

** ログの保存２ [#y65afdea]
 from datetime import date
 import csv
 def debuglog(s=None, data=None):
    d = date.today().isoformat()
    with open('log.csv', 'a') as f:
        writer = csv.writer(f, lineterminator='\n')
        if not s is None:
            writer.writerow([d , s])
        if not data is None:
            writer.writerow([d, data])
        print s, data
** バイナリデータの読み書き [#gd88ed0e]
http://hydrocul.github.io/wiki/numpy/ndarray-io.html

*** np.save, np.load [#f95901d9]
 ndarr1 = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
 np.save('test.npy', ndarr1)
 ndarr2 = np.load('test.npy')

*** np.savez, np.load [#j534b0de]
np.savez を使うと複数のndarrayを名前付きで、しかも圧縮してファイルに保存できる。保存したファイルを読み込むには上と同じく np.load を使う。 np.load はファイルの保存形式を拡張子で判断しているようで、 np.savez での保存形式の場合は拡張子は .npz を使う。

ndarr1 = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
 ndarr2 = np.array([10, 20, 30])
 np.savez('test.npz', x=ndarr1, y=ndarr2)
 ndarr3 = np.load('test.npz')
 print(ndarr3['x'])
 # 出力結果
 # [ 0.   0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1. ]
 print(ndarr3['y'])
 # 出力結果
 # [10 20 30]

*** np.savetxt, np.loadtxt [#v90400eb]
2次元ndarryをテキスト形式で保存・読み込み. 区切り文字は空白がデフォルトなのでcsvにするためにはdelimiterオプションを利用

ndarr1 = np.array([[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]])
 np.savetxt('test.csv', ndarr1, delimiter=',')
 ndarr2 = np.loadtxt('test.csv', delimiter=',')
 print(ndarr2)
 # 出力結果
 # [[ 1.  2.  3.]
 #  [ 4.  5.  6.]
 #  [ 7.  8.  9.]]

** 初めての信号処理 [#a94b925d]

人工信号に対して、FFT、STFT解析を行うサンプルコード

[[→初めての信号処理>https://gist.github.com/peace098beat/b233857e44b1abe3bebf]]

** オーディオのロード [#sac51b33]
+ scipyを使う場合

返り値のdataはnumpy.array()

# オーディオのロード
    from scipy.io.wavfile import read
    fs, data = read(filename)
    data = data[0:fs/2-1]

** Numpy Array (np.ndarray) [#wdb3c394]
 ary = np.array([1,2,3,4,5]) = np.array( (1,2,3,4,5))
 # ゼロベクトル
 ary = np.zeros(3) # 1次元の場合はスカラーでよい
 ary = np.zeros((3, 2)) # N次元の場合はタプルを忘れずに
 # 空配列
 ary = np.empty((Nx,Ny))
 ary = np.empty_like(ary0)
 # 1ベクトル
 ary = np.ones((Nx,Ny))
 ary = np.ones_like(ary0)
 # 単位行列
 ary = identity(N) #NxNの単位行列が生成される

- ndarrayの属性・プロパティ

ary = np.array([1,2,3])
 # データタイプ(np.bool, np.int, np.float, np.complex... np.int32, np.float64...)
  ary.dtype 
 # タイプの変更
 ary = ary.astype(np.float)
 # 配列サイズ
 ary = np.empty((3,5,2), dtype = np.float)
 ary.shape
 >> (3,5,2)
 # 配列次元
 ary = np.empty((3,5,2), dtype = np.float)
 ary.ndim
 >> 3
 # 要素数
 ary = np.empty((3,5,2), dtype = np.float)
 ary.size
 >> 30
 
** Numpy Array (ndarray) 便利機能 [#q7969d99]

[[→The N-dimensional array>http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html]]

#ref(s-ndarray.jpg)
[[→ndarray atribute and method>https://gist.github.com/peace098beat/9eb9270c97060ae0c48f]]

** Numpy dtype [#g7d9e97b]

dtype_range = {np.bool_: (False, True),
               np.bool8: (False, True),
               np.uint8: (0, 255),
               np.uint16: (0, 65535),
               np.int8: (-128, 127),
               np.int16: (-32768, 32767),
               np.int64: (-2**63, 2**63 - 1),
               np.uint64: (0, 2**64 - 1),
               np.int32: (-2**31, 2**31 - 1),
               np.uint32: (0, 2**32 - 1),
               np.float32: (-1, 1),
               np.float64: (-1, 1)}

** オーディオI/O [#ra14cb50]
pythonでオーディオを再生するにはpyaudioがよさそう。pyaudioに再生させるためには、バイナリに戻す必要がある。

[[→numpy.unpackbit>http://docs.scipy.org/doc/numpy/reference/generated/numpy.unpackbits.html]]

>>> a = np.array([[2], [7], [23]], dtype=np.uint8)
 >>> a
 array([[ 2],
       [ 7],
       [23]], dtype=uint8)
 >>> b = np.unpackbits(a, axis=1)
 >>> b
 array([[0, 0, 0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 1, 1, 1],
       [0, 0, 0, 1, 0, 1, 1, 1]], dtype=uint8)

** np.frombufferで24bitファイルを読み込む [#u09b00d3]

[[→Speed up loading 24-bit binary data into 16-bit numpy array>http://stackoverflow.com/questions/22128961/speed-up-loading-24-bit-binary-data-into-16-bit-numpy-array]]

output = np.frombuffer(data,'b').reshape(-1,3)[:,1:].flatten().view('i2')

** sndhdr — サウンドファイルの識別 [#y37af40a]

path = './audio/sin_44100_24bit_stereo_5s.wav'
    print sndhdr.what(path)
    >>> ('wav', 44100, 2, -1, 24)
    print sndhdr.whathdr(path)
    >>>('wav', 44100, 2, -1, 24)

import sndhdr, wave, struct
 if sndhdr.what(fname)[0] != 'wav'
  raise StandardError("file doesn't have wav header")
 try:
  wav = wave.open(fname)
  params = (nchannels,sampwidth,rate,nframes,comp,compname) =  wav.getparams()
  frames = wav.readframes(nframes*nchannels)
 finally:
  wav.close()
 out = struct.unpack_from ("%dh" % nframes*nchannels, frames)

** 配列の拡張(行列の各列ベクトルを標準化) [#t268cbfc]
** 配列の拡張(行列の各列ベクトルを標準化) [#n529e4d8]

X = np.arange(12).reshape((3,4))
 >> array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
 // 平均ベクトル
 m = np.mean(X, axis=0)
 >> array([ 4.,  5.,  6.,  7.])
 // X-m そのまま各行ごとに演算できる
 >> array([
      [-4., -4., -4., -4.],
       [ 0.,  0.,  0.,  0.],
       [ 4.,  4.,  4.,  4.]])
 //標準偏差
 s = np.std(X, axis=0)
 >> array([
       [-0.3381028 , -1.20939778,  0.5124003 ,  1.35738762],
       [-2.12531859,  0.60120223,  1.75273089,  0.4617148 ],
       [-0.14482128, -0.73401208,  0.63781457, -0.33835372]])
 // 標準化
 // (X-m)/s
 >> array([
       [-1.22474487, -1.22474487, -1.22474487, -1.22474487],
       [ 0.        ,  0.        ,  0.        ,  0.        ],
       [ 1.22474487,  1.22474487,  1.22474487,  1.22474487]])

なんて便利なんだ。

** ** 配列の拡張(行列の各行ベクトルを標準化) [#y6f5308c]
** ** 配列の拡張(行列の各行ベクトルを標準化) [#z84ef687]
ポイントはm = np.mean(X, axis=1)[:, np.newaxis]
np.newaxisは軸を追加する。これによりnumpyが自動で軸ごとに演算してくれる。[[ブロードキャスト>http://seesaawiki.jp/met-python/d/array#content_4_12]]というらしい。

X=np.arange(12).reshape(3,4)
 [[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]
 m = np.mean(X, axis=1)
 [ 1.5  5.5  9.5]
 m.T[:,np.newaxis]
  [[ 1.5]
  [ 5.5]
  [ 9.5]]
 X-m
 [[-1.5 -0.5  0.5  1.5]
  [-1.5 -0.5  0.5  1.5]
  [-1.5 -0.5  0.5  1.5]]

** # numpy配列をプリントした場合の表示桁数 [#r9594de7]
 np.set_printoptions(precision=3)

def ndprint(a, format_string='{0:.2f}'):
    """
    ndarrayをprintする関数
    :example: ndprint(x)
    """
    return [format_string.format(v, i) for i, v in enumerate(a)]

def ndprints(s, a, format_string='{0:.2f}'):
    """
    ndarrayをprintする関数
    :example: ndprint(x)
    """
    print s, [format_string.format(v, i) for i, v in enumerate(a)]

** 配列の結合 [#h5fb7178]
 >>> a = np.arange(9).reshape((3,3))
 >>> a
 array([[0, 1, 2],
        [3, 4, 5],
        [6, 7, 8]])
 >>> b = np.arange(8,-1,-1).reshape((3,3))
 >>> b
 array([[8, 7, 6],
        [5, 4, 3],
        [2, 1, 0]])
  
 ## numpy.dstack() で二次元配列を結合して三次元配列にする
 >>> np.dstack((a,b))
 array([[[0, 8],
         [1, 7],
         [2, 6]],
  
        [[3, 5],
         [4, 4],
         [5, 3]],
  
        [[6, 2],
         [7, 1],
         [8, 0]]])
  
 ## numpy.hstack() で列方向に結合
 >>> np.hstack((a,b))
 array([[0, 1, 2, 8, 7, 6],
        [3, 4, 5, 5, 4, 3],
        [6, 7, 8, 2, 1, 0]])
  
 ## numpy.vstack() で行方向に結合
 >>> np.vstack((a,b))
 array([[0, 1, 2],
        [3, 4, 5],
        [6, 7, 8],
        [8, 7, 6],
        [5, 4, 3],
        [2, 1, 0]])
  
 ## numpy.dsplit() で三次元配列を分割
 >>> a = np.arange(16).reshape(2,2,4)
 >>> a
 array([[[ 0,  1,  2,  3],
         [ 4,  5,  6,  7]],
  
        [[ 8,  9, 10, 11],
         [12, 13, 14, 15]]])
  
 >>> np.dsplit(a,2)
 [array([[[ 0,  1],
         [ 4,  5]],
  
        [[ 8,  9],
         [12, 13]]]), array([[[ 2,  3],
         [ 6,  7]],
  
        [[10, 11],
         [14, 15]]])]
  
 ## numpy.hsplit() で列方向に分割
 >>> a = np.arange(16).reshape(4,4)
 >>> a
 array([[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11],
        [12, 13, 14, 15]])
 >>> np.hsplit(a,2)
 [array([[ 0,  1],
        [ 4,  5],
        [ 8,  9],
        [12, 13]]), array([[ 2,  3],
        [ 6,  7],
        [10, 11],
        [14, 15]])]
  
 ## numpy.vsplit() で行方向に分割
 >>> np.vsplit(a,2)
 [array([[0, 1, 2, 3],
        [4, 5, 6, 7]]), array([[ 8,  9, 10, 11],
        [12, 13, 14, 15]])]
  
 ## numpy.transpose() で配列を転置
 >>> a = np.array([[1, 2], [3, 4]])
 >>> a
 array([[1, 2],
        [3, 4]])
 >>> np.transpose(a)
 array([[1, 3],
        [2, 4]])
 ## ndarray.T でも良い
 >>> a.T
 array([[1, 3],
        [2, 4]])
        
 ## numpy.swapaxes() で軸の交換
 >>> a = np.array([[1,2,3]])
 >>> np.swapaxes(a, 0, 1)
 array([[1],
        [2],
        [3]])

タイムスタンプを変更しない