ありとあらゆる音空間の建築設計を手掛けてきました。ソフトウェアの制作を行ってきました。工場、高速道路などの周辺の騒音伝搬を予測するソフト「ファクトリー」です。境界で幾何反射を繰り返すことによって音が伝搬するという考えに基づいています。回折や干渉などの現象を捕えられないという大きな欠点もあります。主な用途が初期反射音の伝搬経路や反射板のカバーエリアの検討などに限られているのが実状でした。これを解決する音場シミュレーションシステムの構築が緊急に求められていました。そこで当社では、近年のコンピュータ技術の発達に伴って注目度が高まっている波動音響理論に基づいた音場解析手法の基礎研究に着手しました。呼ばれるもので、境界要素法の支配方程式に波動方程式を適用して波動音場を解析するものです。境界条件の適用方法の改良による基礎方程式の適用範囲の拡大などの成果が得られたと考えています。今回は、1.についての概略と従来の幾何シミュレーション手法との対比に焦点を絞ってご紹介させていただきます。定常音場及び過渡音場解析の定式化をまとめた論文が数多く発表されています。空間内の唯一の点音源Sが角周波数ωの正弦波で定常に駆動する場合、受音点Pでの定常音場を次のように表すことができます。rはPQ間の距離、k(=ω/c)は波数を表します(cは音速)。境界Bの内側にある場合4π、B上にある場合2π、Bの外側にある場合0となる定数です。各記号については図1を参照してください。この式は、Greenの公式に含まれる2つの連続な関数が波動方程式を満たすと仮定して、音源や受音点などの極となる点を取り除くことにより導かれます。この式が表すように、空間内の任意の点の音場は、音源からの直接音と境界からの寄与を加算することによって求めることができます。(1)式を逆Fourier変換すると次のような過渡音場に関する積分方程式が得られます(tは時刻)。境界積分方程式で使用する各記号さて、これらの積分方程式を計算機上で数値的に解くには、積分の離散化や境界条件の適用などを検討する必要があります。(1)式を用いた定常音場の計算方法を簡単に説明します。(1)式にはΦ(P)、Φ(Q)、及び∂Φ(Q)/∂nQの三つの未知変数が含まれるため、このままでは解くことができません。そこで、まず受音点Pを境界上に置くことにより、境界上の音場に関する積分方程式を用意します。境界値問題の解法を用いることにより、∂Φ(Q)/∂nQをΦ(Q)の関数として表します。Φ(P)∈Φ(Q)、∂Φ(Q)/∂nQ=f[Φ(Q)]となり、式中の未知変数をΦ(Q)の一つのみに変換することができます。このような積分方程式は第二種Fredholm型積分方程式と呼ばれ、一般的な離散化手法によって解くことが可能です。そこで、対象とする周波数に応じた大きさの面積要素に境界を分割して積分の離散化を行い、各要素上では速度ポテンシャルが一定であるものと仮定します。N元の連立一次方程式を構成することができます。得られた値を受音点Pが空間内にある場合の積分方程式に代入すれば、所望の音場解析が完了することになります。過渡音場解析については、大きく分けて次の二つのアプローチがあります。一つは(1)式による定常音場解析を一定の周波数間隔を決めて各周波数毎に行い、その結果を逆Fourier変換する方法です(以下ではIFFT法と呼びます)。一つは(2)式を用いて時間応答を直接計算する方法です(以下では直接法と呼びます)。その結果計算量が多くなるという問題点があります。周波数応答を逆Fourier変換する際に回り込み誤差が生じる危険性もあります。少ない計算量で必要な解析を行えるという長所があります。先に述べた寺井先生や河井先生らが三角波応答を計算する手法を定式化しています。これは、サンプリング周期を指定して各サンプル毎の応答を過去から未来へ逐次計算するものです。計算アルゴリズムがかなり複雑になっています。法線方向微分型の積分方程式を用いることにより解の安定性を向上させることができるとしています。変動の少ない信号の使用や時間軸上での平滑化は、高い周波数ほど計算精度が低下するという欠点を生じると考えられます。法の計算フロー(2)式を用いた過渡音場解析の手順は、先に述べた定常音場の解法とほぼ同様です。最後の∂φ(Q)/∂tについても微分の定義式から近似的にφ(Q)の関数で表すことができます。後は同様に境界の離散化を行い、さらに時間軸の離散化を行うことにより、数値解析可能な境界上の音場を求めるための方程式を導出します。ところで、(2)式の右辺に含まれるすべての未知変数は、左辺のφ(P)に対してr/cだけ過去のものです。このことは、ある時刻からシーケンシャルに数値計算を行う際には、右辺が既知であることを示しています。したがって、(1)式で使用した連立一次方程式を解くステップが不必要となり、記憶容量の少ない計算機でも計算できることになります。双方に長所と短所があり、目的に応じてこれらの方法を使い分ける必要があると思われます。主な設計指標であるC値やD値などの定義に見られるように初期の応答が重要になります。早急な実用化を考えた場合、汎用の小型コンピュータで計算できると好都合です。これらの条件を考慮すると、直接法の採用が現実的であると考えられます。法は前項で紹介した河井先生らの手法が発表されている程度であり、未だに研究段階であると言えます。この手法にはいくつかの問題が残されています。従来の手法よりも簡単で、しかも任意の離散信号を音源に入力できるなどの自由度の高い手法を構築することができました。以下では、解析手法の発案に至った経緯とその解析手順のアウトラインを簡単に説明させていただきます。解析アルゴリズムを検討する際に最も注意したのは、(2)式の右辺に含まれる各関数のr/cの遅延計算をいかに精度良く行うかという点です。実際の計算機を用いた数値計算では、サンプリング定理に基づいた離散信号を扱うことになります。それ以外の場合にはサンプル間に到来する波を見逃さないためにサンプル間の信号を精度良く補間する必要があります。そこで、この問題に対処するために、ディジタル信号処理の概念を積極的に取り入れることにしました。連続信号と離散信号は完全に一対一に対応します。任意の時刻の正確な応答を得ることができる、即ち任意時間の遅延を正確に計算することができます。この関係を利用するには、離散信号の全体が既知であり、有限の長さを持っているという条件が必要になります。そこで、数値計算手順の第一段階である境界上の音場を計算する過程において、(2)式を次のように分割して計算することを考えました。音源に有限の長さの離散信号を入力し、この信号が音源から放射されて境界上の各点に到達するまでの過程を次の(3)式により計算します。(3)ここで、φ[0](P)は境界上の点Pでの直接音による応答を表します。反射音の計算ステップIの計算完了後は、音源はすべての信号を放射してしまったので、無視して良いことになります。次のステップではこれらの境界上の信号が音場内に再放射されて他の境界上の点へ到達するまでの過程を次の(4)式により計算します。(4)ここで、mは反射次数と対応し、ステップIの直後ではm=1となります。したがって、φ[m](P)は境界上の点Pでのm次反射音による応答を表します。これ以降は、mを更新して(4)式による計算を繰り返し行うことにより、各反射次数の境界上の点での応答が求められます。必要な反射次数までステップIIを繰り返し、最後に(5)式で直接音と各次の反射音を加算すれば、境界上の点での応答が得られます。音源からの直接音と境界からの寄与を加算すれば音場解析が完了します。本手法の計算フローこの解析手法の特長をまとめると次のようになります。音源信号にはサンプリング定理に関するもの以外の制約を受けることがなく、任意の有限長の離散信号を入力可能(自由な周波数選択性)。鏡像法や音線法などのように反射次数を単位とした解析結果を導出するため、これらの幾何的な手法との比較が容易。適当な次数で計算を打切り、初期応答だけを計算することが可能(自由な時間選択性)。サンプル間に到来する波についても忠実な推定が可能(各周波数を通して均一な計算精度)。各計算ステップをFFTを利用して周波数領域で行うことが可能(畳み込み演算、遅延計算の高速化)。計算手順のイメージ最後の特長は特筆すべき点であり、ディジタル信号処理の考え方を取り入れたことによる恩恵とも言えるべき長所です。メモリの少ない計算機でも、十分に実用に耐えるプログラムを作成できます。この解析アルゴリズムは並列処理にも適しており、使用する計算機やプログラミング次第で計算の高速化を計ることも可能です。興味のある方はそちらを参照してください。自由空間に浮かぶ剛な矩形板に無指向性点音源から球面波が入射するときの音波の反射性状を対象とした数値計算を行い、実測値との比較を試みました。実験ではこのような音場を無響室内に2枚の矩形板を配置し、音源に12面体スピーカを用いて模擬しました。実験のブロックダイアグラムは図6に示すとおりです。M系列-高速アダマール変換法によりインパルス応答を測定しました。同期加算の回数は15回としています。80Hz付近で無響室の固有モードが顕著であったため、音源への入力信号を100~1000Hzのバンドパスフィルタにより帯域制限しました。数値計算では、音源入力信号として実験に用いたスピーカの特性を測定したものを使用しました。計算でのサンプリング周波数は4000Hzとし、矩形板の表面を40×40個の正方形要素に分割して行いました。図7に示す結果から、各点の実測と計算の時間応答が良く一致していることがわかります。得られた時間応答をFourier変換して周波数特性を求めて見ましたが、その結果からも大きな差は確認できません。残響室内の矩形パルス応答(鏡像法との比較)実験Iでは、平らな矩形板の表面だけを解析対象としましたので、数値計算では一次反射波のみを計算しています。そのため、反射次数の増加に伴う演算誤差の蓄積の度合いや、解の安定性を確認するには至っていません。そこで次のステップとして、剛壁で囲まれた残響室を対象として数値計算を行い、実測値との比較を行いました。音の波動性を考慮した音場解析の有効性を検証しました。使用した測定システムのブロックダイアグラムは図8のようになります。残響室内に設置した12面体スピーカを矩形パルスにより駆動し、サンプリング周波数1000Hzで過渡応答を測定しました。100~300Hzのバンドパスフィルタで帯域制限しています。12面体スピーカの矩形パルス応答を無響室にて測定したものを使用しました。境界は3546個の四角形要素に分割し、境界上の応答を八次反射まで計算しました。これに音源入力信号を畳み込んで矩形パルス応答を算出しました。受音点の位置は図8に示す三つです。図9に示す結果を見ると、本手法の計算結果は実測値に良く一致しています。それ以降に誤差が生じているのが確認できます。対象とした残響室の形状が不整形であるために生じる音線経路の影による誤差が原因として考えられます。反射次数毎に誤差の要因を具体的に考察すると、以下のようになります。音源の位置によってはその淵の部分に鏡面反射点が隠れてしまう場合があります。エッジからの回折波だけが受音点に寄与するため、波動性を考慮しない鏡像法では計算することができません。これによる回折波についても鏡像法では計算できません。生じる回折波も増加していきますが、さらに以前の次数で生じた誤差も反射を重ねる毎に蓄積されていきます。対象とする空間の形状が複雑になり、反射面の大きさが波長に比べて小さくなるに伴って誤差を生じる危険性が高まることになります。これからの音場解析では音の波動性を考慮することが主流となってくるでしょう。本手法の実用化に際して大きな問題が残されています。それは、壁や天井などに取付ける建築材料の吸音特性のデータベースをどのように収集するかです。過渡音場の数値解析では、吸音材の時間応答特性が必要になります。吸音率の1/1あるいは1/3オクターブバンドの周波数特性がほとんどです。このようなデータには反射の際の位相遅れに関する情報は含まれていませんので、時間応答を算出することはほとんど不可能であると言えます。したがって、吸音材の時間応答特性に関するデータベースを新規に構築する必要があり、またそのための測定方法についても検討しなければなりません。良いアイデアが浮かばず、少々停滞気味になっています。ご教授していただければ幸いです。設計施工の効率化を図ると同時に、さらに良い音響空間を実現できるようになります。アクティブノイズコントロールのシミュレーションなどで早急に対象となる音場の特性を把握されたい方に、音場シミュレーションを実施します。ソフトウェアの販売お客様のご要望に応じてプログラムをカスタマイズし、販売させていただきます。市販の境界要素法ソフトウェアをお持ちの方には、そのプログラムとのリンクも検討いたします。これ以外にも色々な御相談に応じます。さらなる検討を進めていく所存です。今後に御期待ください。本研究に関しては、東京電機大学の浜田晴夫助教授の御指導を頂きました。実験に関しては、(株)奥村組筑波研究所の飛松様、木村様の御協力で残響室及び測定システムを貸与して頂きました。以上の方々をはじめ関係者各位に心から感謝いたします。