【Python/PyRoomAcoustics】Pythonで室内音響シミュレーション

# PyRoomAcousticsとは

室内音響のほかビームフォーミングから、音源方向推定、音源分離などあらゆる音響信号処理をサポートする音響系最強Pythonモジュールです。

[【Python/PyRoomAcoustics】ILRMAでブラインド音源分離](https://qiita.com/oozzZZZZ/items/496c44fcab879871a247)

PyRoomAcousticsを使うと簡単に室内音響シミュレーションができます。
シミュレーションはImage Source(鏡像法)で実行され、Image Source(鏡像法)とRay Tracing(光線追跡法)によるハイブリッドなシミュレータも実験的に試せるようです。

1. 部屋の設定
2. マイクの配置
3. 音源の配置
4. シミュレーションの実行

の４ステップで簡単に室内音響シミュレーションが試せます。

今回はJupyter Notebookでのデモです。

```module.py
import pyroomacoustics as pra
import matplotlib.pyplot as plt
from IPython.display import display, Audio
```

# (1)部屋の設計
roomオブジェクトで、二次元または三次元の室を設定できます。

## シューボックス（直方体）型の部屋を作る
一般的な直方体の部屋を作るには`ShoeBox`メソッドを使います。
まず引数について触れた後使用例を示します。

ここでまとめて設定しなくても、あとからでも`room.add_microphone_array(mic_locs)`みたいな感じで設定追加できます。マイクの追加方法については後述します。

```py
import pyroomacoustics as pra

room=pra.room.ShoeBox(
    p, # 2次元[幅、奥行き], ３次元[幅、奥行き、高さ]
    fs=8000, # サンプリング周波数
    t0=0.0, # シミュレーションの開始時間
    absorption=None,  # 平均壁面吸音率（materialsを設定するときは不要）
    max_order=1, # 鏡像法における反射回数の上限 *1
    sigma2_awgn=None, # ガウシアンノイズを乗せる
    sources=None, # 音源データ(リスト)。マイク設定と一緒に使う
    mics=None,  # マイク座標のリスト（リスト）。音源設定と一緒に使う。
    materials=None, # 建材の設定。詳しくはこの後に書きます。
    temperature=None, # 室温(摂氏)。音速に影響。特に書かないと343m/s
    humidity=None, # 湿度(%)。音速に影響。
    air_absorption=False, # 空気吸収。*2
    ray_tracing=False, # レイトレーシングと鏡像法のハイブリッド利用
    use_rand_ism=False, # よくわからん。音像位置がランダムになる？
    max_rand_disp=0.0 # 上のと一緒に使う
)
```

*1 鏡像法の反射回数は`max_order=3`推奨
*2 コンサートホールのような巨大空間では音が空気に吸収され、特に高音から伝わりにくくなる。

よくやる使い方

```py
# 残響時間と部屋の寸法
rt60 = 0.5  # seconds
room_dim = [9, 7.5, 3.5]  # meters ここを二次元にすると二次平面の部屋になります

# Sabineの残響式から壁面の平均吸音率と鏡像法での反射回数の上限を決めます
e_absorption, max_order = pra.inverse_sabine(rt60, room_dim)

# 部屋をつくります
# fsは生成されるインパルス応答のサンプリング周波数です。入力する音源があるならそれに合わせる。
room = pra.ShoeBox(
    room_dim, fs=16000, materials=pra.Material(e_absorption), max_order=max_order
)
```

また次のように建材を設定してあげることもできます。
詳しくはこちらを参考に　[Matelial Database](https://pyroomacoustics.readthedocs.io/en/pypi-release/pyroomacoustics.materials.database.html)
キーワードの代わりに具体的な値を入れてあげても良いです。

```py
import pyroomacoustics as pra
m = pra.make_materials(
    ceiling="hard_surface",
    floor="6mm_carpet",
    east="brickwork",
    west="brickwork",
    north="brickwork",
    south="brickwork",
)
room = pra.ShoeBox(
    [9, 7.5, 3.5], fs=16000, materials=m, max_order=17
)
```

`room.plot()`で部屋を表示できます。

```plot.py
fig, ax = room.plot()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/543519/3f4b2ca5-369b-dd3c-32cf-5b8095352d50.png)

## 任意形状の部屋を作る
`Room.from_corners(corners)`を使えば、座標から室形状を設定してやることもできます。

まず、平面図をつくります。部屋のカドをぐるっと周るイメージで座標設定。

```corner.py
# コーナーの座標を設定します
corners = np.array([[0,0], [0,3], [5,3], [5,1], [3,1], [3,0]]).T  # [x,y]
# 部屋をつくります
room = pra.Room.from_corners(corners, materials=pra.Material(e_absorption), fs=16000, t0=0.0, max_order=1, sigma2_awgn=None, sources=None, mics=None, materials=None, **kwargs)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/543519/99c2067a-70ea-5627-e757-f8c237232297.png)

できたら、`extrude`で高さ情報を追加し、壁をつくります。

```extrude.py
room.extrude(2.)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/543519/8b6c4f7e-2ed8-fc63-0467-c4d70d5e2151.png)

# (2)マイクの設定
作った部屋`room`にマイクアレイを追加していきます。マイクアレイというのは複数のマイクからなる集音システムの意味です。

```mic.py
# マイクの座標を与える
mic_locs = np.c_[
    [6.3, 4.87, 1.2], [6.3, 4.93, 1.2],  # mic 1  # mic 2
]

# room にマイクを追加します
room.add_microphone_array(mic_locs)
```

円形にマイク配置したい場合や直線にマイクを配置したい場合、わざわざ座標を計算しなくても自動でいろんなマイク配置を計算してくれる関数が用意されています。(ただし２Dでの座標しか返してくれないので、高さ方向の座標を追加して上げる必要があります)

## （例）円形マイクアレイ

室内のマイクアレイ中心の配置場所、マイク個数、x軸からの半時計回りの回転、マイク半径を引数に入れると、それぞれのマイクの（x、y）座標を返してくれます。３次元で使いたい場合はz軸方向の座標を追加してやってください。

```circular_2D_array.py
mic_locs = pra.circular_2D_array(center=[2.,2.], M=6, phi0=0, radius=0.1)

>>>array([[2.1, 2.05, 1.95, 1.9, 1.95,2.05],
       [2., 2.08660254, 2.08660254, 2., 1.91339746,1.91339746]])
```

円形配置の他にも
linear_2D_array()
circular_2D_array()
square_2D_array()
poisson_2D_array()
spiral_2D_array()
などいろいろ用意されているので[ここから](https://pyroomacoustics.readthedocs.io/en/pypi-release/pyroomacoustics.beamforming.html#pyroomacoustics.beamforming.linear_2D_array)調べてみてください。

## マイクを一つだけ追加したい時

```mic.py
mic_loc = [1.0,2.0,2.0]
room.add_microphone(mic_loc)
```

## その座標が部屋の中にあるか調べたい時
マイクや音源をおこうとしている座標が部屋の中にちゃんと暑かどうか調べるには次のようにします。
`include_borders`は壁上も含めるかどうか。
入っていれば`True`、そうでなければ`False`を返します。

```check_inside.py
p = [1.,2.5,12.2]
room.is_inside（p、include_borders = True ）
```

# (3)音源の設定
マイクと同様に座標と、それに加えて配置したい音源のデータを与えます。音源データは自作の信号でも良いですし、手持ちのWAVファイルなどから入力することもできます。

```source.py
# wavファイルを読み込んで配置してみます
from scipy.io import wavfile
_, audio1 = wavfile.read('speech1.wav')
_, audio2 = wavfile.read('speech2.wav')
_, audio3 = wavfile.read('speech3.wav')

# 音源ごとに座標情報を与え、`room`に追加していきます。
# オプションで delay を追加することもできます。
room.add_source([2.5, 3.73, 1.76], signal=audio1, delay=1.3)
room.add_source([1.0, 2.34, 2.12], signal=audio2)
room.add_source([3.2, 1.7, 5.2], signal=audio3, delay=2.)
```

マイクと音源を追加してから`room.plot()`すると、図にも反映されます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/543519/f2f54135-f3df-7942-17c1-3cd099a9e97c.png)


# (4)シミュレーションの実行
マイクと音源が配置できたら、いよいよシミュレーションを実行しましょう。
実行はこれだけです。

```simulate.py
room.simulate()
```

実行の際にSN比の影響も加味したい場合には次のようにします。

```simulation.py
# S/N 比
SNR = 90.
room.simulate(snr=SNR)
```

# (5)シミュレーション結果

## シミュレーションされた音を聴く
配置したそれぞれのマイクに到達した音は次のようにして取り出すことができます。
サンプリング周波数は`room.fs`で取得できます。

```result.py
simulation_data = room.mic_array.signals # シミュレーション音源
display(Audio(simulation_data[0,:], rate=room.fs))
```

## インパルス応答を確認する
すべての音源からマイクロホンまでのインパルス応答を確認することができます。インパルス応答は次のように取り出します。

```rir.py
impulse_responses =room.compute_rir（）
display(Audio(impulse_responses[0][0], rate=room.fs))
```

## 室内の残響時間を計算する

```rt60.py
rt60 = room.measure_rt60()
print("残響時間：{}".format(rt60))
```

この関数はシミュレーションされたインパルス応答からだけではなく、手持ちのインパルス応答を使っても計算できます。

```rt60.py
rt60=pra.experimental.measure_rt60(impulse_responses[0][0],fs=rate)
print("残響時間：{}".format(rt60))
```

## 鏡像法を可視化する

```model.py
# compute image sources
room.image_source_model()

# visualize 3D polyhedron room and image sources
fig, ax = room.plot(img_order=3)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/543519/0ea991f7-d446-13c1-c097-7d8a09ba49a1.png)

# 終わりに

まあ、とりあえずこんなもんですかね。他にも無数のオプションがあるので、調べてみてください。
[Docs » Room Simulation](https://pyroomacoustics.readthedocs.io/en/better_doc/pyroomacoustics.room.html)

# 参考
[公式ドキュメンテーション](https://pyroomacoustics.readthedocs.io/en/better_doc/index.html)
[公式GitHub](https://github.com/LCAV/pyroomacoustics)