鏡像音源法による直方体形状の部屋の伝送特性計算法

鏡像音源法による直方体形状の部屋の伝送特性計算法
●原理と先行研究例、問題点の洗い出し反射音を鏡像音源として計算することは誰でもすぐに思いつくことであり、実際に筆者が1997年に計算を始めた時も全く白紙の状態から独自に手法を開発している。(当時はWebで研究論文を検索することは普及していなかったし、通信手段もアナログ回線にModem接続であった。) 但し、これは「車輪の再発明」であり、コンピューターが利用できるようになった初期の頃から様々な研究論文が発表されていたのである。・A simple image method for calculating the distribution of sound pressure level within an enclosure (1972) 　：（詳細不明だが、原点となる研究らしい）・Image Method for efficiently simulating small‐room acoustics (1978) 　：モデル化と実際の計算アルゴリズムが紹介されている。・Rigid sphere room impulse response simulation: Algorithm and applications (2012) 　：鏡像音源法からインパルスレスポンスを求める、さらに発展した手法の解説。壁面の一部でも不整形な部分があると、反射回数が少なくても鏡を通して見える像がゆがみ、複数回の反射による鏡像音源の３次元的な座標を計算するのは非常に困難になる。不整形な反射面を持つ空間の計算には、光源から発射された光を微小角に分割して１本１本の反射して行く先を追跡する、「レイ・トレーシング法」が有効となる。しかし、条件を完全な直方体形状の部屋に限定した場合は、単純な鏡像音源法が最適である。直方体なら反射によって測定点から見える鏡像空間は単純にブロックを隙間なく積み上げた形となり、鏡像による仮想音源はそれぞれの鏡像空間に１個ずつなので、鏡像空間も音源も互いに重なることはない。部屋の寸法をW,H,Dとして、直方体内部の座標はX,Y,Zの各方向に１回反射するたびに座標の向きが反転する。鏡像音源位置は各直方体の原点に対し、音源座標XYZを反射回数に応じて反転するだけで容易に決定できる。伝送周波数特性は、これらの音源が十分小さいと見なせるレベルに減衰するまでの反射回数の全鏡像音源からの音圧合成値を計算することで得られるが、反射回数が増えると鏡像音源の数は反射回数の３乗で増加し、50回で100万個を超える。計算量の制限により反射を全て計算するのは困難で、計算機パワーのなかった昔は十数回の反射までしか計算できなかったが、2000年以降では家庭用のＰＣでも十分な能力を持つようになり、現在では100回の反射回数でも計算可能となった。それぞれの鏡像音源から測定点に到達する音圧は、測定点からの距離による減衰と、壁面での反射損失による減衰により反射回数が多いほど小さくなる。（何枚もの半透明窓を通して見る遠い灯りは、その分暗くなるのと同じ。）時刻０で音源から単位インパルス（距離１ｍで音圧１となるレベル）が放射されたとして、鏡像音源の距離遅延と減衰量を計算し、全音源からの反射波を加算した結果がインパルス応答となる。連続正弦波に対する伝送特性は、このインパルス応答波形を音源波形（測定周波数）に対し畳み込み積分すれば求められる。（周波数特性を求めるには、インパルス応答をＦＦＴしても良い。）最終的な誤差を１％以下とするには、反射減衰量が60dBになるまでの全ての鏡像音源の計算が必要であり、必要な反射回数は壁面の反射率によって左右される。反射率0.8（減衰2dB）では30回だが、0.95（減衰0.45dB）では133回が必要となり、約1300万個の鏡像音源を計算しなければならない。　単純な立方体形状の計算では、問題は「いかに計算を早くするか」だけであり、早くするために発生する誤差を考慮しながら必要な精度と計算量をバランスさせる工夫が必要となる。 ●ＤＣ増幅のメカニズムと回避方法反射音のインパルス応答を計算する際、反射減衰が少なくて非常に大きい反射回数を設定すると、問題が発生する。それは、反射率が１の場合に無限の反射回数を設定するとＤＣ成分が無限大に発散してしまう現象である。これについての詳細な検討結果は、こちら。→（ＤＣ成分発生の解説へのリンク）インパルス入力以外の時間に均一にＤＣ成分を打ち消す信号を設定すれば対策可能であるが、通常の条件ではその影響は無視できる。なお、１ｋＨｚ以上の高い周波数で数秒以上の長距離の伝播では空気自身の損失で距離減衰が大きくなるが、比較的狭い室内での低域（300Hz以下）での特性を論じる場合は無視して全く問題ない。 ●シミュレーションの高速化と精度最適化鏡像音源法では、大きく分けると次の２つの誤差発生要因がある。１、反射回数打切りによる残響音欠落２、インパルス応答計算の際の時間分解能による時間軸誤差それぞれの要因に対し、必要な精度の実現に必要な基本部分の維持と、不要に時間を掛けている部分のそぎ落としが全体最適の条件となる。・反射回数打切りによる誤差量と最適化個々の鏡像音源レベルは、距離減衰と反射減衰の積で決まる。反射回数は時間経過に比例して一様に増加して行くので、反射減衰(dB)は時間に比例する。距離減衰は１個の音源あたりでは1/Lだが、Lが２倍に増えると２倍の距離の面積は４倍（Lの２乗）になるので、時間当たりの音源密度は４倍となり、距離減衰を打ち消して距離に比例して逆に増大する。この時、ＤＣ成分は単純に増加するが、変動(AC)成分はエネルギー加算で√となるので増加せず、一定値となる。つまり反射音構造（インパルス応答）の時間経過は、距離減衰は数の効果で打ち消され、反射減衰のみによって時間経過につれてほぼ一様に低下する。（この値が初期値の-60dBになる所が「残響時間」である。）全体がLogで一様に減衰する場合は、無限小までの積分値は切り捨て時のレベルに単純比例するので、初期値から-40dBに減衰した時点で打ち切れば-∞まで計算した場合との誤差は-40dBとなる。１％精度で計算するには残響音が-40dBになるまでの鏡像音源を全て加算すれば良い。下の図は、天井の低い１２畳の部屋(W=5.4m, D=3.6m, H=2.4m)での1ms毎の反射音レベルの時間構造グラフ（反射率は床以外0.8）の例である。実際の残響音は多数の反射音の重なりなので、切捨てレベルの約10倍から誤差が発生してくる。切り捨て誤差を１％以下にするには、切捨て時反射率を0.001とする必要がある。（誤差が３％を超えると周波数特性の計算結果にリプルが現れてくる）反射回数をＸＹＺの３方向で一定値までとすると、立方体の角の部分はＸＹＺの全ての反射減衰が加わるので１方向のみの減衰量の３倍になる。つまり角の部分は計算不要。「反射減衰が合計40dBまで」で制限すれば、計算する範囲は正方形を斜めに切った四角錐形状になるので、全体の個数は立方体の1/6に削減できる。・時間軸精度による誤差量と改善手法最初にインパルス応答を求めてFFTで周波数特性を求める方法は、計算時間が飛躍的に短縮できる半面、時間軸丸めによる誤差が新たに発生する欠点がある。その誤差量は、ΔP=２πｆΔｔ/２（時間分解能をΔｔとすると、最大誤差はΔｔ/２となるため）であり、周波数に比例したランダムにうねる傾向の誤差となる。例えば、時間軸分解能が0.1msの場合には300Hzでの誤差は最大10％近くになる。300Hzに於いて誤差を１％以下にするには、時間分解能は0.01msにしなければならない。この時のFFTのポイント数は１秒で100,000個にもなってしまうが、13万ﾎﾟｲﾝﾄまで高速にFFT計算できるエンジンを導入すればとりあえずは解決できる。発生する誤差の量は１個１個の反射音レベルに比例するので、直接音と第一反射音が主要因であり、レベルの低い残響音に同じ時間分解能を割り当てるのは非常に無駄が多い。初期反射音のレベルの高い部分は、時間軸丸めなしで積和計算を１個１個直接行う（インパルス応答に含めない）事で計算量を節約すると同時に精度を向上させることが可能となる。初期反射音は時間密度が低いので、むやみに分解能を上げても無音部分が多くなり、これも無駄なのである。反射４回まで（729個）を直接計算した場合、時間軸の丸めを行うのは反射５回以上となり、距離減衰と合わせて誤差を1/20以下に出来る。さらに、５回の反射では0.1ms当りの反射音密度は５個程度になるので誤差は加算平均されてさらに小さくなることが期待できる。つまり、ｲﾝﾊﾟﾙｽ応答で計算するのを反射回数５回以上成分に限定することで時間分解能0.1msのままでも誤差は0.3%以下にまで改善出来、さらに反射遅延100ms以降は時間分解能は0.5msに下げる事で精度を保ったまま計算量をもう１段節約することもできる。（FFTでダイレクトに周波数特性を求めることは出来ないので計算時間が改善するかは微妙だが、時間軸誤差による余計な誤差要因を低減できる効果は大きい）・計算アルゴリズムの高速化反射回数が多くなると、最初にインパルス応答を求めるだけで100万個以上の鏡像音源の座標と測定点からの距離、反射減衰量を計算しなければならない。計算機では、加減算→掛算→割算→階乗（指数）計算の順で速度は遅くなる。計算アルゴリズムで階層を何段にも分けてループ計算させるが、最も深い階層（時間的に90%以上を占める）で計算量を最小限にすると高速化できる。座標軸の前後を反射毎に反転させるために(-1)^nの階乗が最深ループに含まれていると遅くなるのである。さらに、ＸＹＺの３方向の座標・反射率の組合せを計算する前に、ＸＹＺそれぞれの反射回数毎の座標・反射率のテーブルを先に計算しておくと良い。つまり、反射回数100回でも、±方向別で200個のテーブルをＸＹＺで作れば、最後に100万回のループを回す時に計算済みのテーブルを参照するだけになり、高速化できる。３次元的な距離は√(X²+Y²+Z²)で計算されるが、二乗の計算はテーブルの距離設定で先に済ませられるのである。参考として、単発の伝送周波数特性を計算できるExcelシートを公開する。（マクロを含むので使用は自己責任で）　　次に進む　　室内の音響伝達特性検討の目次に戻る　　オーディオのTOPに戻る

　
●原理と先行研究例、問題点の洗い出し

反射音を鏡像音源として計算することは誰でもすぐに思いつくことであり、実際に筆者が1997年に計算を始めた時も全く白紙の状態から独自に手法を開発している。(当時はWebで研究論文を検索することは普及していなかったし、通信手段もアナログ回線にModem接続であった。)

但し、これは「車輪の再発明」であり、コンピューターが利用できるようになった初期の頃から様々な研究論文が発表されていたのである。

・A simple image method for calculating the distribution of sound pressure level within an enclosure (1972)
　：（詳細不明だが、原点となる研究らしい）
・Image Method for efficiently simulating small‐room acoustics (1978)
　：モデル化と実際の計算アルゴリズムが紹介されている。
・Rigid sphere room impulse response simulation: Algorithm and applications (2012)
　：鏡像音源法からインパルスレスポンスを求める、さらに発展した手法の解説。

壁面の一部でも不整形な部分があると、反射回数が少なくても鏡を通して見える像がゆがみ、複数回の反射による鏡像音源の３次元的な座標を計算するのは非常に困難になる。
不整形な反射面を持つ空間の計算には、光源から発射された光を微小角に分割して１本１本の反射して行く先を追跡する、「レイ・トレーシング法」が有効となる。

しかし、条件を完全な直方体形状の部屋に限定した場合は、単純な鏡像音源法が最適である。
直方体なら反射によって測定点から見える鏡像空間は単純にブロックを隙間なく積み上げた形となり、鏡像による仮想音源はそれぞれの鏡像空間に１個ずつなので、鏡像空間も音源も互いに重なることはない。

部屋の寸法をW,H,Dとして、直方体内部の座標はX,Y,Zの各方向に１回反射するたびに座標の向きが反転する。鏡像音源位置は各直方体の原点に対し、音源座標XYZを反射回数に応じて反転するだけで容易に決定できる。

伝送周波数特性は、これらの音源が十分小さいと見なせるレベルに減衰するまでの反射回数の全鏡像音源からの音圧合成値を計算することで得られるが、反射回数が増えると鏡像音源の数は反射回数の３乗で増加し、50回で100万個を超える。
計算量の制限により反射を全て計算するのは困難で、計算機パワーのなかった昔は十数回の反射までしか計算できなかったが、2000年以降では家庭用のＰＣでも十分な能力を持つようになり、現在では100回の反射回数でも計算可能となった。

それぞれの鏡像音源から測定点に到達する音圧は、測定点からの距離による減衰と、壁面での反射損失による減衰により反射回数が多いほど小さくなる。（何枚もの半透明窓を通して見る遠い灯りは、その分暗くなるのと同じ。）
時刻０で音源から単位インパルス（距離１ｍで音圧１となるレベル）が放射されたとして、鏡像音源の距離遅延と減衰量を計算し、全音源からの反射波を加算した結果がインパルス応答となる。

連続正弦波に対する伝送特性は、このインパルス応答波形を音源波形（測定周波数）に対し畳み込み積分すれば求められる。（周波数特性を求めるには、インパルス応答をＦＦＴしても良い。）
最終的な誤差を１％以下とするには、反射減衰量が60dBになるまでの全ての鏡像音源の計算が必要であり、必要な反射回数は壁面の反射率によって左右される。
反射率0.8（減衰2dB）では30回だが、0.95（減衰0.45dB）では133回が必要となり、約1300万個の鏡像音源を計算しなければならない。
　
単純な立方体形状の計算では、問題は「いかに計算を早くするか」だけであり、早くするために発生する誤差を考慮しながら必要な精度と計算量をバランスさせる工夫が必要となる。

●ＤＣ増幅のメカニズムと回避方法

反射音のインパルス応答を計算する際、反射減衰が少なくて非常に大きい反射回数を設定すると、問題が発生する。それは、反射率が１の場合に無限の反射回数を設定するとＤＣ成分が無限大に発散してしまう現象である。

これについての詳細な検討結果は、こちら。→（ＤＣ成分発生の解説へのリンク）

インパルス入力以外の時間に均一にＤＣ成分を打ち消す信号を設定すれば対策可能であるが、通常の条件ではその影響は無視できる。
なお、１ｋＨｚ以上の高い周波数で数秒以上の長距離の伝播では空気自身の損失で距離減衰が大きくなるが、比較的狭い室内での低域（300Hz以下）での特性を論じる場合は無視して全く問題ない。

●シミュレーションの高速化と精度最適化

鏡像音源法では、大きく分けると次の２つの誤差発生要因がある。
１、反射回数打切りによる残響音欠落
２、インパルス応答計算の際の時間分解能による時間軸誤差

それぞれの要因に対し、必要な精度の実現に必要な基本部分の維持と、不要に時間を掛けている部分のそぎ落としが全体最適の条件となる。

・反射回数打切りによる誤差量と最適化

個々の鏡像音源レベルは、距離減衰と反射減衰の積で決まる。反射回数は時間経過に比例して一様に増加して行くので、反射減衰(dB)は時間に比例する。

距離減衰は１個の音源あたりでは1/Lだが、Lが２倍に増えると２倍の距離の面積は４倍（Lの２乗）になるので、時間当たりの音源密度は４倍となり、距離減衰を打ち消して距離に比例して逆に増大する。この時、ＤＣ成分は単純に増加するが、変動(AC)成分はエネルギー加算で√となるので増加せず、一定値となる。

つまり反射音構造（インパルス応答）の時間経過は、距離減衰は数の効果で打ち消され、反射減衰のみによって時間経過につれてほぼ一様に低下する。（この値が初期値の-60dBになる所が「残響時間」である。）

全体がLogで一様に減衰する場合は、無限小までの積分値は切り捨て時のレベルに単純比例するので、初期値から-40dBに減衰した時点で打ち切れば-∞まで計算した場合との誤差は-40dBとなる。１％精度で計算するには残響音が-40dBになるまでの鏡像音源を全て加算すれば良い。

下の図は、天井の低い１２畳の部屋(W=5.4m, D=3.6m, H=2.4m)での1ms毎の反射音レベルの時間構造グラフ（反射率は床以外0.8）の例である。
実際の残響音は多数の反射音の重なりなので、切捨てレベルの約10倍から誤差が発生してくる。切り捨て誤差を１％以下にするには、切捨て時反射率を0.001とする必要がある。（誤差が３％を超えると周波数特性の計算結果にリプルが現れてくる）

反射回数をＸＹＺの３方向で一定値までとすると、立方体の角の部分はＸＹＺの全ての反射減衰が加わるので１方向のみの減衰量の３倍になる。つまり角の部分は計算不要。
「反射減衰が合計40dBまで」で制限すれば、計算する範囲は正方形を斜めに切った四角錐形状になるので、全体の個数は立方体の1/6に削減できる。

・時間軸精度による誤差量と改善手法

最初にインパルス応答を求めてFFTで周波数特性を求める方法は、計算時間が飛躍的に短縮できる半面、時間軸丸めによる誤差が新たに発生する欠点がある。
その誤差量は、ΔP=２πｆΔｔ/２（時間分解能をΔｔとすると、最大誤差はΔｔ/２となるため）であり、周波数に比例したランダムにうねる傾向の誤差となる。

例えば、時間軸分解能が0.1msの場合には300Hzでの誤差は最大10％近くになる。300Hzに於いて誤差を１％以下にするには、時間分解能は0.01msにしなければならない。この時のFFTのポイント数は１秒で100,000個にもなってしまうが、13万ﾎﾟｲﾝﾄまで高速にFFT計算できるエンジンを導入すればとりあえずは解決できる。

発生する誤差の量は１個１個の反射音レベルに比例するので、直接音と第一反射音が主要因であり、レベルの低い残響音に同じ時間分解能を割り当てるのは非常に無駄が多い。
初期反射音のレベルの高い部分は、時間軸丸めなしで積和計算を１個１個直接行う（インパルス応答に含めない）事で計算量を節約すると同時に精度を向上させることが可能となる。初期反射音は時間密度が低いので、むやみに分解能を上げても無音部分が多くなり、これも無駄なのである。

反射４回まで（729個）を直接計算した場合、時間軸の丸めを行うのは反射５回以上となり、距離減衰と合わせて誤差を1/20以下に出来る。さらに、５回の反射では0.1ms当りの反射音密度は５個程度になるので誤差は加算平均されてさらに小さくなることが期待できる。

つまり、ｲﾝﾊﾟﾙｽ応答で計算するのを反射回数５回以上成分に限定することで時間分解能0.1msのままでも誤差は0.3%以下にまで改善出来、さらに反射遅延100ms以降は時間分解能は0.5msに下げる事で精度を保ったまま計算量をもう１段節約することもできる。
（FFTでダイレクトに周波数特性を求めることは出来ないので計算時間が改善するかは微妙だが、時間軸誤差による余計な誤差要因を低減できる効果は大きい）

・計算アルゴリズムの高速化

反射回数が多くなると、最初にインパルス応答を求めるだけで100万個以上の鏡像音源の座標と測定点からの距離、反射減衰量を計算しなければならない。
計算機では、加減算→掛算→割算→階乗（指数）計算の順で速度は遅くなる。

計算アルゴリズムで階層を何段にも分けてループ計算させるが、最も深い階層（時間的に90%以上を占める）で計算量を最小限にすると高速化できる。
座標軸の前後を反射毎に反転させるために(-1)^nの階乗が最深ループに含まれていると遅くなるのである。

さらに、ＸＹＺの３方向の座標・反射率の組合せを計算する前に、ＸＹＺそれぞれの反射回数毎の座標・反射率のテーブルを先に計算しておくと良い。つまり、反射回数100回でも、±方向別で200個のテーブルをＸＹＺで作れば、最後に100万回のループを回す時に計算済みのテーブルを参照するだけになり、高速化できる。

３次元的な距離は√(X²+Y²+Z²)で計算されるが、二乗の計算はテーブルの距離設定で先に済ませられるのである。

参考として、単発の伝送周波数特性を計算できるExcelシートを公開する。（マクロを含むので使用は自己責任で）

　　次に進む　　室内の音響伝達特性検討の目次に戻る　　オーディオのTOPに戻る
　