(57)【要約】【課題】仮想空間内の音響をリアルタイムで再現すること【解決手段】リアルタイム音響再現システムは、室内音響を再現する過程の中で大きく３段階に分けることができる。部屋の形状，環境，条件，位置など様々な情報を収集する情報入力処理、音響伝達関数を用いた音響モデリングを行う変換処理、実際の音声として再生する出力処理（可聴化処理）の３段階である。  音響モデリングを行う変換処理は、部屋２１０の形状に基づいて生成された反射パス２２６，２２７，２２８に対して、反射における吸音や減衰などの損失を含むインパルス応答２３１を算出し、これを伝達関数２３６として用い、入力音声２３８を変換することで音響が再現される。【選択図】図２
2005080124.tif 000002
【特許請求の範囲】
【請求項１】
仮想空間内におけるリアルタイム音響再現システムであって、
少なくとも仮想空間の形状、壁面の材質、温度，湿度の基本情報が格納されている記憶手段と、
前記記憶手段から読み出した基本情報により、仮想空間内の音源位置から動的に与えられた受音点位置までのインパルス応答をフーリエ変換した伝達関数を生成する伝達関数生成手段と、
動的に与えられたデジタル音声データを一定時間切り出してフーリエ変換し、フーリエ変換された周波数帯のデータに対して前記伝達関数を適用し、逆フーリエ変換して、仮想空間内の音声再現データとして出力する可聴化手段と
を備えることを特徴とするリアルタイム音響再現システム。
【請求項２】
請求項１に記載のリアルタイム音響再現システムにおいて、
前記伝達関数生成手段における前記伝達関数生成は、前記仮想空間内の音源位置から受音点位置までのパスを生成し、それぞれのパスごとの室内伝達関数、頭部伝達関数を求めることで行われることを特徴とするリアルタイム音響再現システム。
【請求項３】
請求項２に記載のリアルタイム音響再現システムにおいて、
前記伝達関数生成手段における前記室内伝達関数は、統計学的な推測により求めた高次反射残響の推測値を用いて求められることを特徴とするリアルタイム音響再現システム。
【請求項４】
請求項２又は３に記載のリアルタイム音響再現システムにおいて、
前記音源位置は、基本情報として与えられるとともに、
前記仮想空間をグリッドに分割して、各グリッドに受音点があるとして、グリッドごとの伝達関数を予め計算して伝達関数記憶手段に格納し、
前記伝達関数生成手段は、前記受音点が動的に与えられると、該受音点位置のグリッドに対する伝達関数を前記伝達関数記憶手段から読み出す
ことを特徴とするリアルタイム音響再現システム。
【請求項５】
請求項４に記載のリアルタイム音響再現システムにおいて、
前記グリッドを反射次数ごとに大きさを変えて、反射次数ごとの伝達関数を計算して前記伝達関数記憶手段に格納し、
前記伝達関数生成手段は、受音点が動的に与えられると、該受音点位置の反射次数ごとのグリッドに対する伝達関数を前記伝達関数記憶手段から読み出して合成する
ことを特徴とするリアルタイム音響再現システム。
【請求項６】
請求項４に記載のリアルタイム音響再現システムにおいて、
前記音源位置は、基本情報として与えられるとともに、
前記仮想空間をグリッドに分割し、各グリッドに到達するパスにＩＤをつけてグリッドごとにグリッドパス記憶手段に格納し、
前記伝達関数生成手段は、受音点位置が動的に与えられると、受音点位置のグリッド中のパスを前記グリッドパス記憶手段から読み出し、読み出したパスを用いて、受音点のパスを得る
ことを特徴とするリアルタイム音響再現システム。
【請求項７】
請求項４〜６のいずれかに記載のリアルタイム音響再現システムにおいて、
前記伝達関数生成手段における前記グリッドは、前記インパルス応答の音圧レベルの変化に応じて、グリッドの分割サイズを変化させることを特徴とするリアルタイム音響再現システム。
【請求項８】
請求項４〜６のいずれかに記載のリアルタイム音響再現システムにおいて、
前記伝達関数生成手段における前記グリッドは、前記インパルス応答のベクトル内積に応じて分割のサイズを変化させることを特徴とするリアルタイム音響再現システム。
【請求項９】
請求項１〜８のいずれかに記載のリアルタイム音響再現システムにおいて、
前記可聴化手段は、複数のバッファを有しており、
入力音声を複数のブロックに分けて前記複数のバッファに入力してフーリエ変換し、フーリエ変換された周波数帯のデータに対して、バッファごとに前記伝達関数を適用して逆フーリエ変換し、複数のバッファのデータを足し合わせて仮想空間内の音声再現データとして出力することを特徴とするリアルタイム音響再現システム。
【請求項１０】
請求項１〜９のいずれかに記載のリアルタイム音響再現システムをコンピュータ・システムに構築させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、仮想空間におけるリアルタイム音響再現システムに関するものである。
【背景技術】
【０００２】
現在、計算機上の仮想空間(Virtual Space) を実現するために、いかにリアルに見せることができるのかといった視覚的な面を重視する研究が多い。鏡面反射の計算，ライティング，陰影表示，テクスチャーおよび，バンプマッピングなどの技術は現実に近い空間をつくることを可能にする。しかし、人間に伝わる情報は様々であり、視覚的な情報だけでリアルな仮想空間を作ることができるとは限らない。
例えば、レンガで作られた部屋を木造に変更した場合、視覚的にそれらしいように感じるのだが、音響に変化がなければ当然我々は違和感を覚える。また、両耳に到達する音の時間的および音量的な差により我々は音の到達方向を感覚的に認知でき、音源が後方にあった場合、その位置を予測できる。
我々は、５感のうち、大きな割合で視覚に頼っている。そのため、周りの空間を認知するには視覚的な情報だけで十分であると考えがちである。しかし、全ての感覚機能は同じ神経を伝わり、最終的には唯一の情報として認知される。すなわち、これらの情報が影響し合うことで、より正確に空間を認知できる。
【発明の開示】
【発明が解決しようとする課題】
【０００３】
本発明では、グラフィックスによる情報伝達の支援として、音響をリアルタイムで再現することにより、グラフィックスだけでは表せない空間的な臨場感をユーザに伝達することを目的とする。すなわち、現実に近い音響よりも、仮想空間上での変化を対話的に再現することを重視することで、グラフィックスとの連携を図っている。
特に音響計算は非常に時間がかかり、リアルタイムで実現するのは非常に難しいとされている。そのため、本発明では、特に仮想空間における音響計算の高速化も目的としている。
【課題を解決するための手段】
【０００４】
上記目的を達成するために、本発明は、仮想空間内におけるリアルタイム音響再現システムであって、少なくとも仮想空間の形状、壁面の材質、温度，湿度の基本情報が格納されている記憶手段と、前記記憶手段から読み出した基本情報により、仮想空間内の音源位置から動的に与えられた受音点位置までのインパルス応答をフーリエ変換した伝達関数を生成する伝達関数生成手段と、動的に与えられたデジタル音声データを一定時間切り出してフーリエ変換し、フーリエ変換された周波数帯のデータに対して前記伝達関数を適用し、逆フーリエ変換して、仮想空間内の音声再現データとして出力する可聴化手段とを備えることを特徴とする。
前記伝達関数生成手段における前記伝達関数生成は、前記仮想空間内の音源位置から受音点位置までのパスを生成し、それぞれのパスごとの室内伝達関数、頭部伝達関数を求めることで行ってもよい。
前記室内伝達関数は、統計学的な推測により求めた高次反射残響の推測値を用いて求めてもよい。
【０００５】
前記音源位置は、基本情報として与えられるとともに、前記仮想空間をグリッドに分割して、各グリッドに受音点があるとして、グリッドごとの伝達関数を予め計算して伝達関数記憶手段に格納し、受音点が動的に与えられると、該受音点位置のグリッドに対する伝達関数を前記伝達関数記憶手段から読み出し、該伝達関数を用いて音声再現データを出力してもよい。
前記グリッドを反射次数ごとに大きさを変えて、反射次数ごとの伝達関数を計算して前記伝達関数記憶手段に格納し、受音点が動的に与えられると、該受音点位置の反射次数ごとのグリッドに対する伝達関数を前記伝達関数記憶手段から読み出して合成し、合成した伝達関数を用いて音声再現データを出力してもよい。
【０００６】
前記音源位置は、基本情報として与えられるとともに、前記仮想空間をグリッドに分割し、各グリッドに到達するパスにＩＤをつけてグリッドごとにグリッドパス記憶手段に格納し、受音点位置が動的に与えられると、受音点位置のグリッド中のパスを前記グリッドパス記憶手段から読み出し、読み出したパスを用いて、受音点のパスを得てもよい。
前記グリッドは、前記インパルス応答の音圧レベルの変化に応じて、グリッドの分割サイズを変化させてもよく、又、前記インパルス応答のベクトル内積に応じて分割のサイズを変化させてもよい。
【０００７】
前記可聴化手段は、複数のバッファを有しており、入力音声を複数のブロックに分けて前記複数のバッファに入力してフーリエ変換し、フーリエ変換された周波数帯のデータに対して、バッファごとに前記伝達関数を適用して逆フーリエ変換し、複数のバッファのデータを足し合わせて仮想空間内の音声再現データとして出力してもよい。
上記記載のリアルタイム音響再現システムをコンピュータ・システムに構築させるプログラムも本発明である。
【発明の効果】
【０００８】
本発明では、上記構成により、音響をリアルタイムで再現することにより、グラフィックスだけでは表せない空間的な臨場感をユーザに伝達することができる。
本発明では、音響計算の高速化を図ることにより、リアルタイム性を確保できている。
【発明を実施するための形態】
【０００９】
本発明におけるリアルタイム音響再現システムは、図１及び図２に示す概要のように、室内音響を再現する過程の中で大きく３段階に分けることができる。図１において、部屋の形状，環境，条件，位置など様々な情報を収集する情報入力処理（Ｓ１０２）、音響伝達関数を用いた音響モデリングを行う変換処理（Ｓ１０４）、実際の音声として再生する出力処理（可聴化処理）（Ｓ１０６）の３段階である。
【００１０】
基本的な音響再現システムの概要は、図２に示すように、部屋２１０の形状に基づいて生成された反射パス２２６，２２７，２２８に対して、反射における吸音や減衰などの損失を含むインパルス応答２３１を算出し、これを伝達関数２３６として用い、入力音声２３８を変換することで音響が再現される。
【００１１】
（１）条件・部屋の形状情報の入力（Ｓ１０２）
最初にシミュレーションする部屋２１０の形状を入力する。この段階で３次元ポリゴンモデルおよび各反射面（周囲の壁）の材質（吸音特性）２３２，２３４を定義する。これらは音響再現に必要な基本情報として扱われる。各面の角度によって音が到達する方向が決まり、吸音率によって空間の響き、つまり残響が決まる。
上記の部屋の形状や材質などの、ほぼ固定された情報に対し、音源や受音点の位置または入力音声信号などの動的な情報がある。これらを元に音響を再現する伝達関数２３６を生成する。
ただし、以下の説明では、音源は固定されたものとして述べる。
【００１２】
（２）音源Ｓから受音点Ｒへの伝達関数の生成（Ｓ１０４）
構築された部屋の形状の中で音がどのように反射し、受音点Ｒまで到達するかを幾何音響解析による算出によって求める。計算は、Path Tracing（音線法）を用いて、音源Ｓ２２２から放出される音を線（パス）として扱い、任意の受音点Ｒ２２４に到達する複数の反射音（反射パス）２２６，２２７，２２８を求める。計算結果として、音源Ｓ２２２から受音点Ｒ２２４までの伝搬経路を線で表し、表示したものは音線図という。伝達関数計算に対する条件には、反射する回数（反射次数）や音源の指向性、受音点の大きさ、室温や湿度を設定する。
【００１３】
音線の伝搬経路を計算することで、任意の受聴位置で最初に直接音が到達した後、設定された吸音率の反射面に反射した１次反射音，２次反射音，…，ｎ次反射音が時間遅れで到達する様子を、時系列上で配列するエコーダイアグラム２３１として求めることができる。このエコーダイアグラム２３１に空間の様々な特性を持たせれば、インパルス応答波形が得られる。インパルス応答をフーリエ変換すれば音源Ｓ・受音点Ｒ間の伝達関数が得られ、これが音源Ｓから受音点Ｒへの伝送周波数特性を表す。
エコーダイアグラム２３１は、図２に示すように、音の遅れと大きさの情報だけをもつ波形であるのに対し、インパルス応答は、エネルギー的に平坦な周波数特性をもつ短音（ インパルス）に対する実際の応答であるので、空間の全ての音響情報が含まれると考えられる。
【００１４】
（３）可聴化処理（Ｓ１０６）
インパルス応答に無響室で録音されたデジタルデータである音源（ドライソース）２３８を畳み込む（２３７）ことにより、その空間の響きを持つ音２４２として変換できる。つまり、あたかもその室内で聴いたような感覚が得られる。しかし、インパルス応答の算出および畳み込みを、ユーザの操作に対して対話的に行うことは非常に困難である。本発明では、リアルタイムおよび対話性に適応するため、時間帯域での畳み込みではなく、周波数帯域で処理を行う。この理由は後で詳しく述べる。
本発明の音響再現システムでは、視覚的なグラフィックス表示とリアルタイムで組み合わせることにより、ユーザにより多くの情報を与え、さらに、これらの情報の協調により仮想空間にさらなる現実感を与えることを目指している。以下に、それぞれの処理について詳しく説明する。
【００１５】
＜空間形状の入力（基本・動的情報の入力）＞
本発明の実施形態における音響再現システムでは、仮想空間の環境情報，位置，条件，および音声を入力パラメータとし、音響効果をもつ音を出力している。入力パラメータには２つの種類があり、図３に示すように、あらかじめ設定する必要のある基本情報パラメータ３１０とユーザの操作時にリアルタイムで更新される動的情報パラメータ３３０に分けられる。
基本情報パラメータ３１０は、部屋の３次元形状情報、これに含まれる各面の材質（吸音特性）や、温度および湿度を含む環境情報から構成される。これらの情報は、基本情報データベース３２０に格納され、必要なときに読み出され、これらの情報を基に、これから述べる空間分割やパス情報データベース生成などの高速化のための前処理を行う。音源の位置は本来動的な情報であるが、以下の説明では、音源の位置を固定されたものとして述べる。すなわち、音源位置も基本情報パラメータとして扱っている。
【００１６】
動的情報パラメータ３３０は、受音点の位置や入力音声（音源から発するドライソースのデジタル音声信号）からなる。仮想空間での位置を受音点として、パス情報データベースで限定された反射経路のパスをトレースし、パス総距離やパス毎の減衰率などの情報をリアルタイムで生成・更新する。このデータを基に生成される伝達関数（後で詳しく述べる）を用いて、音源から再生される音声の変換を行う。なお、前にも述べたように、動的パラメータとして、音源位置も含めることができる。
この受音点データは、例えばこのシステムをゲームに適用した場合、仮想空間で行うゲームの進行に従って、ゲーム参加者が操作しているキャラクタの仮想空間内の位置として与えられる。また、音源からの音声データは、ゲームの進行により、その場で発生する音として与えられる。また、音源の位置も動的情報として与えることもできる。
以上のデータの流れは、リアルタイム再現において重要な意味をもたらす。
【００１７】
＜音響伝達関数の生成＞
システムの入力処理の段階を経て、入力された情報を元に、音源から受音点への音の伝達関数を求め、音響のモデリングを行い、最終段階である可聴化処理で実際の音として出力する。ここでは、この伝達関数の生成について述べる。
音源を出るドライソースに対する音源位置から受音点までの音響的な変化は、インパルス応答として表すことができる。これをフーリエ変換すれば周波数帯域上の伝達関数になる。これを利用して入力音声に対する音響を再現する。インパルス応答を生成する手法は様々であり、主に次のような手法が用いられる。
最も精密な方法は、測定によるインパルス応答の算出である。実際に測定を行う部屋でインパルスを発し、その応答を解析する。測定方法は、入力インパルス短音の形や長さ、応答の周波数別の測定とそのときに用いる窓関数の種類などとさまざまである。測定を用いる手法は最も正確に算出できる反面、測定を行う部屋が存在する必要がある。これに伴う、測定の手間と費用が膨大である。
【００１８】
測定によるインパルス応答の計算にかかる手間と費用を削減するため、音響学では計算機を用いたインパルス応答の算出法が多く研究されている。これらの手法は、大きく２種類に分けることができる。波動音響と幾何音響である。
波動音響は、差分法や有限要素法などを用いて、格子に分割された空間での音波の伝播を波動現象として捉える分野であり、測定されたインパルス応答の精度には及ばないが、詳細に求めることができる。しかし、計算する周波数は格子の大きさに依存するため、低周波数に対しては有効であるが、高周波数では、計算量が膨大になる。
幾何音響は、幾何的な計算によるインパルス応答の算出を扱う分野であり、虚像法(Image Source Method)や音線法(Ray Tracing Method)などが用いられる。これらは、上述の２つの方法に比べて圧倒的に速い計算速度をもつ。ゆえに、リアルタイム音響再現において重要な要素であるといえる。まして空間が視覚的に表現されている場合、特に仮想空間の表現に関しては、計算精度よりも計算速度（リアルタイム性）を重視する方が空間情報をユーザに伝えやすい。そのため、本システムでは幾何音響を用いてインパルス応答の計算を行う。
【００１９】
音源から受音点への音響伝達関数(ATF)を求める過程を図４に示す。図４に示されるように、音響伝達関数(ATF)は、パス生成Ｓ４０２，損失特性の算出Ｓ４０４，室内伝達関数(RTF)Ｓ４０８，頭部伝達関数(HRTF)Ｓ４０６の算出などいくつかの段階に分けて行われる。
【００２０】
（パス生成（Ｓ４０２））
音源から周囲の壁に反射して受音点に到達する反射音は、波の物理的な現象によって作られ、図５に示すように、音源から円形に広がる。反射物（図５の壁Ａ）で、波は折り返し、反転した状態でさらに広がる。反射した波形が受音点と接する点が描く音源Ｓ−受音点Ｒ間の軌跡は、幾何音響で用いられる反射パスに相当する。
壁に反射したときの円形をした音波の中心点は、音源点Ｓの壁Ａに対する線対称点（３次元空間の場合は面対称点）ＳＡになる。すなわち、壁Ａの反射音波は、ＳＡの地点から広がる音波に等しい（ 壁Ａでの吸音を除く）。したがって、すべての反射パターンにおける虚音源を求めることによって、すべての反射音を算出することができる。この方法は、幾何音響学でいう虚像法（Image Source Method）であり、ＳＡは虚音源(Image Source) と呼ばれる。
虚像法は、虚音源をあらかじめ求めるが、反射を除く壁との衝突計算が行われないため、複雑な形をした部屋には不向きである。したがって、本システムでは音線法を用いる。
図６は、音源Ｓから壁Ａ，Ｂ，Ｃの順に反射し、受音点Ｒに到達する反射パスの算出方法を示す。音源Ｓの壁Ａに対する虚音源ＳＡを求め、次に虚音源ＳＡの壁Ｂに対する虚音源ＳＡＢを求める。同様に虚音源ＳＡＢＣを求め、受音点Ｒを結ぶ線分ＳＡＢＣＲと壁Ｃとの交点ＰＣが第３反射点になる。同様にＳＡＢＰＣと壁Ｂとの交点である第２反射点ＰＢを求め、さらに、第１反射点ＰＡを同様に求める。音源Ｓとこの３つの点ＰＡ，ＰＢ，ＰＣを介し受音点Ｒを結ぶパスが、反射パスＰＡＴＨＡＢＣである。
図７（ａ）〜（ｆ）に、直接音と、１〜５次反射パスを計算した例を示す。各図には、それぞれに対応するエコーダイヤグラムを添えている。  
【００２１】
（大気伝播及び反射吸音の損失特性（Ｓ４０４））
音響の計算は、Path Tracing によって求められた反射パス別に特性を算出した後、それぞれの遅延時間に伴って推移させられ、最終的に足し合わされる。反射パスｋに代表される反射音の特性は、図８に示すいくつかの要素から求められる。
・音源の指向特性関数ＳＦθ（ω）（図８の（１））
・大気伝播による減衰特性Ａｋ（ω）（図８の（２），（４），（６），（８））
・拡散音場による減衰率Ｅｋ（図８の（２），（４），（６），（８））
・反射吸音による損失特性Ｒｋ（ω）（図８の（３），（５），（７））
・頭部伝達関数による特性ＨＦθ（ω）（図８の（９））
音源には指向特性が設定されており、正面からの角度によって周波数特性が異なる。すなわち、直接音に関して述べると、正面から聴いた音の方が、後方から聞いた音に比べて大きく、高周波成分が多く含まれる。音源を発する反射音（反射パス）は、音源の向きとの角度から、指向特性ＳＦθが決まる。
また、反射音は音源と受音点間の経路でさまざまな損失が生じる。音源から受音点への距離に応じて、空気伝播中に摩擦による減衰Ａｋおよび拡散音場による減衰Ｅｋが生じる。反射がある場合、反射壁面の材質に伴って吸音による損失ｍｒが生じる。音は反射を繰り返し、最終的に受音点に到達したとき、パスｋの反射による損失はＲｋ（ω）になる。ＳＦθ（ω），Ａｋ（ω），Ｒｋ（ω）をそれぞれフィルターとして扱い、音源Ｓ（ドライソース）に適用する。
【００２２】
（音源の指向性による減少）
すべての方向に対して同様な特性をもつ無指向性音源（点音源ともいう）に対して、指向性を持つ音源は方向に伴い周波数特性が異なる。実際、理想的な点音源は非現実的であり、実世界ではあり得ない存在である。すなわち、すべての音源には指向性があり、空間の状態を認知するには重要な要素である。
本システムで音源を設置する際、ユーザが音源の指向特性を設定できる。音源の指向性は、本来、発音体の形に依存する特性を持ち、すべての方向へ発する音の周波数特性は異なるが、特殊なケースを除き、本システムの対象となる音源のほとんどは前方の方が音は大きく、後方が音は小さい。これによって、リスナーは音源の向きを認知する。したがって、図９に示すように、本システムではこの性質を利用し、音源の前方を表すベクトルに対するパスの放出角度からユーザが設定した特性を用いる。
【００２３】
図１０において、音源の指向方向のベクトルをｖ１とし、パストレーシングにより算出された反射パスが、音源をベクトルｖ２方向へ出るとしたとき、ベクトルｖ１，ｖ２の内積を求めることにより、ベクトル間の角度θは次のように算出できる。
【数１】
2005080124.tif 000003
この角度θ（０〜π）方向に相当する周波数特性（図１１参照）は、ユーザが設定でき、再現時に生成された各反射パスの方向から、相当する周波数特性のフィルターＳＦθ（ω）が決まる。
【００２４】
（反射面の吸音特性）
反射パスｋは伝播中に反射吸音による損失Ｒｋ（ω）が含まれる。音源から放出されるすべての反射音は、異なる反射経路を辿って受音点へ到達する。音が壁面に反射する際、図１２に示すように、その一部が反射面に吸収され、壁の材質（および入射角）によってその吸音特性が異なる。図１２（ｂ）は、図１２（ａ）を拡大した模式図である。反射パスが最終的に受音点に到達したときの反射吸音による減少率は、反射回数をｎとして、
【数２】
2005080124.tif 000004
になる（ただし、ｍｒｋｒｅは、ｎ次反射パスｋのｒ反射目での吸音特性）。これをデシベル単位に直すと、損失デシベル数ＭＲｋｒは、
【数３】
2005080124.tif 000005
となる。したがって、
【数４】
2005080124.tif 000006
となる。
【００２５】
（大気伝播による損失特性）
反射吸音による損失以外にも、大気伝播による減衰がある。大気伝播による減衰ｍ（Ｔ，Ｈ，ω）［ｄＢ／ｍ］は、空気に含まれる成分（窒素や酸素など）によって異なるが、成分を固定し、気温Ｔと湿度Ｈから吸音率を求めることができる。本システムでは、米国標準規格ＡＮＳＩ  Ｓ１，２６−１９７８で定義されている計算法を用いて、音が１ｍ進むときの減衰率ｍ（ω，Ｔ，Ｈ）を求め、パスの総距離ｌｋから損失量Ａｋ（ω）を次のように求める。
【数５】
2005080124.tif 000007
同様にデシベル単位で表すと
【数６】
2005080124.tif 000008
になる。
【００２６】
＜音響伝播関数＞
音響伝達関数(ATF: Acoustic Transfer Function) とは、室内における発音体から受聴者の鼓膜に到達する音の伝達関数である。すなわち、音が周囲の壁の反射及び大気伝播による減衰，または、頭部周辺での回折，透過損失，反射など、両耳の鼓膜に到達するまでの伝達特性を指す。本リアルタイム音響再現システムの生成論文では、音響伝達関数を、部屋の形状や環境による特性（反射，減衰など）と受聴者の頭部形状による特性に分け、それぞれ、室内伝達関数(RTF)と頭部伝達関数(HRTF)と呼ぶ。
【００２７】
（室内伝達関数の算出）
室内伝達関数（RTF：Room Transfer Function）は、室内における音源から受音点まで、壁の反射や大気減衰などを含む伝達特性を表す。これを算出するには、反射パスの遅延時間と損失特性の２つの要素が必要になる。損失特性に関しては、前に述べた大気伝播による損失特性および反射面の吸音特性からなる。遅延時間は、反射パスの総距離から音速ｃを用いて算出できる。また、室内伝達関数は１つの反射音だけではなく、複数の反射音からなる。すなわち、すべてのパスの合成関数を求めることになる。
吸音や減衰などの周波数特性を考慮する場合、インパルス応答を算出する際に損失特性情報を含める必要がある。このとき、遅延時間ｄは時間帯域にあるのに対し、損失特性Ｔｌは周波数帯域の関数であるので、一般に、帯域変換を行わなければならない。帯域変換には、フーリエ変換および逆変換などのアルゴリズムを要するため、計算量が多く、無数の反射パスに対して変換を行うと、膨大な時間を要する。その結果、リアルタイム性が失われる。
【００２８】
本システムでは、壁の材質が持つ吸音特性などを考慮した音響再現を可能にするため、周波数帯域で再現を行なう。この手法について以下に述べる。
次に示すフーリエ変換が持つ性質を利用し、フーリエ変換回数を削減することによって、計算量の問題点を解決し、周波数特性を考慮したリアルタイム再現を実現させることができる。
時間帯域関数ｆ（ｔ）のフーリエ変換Ｆ｛ｆ（ｔ）｝をＦ（ω）としたとき、Ｆ（ω）は周波数帯域の関数となる。このとき次の性質が成り立つ。
【数７】
2005080124.tif 000009
これらの性質より、周波数帯域に変換された波形Ｆ（ω）から、
ａ．ｘ秒の時間的な遅延を持つ波形Ｆ（ω）・ｅｉωｘ
ｂ．遅延時間の異なる２つ以上の波形の合成波Ｆ１（ω）＋Ｆ１（ω）＋・・・＋Ｆ１（ω）
を求めることができる。
【００２９】
また、総損失特性αωによる損失を適用した波形はαω・Ｆ（ω）となる。
したがって、音源から、直接音と反射音をまとめて周波数帯域上で計算することが可能となり、フーリエ変換のプロセスは順・逆変換を一回実行するだけで反響（全ての反射音）を再現することができる。音源波形を標本化し、入力信号列をｆ（ｔ）とするとき、そのフーリエ変換Ｆ（ω）は次のようになる。
【数８】
2005080124.tif 000010
求める反響における反射音の数がｎ個，ｋ番目の反射音の遅れをｘｋ，総吸音特性をαωｋとしたとき、入力信号列に対する反響Ｆ（ω）は
【数９】
2005080124.tif 000011
となる。これをまとめると
【数１０】
2005080124.tif 000012
となる。
【００３０】
入力信号列を除く項をＲＴＦ（ω）とすると、ＲＴＦ（ω）は、入力信号列に依存しない係数になる。すなわち、任意の信号列のフーリエ変換とこの係数の積を求めることによって、時間帯域に戻したとき音響が再現される。
【数１１】
2005080124.tif 000013
上記の式ではＲＴＦ（ω）は音響再現に用いられる音響伝達関数となる。音響伝達関数は入力信号に依存しないことから、あらかじめ求めることができ、再現時の負担を軽減する。つまり、リアルタイム音響再現の実現に大きく寄与する。
【００３１】
一般に用いられるフィードバック回路を用いたインパルス応答フィルタ（ＦＩＲフィルタ）の場合、時間帯域上の計算であり、ｎ反射音再現にはｎ個の加算器，乗算器，遅延器が必要となる。よって、再現処理時間は反射音の数に依存する。
これに対して、上式の音響伝達関数ＲＴＦ（ω）を用いた場合、加算，乗算，遅延とも周波数帯域で一つのマトリックスで計算されるため、反射音の数は音響伝達関数生成時のみに関係する。再現時の処理速度は反射音の数に依存しない。
【００３２】
（頭部伝達関数の算出）
聴覚における空間情報の把握は、主に両耳受聴によって行われており、特に音の方向性の知覚は両耳の音圧レベル差，時間差によると言われている。このような音像定位を再現するには、一般に頭部伝達関数が用いられる。
頭部伝達関数(Head-Related Transfer Function：HRTF)は、反射波が全くない空間（自由空間）における、頭部や耳たぶの影響を含め、音源から聴取者の鼓膜までの音響伝達関数であり、両耳の音圧レベル差，時間差の情報を含んでいる。したがって、反射物を含む室内環境においては、虚像法の原理からそれぞれの反射音（虚音源）は独立した音像を持っていると考えられる。このことから、ＨＲＴＦをＲＴＦと分離して扱うことができる。
測定された両耳のＨＲＴＦを音声信号に畳み込むことによって、バイノーラルな信号が得られる。すなわち、立体的な音場を再現することができ、結果的に受聴者は音の方向性を知覚することが可能になる。なお、バイノーラル(binaural)とは、両耳で聴くことを意味するが、立体音響を指す。一般にバイノーラル方式といい、人間の耳の間隔にマイクをセットする録音方法で、３６０度の立体音場をヘッドフォンで再現する技術である。
ＨＲＴＦは人の頭部やダミーヘッドによって実測することができる。しかし、ＨＲＴＦは方向に依存した関数であるので、測定した音源方向に対しては忠実な立体音場を再現することができるが、その他の音源方向については再現に必要なＨＲＴＦを新たに測定する必要がある。実際には、全方位に対してＨＲＴＦを測定することは困難であり、実測値から補間を行う方法，球体・回転楕円体モデルによる合成による方法，などでＨＲＴＦを求める手法が行われている。
【００３３】
図１３に示すのは、水平位置に対して０〜３６０度に対する周波数特性である。受音点向きに対する音源の位置を指すベクトルの角度より、該当する周波数特性を音源に畳み込むことで、音像定位を再現することができ、ユーザに音源の位置を認知させることが可能である。実際、室内において、パス生成により求められた数多くの反射音すべてに対しＨＲＴＦ関数を適用することは、畳み込み演算量を考えると現実的ではない。そのため、通常、低次反射のみに対して計算する。
他の手法として、図１４に示すように、受音点から角度別にエリア分割し、エリア毎に反射音を合成してから処理を行うことで、分割数の数以下の畳み込み演算が保証される。しかし、エリア分割を用いた手法の場合は、当然のことながら、頭部伝達関数は離散化されるため、補間が大きな課題になる。これに対し、後で、頭部伝達関数の単純化手法を提案し、すべての反射音に対して音像定位を再現する手法について述べる。
【００３４】
（後部残響音の予測）
幾何的な音響の予測によって高次反射音を計算する際、膨大な計算時間を要することから、通常、減衰して音圧レベルの低い無数の高次反射音を後部残響音(Late Reverberation)として分けて扱うことが多い。
本システムでは、同様に初期反射(Early Reflections)について、幾何的手法を用いて反射音の遅れや特性を求め、後部残響音を統計学的な予測により静的モデルとして生成する。このときの分岐点を、初期反射と後部残響音が分離しないようにするために、時間ではなく反射回数を基準に分ける。
【００３５】
一般に残響音の予測手法としては、次のようなものが挙げられる。
・初期反射群をフィードバックさせることによる残響生成
・残響式から残響時間を算出し、これに基づいてwhite noise を指数関数と合成することによる残響生成
・1 回反射，2 回反射，…ｎＥ 回反射の初期反射音から統計学的な手法を用いて、ｎＥ＋１以降の反射を推定することによる残響生成
本システムでは、これらのうちの統計学的解法を用いた手法を利用する。
【００３６】
統計学的な分析を行う理由の一つは、エコーダイアグラム波形のこれから述べる特徴にある。図１５に示す部屋において、１〜５次反射まで計算したときのエコーダイアグラムを図１６に示す。
図１６（ａ）の縦軸を対数表示にしたグラフを図１６（ｂ）に示す。図１６（ｂ）のグラフを反射回数別に分けると図１７となる。図１７のグラフから、エコーダイアグラム波形は反射次数ごとに、規則的に分かれるという特徴がみられる。したがって、本システムにおける後部残響音予測モデルは、規則性を見出し、以下の統計学に基づく手法を利用する。
・エコーダイアグラムの回帰分析
・各反射次数パスの分散
・各反射次数パスの重心
【００３７】
（反射次数パス別の回帰分析）
本手法で用いるエコーダイアグラムの回帰分析は、ｒ次反射パスの回帰線を求めることで行う。音のインパルス応答の波形は、伝播中の損失によって、時間に伴い音圧レベルが低くなる傾向がある。すなわち、回帰線を求めることで平均的な減衰率を求めることができる。パスｋの遅れをＸｋ，音圧レベルをｌｏｇＹｋ  とし、回帰直線ｌｏｇＹ＝ａＸ＋ｂは次のように求められる。
まず、それぞれの変数ＸとｌｏｇＹの平均値バーＸとバーｌｏｇＹは、
【数１２】
2005080124.tif 000014
になる。次に、ＸとｌｏｇＹの共変動ΦＸｌｏｇＹ を求める。
【数１３】
2005080124.tif 000015
さらに、Ｚの変動ΦＸを求める。
【数１４】
2005080124.tif 000016
ここで、直線（対数表示）ｌｏｇＹ＝ａＸ＋ｂの各定数ａ，ｂは、
【数１５】
2005080124.tif 000017
になる。両辺の指数を求めることによって、次の回帰曲線を表す式になる。
【数１６】
2005080124.tif 000018
【００３８】
図１７の各反射次数パスに対して回帰分析を行った結果を、図１８に示す。図１８において、それぞれの反射次数パス別に求められた回帰線は、同じ反射回数のパス群がもつ線形的な関係を表し、その係数から時間に対する傾きが得られる。同じ反射次数の系列であることから、この傾きは反射を除く距離に基づく大気減衰や拡散音場などの要素が含まれる。すなわち、音源の指向性，反射時の入射角による吸音率の差，壁の材質による吸音率の差などの要素を除けば、距離に比例する完全な直線（対数）上に並ぶ。いいかえると、回帰線に対するデータの縦方向の分散は、音源の指向性などを含む上記の要素のばらつきを表す指標であるといえる。
このばらつきには、壁の角度，材質のばらつきなどが含まれることで、材質の特性および部屋の形状的な特徴を表し、残響の推定のための重要な要素であると考え、反射次数別にそれぞれのパスの時間，強さに関する分散を考慮する。これに関して次に述べる。
【００３９】
（反射次数パス別の時間／強さにおける分散）
データの分散を計算する際に、分析する次元を定めなければならない。例えば、平面状のデータならば、Ｘ方向やＹ方向の分析などの１次元的な分散と、重心から放射状の分析を行う２次元的な分散とがある。本システムでは、Ｘ方向，Ｙ方向に分けて１次元的な分散の解析を行う。その理由は、時間的な（横方向）分散と音の強さの（縦方向）分散は独立しているからである。前述の通り、縦方向の分散は部屋のさまざまな特徴が含まれており、これを２次元的に分析すると、その情報は失われるため、縦方向，横方向に分けて分析する。なお、時間と音の強さは直線的な関係をもち、回帰分析によって求めることができる。
【００４０】
２次元グラフにプロットされた反射パスｋのデータ（Ｘｋ，Ｙｋ）に対する縦方向の分散は、以下のように求める。
まず、各測定値の平均値からのずれの総量が０にならないようにするために、平均値からの変位の二乗和をとる。要素数をＫ，音圧レベルをＸｋ（ｋ＝１，２，…，Ｋ）とすると、音圧レベルの分散ＶＹは次式で定義される。
【数１７】
2005080124.tif 000019
すなわち、ｌｏｇＹの変動ΦｌｏｇＹをＫで割ったものである。
しかし、上記に示す分散は平均値バーｌｏｇＹからの分散であり、求めたいのは回帰線からの分散である。したがって、バーｌｏｇＹに回帰線関数のＸｋに対するｌｏｇＹを代入すればよい。よって、回帰線からの分散Ｖ’Ｙは次式のようになる。
【数１８】
2005080124.tif 000020
通常の縦方向分散ＶＹは重心に対する変動を元に計算する。すなわち、（対数表示の）直線ｌｏｇＹ＝￣Ｘ（￣Ｘは定数）からの分散を表す。一方、Ｖ’Ｙは（対数表示の）直線ｌｏｇＹ＝ａＸ＋ｂからの分散であり、前述の通り空間形状、音源の指向性、壁の材質などの特性をもっている。この分散を残響の予測に利用することで、幾何学的に求められた初期反射の上記の性質を継承させることが可能である。
【００４１】
以上、縦方向（音圧レベル）の分散について述べた。次に、時間的な分散について述べる。音は高次反射になるにつれ、その遅れ時間のばらつきも大きくなる（図１７参照）、したがって横方向の分散も解析する必要がある。上述のＶＹの式と同様に、横方向の分散ＶＸは次式のように求められる。
【数１９】
2005080124.tif 000021
図１７の各反射次数パスに対して、それぞれの縦横方向の分散ＶＹ，ＶＸ，および回帰線からの縦方向の分散Ｖ’Ｙを求めた結果を表１に示す。
[表１]反射パスの反射次数別音圧レベルの分散
2005080124.tif 000022
【００４２】
（反射次数別パスデータ（時間／強さ）の重心）
図１７のグラフを反射次数別に眺めると、データが規則的に移動していることがわかる。これは、各反射次数に対して平均的なパスの距離が著しく変化するからである。また、大気伝播による損失よりも反射時の損失が大きく、反射次数に対する音圧レベルも大きく変化する。その結果、反射次数が増加するにつれてデータはグラフの右下の方に移動する。
このデータの移動に基づいて、幾何学的に計算された１〜ｎＥ次反射パスの重心からｎＥ＋１次反射パスの重心を推定する。重心は、それぞれのデータ（Ｘｋ，ｌｏｇＹｋ）の平均値である。ｎＥ＋１次反射パス群の重心は、各反射次数に対する重心の近似線（つまりその回帰線）及び平均移動量を求めることによって、推定することができる（図１９参照）。
【００４３】
＜分散の適用＞
最後に、分散の適用手法について述べる。
これまで説明した計算によって、反射次数別に求められた分散Ｖ’ＹとＶＸから、分散の比ＰＹ’，ＰＸ  を算出し、これをもとに、すでに計算されたデータを移動、分散させる。なお、縦と横方向の分散を独立に適用するために次のように計算している。
まず、各データを回帰線に対して縦方向にＰＹだけ分散させる。次に、ＰＹ’を保存するため、回帰線方向に限定した上で重心に対する分散ＰＸを適用する。すなわち、図１９に示すとおり、点Ｕを縦方向に分散させた点をＵＹ，ＵＹを横方向に分散させた点をＵＹＸとしたとき、ＵＸを通る回帰線に平行な直線上の点Ｕ’ＹＸになるように補正を行う。これによって、縦方向の分散が保存される。
これまでの計算の手順を、以下にまとめる。（ｎＥまでのデータは幾何的に算出されたとする。）
１）  １〜ｎＥ次反射音群の重心を計算
２）  ｎＥ＋１次反射音群の重心を推定
３）  １〜ｎＥ次反射音群の回帰線を計算
４）  ｎＥ＋１次反射音群の回帰線を推定
５）  １〜ｎＥ次反射音群の分散および分散比計算
６）  ｎＥ＋１次反射音群の分散を推定
７）  １〜ｎＥ次反射音群の重心を移動し、分散の調整（および補正）を行う
【００４４】
以下に、本システムによる残響予測の実行例を示す。
図２１（ａ）では、１〜２次反射まで計算した応答に対し、３，４次反射の応答を推定した。また、図２１（ｂ）では、１〜３次反射まで計算した応答に対し、４次反射の応答を推定する。図２１（ｃ）では、４次反射まで計算している。
結果的に、左側の３つのインパルス応答は、１〜４次反射の応答であり、推定部分を除いて述べると、上から２反射まで、３反射まで、４反射まで幾何学的に計算した結果を示す。これから、４反射までの推定による残響予測は、計算による残響結果と同様であり、推定による残響予測を用いることができる。
なお、１５次反射、２５次反射の推定値の計算結果から、１５次反射と２５次反射はほとんど変わらなく、これ以降の計算は（元になった部屋において）不必要であることを意味する。
【００４５】
＜システムへの実装＞
これまでに説明した原理をまとめ、音響再現システムへの実装について述べる。幾何音響に基づく室内音響の再現は、直接音及び反射音を求めることによって行われる。前述の通り、室内の音響を音響伝達関数としてあらかじめ計算しておくことで、任意の音声に畳み込むことで音響を再現できる。
室内における音響を伝達関数として表現する際、その係数を求める必要がある。すなわち、ある位置（音源および受音点）における音響伝達関数を保存することにより、任意の入力信号を幾何的な演算と切り離して、即座に変換することができる。よって、音源を固定し、室内での受音点の各位置の音響伝達関数を保存すれば、受音点の位置に相当する音響伝達関数を用いて変換ができ、移動しつつ音響の変化を再現できる。しかし、室内での全ての点における音響伝達関数を保存することはメモリ量を考えても非現実的である。ここでは、空間全体をグリッド分割し、各セルの音響伝達関数を保存することによって、移動するユーザ（受音点）の位置によって適用する伝達関数を使い分け、リアルタイム再現を行っている（図２２参照）。
【００４６】
図２２に示すようなグリッド分割を用いる場合、グリッドの分割サイズが大きな問題になる。また、音響伝達関数を保存する場合、グリッド内のどの点を代表点として計算するのかを考慮する必要がある。
グリッド分割の際、セルサイズが大きくなるにつれ、エッジでの変化が著しいため、本システムでは、パスの反射次数が高いほど、反射音の明瞭度が低くなる性質を利用して反射次数別に分割のサイズを変更している。図２３では、左から１次反射，２次反射，３次反射パス計算に用いるグリッドを示す。図２３に示すように、反射の次数が高くなるにつれて、使用するグリッドが大きくなっている。受音点位置のそれぞれの反射次数ごとのグリッドに保存されている音響伝達関数を読み出し、合成して音響伝達関数として使用する。
音響伝達関数の計算における受音点の位置に関しては、現時点ではセルの中央点を計算の対象としているが、これ以外にさまざまな方法が考えられる。その１つとしては、セル内の４点を計算し、平均を取る方法が挙げられる。この場合は周りのセルに近い点で計算するため、エッジでの急激な変化を多少抑えることができると考えられる。
【００４７】
これまで、グリッド分割による再現を行ってきたが、この際、いくつかの問題が生じる。グリッドの形状及び大きさによって不自然な音響の変化が目立つなど、グリッド間での連続性が保たれないことがしばしばある。この解決法として、セル内での伝達関数の動的化を考える。セル内での受音点の移動に伴い、パス距離を更新することにより、セルの境界部分での誤差を和らげることができ、音響の大きな変化が避けられる。これに関しては、次に詳しく述べる。
【００４８】
（音響伝達関数の動的化）
前述の通り、グリッド分割において各セルに音響伝達関数を保存する場合、セルのエッジ部分で不具合が生じる。これは、セルの領域内では音響が変化しないため、その分セル同士の境界で、反射音の音圧，遅延および，周波数特性の変化が大きくなる。
これを解決するため、本システムでは音響伝達関数の動的化を行った。セル内の領域を代表点（受音点）に近い領域として考え、計算された反射パスの距離は変化するものの、反射経路はさほど変わらない。これを利用して、図２４に示すように、今までセルの代表点の音響伝達関数の保存を行っていたのに対し、各セル全体を受音点として、このセルに到達するすべてのパスのＩＤをつけ、音響伝達関数ではなくＩＤをつけたパスリストをディスクやメモリ等に保存する。ここで用いるＩＤはパスの反射経路である（１つの反射経路には唯一のパスが存在する）。これによって、使用メモリ量も抑えることができる。
なお、図２４に示したＩＤリストの例は、各壁にＡ〜Ｍの符号をつけ、反射する壁の符号を利用して、各パスのＩＤとしている。例えば、ＡＤは、音源Ｓ→壁Ａ→壁Ｄ→受音セルのパスを示している。
【００４９】
各セルのパスリストの保存後、受音点の位置から、セルのパスリストを探索し、パスと壁面の衝突判定を行いつつ、動的に音響伝達関数を求める。したがって、パスの生成／探索を前処理として行うことで、リストに含まれるパスのみに対して更新を行うことができ、無駄な処理を省くことができる。また、Beam Tracingとは異なり、領域が長方形（長方体）であるため、現在位置からセルを判定するのが容易である。音響伝達関数の動的化によって必要な情報のみを更新できるため、リアルタイム再現には欠かせない要素であるが、図２５に示すような小範囲で音響が変化する部屋の形状においては、リストが大きくなるため、効率が低下する。これは、後で述べる分割法を用いることによって解決できる。
音響伝達関数を用いた再現のプロセスは３段階に分けることができる。はじめに幾何的な計算，次に音響伝達関数の算出，最後に伝達関数による再現である。上記の動的生成は、このうちの幾何的な計算を省略し、保存された幾何的な情報に基づき反射パスを更新する手法で、音響伝達関数の算出および最終的な再現は動的に行う。すなわち、反射次数が高く、パスの数が多ければ当然計算量も増加する。したがって、高次数反射に対して行うとリアルタイム性を失う恐れがある。
【００５０】
本システムでは、リアルタイム性を考慮し、音響伝達関数の算出処理をパスの反射回数に応じて、静的または動的計算に分けて処理する。反射回数が多い場合は、静的生成で音響伝達関数を保存するのに対し、反射回数が少ない場合は、動的生成でパスＩＤのリストを保存する。静的または動的計算から得られた音響伝達関数に対して音声信号を入力し、合成信号を求めることによって、高次数反射にも対応する。
【００５１】
＜可聴化システム(AURALIZATION SYSTEM)＞
本音響再現システムの最終段階は可聴化処理である。可聴化処理では、入力された音声に対し音響的な変換をした後、音声としての出力を行う。無限長の入力音声に対する変換処理は、ブロックに区切って行う必要がある。つまり、入力音声から一つのブロックを読み込み、変換した結果を再生バッファに送って再生する。このバッファ制御システムについて述べる。
【００５２】
（マルチバッファ・システム）
本システムでは、フーリエ変換による帯域変換を用いて、周波数帯域で音響伝達関数を適用している。周波数帯域で処理を行っている理由は、残響や拡散などの効果を適用しやすく、これによって、リアルタイムの自由度が向上するからである。
フーリエ変換を用いた手法のひとつの問題点として、残響時間がほぼ無限であるフィードバックを用いた手法に対して、残響時間が制限されることが挙げられる。無限長のフーリエ変換は不可能であり、入力信号を一定時間のブロックに区切って処理を行う必要がある。このブロックに対するフーリエ変換のサイズはサンプリング定理によりその２倍以上とる必要がある。図２６では、音声からデータを切り取り、帯域変換及び逆変換を行った例を示す。このとき、切り取られたブロックがＦＦＴバッファサイズの半分を超えたサイズである場合、エイリアシングが生じる。
【００５３】
しかしながら、本システムでは音響伝達関数を用いて、特性フィルターだけではなく波形の時間的な推移を行っていることが一つの特徴でもある。すなわち、ＦＦＴバッファに格納されたデータを伝達関数によって変換を行ったとき、信号読み取りブロックサイズをＦＦＴバッファサイズの半分に設定した場合、推移した音声はＦＦＴバッファサイズを越える。よって、図２７に示すように、エイリアシングが生じることは明らかであり、音声の復元が不可能になる。これを避けるため、最大の推移時間を考慮した上で切り取りブロックのサイズを定める必要がある。
したがって、ＦＦＴバッファにおける最大遅延時間Δｔは
【数２０】
2005080124.tif 000023
になる（図２８参照）。残響時間を長くするためには、このΔｔを大きく取る必要があり、これに比例してＦＦＴバッファサイズも増大し、膨大な計算量になる。
本システムでは、この問題を改善するために、マルチバッファシステムを採用した。計算量を考えると、いくつかの小さいＦＦＴバッファに分けて計算した方がよい。そこで、インパルス応答を一定時間に区切り、これらを小さいＦＦＴバッファ複数個に設定して音響伝達関数ＡＴＦを求める。これによって、計算機の処理能力に応じた任意の残響時間を計算することができる。
【００５４】
このシステムの概略を図２９に示す。図２９に示すシステムでは、入力信号のサンプリングレートを１６３８４Ｈｚであるとして述べる。４個の１６３８４サンプルのＦＦＴバッファ５０２，５０４，５０６，５０８を用いたシステムを表している。
まず、入力音声をブロックに区切り、それぞれの音響伝達関数ＡＴＦによって変換する。この図２９の例では、ブロックサイズと遅延時間が最大となるように、ブロックサイズ（samples）＝１ＦＦＴバッファの最大遅延（samples）としている。そのため、４個のＦＦＴバッファを使った場合の最大の遅延は４０９６×４＝１６３８４samplesになる。最後に、４個のＦＦＴバッファ５０２，５０４，５０６，５０８の足し合わせたデータを再生バッファ５１０に送り、再生処理を行う。
【００５５】
（再生バッファシステム）
本システムの最終プロセスとして、計算された出力波形を実際の音として出力する。マルチバッファから出力された波形の管理について次に述べる。
再生バッファも複数のストリーミング・バッファから構成される。ＦＦＴバッファの数をＮＢ，そのサイズをｌＢとして、再生バッファの数はＮＢ＋２，サイズはｌＢ／４になる。図２９では、１つのＦＦＴバッファに含まれる音声のサイズはｌＢ／２であり、このうち半分はとなりのバッファと共有するため、計算されたデータを格納するためにはＮＢ＋１個のバッファが必要となる。計算中にこれらのバッファは全てロックされるため、再生を可能にするためには、フリーとなる一つの余分なバッファが必要である。したがって、ＮＢ＋２個の再生バッファが必要になる。図２９のシステムの場合、４つのＦＦＴバッファ５０２，５０４，５０６，５０８で構成されるため、図３０に示すように、６つ分の再生バッファ５１０が必要になる。
【００５６】
再生ポインタを除く残りの５つのバッファはロックされ、計算されたデータが格納される。これらのバッファはループ構造をしており、末端バッファを終了すると、再生ポインタは先頭のバッファを指す(図３０参照）。これによって、無限長の波形を、バッファ格納処理（ロック）と干渉することなく再生できる。データ格納の流れを、図３１に示す。この図示した処理を制御するために、再生ポインタがバッファの境界を指したとき、イベントを発生させ、計算を開始する。すなわち、リアルタイム再現を可能にするためには、バッファのサイズを時間単位で表した時間（図２９の場合は１／４ｓｅｃ）以内に音響再現処理を終了させる必要がある。
【００５７】
＜システムの応用＞
以下に、いくつかの点で、システム全体の効率化を実現するための手法について述べる。
（時間帯域での室内伝達関数の生成）
前に述べたシステムでは、反射パス情報リストを分割されたセルに対して保存することによって、パストレース・プロセスを省略でき、動的に室内伝達関数の計算が可能になった。このとき、算出処理は周波数帯域で行われ、室内伝達関数は次の式のように表される。
【数２１】
2005080124.tif 000024
（ただし、αｋ（ω）とｘｋは、それぞれｋ番目のパスの損失特性と遅延である。）
しかし、上述の式では、乗算回数はｋ×Ｎになる（Ｎはサンプル数）。ゆえに、最も計算時間を要するパストレース処理を省略したとしても、ＲＴＦ生成の処理時間の問題が残る。よって、ＲＴＦ算出の高速化が求められる。ここでは、ＲＴＦに関して時間帯域での算出を行うことによる高速化について述べる。
【００５８】
室内伝達関数の算出処理を周波数帯域で行う理由は、反射パスの損失特性（周波数特性）が周波数帯域の関数であるからである。したがって、それぞれの反射パスの損失特性αｋ（ω）に遅延ｅｉωｘｋをかけ、全反射パスの合成関数を求めればＲＴＦ（ω）が得られる。
周波数特性は、計測データとして１２８Ｈｚ，２５６Ｈｚ，５１２Ｈｚ，１ＫＨｚ，２ＫＨｚ，４ＫＨｚ，８ＫＨｚの損失率で表される。これを特性パラメータと呼ぶ。特性パラメータを線形または，スプライン関数で補間して使用する。８点を操作点として、１６３８４点の周波数帯域上の関数として表わされ、これが損失特性として用いる。図３２は、表２の材料に対する周波数ごとの吸音特性（実測データ）に基き、補間された特性関数（図３２（ａ）：線形補間，図３２（ｂ）：B-Spline補間）を示している。
[表２]建築材質の吸音率（特性パラメータ）
2005080124.tif 000025
ここで、補間された周波数特性関数を三角関数で表現することを考える。補間を行ったあと（図３３の周波数特性の線形補間例参照）、生成された関数に帯域変換をかける。その結果を図３４に示す。図３４に示すように、周波数成分ごとの三角関数を得ることができる。フーリエ変換の可逆性より、図３４の関数は時間帯域の関数であるといえる。
【００５９】
次に、周波数特性を表す操作点の間隔を考える。波形の特徴的な形状を崩さないためには、最低限として、１００サンプル長の波長を持つ三角関数が必要である。つまり、１６３８４／１２８＝１２８Ｈｚ以下の周波数成分が必要である。したがって、帯域変換によって得られた成分のうち、１〜１２８Ｈｚ成分以外の成分を放棄し（ローパス・フィルタを用いる）、逆変換をかければ、形は特徴を失わずにほとんど変わらない。なお、音響的な変化として、周波数特性の詳細な精度よりも各周波数の成分比率が重要となる。
図３５の周波数特性対数表示の図で示すように、１２８Ｈｚ以上の成分をカットした場合、元の波形と比較しても波形の特徴は変わらず、ほとんど一致する波形を生成できることが分かる。同様に、６４Ｈｚ以上の成分をカットした場合、周波数成分のグラフの１２８Ｈｚ付近では形が崩れてきている事が分かる（図３５を参照）。ゆえに、最低１２８Ｈｚ以下の波形が特徴的な形を復元するのに必要であると考えられる。実際に、図３６の周波数特性波形を表現するのに必要となるデータの図では、１２８Ｈｚ以上の成分はほとんどない。
【００６０】
このように、周波数特性が簡略化され、時間帯域上の１２８サンプルサイズの波形で表現されることによって、次のことがいえる。
フーリエ変換の畳み込みの性質から、周波数帯域関数の乗算は時間帯域関数の畳み込み演算と等価である。すなわち、図３６に示す１２８サンプルの関数に遅延関数を畳み込み、合成関数を求めれば、時間帯域上の室内伝達関数、つまり室内インパルス応答（RIR:Room Impulse Response）が得られる。
ここでいう遅延関数は、前式のｅｉωｔを時間帯域に直したものである。これは、ｔだけ遅らせた単位パルスになる。したがって、上記に述べた畳み込み演算は１２８サンプルのシフト演算ということになる。以前のシステムでは、１６３８４×ｋ（ｋはパス数）の計算量を要するのに対し、時間帯域で計算および周波数特性の簡略化を行うことによって、１２８サンプルのシフト演算のみで計算でき、なおかつ、周波数特性が考慮される。ゆえに、以前より明らかに高速化されることがわかる。
【００６１】
以上、時間帯域でＲＴＦの算出を行う手法による高速化について述べてきた。上述の手法では、周波数帯域関数である損失特性を帯域変換によって時間帯域の関数に直し、後部（１２８サンプル以降）をカットしても周波数帯域上の関数にはさほどの影響がないことがわかった。これによって、小さいサイズの配列シフト演算でＲＴＦを求めることができることから高速化が実現された。しかし、上記の事柄は、反射音が一つの反射経路に対し固定された特性をもつことを前提としている。すなわち、反射吸音による損失計算に対しては有効であるが、大気伝播による損失計算のような、反射パスの距離に依存した特性の適用に対しては後で述べているとおり、逆効果である。
このことから、上記の手法は一見否定されるように思えるが、周波数特性を複数の窓関数を用いて表現することにより、パスの距離に依存した特性に対しても有効となり、なおかつ、更なる高速化が可能になる。これについて次に述べる。
【００６２】
（窓関数を利用した周波数特性の表現）
前述の通り、周波数特性を時間帯域で計算する場合の一つの難点として、音の大気伝播による損失を計算する際に、上述の式に示す周波数特性は距離に依存した特性になるため、前処理として計算することが難しいことが挙げられる。その理由は減衰の適用方法にある。大気伝播による損失は通常、単位距離に対する減衰率として扱われる。したがって、距離ｌに対する減衰率ｍはｍｌになる。ここで、大気伝播による損失は周波数特性をもっていることから、波形ｆ（ｔ）の減衰を求めるためには、ｆ（ｔ）のフーリエ変換Ｆ（ω）に周波数成分別の減衰率Ｍ（ω）のべき乗Ｍ（ω）ｌをかける必要がある。周波数帯域でのべき乗は時間帯域との互換性はないため、次式のように距離に比例した畳み込み演算回数を行う必要がある。
【数２２】
2005080124.tif 000026
その結果、時間帯域での演算は逆に時間がかかり高速化の意味がなくなる。本システムでは、これを解決するために周波数特性曲線の制御点を独立させることで、周波数特性関数Ｍ（ω）と距離ｌを定数化する手法を採用する。
【００６３】
前述では、図３２（ａ）に示すとおり、測定された吸音特性パラメータを通過点とするB-Spline補間を用いて、周波数特性関数を前処理として生成する手法について述べた。この手法を用いると、時間帯域での計算を行った場合には非常に時間がかかる。これは、曲線補間を行った時点で全ての特性パラメータが同じ関数上に配置されることが原因である。したがって、本システムでは曲線補間処理を前処理ではなく、リアルタイムで行うことで、特性パラメータ同士を反射音の再現の段階まで独立な状態に保つことを考えた。その結果、上式に示した畳み込みの問題が解決される。
しかし、B-Spline補間などの逆変換を用いた補間手法は時間がかかるという問題がある。また、周波数特性関数は時間帯域の関数ではないため、時間帯域での計算においては曲線補間を行ったあと帯域変換（逆フーリエ変換）を行う必要がある。これらの問題点を解決するため、窓関数を利用した時間帯域で周波数帯域上の曲線補間を行う手法を用いている。
曲線補間に用いられるB-Spline  曲線やBezier  曲線は重心結合（Barycentric combination）に基づいている。重心結合は複数の関数（基底関数）によって構成されており、以下の式で定義されている。
【数２３】
2005080124.tif 000027
ただし、Ｐｉはベクトル、Ｂｉ（ｔ）は基底関数、つまりスカラーである。
【００６４】
基底関数の和が１であることから、周波数特性で考えれば無特性状態、すなわち、あらゆる信号にかけても変わらない状態を表せる。基底関数を高さ１の関数として、周波数特性の各特性パラメータに対し基底関数を一つ用意することで、それぞれの基底関数を定数倍することによって、周波数特性を線形結合で表すことができる。この高さ１の関数を窓関数として扱い、時間帯域で各定数値を変更することで周波数特性を操作できる。すなわち、時間帯域で周波数帯域の曲線（周波数特性）を、曲線補間手法を用いることなく自由に操作できる。
これにより、距離ｌは定数化し、周波数特性はｌに依存しなくなる。
吸音特性の特性パラメータＭ１，Ｍ２，…，Ｍｎ（ｎは特性パラメータの数）を、前処理として補間した場合の曲線をＭ（ω）とすると、周波数特性Ｕ（ω）は
【数２４】
2005080124.tif 000028
になる。一方、線形結合によって周波数特性パラメータを表すときの各高さ１の関数をそれぞれＢ１(ω），Ｂ２(ω），…，Ｂｎ(ω）としたときの周波数特性Ｕ’（ω）は
【数２５】
2005080124.tif 000029
になる。
【００６５】
関数Ｕ’（ω）とＵ（ω）は必ずしも一致しないが、同じ特性パラメータを通る。もともとＵ（ω）は正しい周波数特性ではなく、補間によって生成されたものであり、点Ｍ１，Ｍ２，…，Ｍｎを通る曲線を求めることが本来の目的である。したがって、Ｕ’（ω）は別の補間手法であると考えれば、望む結果であることがわかる。
図３７では、２つの手法の比較を図示している。図３７の窓関数を用いた場合、基底関数の定数倍したものをそれぞれ別々に帯域変換を行っている。フーリエ変換の線形性から、この定数の処理は時間帯域で行っても同じであるから、Ｍ１，Ｍ２，…，Ｍｎをｌ乗すればそれぞれの窓関数（時間帯域上の基底関数）にかけられる定数になる。このことから、以前に述べたｌ乗の問題は解決されることがわかる。
【００６６】
（頭部伝達関数の単純化手法を利用した実装）
リアルタイム音響再現において、音場予測手法として幾何音響を用いる利点は、反射パス（反射音）の到来方向がわかることである。反射パスの到来角度からＨＲＴＦ（頭部伝達関数）（図１３）を適用することにより、反射音の音像定位まで再現することができる。すなわち、立体的な音場を再現することが可能になる。ＨＲＴＦは反射物が存在しない自由空間における伝達関数を表し、音源から到来する反射音をそれぞれ独立に扱う幾何音響と合成すれば、空間形状に基づいた３次元音場情報をリアルタイムでリスナーに伝えることができ、グラフィックスだけでは表せないような臨場感をユーザに与えることができる。
本システムではＨＲＴＦを用いて音場を再現している。仮想空間でのユーザの指向方向に対する音源の位置を指すベクトルの角度および、反射パスが受音点に到来する角度から、あらかじめメモリに読み込まれたダミーヘッドの頭部伝達関数データを元に音像定位を再現できる。
さて、角度に該当する伝達関数を入力音声に畳み込むことによって音像を再現できるが、パスによって到来角度が様々である。そのため、パス数の分だけ畳み込みが必要になる。頭部伝達関数のインパルス応答は５〜６ｍｓ以内でサイズが小さいとはいえ、１００〜１０００本にもわたるパスに適用すれば、当然リアルタイム性が失われる恐れがある。そのため、本システムではＨＲＴＦの適用レベルを設定している。適用レベル以外の反射パスの扱いについては、以降に述べる。
【００６７】
図３８では、ＨＲＴＦを適用するバイノーラル化モジュールの流れを示す。各反射音に対し、吸音・減衰・遅延等の処理よりパスの伝達関数（以降ＰＴＦと呼ぶ）を生成後（Ｓ６０４）、ＨＲＴＦ適用レベルの判定および角度の算出を行い（Ｓ６０８）、ＰＴＦに畳み込む（Ｓ６１０）。これにより、生成されたインパルス応答が定位を持った反射音の音響係数になる（Ｓ６１２）。最終的にすべての係数をＨ’１〜Ｈ’ｎの和を取る（Ｓ６１４）ことで、バイノーラルな音響伝達関数を算出できる（Ｓ６１６）。
ＨＲＴＦ適用レベルに除外されたパスは音像定位計算を省略するわけではない。本システムでは、数多くの反射パスに対してＨＲＴＦを適用するため、ＨＲＴＦ関数の簡略化を提案し、減衰して変化が目立たない高次反射パスに適用した。この原理について以下に述べる。
【００６８】
図３９に示す時間帯域の頭部伝達関数を３軸方向から眺めると、図４０（ａ）〜（ｃ）になる。このいくつかの特徴的な形状を次のようにまとめた。
（１）角度によって、インパルス応答の全体の遅れが三角関数的に変化している（図４０（ａ））
（２）角度によって、全体的に音量が変化している（図４０（ｂ））
（３）角度ごとのインパルス応答が似たような形状をしている（図４０（ｃ））
音響の変化を認知するとき、遅延，音量は最も影響する要素であり、詳細な形状（特性）よりも優先すべきであると考え、図４０で見られるインパルス応答の全体的な遅れと音量に着目した。
これが立体的な音場を感知するための重要な要素であり、どんな頭部伝達関数においても共通な特徴であると言える。
【００６９】
ＨＲＴＦの頭部形状に依存しない遅延時間である（１）に関しては、次のように解釈している。ＨＲＴＦデータの測定はダミーヘッドの周囲円上に音源を配置しており、右耳に配置された受音器は円の中央にないため、当然この円上の音源からの距離が異なる（図４１）。その結果、図４０（ａ）のように到達する音の遅延にも変化が出てくる。角度に対する音源−受音点の距離を示す関数ｌｓｒは次のようになる。
【数２６】
2005080124.tif 000030
音源：（cosθ, sinθ）  受音点：左（０，１０）  右（０，−１０）
この関数は図４２のようになる。これを比較すると、測定データと同じ形であることが分かる。
【００７０】
（２）に関しては、次のように解釈できる。図４１に示すとおり、真右の音源から右耳に到達する場合は頭の影響がほとんど無いのに対し、真左の音源から右耳の音に到達する音は頭を真っ向にすり抜けるため、その損失は大きい。つまり、回折を除き、透過による損失が大きく影響していることが分かる。ここで、図４１のように頭を円として考え、音源から受音点（右耳の鼓膜）へ直線を引き、円内部に含む線分の距離は透過損失に比例するので、これに従って式にまとめると以下のようになる。
まず、図４３に示している音源−受音点を通る直線は以下の式となる。
【数２７】
2005080124.tif 000031
この直線と頭に相当する半径１０の円との交点は、
【数２８】
2005080124.tif 000032
である。この２点のうち、音源に近い方をｑ（ｑｘ，ｑｙ）として、音源−受音点を結ぶ線分から頭部（半径１０の円）の内部を通る線分の長さｌｈは、
【数２９】
2005080124.tif 000033
である。このとき、交点は必ず２つ以上存在するとしている。この結果を図４４の平均的損失量として示す。
【００７１】
以上のことをまとめて，計算により生成されたＨＲＴＦは図４５のようになる。図４５（ａ）では、左鼓膜に到達する代表的なＨＲＴＦとして、遅延と損失の情報を表している。ここで、左鼓膜に対して頭部の影響が最も小さくなる真左の音源からのインパルス応答ＲＩＲ（図４５（ｂ））を考え、これをＨＲＴＦに畳み込むことにより左鼓膜へ到達したインパルス応答が求まる（図４５（ｃ））。この計算によって求められた結果（図４５（ｃ））は、実際に測定されたＨＲＴＦに対応するインパルス応答（図４５（ｄ））と同じ形として得られている。
【００７２】
次に、システムへの実装について述べる。伝達関数が単純化されたとはいえ、すべてのパスに対して畳み込むことは精度が落ちるだけで、計算量的には全く無意味である。重要なのは、反射パスの到来角度から遅延時間と損失量がわかることである。これによって、前のフローチャート（図３８参照）を図４６のフローチャートように変更することができる。図４６と図３８のフローチャートを比較すると、ｎ本の反射パスの場合、畳み込み演算の回数はｎ回であったのに対し、遅延・損失量の訂正と１回の畳み込みで済ますことができる。これについて以下に詳しく述べる。
【００７３】
すでに述べたとおり、音響伝達関数ＡＴＦが室内伝達関数ＲＴＦと頭部伝達関数ＨＲＴＦによって構成されている。ここでＲＴＦは音源から受音点へ到達する様々な反射パスの伝達関数ＰＴＦの合成関数であり、時間帯域上では単純な和として表される。一方、ＨＲＴＦは、図１３に示すような角度に依存する伝達関数となっており、これを時間帯域へ変換することにより、頭部インパルス応答（HRIR:Head-Related Impulse Response）が得られる。これらのことから、lossｋをパスｋの総合的な損失，Δｔｋをパスｋの到達遅延時間，φｋをパスの到来角度，ｓｒｃ（ｔ）をドライソースの信号系列としたとき、鼓膜で聞き取られる音響を表すＡＴＦに対応するインパルス応答ＲＩＲ（ｔ）（室内インパルス応答）は次式のように表せる。
【数３０】
2005080124.tif 000034
ＨＲＴＦを単純化し、代表的なインパルス応答ｒＩＲ（representative Impulse Response）とその遅延（前述の式のｌｓｒ）及び損失（前述の式のｌｈ）で表すこととしたとき、単位振幅のインパルス応答ｒＩＲに遅延と振幅を掛け合わせることで、擬似的な頭部伝達関数が得られる。これを用いて上式を表すと、次式のようになる。
【数３１】
2005080124.tif 000035
この式は図３８のフローチャートに対応しており、数式上の変形によって次式のように書き直すことができる。
【数３２】
2005080124.tif 000036
上式は、図４６のフローチャートを表している。
ＨＲＴＦのデータは個人差があり、説明した計算法は決して厳密ではない。そのため、直接音や初期反射に対して測定データの畳み込みを用い、減衰した音に対しては簡略された手法を適用している。
【００７４】
（音響の変化に基づいた空間分割）
前述では、グリッド分割による反射パス情報リストの保存について述べた。空間分割は、通常、位置判定が容易であることから、グリッド分割が用いられることが多い。しかし、音響再現においては、受音点の場所によっては音響が激しく変化し、セルの大きさはこれに対応できない場合が発生する。この一例を図２５に示した。ここでは、空間をグリッド分割したときに同じセルに到達する反射パスが大幅に変化し、セルに保存される反射パス情報リストが大きくなる例を示している。
本システムでは、この問題を解決するために音響の変化に基づく分割法を考慮する必要があると考え、コンピュータ・グラフィックスで形状表現法としてよく用いられる四分木を利用する。
音響の変化に関連する要素は様々であるが、主に反射音の遅れ，音圧レベル，周波数特性は大きく影響している。これらを考慮した上で分割を行う手法を提案する。本システムでは四分木分割を用いることによって、これらの問題を解決した。この適用方法について以下に述べる。なお、本システムに用いられる分割手法は、２次元分割を行っている。これは、ユーザは平面上にしか移動していないことを前提としているからである。３次元への拡張を行う場合、四分木は八分木(Octo-tree)になる。
【００７５】
本システムでは、音響の変化が著しい場所では分割を詳細にし、変化が微小である場所では分割を粗くすることによって、メモリを効率的に使っている。これに適した分割法として、四分木分割が挙げられる。
四分木（Quad-tree）とはピクセルの集合の２次元図形を表す手法である。平面全体をそれぞれ縦横半分に分割すると４つの領域に分割される。各領域に図形がすべて占めるかあるいは全くない場合は分割を中止し、それ以外の場合は各領域をさらに４つの領域に分割していく。このような再帰的手法により全ての図形は図４７のような四分木のツリー構造で表現できる。図形の形状にもよるが、四分木表現によりデータ量を縮小させることができる。
音響の比較対照はさまざまであるが、反射パス情報リストの保存が本分割の目的であることから、本システムではインパルス応答に着目する。上記に述べた分割手法の四分木の分割条件を、一定以上の音響の変化量に置き換えれば、図４８に示すように空間を分割できる。
インパルス応答といっても、これ自体には音響の変化に関わる要素がたくさん含まれる。したがって、本手法で用いる比較の要素を明確にする必要がある。比較の対象とする要素としては、音圧レベルの変化，インパルス応答の多次元ベクトル内積がある。以下にそれぞれを詳しく説明する。
【００７６】
（音圧レベルの解析）
インパルス応答から、パルスに対する音圧レベルを図ることができる。これはつまり単位的な指標として考えられる。高さ１のインパルスの積分（面積）は１として（一般にδ関数という）、これに対するインパルス応答を積分すれば、指定時間内（積分の範囲）のエネルギー量を図ることができる。これを正規化することで、指定時間の平均音圧レベルとして考えることができる。したがって、位置ｐｉ，ｊにおけるインパルス応答の系列をＡＴＦｉ，ｊ（ｔ）として、平均音圧レベルＰＬｉ，ｊは次のようになる。
【数３３】
2005080124.tif 000037
空間を細分割し、すべての格子に対する平均音圧レベルを求めた例を、図４９に示す。
ここで、空間における音圧の変化を比較する際のオペレータについて述べる。本手法では、画像処理のフィルタリングに用いられる１次微分（ソーベルオペレータ，Sobel Operator)を利用する。分割された各格子に対して平均音圧レベルを算出することで、２次元のグレースケール画像と同等に扱うことができる。
画像のおける明るさを平均音圧レベルに置き換えると、空間全体の音圧の変化量を算出できる。図４９の平均音圧レベルから１次微分を求めたものを図５０に示す。
差分を用いた音圧の変化の解析方法について述べてきた。この結果を空間分割のための情報として利用する。しかし、平均音圧レベルの情報だけでは音の遅れや特性等の変化が検出できない。したがって、本システムでは次に述べるアルゴリズムも用いる。
【００７７】
（インパルス応答のベクトル内積）
インパルス応答はパルス音の集合であり、それぞれにパルス（応答）は反射パスの距離に応じた遅れ、または、伝達の経路での損失特性をもっている。したがって、音響の変化を比較する際にこの２つの要素を考慮する必要がある。例えば、図５１（ａ）に示すような例では、比較対照の波形の1 つのパルスが微小な遅延をもつ。この場合は、我々の聴覚にはその変化はさほど分からない。一方、図５１（ｂ）のように、大幅に遅れたときは音が分離して聞こえるため、音響に大きな変化があらわれる。上記の手法（音圧レベル）の場合には、いずれも同等に扱われるため、検出することができない。逆に、線形的に比較した場合は時系列の各値は独立しているため、図５１（ａ）の場合でも大きな変化が現れる。
上記の変化の差を考慮するために、インパルス応答を多次元ベクトルとして扱い、その内積を求めることによって、ベクトル間の角度を算出できる。ｎ次元ベクトルの場合、同じ高さのパルス音が時系列上でシフトしたとき、このベクトルは多ｎ次元空間上の球面上を移動する。すなわち、遅延時間差が小さければ、球面上での移動距離が小さく、結果的にベクトル同士の角度も小さくなる。
このアルゴリズムを用いて、測定地点の周囲の４点もしくは８点と比較し、平均を取ることによって空間全体の変化を解析している。図５０と同じ部屋における、この手法による解析結果を図５２に示す。
上記の２つのアルゴリズムに基づいて空間分割を行う。四分木の上位階層を開始し、ボクセル内に含まれる格子の上記のアルゴリズムによって算出されたデータの平均値を求める。この値が指定値（適用度）を越えた場合のみに再分割する。分割されたボクセルに対し、さらにこの処理を繰り返す。この適用度を変更して実行した分割例を図５３に示す。
【００７８】
この空間分割法を用いることによって，次の利点が得られる。
・音響の変化位置の制度が向上する
空間を移動する場合、セルの境界部分では多少の変化が生じる。この変化位置は空間の分割の詳細度に依存しており、グリッドのサイズが大きければ、その誤差も大きくなる。四分木を用いることによって、変化が著しい場所での分割が詳細になり、精度が向上する。
・ひとつのセルに保存されるパス情報リストのサイズが縮小される
本システムでは、音源からセルに到達するすべてのパスＩＤをリストとして保存する。受音セルのサイズが大きければ、これに伴い到達するパス数が多い。本システムでは、リストに含まれるパスの壁との衝突判定を行っており、音響の変化が著しい場所では、パスが出現と消滅が多いと考えられるため、無駄に保存されているパスが多いことになる。
したがって、変化が著しい場所での細分割によって、衝突判定の回数を減らすことができる上、メモリを効率的に使うことができる。
・静的モデルにおけるセルでの境界で起こる変化が小さくなる
残響など静的モデルとして保存している伝達関数に関して、セル内は同じ伝達関数を用いるため音響は変化しない。逆に、セル境界部分ではその分だけ変化が大きい。これはセルのサイズに伴うため、音響事態の変化が著しい場所での詳細化によって大いに向上する。
以上、音響の変化に基づく分割手法について述べた。
【図面の簡単な説明】
【００７９】
【図１】本実施形態の音響再現システムのフローチャートを示す図である。
【図２】本実施形態における音響再現を説明する概要図である。
【図３】本システムにおけるデータのやり取りを示す図である。
【図４】音響伝達関数生成処理を示すフローチャートである。
【図５】音の反射を説明する図である。
【図６】反射パスの幾何的な算出手法を説明する図である。
【図７】１〜５次反射パスを算出した例を示す図である。
【図８】音の損失を説明するための図である。
【図９】音源の指向性を説明する図である。
【図１０】音響の指向性に対する放出パスの角度を説明する図である。
【図１１】音源の角度に対する特性を説明する図である。
【図１２】反射吸音特性を説明する図である。
【図１３】頭部伝達関数を示す図である。
【図１４】頭部伝達関数適用手法の例を示す図である。
【図１５】対象空間及び音源・受音点の位置を示す図である。
【図１６】エコーダイヤグラムを示す図である。
【図１７】パスの反射次数別の応答（対数表示）を示す図である。
【図１８】パスの反射次数別の回帰分析を示す図である。
【図１９】９次反射音群の重心の推定を説明する図である。
【図２０】分散の適用を説明する図である。
【図２１】２，３次反射パスからの４次反射パスの推定を示す図である。
【図２２】空間のグリッド分割を説明する図である。
【図２３】反射次数別の空間分割を示す図である。
【図２４】セルのパス情報リストを示す図である。
【図２５】グリッド分割の問題を説明する図である。
【図２６】サンプリング定理を説明する図である。
【図２７】エイリアシングを示す図である。
【図２８】切り取りブロックサイズとバッファ最大遅延可能時間を説明する図である。
【００８０】
【図２９】マルチバッファシステムを説明する図である。
【図３０】再生バッファのループ構造を説明する図である。
【図３１】データ格納の流れを説明する図である。
【図３２】材質の吸音特性を示す図である。
【図３３】周波数特性の線形補間の例を示す図である。
【図３４】図３３の周波数特性を帯域変換した結果を示す図である。
【図３５】周波数特性波形の高周波成分をカットした結果を対数表示で示す図である。
【図３６】サンプリングした周波数特性波形の表現に必要になるデータを示す図である。
【図３７】周波数特性の周波数帯域／時間帯域での計算の手法を示す図である。
【図３８】ＨＲＴＦ適用システムのフローチャートを示す図である。
【図３９】時間帯域におけるＨＲＴＦ関数を示す図である。
【図４０】図３９のＨＲＴＦ関数を３方向から見たグラフを示す。
【図４１】ＨＲＴＦの頭部形状に依存しない遅延時間を説明する図である。
【図４２】計算によって求められた遅延時間を示す図である。
【図４３】頭部による平均的な吸音量を説明する図である。
【図４４】計算によって求められた平均的な損失量を示す図である。
【図４５】本手法に基づいた計算によるＨＲＴＦと測定のＨＲＴＦの比較（右耳）を示す図である。
【図４６】簡略化されたＨＲＴＦ適用システムのフローチャートを示す図である。
【図４７】四分木及び四分木分割を説明する図である。
【図４８】四分木分割を利用した例を示す図である。
【図４９】空間における音圧の変化を示す図である。
【図５０】ソーベルオペレータによる、空間における音圧の変化を示す図である。
【図５１】インパルス応答の比較を示す図である。
【図５２】インパルス応答のベクトル内積による、空間における音圧の変化を示す図である。
【図５３】ベクトル内積を比較アルゴリズムとして用いて分割を行った様子を適用度別に示した図である。
【図１】
2005080124.tif 000038
【図３】
2005080124.tif 000039
【図４】
2005080124.tif 000040
【図５】
2005080124.tif 000041
【図６】
2005080124.tif 000042
【図４５】
2005080124.tif 000043
【図２】
2005080124.tif 000044
【図７】
2005080124.tif 000045
【図８】
2005080124.tif 000046
【図９】
2005080124.tif 000047
【図１０】
2005080124.tif 000048
【図１１】
2005080124.tif 000049
【図１２】
2005080124.tif 000050
【図１３】
2005080124.tif 000051
【図１４】
2005080124.tif 000052
【図１５】
2005080124.tif 000053
【図１６】
2005080124.tif 000054
【図１７】
2005080124.tif 000055
【図１８】
2005080124.tif 000056
【図１９】
2005080124.tif 000057
【図２０】
2005080124.tif 000058
【図２１】
2005080124.tif 000059
【図２２】
2005080124.tif 000060
【図２３】
2005080124.tif 000061
【図２４】
2005080124.tif 000062
【図２５】
2005080124.tif 000063
【図２６】
2005080124.tif 000064
【図２７】
2005080124.tif 000065
【図２８】
2005080124.tif 000066
【図２９】
2005080124.tif 000067
【図３０】
2005080124.tif 000068
【図３１】
2005080124.tif 000069
【図３２】
2005080124.tif 000070
【図３３】
2005080124.tif 000071
【図３４】
2005080124.tif 000072
【図３５】
2005080124.tif 000073
【図３６】
2005080124.tif 000074
【図３７】
2005080124.tif 000075
【図３８】
2005080124.tif 000076
【図３９】
2005080124.tif 000077
【図４０】
2005080124.tif 000078
【図４１】
2005080124.tif 000079
【図４２】
2005080124.tif 000080
【図４３】
2005080124.tif 000081
【図４４】
2005080124.tif 000082
【図４６】
2005080124.tif 000083
【図４７】
2005080124.tif 000084
【図４８】
2005080124.tif 000085
【図４９】
2005080124.tif 000086
【図５０】
2005080124.tif 000087
【図５１】
2005080124.tif 000088
【図５２】
2005080124.tif 000089
【図５３】
2005080124.tif 000090