畳み込みは, 信号処理や画像処理, 統計学など, さまざまな分野で応用されている便利な数学演算です. 音響エンジニアは, 必要な情報を抽出したり, 音をより良く解釈するために, 畳み込みを使って音響信号を処理することがよくあります. このブログでは, COMSOL Multiphysics®ソフトウェアで畳み込みを実行するための3つの手法をご紹介します. 部屋のインパルス応答 (IR) のローパスフィルタリングにこれらの手法を適用することによる畳み込みの実装について説明し, 室内音響の可聴化の一例をお見せします.
畳み込み: 定義と手法
物理的には, 畳み込みによって, 時間領域, 周波数領域, または空間領域で表現される2つの信号間のオーバーラップ量が得られます. 時間領域の信号の場合, 数学的には次のように定義されます.
ここで, tと\tauはそれぞれ, 時間積分に使われる時間変数とダミー変数であり, \astは畳み込み演算子を表します.
各 tにおいて, この方程式は, ある関数f(\tau)の元の形式と, tによって反射されシフトされたもう一方の関数g(t-\tau)との積の時間積分を計算します. この演算は可換です. すなわち, どの関数が反射され, シフトされても結果は同じです. 積分はシフト (t)のすべての値に対して行われ, tの関数としての畳み込み結果を与えます.
この積分形式は, ここでは畳み込み積分と呼ばれ, 滑らかで連続的な関数の処理に適しています. デジタル化された音響信号のような離散的なデータの場合, この方法は非常に負荷の高い数値積分スキームを必要とし, 計算負荷が大きくなることが多くあります. これを避けるために, 離散信号の処理には次のような離散畳み込みを使用できます:
ここで, \Delta tはサンプリング間隔です.
積分を離散サンプルの和に近似することで, 離散畳み込みは畳み込み積分よりも高速に計算できます.
2つの信号のフーリエ変換が存在する場合, 畳み込み定理に基づく高速フーリエ変換 (FFT) を用いて離散畳み込みをより効率的に評価することができます:
ここで, \mathcal{F}と\mathcal{F}^{-1}はそれぞれフーリエ変換と逆フーリエ変換の演算子です. 畳み込み定理は, 時間ドメインにおける2つの関数の畳み込みのフーリエ変換が, 信号（周波数ドメインでの信号）のフーリエ変換の積と等価であることを利用しています. 最新のリアルタイム畳み込み技術は, その優れた効率性から, 一般的にFFTを使用しています.
室内インパルス応答のローパスフィルタリング
ここでは, COMSOL Multiphysics®の畳み込み積分, 離散畳み込み, 畳み込み定理を使った畳み込みの実装方法を紹介します. チュートリアルモデル小規模コンサートホール音響から得られたインパルス応答 (IR) にローパスフィルタリングを適用する例を通して, これらの方法の重要なステップをご説明します. (ローパスフィルタリングは, サウンドデザインにおいて低音だけを抽出するために広く使われています.) このIRデータをロードして, 補間関数に保存します (ここでは補間1). 以下にデータのプロットを示します.
室内IRのプロット.
データの継続時間は2秒, サンプリング周波数は44,100Hzです.
ローパスフィルターには, 4次バターワースフィルターを使用します. フィルターの伝達関数TF(f)は以下のように定義されます:
と
ここでfと\omegaはそれぞれ周波数と角周波数を表します. \omega_cはカットオフ周波数での角周波数です. \Piは, シーケンスの積を表す積演算子です.
この例では, カットオフ周波数2kHzの4次フィルターを解析関数 (この場合は解析1) で定義しています. 以下の画像は, 関数の定義方法, 関数の実数成分と虚数成分に対する関数の周波数応答, およびこのような関数で定義されたフィルターのゲイン特性を示しています.
4次Butterworth フィルターの伝達関数は, 解析関数を用いて定義されます.
実数成分と虚数成分に対する伝達関数の周波数応答.
ローパスフィルターのゲイン特性.
ゲイン特性は, ローパスフィルターが2kHzのカットオフ周波数まで平坦な通過帯域を持つことを示しています (この時点で, フィルターは入力電力を半分または3dB減衰させます). カットオフ周波数を超えると, フィルターゲインは1オクターブあたり-24dBの割合で減少します.
では, COMSOL Multiphysics®で畳み込みを行う3つの手法を見ていきましょう.
手法1: 畳み込み積分の直接処理
まず, 畳み込み積分の直接処理から始めます. ここでは, 入力として2つの時間ドメイン信号が必要です. 1つ目の信号である室内IRについては, すでに補間1として定義されています. 2つ目の信号であるローパスフィルターは, 周波数ドメインで定義されています. この信号を逆離散フーリエ変換 (IDFT) して時間ドメインに変換する必要があります. これは以下の手順で行われます.
ステップ1
Grid 1Dデータセットを作成し, Grid 1D_fと名前を付けます. これにより, 室内IRデータの周波数スペクトルに対応する周波数範囲(-22050Hz, 22050Hz)で指定された周波数グリッドが生成されます. これは, 周波数領域で信号をプロットする際に使用されます.
Grid 1D_f データセットの設定.
ステップ2
Grid 1D_fデータセットで1Dプロットグループ機能の関数プロットを使ってローパスフィルターのIDFTを実行します. 設定ウィンドウの出力セクションで, ディスプレイリストから離散フーリエ変換を選択し, 表示リストから実部を選択し, 逆変換をチェックします.
ローパスフィルターのIDFTの設定とプロット.
ステップ3
プロットデータをテーブル1に追加します.
ステップ4
テーブル1に格納されているデータを使用して補間関数 (補間2) を定義します.
補間2の設定.
これで2つの時間ドメイン信号の畳み込みの準備ができました. 定義では無限の時間間隔-\infty, \inftyが仮定されていますが, その時間範囲外では両方の入力信号がゼロに設定されるため, 積分は0–2秒の有限の時間間隔で計算するだけで済みます. 畳み込み積分の手順を以下に示します.
ステップ1
ジオメトリノードで区間を定義し, その開始値と終了値が積分区間 (0–2秒) に対応するようにします.
間隔を作ります.
ステップ2
区間上で積分演算子intop1を定義します.
積分演算子の定義.
ステップ3
元の室内IRデータのサンプリング間隔に等しい大きさの均一なラインメッシュで区間を離散化します.
間隔の離散化.
ステップ4
補間関数と intop1演算子を結果セクションで使用できるように, 定常スタディを実行します.
ステップ5
Grid 1Dデータセットを作成し, Grid 1D_tと名前を付けます. これにより, サンプリング周波数44,100Hzで, 0–2秒の範囲の時間信号を定義するための時間グリッドが生成されます.
Grid 1D_t データセットの設定.
ステップ6
Grid 1D_tをソースデータセットとして, 1Dプロットで畳み込み積分を実行します.
dest 演算子を使った畳み込み積分の式.
ここで, IRとFilter_IDFTは, 補間1と補間2で定義された室内IRとローパスフィルターIRの補間関数です.
dest演算子は, ソース点ではなく, 行先点での関数の評価を強制することにご注意ください.
手法2: 標準的な離散畳み込み
離散畳み込みはよく使用される手法です. COMSOL Multiphysics®で離散畳み込みを実行するには, アプリケーションビルダーの使用が不可欠です. アプリケーションビルダーのメソッドエディターを使用すると, モデルツリー内の操作を自動化または拡張するために実行できるメソッドを, Java®コードを使って作成することができます. このブログでは, テーブルに格納されたデータで畳み込みを実行する方法を例示します.
以下の画像は, 実装されたコードと設定を示しており, 以下の表は, コードで定義された変数を示しています.
テーブルに格納されたデータを使って離散畳み込みを実行するJavaメソッド.
| 名前 | Type | 説明 |
|---|---|---|
| a | 倍精度 (配列2D) | 1番目の入力テーブルのデータ |
| b | 倍精度 (配列2D) | 2番目の入力テーブルのデータ |
| c | 倍精度 (配列2D) | 畳み込み結果 |
| ndata1 | 整数 (スカラー) | aの長さ |
| ndata2 | 整数 (スカラー) | bの長さ |
| ndata | 整数 (スカラー) | cの長さ |
| dt | 倍精度 (スカラー) | サンプリング間隔 |
| js | 整数 (スカラー) | 合計指数の開始値 |
| je | 整数 (スカラー) | 合計指数の終了値 |
Java®プログラムで定義された変数.
このコードは, 2つの入力テーブルに格納されたデータに対して離散畳み込みを実行し, 結果を出力テーブルにエクスポートします. コードをよりよく理解するために, 以下のメモをご参照ください:
- 2–4行目: 1番目の入力テーブルからデータとデータ長を抽出する
- 6–8行目: 2番目の入力テーブルからデータとデータ長を抽出する
- 11–13行目: 結果の長さを定義し, 畳み込み結果を格納する倍精度配列を作成し, サンプリング間隔を定義する
- 18–28行目: forループを実行し, 最初の時間ステップから最後の時間ステップまで, 離散畳み込みを開始する
- 30–33行目: 結果を, 離散畳み込み結果とラベル付けされた出力テーブルにエクスポートする
結果データの長さは2つの入力テーブルの長さの和から1を引いたものであり, 和のインデックスの開始値 (js) と終了値 (je) は, jsが2番目のテーブルの長さより小さく, jeが1番目のテーブルの長さより小さくなるように定義されていることにご注意ください (それぞれコードの22行目と23行目を参照).
プログラムを実行するには, 2つの時間信号のプロットをテーブル格納する必要があります. ローパスフィルターのIDFTのプロットデータはテーブル1に格納されています. 室内IRについては, Grid 1D_tデータセットを使用して1Dプロットグループ機能で室内IRの関数プロットを作成し, プロットデータをテーブル2にコピーします.
テーブルのタグ名は, モデルツリーに追加されたメソッド呼出機能の設定ウィンドウで入力できます. すべての設定が完了したら, メソッド呼出設定ウィンドウの 実行ボタンをクリックして, 離散畳み込みを実行できます.
メソッド呼出 機能の 設定 ウィンドウでテーブルタグを入力.
方法3: 畳み込み定理に基づく離散畳み込み
最後の手法は, 畳み込み定理に基づいて (すなわちDFTを使用して) 畳み込みを実行します. この場合, 室内IRの DFT とローパスフィルターの伝達関数が使用されます. その方法は以下の通りです:
ステップ1
室内 IR の DFT の実数部と虚数部の両方を, 同じ1Dプロットグループ内の2つの関数プロット (それぞれ1つずつ) にプロットします. 設定で, 表示リストから離散フーリエ変換を選択します. 次に, 表示リストから, 室内IRの DFT の実数部と虚数部をプロットするために, それぞれ実数部または虚数部を選択します. チェックボックスはデフォルト設定です.
室内IRの DFT を1Dプロットしたもの.
ステップ2
室内 IR の DFT の実数部と虚数部のプロットデータをそれぞれテーブル3とテーブル4にコピーします.
ステップ3
テーブル3とテーブル4にそれぞれ保存されたデータを使って, 補間3と補間4を定義します.
ステップ4
室内 IR とローパスフィルターの DFT の点積の IDFT を計算し, 畳み込みを実行します. 点ごとの乗算と IDFT は, Grid 1D_fデータセットを使用して, 単一の関数プロットで実行されます.
室内 IR の DFT とローパスフィルターの積の IDFT. 逆変換を実行するために, 設定 ウィンドウの逆変換チェックボックスがオンになっています.
ここで, real_IRとimag_IRは, 室内IRのDFTの実数部と虚数部であり, それぞれ補間3 と補間4で定義されます. Butterworthは, 解析1で定義されたローパスフィルターの伝達関数です.
畳み込み定理手法では循環畳み込みが実行されることにご注意ください. つまり, グリッドの間隔の長さが不十分な場合, 循環的な重なりが生じます.
ローパスフィルタリングの結果
下のグラフは, 3つの手法で計算された室内IRをローパスフィルターしたもので, 3つの手法がよく一致しているのがわかります.
畳み込み積分, 離散畳み込み, および畳み込み定理によって計算された畳み込み結果の比較.
畳み込みの結果は0.02–0.07秒の範囲になります.
次の図は, フィルタリング前後の室内IRのスペクトルの比較です. フィルタリング後のスペクトルは, ローパスフィルターのカットオフ周波数である2kHzより上で減少していることが確認できます. この結果は, 畳み込み実装の成功を裏付けています.
フィルタリング前後の室内インパルス応答のスペクトル.
可聴化への応用
可聴化の例に移りましょう. この例では, 室内 IR と無響室で録音された音声の畳み込みを取り上げます. 線形時不変システムの仮定の下で, 任意の入力からの応答はIRと入力信号の畳み込みによって評価することができます. この理論に基づき, 音響学者は計算手法でシミュレートした IR と無響音の畳み込みによって音場を聴感的に評価します. このシミュレーションのプロセスを可聴化と呼び, 音声データの作成から聴取に至るまでの過程に及びます. このモデル例 (次のセクションのボタンからアクセス可能) では, 前述の標準的な離散畳み込み法を用いて, 小規模コンサートホールモデルでトランペットの音の可聴化を行っています. 畳み込みの結果は, WAV オーディオファイルとしてエクスポートして, オーディオプレーヤーまたはメディアプレーヤーで再生することもできます. 以下では, 無響室でのトランペットの音とリバーブのかかったトランペットの音を聴き比べることができます.
無響室で演奏したときのトランペットの音.
可聴化を使い, 小さなコンサートホールで演奏したときのトランペットの音.
上記の例はモノラル音であり, 私たちが通常耳にするバイノーラル音とは異なります. しかし, 頭部伝達関数やアンビソニックスのような音場再生技術を組み合わせることで, よりリアルな音を生成することができます (例えば, 文献1参照). このような例については, 今後のブログでご紹介する予定です.
次のステップ
このブログで紹介したモデルをさらに詳しく見たい場合は, 下のボタンをクリックすると, アプリケーションギャラリに移動できます.
さらに詳しい文献
音響モデリングにおけるデータ処理については, 以下のブログで詳しく説明しています:
参考文献
- M. Vorländer, Auralization: Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality, Springer Science & Business Media: Berlin/Heidelberg, 2007.
無響音はThe Open AIR Libraryより,CC BY 4.0で提供されています.
Oracle および Java は, Oracle および/またはその関連会社の登録商標です.
コメント (0)