はじめに
1. 鏡像法とは
- 1.1 計算のイメージ
- 1.2 CPUでのシミュレーションの限界と、GPUの利点
2. gpuRIRとは
3. 環境構築
- 3.1 必要なソフトウェア
- 3.2 Windows11での利用方法
4. gpuRIRの基本的な使い方
4.4 反射係数とシミュレーションパラメータの計算
- 4.5 インパルス応答の生成
- 4.6 結果のプロット
まとめ
関連記事：PythonでCPUを使った鏡像法によるRIR生成

はじめに

この記事では、Python用ライブラリ「gpuRIR」を使って、GPUで高速に室内インパルス応答（RIR: Room Impulse Response）を生成する方法を紹介します。

鏡像法でよりリッチなRIR（例：長い残響）を生成しようとすると、CPUでは計算時間が長くなってしまいます。

そこでgpuRIRは、鏡像法をGPUで実行することで、従来のCPUでの処理よりも大幅に高速化できます。室内音響シミュレーションは、バーチャルリアリティやゲーム開発、音響機器の設計において不可欠です。gpuRIRの使い方を理解することで、これらの分野で効率的に音響シミュレーションや残響生成用のフィルタ作成を行うことができます。具体的には、コンボリューションリバーブ（畳み込みリバーブ）のフィルタを生成することができます。

なお、GPU非搭載PCでPythonを使った鏡像法によるRIR生成は、以下の記事を参考にしてください。 www.wizard-notes.com

1. 鏡像法とは

reuk.github.io

鏡像法（Image Source Method）は、室内で音が壁や床、天井に反射する様子をシミュレートするための手法です。音源が壁に反射すると、あたかも壁の向こうに「鏡像」のように音源があると仮定します。この仮想的な音源を使って、反射音がリスナーにどのように届くかを計算します。これを全ての壁、天井、床について繰り返すことで、多重反射による残響をシミュレートできます。

1.1 計算のイメージ

https://reuk.github.io/wayverb/image_source.html Figure1 より引用

部屋の壁を鏡 (boundary) と考えて、音源 (source) が壁に映るように仮想音源 (image source) を生成します。
仮想音源からリスナー（マイク, receiver）までの距離と音の減衰を計算することで、音の伝播や反射の影響を表現します。

1.2 CPUでのシミュレーションの限界と、GPUの利点

鏡像法を使ったシミュレーションでは、多数の仮想音源（反射音）の計算に時間がかかります。

RIRを生成する際、以下の要因が計算量を増加させます。

部屋のサイズ
- 部屋が広いほど、音が壁に反射してリスナーに届くまでの距離が長くなります。その結果、計算する必要のある仮想音源（鏡像音源）の数が増加します。
- 広い部屋では、遠くにある反射音まで含めてシミュレーションする必要があり、これが計算量を増やします。
反射回数（仮想音源の数）
- 多重反射の考慮: 鏡像法では、壁、床、天井などに対する反射ごとに仮想音源を生成します。例えば、1回反射だけでなく、2回、3回と多重反射を考慮する場合、仮想音源の数が急激に増加します。
- 反射回数が増えると、仮想音源の位置やそれに対応する伝播時間、減衰の計算も複雑になるため、計算量が増えます。
音源・マイクの数
- 音源やマイクの数が多いほど、それぞれの音源からマイクへの伝播経路を計算する必要があり、計算量が増加します。例えば、部屋内に2つの音源と3つのマイクがある場合、それぞれの組み合わせに対してRIRを計算する必要があります。
- 音源とマイクの数が増えると、必要な計算量が線形的に増加します。
サンプリング周波数
- 高いサンプリング周波数では、1秒間に生成するサンプルの数が増加します。これにより、インパルス応答の長さ（時間方向のデータ量）が増え、各仮想音源からリスナーまでの伝播を高精度にシミュレートするための計算量が増えます。
残響時間（T60）
- 長い残響時間では、音が部屋の中で減衰するまでの時間が長くなるため、シミュレーションの範囲も広がります。これにより、多くの反射を考慮してRIRを生成する必要が生じ、計算量が増加します。
吸音特性の複雑さ
- 部屋の壁や天井、床の吸音特性が複雑な場合、音の減衰や拡散をより詳細に計算する必要があります。各面の異なる吸音特性を考慮して、音の反射と減衰をモデル化するために多くの計算を行う必要があります。

これらの要素により、反射音の数と各反射音に関する伝播時間、減衰量、位相の計算が増えるため、計算量が大幅に増加します。

CPUでは、並列処理の数に限りがあるため、このような大規模なシミュレーションには時間がかかることが問題となります。そこで、gpuRIRのようなGPUを活用したライブラリを使うと、これらの計算を並列化して高速に処理できるため、複雑なRIRを生成する際に非常に有効です。

2. gpuRIRとは

gpuRIRは、Python環境でGPUを利用して室内音響シミュレーションを行うためのライブラリです。具体的には、鏡像法を用いて、室内の壁や床などでの音の反射をシミュレートし、インパルス応答を生成します。

公式リポジトリには、gpuRIRのインストール手順や基本的な使い方が記載されています。本記事では、それを参考に導入方法を説明します。

github.com

論文はこちらです：https://arxiv.org/pdf/1810.11359

Diaz-Guerra, D., Miguel, A., & Beltran, J. R. (2021). gpuRIR: A python library for room impulse response simulation with GPU acceleration. Multimedia Tools and Applications, 80(4), 5653-5671.

3. 環境構築

3.1 必要なソフトウェア

CUDA対応のGPUドライバ
Python（3.x）
gpuRIRのインストールは以下のコマンドで行います。

pip install https://github.com/DavidDiazGuerra/gpuRIR/zipball/master

OSは，GNU/Linux システム (Ubuntu および centOS) および Windows 10 で動作確認済みとなっています。

注意点: CUDA対応のGPUとそのドライバが必要です。また、gpuRIRは特定のPythonおよびCUDAのバージョンに依存するため、公式ドキュメントの互換性情報を確認してください

3.2 Windows11での利用方法

Windowsユーザーの場合、前回の記事で説明したように、WSL2上でCUDAとPython環境を構築する必要があります。特にgpuRIRはGPUを利用するため、WSL2でのCUDA環境設定が重要です。

www.wizard-notes.com

4. gpuRIRの基本的な使い方

以下に、部屋のサイズ、音源とマイクの位置、吸音係数などを指定してインパルス応答を生成する基本的なコード例を紹介します。

解説するコードは、公式の examples.py です。

4.1 ライブラリのインポート、初期設定

import numpy as np
import numpy.matlib
import matplotlib
import matplotlib.pyplot as plt
from math import ceil

import gpuRIR
gpuRIR.activateMixedPrecision(False)  # 混合精度計算を無効にします
gpuRIR.activateLUT(True)  # ルックアップテーブルを有効にします

activateMixedPrecision(False):
- 混合精度の計算を無効にして、計算の精度を維持します。
activateLUT(True):
- ルックアップテーブル（LUT）を使用して計算の高速化を行います。

4.2 シミュレーションのパラメータ設定

room_sz = [3, 3, 2.5]  # 部屋のサイズ [m]
nb_src = 2  # 音源の数
pos_src = np.array([[1, 2.9, 0.5], [1, 2, 0.5]])  # 音源の位置 [m]
nb_rcv = 3  # 受音点（マイク）の数
pos_rcv = np.array([[0.5, 1, 0.5], [1, 1, 0.5], [1.5, 1, 0.5]])  # 受音点の位置 [m]
orV_rcv = np.matlib.repmat(np.array([0, 1, 0]), nb_rcv, 1)  # 受音点が向いている方向のベクトル
mic_pattern = "card"  # マイクの指向性パターン

音源とマイクの位置関係は、天井から見下ろすと以下のようになっています

fig, ax = plt.subplots(figsize=(6, 6))
for k in range(nb_src):
    ax.plot(pos_src[k,0], pos_src[k,1], label=f"src {k}", marker="o")
    
for k in range(nb_rcv):
    ax.plot(pos_rcv[k,0], pos_rcv[k,1], label=f"rcv {k}", marker="*")
    
ax.set_xlabel("x")
ax.set_ylabel("y")
plt.xlim([0,room_sz[0]])
plt.ylim([0,room_sz[1]])
plt.legend()
plt.show()

4.3 室内音響特性の設定

abs_weights = [0.9]*5 + [0.5]  # 壁の吸音係数の割合
T60 = 1.0  # 残響時間（60dB減衰するまでの時間） [s]
att_diff = 15.0  # 拡散反射モデルを使用し始めるときの減衰量 [dB]
att_max = 60.0  # シミュレーション終了時の減衰量 [dB]
fs = 16000.0  # サンプリング周波数 [Hz]

abs_weights:
- 部屋の壁の吸音係数を設定。[0.0, 1.0] の間で値を設定
- [0.9]*5 + [0.5] は、6つの壁（4つの側面、天井、床）に対する吸音特性を表しており、最初の5つの壁は0.9、最後の1つは0.5の吸音係数を持ちます。
T60:
- 音が60dB減衰するまでの時間（残響時間）。
att_diff, att_max:
- シミュレーションの際の拡散反射の開始タイミングや終了時の減衰量を設定。

4.4 反射係数とシミュレーションパラメータの計算

beta = gpuRIR.beta_SabineEstimation(room_sz, T60, abs_weights=abs_weights)  # 反射係数の推定
Tdiff = gpuRIR.att2t_SabineEstimator(att_diff, T60)  # 拡散反射モデルを開始する時間 [s]
Tmax = gpuRIR.att2t_SabineEstimator(att_max, T60)  # シミュレーション終了時間 [s]
nb_img = gpuRIR.t2n(Tdiff, room_sz)  # 各次元における鏡像音源の数

各関数の説明：

beta_SabineEstimation():
- サビンの公式を使用して、部屋の反射係数を推定します。
att2t_SabineEstimator():
- 減衰量に基づき、拡散反射モデルの開始時間やシミュレーションの終了時間を計算します。
t2n():
- 時間と部屋のサイズから、各次元における鏡像音源の数を計算します。

各関数の返り値：

beta
- 6要素の array_like
- 各壁の反射係数を格納
Tdiff :
- float、opitional
- 鏡像法が拡散残響モデルに置き換えられる時間 (秒単位)。
- デフォルトは Tmax (完全な ISM シミュレーション)
Tmax :
- float
- RIR の長さ (秒単位)。
nb_img :
- 3 つの整数要素を持つ array_like
- 各次元でシミュレートする仮想音源の数

4.5 インパルス応答の生成

RIRs = gpuRIR.simulateRIR(room_sz, beta, pos_src, pos_rcv, nb_img, Tmax, fs, Tdiff=Tdiff, orV_rcv=orV_rcv, mic_pattern=mic_pattern)

simulateRIR() の引数として，部屋のサイズ、反射係数、音源と受音点の位置、鏡像音源の数、シミュレーション時間などを設定し、インパルス応答を計算します

ここで、RIRの形状は以下のようになっています：（音源、受音点、時間長）

RIRs.shape  # (nb_src, nb_rcv, num_samples)

4.6 結果のプロット

t = np.arange(int(ceil(Tmax * fs))) / fs
plt.plot(t, RIRs.reshape(nb_src * nb_rcv, -1).transpose())
plt.show()

np.arangeでシミュレーション時間に基づいた時間軸を生成しています。

しかし、このプロットだとよくわからないので、以下のコードで音源ごとにプロットしてみます。また、プロットする時間長を変えてみます。

t = np.arange(int(ceil(Tmax * fs))) / fs
length = 160*3
plt.figure(figsize=(8,8))
for k_src in range(nb_src):
    print(k_src, nb_src)
    plt.subplot(nb_src, 1, k_src+1)
    for k_rcv in range(nb_rcv):
        plt.plot(t[0:length], RIRs[k_src,k_rcv,0:length].transpose(), label=f"src={k_src}, rcv={k_rcv}", alpha=0.5)
    plt.legend() 
plt.show()