2019/3/24にmatplotlibのアニメーション機能に関する発表を行いました

2019/3/24にmatplotlibのアニメーションの機能に関する発表を行いました

2019/3/24に「はんなりPython+PyData Osakaの可視化特集会 」のイベントにて標題の発表を行いました。 イベントページはこちらです。 発表資料はいつものように下記のnotebookをRISEでスライドにしたものを使用しました。

In [8]:
from IPython.display import IFrame
IFrame("https://www.hiromasa.info/slide/15.slides.html", "100%", "450px")
Out[8]:

matplotlibのアニメーション機能の紹介

自己紹介

  • 大橋宏正(@wrist)
  • PyData Osaka オーガナイザ
  • メーカー職(音響信号処理の業務に従事)

概要

  • matplotlibのアニメーション機能を紹介
  • 単純な方法の紹介(plt.pause)
  • Animationクラスを用いた方法の紹介
    • ArtistAnimation, FuncAnimation
    • TimedAnimation
  • 動画書き出し

なぜアニメーション描画を行いたいか

  • センサーデータの時系列変化を可視化したい
    • ex. PCM信号
    • c/c++バイナリが吐いたdumpデータの描画
In [1]:
import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt
In [2]:
print(np.__version__)
print(sp.__version__)
print(mpl.__version__)
1.16.2
1.2.1
3.0.3
In [3]:
# osxの場合
%matplotlib osx

osxの場合の制約

  • アニメーション描画
    • backendに"osx"の指定が必要
  • "osx"の指定
    • FrameworkとしてビルドされたPythonが必要
    • 仮想環境(conda含む)のPythonはFrameworkビルドではない

osx(conda)の場合の対処方法

  • 公式FAQでの対処法
    • conda install python.app
    • pythonの代わりにpythonwを実行
    • pythonwを使うkernelを作成すればjupyter notebook上でも実行可能

pythonwを使用するカーネルの作成方法

  • 要ipykernel
$ ipython kernel install --user --name pythonw3 --display-name pythonw3
$ jupyter kernelspec list  # pythonw3カーネルのパスを確認
$ vim ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3/kernel.json # kernel.json記載のpathをpythonwに変更

実際の実行結果

$ jupyter kernelspec list
Available kernels:
  julia-0.4    /Users/wrist/Library/Jupyter/kernels/julia-0.4
  julia_0.3    /Users/wrist/Library/Jupyter/kernels/julia_0.3
  python3      /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/python3

$ ipython kernel install --user --name pythonw3 --display-name pythonw3
Installed kernelspec pythonw3 in /Users/wrist/Library/Jupyter/kernels/pythonw3

$ jupyter kernelspec list
Available kernels:
  julia-0.4    /Users/wrist/Library/Jupyter/kernels/julia-0.4
  julia_0.3    /Users/wrist/Library/Jupyter/kernels/julia_0.3
  pythonw3     /Users/wrist/Library/Jupyter/kernels/pythonw3
  python3      /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/python3

# condaのpython3と同じディレクトリに移動
$ mv ~/Library/Jupyter/kernels/pythonw3 ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/
/Users/wrist/Library/Jupyter/kernels/pythonw3 -> /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3

$ jupyter kernelspec list
Available kernels:
  julia-0.4    /Users/wrist/Library/Jupyter/kernels/julia-0.4
  julia_0.3    /Users/wrist/Library/Jupyter/kernels/julia_0.3
  python3      /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/python3
  pythonw3     /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3

# kernel.jsonに記載されているpathをpythonwに書き直す
$ vim ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3/kernel.json

最も単純なアニメーション方法

  • plt.pause(second)を使って描画を繰り返す
    • plt.plot(),plt.pause(sec), plt.cla()を繰り返す
In [4]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)

fig = plt.figure(1)
ax = fig.add_subplot(111)
In [5]:
# アニメーションしないバージョン
for _ in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx;
plt.show()
In [6]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)

fig = plt.figure(1)
ax = fig.add_subplot(111)
In [7]:
# アニメーションするバージョン
for n in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx;
    plt.pause(0.1); plt.cla();  # 0.1s待ってclear
plt.show()
In [8]:
# アニメーションするバージョン
for n in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx;
    plt.pause(0.1); ax.clear(); # ax.clearでも良い
plt.show()

plotのコストを下げる

  • ax.plotは毎回プロットし直すためコストが大きい
  • linesを取得して毎回set_dataでデータを設定
In [9]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)

fig = plt.figure(1)
ax = fig.add_subplot(111)
In [10]:
lines = ax.plot(xs, np.sin(xs))
for n in range(100):
    lines[0].set_data(xs, np.sin(xs)); xs += dx;
    ax.set_xlim((np.min(xs), np.max(xs)))  # 横軸範囲の再設定
    plt.pause(0.1)  # clear不要
plt.show()

plt.ion(), plt.draw()を使う方法

  • 昔の記事でよく出てくる方法
    • plt.ion() -> plt.show()
    • ax.plot() -> plt.draw() -> plt.cla()を繰り返す
    • インタラクティブモードをonにした上で再描画を都度命令
  • windowsでは過去は確かにこれでいけた(現状は未確認)
  • 手元のosxだとplt.pauseを入れないと動作しない
    • plt.ion()しておく必要は特にない?
In [11]:
print(plt.isinteractive())
True
In [12]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0; xs = np.arange(-np.pi, np.pi, dx)
fig = plt.figure(1); ax = fig.add_subplot(111)
In [13]:
plt.ion(); plt.show();
for n in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx
    plt.draw(); plt.pause(0.1); plt.cla();
In [14]:
print(plt.isinteractive())
True
In [15]:
plt.ioff()
In [16]:
print(plt.isinteractive())
False

matplotlibのアニメーション描画クラス

In [6]:
from IPython.display import Image
Image("./15/animation_class.png")
Out[6]:

ArtistAnimation

  • 描画するプロットのリストを予め作成し与える
  • アニメーション描画クラスの中でも単純なクラス

class matplotlib.animation.ArtistAnimation(fig, artists, *args, **kwargs)

  • artists
    • 各フレームで描画するArtistコレクションのlist
In [17]:
import matplotlib.animation as anim
In [18]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
interval_ms = 100
images = []
In [19]:
fig = plt.figure(1); ax = fig.add_subplot(111)
for _ in range(100):
    images.append(ax.plot(xs , np.sin(xs))); xs += dx;
ani = anim.ArtistAnimation(fig, images, interval=interval_ms, repeat=False)    
plt.show()

問題点

  • 範囲制御ができていない
  • 理由
    • xsの範囲が都度増加しているため
  • 対策
    • xsを固定
    • x軸レンジが固定されてしまう
In [20]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
interval_ms = 100
images = []
In [21]:
fig = plt.figure(1); ax = fig.add_subplot(111)
for n in range(100):
    images.append(ax.plot(xs , np.sin(xs + n*dx)))
ani = anim.ArtistAnimation(fig, images, interval=interval_ms, repeat=False)    
plt.show()

FuncAnimation

  • 描画関数とデータ配列を与えてアニメーション
  • データ配列はgeneratorなどでも良い

class matplotlib.animation.FuncAnimation(fig, func, frames=None, init_func=None, fargs=None, save_count=None, **kwargs)

  • func : callable
    • def func(frame, *fargs) -> iterable_of_artists
    • framesで指定されたデータが渡される
  • init_func : callable, optional
    • def init_func() -> iterable_of_artists
  • frames : iterable, int, generator function, or None, optional
    • iterable: イテレーション可能なデータから描画データを供給
    • integer: 連番range(frames)を供給
    • generator function: def gen_function() -> objを指定
    • None: itertools.countを供給
In [22]:
import matplotlib.animation as anim
In [23]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
ys = np.sin(xs)
interval_ms = 100
In [24]:
# 呼び出しごとにデータを生成するジェネレータ
def gen_data():
    counter = 0
    dx = np.pi / 10.0
    xs = np.arange(-np.pi, np.pi, dx)
    while counter < 100:
        ys = np.sin(xs); xs += dx;
        counter += 1
        yield (xs, ys)
In [25]:
fig = plt.figure(1); ax = fig.add_subplot(111)
gline2d = ax.plot(xs, ys)[0]  # これを上書き

def update_data(data):
    xs, ys = data
    ax.set_xlim((np.min(xs), np.max(xs)))
    gline2d.set_data(xs, ys)
    return gline2d
    
ani = anim.FuncAnimation(fig, update_data, gen_data,
                         interval=interval_ms, repeat=False)
plt.show()
  • update_dataでx軸レンジを都度再設定できる

FuncAnimationを用いたアニメーション描画クラス

  • 先の例はグローバル変数を直接上書きしているため気持ち悪い
  • クラス化することで状態を同時に保持可能
  • generatorや描画関数はインスタンスメソッドにできる
In [26]:
class MyAnimation:
    def __init__(self):
        plt.close("all")
        self.fig = plt.figure(1)
        self.ax = self.fig.add_subplot(111)
        self.dx = np.pi / 10.0
        self.xs = np.arange(-np.pi, np.pi, self.dx)
        self.ys = np.sin(self.xs)
        self.interval_ms = 100
        self.line2d = self.ax.plot(self.xs, self.ys)[0]
        self.counter = 0
    
    def gen(self):
        while self.counter < 100:
            self.ys = np.sin(self.xs)
            self.xs += self.dx;
            self.counter += 1
            yield (self.xs, self.ys)
            
    def update_data(self, data):
        xs, ys = data
        self.ax.set_xlim((np.min(xs), np.max(xs)))
        self.line2d.set_data(xs, ys)
        return self.line2d

    def animate(self):
        ani = anim.FuncAnimation(self.fig, self.update_data, self.gen,
                         interval=self.interval_ms, repeat=False)
        plt.show()
In [27]:
MyAnimation().animate()

TimedAnimation

  • ArtistAnimationFuncAnimationが継承しているクラス
  • 自作クラスで継承し下記メソッドを定義して使用
    • self._draw_frame(self, frame)
    • self.new_frame_seq(self)
    • self._init_draw(self)
In [28]:
class MyTimedAnimation(anim.TimedAnimation):
    def __init__(self):
        plt.close("all")
        self.fig = plt.figure(1)
        self.ax = self.fig.add_subplot(111)
        self.dx = np.pi / 10.0
        self.xs = np.arange(-np.pi, np.pi, self.dx)
        self.ys = np.sin(self.xs)
        self.interval_ms = 100
        self.line2d = self.ax.plot(self.xs, self.ys)[0]
        self.counter = 0

        plt.show()
        anim.TimedAnimation.__init__(self, self.fig, self.interval_ms)
    
    def new_frame_seq(self):
        while self.counter < 100:
            self.ys = np.sin(self.xs)
            self.xs += self.dx;
            self.counter += 1
            yield (self.xs, self.ys)
            
    def _draw_frame(self, data):
        xs, ys = data
        self.ax.set_xlim((np.min(xs), np.max(xs)))
        self.line2d.set_data(xs, ys)
        return self.line2d

    def _init_draw(self):
        pass
In [29]:
MyTimedAnimation()
Out[29]:
<__main__.MyTimedAnimation at 0x11bb5dd30>

動画書き出し

  • saveメソッド
    • アニメーションを動画としてファイル保存可能
    • save(filename, writer="ffmpeg")
      • writerは文字列 or MovieWriterを継承したクラスを指定
  • to_html5_videoメソッド
    • videoタグを生成

save(filename, writer=None, fps=None, dpi=None, codec=None, bitrate=None, extra_args=None, metadata=None, extra_anim=None, savefig_kwargs=None)

In [7]:
Image("./15/animation_inheritance.png")
Out[7]:
In [30]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
interval_ms = 100
images = []
In [31]:
fig = plt.figure(1); ax = fig.add_subplot(111)
for _ in range(100):
    images.append(ax.plot(xs , np.sin(xs))); xs += dx;
ani = anim.ArtistAnimation(fig, images, interval=interval_ms, repeat=False)    
ani.save("./test.gif", "imagemagick")
In [32]:
!ls
animation_class.webp
animation_inheritance.webp
pydata_osaka_20190323_matplotlib_animation.ipynb
test.gif
In [4]:
Image("./15/test.gif")
Out[4]:
<IPython.core.display.Image object>
In [33]:
# 要ffmpeg
from IPython.display import HTML
HTML(ani.to_html5_video())
Out[33]:

まとめ

  • matplotlibのアニメーション機能を紹介
    • 単純な方法
    • ArtistAnimation, FuncAnimation, TimedAnimation
    • 動画書き出しの方法

2018/12/21にPyData Osakaの梅キャンPython勉強会コラボ回でpybind11について発表しました

2018/12/21(金)にPyData Osakaの梅キャンPython勉強会コラボ回でpybind11について発表しました

ブログ記事にするのが遅くなりましたが2018/12/21に【大阪工業大学】特別回 梅キャンPython勉強会【梅田キャンパス】 でpybind11に関する発表を行いました。 下記はその際の資料になります。

In [1]:
from IPython.display import IFrame
IFrame("https://www.hiromasa.info/slide/14.slides.html", "100%", "450px")
Out[1]:

pybind11の紹介

大橋宏正(PyDataオーガナイザー)

自己紹介

  • 大橋 宏正(@wrist)
  • PyDataオーガナイザー
  • 某メーカー勤務
    • 音響信号処理屋

pybind11とは

他のライブラリとの比較

  • Python C API
    • 自分で参照カウント制御が必要
    • 学習障壁が高い
  • ctypes
    • dllがあればビルド不要で使える
    • 返り値や引数の型を厳密に設定する必要有
  • cython
    • pyxの独自記法を覚える必要
    • 一旦pyxファイルをcにコンパイル
      • 更にcファイルをコンパイル
  • Boost.python
    • Boostが必要
    • 便利だが巨大でありコンパイルに時間がかかる

インストール

In [1]:
!pip install pybind11
Requirement already satisfied: pybind11 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (2.2.4)

ipybind

ipybindのインストール

  • リポジトリをクローン
  • python setup.py build
  • python setup.py install

使ってみる

  • ipybind経由でpybind11を使用
  • 最初にload_extでipybindを有効化
In [2]:
%load_ext ipybind

pybind11の基本

  • PYBIND_MODULE(モジュール名, m)
  • m.def(公開名, 関数へのポインタ, 説明)
In [3]:
%%pybind11

int add(int x, int y){
    return x + y;
}

PYBIND11_MODULE(myadd, m){
    m.def("c_add", &add, "Add two integers.");
}
In [4]:
c_add(1, 2)
Out[4]:
3
In [5]:
help(c_add)
Help on built-in function c_add in module pybind11_3e2e466:

c_add(...) method of builtins.PyCapsule instance
    c_add(arg0: int, arg1: int) -> int
    
    Add two integers.

In [6]:
# c_add("foo", "bar")
In [7]:
%%pybind11

PYBIND11_MODULE(myadd, m){
    m.def("lambda_add", [](int a, int b){ return a + b;},
         "Add two integers.");
}
In [8]:
lambda_add(1, 2)
Out[8]:
3

その他のコンパイル方法

  • 手動ビルド
  • distutils
    • setup.pyを書く
  • cmake
    • pybind11_add_module(myadd myadd.cpp)

単純なクラスをラップ

  • c++のクラスを定義
  • pybind11でラップ
    • py::classs_<C++クラス名>(m, 公開クラス名)
In [9]:
%%pybind11

namespace py = pybind11;

class Coordinate {
public:
    int x_;
    int y_;
    Coordinate(){ x_ = 0; y_ = 0; }
    Coordinate(int x, int y){ x_ = x; y_ = y; }
    float Norm(){ return sqrt(x_*x_ + y_*y_); }
};

PYBIND11_MODULE(coord, m) {
    py::class_<Coordinate>(m, "Coordinate")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate::x_)
        .def_readonly("y", &Coordinate::y_)
        .def_property_readonly("norm", &Coordinate::Norm);
}

コンストラクタの定義

  • ... .def(py::init<オーバーロード引数>())
Coordinate(){ x_ = 0; y_ = 0; }
    Coordinate(int x, int y){ x_ = x; y_ = y; }

は下記のように記載

py::class_<Coordinate>(m, "Coordinate")
        .def(py::init<>())
        .def(py::init<int, int>())

プロパティの定義

  • ... .def_readonly(公開変数名, メンバへのポインタ)
py::class_<Coordinate>(m, "Coordinate")
        .def_readonly("x", &Coordinate::x_)
        .def_readonly("y", &Coordinate::y_)

computed propetyの定義

  • ... .def_propety_readonly(公開変数名, 計算関数へのポインタ)
float Norm(){ return sqrt(x_*x_ + y_*y_); }

.def_property_readonly("norm", &Coordinate::Norm);
In [10]:
point = Coordinate(3, 4)
point
Out[10]:
<pybind11_c19ca2a.Coordinate at 0x10ba96df8>
In [11]:
print(point.x, point.y)
3 4
In [12]:
point.norm
Out[12]:
5.0

さらなる例

  • 演算子オーバーロード、python特殊メソッド
  • docstring、引数名、キーワード引数
  • *arg, **kwarg

演算子オーバーロード

  • 等値演算子
    • cppだとbool operator==(const Cls& rhs)
    • pythonだとCls.__eq__(self, other)
    • これらを対応付ける
In [13]:
int.__eq__(1, 1)
Out[13]:
True
In [14]:
%%pybind11

namespace py = pybind11;

class Coordinate2 {
public:
    int x_; int y_;
    Coordinate2(){ x_ = 0; y_ = 0; }
    Coordinate2(int x, int y){ x_ = x; y_ = y; }
    bool operator==(const Coordinate2& rhs){
        return (x_ == rhs.x_) && (y_ == rhs.y_);
    }
};

PYBIND11_MODULE(coord2, m) {
    py::class_<Coordinate2>(m, "Coordinate2")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate2::x_)
        .def_readonly("y", &Coordinate2::y_)
        .def("__eq__", &Coordinate2::operator==);
}
In [15]:
point1 = Coordinate2(3, 4)
point2 = Coordinate2(3, 4)
point3 = Coordinate2(1, 2)
In [16]:
point1 == point2
Out[16]:
True
In [17]:
point1 == point3
Out[17]:
False

演算子の代替記法

.def("__eq__", &Coordinate2::operator==);

はc++のクラスでoperator==が定義されていれば下記のようにも書ける

.def("__eq__", 
             [](const Coordinate2& self,
                const Coordinate2& other){
                 return self == other; 
             });

Python特殊メソッド

  • __repr__の定義を考える
  • 先の__eq__と同様に第一引数にselfを取るメソッドを定義
In [18]:
%%pybind11

namespace py = pybind11;

class Coordinate3 {
public:
    int x_; int y_;
    Coordinate3(){ x_ = 0; y_ = 0; }
    Coordinate3(int x, int y){ x_ = x; y_ = y; }
};

PYBIND11_MODULE(coord3, m) {
    py::class_<Coordinate3>(m, "Coordinate3")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate3::x_)
        .def_readonly("y", &Coordinate3::y_)
        .def("__repr__" , [](const Coordinate3& self){
            return std::string() + "(x, y) = (" + std::to_string(self.x_) + ", " + std::to_string(self.y_) + ")";
        });
} 
.def("__repr__", [](const Coordinate3& self){
    return std::string() + 
        "(x, y) = (" +
        std::to_string(self.x_) +
        ", " +
        std::to_string(self.y_) + ")";
});
In [19]:
point = Coordinate3(1, 2)
point
Out[19]:
(x, y) = (1, 2)

docstirng、引数名、キーワード引数

  • docstring
    • m.defの第三引数に文字列を渡す
  • 引数名
    • 第四引数以降にpy::arg("...")を渡す
  • キーワード引数
    • py::arg("...") = デフォルト値とすれば良い
In [20]:
%%pybind11

namespace py = pybind11;

class Coordinate4 {
public:
    int x_; int y_;
    Coordinate4(){ x_ = 0; y_ = 0; }
    Coordinate4(int x, int y){ x_ = x; y_ = y; }
    float DistanceTo(const Coordinate4& target){ 
        int dx = target.x_ - x_, dy = target.y_ - y_;
        return sqrt( dx*dx + dy*dy);
    }
};

PYBIND11_MODULE(coord4, m) {
    py::class_<Coordinate4>(m, "Coordinate4")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate4::x_)
        .def_readonly("y", &Coordinate4::y_)
        .def("distance_to", &Coordinate4::DistanceTo,
             "distance to target",
             py::arg("target"));
}
float DistanceTo(const Coordinate4& target){ 
        int dx = target.x_ - x_, dy = target.y_ - y_;
        return sqrt( dx*dx + dy*dy);
    }

を下記のように設定

.def("distance_to", &Coordinate4::DistanceTo,
             "distance to target",
             py::arg("target"));
In [24]:
point = Coordinate4(1, 2)
target = Coordinate4(3, 4)
point.distance_to(target)
Out[24]:
2.8284270763397217
In [40]:
point.distance_to?

*args, **kwarg

  • py::args (py::tupleのサブクラス)
  • py::kwargs (py::dictのサブクラス)
m.def("count_args", [](py::args a, py::kwargs kw) {
    py::print(a.size(), "args,", kw.size(), "kwargs");
});

その他の話題

  • C++とPython間のオブジェクトの取扱
  • 関数とコールバック
  • Numpy対応

C++とPython間のオブジェクトの取扱

  • C++とPythonのオブジェクトを混在してラッパーを定義可能
  • 3種類の方法
    • Python側でC++の型のインスタンスを引数に取る
    • C++側でPythonの型のインスタンスを引数に取る
    • C++の型とPythonの型を型変換する

Python側でC++の型のインスタンスを引数に取る

py::class_<Foo>(m, "Foo");
m.def("f1", [](const Foo& foo){ ...}

C++側でPythonの型のインスタンスを引数に取る

m.def("f2", [](py::list list){ ... }

C++の型とPythonの型を型変換する

m.def("f3", [](int x) { ... });
m.def("f4", [](const std::string& s) { ... });
m.def("f5", [](const std::vector<int>& v) { ... });

サポートされている型変換

  • スカラー値
    • integer types, float, double, bool, char
  • 文字列
    • std::string, const char *
  • タプル
    • std::pair<F, S>, std::tuple<...>
  • シーケンス
    • std::vector, std::list, std::array<T, n>
  • マップ
    • std::map<K, V>, std::unordered_map<K, V>
  • 集合
    • std::set, std::unordered_set
  • 関数
    • std::function<...>
  • Date/time
    • std::chrono::duration, std::chrono::time_point
  • Optional
    • std::optional, std::experimental::optional

関数とコールバック

  • std::functionでpythonの関数を受け取れる
In [38]:
%%pybind11

#include <pybind11/functional.h>

PYBIND11_MODULE(callback_test, m) {
    m.def("for_even",
          [](int n, std::function<void(int)> f) {
              for (int i = 0; i < n; ++i){ if (i % 2 == 0) f(i); }
          }
         );
}
In [39]:
def py_callback_func(x):
    print("called {0}".format(x))
    
for_even(10, py_callback_func)
called 0
called 2
called 4
called 6
called 8

numpy対応

  • #include <pybind11/numpy.h>
  • py::array_t<type>
    • ndarrayを受け取ることが可能
  • arr.ndim, arr.shape(n), arr.size()
  • バッファへのポインタ
    • arr.data(), arr.mutable_data()
    • インデックス指定も可能
  • 要素への直接アクセス
    • arr.unchecked(), arr.mutable_unchecked()
    • ref(i, j, k)としてアクセス可能
In [31]:
%%pybind11

#include <pybind11/numpy.h>

namespace py = pybind11;

PYBIND11_MODULE(numpy_bind, m) {
    m.def("add_2d",
          [](py::array_t<double> x, py::array_t<double> y){
              auto r = x.mutable_unchecked<2>();
              auto c = y.unchecked<2>();

            for(ssize_t i = 0; i < r.shape(0); i++){
                for(ssize_t j = 0; j < r.shape(1); j++){
                    r(i, j) += c(i, j);
                }
            }
          },  py::arg().noconvert(), py::arg().noconvert());
}
In [33]:
import numpy as np
In [34]:
x = np.array([[1.0, 2.0], [3.0, 4.0]])
y = np.array([[2.0, 3.0], [4.0, 5.0]])
add_2d(x, y)
x
Out[34]:
array([[3., 5.],
       [7., 9.]])

まとめ

  • Python向けのC拡張モジュールをC++11で作成するためのライブラリpybind11について紹介
  • ipybindを用いたJupyter Notebookからの実行方法の紹介
  • numpy連携を紹介

Pyroomacousticsの紹介

Pyroomacousticsの紹介

この記事はQiita Python Advent Calendar 2018の9日目の記事です。

Pyroomacousticsとは

Pyroomacousticsは音響アレイ信号処理に関するアルゴリズムの迅速な開発や評価を目的としたライブラリです。ライブラリは主に下記の3つの要素から成り立っています。

  • 2次元または3次元空間に配置された複数の音源と複数のマイクロホンに関するシミュレーション環境を迅速に構築するためのオブジェクト指向インタフェースの提供
  • 一般的な多面体形状の室内に対して室内インパルス応答を効率的に生成するための鏡像法の高速なC実装と、音源と受音点の間の伝搬のシミュレーション環境の提供
  • STFT, ビームフォーミング、方向推定(音源定位)、適応フィルタ、音源分離に対する一般的なアルゴリズムのリファレンス実装の提供

ここではデモ用notebookを参考に、使用方法や何ができるかについて説明していきます。 なおドキュメントはhttps://pyroomacoustics.readthedocs.io/en/pypi-release/index.html にあります。

インストール方法

pipで下記のようにインストールできます。また別途pipでcookiecutterをインストールしておくと、後述のようにpyroomacousitcs用のシミュレーションスクリプトを対話的に生成することが可能であるため便利です。

In [1]:
!pip install -U pyroomacoustics
Requirement already up-to-date: pyroomacoustics in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.1.20)
Requirement already satisfied, skipping upgrade: numpy in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from pyroomacoustics) (1.14.2)
Requirement already satisfied, skipping upgrade: scipy>=0.18.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from pyroomacoustics) (1.1.0)
Requirement already satisfied, skipping upgrade: Cython in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from pyroomacoustics) (0.29)
In [2]:
!pip install cookiecutter
Requirement already satisfied: cookiecutter in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (1.6.0)
Requirement already satisfied: jinja2>=2.7 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (2.10)
Requirement already satisfied: poyo>=0.1.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.4.2)
Requirement already satisfied: jinja2-time>=0.1.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.2.0)
Requirement already satisfied: requests>=2.18.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (2.19.1)
Requirement already satisfied: future>=0.15.2 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.16.0)
Requirement already satisfied: click>=5.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (7.0)
Requirement already satisfied: whichcraft>=0.4.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.5.2)
Requirement already satisfied: binaryornot>=0.2.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.4.4)
Requirement already satisfied: MarkupSafe>=0.23 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from jinja2>=2.7->cookiecutter) (1.0)
Requirement already satisfied: arrow in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from jinja2-time>=0.1.0->cookiecutter) (0.12.1)
Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (3.0.4)
Requirement already satisfied: idna<2.8,>=2.5 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (2.7)
Requirement already satisfied: certifi>=2017.4.17 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (2018.10.15)
Requirement already satisfied: urllib3<1.24,>=1.21.1 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (1.22)
Requirement already satisfied: python-dateutil in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from arrow->jinja2-time>=0.1.0->cookiecutter) (2.7.5)
Requirement already satisfied: six>=1.5 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from python-dateutil->arrow->jinja2-time>=0.1.0->cookiecutter) (1.11.0)

Pyroomacousitcsの使用方法

まず下記ライブラリをインポートしておきます。

In [3]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.signal import fftconvolve
import IPython
import pyroomacoustics as pra

roomオブジェクトの生成

下記のように2次元座標を要素に持つndarrayを作成し、pra.Room.from_cornersに与えることで2次元または3次元の室内空間を表すオブジェクトを生成することが可能です。またroomオブジェクトのplotメソッドで空間を描画することが可能です。座標の単位はmです。

In [4]:
corners = np.array([[0,0], [0,3], [5,3], [5,1], [3,1], [3,0]]).T  # [x,y]
room = pra.Room.from_corners(corners)

fig, ax = room.plot()
ax.set_xlim([-1, 6])
ax.set_ylim([-1, 4]);

3次元空間については2次元空間をextrudeメソッドで指定した高さまで持ち上げます。

In [5]:
room = pra.Room.from_corners(corners)
room.extrude(2.)

fig, ax = room.plot()
ax.set_xlim([0, 5])
ax.set_ylim([0, 3])
ax.set_zlim([0, 2]);

なお、単純な正六面体の部屋であればShoeBoxメソッドを使ったほうが早いです。

In [6]:
room = pra.ShoeBox([2.0, 2.5])
fig, ax = room.plot()
ax.set_xlim([0, 3])
ax.set_ylim([0, 3])
Out[6]:
(0, 3)
In [7]:
room = pra.ShoeBox([2.0, 2.5, 2.0])
fig, ax = room.plot()
ax.set_xlim([0, 3])
ax.set_ylim([0, 3])
ax.set_zlim([0, 3]);

音源とマイクロホンの追加

roomオブジェクトに対してadd_sourceメソッドを用いることで指定した座標に音源を追加することが可能です。音源はndarrayである必要があるので、ここではscipy.ioのwavfileを用いて読み込んだ信号を用います

音源はデモ用ノートブックではCMU Arctic Databaseから取得した音声ファイルを用いています。 ここからダウンロードしたtar.gzを解凍したその中のwavディレクトリに入っています。後で気づきましたが、pyroomacoustics自体がデータセットのラッパーも提供しています。

In [8]:
# specify signal source
fs, signal = wavfile.read("work/arctic_a0010.wav")

# add source to 2D room
room = pra.Room.from_corners(corners, fs=fs)
room.add_source([1.,1.], signal=signal)

fig, ax = room.plot()

黒点が音源位置を表しています。追加した音源はこんな感じです。

In [9]:
IPython.display.Audio(signal, rate=fs)
Out[9]:

音源を追加した時点でroomオブジェクトにもfsが設定される模様です。試したところデフォルトでは8kHzが指定されているようでした。

In [10]:
room.fs
Out[10]:
16000

この空間に更に環状のマイクロホンアレイを追加します。 ここではpra.circular_2D_arrayを用いて6素子からなる中心座標[2.0, 2.0]、半径0.1の環状アレイの各マイクロホンの座標を計算します。 各マイクロホンアレイの座標をfsとともにpra.MicrophoneArrayに渡すことでマイクロホンアレイオブジェクトを生成し、これをroomオブジェクトのadd_microphone_arrayメソッドにわたすことでマイクロホンアレイを室内空間へと追加します。

In [11]:
R = pra.circular_2D_array(center=[2.,2.], M=6, phi0=0, radius=0.1)
room.add_microphone_array(pra.MicrophoneArray(R, room.fs))

fig, ax = room.plot()

xがマイクロホンを表します。

鏡像法を用いた室内インパルス応答の生成

音源とマイクロホンアレイが追加された室内では、音源が各マイクロホンへとどのように伝搬するかを表すインパルス応答を生成することができます。 ここではこのインパルス応答を鏡像法に基づいて生成します。 鏡像法とはざっくり言うと室内の壁面での反射波は壁面に対して音源と対象な位置に鏡像が存在すると考え、音源位置から放射された信号が壁面に伝達する度に鏡像を計算して、最終的には実音源とこれらの鏡像からマイクロホンへと到来する波の重ね合わせがマイクロホンにおける観測信号となるシミュレーション手法です。反射回数を何回まで考慮するかや、壁面での吸音率、取り扱うサンプリングレートなどによって鏡像法で推定されるインパルス応答の精度が変わってくることが知られています。

ここでは3次元の室内空間を生成し、[1.0, 1.0, 0.5]の位置に音源を追加、[3.5, 2.0, 0.5]と[3.5, 2.0, 0.5]の位置にマイクロホンアレイを追加します。 更に、roomオブジェクトの生成時に最大何次の反射までを考慮するかを表すmax_orderと壁面での反射時の吸音率を表すabsorptionが指定されています。max_orderが8であるため8次の反射までが考慮され、吸音率が0.2であるため反射率が0.8、すなわち壁面での反射時には音波が0.8倍されることを意味しています。

この状態でroomオブジェクトに対してimage_source_modelメソッドを呼ぶと、鏡像法に基づき室内の反射を模擬することが可能です。use_libroomのオプションは鏡像法のC実装を用いるかどうかを制御するオプションのようです。

なお、以下では3次反射までを描画しています。

In [12]:
# specify signal source
fs, signal = wavfile.read("work/arctic_a0010.wav")

# set max_order to a low value for a quick (but less accurate) RIR
room = pra.Room.from_corners(corners, fs=fs, max_order=8, absorption=0.2)
room.extrude(2.)

# add source and set the signal to WAV file content
room.add_source([1., 1., 0.5], signal=signal)

# add two-microphone array
R = np.array([[3.5, 3.6], [2., 2.], [0.5,  0.5]])  # [[x], [y], [z]]
room.add_microphone_array(pra.MicrophoneArray(R, room.fs))

# compute image sources
room.image_source_model(use_libroom=True)

# visualize 3D polyhedron room and image sources
fig, ax = room.plot(img_order=3)
fig.set_size_inches(18.5, 10.5)

一度鏡像法のシミュレーションを実行すると、下記のように室内インパルス応答を描画することが可能になります。

In [13]:
room.plot_rir()
fig = plt.gcf()
fig.set_size_inches(20, 10)

なおこの室内インパルス応答(RIR)はroomオブジェクトの下にrirとして格納されています。

In [14]:
room.rir
Out[14]:
[[array([0., 0., 0., ..., 0., 0., 0.])],
 [array([0., 0., 0., ..., 0., 0., 0.])]]

pyroomacousticsではこの室内インパルス応答を音源に対して畳み込んだ場合のシミュレーションを行い、マイクロホンでの観測音源を生成することができます。 roomオブジェクトに対してsimulateメソッドを呼ぶことで、室内インパルス応答を畳み込んだ観測信号がroom.mic_array.signalsとして生成されます。ここではマイクロホンが2つ存在しているため、2つの観測信号が生成されています。

In [15]:
room.simulate()
print(room.mic_array.signals.shape)
(2, 55372)

実際に聞いてみて元音源と比較することが可能です。

In [16]:
# original signal
print("Original WAV:")
IPython.display.Audio(signal, rate=fs)
Original WAV:
Out[16]:
In [17]:
print("Simulated propagation to first mic:")
IPython.display.Audio(room.mic_array.signals[0,:], rate=fs)
Simulated propagation to first mic:
Out[17]:

cookiecutterを用いたシミュレーションスクリプトの雛形作成

冒頭でpipでインストールしたcookiecutterを使うとここまでの内容を実行するためのシミュレーションスクリプトの雛形を下記のように対話的に作成することができます。 ここで$はシェルのプロンプトを表します。全質問に対して何も入力していないためデフォルトの項目([]内の値)が選択されています。

$ cookiecutter gh:fakufaku/cookiecutter-pyroomacoustics-sim
You've downloaded /Users/wrist/.cookiecutters/cookiecutter-pyroomacoustics-sim before.
Is it okay to delete and re-download it? [yes]:
script_name [my_room_simulation]:
samplerate [16000]:
absorption [0.25]:
max_order [17]:
Select room_type:
1 - shoebox
2 - polygonal
Choose from 1, 2 (1, 2) [1]:
Select dimension:
1 - 3D
2 - 2D
Choose from 1, 2 (1, 2) [1]:

これによって生成されたスクリプトの中身は下記のようになっています。

In [18]:
!cat work/my_room_simulation.py | pygmentize
import numpy as np
import pyroomacoustics as pra
import matplotlib.pyplot as plt

# Simulation parameters
fs = 16000
absorption = 0.25
max_order = 17

# Geometry of the room and location of sources and microphones
room_dim = np.array([10, 7.5, 3])
source_loc = np.array([2.51, 3.57, 1.7])
mic_loc = np.c_[[7, 6.1, 1.3],[6.9, 6.1, 1.3]]

# Create the room itself
room = pra.ShoeBox(room_dim, fs=fs, absorption=absorption, max_order=max_order)

# Place a source of white noise playing for 5 s
source_signal = np.random.randn(fs * 5)
room.add_source(source_loc, signal=source_signal)

# Place the microphone array
room.add_microphone_array(
        pra.MicrophoneArray(mic_loc, fs=room.fs)
        )

# Now the setup is finished, run the simulation
room.simulate()

# As an example, we plot the simulated signals, the RIRs, and the room and a
# few images sources

# The microphone signal are in the rows of `room.mic_array.signals`
mic_signals = room.mic_array.signals
plt.figure()
plt.subplot(1,2,1)
plt.plot(np.arange(mic_signals.shape[1]) / fs, mic_signals[0])
plt.title('Microphone 0 signal')
plt.xlabel('Time [s]')
plt.subplot(1,2,2)
plt.plot(np.arange(mic_signals.shape[1]) / fs, mic_signals[1])
plt.title('Microphone 1 signal')
plt.xlabel('Time [s]')
plt.tight_layout()

# Plot the room and the image sources
room.plot(img_order=4)
plt.title('The room with 4 generations of image sources')

# Plot the impulse responses
plt.figure()
room.plot_rir()

plt.show()

音源としてnp.random.randnによるホワイトノイズが指定されていることが分かります。 これを実行すると下記のように観測信号、室内環境、室内インパルス応答がプロットされます。 分析の取っ掛かりに便利そうですね。

In [19]:
%run work/my_room_simulation.py

ビームフォーミング

Pyroomacousticsはビームフォーミングのためのリファレンス実装を提供しています。 時間領域および周波数領域のどちらにおいても固定係数のビームフォーミングが使用可能です。 Pyroomacousticsが提供しているのは古典的な遅延和法(Delay and Sum;DAS)と最小分散法(MVDR)に基づくビームフォーマです。

下記では遅延和法を使用しています。 使用しているノイズはGoogleのSpeech Command datasetに含まれるものであり、datasetを解凍したディレクトリの_background_noise_ディレクトリ以下に入ってました。読み込みの際にwarningが出ますがコメントにある通り問題ないようですね。 音声信号を[1, 4.5]、妨害音(ノイズ)を[3.5, 3]の位置に配置し、マイクロホンアレイを[2, 1.5]位置を中心とした半径0.0375[m]=3.75[cm]の円周上および中心に配置しています。これはAmazon Echo相当の配置だとコメントに書いてありますね。ここではマイクの座標情報echopra.MicrophoneArrayではなくpra.Beamformerにfs, fft長、秒数で表されたフィルタサイズとともに渡しています。 このマイクオブジェクトをルームオブジェクトroom_bfadd_microphone_arrayした上で、マイクロホンアレイオブジェクトmicsrake_delay_and_sum_weightsを音源オブジェクトroom_bf.sources[0][:1]を与えて呼ぶことで遅延和法の固定係数を計算しているようです。

In [20]:
Lg_t = 0.100                # filter size in seconds
Lg = np.ceil(Lg_t*fs)       # in samples

# specify signal and noise source
fs, signal = wavfile.read("work/arctic_a0010.wav")
fs, noise = wavfile.read("work/exercise_bike.wav")  # may spit out a warning when reading but it's alright!

# Create 4x6 shoebox room with source and interferer and simulate
room_bf = pra.ShoeBox([4,6], fs=fs, max_order=12)
source = np.array([1, 4.5])
interferer = np.array([3.5, 3.])
room_bf.add_source(source, delay=0., signal=signal)
room_bf.add_source(interferer, delay=0., signal=noise[:len(signal)])

# Create geometry equivalent to Amazon Echo
center = [2, 1.5]; radius = 37.5e-3
fft_len = 512
echo = pra.circular_2D_array(center=center, M=6, phi0=0, radius=radius)
echo = np.concatenate((echo, np.array(center, ndmin=2).T), axis=1)
mics = pra.Beamformer(echo, room_bf.fs, N=fft_len, Lg=Lg)
room_bf.add_microphone_array(mics)

# Compute DAS weights
mics.rake_delay_and_sum_weights(room_bf.sources[0][:1])

# plot the room and resulting beamformer before simulation
fig, ax = room_bf.plot(freq=[500, 1000, 2000, 4000], img_order=0)
ax.legend(['500', '1000', '2000', '4000'])
fig.set_size_inches(20, 8)
/Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages/scipy/io/wavfile.py:273: WavFileWarning: Chunk (non-data) not understood, skipping it.
  WavFileWarning)

room_bf.plotの引数にfreqを与えると周波数に応じたpolar patternも表示してくれるようです。 なお遅延和法で推定された固定係数はmics.weightsに格納されているようです。

In [21]:
mics.weights.shape
Out[21]:
(7, 257)

shapeが(7, 257)なので、マイクロホンアレイの7本のマイクそれぞれに対して257ptの固定係数、即ちFFT長/2+1個の周波数領域の係数が入っているようです。一般に遅延和法は音源方向から各マイクロホンへの到来時間が同一になるように遅延を補正して重畳させることで強調する手法なので、この固定係数は一定ゲインで位相差のみを制御するようなものになっていることが想像されます。

In [22]:
mics.weights
Out[22]:
array([[ 3.58128812e-03+0.j        , -8.73334814e-04-0.00347317j,
        -3.15534434e-03+0.00169394j, ..., -3.46061631e-03-0.00092182j,
        -5.00866112e-05+0.00358094j,  3.48504461e-03-0.00082467j],
       [ 3.62330738e-03+0.j        , -8.09341157e-04-0.00353176j,
        -3.26174098e-03+0.00157778j, ..., -2.83947598e-03+0.00225072j,
         2.82810664e-03+0.00226499j,  1.57604295e-03-0.00326258j],
       [ 3.63718941e-03+0.j        , -7.88120007e-04-0.00355078j,
        -3.29564374e-03+0.00153879j, ..., -1.74332143e-03-0.00319217j,
        -2.73858467e-03+0.0023936j ,  2.93013528e-03+0.00215487j],
       ...,
       [ 3.56674207e-03+0.j        , -8.95397732e-04-0.00345252j,
        -3.11717936e-03+0.00173345j, ...,  1.82547501e-04+0.00356207j,
         3.40217087e-03-0.00107093j, -1.89071480e-03-0.00302438j],
       [ 3.55364827e-03+0.j        , -9.15215577e-04-0.00343377j,
        -3.08223439e-03+0.00176869j, ...,  3.49104208e-03-0.00066411j,
        -1.54079876e-03-0.00320224j, -2.69739964e-03+0.00231354j],
       [ 3.59494372e-03+0.j        , -8.52579157e-04-0.00349238j,
        -3.19054728e-03+0.00165651j, ...,  1.54024356e-03-0.00324827j,
        -3.52088528e-03-0.00072594j,  1.29787261e-04+0.0035926j ]])
In [23]:
np.abs(mics.weights)
Out[23]:
array([[0.00358129, 0.00358129, 0.00358129, ..., 0.00358129, 0.00358129,
        0.00358129],
       [0.00362331, 0.00362331, 0.00362331, ..., 0.00362331, 0.00362331,
        0.00362331],
       [0.00363719, 0.00363719, 0.00363719, ..., 0.00363719, 0.00363719,
        0.00363719],
       ...,
       [0.00356674, 0.00356674, 0.00356674, ..., 0.00356674, 0.00356674,
        0.00356674],
       [0.00355365, 0.00355365, 0.00355365, ..., 0.00355365, 0.00355365,
        0.00355365],
       [0.00359494, 0.00359494, 0.00359494, ..., 0.00359494, 0.00359494,
        0.00359494]])
In [24]:
np.angle(mics.weights)
Out[24]:
array([[ 0.        , -1.8171409 ,  2.6489035 , ..., -2.88126198,
         1.58478242, -0.23235848],
       [ 0.        , -1.79606764,  2.69105004, ...,  2.47134787,
         0.67528024, -1.1207874 ],
       [ 0.        , -1.7892126 ,  2.70476011, ..., -2.07065847,
         2.42331424,  0.63410164],
       ...,
       [ 0.        , -1.82455165,  2.63408201, ...,  1.5195935 ,
        -0.30495815, -2.1295098 ],
       [ 0.        , -1.83127441,  2.6206365 , ..., -0.18798625,
        -2.01926066,  2.43265024],
       [ 0.        , -1.81023839,  2.66270853, ..., -1.12802278,
        -2.93826117,  1.53468576]])

上記のように各周波数ビンにおける振幅値は一定ですが、位相についてはビンごとに異なっていることが分かりました。

続いてこの室内において鏡像法により推定された室内インパルス応答を確認した上で、観測信号が遅延和法によってどのように変化したかを聞いてみます。

In [25]:
# RIRを計算
room_bf.compute_rir()
# 畳み込みを実行し観測信号を生成
room_bf.simulate()
In [26]:
# room_bf.plot_rir()が見にくいので自分で描画
fig = plt.figure(1, figsize=(12,8))

for i, rir_s in enumerate(room_bf.rir):
    for j, rir in enumerate(rir_s):
        ax = fig.add_subplot(7, 2, 2*i + j + 1)
        ax.plot(rir,label="source {0} to mic {1}".format(j, i))
        ax.legend(loc=0)
In [27]:
print("Center Mic:")
IPython.display.Audio(room_bf.mic_array.signals[-1,:], rate=fs)  # -1はcenter mic
Center Mic:
Out[27]:
In [28]:
signal_das = mics.process(FD=False)
print("DAS Beamformed Signal:")
IPython.display.Audio(signal_das, rate=fs)
DAS Beamformed Signal:
Out[28]:

遅延和法による出力ではより鋭い指向性を持つ高域のノイズが顕著に抑圧されていることを確認できます。

その他

時間がなくなってしまったためここでは紹介しきれませんでしたが、pyroomacousticsには他にも音源定位、適応フィルタ、音源分離アルゴリズムのリファレンス実装が用意されています。これらについては改めて今後別の記事で紹介したいと思います。 特に音源分離アルゴリズムとして補助関数法を用いたIVAやILRMAといった最近のアルゴリズムが実装されている点が非常に熱いです。

まとめ

室内形状から室内インパルス応答を算出しシミュレーションが可能なpyroomacousticsを紹介し、シミュレーションの実行方法とビームフォーミングの紹介を行いました。実際にインパルス応答を収録することなくシミュレーションを手軽に行うことができるため個人的にも非常に便利に使用できそうだと感じました。