Travis CIによるNikolaブログ構築の自動化

Travis CIによるNikolaブログ構築の自動化

Nikolaでブログを構築するための方法は過去の記事(1, 2, 3)に書いていますが、手元でのビルドのためにpythonやnikolaがインストールされたPCが必要になるため、ブログ記事を書くための環境が限定されてしまうという問題がありました。 そこで、この記事ではNikola公式の記事を参考に、Travis CIを用いることでsrcブランチにブログ記事をコミットするだけでTravis CI側で自動的にブログを構築し、masterブランチにプッシュしてくれる仕組みを構築する過程を紹介します。 なお、以下の作業はすべてsrcブランチ上で行います。

conf.pyの編集

まず、nikola github_deployコマンド実行時に、デプロイと同時にsrcブランチもコミットするかどうかを制御するためのオプションをFalseに設定します。

GITHUB_COMMIT_SOURCE = False

これはTravis CIではソースブランチへのコミットをトリガーにnikola build && nikola github_deployを実行するため、nikola github_deployによってsrcブランチがコミットされると再度それがトリガーとなって無限ループに陥ってしまうことを防ぐためです。

.travis.ymlの作成

nikolaブログのルートディレクトリ(conf.pyが置いてあるディレクトリ)に以下のような.travis.ymlファイルを作成します。 これはNikola公式の記事に記載のtravis.ymlを改変したものになります。

※2019/4/8追記, travis CIでのビルドに失敗する対策のissueの内容を反映

language: python
cache: apt
sudo: false
addons:
  apt:
    packages:
    - language-pack-ja-base
    - language-pack-ja
branches:
  only:
  - src
python:
- 3.6
before_install:
- git config --global user.name 'USERNAME'
- git config --global user.email 'travis@invalid'
- git config --global push.default 'simple'
- pip install --upgrade pip wheel
- echo -e 'Host github.com\n    StrictHostKeyChecking no' >> ~/.ssh/config
- eval "$(ssh-agent -s)"
- chmod 600 id_rsa
- ssh-add id_rsa
- git remote rm origin
- git remote add origin git@github.com:USERNAME/REPO.git
- git fetch origin master
- git branch master FETCH_HEAD
install:
- pip install 'ghp-import2'
- pip install 'webassets'
- pip install -U --upgrade-strategy=eager 'Nikola[extras]'
script:
- nikola build && nikola github_deploy -m 'Nikola auto deploy [ci skip]'
notifications:
  email:
    on_success: change
    on_failure: always

ここで、

- git config --global user.name 'USERNAME'
- git config --global user.email 'travis@invalid'

の行は適切なユーザー名、メールアドレスに、

- git remote add origin git@github.com:USERNAME/REPO.git

の行は適切なユーザー名およびリポジトリ名に変更する必要があります。

SSH鍵の生成

まず、.gitignoreファイルにid_rsaid_rsa.pubを無視する設定を追記した上でSSH鍵を生成します。

$ echo id_rsa >> .gitignore
$ echo id_rsa.pub >> .gitignore
$ ssh-keygen -C TravisCI -f id_rsa -N ''

上記を実行すると秘密鍵id_rsaおよび公開鍵id_rsa.pubが作成されます。 念のため、.gitignoreの無視設定が合っているかを確認するためにgit statusを実行してもaddされる対象となっていないことを確認しておきましょう。 なおssh-keygen-Cはコメント、-fは鍵名、-Nはパスフレーズの指定(ここでは空文字なので指定なし)です。

公開鍵のgithubリポジトリへの登録

生成した公開鍵はgithubリポジトリに登録しておく必要があります。 リポジトリページ -> Settings -> Deploy Keys -> Add deploy key からTitleをTravis CIとし、Keyにidrsa.pubの中身をコピペして登録しましょう。 また、Allow write accessはチェックしておく必要があります。 これらの作業を忘れるとTravis CIの自動ビルド時にアクセスエラーが発生します。

travis encrypt-fileコマンドによる秘密鍵の暗号化

生成した秘密鍵id_rsatravis encrypt-fileによってid_rsa.encに共通鍵による暗号化を施した上でリポジトリに追加します。 この暗号化されたid_rsa.encはTravis CIでの自動ビルド時に共通鍵によって復号化され、Travis CIのgithubへのアクセスに使用されます。

この作業を実行するためにはtravisコマンドが必要となりますが、これはrubyのgemとして配布されているため、インストールされていない場合は下記にようにgemでインストールする必要があります。

$ gem install --user-install travis

なおtravisコマンドがインストールされる場所にPATHが通っていない場合はPATHに追加するか、下記一連のコマンドをフルパスで実行する必要があります。

travisコマンドが実行可能になったら、更に下記を実行します。

$ travis login
$ travis enable
$ travis encrypt-file id_rsa --add

travis loginを実行するとgithubアカウントのユーザ名、パスワードを求められるため入力してください。 この上でtravis enableを実行すると自動ビルドを有効化するリポジトリが正しいか確認されるため、正しければyesと入力しましょう。 更にtravis encrypt-file id_rsa --addを実行すると、秘密鍵id_rsaが暗号化されてid_rsa.encが生成されます。 このid_rsa.enc.gitignoreに追加されていないため、git addによってgitの管理下に置かれることになります。 なお引数の--addを付けておくと、.travis.ymlファイルにid_rsa.encの復号化を行うための下記のようなopensslのコマンドを追加してくれます。

before_install:
- openssl aes-256-cbc -K $encrypted_XXXXXX_key -iv $encrypted_XXXXXX_iv
  -in id_rsa.enc -out id_rsa -d

-Kは共通鍵、-ivは初期ベクトルの指定であり、指定されている値はTravis CI側で環境変数として設定されています。 (ブラウザからTravis CIの設定を見ると確認することができます)

srcブランチへの各種ファイルのadd, commitおよびgithubへのpush

以上の作業によりconf.py.gitignore.travis.ymlid_rsa.encの4つのファイルが編集・生成されているため、これをsrcブランチにaddした上でcommitし、更にremoteとなっているgithubにpushします。

$ git add .
$ git commit -am "Automate builds with Travis CI"
$ git push origin src

これによりTravis CIの自動ビルドが実行されるはずなので、あとはブラウザからTravis CIのページを確認し、ビルドが通っているかを確認するのが良いでしょう。

上記設定以降の記事の追加方法

上記までで設定した方法によってsrcブランチが変更される度にTravis CIがnikola buildおよびnikola github_deployを行ってくれるようになったため、postsディレクトリ以下に新規記事を追加したら、あとはこの新規記事をsrcブランチにgit add posts/XX.mdのようにgit addした上でgit commitして、更にgithubにgit push origin srcすれば自動的にブログがビルドされます。

参考

2018/5/19にPycon Osakaで'librosaで始める音楽情報処理'というタイトルで発表しました

2018/5/19にPyCon mini Osakaで「librosaで始める音楽情報検索」というタイトルで発表しました

公開が遅くなりましたが先月開催されたPyCon mini Osakaで発表した内容を公開します。この記事もAzure Notebooksで既に公開済ではありましたが、nikolaのipynb公開機能を使ってブログポストとして投稿しておきたいと思います。

librosaで始める音楽情報検索

  • 大橋 宏正(twitter: @wrist)

この資料

自己紹介

  • 大橋宏正(@wrist)
    • 某メーカで音響信号処理の業務に従事
    • 2013年ごろから大阪PRML読書会を主催(現在休止中)
    • PyData Osakaオーガナイザ

PyCon mini Osakaのテーマ

  • 「Pythonでなにしよう?」
    • Webアプリ
    • 科学技術計算
  • Pythonで音楽の分析を行いたい!

この発表

  • Pythonで音楽情報検索(MIR)を行うためのライブラリlibrosaの紹介

音楽情報検索(MIR)とは?

  • Music Information Retrievalの略
  • 音楽情報に関する処理全般
    • ビートトラッキング、楽器音分離、音楽推薦など
  • コミュニティ
    • ismir.net
    • 音楽情報科学研究会(www.sigmus.jp)

MIRをPythonで行う

  • Pythonで音楽データの処理が必要
  • 読み書き
    • wave, audioread, pysoundfile, scipy.io, ...
  • 線形代数、信号処理、フーリエ変換
    • numpy, scipy(signal, fftpack, ...)
  • 可視化
    • matplotlib, seaborn
  • これらを組み合わせることでさまざまな処理を実現

参考

音響信号処理についてはQiitaに過去に諸々記事を書いてます

librosa

librosaの特徴

  • MIRに必要となる典型的な処理が多く実装
  • 低い学習障壁
  • 一貫性、後方互換性、モジュール性
In [1]:
# pipでインストール
!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: 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: 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: 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: 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: 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: 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: 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: 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)

librosaの具体的な紹介

  • 適宜デモを交えながらlibrosaを紹介
  • 発表の流れ
    • 音データの短時間周波数分析
    • 特徴量抽出
    • オンセット検出、ビートトラッキング
In [2]:
# 本発表では以下を仮定
import numpy as np
import scipy as sp
import scipy.signal as sg
import scipy.fftpack as fft
import matplotlib as mpl
import matplotlib.pyplot as plt
import IPython.display
In [3]:
# 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
In [4]:
# librosa以下の各モジュールをimport
import librosa
import librosa.display
import librosa.feature
import librosa.effects
import librosa.onset
import librosa.beat

サウンドデータの読み込み

  • librosa.load
    • ファイルパスを与えることでデータを読み込み
    • デフォルトでモノラル化、サンプリングレート22050Hzとなる
    • バックエンドにaudioreadライブラリを使用
  • librosa.display.waveplotで波形データの可視化
In [5]:
# デモ音源のロード&再生
# y: モノラル音声データ, sr: サンプリングレート
y, sr = librosa.load(librosa.util.example_audio_file())
IPython.display.Audio(data=y, rate=sr)
Out[5]:
In [6]:
# 波形の描画
# 縦軸レンジを自動的に指定し、横軸を秒数表示にしてくれる
librosa.display.waveplot(y, sr)
Out[6]:
<matplotlib.collections.PolyCollection at 0x115239668>

サウンドデータのデジタル表現

  • 現実の音はアナログデータ
    • 時間方向にも振幅方向にも連続
  • コンピュータで音を取り扱うためには離散化が必要
    • 時間方向: サンプリング
    • 振幅方向: 量子化
In [7]:
b_pt = int(10.0*sr); e_pt = int((10.0+0.002)*sr)
fig = plt.figure(1, figsize=(14,6)); ax = fig.add_subplot(111)
ax.plot(y[b_pt:e_pt]); ax.plot(y[b_pt:e_pt], "o"); ax.grid();
ax.set_xlabel("time [pt]"); ax.set_ylabel("amplitude")
Out[7]:
Text(0,0.5,'amplitude')
In [8]:
# サンプリングレート: 1秒辺りのサンプル数
print("sampling rate: {0}".format(sr))
print("length: {0}[pt]={1}[s]".format(y.shape[0], y.shape[0]/sr))
sampling rate: 22050
length: 1355168[pt]=61.45886621315193[s]

音声データの周波数分析

  • 音声データの長さはデータにより異なる
  • 固定長の短時間波形に切り出した上で処理を行う
  • 短時間波形に対し離散フーリエ変換を適用することで周波数分析を行う
In [9]:
# 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
In [10]:
def plot_overlap_waves():
    import matplotlib.gridspec as gridspec

    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    gs = gridspec.GridSpec(3, nframes + 1)
    
    ax_wav = fig.add_subplot(gs[0, :])
    ax_wav.plot(np.arange(half_width*(nframes+1))/sr, 
                y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)/sr])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    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)[:nframes]):
        y_win = win * y_cut
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width)/sr, 0.2 * win)
        ax = None
        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:])
            ax.set_xlim([0, frame_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_xlim([0, frame_width])

        ax.tick_params(labelleft="off",left="off") # y軸の削除
        ax.set_ylim([-ymax, ymax]); ax.grid()

    fig.subplots_adjust(wspace=0.0)
In [11]:
# 周波数分析時には通常窓関数を乗算した短時間波形を使用
# 切り出しはここでは1/2区間をオーバーラップさせている
plot_overlap_waves()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

用語の整理

  • sr = sampling rate [Hz = pt/s]
    • 22050 [Hz]
  • frame = 各短時間波形
  • n_fft = 各フレームに含まれるサンプル数
    • 2048 [pt]
  • hop_length = フレーム間のサンプル数
    • 512 [pt]
In [12]:
def plot_overlap_freqs():
    import matplotlib.gridspec as gridspec

    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    gs = gridspec.GridSpec(3, nframes + 1)
    
    ax_wav = fig.add_subplot(gs[0, :])
    ax_wav.plot(y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)
    w = np.linspace(0.0, 2.0 * np.pi, frame_width, endpoint=False)
    freqs = w[:frame_width//2+1] * (sr / 2) / np.pi
    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
        y_win = win * y_cut
        Y_abs = np.abs(fft.fft(y_win, frame_width)[:frame_width//2+1])
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width), 0.2 * win)
        ax = None
        if i % 2 == 0:
            ax = fig.add_subplot(gs[1, i:i+2])
            ax.semilogy(Y_abs, freqs)
        else:
            ax = fig.add_subplot(gs[2, i:i+2])
            ax.semilogy(Y_abs, freqs)

        if i > 1:
            ax.tick_params(labelleft="off",left="off") # y軸の削除
        else:
            ax.tick_params(labelbottom="off",bottom="off") # x軸の削除
            ax.set_ylabel("frequency [Hz]")

        ax.set_xlabel("amplitude")
        ax.set_ylim([10, sr/2])
        ax.grid()

    fig.subplots_adjust(wspace=0.0)
In [13]:
# 各短時間波形を周波数分析したグラフを描画
# 縦軸が周波数、横軸が振幅スペクトル
plot_overlap_freqs()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [14]:
def plot_spectrogram():
    frame_width = 2048
    half_width = int(frame_width // 2)
    nframes = 20

    fig = plt.figure(1, figsize=(12,8))
    
    ax_wav = fig.add_subplot(211)
    ax_wav.plot(np.arange(half_width*(nframes+1))/sr,
                y[:half_width*(nframes+1)])
    ymax = np.max(y[:half_width*(nframes+1)])
    ax_wav.set_xlim([0, half_width*(nframes+1)/sr])
    ax_wav.set_ylim([-ymax, ymax]); ax_wav.grid()

    win = sg.get_window("hanning", int(frame_width))
    xs = np.arange(frame_width)
    w = np.linspace(0.0, 2.0 * np.pi, frame_width, endpoint=False)
    freqs = w[:frame_width//2+1] * (sr / 2) / np.pi

    Zs = []
    for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
        y_win = win * y_cut
        Y_abs = np.abs(fft.fft(y_win, frame_width)[:frame_width//2+1])
        ax_wav.plot(np.arange(half_width*i, half_width*i + frame_width), 0.2 * win)
        Zs.append(20.0*np.log10(Y_abs))
    
    ax = fig.add_subplot(212)
    Xs = np.arange((nframes+1)) / sr / 2
    pcol = ax.pcolormesh(Xs, freqs, np.array(Zs).T, cmap='magma')
    cax = fig.colorbar(pcol)
    cax.set_clim(-80, 20)

    fig.subplots_adjust(wspace=0.0)
In [15]:
# 短時間波形を周波数分析したものをスペクトログラムとして表示
plot_spectrogram()

librosaにおける短時間周波数分析

  • 短時間離散フーリエ変換(Short time fourie transform)
    • librosa.stft
  • Constant-Q Transform(CQT)
    • librosa.cqt
  • 分析結果の可視化
    • librosa.specshow
In [16]:
# librosaでは短時間波形分析をメソッド一発で実行可能
D = librosa.stft(y)
In [17]:
# D[k, l]
# k: 周波数ビン  l: フレームインデックス
# kは実信号に対する2048ptフーリエ変換の結果なので1025ptのみ保持
D.shape
Out[17]:
(1025, 2647)

離散フーリエ変換$D[k,l]$の算出法

$$D[k, l]=\sum_{n=0}^{{\textrm N_{\rm fft}-1}} w[n]\cdot x\left[\left(l\cdot {\rm M}-\frac{N_{\rm fft}}{2}\right)+n\right]\exp\left\{-j\frac{2\pi k n}{N_{\rm fft}}\right\}$$
* $x[n]$は入力信号、$w[n]$は窓関数
* 周波数ビン$k$、フレームインデックス$l$、$M$はhop length
* $lM$を中心に前後$N/2$点に渡り離散フーリエ変換を実行
In [18]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_power = librosa.amplitude_to_db(np.abs(D), ref=np.max)
librosa.display.specshow(log_power, x_axis="time", y_axis="linear")
plt.colorbar()
Out[18]:
<matplotlib.colorbar.Colorbar at 0x115a2aba8>
In [19]:
# 縦軸を対数表示
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_power = librosa.amplitude_to_db(np.abs(D), ref=np.max)
librosa.display.specshow(log_power, x_axis="time", y_axis="log")
plt.colorbar()
Out[19]:
<matplotlib.colorbar.Colorbar at 0x1c3344cef0>

Constant-Q変換

$$C[k]= \frac{1}{N[k]} \sum_{n=0}^{{\textrm N[k]-1}} w[n, k]\cdot x\left[n\right]\exp\left\{-j\frac{2\pi Q[k] n}{N[k]}\right\}$$
  • 特徴
    • 分析する周波数間隔を対数的に配置
    • 低域は細かく分析し、高域は荒く分析
    • 音楽信号の分析に適する
In [20]:
C = librosa.cqt(y, sr)
In [21]:
# STFTと比較して縦軸の周波数データ数は圧倒的に少ない
# STFTは周波数間隔が等幅なのに対してCQTだと対数間隔となるため
print(D.shape)
print(C.shape)
(1025, 2647)
(84, 2647)
In [22]:
fmin = 32.70  # C1に相当するHz
fmax = fmin * ((2**(1/12)) ** 84)
print(fmax)
print(librosa.hz_to_note(fmin))
print(librosa.hz_to_note(fmax))
4185.600000000015
C1
C8
In [23]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
plt.colorbar()
Out[23]:
<matplotlib.colorbar.Colorbar at 0x115933cf8>
In [24]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_note")
plt.colorbar()
Out[24]:
<matplotlib.colorbar.Colorbar at 0x11589d518>
In [25]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

# デフォルトだと80dBまでだがtop_dbを付けるとレンジが縮まる
log_cqt_power = librosa.amplitude_to_db(np.abs(C), ref=np.max, top_db=40)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_note")
plt.colorbar()
Out[25]:
<matplotlib.colorbar.Colorbar at 0x1c1f823518>

クロマベクトル

  • 12音階それぞれの強さを要素に持つベクトル
  • CQTで算出した対数周波数帯域ごとの強度を元に計算
  • librosa.feature.chroma_cqt
In [26]:
chroma = librosa.feature.chroma_cqt(C=C, sr=sr)

fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

librosa.display.specshow(chroma, x_axis="time", y_axis="chroma")
plt.colorbar()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/librosa/display.py:657: UserWarning: Trying to display complex-valued input. Showing magnitude instead.
  warnings.warn('Trying to display complex-valued input. '
Out[26]:
<matplotlib.colorbar.Colorbar at 0x1159a3128>
In [27]:
# STFTの結果Dからクロマベクトルを算出
chroma_stft = librosa.feature.chroma_stft(S=D, sr=sr)

fig = plt.figure(1, figsize=(12,4))
ax = fig.add_subplot(2,1,1)
librosa.display.specshow(chroma_stft, y_axis='chroma')
plt.title('chroma_stft'); plt.colorbar()
ax= fig.add_subplot(2,1,2)
librosa.display.specshow(chroma, y_axis='chroma', x_axis='time')
plt.title('chroma_cqt'); plt.colorbar()
plt.tight_layout()
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/librosa/display.py:657: UserWarning: Trying to display complex-valued input. Showing magnitude instead.
  warnings.warn('Trying to display complex-valued input. '

その他特徴量抽出

  • メルスペクトル、MFCC、tonnetzなども算出可能
  • $\Delta$特徴量の計算も可能
In [28]:
M = librosa.feature.melspectrogram(y=y, sr=sr)
mfcc = librosa.feature.mfcc(y=y, sr=sr)
tonnetz = librosa.feature.tonnetz(y=y, sr=sr)

librosa.effectsモジュール

  • 音源分離
    • librosa.hpss
      • 調波成分と打音成分に分解
  • 時間伸縮、ピッチシフト
    • librosa.time_stretch
    • librosa.pitch_shift
In [29]:
y_harmonic, y_percussive = librosa.effects.hpss(y)
In [30]:
IPython.display.Audio(y, rate=sr)
Out[30]:
In [31]:
IPython.display.Audio(y_harmonic, rate=sr)
Out[31]:
In [32]:
IPython.display.Audio(y_percussive, rate=sr)
Out[32]:
In [33]:
# それぞれに対するCQTを求める
C      = librosa.cqt(y=y, sr=sr)
C_harm = librosa.cqt(y=y_harmonic, sr=sr)
C_perc = librosa.cqt(y=y_percussive, sr=sr)
In [34]:
fig = plt.figure(1, figsize=(12,4)); ax = fig.add_subplot(1,1,1)

log_cqt_power= librosa.amplitude_to_db(np.abs(C), ref=np.max)
log_cqt_harm = librosa.amplitude_to_db(np.abs(C_harm), ref=np.max(np.abs(C)))
log_cqt_perc = librosa.amplitude_to_db(np.abs(C_perc), ref=np.max(np.abs(C)))

ax = fig.add_subplot(3,1,1); librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
ax = fig.add_subplot(3,1,2); librosa.display.specshow(log_cqt_harm , x_axis="time", y_axis="cqt_hz")
ax = fig.add_subplot(3,1,3); librosa.display.specshow(log_cqt_perc , x_axis="time", y_axis="cqt_hz")
Out[34]:
<matplotlib.axes._subplots.AxesSubplot at 0x1139ae2e8>
In [35]:
# 2倍の速度に時間伸縮
y_fast = librosa.effects.time_stretch(y, 2.0)
IPython.display.Audio(y_fast, rate=sr)
Out[35]:
In [36]:
# 3度上にピッチシフト
y_third = librosa.effects.pitch_shift(y, sr, n_steps=4)
IPython.display.Audio(y_third, rate=sr)
Out[36]:
In [37]:
# 減五度にピッチシフト
y_tritone = librosa.effects.pitch_shift(y, sr, n_steps=-6)
IPython.display.Audio(y_tritone, rate=sr)
Out[37]:

onset検出とビートトラッキング

  • librosa.onsetモジュール、librosa.beatモジュールを使用
    • 出音の強度をonset.onset_strengthで算出
    • 強度を元にonset箇所をonset.onset_detectで検出
    • 強度を元にbeat.beat_trackでテンポとビートを推定
In [38]:
onset_envelope = librosa.onset.onset_strength(y, sr)
onsets = librosa.onset.onset_detect(onset_envelope=onset_envelope)
In [39]:
fig = plt.figure(1, figsize=(12,8))
ax = fig.add_subplot(2, 1, 1)
ax.plot(onset_envelope, label="onset strength")
ax.vlines(onsets, 0, np.max(onset_envelope), color="red", alpha=0.25, label="onsets")
ax.set_xlim([0, len(onset_envelope)]); ax.legend(loc=1)
ax = fig.add_subplot(2, 1, 2)
librosa.display.waveplot(y, sr)
Out[39]:
<matplotlib.collections.PolyCollection at 0x113aa6c88>
In [40]:
tempo, beats = librosa.beat.beat_track(onset_envelope=onset_envelope)
In [41]:
print(tempo)
129.19921875
In [42]:
fig = plt.figure(1, figsize=(12,8))
ax = fig.add_subplot(1, 1, 1)
ax.plot(onset_envelope, label="onset strength")
ax.vlines(onsets, 0, 0.8 * np.max(onset_envelope), color="red", alpha=0.5, linestyle="dotted", label="onsets")
ax.vlines(beats, 0, np.max(onset_envelope), color="black", alpha=0.25, label="beats")
ax.set_xlim([0, len(onset_envelope)]); ax.legend(loc=1)
Out[42]:
<matplotlib.legend.Legend at 0x114fc1128>

クリック音を生成し同時に再生

  • MIRの評価に用いられるmir_evalライブラリを使用
In [43]:
!pip install mir_eval
Requirement already satisfied: mir_eval in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.4)
Requirement already satisfied: six in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.11.0)
Requirement already satisfied: numpy>=1.7.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.14.2)
Requirement already satisfied: scipy>=0.9.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (1.1.0)
Requirement already satisfied: future in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from mir_eval) (0.16.0)
In [44]:
# http://colinraffel.com/publications/ismir2014mir_eval.pdf
import mir_eval
In [45]:
beat_times = librosa.frames_to_time(beats)
y_click = mir_eval.sonify.clicks(beat_times, sr, length=len(y))
In [46]:
IPython.display.Audio(y_click, rate=sr)
Out[46]:
In [47]:
IPython.display.Audio(y + y_click, rate=sr)
Out[47]:

特徴量のビート同期

  • librosa.util.syncを用いるとビートとその他特徴量を同期させることができる
    • 各ビートの間の成分を平均化したり中央値を取ったりできる
In [48]:
tempo, beat_frames = librosa.beat.beat_track(y=y_percussive, sr=sr)
# 13次元のMFCCおよびΔを抽出
mfcc = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)
mfcc_delta = librosa.feature.delta(mfcc)
In [49]:
print(beat_frames.shape)
print(mfcc.shape)
print(mfcc_delta.shape)
(98,)
(13, 2647)
(13, 2647)
In [50]:
# ビート間の平均値を用いてMFCCとΔMFCCをビートに同期
beat_mfcc_delta = librosa.util.sync(np.vstack([mfcc, mfcc_delta]),
                                    beat_frames)
In [51]:
beat_mfcc_delta.shape
Out[51]:
(26, 99)
In [52]:
# 調波成分からクロマベクトルを推定
chromagram = librosa.feature.chroma_cqt(y=y_harmonic, sr=sr)

# ビート間の中央値を用いてクロマベクトルをビートに同期
beat_chroma = librosa.util.sync(chromagram, beat_frames, aggregate=np.median)

# ビート同期特徴量(クロマベクトル、MFCC、ΔMFCC)を結合
beat_features = np.vstack([beat_chroma, beat_mfcc_delta])
In [53]:
# MFCC 13dim, ΔMFCC 13dim, クロマベクトル(12音階; 12dim)
beat_features.shape
Out[53]:
(38, 99)
In [54]:
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(211)
librosa.display.specshow(log_cqt_power, x_axis="time", y_axis="cqt_hz")
log_beat_features = librosa.amplitude_to_db(np.abs(beat_features), ref=np.max)
ax = fig.add_subplot(212)
ax.pcolormesh(log_beat_features)
Out[54]:
<matplotlib.collections.QuadMesh at 0x113a0ac50>

まとめ

  • librosaを用いた音楽情報検索(MIR)について説明
    • 特徴量抽出、可視化、音源分離、ビート推定
  • 音楽データに対する探索的分析の手段の一つ

以下補足およびおまけ

librosa.decomposeによる音源分離

  • NMFによるモノラル音源分離を行う
    • スペクトログラムを行列と見てscikit-learnのNMFに投げる
    • アクティベーション行列と基底行列に分解
    • 教師なし音源分離が可能

NMF

  • 非負値行列$\textbf{X}$を$\textbf{X} \simeq \textbf{WH} $と分解する手法
    • Hを基底行列、Wをアクティベーションと呼ぶ
  • scikit-learnのNMFでは下記目的関数を最小化
$$0.5 * ||\textbf{X} - \textbf{WH}||_{Fro}^2 + \alpha l1_{ratio} ||vec(W)||_1 + \alpha l1_{ratio} ||vec(H)||_1 + 0.5 \alpha (1 - l1_{ratio}) ||W||_{Fro}^2 + 0.5 \alpha (1 - l1_{ratio}) ||H||_{Fro}^2 $$
  • 通常最適化には補助関数法が使われる
In [55]:
D = librosa.stft(y)
S, phase = librosa.magphase(D)
comps, acts = librosa.decompose.decompose(S, n_components=8, sort=True)
S_approx = comps @ acts
In [56]:
print(comps.shape)
print(acts.shape)
print(S.shape, S_approx.shape)
(1025, 8)
(8, 2647)
(1025, 2647) (1025, 2647)
In [57]:
def plot_decomposed(S, comps, acts, S_approx):
    fig = plt.figure(figsize=(10,8))
    ax = fig.add_subplot(311)
    librosa.display.specshow(librosa.amplitude_to_db(S, ref=np.max), y_axis='log', x_axis='time')
    ax.set_title('Input spectrogram'); plt.colorbar(format='%+2.0f dB')

    ax = fig.add_subplot(323)
    librosa.display.specshow(librosa.amplitude_to_db(comps, ref=np.max), y_axis='log')
    plt.colorbar(format='%+2.0f dB'); ax.set_title('Components')

    ax = fig.add_subplot(324)
    librosa.display.specshow(acts, x_axis='time')
    ax.set_ylabel('Components'); ax.set_title('Activations')
    plt.colorbar()

    ax = fig.add_subplot(313)
    librosa.display.specshow(librosa.amplitude_to_db(S_approx, ref=np.max), y_axis='log', x_axis='time')
    plt.colorbar(format='%+2.0f dB')
    ax.set_title('Reconstructed spectrogram')

    fig.tight_layout()
In [58]:
plot_decomposed(S, comps, acts, S_approx)
In [59]:
S_partly = comps[:, :3] @ acts[:3, :]
In [60]:
plot_decomposed(S, comps[:, :3], acts[:3, :], S_partly)
In [61]:
y_partly = librosa.istft(S_partly * phase)
IPython.display.Audio(y_partly, rate=sr)
Out[61]:

decomposeの用途

  • 視聴用途の音源分離には(このままでは)向かない
    • 振幅スペクトルのみを分解するが、位相情報は元信号の元を使うため歪が大きい
    • 基底数を増加した上で適切なクラスタリングを施して使用したり、事前学習した楽器別の基底などを用いて分離を行うなどの工夫が必要
  • 音響イベント検出のための特徴量に用いることも可能

2018/2/28にPydata Osakaで'SciPyの概要と各モジュールの紹介'というタイトルで発表しました

2018/2/28にPydata Osakaで発表しました

非常に今更ですが2018/2/28に"SciPyの概要と各モジュールの紹介"というタイトルで発表しました。当日の資料は下記途中のセルにもあるようにAzure Notebooksでも公開していましたが、nikolaではipynbファイルを記事としても公開可能であるため今更にはなりますが公開いたしました。 以下当日の資料になります。

SciPyの概要と各モジュールの紹介

自己紹介

  • 大橋宏正
    • 某メーカで音響エンジニア
    • scipy.signalとscipy.fftpackを仕事でよく使う
  • twitter: @wrist

SciPyとは

  • NumPyの上に作られた数学的なアルゴリズムや便利な関数の集合
  • 機能ごとに分けられたサブモジュール群から成る

本発表で紹介するサブモジュール(1)

  • io
    • 様々なファイル入出力
  • special
    • 特殊関数
  • fftpack
    • FFTルーチン
  • interpolate
    • 補間と平滑化
  • linalg
    • 線形代数

サブモジュール一覧(2)

  • optimize
    • 最適化および求根ルーチン
  • signal
    • 信号処理
  • spatial
    • 空間データ構造およびアルゴリズム
  • stats
    • 確率分布および密度関数
  • misc
    • 雑多な諸々

本日説明しないモジュール

  • cluster
    • クラスタリングアルゴリズム
  • constants
    • 科学定数、数学定数
  • odr
    • 直交距離回帰を行うためのFortranライブラリODRPACKのラッパー
  • ndimage
    • N次元画像処理
  • sparse
    • 疎行列および関連ルーチン
  • integrate
    • 積分および常微分方程式のソルバ

SciPy 1.0

History

  • 2011: GitHubに移行
  • 2011: Python3サポート
  • 2012: sparse graphモジュールと統合的な最適化インタフェースの追加
  • 2012: scipy.maxentropyの削除
  • 2013: TravisCIでの継続的インテグレーションを開始
  • 2015: BLAS/LAPACKに対するCythonインタフェースとベンチマークスイートの追加
  • 2017: scipy.LowLevelCallableを用いた統一的なC APIの追加; およびscipy.weaveの削除
  • 2017: SciPy 1.0のリリース
    • version 0.1から実に16年

何故Version 1.0か?(意訳)

  • 大枠の背景
    • Scipyはプロダクション環境で使われるほどに成熟し安定したことを反映
    • 技術的な目標(WindowsでのWheel配布、CI)と組織的な目標(組織構造、行動規範、ロードマップ)の両方を近年達成
  • 1.0の意味
    • オープンソースプロジェクトにとって"1.0"と言うと100%完成したという印象があるが、実際は完璧ではなくいくつもの修正が必要な点があることを自分たちとしても認めている
    • そうだとしてもSciPyは究極的にユーザーにとって便利なものであり、平均して高いクオリティのコードとドキュメントがあり、1.0のラベルが示すであろう安定性と後方互換性の保証がある。

SciPyを学ぶ際のリソース

scipy lecture note: 1.5 Scipy: high-level scientific computing

  • http://www.scipy-lectures.org/intro/scipy.html
  • scipyはさまざまな科学技術計算で共通する問題に対するtoolbox
    • 補間、積分、最適化、画像処理、統計、特殊関数、など
  • GSLやmatlabと比較されうるもの

本発表

  • tutorialとscipy lecture noteからコード等を抜粋して紹介
  • scipyを使った計算の簡単なデモ
In [1]:
# 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)
In [2]:
# 本発表では以下を仮定
import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt
In [3]:
# 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の代わりに@演算子を使用

In [4]:
# 本発表では以下を仮定
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

scipy.io

  • io
    • さまざまな形式のファイルIOを提供するモジュール
      • matlabやoctaveのmatファイル
      • netcdf, fortranのunformatted seaquential file, matrix market, harwell-boeing, IDL, Arff形式
    • 何故かwavファイルのIOもできるが24bit音源が読めないのでハイレゾ音源が扱えない
In [5]:
a = np.ones((3, 3))
spio.savemat('file.mat', {'a': a}) # savematには辞書を渡す
data = spio.loadmat('file.mat')
data['a']
Out[5]:
array([[1., 1., 1.],
       [1., 1., 1.],
       [1., 1., 1.]])
In [6]:
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

  • さまざまな特殊関数が含まれている
    • docstringに色々書かれている
  • よく使われると思われるもの
    • ベッセル関数 scipy.special.jn()
    • 楕円関数 scipy.special.ellpj()
    • ガンマ関数 scipy.special.gamma()、対数ガンマ関数 scipy.special.gammaln()
    • 誤差関数 scipy.special.erf()
In [7]:
special?

scipy.special.cython_special

  • cythonバインディングも存在
    • cimport scipy.special.cython_specialでcythonコード上でcimport
  • 型が付いている
In [11]:
%load_ext cython
The cython extension is already loaded. To reload it, use:
  %reload_ext cython
In [12]:
%%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)
Out[12]:
Cython: _cython_magic_684790e4af39b5c80e39e5974eacdfbe.pyx

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.special.sici
    • $\int_0^x \frac{\sin t}{t} dt$と$\gamma + \log (x) + \int_0^x \frac{\cos t - 1}{t} dt$を計算

Tutorial記載の他の内容

  • cythonのprangeを使うことでGILを外して並列に回す例
  • specialに含まれないがnumpy, scipyで容易に実装可能な関数の紹介
    • 2値エントロピー関数、ヘビサイド関数、ステップ関数、ランプ関数

scipy.linalg

  • 線形代数関連の操作を扱うモジュール

numpy.linalgとの違い

  • scipy.linalgnumpy.linalgに存在しない機能を持つ
    • LU分解、シュール分解
    • 様々な方法による擬似逆行列の計算
    • 行列の対数のような行列操作
    • 追加の引数を持つもの(scipy.linalg.eig()は一般化固有値問題を解くために第二の行列を引数に取ることが可能)

numpy.matrixと2次元numpy.ndarray

  • numpy.matrix
    • ndarrayよりも行列操作に便利なインタフェースを持ったクラス
    • 1次元matrixの転置は転置となり、matrixオブジェクト同士の*による乗算は行列の乗算となる
  • numpy.ndarray
    • 1次元ndarrayは.Tによる転置操作は転置にならず、*による乗算は要素同士の乗算となる
  • scipy.linalgの操作はどちらのオブジェクトに対しても同一の操作を実現

行列式の計算

  • linalg.detで計算可能
  • 行列式は余因子展開によって再帰的に求めることが可能
$$ |\textbf{A}|=\sum_j (-1)^{i+j} a_{i,j} M_{i,j} $$
  • $ M_{i,j} = | \textbf{A}_{ij}| $はAの(i,j)小行列式(第i行と第j列を取り除いて得られる行列)$\textbf{A}_{ij}$の行列式
In [13]:
# 行列式の計算
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を用いる
In [14]:
# 逆行列の計算
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
In [15]:
# 特異行列(行列式が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を用いる
  • 下記線形方程式を解くことを考える
$$ \begin{aligned} x + 2y &=& 5 \\ 3x + 4y &=& 6 \\ \end{aligned} $$
  • $\textbf{A}\textbf{x}=\textbf{b}$とみなすと$\textbf{x}=\textbf{A}^{-1}\textbf{b}$と計算しても解けるが、速度の観点と数値的な安定性からlinalg.solveを用いた方が良い
In [16]:
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)で行列やベクトルの様々なノルムを計算可能
  • ベクトルに対してはordパラメータで制御
$$\begin{aligned} || \textbf{x} || = \begin{cases} \max |x_i| & \textrm{ord} = \inf \\ \min |x_i| & \textrm{ord} = -\inf \\ \left( \sum_i |x_i|^{\textrm{ord}} \right)^{1/\textrm{ord}} & |\textrm{ord}| < \infty. \end{cases} \end{aligned} $$
  • 行列に対しては$\pm2, \pm1, \pm\inf$, 'fro'(または'f')が指定可能
$$\begin{aligned}|| \textbf{A} || = \begin{cases} \max_i \sum_j |a_{ij}| & \textrm{ord} = \inf \\ \min_i \sum_j |a_{ij}| & \textrm{ord} = -\inf \\ \max_j \sum_i |a_{ij}| & \textrm{ord} = 1 \\ \min_j \sum_i |a_{ij}| & \textrm{ord} = -1 \\ \max \sigma_i & \textrm{ord} = 2 \\ \min \sigma_i & \textrm{ord} = -2 \\ \sqrt{\textrm{trace} \left(\textbf{A}^H\textbf{A}\right)} = \sqrt{\sum_{i,j} |a_{ij}|^2} & \textrm{ord} = \textrm{'fro'} \end{cases}\end{aligned} $$
In [17]:
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
    • 特異値分解に基づく擬似逆行列

目標変数$y_i$は説明変数$\textbf{x}_i$, パラメータ$c_j$, モデル関数$f_j (\textbf{x}_i )$を用いて下記のように書けるとする

$$ y_i = \sum_j {c_j f_j (\textbf{x}_i)} + \epsilon_i $$

ここで$\epsilon$はデータにおける不確かさを表す。 パラメータ$c_j$を得るために、以下の式で表される二乗誤差をJパラメータ$c_n$で微分して0とおくことで最小化する。

$$ J(\textbf{c}) = \sum_i \left| y_i - \sum_j { c_j f_j (\textbf{x}_i) }\right|^2 $$$$ \frac{\partial J}{\partial c_n^*} = 0 = \sum_i \left( y_i - \sum_j c_j f_j (x_i)\right)(-f_n^*(x_i)) $$

即ち $$ \sum_j c_j \sum_i f_j(x_i)f_n^*(x_i) = \sum_i y_i f_n^* (x_i) $$ $$ \textbf{A}^H\textbf{A}\textbf{c}=\textbf{A}^H\textbf{y}$$ ここで $$\left\{ \textbf{A} \right\}_{ij} = f_j (x_i). $$ $\textbf{A}^H\textbf{A}$に逆行列が存在する場合、 $$ \textbf{c} = \left(\textbf{A}^H\textbf{A}\right)^{-1} \textbf{A}^H\textbf{y} = \textbf{A}^\dagger \textbf{y}$$ $\textbf{A}^\dagger$は擬似逆行列と呼ばれる

疑似逆行列

  • $\textbf{A}$をM行N列の行列とする
  • M > Nのとき
$$ \textbf{A}^\dagger = \left(\textbf{A}^H\textbf{A}\right)^{-1} \textbf{A}^H $$
  • M < Nのとき
$$ \textbf{A}^\# = \textbf{A}^H \left( \textbf{A}\textbf{A}^H\right)^{-1}$$
  • M = Nのとき
$$ \textbf{A}^\dagger = \textbf{A}^\# = \textbf{A}^{-1} $$
  • $x_i$=0.1i for i = 1...10に対し、下記のモデルによって生成される$y_i$を考える。 $$ y_i = 5e^{-x_i} + 4x_i $$ 観測値として下記が得られるとき、 $$ z_i = y_i + \epsilon = c_1 e^{-x_i} + c_2 x_i + \epsilon$$ $z_i$と$[e^{-x_i}, x_i]$を用いてモデルパラメータc1とc2を最小二乗法で推定したい
In [18]:
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.        ]]
In [19]:
# 計画行列と観測値を与えてモデルパラメータを推定
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.linalgeigs, eighs
In [20]:
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]]

特異値分解

  • 擬似逆行列を求めるときや次元削減に使われる
$$ M = U\Sigma V^* $$
  • Mはrank rのm行n列
    • Uはm行m列のユニタリ行列
    • $\Sigma$は対角方向に特異値が並んだm行n列行列
    • $V^*$はn行n列ユニタリ行列Vの随伴行列
  • (m, n) = (m, m) @ (m, n) @ (n, n)
In [21]:
# 正方行列に対する特異値分解
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

  • 高速フーリエ変換、離散コサイン変換などができる
  • scipy.fftpack.fftでFFTを実行
  • scipy.fftpack.fftfreqで周波数軸の値を取得
In [22]:
import IPython.display
In [23]:
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)
Out[23]:
In [24]:
# 窓関数を取得
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
In [25]:
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

  • 曲線の補間が可能
  • scipy.interpolate.interp1dで補間関数を生成
  • 生成した補間関数に補間点を与えると補間値を算出
In [26]:
# 観測値
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)
In [27]:
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()
Out[27]:
<matplotlib.legend.Legend at 0x1d1fb80080>

scipy.signal

  • 信号処理のためのライブラリ
    • フィルタ設計・フィルタリング
    • 畳み込み・時間相関計算
    • 窓関数

ホワイトノイズに対するローパスフィルタの適用

  • numpy.random.randnでホワイトノイズを生成
  • scipy.signal.firwinでローパスフィルタを設計
  • scipy.signal.lfilterでローパスフィルタをホワイトノイズに対して適用
In [28]:
# ホワイトノイズの生成
time_s = 5.0
ys = 0.1 * np.random.randn(int(time_s * fs))
IPython.display.Audio(ys, rate=fs)
Out[28]:
In [29]:
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)
Out[29]:
In [30]:
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]")
Out[30]:
Text(0,0.5,'power [dB]')

参考

音響信号処理についてはQiitaに過去に諸々記事を書いてます

scipy.optimize

  • 各種最適化アルゴリズムが実装
  • scipy.optimize.curve_fitで曲線フィッティング
  • scipy.optimize.minimizeで最小値を見つける
  • scipy.optimize.rootで方程式の根を見つける
In [31]:
# ノイズの乗った正弦波を観測値とする
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]]
In [32]:
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]))
Out[32]:
[<matplotlib.lines.Line2D at 0x1d1fced470>]
In [33]:
# 最小値を見つけたい対象
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))
Out[33]:
[<matplotlib.lines.Line2D at 0x1d1fd7af60>]
In [34]:
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]
In [35]:
# nfev(評価回数)が12に減る
opt.minimize(f, x0=0, method="L-BFGS-B")
Out[35]:
      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])
In [36]:
# 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]
In [37]:
# `scipy.optimize.basinhopping`を用いて大域解を探す
opt.basinhopping(f, 1.0)
Out[37]:
                        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])
In [38]:
# 探索範囲に制約を与えることも可能
res = opt.minimize(f, x0=1,
                   bounds=((0, 10), ))
res.x    
Out[38]:
array([0.])
In [39]:
# 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.]
In [40]:
# もうひとつの根を探す
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]

scipy.stat

  • 統計や確率分布、統計的検定のためのモジュール
  • サンプル点の正規分布へのフィッティング等を紹介
In [41]:
# サンプリングした点のヒストグラム(頻度分布)と確率分布を同時に描画
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()
Out[41]:
<matplotlib.legend.Legend at 0x1d1fc61048>
In [42]:
# 観測値を確率分布にフィッティング
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
Out[42]:
[<matplotlib.lines.Line2D at 0x1d1fad3048>]
In [43]:
# 平均、中央値、四分位数
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

  • 空間データや距離を扱うためのモジュール
  • 距離行列を計算するscipy.spatial.distancecdist, pdistのみ紹介
In [44]:
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
Out[44]:
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.        ]])
In [45]:
# 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.        ]]

scipy.misc

  • deprecatedな関数を除くと4つのみ
    • face
    • central_diff_weight
    • derivative
    • ascent
In [46]:
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
In [47]:
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()
In [48]:
# グレースケールのデモ用画像
ascent = scipy.misc.ascent()
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.imshow(ascent, cmap=cm.gray)
plt.show()
In [49]:
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)
Out[49]:
4.999999999921734
In [50]:
print((f(1.0 + 1e-6) - f(1.0 - 1e-6)) / (2.0 * 1.e-6))
4.999999999921734

デモ

  • scipy.fftpackとscipy.signalを使ったNMFによるモノラル音源分離のデモ
    • NMFはscikit-learnの実装を使用
  • 音楽情報処理ライブラリlibrosaのデモ音源を使用
  • librosa.decompose.decomposeと同等の処理

NMF

  • 非負値行列$\textbf{X}$を$\textbf{X} \simeq \textbf{WH} $と分解する手法
    • Hを基底行列、Wをアクティベーションと呼ぶ
  • scikit-learnのNMFでは下記目的関数を最小化
$$0.5 * ||\textbf{X} - \textbf{WH}||_{Fro}^2 + \alpha l1_{ratio} ||vec(W)||_1 + \alpha l1_{ratio} ||vec(H)||_1 + 0.5 \alpha (1 - l1_{ratio}) ||W||_{Fro}^2 + 0.5 \alpha (1 - l1_{ratio}) ||H||_{Fro}^2 $$
  • 通常最適化には補助関数法が使われる

音響信号における非負行列

  • 音を短時間波形に切り出した各断片をFFT
  • 断片から求めた振幅スペクトルを時間方向に並べたスペクトログラムを行列とみなす
In [51]:
!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)
In [52]:
import librosa
import librosa.display
In [53]:
y, sr = librosa.load(librosa.util.example_audio_file())
IPython.display.Audio(data=y, rate=sr)
Out[53]:
In [54]:
np.pad(np.arange(5), 2, "reflect")
Out[54]:
array([2, 1, 0, 1, 2, 3, 4, 3, 2])
In [55]:
# 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
In [56]:
# 信号の短時間波形への切り出し
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()
In [57]:
# 切り出し間隔を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()
In [58]:
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()
In [59]:
# 更に各短時間波形に窓掛けを行った場合
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()
In [60]:
# 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
In [61]:
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")
Out[61]:
<matplotlib.axes._subplots.AxesSubplot at 0x1d216b7320>
In [62]:
# 得られたスペクトログラム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
In [63]:
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')
Out[63]:
Text(0.5,1,'Components')
In [64]:
# 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
In [65]:
# 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
Out[65]:
In [66]:
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
Out[66]:

まとめ

  • SciPyの各モジュールとコード片を紹介
  • fftpackとsignalで音の信号処理ができる

個人的によく使う関数

  • scipy.ndimage.interpolation.shift
    • 配列を単純にシフトし、要素が無くなったところに指定した値を埋める
    • numpyのrollだと要素が回転しまう
In [67]:
from scipy.ndimage.interpolation import shift
xs = np.array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.]) 
print(shift(xs, -3, cval=0.0))
print(np.roll(xs, -3))
[3. 4. 5. 6. 7. 8. 9. 0. 0. 0.]
[3. 4. 5. 6. 7. 8. 9. 0. 1. 2.]