音場測定システム
【課題】多数のスピーカより同時に音響信号を出力させると共に複数のマイクで収録することにより、全てのスピーカおよびマイクに対応する伝達関数を同時に測定すること。
【解決手段】 伝達関数算出手段１４は、複数のスピーカ９より出力される評価音と複数のマイク１０により収録される収録音とに基づいて、時間領域における反復解法を適用することにより、伝達関数の算出処理を実行し、反復解法の適用に生じる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に時間領域に演算結果を逆変換させる演算方法を用いる。評価音生成手段１２は、評価音におけるスピーカ９毎の畳み込み行列を行列要素とし、各スピーカ９に対応させて行列要素を整列させた行列Ａと、行列Ａの転置行列Ａｔとの積ＡｔＡの逆行列が存在し得るように評価音を生成する。
【解決手段】 伝達関数算出手段１４は、複数のスピーカ９より出力される評価音と複数のマイク１０により収録される収録音とに基づいて、時間領域における反復解法を適用することにより、伝達関数の算出処理を実行し、反復解法の適用に生じる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に時間領域に演算結果を逆変換させる演算方法を用いる。評価音生成手段１２は、評価音におけるスピーカ９毎の畳み込み行列を行列要素とし、各スピーカ９に対応させて行列要素を整列させた行列Ａと、行列Ａの転置行列Ａｔとの積ＡｔＡの逆行列が存在し得るように評価音を生成する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音場測定システムに関し、より詳細には、測定音場において複数のスピーカ（音源）から一斉に評価音を出力すると共に、複数のマイクで複数のスピーカより出力された評価音を収録することにより、各スピーカから各マイクに対応するそれぞれの伝達関数を一度の測定で算出することが可能な音場測定システムに関する。
【背景技術】
【０００２】
従来より、ある空間（原音場）の音場特性を収録（測定）した後、その空間とは異なる空間（再生音場）において、測定した原音場の音場環境を仮想的に再現する音場再現方法が提案されている（例えば、特許文献１参照）。
【０００３】
上述した音場再現方法では、まず、ある空間における各スピーカから各マイクへの全ての伝達関数を原音場の音響特性として求めると共に、異なる空間（再生音場）における各スピーカから各マイクへの全ての伝達関数を再生音場の音響特性として求め、この２つの音響特性に基づいて音場再現を行うための音場再現フィルタを算出する。次に、算出された音場再現フィルタを用いて、再生音場において出力される音響信号に対して音響処理（例えば、ＦＩＲ（Finite Impulse Response）フィルタを用いた畳み込み演算処理）を行うことによって、再生音場へ出力される出力音に音響処理を施す。このようにして音響処理が施された音響信号を再生音場で出力することによって、聴取者はあたかも原音場で聴取しているかのような感覚を再生音場で得ることが可能となる。
【０００４】
具体的に、既知の入力信号（入力波形）ｘ（ｔ）を入力（スピーカから出力）させ、この入力に対応する出力信号（出力波形）ｙ（ｔ）を測定する（マイクで収録する）ことにより原音場あるいは再生音場における伝達関数ｈ（ｔ）を求める場合について説明する。
【０００５】
上述した入力信号ｘ（ｔ）、伝達関数ｈ（ｔ）、出力信号ｙ（ｔ）のそれぞれをＺ変換（離散関数のフーリエ変換）することにより、Ｘ（ｚ）、Ｈ（ｚ）、Ｙ（ｚ）が求められるとすると、Ｘ（ｚ）、Ｈ（ｚ）、Ｙ（ｚ）のそれぞれは、
【数１】
【数２】
【数３】
と表される。
【０００６】
このとき、システムの伝達関数Ｈ（ｚ）は既知の入力信号Ｘ（ｚ）と出力信号Ｙ（ｚ）とを用いて
【数４】
を計算することによって求めることができる。
【０００７】
式４に示す式は、具体的に以下の２通りの計算方法を用いることにより、求めることができる。
【０００８】
（１）時間領域で計算する方法
時間領域で計算する方法とは、入力信号Ｘ（ｚ）の逆特性である１／Ｘ（ｚ）を求め、この１／Ｘ（ｚ）の特性を持つフィルタを出力信号Ｙ（ｚ）に畳み込むことにより（畳み込み演算を行うことによって）、伝達関数Ｈ（ｚ）を計算する方法を意味する。この方法を用いて伝達関数ｈ（ｔ）を求めるためには、１／Ｘ（ｚ）の特性を持つフィルタが既知であることが必要とされる。例えば、入力信号が時間伸張パルス（ＴＳＰ：Time stretched pulse）の場合は、入力信号の逆特性を持つフィルタの時間応答が既知となるため、この方法を用いることができる。
【０００９】
（２）周波数領域で計算する方法
周波数領域で計算する方法とは、式４をフーリエ変換し、周波数領域上で周波数成分毎に、
【数５】
を計算して伝達関数Ｈ（ω）の周波数特性を求める方法を意味する。時間応答に該当する伝達関数ｈ（ｔ）は、式５により求められる伝達関数Ｈ（ω）を逆フーリエ変換することにより求めることが可能である。
【００１０】
上述した（１）の方法および（２）の方法のいずれの方法を用いて伝達関数ｈ（ｔ）を求める場合であっても、Ｓ個の音源（スピーカ）より出力される音響信号をＸ（ｚ）（音源１から音源Ｓまでの各音源より出力される信号をＸ１（ｚ）〜ＸＳ（ｚ）とする）、Ｍ個のマイクにより録音される音響信号をＹ（ｚ）（マイク１からマイクＭまでの各マイクにおいて録音される信号をＹ１（ｚ）〜ＹＭ（ｚ）とする）とし、それぞれの音響信号Ｘ（ｚ）および音響信号Ｙ（ｚ）に基づいて求められる伝達関数をＨ（ｚ）（音源ｊより出力される音響信号に対するマイクｉの伝達関数をＨｉｊ（ｚ）とする）とすると、音響信号の入出力関係は、次の式で表現することができる。
【００１１】
Ｙ（ｚ）＝Ｈ（ｚ）Ｘ（ｚ） ・・・式６
【数６】
【数７】
【数８】
【００１２】
入力データとして、ある音源ｊについてのみ入力信号を加えた状態を測定すると、音響信号の入出力関係は、次の式で表現することができる。
【数９】
得られた測定データから音源ｊより出力された音響信号に対する伝達関数Ｈｉｊ（ｚ）は、
【数１０】
により求めることができる。
【００１３】
各音源（ｊ＝１からｊ＝ＳまでのＳ個の音源）に対応する伝達関数Ｈｉｊ（ｚ）を順番に測定することによって、音場全体の伝達関数ｈ（ｔ）を測定することが可能となる。
【先行技術文献】
【特許文献】
【００１４】
【特許文献１】特開２００７−４１１６４号公報
【発明の概要】
【発明が解決しようとする課題】
【００１５】
ところで、従来の伝達関数の測定方法では、１つの音源より出力された音響信号を、複数のマイクで同時に収録する作業に基づいて、スピーカ（音源）毎の伝達関数を求める。このように、スピーカ（音源）毎に順番に伝達関数を求めることによって、全ての音源に対応する伝達関数を求める方法を採用している。
【００１６】
このため、例えば、測定用のマイクを人間の耳位置に取り付けて、複数のスピーカ（音源）とマイクとの間の伝達関数を測定する場合、聴取者（人間）が姿勢を完全に静止させようとしても、姿勢変化を完全に抑制することが現実には困難であった。従って、スピーカ毎に１つずつ音響信号の出力を行っている間に、マイクの位置が変化してしまう恐れが高く、複数の音源に対して完全に同じ測定条件（この場合はマイク位置が変動しないこと）で測定データを得ることが極めて困難であるという問題があった。
【００１７】
さらに、非常に多くのスピーカを使用してスピーカ毎の伝達関数を測定する場合には、全てのスピーカの伝達関数を測定するための時間が、スピーカの設置数に比例して極めて長くなってしまうという問題があった。例えば、１つのスピーカについての測定時間が数秒の時間を要するとしても、スピーカの設置数が６０以上の場合には、総測定時間が数分単位で必要になる。
【００１８】
また、多数のスピーカ（音源）より同時に多数の音響信号を出力させ、出力された音響信号を、複数のマイクで同時に録音することにより、全てのスピーカ（音源）およびマイクに対応する伝達関数を同時に求める方法も考えられるが、従来の伝達関数の測定方法では、複数のスピーカ（音源）で出力され、複数のマイクにおいて収録された音響信号に基づいて、スピーカ（音源）およびマイクに対応する伝達関数成分を、スピーカ（音源）毎に分解する方法が確立されていなかった。このため、多数のスピーカ（音源）より同時に多数の音響信号を出力させ、出力された音響信号を、複数のマイクで同時に録音することにより、全てのスピーカ（音源）およびマイクに対応する伝達関数を同時に高精度に測定する方法を用いることができないという問題があった。
【００１９】
本発明は、上記問題に鑑みて成されたものであり、多数のスピーカ（音源）より同時に音響信号を出力させ、出力された音響信号を、複数のマイクで同時に収録することにより、全てのスピーカ（音源）およびマイクに対応する伝達関数を同時かつ高精度に測定することが可能な音場測定システムを提供することを課題とする。
【課題を解決するための手段】
【００２０】
上記課題を解決するために、本発明に係る音場測定システムは、特定の音響空間に設定される複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録することによって、各スピーカから各マイクに対応する伝達関数を算出するための音場測定システムであって、前記スピーカより出力するための評価音の生成を行う評価音生成手段と、前記マイクを用いて収録された収録音と前記評価音生成手段により生成された評価音とに基づいて、各スピーカから各マイクに対応する伝達関数を算出する伝達関数算出手段とを備え、前記伝達関数算出手段は、前記収録音と前記評価音とに基づいて、時間領域における反復解法を適用することにより、前記伝達関数の算出処理を実行し、前記反復解法の適用において用いる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させる演算方法を用い、前記評価音生成手段は、前記評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第１行列と、該第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成することを特徴とする。
【００２１】
このように、本発明に係る音場測定システムでは、反復解法を用いて伝達関数を算出する場合において、伝達関数算出手段が、相関演算処理と畳み込み演算処理とを行う場合にのみ、周波数領域における演算結果を利用して時間領域における解を求めるので、演算処理の負担軽減と演算処理に利用されるメモリ容量の低減とを図ることが可能となる。
【００２２】
さらに、評価音生成手段が、評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する行列要素を備えた第１行列と、第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成するため、伝達関数算出手段が、複数のスピーカより出力された評価音からなる複数要素の入力行列と、複数のマイクより収録された収録音からなる複数要素の出力行列とに基づいて各伝達関数を求める場合において、相関演算処理と畳み込み演算処理とを行う際に対応する伝達関数の値を一意的に求めることが可能となる。
【００２３】
このため、複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクで評価音を同時に収録する場合であっても、伝達関数算出手段において、各スピーカから各マイクに対応する伝達関数をそれぞれ算出することが可能となる。
【００２４】
また、上述した音場測定システムにおいて、前記伝達関数算出手段は、前記伝達関数の算出処理において、前記伝達関数を求めるための計算式に対して時間窓を構成する対角行列を掛け合わせ、該対角行列の要素を変更することにより、前記複数のスピーカから連続して評価音を出力すると共に前記複数のマイクから連続して収録音の収録を行った場合における任意の時間範囲の伝達関数を算出するものであってもよい。
【００２５】
このように、本発明に係る音場測定システムでは、長時間連続して同時に複数のスピーカから評価音を出力すると共に、スピーカから出力された評価音を連続してマイクで同時に測定した場合であっても、対角行列要素による時間重みを適用して任意の時間範囲における伝達関数を算出することができる。このため、連続的に収録音の収録位置が音響空間において変更される場合であっても、時間変化に応じた適切な伝達関数を求めることが可能となる。
【００２６】
さらに、上述した音場測定システムにおいて、前記評価音生成手段は、前記評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、前記マイクにより収録される収録音の低域周波数帯域におけるＳ／Ｎの改善を図るものであってもよい。
【００２７】
一般的に、特定の音響空間においてスピーカより出力される音をマイクを用いて収録する場合には、マイクで収録された収録音に対してノイズ成分が含まれる場合がある。特に、マイクを使用して音響特性を測定する場合には、低域周波数であるほど、レベルの高い環境ノイズが含まれる傾向がある。
【００２８】
このため、本発明に係る音場測定システムでは、評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、低域周波数帯域に含まれるノイズの出力値を評価音の出力値に対して相対的に低くすることができ、収録音の低域周波数帯域におけるＳ／Ｎの改善を図ることが可能となる。
【００２９】
また、上述した音場測定システムにおいて、前記伝達関数算出手段は、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させて求められた前記相関演算処理または前記畳み込み演算処理の演算結果のうち、前記時間領域のみで前記相関演算処理または前記畳み込み演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、前記時間領域における前記伝達関数の算出処理を行うものであってもよい。
【００３０】
このように、上述した音場測定システムにおいて、伝達関数算出手段が、周波数領域から時間領域に逆変換させた相関演算処理または前記畳み込み演算処理の演算結果のうち、時間領域のみで演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、伝達関数の算出処理を実行することによって、最終的に算出される伝達関数の値に誤差等が含まれてしまうことを効果的に防止することが可能となる。
【００３１】
また、上述した音場測定システムにおいて、前記誤差が生じない範囲は、前記評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定されるものであってもよい。
【００３２】
このように、誤差が生じない範囲が、評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定される場合には、算出された伝達関数の性能（品質）を十分に確保しつつ、最終的に算出される伝達関数の値に誤差等が含まれてしまうことを抑制することが可能となる。
【発明の効果】
【００３３】
本発明に係る音場測定システムによれば、評価音生成手段が、評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第１行列と、該第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成するため、伝達関数算出手段が、複数のスピーカより出力された評価音からなる複数要素の入力行列と、複数のマイクより収録された収録音からなる複数要素の出力行列とに基づいて各伝達関数を求める場合において、相関演算処理と畳み込み演算処理とを行う際に対応する伝達関数の値を一意的に求めることが可能となる。
【００３４】
このため、複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録する場合であっても、伝達関数算出手段において、各スピーカから各マイクに対応する伝達関数をそれぞれ算出することが可能となる。
【図面の簡単な説明】
【００３５】
【図１】実施の形態１に係る音場測定システムの概略構成を示すブロック図である。
【図２】実施の形態１に係る音場再現システムにおける、評価音の信号（入力信号）Ｘ（ｚ）と、収録音の信号（測定信号）Ｙ（ｚ）と、伝達関数Ｈ（ｚ）との関係を示した図である。
【図３】実施の形態１に係る伝達関数算出部において行われる、Ｄ＝ＡＱを解くための数学計算方式の処理内容を示したフローチャートである。
【図４】実施の形態１に係る伝達関数算出部において行われる、方程式を組み立てる計算処理内容を示したフローチャートである。
【図５】実施の形態１に係る伝達関数算出部において行われる、方程式を解く処理内容を示したフローチャートの前半図である。
【図６】実施の形態１に係る伝達関数算出部において行われる、方程式を解く処理内容を示したフローチャートの後半図である。
【図７】実施の形態１に係る伝達関数算出部において行われる、行列ＲとベクトルＶとの演算処理内容を示したフローチャートである。
【図８】実施の形態２に係る伝達関数算出部において行われる、行列ＡとベクトルＶとの積の計算処理内容を示したフローチャートである。
【図９】実施の形態２に係る伝達関数算出部において行われる、転置行列ＡｔとベクトルＶとの積の計算処理内容を示したフローチャートである。
【図１０】実施の形態２に係る伝達関数算出部において行われる、連続測定した測定データの一部を用いる伝達関数の計算処理内容を示したフローチャートである。
【図１１】実施の形態２に係る伝達関数算出部において行われるＡｔ（Ｗ２Ｄ）の計算処理内容を示したフローチャートの後半図である。
【図１２】実施の形態２に係る伝達関数算出部において行われる、Ａｔ（Ｗ２（Ａσμ））の計算処理内容を示したフローチャートである。
【図１３】実施の形態３に係る音場測定システムの音場空間においてノイズが存在する場合における測定データＹ（ｚ）を、入力信号Ｘ（ｚ）と、伝達関数Ｈ（ｚ）と、ノイズ成分Ｎ（ｚ）との関係を示したブロック図である。
【図１４】実施の形態３に係る音場測定システムの入力信号の周波数特性と、ノイズの周波数特性とを示した図であり、（ａ）は、入力信号の周波数特性が平坦な場合示しており、（ｂ）は、入力信号の周波数特性を変えることにより、低域周波数成分の信号パワーが大きくなるようにした場合を示している。
【発明を実施するための形態】
【００３６】
以下、本発明に係る音場測定システムについて、図面を用いて詳細に説明を行う。
【００３７】
［実施の形態１］
図１は、本実施の形態１に係る音場測定システムの概略構成を示した図である。音場測定システム１は、再現を望む音響空間（原音場）あるいは、原音場の音響環境の再現を行う音響空間（再生音場）において設置されるシステムである。
【００３８】
音場測定システム１は、図１に示すように、音響空間２に設置されるスピーカ部３と、スピーカ部３に対して音声出力を行う音声出力装置４と、音声出力装置４に対して評価音を出力するコンピュータ５と、音響空間２に設置されるマイクアレイ７と、マイクアレイ７により収録された収録音の入力を行う音声入力装置８とによって概略構成されている。
【００３９】
スピーカ部３は、図１に示すように、マイクアレイ７を中心として等距離を保つように円形に配置された複数（例えばＳ個）のスピーカ９で構成されている。このスピーカ９は、音声出力装置４より出力された評価音（音響信号）を出力することが可能となっている。
【００４０】
なお、図１に示すスピーカ部３の配置は、再生音場における一般的なスピーカの配置状態として実験室を用いた場合を一例として示している。例えば、本実施の形態１に係る音場測定システム１を原音場に設置する場合には、その原音場の状況に応じてスピーカの設置位置およびスピーカの設置個数が決定される。例えば、原音場として車両の室内空間を想定する場合には、車室の四隅に対して４個のスピーカが均等に設置される場合などが考えられる。
【００４１】
音声出力装置４は、コンピュータ５より出力された評価音を、スピーカ部３のＳ個のスピーカ９に対して同時に出力する役割を有している。ここで評価音とは、マイクアレイ７および音声入力装置８を介して収録された収録音に基づいて音響空間２の伝達関数（音響特性）を算出するために用いられる出力音であり、一般的には、収録された収録音をインパルス応答として測定することに適した音（例えば、Ｍ系列信号や時間伸張パルス（ＴＳＰ：Time Stretched Pulse））が該当する。なお、本実施の形態において用いられる評価音については、後述する。
【００４２】
マイクアレイ７は、円形に配設された複数（例えばＭ個）のマイク１０の一群により構成されている。マイクアレイ７は、各スピーカ９から等距離となる音響空間２の中心位置に設置されている。音声入力装置８は、マイクアレイ７によって測定されたスピーカ部３の出力音をコンピュータ５に出力する役割を有している。
【００４３】
コンピュータ５は、上述した評価音を生成すると共に、マイクアレイ７により収録された収録音に基づいて音響空間２における伝達関数（伝達関数マトリックス、音響特性）を算出する機能を有している。
【００４４】
コンピュータ５は、評価音生成部（評価音生成手段）１２と、収録音記録部１３と、伝達関数算出部（伝達関数算出手段）１４と、伝達関数記録部１５と、制御部１６とによって概略構成されている。評価音生成部１２は、上述したインパルス応答を収録するための評価音を生成する機能を有している。評価音生成部１２は、制御部１６の制御命令に応じて評価音を生成し、制御部１６を介して音声出力装置４に評価音を出力する。
【００４５】
収録音記録部１３は、マイクアレイ７および音声入力装置８を介して収録されたマイク毎の収録音を記録する機能を有している。伝達関数算出部１４は、音声出力装置４により出力された評価音と、収録音記録部１３に記録された収録音とに基づいて、スピーカ部３のＳ個のスピーカ９とマイクアレイ７のＭ個のマイク１０との間におけるインパルス応答をそれぞれ求めることにより、各スピーカと各マイクとに対応した伝達関数を算出する。なお、各スピーカと各マイクとに対応した伝達関数の要素の集合により、音響空間２における全体の音響特性が求められる。
【００４６】
伝達関数算出部１４は、伝達関数の算出処理を行うＣＰＵ（Central Processing Unit）と、伝達関数を算出するためのプログラム等が記録されるＲＯＭ（Read Only Memory）、ワークエリア等に使用されるＲＡＭ（Random Access Memory）等の一般的な演算処理器によって構成されている。
【００４７】
伝達関数記録部１５は、伝達関数算出部１４により求められた伝達関数を記録する機能を有している。伝達関数記録部１５に記録される伝達関数は、各スピーカと各マイクとに対応付けてそれぞれ記録される。
【００４８】
制御部１６は、音響空間２における伝達関数（音響特性）を求めるために必要な処理を行う役割を有している。制御部１６は、上述したように評価音生成部１２に対して評価音を生成させると共に、生成された評価音を音声出力装置４に出力する役割を有している。また、制御部１６は、音声入力装置８より入力される収録音を収録音記録部１３に記録させると共に、伝達関数算出部１４に伝達関数を算出させて、伝達関数記録部１５に記録させる役割を有している。
【００４９】
なお、本実施の形態１に係る音場測定システム１では、評価音生成部１２において生成された評価音が、スピーカ部３の全て（Ｓ個）のスピーカ９より出力される。マイクアレイ７では、全て（Ｓ個）のスピーカより出力された出力音を全て（Ｍ個）のマイク１０で同時に収録し、収録された収録音は、収録音記録部１３においてマイク毎に記録される。
【００５０】
次に、音声出力装置４より出力される評価音と、収録音記録部１３に記録される収録音とに基づいて、伝達関数算出部１４で音響空間２の伝達関数を算出する方法について説明する。
【００５１】
まず、スピーカ部３より出力される評価音を測定時の既知の入力信号波形（評価音信号、入力信号）をｘｊ（ｔ）（ｊはスピーカを示し、ｊ＝１，２，・・・，Ｓ）とし、マイクアレイ７において収録される収録音を既知の出力信号波形（収録音信号、出力信号）をｙｉ（ｔ）（ｉはマイクを示し、ｉ＝１，２，・・・，Ｍ）とし、それぞれをｚ変換で表現したものを、Ｘ（ｚ）、Ｙ（ｚ）とする。また、入力信号波形ｘｊ（ｔ）および出力信号波形ｙｉ（ｔ）により求められる伝達関数をｈｉｊ（ｔ）とし、ｚ変換で表現されたＸ（ｚ）およびＹ（ｚ）に基づいて求められる未知の伝達関数マトリックスをＨ（ｚ）とする。
【００５２】
音場測定システム１における測定時間を有限とし、入力信号として有限の時間長さＮｘの信号を用いると、入力信号波形ｘｊ（ｔ）の値が０でない時間範囲は、０≦ｔ≦Ｎｘ−１となる。
【００５３】
また、伝達関数ｈｉｊ（ｔ）の時間応答を有限の応答時間長さＬで近似できると仮定し、ｈｉｊ（ｔ）の値が０でない時間範囲を、０≦ｔ≦Ｌ−１であるとする。このとき、出力信号波形ｙｉ（ｔ）の値が入力信号波形ｘｊ（ｔ）に依存する時間応答も有限となり、その時間は、０≦ｔ≦Ｎｘ＋Ｌ−２（入力信号波形ｘｊ（ｔ）の０でない時間範囲（Ｎｘ−１）とｈｉｊ（ｔ）の値が０でない時間範囲（Ｌ−１）とを加えた時間範囲）の範囲となる。
【００５４】
従って、音場測定システム１において測定する出力信号波形ｙｉ（ｔ）の測定サンプル数をＮｙとすると、測定サンプル数Ｎｙは、出力信号波形ｙｉ（ｔ）の応答時間に基づいてＮｘ＋Ｌ−１となるように選べばよくなる。測定サンプル数ＮｙとしてＮｙ＝Ｎｘ＋Ｌ−１を選ぶことにより、出力信号波形ｙｉ（ｔ）の応答時間に対応する値が０以上となる測定サンプルを求めることが可能となる。
【００５５】
複数のスピーカ９と複数のマイク１０とを備えた音場測定システム１における測定系の入出力関係は、以下の関係式で表すことができる。
Ｙ（ｚ）＝Ｈ（ｚ）Ｘ（ｚ） ・・・式１２
【００５６】
但し、Ｘ（ｚ）は、
【数１１】
で示され、Ｘ（ｚ）の要素Ｘｊ（ｚ）（ｊ＝１，２，・・・，Ｓ）は、
【数１２】
を表している。
【００５７】
また、Ｙ（ｚ）は、
【数１３】
で示され、Ｙ（ｚ）の要素Ｙｉ（ｚ）（ｉ＝１，２，・・・，Ｍ）は、
【数１４】
を表している。
【００５８】
また、Ｈ（ｚ）は、
【数１５】
で示され、Ｈ（ｚ）の要素Ｈij（Ｚ）（ｉ＝１，２，・・，Ｍ、ｊ＝１，２，・・，Ｓ）は、
【数１６】
を表している。
【００５９】
式１２の両辺の転置をとると、
Ｙｔ（ｚ）＝Ｘｔ（ｚ）Ｈｔ（ｚ） ・・・式１９
となり、Ｘｔ（ｚ）は、
Ｘｔ（ｚ）＝（Ｘ１（ｚ） Ｘ２（ｚ） ・・・ ＸＳ（ｚ））
・・・ 式２０
Ｙｔ（ｚ）は、
Ｙｔ（ｚ）＝（Ｙ１（ｚ） Ｙ２（ｚ） ・・・ ＹＭ（ｚ））
・・・ 式２１
Ｈｔ（ｚ）は、
【数１７】
となる。
【００６０】
ここで、ｚ変換同士の積は、畳み込み演算マトリックスとベクトルの積として表すことができ、式１９は、等価な連立一次方程式である
Ｄ＝ＡＱ ・・・式２３
で表すことができる。但し、
【数１８】
【数１９】
【数２０】
【数２１】
で表される。
従って、伝達関数Ｈ（ｚ）を求めることは、上述した連立一次方程式Ｄ＝ＡＱを解く問題に帰着することができる。
【００６１】
音場測定システム１において、図２に示すように、測定対象である伝達関数Ｈ（ｚ）が線形の伝達関数モデルである限り、Ｄ＝ＡＱを満たす解は、必ず存在する（但し、等号が成立するのは、測定データ中にノイズがない場合となる）。このため、Ｄ＝ＡＱの解を求めることによって伝達関数Ｈ（ｚ）を算出できるようにするためには、Ｄ＝ＡＱを満たす解が１つだけ存在するように方程式を構築する必要がある。
【００６２】
具体的に、Ｄ＝ＡＱを満たす解が１つだけ存在するように方程式を構築するためには、ＡｔＡの逆行列が存在するように入力信号波形ｘ（ｔ）の組を与えることが必要とされる。つまり、入力信号波形（評価音信号、入力信号）ｘｊ（ｔ）の畳み込み行列（式２６）を行列要素とし、各スピーカに対応する行列要素（式２６）を備えた行列Ａ（式２５：第１行列）と、この行列Ａ（第１行列）の転置行列を成す行列Ａｔ（第２行列）との積であるＡｔＡの逆行列が存在し得るように入力信号波形（評価信号、評価音）ｘｊ（ｔ）が設定されればよい。
【００６３】
ＡｔＡの逆行列が存在する場合、Ｑは、Ｄ＝ＡＱの最小二乗解として求められ、次の連立一次方程式
（ＡｔＤ）＝（ＡｔＡ）Ｑ ・・・式２８
の解として、次式で求めることができる。
Ｑ＝（ＡｔＡ）−１（ＡｔＤ） ・・・式２９
【００６４】
上述した計算原理を用いることにより、複数の音源を同時に再生した場合の測定データを用いて伝達関数マトリックスのすべての成分を正確に計算することが可能となる。このため、１回の測定ですべてのスピーカ・マイク間の伝達関数を求めることが可能になる。
【００６５】
具体的な数値解を得るためには、式２３および式２８を具体的に解く数値計算アルゴリズムおよび入力信号の与え方が重要になる。一般に式２３および式２８で示される方程式は、計算規模（演算量）が極めて大きく、汎用の行列解法を用いる場合には、極めて処理能力の高い計算環境を用意する必要があるため、方程式の演算を実行することが困難である。
【００６６】
具体的な数値を挙げて説明すると、式２８の係数行列ＡｔＡの計算規模は、音源数（スピーカ９数）をＳ、伝達関数の時間応答長さをＬとした場合、Ｓ×Ｌ行、Ｓ×Ｌ列の行列を計算するための計算規模が必要となる。伝達関数を求める場合、時間応答長さＬは、通常、数千〜数万のオーダーで設定されるため、仮に時間応答長さＬを２１５＝３２７６８に設定し、データ精度を単精度実数（４ｂｙｔｅ）とすると、係数行列の演算処理に用いられる記憶容量だけでも、Ｓ×Ｓ×４Ｇｂｙｔｅ程度のメモリ容量が必要となる。従って、１６音源を一度に測定したい場合には、１０００Ｇｂｙｔｅ程度のメモリ容量が係数記憶のためだけに必要になる。さらに、このような大規模の方程式を解くための演算コストや計算時間も膨大なものとなってしまう。
【００６７】
そのため、本実施の形態１に係る伝達関数算出部１４では、畳み込み問題の性質に特化して効率的な計算を行う計算アルゴリズムを用いている。図３〜図６は、本実施の形態１における計算アルゴリズムを示したフローチャートである。
【００６８】
図３は、Ｄ＝ＡＱを解くための数学計算方式の概要を示した図である。図３に示すように、まず、伝達関数算出部１４は、方程式Ｄ＝ＡＱを最小二乗法で解くために、Ｐ＝ＲＱに基づいて
Ｐ＝ＡｔＡ
Ｒ＝ＡｔＤ
とおいて、ＰおよびＲの計算を行い（ステップＳ．１）、この計算式に基づいて連立一次方程式Ｐ＝ＲＱを解くことにより、時間領域における伝達関数マトリックス（行列Ｑ）（Ｑ＝Ｒ−１Ｐ）を求める（ステップＳ．２）
【００６９】
以下、伝達関数算出部１４が、時間領域において伝達関数マトリックス（行列Ｑ）を求める処理を説明する。伝達関数算出部１４は、時間領域において
Ｑ＝（ＡｔＡ）−１（ＡｔＤ） ・・・式２９
の最適解を求めるために、反復解法を用いるものとし、特に本実施の形態１では反復解法の解法例として共役勾配法を用いて演算を行うものとする。
【００７０】
伝達関数マトリックス（行列Ｑ）を算出するためには、上述したように、
Ｑ＝（ＡｔＡ）−１（ＡｔＤ） ・・・式２９
を解く必要がある。ここで、反復解法を用いて伝達関数マトリックス（行列Ｑ）を算出する場合、反復計算の処理において必要とされる行列形式は、規則的な畳み込み行列になるという特徴がある。このため、本実施の形態１に係る伝達関数算出部１４では、この畳み込み行列を利用した特徴的な演算処理を行うことによって、計算過程に必要なメモリ量をフィルタ段数Ｌに比例するオーダーまで削減することが可能となる。このように、連立一次方程式を解く大きな枠組みとして、反復法（共役勾配法）を用いることによって、係数行列の規則性を保ち、メモリの使用量を抑えたままでの計算を可能にすることができる。
【００７１】
また、本実施の形態１に係る計算アルゴリズムでは、係数行列ＡｔＡが相関行列となるため、方程式の組み立て処理の段階で相関行列の規則性を利用することにより、計算に必要な記憶容量をＬ×Ｌに比例するオーダーから、Ｌに比例するオーダーへと大幅に削減することができる。
【００７２】
従って、方程式の組み立て処理と反復計算の処理とに用いられる主要な行列演算を、高速フーリエ変換（ＦＦＴ）を利用した相関演算処理と畳み込み演算処理とで演算することによって、計算量を大幅に削減することができ、音場測定システムの規模が大きくなっても、一般的なコンピュータを用いて演算処理を行うことが可能となる。
【００７３】
なお、上述する計算アルゴリズムにおいて入力信号波形として入力される評価音（入力信号）には、連立一次方程式
（ＡｔＤ）＝（ＡｔＡ）Ｑ ・・・式２８
の係数行列ＡｔＡの逆行列が存在するような信号系列を選ぶことが必要である。
【００７４】
このように、係数行列ＡｔＡの逆行列が存在するような信号系列を選ぶことによって、原理上、計算アルゴリズムにおいて解の精度が劣化することを防止できる。なお、有限桁数の数値計算で有限の計算時間で解を求める関係上、係数行列が単位行列に近くなるように入力信号を選ぶことが望ましい。係数行列が単位行列に近いほど、連立一次方程式を反復解法で解く場合の収束回数が少なくなるので、有限の反復回数で計算処理を打ち切った場合に得られる解の精度を向上させることができる。
【００７５】
このような入力信号として、正規分布や一様分布などの自己相関特性が鋭くて（自己相関波形が時刻０でのみ大きな値を持ち、その他の時刻では０に近い値を持つもの）、相互相関の低い信号を用いることが好ましい。従って、評価音生成部１２では、このような条件を満たし得る入力信号（評価音）の生成を行う。
【００７６】
図４は、図３に示したステップＳ．１の処理内容である方程式を組み立てるための計算処理を示したフローチャートである。また、図５および図６は、図３に示したステップＳ．２の処理内容である方程式を解く計算処理を示したフローチャートであり、共役勾配法を使用して伝達関数マトリックス（行列Ｑ）を算出する処理を示している。
【００７７】
フローチャートに従って時間領域だけで演算処理を行う場合には、次述する（１）〜（４）の演算処理において処理負担が増大するという問題があった。
（１）行列Ｒ：Ｒ＝ＡtＡの演算処理
（２）行列Ｐ：Ｐ＝ＡtＤの演算処理
（３）初期値計算におけるＲとベクトルとの積Ｒｑ0の演算処理
（４）反復計算におけるＲとベクトルとの積Ｒσμの演算処理
また、行列Ｒの演算に要する記憶容量も増大するという問題があった。
【００７８】
本実施の形態１に係る計算アルゴリズムでは、上述したように、（１）〜（４）の演算を行う場合に、高速フーリエ変換（ＦＦＴ）を用いて周波数領域における演算処理を行い、演算処理された演算結果を逆高速フーリエ変換（ＩＦＦＴ）によって時間領域の演算結果に逆変換する（時間領域の演算結果に戻す）ことによって、演算処理負担の軽減を図っている。
【００７９】
まず、図４に示すフローチャートに基づいて、初期演算として行列Ｒと行列Ｐとを計算する場合について説明する。
【００８０】
行列Ｒは、上述したようにＲ＝ＡtＡで示される。
ここで行列Ｒの部分行列Ｒi,jは、相関行列であって規則正しいパターンを持ち、斜め方向の要素の値が同一の値になるという特徴を有している。この部分行列Ｒi,jでは、２Ｌ−１個の異なる値を備えている。このため、係数行列Ｒの演算を行う場合には、重複する値をメモリ等の記憶領域に記憶せず、異なる値のみをベクトルとして記憶することによって、伝達関数算出部１４の演算処理で使用されるメモリ使用量をＬのオーダーで削減することが可能となる。
【００８１】
また、行列Ｒi,jの各要素ｒi,jは、伝達関数同士の相関演算の形を含んでいる。このため、相関演算処理を行う場合、高速フーリエ変換（ＦＦＴ）を用いてｒi,jの値を周波数領域の値に変換して解を求めた後に、求められた解を逆フーリエ変換（ＩＦＦＴ）により逆変換して、時間領域におけるｒi,jの相関演算処理の解を求める。このように周波数領域において演算処理を行うことによって、処理負担を軽減させることができる。
【００８２】
また、行列Ｐは、上述したようにＰ＝ＡtＤで示され、行列Ｐi,jの各要素ｐi,jも、行列Ｒと同様に、伝達関数同士の相関演算の形を含んでいるため、ｐi,jに対して高速フーリエ変換を適用し、ｐi,jの値を周波数領域の値に変換して解を求め、求められた解を逆フーリエ変換して時間領域の解に逆変換することにより値を求めることができる。このように周波数領域の解を時間領域の解に逆変換する処理を用いることにより、演算処理の負担を軽減させることができる。
【００８３】
ここで、伝達関数算出部１４は、ＲおよびＰの高速フーリエ変換を行う場合における段数の数Ｎcorrとして、Ｎcorr≧Ｎｙ＋Ｌ−１を満たす値、例えば、Ｎcorr＝Ｎｙ＋Ｌを設定する（ステップＳ．２１）。
【００８４】
そして、伝達関数算出部１４は、入力信号波形であるｘｊ（ｔ）（但しｊ＝１，２，・・・，Ｓ）および出力信号波形であるｙｉ（ｔ）（但しｉ＝１，２，・・・，Ｍ）に対してＮcorr点の高速フーリエ変換を適用する（ステップＳ．２２）。
【００８５】
なお、ｘｊ（ｔ）およびｙｉ（ｔ）における時系列のデータ数がＮcorrに満たない場合、伝達関数算出部１４は、各データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加してデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。このように調整を行うことにより、共役勾配法を用いて求められる伝達関数に誤差が生じてしまうことを防止することができる。
【００８６】
高速フーリエ変換を適用することにより、
Ａ（ω）≡Ｘｔ（ω）＝（Ｘ１（ω） Ｘ２（ω） …ＸＳ（ω））
・・・式３０
Ｄ（ω）≡Ｙｔ（ω）＝（Ｙ１（ω） Ｙ２（ω） …ＹＭ（ω））
・・・式３１
となる。
【００８７】
次に、伝達関数算出部１４は、求められたＡ（ω）およびＤ（ω）を用いて、Ｒ（ω）およびＰ（ω）の周波数領域における相関演算を実施する（ステップＳ．２３）。
式３０および式３１より、Ｒ（ω）およびＰ（ω）は、
Ｒ（ω）＝Ａ＊（ω）Ａ（ω） ・・・式３２
Ｐ（ω）＝Ａ＊（ω）Ｄ（ω） ・・・式３３
で示すことができる。ここで、Ａ＊（ω）は、Ａ（ω）の共役転置行列を示しており、
Ａ＊（ω）は、
【数２２】
で表される。
【００８８】
また、Ｒ（ω）の要素ｒi,j（ω）は、
ｒi,j（ω）＝conj（Ｘｉ（ω））Ｘｊ（ω）
で表すことができ、Ｐ（ω）の要素ｐi,j（ω）は、
ｐi,j（ω）＝conj（Ｘｉ（ω））Ｙｊ（ω）
で表すことができる。
なお、conj（Ｘ）はＸの複素共役を示し、conj（Ｙ）はＹの複素共役を示している。
【００８９】
このように、フーリエ変換されたｒi,j（ω）およびｐi,j（ω）の値を周波数成分の掛け算を用いて求めた後に、求めた解を逆フーリエ変換して時間領域の時系列データに変換することによって、時間領域において演算を行う場合に比べて迅速かつ容易に解を求めることができる。
【００９０】
伝達関数算出部１４は、求められたｒi,j（ω）の値を逆フーリエ変換（ＩＦＦＴ）することにより、時系列データであるｒi,j（ｔ）に変換する（ステップＳ．２４）。この場合、伝達関数算出部１４は、逆フーリエ変換によって求められた解が正確な値となるように（時間領域のみの演算により求められる解とに誤差が生じないように）、−Ｌ＋１≦ｔ≦Ｌ−１の範囲のデータのみを用いて解を求める。つまり、誤差の生じない範囲は、出力音の残響時間等に応じて予め設定される有限の応答時間長さＬに応じて決定される。
【００９１】
伝達関数算出部１４は、−Ｌ＋１≦ｔ≦Ｌ−１の範囲のデータのうち、ｔ＜０のときに、ｒi,j（ｔ＋Ｎcorr）を負の時刻のｒi,j（ｔ）として所定の記憶手段に記憶し、ｔ≧０のときに、ｒi,j（ｔ）をそのまま記憶する。
【００９２】
一方で、伝達関数算出部１４は、求められたｐi,j（ω）の値を逆フーリエ変換（ＩＦＦＴ）することにより、時系列データであるｐi,j（ｔ）に変換する（ステップＳ．２５）。この場合、伝達関数算出部１４は、逆フーリエ変換によって求められた解が正確な値となるように、０≦ｔ≦Ｌ−１の範囲のデータのみを用いて解を求め、記憶する。
【００９３】
次に、図５および図６に示すフローチャートに基づいて、伝達関数算出部１４の具体的な伝達関数の演算処理について説明する。
【００９４】
伝達関数算出部１４は、伝達関数の具体的な演算処理においても、行列Ｒi,jの各要素ｒi,jに含まれる伝達関数同士の相関演算の形を利用する。伝達関数算出部１４は、相関演算処理を行う場合において、高速フーリエ変換（ＦＦＴ）を用いてｒi,jの値を周波数領域の値に変換して解を求めた後に、求められたｒi,jを周波数領域において演算処理し、演算処理された演算結果を逆フーリエ変換（ＩＦＦＴ）を用いて逆変換して、時間領域におけるｒi,jの相関演算処理の解を求める。このように周波数領域において演算処理を行うことによって、処理負担を軽減させることができる。
【００９５】
まず、伝達関数算出部１４は、Ｒに対して高速フーリエ変換を行う場合の段数をＮＣＧとして、ＮＣＧ≧２Ｌ−１を満たす値、例えば、ＮＣＧ＝２Ｌを設定する（ステップＳ．３１）。
【００９６】
そして、伝達関数算出部１４は、図４のステップＳ．２４において求められたｒi,j（ｔ）を並べ直したベクトルｒ’i,j（ｔ） ０≦ｔ≦ＮＣＧ−１を定義する（ステップＳ．３２）。具体的にｒ’i,j（ｔ）は、ｔの値に応じて、下記のように示される。
【数２３】
【００９７】
伝達関数算出部１４は、ｒ’i,j（ｔ）のＮＣＧ点における高速フーリエ変換（ＦＦＴ）を行うことにより、周波数領域におけるＲＣＧ（ω）のｉｊ成分を計算する（ステップＳ．３３）。
【００９８】
次に、伝達関数算出部１４は、初期値としてｊに１を代入し、μに０を代入する（ステップＳ．３４）、ここでｊは、音響空間２におけるスピーカ９の数を示している。またμは、行列Ｑのｊ列目の値を算出するために実行された演算回数（反復回数）を示すパラメータであり、具体的には、ステップＳ．３６〜ステップＳ．３８の処理を繰り返し実行した回数を示している。
【００９９】
次に伝達関数算出部１４は、初期設定を行う（ステップＳ．３５）。初期設定では、Ｐ，Ｑのｊ列目の値をｐ，ｑとおき、適当な初期値ｑ0を用いて、
【数２４】
の初期値であるε0の値
ε0＝ｐ−Ｒｑ0 ・・・式３７
を求める。なお、εμは、反復回数μでの暫定的な解をｑμとしたときの、誤差ベクトルを示したものである。
【０１００】
ここで。ε0の算出に用いられる行列Ｒと任意のベクトルＶとの演算処理（例えば、Ｒｑ0の演算）について、図７に示すフローチャートに沿って説明する。
【０１０１】
係数行列Ｒに対する積算対象となる任意のベクトルＶを
【数２５】
とする。ベクトルＶは、長さＬのサブベクトルをＳ個並べたベクトルとして定義される（ステップＳ．５１）。
【０１０２】
このようにして示されるベクトルＶのそれぞれを、ステップＳ．３１において設定したＮＣＧ点で高速フーリエ変換（ＦＦＴ）すると、
【数２６】
を求めることができる（ステップＳ．５２）。
【０１０３】
このようにして求められたＶ（ω）とステップＳ．３３において求められたＲＣＧ（ω）との積を、周波数領域において算出すると、
【数２７】
を求めることができる（ステップＳ．５３）。このようにして算出されたＲＣＧ（ω）Ｖ（ω）の算出結果を時間領域に変換（逆フーリエ変換）することにより
【数２８】
と示すことができ（ステップＳ．５４）、この（ＲＶ）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＲＶを求めることができる（ステップＳ．５５）。なお、０≦ｔ≦Ｌ−１の範囲のデータのみを用いて計算を行うのは、逆フーリエ変換によって求められた解が正確な値となるように（時間領域のみの演算により求められる解とに誤差が生じないように）するためのである。
【０１０４】
例えば、Ｒｑ0の演算を行う場合には、積算対象となるベクトルｑ0をベクトルＶと置き換えて、図７に示すフローチャートの計算処理を行うことにより、周波数領域における演算処理を用いてＲｑ0の演算を行うことができる。
【０１０５】
具体的に伝達関数算出部１４は、ｑ0をＮＣＧ点で高速フーリエ変換（ＦＦＴ）してｑ0（ω）を計算する。伝達関数算出部１４では、ステップＳ．３３において計算されたＲＣＧ（ω）と、ｑ0（ω）とを周波数領域において積算することにより、ＲＣＧ（ω）ｑ0（ω）を算出し、算出されたＲＣＧ（ω）ｑ0（ω）の算出結果を時間領域に変換（逆フーリエ変換）することにより（Ｒｑ0）’である
（Ｒｑ0）’＝ＩＦＦＴ（ＲＣＧ（ω）ｑ0（ω）） ・・・式４２
を求める。そして、伝達関数算出部１４は、求められた（Ｒｑ0）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＲｑ0を求めることができる。このように、伝達関数算出部１４において、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部１４における処理負担を軽減させることができる。
【０１０６】
また、伝達関数算出部１４は、
【数２９】
の初期値であるσ0の値
【数３０】
の値を算出する（ステップＳ．３５）。ここで、σμは、後述する演算処理におけるＲσμの計算に用いられる。
【０１０７】
さらに、伝達関数算出部１４は、共役勾配法の収束判定に用いられるパラメータをηとして設定する（ステップＳ．３５）。ηは収束を判断するための閾値パラメータとしての定数であって、このパラメータに基づいて誤差が収束判定条件を満たす値となった場合（次式４５を満たす場合）、伝達関数算出部１４は、解が収束したと判断して反復計算を終了する。
【０１０８】
そして、伝達関数算出部１４は、収束条件の判断処理（ステップＳ．３６）を行う。収束条件の判断式は、
【数３１】
で表される。伝達関数算出部１４は、この判断式に基づいて解の収束を判断する。
【０１０９】
収束条件を満たさなかった場合（ステップＳ．３６においてＮｏの場合）、伝達関数算出部１４は、Ｒσμの計算を行う（ステップＳ．３７）
【０１１０】
Ｒσμの計算を行う場合にも、図７に示すフローチャートを用いて説明したように、行列Ｒとベクトルとの積を求める計算に該当する畳み込み演算処理を用いて計算を行うことが可能である。
【０１１１】
具体的に伝達関数算出部１４は、σμをＮＣＧ点で高速フーリエ変換（ＦＦＴ）してσμ（ω）を計算する。伝達関数算出部１４では、ステップＳ．３３において計算されたＲＣＧ（ω）と、σμ（ω）とを周波数領域において積算することにより、ＲＣＧ（ω）σμ（ω）を算出し、算出されたＲＣＧ（ω）σμ（ω）の算出結果を時間領域に変換（逆フーリエ変換）することにより、
（Ｒσμ）’＝ＩＦＦＴ（ＲＣＧ（ω）σμ（ω）） ・・・式４６
を求める。
【０１１２】
この場合においても、伝達関数算出部１４は、時間領域において解に誤差が生じないように、０≦ｔ≦Ｌ−１の範囲のデータのみを用いて解を求める。従って、伝達関数算出部１４は、求められた（Ｒσμ）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＲｑ0を求める。このように、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部１４における処理負担を軽減させることができる。
【０１１３】
そして、伝達関数算出部１４は、求められたＲσμの値と、εμ、σμとを用いて、
【数３２】
【数３３】
【数３４】
を算出し、さらに、μに１を加える（μ＝μ＋１）処理を行う（ステップＳ．３８）。
【０１１４】
そして、伝達関数算出部１４は、ステップＳ．３８において算出されたεμとＰとを用いて、収束条件の判断処理（ステップＳ．３６）を繰り返し実行する。
【０１１５】
収束条件を満たす場合（ステップＳ．３６においてＹｅｓの場合）、伝達関数算出部１４は、ステップＳ．３８において算出されたｑμの値を、行列Ｑにおけるｊ列目の解（伝達関数）とする（ステップＳ．３９）。
【０１１６】
その後、伝達関数算出部１４は、全てのｊについて、つまり全てのスピーカ９に対応する各マイク１０の伝達関数マトリックス（行列Ｑ）における全ての伝達関数を求めたか否かを判断し（ステップＳ．４０）、全ての列において伝達関数を求めていない場合（ステップＳ．４０においてＮｏの場合）には、ｊに１を加え、μに０を代入した後に（ステップＳ．４１）、処理をステップＳ．３５に移行して上述した処理を繰り返し実行する。
【０１１７】
全ての列において伝達関数を求めた場合（ステップＳ．４０においてＹｅｓの場合）、伝達関数算出部１４は、全ての列の伝達関数の値を求めたものと判断して、方程式を解く計算処理（伝達関数マトリックス（行列Ｑ）の算出処理）を終了する。
【０１１８】
このようにして伝達関数算出部１４が、複数のスピーカ９より出力される入力信号波形に関するデータｘｊ（ｔ）と、複数のマイクにより収録される出力信号波形に関するデータｙｉ（ｔ）とに基づいて伝達関数の算出を行うことによって、同時に複数のスピーカ９から評価音を出力すると共に、スピーカ９から出力された評価音をマイク１０で同時に測定することにより、各スピーカ９から各マイク１０に対応する伝達関数を同時に算出することが可能となる。このため、従来のように、評価音をそれぞれのスピーカ９より順番に出力することにより、出力先となるスピーカ９と収録先となるマイク１０とを対応させて伝達関数を算出する必要がなくなり、伝達関数の測定時間の短縮化を図ることが可能となる。
【０１１９】
また、一度に全てのスピーカから評価音を出力させて、同時に全てのマイクで収録を行うことにより伝達関数の算出を行うため、従来のように集音時におけるマイク位置の変動などによる影響を受けにくくなり、測定精度の向上を図ることが可能となる。
【０１２０】
［実施の形態２］
次に、実施の形態２に係る音場測定システムについて説明を行う。
【０１２１】
実施の形態１に係る音場測定システム１では、複数のスピーカ９から同時に評価音を出力すると共に、複数のマイク１０から同時に評価音の収録を行うことにより、一度の測定でそれぞれのマイク１０とスピーカ９とに対応させた伝達関数を求める方法について説明を行った。実施の形態１に示した方法を用いることにより、一度の測定で、各スピーカ９から各マイク１０に対応する伝達関数を求めることができるので、マイク位置の変動（聴取者の耳位置のぐらつき等）があっても高い検出精度で伝達関数の測定を行うことが可能となる。
【０１２２】
一方で、スピーカ９からマイク１０までの伝達関数は、そのスピーカ９に対するマイク１０の設置位置に応じて変化する。従って、特定の空間における伝達関数を算出する場合、通常は、聴取者の聴取位置を予め決定し、決定された聴取位置で最適な音響効果を得ることができるように、決定位置にマイク１０を設置して収録作業を行うことが一般的である。
【０１２３】
一方で、再現音場において聴取者の聴取位置（つまりマイク１０の設置位置）が固定されておらず、移動（変化）する場合などにおいては、その変動位置に応じて、個別に伝達関数の測定を行う必要があった。特に、従来の方法では、測定位置において各スピーカ９より順番に評価音を出力させることによって、各スピーカ９から各マイク１０に対応する伝達関数を算出する必要が有るため、連続的に測定位置を変更させて伝達関数の測定を行うことが困難であった。
【０１２４】
しかしながら、実施の形態１に示した方法を用いることにより、一度の測定で、各スピーカから各マイクに対応する全ての伝達関数を求めることができるため、測定位置を連続的に移動させることにより、より具体的には、長時間の連続した測定データの中から、一部のデータを用いて伝達関数を求めることにより、伝達関数の時間的変化を求めることが可能となる。実施の形態２では、このような時間的な伝達関数の変化を求めるために、時間重みを適用した場合における伝達関数の算出方法について説明する。
【０１２５】
なお、実施の形態２に係る音場測定システムは、実施の形態１において図１を用いて説明した音場測定システム１と同一の構成であり、同一の機能を有するシステムを用いることができるため、実施の形態２において同一符号を用いると共に、その詳細な説明は省略する。また、時間重みを適用した具体的な伝達関数の算出は、実施の形態１に示した伝達関数算出部１４で行われるものとし、その詳細な構成についての説明は省略する。以下、伝達関数算出部１４における処理内容について説明を行う。
【０１２６】
測定データの全てを用いた関係式は、式１９に示すように
Ｙｔ（ｚ）＝Ｘｔ（ｚ）Ｈｔ（ｚ）
で表すことができ、式２３に示すように、この式と等価な連立一次方程式を示す
Ｄ＝ＡＱ
を解くことにより、伝達関数Ｈｔ（ｚ）を求めることができる。
【０１２７】
このようにして求められる伝達関数が、長時間の連続した測定で得られた測定データの一部を用いることにより求められる場合には、長時間の測定により得られた測定データのうち、任意の時間範囲を解析する処理が必要とされる。このため、本実施の形態２に係る伝達関数算出部１４では、時間窓を表す対角行列Ｗを用いることにより、式２３に示す方程式を求めて、任意の時間範囲の伝達関数を求める。
【０１２８】
具体的には、式２３に示す連立一次方程式の両辺に対して対角行列Ｗ
【数３５】
を掛けて、
ＷＤ＝（ＷＡ）Ｑ ・・・式５１
を解くことにより、所定の時間範囲に対応する伝達関数を求めることができる。
なお、この対角行列Ｗは、対角成分が０である行は、伝達関数を求めるための計算において無視され、その行の測定データ（その時刻の測定データ）が使用されないことを意味する。
【０１２９】
ＷＤ＝（ＷＡ）Ｑ ・・・式５１
の最小二乗解は、
ＡｔＷ２Ｄ＝（ＡｔＷ２Ａ）Ｑ ・・・式５２
の解として、
Ｑ＝（ＡｔＷ２Ａ）−１（ＡｔＷ２Ｄ） ・・・式５３
と表すことができる。
この式５３に示す連立方程式の数値解を得るための計算手順を図８〜図１２に示す。
【０１３０】
図８〜図１２の計算手順は、連立一次方程式を解く大きな枠組みとして共役勾配法を採用し、個々の計算過程で必要な演算を畳み込み問題に特化して効率化する方法を用いたものとなっている。共役勾配法の実施で問題となる係数行列とベクトルの演算については、図８に示す畳み込み行列Ａと任意のベクトルＶとの積の計算と、図９に示す行列Ａの転置行列Ａｔと任意のベクトルＶとの積の計算を用いることにより、高速フーリエ変換（ＦＦＴ）を利用した畳み込み演算または相関演算で効率的に計算する方法を採用している。
まず、図８に示す、行列ＡとベクトルＶとの積の計算について説明する。
【０１３１】
行列Ａに対する積算対象となる任意のベクトルＶを
【数３６】
とする。ベクトルＶは、長さＬのサブベクトルをＳ個並べたベクトルとして定義される（ステップＳ．６１）。
【０１３２】
このようにして示されるＶのそれぞれに対して、Ｎconv≧Ｎｘ＋Ｌ−１の条件を満たすＮconvを決定する（ステップＳ．６２）。本実施の形態２では、例えば、Ｎconv＝Ｎｘ＋Ｌを用いる。
【０１３３】
次に、伝達関数算出部１４は、入力信号波形であるｘｊ（ｔ）（但しｊ＝１，２，・・・，Ｓ）に対してＮconv点の高速フーリエ変換を適用する（ステップＳ．６３）。
【０１３４】
なお、ｘｊ（ｔ）における時系列のデータ数がＮconvに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加してデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められる時間領域の演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１３５】
高速フーリエ変換を適用することにより、
Ａ（ω）≡Ｘｔ（ω）＝（Ｘ１（ω） Ｘ２（ω） …ＸＳ（ω））
・・・式５５
となる。
【０１３６】
一方で、伝達関数算出部１４は、ステップ６１において定義されたベクトルＶに対して、Ｎconv点の高速フーリエ変換を適用する（ステップＳ．６４）。
【０１３７】
なお、ベクトルＶの構成要素であるｖｊ（ｔ）（ｊ＝１，２，・・・，Ｓ）における時系列のデータ数がＮconvに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用して求められる時間領域の演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１３８】
ステップＳ．６４においてＮconv点で高速フーリエ変換（ＦＦＴ）すると、
【数３７】
を求めることができる（ステップＳ．６４）。
【０１３９】
このようにして求められたＶ（ω）とステップＳ．６３において求められたＡ（ω）との積を、周波数領域において算出する（周波数領域において畳み込み算出を実施する）と、
【数３８】
を求めることができる（ステップＳ．６５）。
【０１４０】
そして、伝達関数算出部１４は、このようにしても求められた（ＡＶ）（ω）の演算結果に対して逆フーリエ変換を適用することにより、時間領域におけるデータ（時系列データ）ＡＶに変換を行う（ステップＳ．６６）。
【０１４１】
このように、行列ＡおよびベクトルＶをそれぞれ高速フーリエ変換することにより周波数領域におけるデータの変換を行い、変換された周波数領域において畳み込み演算を行った後に、演算結果を逆フーリエ変換により時間領域に戻して演算を行うことにより、伝達関数算出部１４における演算処理の負担を軽減させることができる。
【０１４２】
次に、図９に示す、転置行列ＡｔとベクトルＶとの積の計算について説明する。
【０１４３】
転置行列Ａｔに対する積算対象となる任意のベクトルＶを
【数３９】
とする。ベクトルＶは、長さＬのサブベクトルをＳ個並べたベクトルとして定義される（ステップＳ．７１）。
【０１４４】
このようにして示されるＶのそれぞれに対して、Ｎcorr≧Ｎｘ＋Ｌ−１の条件を満たすＮcorrを決定する（ステップＳ．７２）。本実施の形態２では、例えば、Ｎcorr＝Ｎｘ＋Ｌを用いる。
【０１４５】
次に、伝達関数算出部１４は、入力信号波形であるｘｊ（ｔ）（但しｊ＝１，２，・・・，Ｓ）に対してＮcorr点の高速フーリエ変換を適用する（ステップＳ．７３）。
【０１４６】
なお、ｘｊ（ｔ）における時系列のデータ数がＮcorrに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１４７】
高速フーリエ変換を適用することにより、
Ａ（ω）≡Ｘｔ（ω）＝（Ｘ１（ω） Ｘ２（ω） …ＸＳ（ω））
・・・式５９
となる。
【０１４８】
一方で、伝達関数算出部１４は、ステップ７１において定義されたベクトルＶに対して、Ｎcorr点の高速フーリエ変換を適用する（ステップＳ．７４）。
なお、ベクトルＶの構成要素であるｖｊ（ｔ）（ｊ＝１，２，・・・，Ｓ）における時系列のデータ数がＮcorrに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１４９】
ステップＳ．７４においてＮcorr点で高速フーリエ変換（ＦＦＴ）すると、
【数４０】
を求めることができる。
【０１５０】
このようにして求められたＶ（ω）と転置行列Ａｔとの積を、周波数領域において算出する（周波数領域において畳み込み算出を実施する）と、
【数４１】
を求めることができる（ステップＳ．７５）。
【０１５１】
ここで、Ａ＊（ω）は、Ａ（ω）の共役転置行列を意味しており、
【数４２】
で表される。
【０１５２】
そして、伝達関数算出部１４は、上述した式により算出された（ＡｔＶ）（ω）の演算結果に対して逆フーリエ変換を適用することにより、時間領域におけるデータ（時系列データ）（ＡｔＶ）’に変換を行う（ステップＳ．７６）。
【０１５３】
伝達関数算出部１４は、この（ＡｔＶ）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＡｔＶを求める（ステップＳ．７６）。このように、０≦ｔ≦Ｌ−１の範囲の時間窓を利用して、対応するデータのみを用いて計算を行うことによって、逆フーリエ変換によって求められた解が正確な値となるように（時間領域のみの演算により求められる解とに誤差が生じないように）することが可能となる。
【０１５４】
上述した図８および図９に示す処理手順に基づいて演算処理を行うことによって、伝達関数算出部１４における演算処理負担の低減を図ることが可能となる。
【０１５５】
次に、図１０〜図１２を用いて、長時間の連続した測定で得られた測定データの一部を用いることにより伝達関数を求める処理内容について説明する。
【０１５６】
上述したように、長時間の連続した測定で得られた測定データの一部を用いることにより伝達関数を求めるためには、
ＡｔＷ２Ｄ＝（ＡｔＷ２Ａ）Ｑ ・・・式５２
の解である、
Ｑ＝（ＡｔＷ２Ａ）−１（ＡｔＷ２Ｄ） ・・・式５３
を求めることが必要となる。
【０１５７】
伝達関数算出部１４は、式５２に示す方程式を最小二乗法で解くために、式５２の左辺をＰ、右辺の係数行列をＲとし、Ｐ＝ＲＱに基づいて
Ｐ＝Ａｔ（Ｗ２Ｄ）
Ｒ＝Ａｔ（Ｗ２Ａ）
とおいて、ＰおよびＲの計算を行い、この計算式に基づいて連立一次方程式Ｐ＝ＲＱを解くことにより、時間領域における行列Ｑ：Ｑ＝Ｒ−１Ｐを求める。
【０１５８】
なお、
Ｑ＝（ＡｔＷ２Ａ）−１（ＡｔＷ２Ｄ） ・・・式５３
の連立一次方程式を求める場合、Ｗが単位行列でない場合には、式５３の係数行列であるＡｔＷ２Ａの規則性が崩れてしまう。このため、伝達関数算出部１４では、式５３の演算を行う際に、ＡｔＷ２Ａを相関関数のベクトルとして記憶させて演算を行うことができない。
【０１５９】
簡単な数値例を挙げて説明すると、
ｘ１（ｔ）＝１，２，３，４，０・・・
ｘ２（ｔ）＝１１，１２，１３，１４，０・・・
ｗ（ｔ）＝０，０，０，１，１，１，０・・・
とし、
【数４３】
【数４４】
とすると、時間重みのない場合（つまりＷを考慮する必要がない場合）には、ＡｔＡが、
【数４５】
として求められる。このＡｔＡは、行列全体が対称行列であることに加え、（）内の個々のサブマトリックスの要素が斜め方向に同一値となる。このような行列の場合には、行列の全ての要素を個別に記憶しておく必要がないため、行列の規模が大きくなった場合であっても規則性を利用して記憶容量（メモリ容量）を節約することが可能である。
【０１６０】
しかしながら、時間重みのある場合には、
【数４６】
となり、サブマトリックスの要素が斜め方向に同一値を持つという規則性が崩れてしまう。このように行列の規模が大きくなった場合には、係数行列を記憶するために膨大な記憶容量（メモリ容量）が必要となるため、係数行列ＡｔＷ２Ａ（＝（ＷＡ）ｔＷＡ）を予め計算して変数として記憶することが困難となる。
【０１６１】
このため、本実施の形態２に係る伝達関数算出部１４では、伝達関数の計算処理において、係数行列ＡｔＷ２Ａそのものを求めることではなく、（ＡｔＷ２Ａ）ｖの形式の行列のベクトル積の演算結果を求める方法を用いている。（ＡｔＷ２Ａ）ｖをＡｔ（Ｗ２（Ａｖ））とする計算は、後述する図１２の計算手順と、図８および図９に示した一連の計算手順を用いることにより、畳み込み演算、時間重み、相関演算を順番に実施することで効率的に行うことが可能になる。従って、実施の形態２に示す数値計算方法を用いることによって、係数行列の記憶に関する問題を回避することができる。以下、上述した計算方法を用いて演算を行う手順について説明する。
【０１６２】
まず、伝達関数算出部１４は、図１０に示すように、Ｐ＝Ａｔ（Ｗ２Ｄ）の計算を行う（ステップＳ．８１）。
【０１６３】
具体的に、伝達関数算出部１４は、図１１に示すように、初期値としてｉに１を代入する（ステップＳ．９１）、ここでｉは、音響空間２におけるマイク１０の数を示しており、１に設定されたｉがマイク１０の設置数Ｍに該当するまでの間、伝達関数算出部１４はステップＳ．９２〜ステップＳ．９５の処理を繰り返し実行する。
【０１６４】
次に伝達関数算出部１４は、行列Ｄのｉ列目の要素ｄｉをとして抽出する（ステップＳ．９２）。そして、伝達関数算出部１４は、対角行列Ｗとｄｉとを用いて、Ｗ２ｄｉの計算を行う（ステップＳ．９３）。具体的に、伝達関数算出部１４は、ｗ（ｔ）ｙｉ（ｔ）の計算を行う。
【０１６５】
その後、伝達関数算出部１４は、行列Ｄのｉ列目の要素ｄｉを用いて求められたＷ２ｄｉの計算結果をベクトルＶとして、既に図９を用いて説明したように、転置行列ＡｔとベクトルＶとの演算方法を用いることにより、Ａｔ（Ｗ２ｄｉ）の計算を行う（ステップＳ．９４）。
【０１６６】
そして、伝達関数算出部１４は、全てのｉに対してＡｔ（Ｗ２ｄｉ）の計算を行ったか否か（つまり、ｉ＝Ｍであるか否か）の判断を行う（ステップＳ．９５）。全てのｉに対して計算を行っていないと判断した場合（ステップＳ．９５においてＮｏの場合）、伝達関数算出部１４は、ｉに１を加えた後に、上述したステップＳ．９２〜ステップＳ．９５の処理を繰り返し実行する（ステップＳ．９６）。
【０１６７】
一方で、全てのｉに対して計算を行ったと判断した場合（ステップＳ．９５においてＹｅｓの場合）、伝達関数算出部１４は、Ｄの全ての列に対応するｉのＡｔ（Ｗ２ｄｉ）を算出したものと判断し、つまり、Ａｔ（Ｗ２Ｄ）（＝Ｐ）を算出したものと判断して、処理を図１０のステップＳ．８１に戻す。
【０１６８】
次に、伝達関数算出部１４は、初期値としてｊに１を代入し、μに０を代入する（ステップＳ．８２）、ここでｊは、音響空間２におけるスピーカ９の数を示している。またμは、行列Ｑのｊ列目の値を算出するために実行された演算回数（反復回数）を示すパラメータであり、具体的には、ステップＳ．８３〜ステップＳ．８８の処理を繰り返し実行した回数を示している。
【０１６９】
次に伝達関数算出部１４は、初期設定を行う（ステップＳ．８３）。初期設定では、Ｐ，Ｑのｊ列目の値をｐ，ｑとおき、適当な初期値ｑ0を用いて、
ε0＝ｐ−Ａｔ（Ｗ２（Ａｑ0）） ・・・式６７
を求める。なお、εμは、反復回数μでの暫定的な解をｑμとしたときの、誤差ベクトルを示したものである。
【０１７０】
ε0＝ｐ−Ａｔ（Ｗ２（Ａｑ0）） ・・・式６７
におけるＡｔ（Ｗ２（Ａｑ0））の計算は、後述するステップＳ．８５に示す
Ｒσμ＝Ａｔ（Ｗ２（Ａσμ））
のσμに対して初期値としてｑ0を設定することによって求めることができる。
Ｒσμ＝Ａｔ（Ｗ２（Ａσμ））
の計算方法については、後述する。
なお、ｐ、ｑ、σμは、長さＬのベクトルをＳ個集めた長さＳＬのベクトルとして定義される。
【０１７１】
また、伝達関数算出部１４は、求められたε0を用いて
【数４７】
の初期値であるσ0の値
【数４８】
の値を算出する。
【０１７２】
さらに、伝達関数算出部１４は、共役勾配法の収束判定に用いられるパラメータをηとして設定する。ηは収束を判断するための閾値パラメータとしての定数であって、このパラメータに基づいて誤差が収束判定条件を満たす値となった場合（次述するステップＳ．８４においてＹｅｓの場合）、伝達関数算出部１４は、解が収束したと判断して反復計算を終了する。
【０１７３】
そして、伝達関数算出部１４は、収束条件の判断処理（ステップＳ．８４）を行う。収束条件の判断式は、
【数４９】
で表される。伝達関数算出部１４は、この判断式に基づいて解の収束を判断する。
【０１７４】
収束条件を満たさなかった場合（ステップＳ．８４においてＮｏの場合）、伝達関数算出部１４は、Ｒσμの計算を行う（ステップＳ．８５）。
【０１７５】
図１２は、Ｒσμの計算方法、つまり、Ａｔ（Ｗ２（Ａσμ））を計算する方法を示したフローチャートである。
【０１７６】
伝達関数算出部１４は、まず、図８に示した行列ＡとベクトルＶとの積の計算方法において、ベクトルＶにベクトルσμを用いることにより、Ａσμを計算する（ステップＳ．１０１）。次に、伝達関数算出部１４は、対角行列Ｗの積であるＷ２にＡσμを積算したＷ２（Ａσμ）の計算を行う（ステップＳ．１０２）。
【０１７７】
そして、伝達関数算出部１４は、図９に示した転置行列ＡｔとベクトルＶとの積の計算方法において、ベクトルＶにベクトルＷ２（Ａσμ）を用いることにより、Ａｔ（Ｗ２（Ａσμ））を計算する（ステップＳ．１０３）。
【０１７８】
このように、Ａｔ（Ｗ２（Ａσμ））の計算を行う場合にも、伝達関数算出部１４において、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部１４における処理負担を軽減させることができる。
【０１７９】
そして、伝達関数算出部１４は、求められたＲσμの値と、εμ、σμとを用いて、
【数５０】
【数５１】
【数５２】
を算出し、さらに、μに１を加える（μ＝μ＋１）処理を行う（ステップＳ．８６）。
【０１８０】
そして、伝達関数算出部１４は、ステップＳ．８６において算出されたεμとＰとを用いて、収束条件の判断処理（ステップＳ．８４）を繰り返し実行する。
【０１８１】
収束条件を満たす場合（ステップＳ．８４においてＹｅｓの場合）、伝達関数算出部１４は、ステップＳ．８６において算出されたｑμの値を、行列Ｑにおけるｊ列目の解（伝達関数）とする（ステップＳ．８７）。
【０１８２】
その後、伝達関数算出部１４は、全ての列ｊについて、つまり全てのスピーカに対応する各マイクの伝達関数（伝達関数マトリックス（行列Ｑ）における全ての伝達関数）を求めたか否かを判断し（ステップＳ．８８）、全ての列において伝達関数を求めていない場合（ステップＳ．８８においてＮｏの場合）には、ｊに１を加え、μに０を代入した後に（ステップＳ．８９）、処理をステップＳ．８３に移行して上述した処理を繰り返し実行する。
【０１８３】
全ての列において伝達関数を求めた場合（ステップＳ．８８においてＹｅｓの場合）、伝達関数算出部１４は、全ての列の伝達関数の値を求めたものと判断して、伝達関数の計算処理（伝達関数マトリックス（行列Ｑ）の算出処理）を終了する。
【０１８４】
このように、実施の形態２に係る伝達関数算出部１４では、長時間連続して同時に複数のスピーカ９から評価音を出力すると共に、スピーカ９から出力された評価音を連続してマイク１０で同時に測定した場合であっても、時間重みを適用して伝達関数を算出することができる。このため、連続的に測定位置が変更される場合であっても時間変化に応じた適切な伝達関数を求めることが可能となる。
【０１８５】
なお、実施の形態２では、上述したように、長時間連続して同時に複数のスピーカ９から評価音を同時に出力すると共に、スピーカ９から出力された評価音を連続してマイク１０で同時に測定した場合について説明を行ったが、実施の形態２に係る発明は、複数のスピーカ９から同時に評価音を出力する場合に限定されるものではない。例えば、従来の方法を用いて伝達関数を求める場合のように、各スピーカから順番に評価音を出力することにより該当するスピーカとマイクとに対応する伝達関数を算出する場合であっても、各評価音の出力を長時間連続して行うことにより、該当するスピーカと各マイクとに対応する伝達関数の算出を、時間重みを適用して求めることが可能となる。
【０１８６】
[実施の形態３]
実施の形態１に示す方法を用いることにより、複数のスピーカから同時に評価音を出力すると共に、複数のマイクから同時に評価音の収録を行うことによって、一度の測定でそれぞれのマイクとスピーカとに対応させた伝達関数を求めることが可能となる。また、実施の形態２に示す方法を用いることにより、時間重みを適用して伝達関数を算出することができるので、連続的に測定位置が変更される場合であっても時間変化に応じた適切な伝達関数を求めることが可能となる。
【０１８７】
しかしながら、音響空間２にノイズが存在する場合、マイク１０で測定された収録音（出力信号、測定データ）Ｙには、評価音（入力信号）Ｘの成分（寄与成分）に加えて、ノイズ成分が含まれることになる。
【０１８８】
例えば、図１３に示すように、測定データＹ（ｚ）を、入力信号Ｘ（ｚ）と、伝達関数Ｈ（ｚ）と、ノイズ成分Ｎ（ｚ）とを用いて表すと、
【数５３】
と表される。
【０１８９】
このように、測定データＹ（ｚ）に、入力信号Ｘ（ｚ）の寄与成分と、ノイズ成分Ｎ（ｚ）とが含まれている場合において、ノイズの影響を低減させるためには、ノイズ信号のパワーに対して入力信号のパワーを相対的に大きくすることにより、測定Ｓ／Ｎ（Signal to Noise ratio）を大きくすることが必要となる。
【０１９０】
ところで、マイク１０を使用して、音響空間２における音響特性を測定する場合には、低域周波数であるほど、レベルの高い環境ノイズが含まれる傾向があり、高い周波数に比例して低域周波数での測定Ｓ／Ｎが劣化する傾向がある。このため、スピーカ９より出力される評価音（入力信号）の低域周波数成分のレベルを大きくすることにより、収録音（測定データ）に含まれる低い周波数成分のＳ／Ｎを改善することが考えられる。
【０１９１】
従って、本実施の形態３に係る伝達関数算出部１４では、入力信号として低域周波数成分のレベルの高い信号を用いることによって、測定データに含まれる低い周波数成分のＳ／Ｎを改善する方法を採用する。
【０１９２】
図１４（ａ）（ｂ）は、入力信号の周波数特性をそれぞれ変えた場合における入力信号の信号パワー値とマイクで集音されるノイズの信号パワーとの対応を示した図である。（ａ）は、入力信号の周波数特性が平坦な場合（つまり、低域周波数成分の信号パワーの補強を行わない場合）を示しており、（ｂ）は、入力信号の周波数特性を変えることにより、低域周波数成分の信号パワーが大きくなるようにした場合を示している。また、（ａ）（ｂ）には、周波数帯域に応じて信号パワーが変化するノイズの状態を対比して示している。
【０１９３】
図１４（ａ）に示す入力信号とノイズとの信号のパワー比をとると、低域周波数帯域において、ノイズのパワー値と入力信号のパワー値とが近い値となるため、Ｓ／Ｎが小さな値になってしまう。一方で、図１４（ｂ）に示す入力信号とノイズとの信号のパワー比をとると、入力信号の信号パワーが周波数の減少にともない増加する状態を示しているため、低域周波数帯域であっても、ノイズのパワー値と入力信号のパワー値との間隔が一定に保たれることになり、Ｓ／Ｎの減少を抑制することが可能となる。
【０１９４】
しかしながら、低域周波数成分の信号パワーが大きな入力信号を用いると、自己相関特性が劣化して、連立一次方程式の係数行列（ＡｔＡ）の性質が悪くなるため、連立方程式を反復解法で解くときの収束性が劣化してしまう。このように、解の収束性が悪くなると、反復解法に必要な収束回数が多く必要になり、有限の反復回数で得られる解の精度が劣化してしまう。そこで、実施の形態３に係る伝達関数算出部１４では、収束性を改善するために、入力信号（評価音）に対して前処理を行うことにより、収束性の劣化を回避して解の精度を保つ方法が用いられている。
【０１９５】
下記の式７５は、低域成分が高いパワーとなる試験信号を用いる場合において、収束性の改善を実現することが可能な数式を示している。実施の形態３の入力信号Ｘには、予めマイク１０に加えられるノイズの周波数特性を考慮して周波数特性が平坦ではない信号を使用しているものとする。このため、
Ｙｔ（ｚ）＝Ｘｔ（ｚ）Ｈｔ（ｚ）
を解く代わりに、この式の両辺にｇ（ｚ）を畳み込む前処理を実施して、
（ｇ（ｚ）Ｙｔ（ｚ））＝（ｇ（ｚ）Ｘｔ（ｚ））Ｈｔ（ｚ） ・・・式７５
を解く。
【０１９６】
ここで、
ｇ（ｚ）Ｘｔ（ｚ）
＝（ｇ（ｚ）Ｘ１（ｚ） ｇ（ｚ）Ｘ２（ｚ） … ｇ（ｚ）ＸＳ（ｚ））
・・・ 式７６
ｇ（ｚ）Ｙｔ（ｚ）
＝（ｇ（ｚ）Ｙ１（ｚ） ｇ（ｚ）Ｙ２（ｚ） … ｇ（ｚ）ＹＭ（ｚ））
・・・ 式７７
である。
【０１９７】
式７５を解くことは、これを代数方程式で表した
ＧＤ＝（ＧＡ）Ｑ ・・・ 式７８
を解いて、
Ｑ＝（（ＧＡ）ｔＧＡ）−１（（ＧＡ）ｔＤ）・・・ 式７９
を得ることと同じである。
【０１９８】
ここで、Ｇはｇ（ｚ）の畳み込みマトリックスを示しており、
【数５４】
で表される。
【０１９９】
式７９に示す式の係数行列（ＧＡ）ｔ（ＧＡ）は、
（ＧＡ）ｔ（ＧＡ）＝Ａｔ（ＧｔＧ）Ａ
が単位行列に近づくようなＧを与えることで、解の収束性を改善することができる。例えば、前処理に用いるフィルタとして、入力信号Ｘの周波数特性とほぼ逆の周波数特性を持つ最小位相フィルタを利用することができる。
【０２００】
一般的に、周波数特性がフラットな信号の相関特性は、周波数特性がフラットでない信号より鋭い相関ピークを持ち、相関特性がよい。また最小位相フィルタを用いることにより前処理フィルタの段数を必要最小限にし、前処理に必要な計算コストの低減を図ることが可能となる。
【０２０１】
以上、本発明に係る音場測定システムについて、図面を用いて詳細に説明を行ったが、本発明に係る音場測定システムは、上述したものに限定されるわけではない。当業者であれば、特許請求の範囲に記載された範疇内において、各種の変更例または修正例に想到しうることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
【０２０２】
なお、上述した音場測定システムを用いて、多数のスピーカから一度に評価音を出力し、多数のマイクにおいて同時に評価音の収録を行うことによって、各スピーカから各マイクに対応する伝達関数を一度に算出する方法が、従来のように各スピーカより順番に評価音を出力し、多数のマイクで評価音の収録を行うことにより伝達関数を算出する場合に比べて、測定時間の短縮化を図ることが可能であることを、具体的に説明する。
【０２０３】
スピーカ（音源）の数をＳとし、伝達関数を測定するための応答時間がＬで打ち切られると仮定する。従来のように、各スピーカより順番に評価音を出力し、多数のマイクにおいて同時に評価音の収録を行うことにより伝達関数を算出する場合において、１つのスピーカ（音源）より出力される評価音をマイクで収録して測定を行う時間は、概ね
α×Ｌ
程度と見積もることができる。
【０２０４】
ここでαは測定条件（積分処理に要する時間や応答時間等）を考慮した定数であるが、測定データの積分時間や平均回数を大きくしてノイズ環境下での測定Ｓ／Ｎを確保するために、αは１より大きな値とする。全てのスピーカから出力される評価音を測定するためには、
α×Ｌ×Ｓ
程度の時間が必要とされる。
【０２０５】
一方で、本実施の形態１〜実施の形態３に示した音場測定システムを用いる場合における測定時間を求めるためには、
１）ノイズ環境下において測定Ｓ／Ｎを確保するために必要な測定時間
２）方程式の係数行列が逆行列を持つために原理上必要となる入力信号の長さ
を考慮する必要がある。
【０２０６】
まず１つ目の「ノイズ環境下での測定Ｓ／Ｎを確保するために必要な測定時間」について考える。本実施の形態に示した測定方式は、他のスピーカ（音源）から再生される音が測定ノイズとして加算されることがない方式である。このため、測定Ｓ／Ｎを確保するために必要な測定時間は、本実施の形態に示したように複数のスピーカから同時に評価音を出力する方式であっても、従来の方式により１つのスピーカで出力される音を測定する方式であっても、一度の測定で必要とされる測定時間は同等である。従って、測定時間は音源の数に依存せず、
α×Ｌ
と見積もることができる。
【０２０７】
次に２つ目の「方程式の係数行列が逆行列を持つために原理上必要となる入力信号の長さ」は
β×Ｌ×Ｓ
と見積もることができる。
なお、βの値は、入力信号（評価音）の長さとして演算処理上必要な値を示しており、応答時間Ｌに比例した値となる。このβの値は、測定環境のノイズの大きさには無関係であり、β＝２程度で十分であると判断できる。
【０２０８】
以上２つの条件を考慮すると、本実施の形態に示した測定方法を用いる際に必要とされる測定時間は、
ＭＡＸ（α×Ｌ，β×Ｌ×Ｓ）
（：（α×Ｌ）と（β×Ｌ×Ｓ）とを比較して大きな値を有効とする）
と見積もることができる。
【０２０９】
従って、従来方式を用いて測定する時間と本提案方式を用いて測定する時間との比は以下のようになる。
【数５５】
ノイズ環境下においては、α＞βとなるため、本提案方式の測定時間の方が従来方式の測定時間より短くなる。
【符号の説明】
【０２１０】
１ …音場測定システム
２ …音響空間
３ …スピーカ部
４ …音声出力装置
５ …コンピュータ
７ …マイクアレイ
８ …音声入力装置
９ …スピーカ
１０ …マイク
１２ …評価音生成部（評価音生成手段）
１３ …収録音記録部
１４ …伝達関数算出部（伝達関数算出手段）
１５ …伝達関数記録部
１６ …制御部
【技術分野】
【０００１】
本発明は、音場測定システムに関し、より詳細には、測定音場において複数のスピーカ（音源）から一斉に評価音を出力すると共に、複数のマイクで複数のスピーカより出力された評価音を収録することにより、各スピーカから各マイクに対応するそれぞれの伝達関数を一度の測定で算出することが可能な音場測定システムに関する。
【背景技術】
【０００２】
従来より、ある空間（原音場）の音場特性を収録（測定）した後、その空間とは異なる空間（再生音場）において、測定した原音場の音場環境を仮想的に再現する音場再現方法が提案されている（例えば、特許文献１参照）。
【０００３】
上述した音場再現方法では、まず、ある空間における各スピーカから各マイクへの全ての伝達関数を原音場の音響特性として求めると共に、異なる空間（再生音場）における各スピーカから各マイクへの全ての伝達関数を再生音場の音響特性として求め、この２つの音響特性に基づいて音場再現を行うための音場再現フィルタを算出する。次に、算出された音場再現フィルタを用いて、再生音場において出力される音響信号に対して音響処理（例えば、ＦＩＲ（Finite Impulse Response）フィルタを用いた畳み込み演算処理）を行うことによって、再生音場へ出力される出力音に音響処理を施す。このようにして音響処理が施された音響信号を再生音場で出力することによって、聴取者はあたかも原音場で聴取しているかのような感覚を再生音場で得ることが可能となる。
【０００４】
具体的に、既知の入力信号（入力波形）ｘ（ｔ）を入力（スピーカから出力）させ、この入力に対応する出力信号（出力波形）ｙ（ｔ）を測定する（マイクで収録する）ことにより原音場あるいは再生音場における伝達関数ｈ（ｔ）を求める場合について説明する。
【０００５】
上述した入力信号ｘ（ｔ）、伝達関数ｈ（ｔ）、出力信号ｙ（ｔ）のそれぞれをＺ変換（離散関数のフーリエ変換）することにより、Ｘ（ｚ）、Ｈ（ｚ）、Ｙ（ｚ）が求められるとすると、Ｘ（ｚ）、Ｈ（ｚ）、Ｙ（ｚ）のそれぞれは、
【数１】
【数２】
【数３】
と表される。
【０００６】
このとき、システムの伝達関数Ｈ（ｚ）は既知の入力信号Ｘ（ｚ）と出力信号Ｙ（ｚ）とを用いて
【数４】
を計算することによって求めることができる。
【０００７】
式４に示す式は、具体的に以下の２通りの計算方法を用いることにより、求めることができる。
【０００８】
（１）時間領域で計算する方法
時間領域で計算する方法とは、入力信号Ｘ（ｚ）の逆特性である１／Ｘ（ｚ）を求め、この１／Ｘ（ｚ）の特性を持つフィルタを出力信号Ｙ（ｚ）に畳み込むことにより（畳み込み演算を行うことによって）、伝達関数Ｈ（ｚ）を計算する方法を意味する。この方法を用いて伝達関数ｈ（ｔ）を求めるためには、１／Ｘ（ｚ）の特性を持つフィルタが既知であることが必要とされる。例えば、入力信号が時間伸張パルス（ＴＳＰ：Time stretched pulse）の場合は、入力信号の逆特性を持つフィルタの時間応答が既知となるため、この方法を用いることができる。
【０００９】
（２）周波数領域で計算する方法
周波数領域で計算する方法とは、式４をフーリエ変換し、周波数領域上で周波数成分毎に、
【数５】
を計算して伝達関数Ｈ（ω）の周波数特性を求める方法を意味する。時間応答に該当する伝達関数ｈ（ｔ）は、式５により求められる伝達関数Ｈ（ω）を逆フーリエ変換することにより求めることが可能である。
【００１０】
上述した（１）の方法および（２）の方法のいずれの方法を用いて伝達関数ｈ（ｔ）を求める場合であっても、Ｓ個の音源（スピーカ）より出力される音響信号をＸ（ｚ）（音源１から音源Ｓまでの各音源より出力される信号をＸ１（ｚ）〜ＸＳ（ｚ）とする）、Ｍ個のマイクにより録音される音響信号をＹ（ｚ）（マイク１からマイクＭまでの各マイクにおいて録音される信号をＹ１（ｚ）〜ＹＭ（ｚ）とする）とし、それぞれの音響信号Ｘ（ｚ）および音響信号Ｙ（ｚ）に基づいて求められる伝達関数をＨ（ｚ）（音源ｊより出力される音響信号に対するマイクｉの伝達関数をＨｉｊ（ｚ）とする）とすると、音響信号の入出力関係は、次の式で表現することができる。
【００１１】
Ｙ（ｚ）＝Ｈ（ｚ）Ｘ（ｚ） ・・・式６
【数６】
【数７】
【数８】
【００１２】
入力データとして、ある音源ｊについてのみ入力信号を加えた状態を測定すると、音響信号の入出力関係は、次の式で表現することができる。
【数９】
得られた測定データから音源ｊより出力された音響信号に対する伝達関数Ｈｉｊ（ｚ）は、
【数１０】
により求めることができる。
【００１３】
各音源（ｊ＝１からｊ＝ＳまでのＳ個の音源）に対応する伝達関数Ｈｉｊ（ｚ）を順番に測定することによって、音場全体の伝達関数ｈ（ｔ）を測定することが可能となる。
【先行技術文献】
【特許文献】
【００１４】
【特許文献１】特開２００７−４１１６４号公報
【発明の概要】
【発明が解決しようとする課題】
【００１５】
ところで、従来の伝達関数の測定方法では、１つの音源より出力された音響信号を、複数のマイクで同時に収録する作業に基づいて、スピーカ（音源）毎の伝達関数を求める。このように、スピーカ（音源）毎に順番に伝達関数を求めることによって、全ての音源に対応する伝達関数を求める方法を採用している。
【００１６】
このため、例えば、測定用のマイクを人間の耳位置に取り付けて、複数のスピーカ（音源）とマイクとの間の伝達関数を測定する場合、聴取者（人間）が姿勢を完全に静止させようとしても、姿勢変化を完全に抑制することが現実には困難であった。従って、スピーカ毎に１つずつ音響信号の出力を行っている間に、マイクの位置が変化してしまう恐れが高く、複数の音源に対して完全に同じ測定条件（この場合はマイク位置が変動しないこと）で測定データを得ることが極めて困難であるという問題があった。
【００１７】
さらに、非常に多くのスピーカを使用してスピーカ毎の伝達関数を測定する場合には、全てのスピーカの伝達関数を測定するための時間が、スピーカの設置数に比例して極めて長くなってしまうという問題があった。例えば、１つのスピーカについての測定時間が数秒の時間を要するとしても、スピーカの設置数が６０以上の場合には、総測定時間が数分単位で必要になる。
【００１８】
また、多数のスピーカ（音源）より同時に多数の音響信号を出力させ、出力された音響信号を、複数のマイクで同時に録音することにより、全てのスピーカ（音源）およびマイクに対応する伝達関数を同時に求める方法も考えられるが、従来の伝達関数の測定方法では、複数のスピーカ（音源）で出力され、複数のマイクにおいて収録された音響信号に基づいて、スピーカ（音源）およびマイクに対応する伝達関数成分を、スピーカ（音源）毎に分解する方法が確立されていなかった。このため、多数のスピーカ（音源）より同時に多数の音響信号を出力させ、出力された音響信号を、複数のマイクで同時に録音することにより、全てのスピーカ（音源）およびマイクに対応する伝達関数を同時に高精度に測定する方法を用いることができないという問題があった。
【００１９】
本発明は、上記問題に鑑みて成されたものであり、多数のスピーカ（音源）より同時に音響信号を出力させ、出力された音響信号を、複数のマイクで同時に収録することにより、全てのスピーカ（音源）およびマイクに対応する伝達関数を同時かつ高精度に測定することが可能な音場測定システムを提供することを課題とする。
【課題を解決するための手段】
【００２０】
上記課題を解決するために、本発明に係る音場測定システムは、特定の音響空間に設定される複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録することによって、各スピーカから各マイクに対応する伝達関数を算出するための音場測定システムであって、前記スピーカより出力するための評価音の生成を行う評価音生成手段と、前記マイクを用いて収録された収録音と前記評価音生成手段により生成された評価音とに基づいて、各スピーカから各マイクに対応する伝達関数を算出する伝達関数算出手段とを備え、前記伝達関数算出手段は、前記収録音と前記評価音とに基づいて、時間領域における反復解法を適用することにより、前記伝達関数の算出処理を実行し、前記反復解法の適用において用いる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させる演算方法を用い、前記評価音生成手段は、前記評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第１行列と、該第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成することを特徴とする。
【００２１】
このように、本発明に係る音場測定システムでは、反復解法を用いて伝達関数を算出する場合において、伝達関数算出手段が、相関演算処理と畳み込み演算処理とを行う場合にのみ、周波数領域における演算結果を利用して時間領域における解を求めるので、演算処理の負担軽減と演算処理に利用されるメモリ容量の低減とを図ることが可能となる。
【００２２】
さらに、評価音生成手段が、評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する行列要素を備えた第１行列と、第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成するため、伝達関数算出手段が、複数のスピーカより出力された評価音からなる複数要素の入力行列と、複数のマイクより収録された収録音からなる複数要素の出力行列とに基づいて各伝達関数を求める場合において、相関演算処理と畳み込み演算処理とを行う際に対応する伝達関数の値を一意的に求めることが可能となる。
【００２３】
このため、複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクで評価音を同時に収録する場合であっても、伝達関数算出手段において、各スピーカから各マイクに対応する伝達関数をそれぞれ算出することが可能となる。
【００２４】
また、上述した音場測定システムにおいて、前記伝達関数算出手段は、前記伝達関数の算出処理において、前記伝達関数を求めるための計算式に対して時間窓を構成する対角行列を掛け合わせ、該対角行列の要素を変更することにより、前記複数のスピーカから連続して評価音を出力すると共に前記複数のマイクから連続して収録音の収録を行った場合における任意の時間範囲の伝達関数を算出するものであってもよい。
【００２５】
このように、本発明に係る音場測定システムでは、長時間連続して同時に複数のスピーカから評価音を出力すると共に、スピーカから出力された評価音を連続してマイクで同時に測定した場合であっても、対角行列要素による時間重みを適用して任意の時間範囲における伝達関数を算出することができる。このため、連続的に収録音の収録位置が音響空間において変更される場合であっても、時間変化に応じた適切な伝達関数を求めることが可能となる。
【００２６】
さらに、上述した音場測定システムにおいて、前記評価音生成手段は、前記評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、前記マイクにより収録される収録音の低域周波数帯域におけるＳ／Ｎの改善を図るものであってもよい。
【００２７】
一般的に、特定の音響空間においてスピーカより出力される音をマイクを用いて収録する場合には、マイクで収録された収録音に対してノイズ成分が含まれる場合がある。特に、マイクを使用して音響特性を測定する場合には、低域周波数であるほど、レベルの高い環境ノイズが含まれる傾向がある。
【００２８】
このため、本発明に係る音場測定システムでは、評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、低域周波数帯域に含まれるノイズの出力値を評価音の出力値に対して相対的に低くすることができ、収録音の低域周波数帯域におけるＳ／Ｎの改善を図ることが可能となる。
【００２９】
また、上述した音場測定システムにおいて、前記伝達関数算出手段は、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させて求められた前記相関演算処理または前記畳み込み演算処理の演算結果のうち、前記時間領域のみで前記相関演算処理または前記畳み込み演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、前記時間領域における前記伝達関数の算出処理を行うものであってもよい。
【００３０】
このように、上述した音場測定システムにおいて、伝達関数算出手段が、周波数領域から時間領域に逆変換させた相関演算処理または前記畳み込み演算処理の演算結果のうち、時間領域のみで演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、伝達関数の算出処理を実行することによって、最終的に算出される伝達関数の値に誤差等が含まれてしまうことを効果的に防止することが可能となる。
【００３１】
また、上述した音場測定システムにおいて、前記誤差が生じない範囲は、前記評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定されるものであってもよい。
【００３２】
このように、誤差が生じない範囲が、評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定される場合には、算出された伝達関数の性能（品質）を十分に確保しつつ、最終的に算出される伝達関数の値に誤差等が含まれてしまうことを抑制することが可能となる。
【発明の効果】
【００３３】
本発明に係る音場測定システムによれば、評価音生成手段が、評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第１行列と、該第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成するため、伝達関数算出手段が、複数のスピーカより出力された評価音からなる複数要素の入力行列と、複数のマイクより収録された収録音からなる複数要素の出力行列とに基づいて各伝達関数を求める場合において、相関演算処理と畳み込み演算処理とを行う際に対応する伝達関数の値を一意的に求めることが可能となる。
【００３４】
このため、複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録する場合であっても、伝達関数算出手段において、各スピーカから各マイクに対応する伝達関数をそれぞれ算出することが可能となる。
【図面の簡単な説明】
【００３５】
【図１】実施の形態１に係る音場測定システムの概略構成を示すブロック図である。
【図２】実施の形態１に係る音場再現システムにおける、評価音の信号（入力信号）Ｘ（ｚ）と、収録音の信号（測定信号）Ｙ（ｚ）と、伝達関数Ｈ（ｚ）との関係を示した図である。
【図３】実施の形態１に係る伝達関数算出部において行われる、Ｄ＝ＡＱを解くための数学計算方式の処理内容を示したフローチャートである。
【図４】実施の形態１に係る伝達関数算出部において行われる、方程式を組み立てる計算処理内容を示したフローチャートである。
【図５】実施の形態１に係る伝達関数算出部において行われる、方程式を解く処理内容を示したフローチャートの前半図である。
【図６】実施の形態１に係る伝達関数算出部において行われる、方程式を解く処理内容を示したフローチャートの後半図である。
【図７】実施の形態１に係る伝達関数算出部において行われる、行列ＲとベクトルＶとの演算処理内容を示したフローチャートである。
【図８】実施の形態２に係る伝達関数算出部において行われる、行列ＡとベクトルＶとの積の計算処理内容を示したフローチャートである。
【図９】実施の形態２に係る伝達関数算出部において行われる、転置行列ＡｔとベクトルＶとの積の計算処理内容を示したフローチャートである。
【図１０】実施の形態２に係る伝達関数算出部において行われる、連続測定した測定データの一部を用いる伝達関数の計算処理内容を示したフローチャートである。
【図１１】実施の形態２に係る伝達関数算出部において行われるＡｔ（Ｗ２Ｄ）の計算処理内容を示したフローチャートの後半図である。
【図１２】実施の形態２に係る伝達関数算出部において行われる、Ａｔ（Ｗ２（Ａσμ））の計算処理内容を示したフローチャートである。
【図１３】実施の形態３に係る音場測定システムの音場空間においてノイズが存在する場合における測定データＹ（ｚ）を、入力信号Ｘ（ｚ）と、伝達関数Ｈ（ｚ）と、ノイズ成分Ｎ（ｚ）との関係を示したブロック図である。
【図１４】実施の形態３に係る音場測定システムの入力信号の周波数特性と、ノイズの周波数特性とを示した図であり、（ａ）は、入力信号の周波数特性が平坦な場合示しており、（ｂ）は、入力信号の周波数特性を変えることにより、低域周波数成分の信号パワーが大きくなるようにした場合を示している。
【発明を実施するための形態】
【００３６】
以下、本発明に係る音場測定システムについて、図面を用いて詳細に説明を行う。
【００３７】
［実施の形態１］
図１は、本実施の形態１に係る音場測定システムの概略構成を示した図である。音場測定システム１は、再現を望む音響空間（原音場）あるいは、原音場の音響環境の再現を行う音響空間（再生音場）において設置されるシステムである。
【００３８】
音場測定システム１は、図１に示すように、音響空間２に設置されるスピーカ部３と、スピーカ部３に対して音声出力を行う音声出力装置４と、音声出力装置４に対して評価音を出力するコンピュータ５と、音響空間２に設置されるマイクアレイ７と、マイクアレイ７により収録された収録音の入力を行う音声入力装置８とによって概略構成されている。
【００３９】
スピーカ部３は、図１に示すように、マイクアレイ７を中心として等距離を保つように円形に配置された複数（例えばＳ個）のスピーカ９で構成されている。このスピーカ９は、音声出力装置４より出力された評価音（音響信号）を出力することが可能となっている。
【００４０】
なお、図１に示すスピーカ部３の配置は、再生音場における一般的なスピーカの配置状態として実験室を用いた場合を一例として示している。例えば、本実施の形態１に係る音場測定システム１を原音場に設置する場合には、その原音場の状況に応じてスピーカの設置位置およびスピーカの設置個数が決定される。例えば、原音場として車両の室内空間を想定する場合には、車室の四隅に対して４個のスピーカが均等に設置される場合などが考えられる。
【００４１】
音声出力装置４は、コンピュータ５より出力された評価音を、スピーカ部３のＳ個のスピーカ９に対して同時に出力する役割を有している。ここで評価音とは、マイクアレイ７および音声入力装置８を介して収録された収録音に基づいて音響空間２の伝達関数（音響特性）を算出するために用いられる出力音であり、一般的には、収録された収録音をインパルス応答として測定することに適した音（例えば、Ｍ系列信号や時間伸張パルス（ＴＳＰ：Time Stretched Pulse））が該当する。なお、本実施の形態において用いられる評価音については、後述する。
【００４２】
マイクアレイ７は、円形に配設された複数（例えばＭ個）のマイク１０の一群により構成されている。マイクアレイ７は、各スピーカ９から等距離となる音響空間２の中心位置に設置されている。音声入力装置８は、マイクアレイ７によって測定されたスピーカ部３の出力音をコンピュータ５に出力する役割を有している。
【００４３】
コンピュータ５は、上述した評価音を生成すると共に、マイクアレイ７により収録された収録音に基づいて音響空間２における伝達関数（伝達関数マトリックス、音響特性）を算出する機能を有している。
【００４４】
コンピュータ５は、評価音生成部（評価音生成手段）１２と、収録音記録部１３と、伝達関数算出部（伝達関数算出手段）１４と、伝達関数記録部１５と、制御部１６とによって概略構成されている。評価音生成部１２は、上述したインパルス応答を収録するための評価音を生成する機能を有している。評価音生成部１２は、制御部１６の制御命令に応じて評価音を生成し、制御部１６を介して音声出力装置４に評価音を出力する。
【００４５】
収録音記録部１３は、マイクアレイ７および音声入力装置８を介して収録されたマイク毎の収録音を記録する機能を有している。伝達関数算出部１４は、音声出力装置４により出力された評価音と、収録音記録部１３に記録された収録音とに基づいて、スピーカ部３のＳ個のスピーカ９とマイクアレイ７のＭ個のマイク１０との間におけるインパルス応答をそれぞれ求めることにより、各スピーカと各マイクとに対応した伝達関数を算出する。なお、各スピーカと各マイクとに対応した伝達関数の要素の集合により、音響空間２における全体の音響特性が求められる。
【００４６】
伝達関数算出部１４は、伝達関数の算出処理を行うＣＰＵ（Central Processing Unit）と、伝達関数を算出するためのプログラム等が記録されるＲＯＭ（Read Only Memory）、ワークエリア等に使用されるＲＡＭ（Random Access Memory）等の一般的な演算処理器によって構成されている。
【００４７】
伝達関数記録部１５は、伝達関数算出部１４により求められた伝達関数を記録する機能を有している。伝達関数記録部１５に記録される伝達関数は、各スピーカと各マイクとに対応付けてそれぞれ記録される。
【００４８】
制御部１６は、音響空間２における伝達関数（音響特性）を求めるために必要な処理を行う役割を有している。制御部１６は、上述したように評価音生成部１２に対して評価音を生成させると共に、生成された評価音を音声出力装置４に出力する役割を有している。また、制御部１６は、音声入力装置８より入力される収録音を収録音記録部１３に記録させると共に、伝達関数算出部１４に伝達関数を算出させて、伝達関数記録部１５に記録させる役割を有している。
【００４９】
なお、本実施の形態１に係る音場測定システム１では、評価音生成部１２において生成された評価音が、スピーカ部３の全て（Ｓ個）のスピーカ９より出力される。マイクアレイ７では、全て（Ｓ個）のスピーカより出力された出力音を全て（Ｍ個）のマイク１０で同時に収録し、収録された収録音は、収録音記録部１３においてマイク毎に記録される。
【００５０】
次に、音声出力装置４より出力される評価音と、収録音記録部１３に記録される収録音とに基づいて、伝達関数算出部１４で音響空間２の伝達関数を算出する方法について説明する。
【００５１】
まず、スピーカ部３より出力される評価音を測定時の既知の入力信号波形（評価音信号、入力信号）をｘｊ（ｔ）（ｊはスピーカを示し、ｊ＝１，２，・・・，Ｓ）とし、マイクアレイ７において収録される収録音を既知の出力信号波形（収録音信号、出力信号）をｙｉ（ｔ）（ｉはマイクを示し、ｉ＝１，２，・・・，Ｍ）とし、それぞれをｚ変換で表現したものを、Ｘ（ｚ）、Ｙ（ｚ）とする。また、入力信号波形ｘｊ（ｔ）および出力信号波形ｙｉ（ｔ）により求められる伝達関数をｈｉｊ（ｔ）とし、ｚ変換で表現されたＸ（ｚ）およびＹ（ｚ）に基づいて求められる未知の伝達関数マトリックスをＨ（ｚ）とする。
【００５２】
音場測定システム１における測定時間を有限とし、入力信号として有限の時間長さＮｘの信号を用いると、入力信号波形ｘｊ（ｔ）の値が０でない時間範囲は、０≦ｔ≦Ｎｘ−１となる。
【００５３】
また、伝達関数ｈｉｊ（ｔ）の時間応答を有限の応答時間長さＬで近似できると仮定し、ｈｉｊ（ｔ）の値が０でない時間範囲を、０≦ｔ≦Ｌ−１であるとする。このとき、出力信号波形ｙｉ（ｔ）の値が入力信号波形ｘｊ（ｔ）に依存する時間応答も有限となり、その時間は、０≦ｔ≦Ｎｘ＋Ｌ−２（入力信号波形ｘｊ（ｔ）の０でない時間範囲（Ｎｘ−１）とｈｉｊ（ｔ）の値が０でない時間範囲（Ｌ−１）とを加えた時間範囲）の範囲となる。
【００５４】
従って、音場測定システム１において測定する出力信号波形ｙｉ（ｔ）の測定サンプル数をＮｙとすると、測定サンプル数Ｎｙは、出力信号波形ｙｉ（ｔ）の応答時間に基づいてＮｘ＋Ｌ−１となるように選べばよくなる。測定サンプル数ＮｙとしてＮｙ＝Ｎｘ＋Ｌ−１を選ぶことにより、出力信号波形ｙｉ（ｔ）の応答時間に対応する値が０以上となる測定サンプルを求めることが可能となる。
【００５５】
複数のスピーカ９と複数のマイク１０とを備えた音場測定システム１における測定系の入出力関係は、以下の関係式で表すことができる。
Ｙ（ｚ）＝Ｈ（ｚ）Ｘ（ｚ） ・・・式１２
【００５６】
但し、Ｘ（ｚ）は、
【数１１】
で示され、Ｘ（ｚ）の要素Ｘｊ（ｚ）（ｊ＝１，２，・・・，Ｓ）は、
【数１２】
を表している。
【００５７】
また、Ｙ（ｚ）は、
【数１３】
で示され、Ｙ（ｚ）の要素Ｙｉ（ｚ）（ｉ＝１，２，・・・，Ｍ）は、
【数１４】
を表している。
【００５８】
また、Ｈ（ｚ）は、
【数１５】
で示され、Ｈ（ｚ）の要素Ｈij（Ｚ）（ｉ＝１，２，・・，Ｍ、ｊ＝１，２，・・，Ｓ）は、
【数１６】
を表している。
【００５９】
式１２の両辺の転置をとると、
Ｙｔ（ｚ）＝Ｘｔ（ｚ）Ｈｔ（ｚ） ・・・式１９
となり、Ｘｔ（ｚ）は、
Ｘｔ（ｚ）＝（Ｘ１（ｚ） Ｘ２（ｚ） ・・・ ＸＳ（ｚ））
・・・ 式２０
Ｙｔ（ｚ）は、
Ｙｔ（ｚ）＝（Ｙ１（ｚ） Ｙ２（ｚ） ・・・ ＹＭ（ｚ））
・・・ 式２１
Ｈｔ（ｚ）は、
【数１７】
となる。
【００６０】
ここで、ｚ変換同士の積は、畳み込み演算マトリックスとベクトルの積として表すことができ、式１９は、等価な連立一次方程式である
Ｄ＝ＡＱ ・・・式２３
で表すことができる。但し、
【数１８】
【数１９】
【数２０】
【数２１】
で表される。
従って、伝達関数Ｈ（ｚ）を求めることは、上述した連立一次方程式Ｄ＝ＡＱを解く問題に帰着することができる。
【００６１】
音場測定システム１において、図２に示すように、測定対象である伝達関数Ｈ（ｚ）が線形の伝達関数モデルである限り、Ｄ＝ＡＱを満たす解は、必ず存在する（但し、等号が成立するのは、測定データ中にノイズがない場合となる）。このため、Ｄ＝ＡＱの解を求めることによって伝達関数Ｈ（ｚ）を算出できるようにするためには、Ｄ＝ＡＱを満たす解が１つだけ存在するように方程式を構築する必要がある。
【００６２】
具体的に、Ｄ＝ＡＱを満たす解が１つだけ存在するように方程式を構築するためには、ＡｔＡの逆行列が存在するように入力信号波形ｘ（ｔ）の組を与えることが必要とされる。つまり、入力信号波形（評価音信号、入力信号）ｘｊ（ｔ）の畳み込み行列（式２６）を行列要素とし、各スピーカに対応する行列要素（式２６）を備えた行列Ａ（式２５：第１行列）と、この行列Ａ（第１行列）の転置行列を成す行列Ａｔ（第２行列）との積であるＡｔＡの逆行列が存在し得るように入力信号波形（評価信号、評価音）ｘｊ（ｔ）が設定されればよい。
【００６３】
ＡｔＡの逆行列が存在する場合、Ｑは、Ｄ＝ＡＱの最小二乗解として求められ、次の連立一次方程式
（ＡｔＤ）＝（ＡｔＡ）Ｑ ・・・式２８
の解として、次式で求めることができる。
Ｑ＝（ＡｔＡ）−１（ＡｔＤ） ・・・式２９
【００６４】
上述した計算原理を用いることにより、複数の音源を同時に再生した場合の測定データを用いて伝達関数マトリックスのすべての成分を正確に計算することが可能となる。このため、１回の測定ですべてのスピーカ・マイク間の伝達関数を求めることが可能になる。
【００６５】
具体的な数値解を得るためには、式２３および式２８を具体的に解く数値計算アルゴリズムおよび入力信号の与え方が重要になる。一般に式２３および式２８で示される方程式は、計算規模（演算量）が極めて大きく、汎用の行列解法を用いる場合には、極めて処理能力の高い計算環境を用意する必要があるため、方程式の演算を実行することが困難である。
【００６６】
具体的な数値を挙げて説明すると、式２８の係数行列ＡｔＡの計算規模は、音源数（スピーカ９数）をＳ、伝達関数の時間応答長さをＬとした場合、Ｓ×Ｌ行、Ｓ×Ｌ列の行列を計算するための計算規模が必要となる。伝達関数を求める場合、時間応答長さＬは、通常、数千〜数万のオーダーで設定されるため、仮に時間応答長さＬを２１５＝３２７６８に設定し、データ精度を単精度実数（４ｂｙｔｅ）とすると、係数行列の演算処理に用いられる記憶容量だけでも、Ｓ×Ｓ×４Ｇｂｙｔｅ程度のメモリ容量が必要となる。従って、１６音源を一度に測定したい場合には、１０００Ｇｂｙｔｅ程度のメモリ容量が係数記憶のためだけに必要になる。さらに、このような大規模の方程式を解くための演算コストや計算時間も膨大なものとなってしまう。
【００６７】
そのため、本実施の形態１に係る伝達関数算出部１４では、畳み込み問題の性質に特化して効率的な計算を行う計算アルゴリズムを用いている。図３〜図６は、本実施の形態１における計算アルゴリズムを示したフローチャートである。
【００６８】
図３は、Ｄ＝ＡＱを解くための数学計算方式の概要を示した図である。図３に示すように、まず、伝達関数算出部１４は、方程式Ｄ＝ＡＱを最小二乗法で解くために、Ｐ＝ＲＱに基づいて
Ｐ＝ＡｔＡ
Ｒ＝ＡｔＤ
とおいて、ＰおよびＲの計算を行い（ステップＳ．１）、この計算式に基づいて連立一次方程式Ｐ＝ＲＱを解くことにより、時間領域における伝達関数マトリックス（行列Ｑ）（Ｑ＝Ｒ−１Ｐ）を求める（ステップＳ．２）
【００６９】
以下、伝達関数算出部１４が、時間領域において伝達関数マトリックス（行列Ｑ）を求める処理を説明する。伝達関数算出部１４は、時間領域において
Ｑ＝（ＡｔＡ）−１（ＡｔＤ） ・・・式２９
の最適解を求めるために、反復解法を用いるものとし、特に本実施の形態１では反復解法の解法例として共役勾配法を用いて演算を行うものとする。
【００７０】
伝達関数マトリックス（行列Ｑ）を算出するためには、上述したように、
Ｑ＝（ＡｔＡ）−１（ＡｔＤ） ・・・式２９
を解く必要がある。ここで、反復解法を用いて伝達関数マトリックス（行列Ｑ）を算出する場合、反復計算の処理において必要とされる行列形式は、規則的な畳み込み行列になるという特徴がある。このため、本実施の形態１に係る伝達関数算出部１４では、この畳み込み行列を利用した特徴的な演算処理を行うことによって、計算過程に必要なメモリ量をフィルタ段数Ｌに比例するオーダーまで削減することが可能となる。このように、連立一次方程式を解く大きな枠組みとして、反復法（共役勾配法）を用いることによって、係数行列の規則性を保ち、メモリの使用量を抑えたままでの計算を可能にすることができる。
【００７１】
また、本実施の形態１に係る計算アルゴリズムでは、係数行列ＡｔＡが相関行列となるため、方程式の組み立て処理の段階で相関行列の規則性を利用することにより、計算に必要な記憶容量をＬ×Ｌに比例するオーダーから、Ｌに比例するオーダーへと大幅に削減することができる。
【００７２】
従って、方程式の組み立て処理と反復計算の処理とに用いられる主要な行列演算を、高速フーリエ変換（ＦＦＴ）を利用した相関演算処理と畳み込み演算処理とで演算することによって、計算量を大幅に削減することができ、音場測定システムの規模が大きくなっても、一般的なコンピュータを用いて演算処理を行うことが可能となる。
【００７３】
なお、上述する計算アルゴリズムにおいて入力信号波形として入力される評価音（入力信号）には、連立一次方程式
（ＡｔＤ）＝（ＡｔＡ）Ｑ ・・・式２８
の係数行列ＡｔＡの逆行列が存在するような信号系列を選ぶことが必要である。
【００７４】
このように、係数行列ＡｔＡの逆行列が存在するような信号系列を選ぶことによって、原理上、計算アルゴリズムにおいて解の精度が劣化することを防止できる。なお、有限桁数の数値計算で有限の計算時間で解を求める関係上、係数行列が単位行列に近くなるように入力信号を選ぶことが望ましい。係数行列が単位行列に近いほど、連立一次方程式を反復解法で解く場合の収束回数が少なくなるので、有限の反復回数で計算処理を打ち切った場合に得られる解の精度を向上させることができる。
【００７５】
このような入力信号として、正規分布や一様分布などの自己相関特性が鋭くて（自己相関波形が時刻０でのみ大きな値を持ち、その他の時刻では０に近い値を持つもの）、相互相関の低い信号を用いることが好ましい。従って、評価音生成部１２では、このような条件を満たし得る入力信号（評価音）の生成を行う。
【００７６】
図４は、図３に示したステップＳ．１の処理内容である方程式を組み立てるための計算処理を示したフローチャートである。また、図５および図６は、図３に示したステップＳ．２の処理内容である方程式を解く計算処理を示したフローチャートであり、共役勾配法を使用して伝達関数マトリックス（行列Ｑ）を算出する処理を示している。
【００７７】
フローチャートに従って時間領域だけで演算処理を行う場合には、次述する（１）〜（４）の演算処理において処理負担が増大するという問題があった。
（１）行列Ｒ：Ｒ＝ＡtＡの演算処理
（２）行列Ｐ：Ｐ＝ＡtＤの演算処理
（３）初期値計算におけるＲとベクトルとの積Ｒｑ0の演算処理
（４）反復計算におけるＲとベクトルとの積Ｒσμの演算処理
また、行列Ｒの演算に要する記憶容量も増大するという問題があった。
【００７８】
本実施の形態１に係る計算アルゴリズムでは、上述したように、（１）〜（４）の演算を行う場合に、高速フーリエ変換（ＦＦＴ）を用いて周波数領域における演算処理を行い、演算処理された演算結果を逆高速フーリエ変換（ＩＦＦＴ）によって時間領域の演算結果に逆変換する（時間領域の演算結果に戻す）ことによって、演算処理負担の軽減を図っている。
【００７９】
まず、図４に示すフローチャートに基づいて、初期演算として行列Ｒと行列Ｐとを計算する場合について説明する。
【００８０】
行列Ｒは、上述したようにＲ＝ＡtＡで示される。
ここで行列Ｒの部分行列Ｒi,jは、相関行列であって規則正しいパターンを持ち、斜め方向の要素の値が同一の値になるという特徴を有している。この部分行列Ｒi,jでは、２Ｌ−１個の異なる値を備えている。このため、係数行列Ｒの演算を行う場合には、重複する値をメモリ等の記憶領域に記憶せず、異なる値のみをベクトルとして記憶することによって、伝達関数算出部１４の演算処理で使用されるメモリ使用量をＬのオーダーで削減することが可能となる。
【００８１】
また、行列Ｒi,jの各要素ｒi,jは、伝達関数同士の相関演算の形を含んでいる。このため、相関演算処理を行う場合、高速フーリエ変換（ＦＦＴ）を用いてｒi,jの値を周波数領域の値に変換して解を求めた後に、求められた解を逆フーリエ変換（ＩＦＦＴ）により逆変換して、時間領域におけるｒi,jの相関演算処理の解を求める。このように周波数領域において演算処理を行うことによって、処理負担を軽減させることができる。
【００８２】
また、行列Ｐは、上述したようにＰ＝ＡtＤで示され、行列Ｐi,jの各要素ｐi,jも、行列Ｒと同様に、伝達関数同士の相関演算の形を含んでいるため、ｐi,jに対して高速フーリエ変換を適用し、ｐi,jの値を周波数領域の値に変換して解を求め、求められた解を逆フーリエ変換して時間領域の解に逆変換することにより値を求めることができる。このように周波数領域の解を時間領域の解に逆変換する処理を用いることにより、演算処理の負担を軽減させることができる。
【００８３】
ここで、伝達関数算出部１４は、ＲおよびＰの高速フーリエ変換を行う場合における段数の数Ｎcorrとして、Ｎcorr≧Ｎｙ＋Ｌ−１を満たす値、例えば、Ｎcorr＝Ｎｙ＋Ｌを設定する（ステップＳ．２１）。
【００８４】
そして、伝達関数算出部１４は、入力信号波形であるｘｊ（ｔ）（但しｊ＝１，２，・・・，Ｓ）および出力信号波形であるｙｉ（ｔ）（但しｉ＝１，２，・・・，Ｍ）に対してＮcorr点の高速フーリエ変換を適用する（ステップＳ．２２）。
【００８５】
なお、ｘｊ（ｔ）およびｙｉ（ｔ）における時系列のデータ数がＮcorrに満たない場合、伝達関数算出部１４は、各データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加してデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。このように調整を行うことにより、共役勾配法を用いて求められる伝達関数に誤差が生じてしまうことを防止することができる。
【００８６】
高速フーリエ変換を適用することにより、
Ａ（ω）≡Ｘｔ（ω）＝（Ｘ１（ω） Ｘ２（ω） …ＸＳ（ω））
・・・式３０
Ｄ（ω）≡Ｙｔ（ω）＝（Ｙ１（ω） Ｙ２（ω） …ＹＭ（ω））
・・・式３１
となる。
【００８７】
次に、伝達関数算出部１４は、求められたＡ（ω）およびＤ（ω）を用いて、Ｒ（ω）およびＰ（ω）の周波数領域における相関演算を実施する（ステップＳ．２３）。
式３０および式３１より、Ｒ（ω）およびＰ（ω）は、
Ｒ（ω）＝Ａ＊（ω）Ａ（ω） ・・・式３２
Ｐ（ω）＝Ａ＊（ω）Ｄ（ω） ・・・式３３
で示すことができる。ここで、Ａ＊（ω）は、Ａ（ω）の共役転置行列を示しており、
Ａ＊（ω）は、
【数２２】
で表される。
【００８８】
また、Ｒ（ω）の要素ｒi,j（ω）は、
ｒi,j（ω）＝conj（Ｘｉ（ω））Ｘｊ（ω）
で表すことができ、Ｐ（ω）の要素ｐi,j（ω）は、
ｐi,j（ω）＝conj（Ｘｉ（ω））Ｙｊ（ω）
で表すことができる。
なお、conj（Ｘ）はＸの複素共役を示し、conj（Ｙ）はＹの複素共役を示している。
【００８９】
このように、フーリエ変換されたｒi,j（ω）およびｐi,j（ω）の値を周波数成分の掛け算を用いて求めた後に、求めた解を逆フーリエ変換して時間領域の時系列データに変換することによって、時間領域において演算を行う場合に比べて迅速かつ容易に解を求めることができる。
【００９０】
伝達関数算出部１４は、求められたｒi,j（ω）の値を逆フーリエ変換（ＩＦＦＴ）することにより、時系列データであるｒi,j（ｔ）に変換する（ステップＳ．２４）。この場合、伝達関数算出部１４は、逆フーリエ変換によって求められた解が正確な値となるように（時間領域のみの演算により求められる解とに誤差が生じないように）、−Ｌ＋１≦ｔ≦Ｌ−１の範囲のデータのみを用いて解を求める。つまり、誤差の生じない範囲は、出力音の残響時間等に応じて予め設定される有限の応答時間長さＬに応じて決定される。
【００９１】
伝達関数算出部１４は、−Ｌ＋１≦ｔ≦Ｌ−１の範囲のデータのうち、ｔ＜０のときに、ｒi,j（ｔ＋Ｎcorr）を負の時刻のｒi,j（ｔ）として所定の記憶手段に記憶し、ｔ≧０のときに、ｒi,j（ｔ）をそのまま記憶する。
【００９２】
一方で、伝達関数算出部１４は、求められたｐi,j（ω）の値を逆フーリエ変換（ＩＦＦＴ）することにより、時系列データであるｐi,j（ｔ）に変換する（ステップＳ．２５）。この場合、伝達関数算出部１４は、逆フーリエ変換によって求められた解が正確な値となるように、０≦ｔ≦Ｌ−１の範囲のデータのみを用いて解を求め、記憶する。
【００９３】
次に、図５および図６に示すフローチャートに基づいて、伝達関数算出部１４の具体的な伝達関数の演算処理について説明する。
【００９４】
伝達関数算出部１４は、伝達関数の具体的な演算処理においても、行列Ｒi,jの各要素ｒi,jに含まれる伝達関数同士の相関演算の形を利用する。伝達関数算出部１４は、相関演算処理を行う場合において、高速フーリエ変換（ＦＦＴ）を用いてｒi,jの値を周波数領域の値に変換して解を求めた後に、求められたｒi,jを周波数領域において演算処理し、演算処理された演算結果を逆フーリエ変換（ＩＦＦＴ）を用いて逆変換して、時間領域におけるｒi,jの相関演算処理の解を求める。このように周波数領域において演算処理を行うことによって、処理負担を軽減させることができる。
【００９５】
まず、伝達関数算出部１４は、Ｒに対して高速フーリエ変換を行う場合の段数をＮＣＧとして、ＮＣＧ≧２Ｌ−１を満たす値、例えば、ＮＣＧ＝２Ｌを設定する（ステップＳ．３１）。
【００９６】
そして、伝達関数算出部１４は、図４のステップＳ．２４において求められたｒi,j（ｔ）を並べ直したベクトルｒ’i,j（ｔ） ０≦ｔ≦ＮＣＧ−１を定義する（ステップＳ．３２）。具体的にｒ’i,j（ｔ）は、ｔの値に応じて、下記のように示される。
【数２３】
【００９７】
伝達関数算出部１４は、ｒ’i,j（ｔ）のＮＣＧ点における高速フーリエ変換（ＦＦＴ）を行うことにより、周波数領域におけるＲＣＧ（ω）のｉｊ成分を計算する（ステップＳ．３３）。
【００９８】
次に、伝達関数算出部１４は、初期値としてｊに１を代入し、μに０を代入する（ステップＳ．３４）、ここでｊは、音響空間２におけるスピーカ９の数を示している。またμは、行列Ｑのｊ列目の値を算出するために実行された演算回数（反復回数）を示すパラメータであり、具体的には、ステップＳ．３６〜ステップＳ．３８の処理を繰り返し実行した回数を示している。
【００９９】
次に伝達関数算出部１４は、初期設定を行う（ステップＳ．３５）。初期設定では、Ｐ，Ｑのｊ列目の値をｐ，ｑとおき、適当な初期値ｑ0を用いて、
【数２４】
の初期値であるε0の値
ε0＝ｐ−Ｒｑ0 ・・・式３７
を求める。なお、εμは、反復回数μでの暫定的な解をｑμとしたときの、誤差ベクトルを示したものである。
【０１００】
ここで。ε0の算出に用いられる行列Ｒと任意のベクトルＶとの演算処理（例えば、Ｒｑ0の演算）について、図７に示すフローチャートに沿って説明する。
【０１０１】
係数行列Ｒに対する積算対象となる任意のベクトルＶを
【数２５】
とする。ベクトルＶは、長さＬのサブベクトルをＳ個並べたベクトルとして定義される（ステップＳ．５１）。
【０１０２】
このようにして示されるベクトルＶのそれぞれを、ステップＳ．３１において設定したＮＣＧ点で高速フーリエ変換（ＦＦＴ）すると、
【数２６】
を求めることができる（ステップＳ．５２）。
【０１０３】
このようにして求められたＶ（ω）とステップＳ．３３において求められたＲＣＧ（ω）との積を、周波数領域において算出すると、
【数２７】
を求めることができる（ステップＳ．５３）。このようにして算出されたＲＣＧ（ω）Ｖ（ω）の算出結果を時間領域に変換（逆フーリエ変換）することにより
【数２８】
と示すことができ（ステップＳ．５４）、この（ＲＶ）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＲＶを求めることができる（ステップＳ．５５）。なお、０≦ｔ≦Ｌ−１の範囲のデータのみを用いて計算を行うのは、逆フーリエ変換によって求められた解が正確な値となるように（時間領域のみの演算により求められる解とに誤差が生じないように）するためのである。
【０１０４】
例えば、Ｒｑ0の演算を行う場合には、積算対象となるベクトルｑ0をベクトルＶと置き換えて、図７に示すフローチャートの計算処理を行うことにより、周波数領域における演算処理を用いてＲｑ0の演算を行うことができる。
【０１０５】
具体的に伝達関数算出部１４は、ｑ0をＮＣＧ点で高速フーリエ変換（ＦＦＴ）してｑ0（ω）を計算する。伝達関数算出部１４では、ステップＳ．３３において計算されたＲＣＧ（ω）と、ｑ0（ω）とを周波数領域において積算することにより、ＲＣＧ（ω）ｑ0（ω）を算出し、算出されたＲＣＧ（ω）ｑ0（ω）の算出結果を時間領域に変換（逆フーリエ変換）することにより（Ｒｑ0）’である
（Ｒｑ0）’＝ＩＦＦＴ（ＲＣＧ（ω）ｑ0（ω）） ・・・式４２
を求める。そして、伝達関数算出部１４は、求められた（Ｒｑ0）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＲｑ0を求めることができる。このように、伝達関数算出部１４において、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部１４における処理負担を軽減させることができる。
【０１０６】
また、伝達関数算出部１４は、
【数２９】
の初期値であるσ0の値
【数３０】
の値を算出する（ステップＳ．３５）。ここで、σμは、後述する演算処理におけるＲσμの計算に用いられる。
【０１０７】
さらに、伝達関数算出部１４は、共役勾配法の収束判定に用いられるパラメータをηとして設定する（ステップＳ．３５）。ηは収束を判断するための閾値パラメータとしての定数であって、このパラメータに基づいて誤差が収束判定条件を満たす値となった場合（次式４５を満たす場合）、伝達関数算出部１４は、解が収束したと判断して反復計算を終了する。
【０１０８】
そして、伝達関数算出部１４は、収束条件の判断処理（ステップＳ．３６）を行う。収束条件の判断式は、
【数３１】
で表される。伝達関数算出部１４は、この判断式に基づいて解の収束を判断する。
【０１０９】
収束条件を満たさなかった場合（ステップＳ．３６においてＮｏの場合）、伝達関数算出部１４は、Ｒσμの計算を行う（ステップＳ．３７）
【０１１０】
Ｒσμの計算を行う場合にも、図７に示すフローチャートを用いて説明したように、行列Ｒとベクトルとの積を求める計算に該当する畳み込み演算処理を用いて計算を行うことが可能である。
【０１１１】
具体的に伝達関数算出部１４は、σμをＮＣＧ点で高速フーリエ変換（ＦＦＴ）してσμ（ω）を計算する。伝達関数算出部１４では、ステップＳ．３３において計算されたＲＣＧ（ω）と、σμ（ω）とを周波数領域において積算することにより、ＲＣＧ（ω）σμ（ω）を算出し、算出されたＲＣＧ（ω）σμ（ω）の算出結果を時間領域に変換（逆フーリエ変換）することにより、
（Ｒσμ）’＝ＩＦＦＴ（ＲＣＧ（ω）σμ（ω）） ・・・式４６
を求める。
【０１１２】
この場合においても、伝達関数算出部１４は、時間領域において解に誤差が生じないように、０≦ｔ≦Ｌ−１の範囲のデータのみを用いて解を求める。従って、伝達関数算出部１４は、求められた（Ｒσμ）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＲｑ0を求める。このように、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部１４における処理負担を軽減させることができる。
【０１１３】
そして、伝達関数算出部１４は、求められたＲσμの値と、εμ、σμとを用いて、
【数３２】
【数３３】
【数３４】
を算出し、さらに、μに１を加える（μ＝μ＋１）処理を行う（ステップＳ．３８）。
【０１１４】
そして、伝達関数算出部１４は、ステップＳ．３８において算出されたεμとＰとを用いて、収束条件の判断処理（ステップＳ．３６）を繰り返し実行する。
【０１１５】
収束条件を満たす場合（ステップＳ．３６においてＹｅｓの場合）、伝達関数算出部１４は、ステップＳ．３８において算出されたｑμの値を、行列Ｑにおけるｊ列目の解（伝達関数）とする（ステップＳ．３９）。
【０１１６】
その後、伝達関数算出部１４は、全てのｊについて、つまり全てのスピーカ９に対応する各マイク１０の伝達関数マトリックス（行列Ｑ）における全ての伝達関数を求めたか否かを判断し（ステップＳ．４０）、全ての列において伝達関数を求めていない場合（ステップＳ．４０においてＮｏの場合）には、ｊに１を加え、μに０を代入した後に（ステップＳ．４１）、処理をステップＳ．３５に移行して上述した処理を繰り返し実行する。
【０１１７】
全ての列において伝達関数を求めた場合（ステップＳ．４０においてＹｅｓの場合）、伝達関数算出部１４は、全ての列の伝達関数の値を求めたものと判断して、方程式を解く計算処理（伝達関数マトリックス（行列Ｑ）の算出処理）を終了する。
【０１１８】
このようにして伝達関数算出部１４が、複数のスピーカ９より出力される入力信号波形に関するデータｘｊ（ｔ）と、複数のマイクにより収録される出力信号波形に関するデータｙｉ（ｔ）とに基づいて伝達関数の算出を行うことによって、同時に複数のスピーカ９から評価音を出力すると共に、スピーカ９から出力された評価音をマイク１０で同時に測定することにより、各スピーカ９から各マイク１０に対応する伝達関数を同時に算出することが可能となる。このため、従来のように、評価音をそれぞれのスピーカ９より順番に出力することにより、出力先となるスピーカ９と収録先となるマイク１０とを対応させて伝達関数を算出する必要がなくなり、伝達関数の測定時間の短縮化を図ることが可能となる。
【０１１９】
また、一度に全てのスピーカから評価音を出力させて、同時に全てのマイクで収録を行うことにより伝達関数の算出を行うため、従来のように集音時におけるマイク位置の変動などによる影響を受けにくくなり、測定精度の向上を図ることが可能となる。
【０１２０】
［実施の形態２］
次に、実施の形態２に係る音場測定システムについて説明を行う。
【０１２１】
実施の形態１に係る音場測定システム１では、複数のスピーカ９から同時に評価音を出力すると共に、複数のマイク１０から同時に評価音の収録を行うことにより、一度の測定でそれぞれのマイク１０とスピーカ９とに対応させた伝達関数を求める方法について説明を行った。実施の形態１に示した方法を用いることにより、一度の測定で、各スピーカ９から各マイク１０に対応する伝達関数を求めることができるので、マイク位置の変動（聴取者の耳位置のぐらつき等）があっても高い検出精度で伝達関数の測定を行うことが可能となる。
【０１２２】
一方で、スピーカ９からマイク１０までの伝達関数は、そのスピーカ９に対するマイク１０の設置位置に応じて変化する。従って、特定の空間における伝達関数を算出する場合、通常は、聴取者の聴取位置を予め決定し、決定された聴取位置で最適な音響効果を得ることができるように、決定位置にマイク１０を設置して収録作業を行うことが一般的である。
【０１２３】
一方で、再現音場において聴取者の聴取位置（つまりマイク１０の設置位置）が固定されておらず、移動（変化）する場合などにおいては、その変動位置に応じて、個別に伝達関数の測定を行う必要があった。特に、従来の方法では、測定位置において各スピーカ９より順番に評価音を出力させることによって、各スピーカ９から各マイク１０に対応する伝達関数を算出する必要が有るため、連続的に測定位置を変更させて伝達関数の測定を行うことが困難であった。
【０１２４】
しかしながら、実施の形態１に示した方法を用いることにより、一度の測定で、各スピーカから各マイクに対応する全ての伝達関数を求めることができるため、測定位置を連続的に移動させることにより、より具体的には、長時間の連続した測定データの中から、一部のデータを用いて伝達関数を求めることにより、伝達関数の時間的変化を求めることが可能となる。実施の形態２では、このような時間的な伝達関数の変化を求めるために、時間重みを適用した場合における伝達関数の算出方法について説明する。
【０１２５】
なお、実施の形態２に係る音場測定システムは、実施の形態１において図１を用いて説明した音場測定システム１と同一の構成であり、同一の機能を有するシステムを用いることができるため、実施の形態２において同一符号を用いると共に、その詳細な説明は省略する。また、時間重みを適用した具体的な伝達関数の算出は、実施の形態１に示した伝達関数算出部１４で行われるものとし、その詳細な構成についての説明は省略する。以下、伝達関数算出部１４における処理内容について説明を行う。
【０１２６】
測定データの全てを用いた関係式は、式１９に示すように
Ｙｔ（ｚ）＝Ｘｔ（ｚ）Ｈｔ（ｚ）
で表すことができ、式２３に示すように、この式と等価な連立一次方程式を示す
Ｄ＝ＡＱ
を解くことにより、伝達関数Ｈｔ（ｚ）を求めることができる。
【０１２７】
このようにして求められる伝達関数が、長時間の連続した測定で得られた測定データの一部を用いることにより求められる場合には、長時間の測定により得られた測定データのうち、任意の時間範囲を解析する処理が必要とされる。このため、本実施の形態２に係る伝達関数算出部１４では、時間窓を表す対角行列Ｗを用いることにより、式２３に示す方程式を求めて、任意の時間範囲の伝達関数を求める。
【０１２８】
具体的には、式２３に示す連立一次方程式の両辺に対して対角行列Ｗ
【数３５】
を掛けて、
ＷＤ＝（ＷＡ）Ｑ ・・・式５１
を解くことにより、所定の時間範囲に対応する伝達関数を求めることができる。
なお、この対角行列Ｗは、対角成分が０である行は、伝達関数を求めるための計算において無視され、その行の測定データ（その時刻の測定データ）が使用されないことを意味する。
【０１２９】
ＷＤ＝（ＷＡ）Ｑ ・・・式５１
の最小二乗解は、
ＡｔＷ２Ｄ＝（ＡｔＷ２Ａ）Ｑ ・・・式５２
の解として、
Ｑ＝（ＡｔＷ２Ａ）−１（ＡｔＷ２Ｄ） ・・・式５３
と表すことができる。
この式５３に示す連立方程式の数値解を得るための計算手順を図８〜図１２に示す。
【０１３０】
図８〜図１２の計算手順は、連立一次方程式を解く大きな枠組みとして共役勾配法を採用し、個々の計算過程で必要な演算を畳み込み問題に特化して効率化する方法を用いたものとなっている。共役勾配法の実施で問題となる係数行列とベクトルの演算については、図８に示す畳み込み行列Ａと任意のベクトルＶとの積の計算と、図９に示す行列Ａの転置行列Ａｔと任意のベクトルＶとの積の計算を用いることにより、高速フーリエ変換（ＦＦＴ）を利用した畳み込み演算または相関演算で効率的に計算する方法を採用している。
まず、図８に示す、行列ＡとベクトルＶとの積の計算について説明する。
【０１３１】
行列Ａに対する積算対象となる任意のベクトルＶを
【数３６】
とする。ベクトルＶは、長さＬのサブベクトルをＳ個並べたベクトルとして定義される（ステップＳ．６１）。
【０１３２】
このようにして示されるＶのそれぞれに対して、Ｎconv≧Ｎｘ＋Ｌ−１の条件を満たすＮconvを決定する（ステップＳ．６２）。本実施の形態２では、例えば、Ｎconv＝Ｎｘ＋Ｌを用いる。
【０１３３】
次に、伝達関数算出部１４は、入力信号波形であるｘｊ（ｔ）（但しｊ＝１，２，・・・，Ｓ）に対してＮconv点の高速フーリエ変換を適用する（ステップＳ．６３）。
【０１３４】
なお、ｘｊ（ｔ）における時系列のデータ数がＮconvに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加してデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められる時間領域の演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１３５】
高速フーリエ変換を適用することにより、
Ａ（ω）≡Ｘｔ（ω）＝（Ｘ１（ω） Ｘ２（ω） …ＸＳ（ω））
・・・式５５
となる。
【０１３６】
一方で、伝達関数算出部１４は、ステップ６１において定義されたベクトルＶに対して、Ｎconv点の高速フーリエ変換を適用する（ステップＳ．６４）。
【０１３７】
なお、ベクトルＶの構成要素であるｖｊ（ｔ）（ｊ＝１，２，・・・，Ｓ）における時系列のデータ数がＮconvに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用して求められる時間領域の演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１３８】
ステップＳ．６４においてＮconv点で高速フーリエ変換（ＦＦＴ）すると、
【数３７】
を求めることができる（ステップＳ．６４）。
【０１３９】
このようにして求められたＶ（ω）とステップＳ．６３において求められたＡ（ω）との積を、周波数領域において算出する（周波数領域において畳み込み算出を実施する）と、
【数３８】
を求めることができる（ステップＳ．６５）。
【０１４０】
そして、伝達関数算出部１４は、このようにしても求められた（ＡＶ）（ω）の演算結果に対して逆フーリエ変換を適用することにより、時間領域におけるデータ（時系列データ）ＡＶに変換を行う（ステップＳ．６６）。
【０１４１】
このように、行列ＡおよびベクトルＶをそれぞれ高速フーリエ変換することにより周波数領域におけるデータの変換を行い、変換された周波数領域において畳み込み演算を行った後に、演算結果を逆フーリエ変換により時間領域に戻して演算を行うことにより、伝達関数算出部１４における演算処理の負担を軽減させることができる。
【０１４２】
次に、図９に示す、転置行列ＡｔとベクトルＶとの積の計算について説明する。
【０１４３】
転置行列Ａｔに対する積算対象となる任意のベクトルＶを
【数３９】
とする。ベクトルＶは、長さＬのサブベクトルをＳ個並べたベクトルとして定義される（ステップＳ．７１）。
【０１４４】
このようにして示されるＶのそれぞれに対して、Ｎcorr≧Ｎｘ＋Ｌ−１の条件を満たすＮcorrを決定する（ステップＳ．７２）。本実施の形態２では、例えば、Ｎcorr＝Ｎｘ＋Ｌを用いる。
【０１４５】
次に、伝達関数算出部１４は、入力信号波形であるｘｊ（ｔ）（但しｊ＝１，２，・・・，Ｓ）に対してＮcorr点の高速フーリエ変換を適用する（ステップＳ．７３）。
【０１４６】
なお、ｘｊ（ｔ）における時系列のデータ数がＮcorrに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１４７】
高速フーリエ変換を適用することにより、
Ａ（ω）≡Ｘｔ（ω）＝（Ｘ１（ω） Ｘ２（ω） …ＸＳ（ω））
・・・式５９
となる。
【０１４８】
一方で、伝達関数算出部１４は、ステップ７１において定義されたベクトルＶに対して、Ｎcorr点の高速フーリエ変換を適用する（ステップＳ．７４）。
なお、ベクトルＶの構成要素であるｖｊ（ｔ）（ｊ＝１，２，・・・，Ｓ）における時系列のデータ数がＮcorrに満たない場合、伝達関数算出部１４は、データの後ろに０を追加することにより高速フーリエ変換を実行する。このように０を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換（ＩＦＦＴ）を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。
【０１４９】
ステップＳ．７４においてＮcorr点で高速フーリエ変換（ＦＦＴ）すると、
【数４０】
を求めることができる。
【０１５０】
このようにして求められたＶ（ω）と転置行列Ａｔとの積を、周波数領域において算出する（周波数領域において畳み込み算出を実施する）と、
【数４１】
を求めることができる（ステップＳ．７５）。
【０１５１】
ここで、Ａ＊（ω）は、Ａ（ω）の共役転置行列を意味しており、
【数４２】
で表される。
【０１５２】
そして、伝達関数算出部１４は、上述した式により算出された（ＡｔＶ）（ω）の演算結果に対して逆フーリエ変換を適用することにより、時間領域におけるデータ（時系列データ）（ＡｔＶ）’に変換を行う（ステップＳ．７６）。
【０１５３】
伝達関数算出部１４は、この（ＡｔＶ）’の内のＳ個の逆フーリエ変換（ＩＦＦＴ）結果に対してそれぞれ０≦ｔ≦Ｌ−１の時間窓をかけることによりＡｔＶを求める（ステップＳ．７６）。このように、０≦ｔ≦Ｌ−１の範囲の時間窓を利用して、対応するデータのみを用いて計算を行うことによって、逆フーリエ変換によって求められた解が正確な値となるように（時間領域のみの演算により求められる解とに誤差が生じないように）することが可能となる。
【０１５４】
上述した図８および図９に示す処理手順に基づいて演算処理を行うことによって、伝達関数算出部１４における演算処理負担の低減を図ることが可能となる。
【０１５５】
次に、図１０〜図１２を用いて、長時間の連続した測定で得られた測定データの一部を用いることにより伝達関数を求める処理内容について説明する。
【０１５６】
上述したように、長時間の連続した測定で得られた測定データの一部を用いることにより伝達関数を求めるためには、
ＡｔＷ２Ｄ＝（ＡｔＷ２Ａ）Ｑ ・・・式５２
の解である、
Ｑ＝（ＡｔＷ２Ａ）−１（ＡｔＷ２Ｄ） ・・・式５３
を求めることが必要となる。
【０１５７】
伝達関数算出部１４は、式５２に示す方程式を最小二乗法で解くために、式５２の左辺をＰ、右辺の係数行列をＲとし、Ｐ＝ＲＱに基づいて
Ｐ＝Ａｔ（Ｗ２Ｄ）
Ｒ＝Ａｔ（Ｗ２Ａ）
とおいて、ＰおよびＲの計算を行い、この計算式に基づいて連立一次方程式Ｐ＝ＲＱを解くことにより、時間領域における行列Ｑ：Ｑ＝Ｒ−１Ｐを求める。
【０１５８】
なお、
Ｑ＝（ＡｔＷ２Ａ）−１（ＡｔＷ２Ｄ） ・・・式５３
の連立一次方程式を求める場合、Ｗが単位行列でない場合には、式５３の係数行列であるＡｔＷ２Ａの規則性が崩れてしまう。このため、伝達関数算出部１４では、式５３の演算を行う際に、ＡｔＷ２Ａを相関関数のベクトルとして記憶させて演算を行うことができない。
【０１５９】
簡単な数値例を挙げて説明すると、
ｘ１（ｔ）＝１，２，３，４，０・・・
ｘ２（ｔ）＝１１，１２，１３，１４，０・・・
ｗ（ｔ）＝０，０，０，１，１，１，０・・・
とし、
【数４３】
【数４４】
とすると、時間重みのない場合（つまりＷを考慮する必要がない場合）には、ＡｔＡが、
【数４５】
として求められる。このＡｔＡは、行列全体が対称行列であることに加え、（）内の個々のサブマトリックスの要素が斜め方向に同一値となる。このような行列の場合には、行列の全ての要素を個別に記憶しておく必要がないため、行列の規模が大きくなった場合であっても規則性を利用して記憶容量（メモリ容量）を節約することが可能である。
【０１６０】
しかしながら、時間重みのある場合には、
【数４６】
となり、サブマトリックスの要素が斜め方向に同一値を持つという規則性が崩れてしまう。このように行列の規模が大きくなった場合には、係数行列を記憶するために膨大な記憶容量（メモリ容量）が必要となるため、係数行列ＡｔＷ２Ａ（＝（ＷＡ）ｔＷＡ）を予め計算して変数として記憶することが困難となる。
【０１６１】
このため、本実施の形態２に係る伝達関数算出部１４では、伝達関数の計算処理において、係数行列ＡｔＷ２Ａそのものを求めることではなく、（ＡｔＷ２Ａ）ｖの形式の行列のベクトル積の演算結果を求める方法を用いている。（ＡｔＷ２Ａ）ｖをＡｔ（Ｗ２（Ａｖ））とする計算は、後述する図１２の計算手順と、図８および図９に示した一連の計算手順を用いることにより、畳み込み演算、時間重み、相関演算を順番に実施することで効率的に行うことが可能になる。従って、実施の形態２に示す数値計算方法を用いることによって、係数行列の記憶に関する問題を回避することができる。以下、上述した計算方法を用いて演算を行う手順について説明する。
【０１６２】
まず、伝達関数算出部１４は、図１０に示すように、Ｐ＝Ａｔ（Ｗ２Ｄ）の計算を行う（ステップＳ．８１）。
【０１６３】
具体的に、伝達関数算出部１４は、図１１に示すように、初期値としてｉに１を代入する（ステップＳ．９１）、ここでｉは、音響空間２におけるマイク１０の数を示しており、１に設定されたｉがマイク１０の設置数Ｍに該当するまでの間、伝達関数算出部１４はステップＳ．９２〜ステップＳ．９５の処理を繰り返し実行する。
【０１６４】
次に伝達関数算出部１４は、行列Ｄのｉ列目の要素ｄｉをとして抽出する（ステップＳ．９２）。そして、伝達関数算出部１４は、対角行列Ｗとｄｉとを用いて、Ｗ２ｄｉの計算を行う（ステップＳ．９３）。具体的に、伝達関数算出部１４は、ｗ（ｔ）ｙｉ（ｔ）の計算を行う。
【０１６５】
その後、伝達関数算出部１４は、行列Ｄのｉ列目の要素ｄｉを用いて求められたＷ２ｄｉの計算結果をベクトルＶとして、既に図９を用いて説明したように、転置行列ＡｔとベクトルＶとの演算方法を用いることにより、Ａｔ（Ｗ２ｄｉ）の計算を行う（ステップＳ．９４）。
【０１６６】
そして、伝達関数算出部１４は、全てのｉに対してＡｔ（Ｗ２ｄｉ）の計算を行ったか否か（つまり、ｉ＝Ｍであるか否か）の判断を行う（ステップＳ．９５）。全てのｉに対して計算を行っていないと判断した場合（ステップＳ．９５においてＮｏの場合）、伝達関数算出部１４は、ｉに１を加えた後に、上述したステップＳ．９２〜ステップＳ．９５の処理を繰り返し実行する（ステップＳ．９６）。
【０１６７】
一方で、全てのｉに対して計算を行ったと判断した場合（ステップＳ．９５においてＹｅｓの場合）、伝達関数算出部１４は、Ｄの全ての列に対応するｉのＡｔ（Ｗ２ｄｉ）を算出したものと判断し、つまり、Ａｔ（Ｗ２Ｄ）（＝Ｐ）を算出したものと判断して、処理を図１０のステップＳ．８１に戻す。
【０１６８】
次に、伝達関数算出部１４は、初期値としてｊに１を代入し、μに０を代入する（ステップＳ．８２）、ここでｊは、音響空間２におけるスピーカ９の数を示している。またμは、行列Ｑのｊ列目の値を算出するために実行された演算回数（反復回数）を示すパラメータであり、具体的には、ステップＳ．８３〜ステップＳ．８８の処理を繰り返し実行した回数を示している。
【０１６９】
次に伝達関数算出部１４は、初期設定を行う（ステップＳ．８３）。初期設定では、Ｐ，Ｑのｊ列目の値をｐ，ｑとおき、適当な初期値ｑ0を用いて、
ε0＝ｐ−Ａｔ（Ｗ２（Ａｑ0）） ・・・式６７
を求める。なお、εμは、反復回数μでの暫定的な解をｑμとしたときの、誤差ベクトルを示したものである。
【０１７０】
ε0＝ｐ−Ａｔ（Ｗ２（Ａｑ0）） ・・・式６７
におけるＡｔ（Ｗ２（Ａｑ0））の計算は、後述するステップＳ．８５に示す
Ｒσμ＝Ａｔ（Ｗ２（Ａσμ））
のσμに対して初期値としてｑ0を設定することによって求めることができる。
Ｒσμ＝Ａｔ（Ｗ２（Ａσμ））
の計算方法については、後述する。
なお、ｐ、ｑ、σμは、長さＬのベクトルをＳ個集めた長さＳＬのベクトルとして定義される。
【０１７１】
また、伝達関数算出部１４は、求められたε0を用いて
【数４７】
の初期値であるσ0の値
【数４８】
の値を算出する。
【０１７２】
さらに、伝達関数算出部１４は、共役勾配法の収束判定に用いられるパラメータをηとして設定する。ηは収束を判断するための閾値パラメータとしての定数であって、このパラメータに基づいて誤差が収束判定条件を満たす値となった場合（次述するステップＳ．８４においてＹｅｓの場合）、伝達関数算出部１４は、解が収束したと判断して反復計算を終了する。
【０１７３】
そして、伝達関数算出部１４は、収束条件の判断処理（ステップＳ．８４）を行う。収束条件の判断式は、
【数４９】
で表される。伝達関数算出部１４は、この判断式に基づいて解の収束を判断する。
【０１７４】
収束条件を満たさなかった場合（ステップＳ．８４においてＮｏの場合）、伝達関数算出部１４は、Ｒσμの計算を行う（ステップＳ．８５）。
【０１７５】
図１２は、Ｒσμの計算方法、つまり、Ａｔ（Ｗ２（Ａσμ））を計算する方法を示したフローチャートである。
【０１７６】
伝達関数算出部１４は、まず、図８に示した行列ＡとベクトルＶとの積の計算方法において、ベクトルＶにベクトルσμを用いることにより、Ａσμを計算する（ステップＳ．１０１）。次に、伝達関数算出部１４は、対角行列Ｗの積であるＷ２にＡσμを積算したＷ２（Ａσμ）の計算を行う（ステップＳ．１０２）。
【０１７７】
そして、伝達関数算出部１４は、図９に示した転置行列ＡｔとベクトルＶとの積の計算方法において、ベクトルＶにベクトルＷ２（Ａσμ）を用いることにより、Ａｔ（Ｗ２（Ａσμ））を計算する（ステップＳ．１０３）。
【０１７８】
このように、Ａｔ（Ｗ２（Ａσμ））の計算を行う場合にも、伝達関数算出部１４において、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部１４における処理負担を軽減させることができる。
【０１７９】
そして、伝達関数算出部１４は、求められたＲσμの値と、εμ、σμとを用いて、
【数５０】
【数５１】
【数５２】
を算出し、さらに、μに１を加える（μ＝μ＋１）処理を行う（ステップＳ．８６）。
【０１８０】
そして、伝達関数算出部１４は、ステップＳ．８６において算出されたεμとＰとを用いて、収束条件の判断処理（ステップＳ．８４）を繰り返し実行する。
【０１８１】
収束条件を満たす場合（ステップＳ．８４においてＹｅｓの場合）、伝達関数算出部１４は、ステップＳ．８６において算出されたｑμの値を、行列Ｑにおけるｊ列目の解（伝達関数）とする（ステップＳ．８７）。
【０１８２】
その後、伝達関数算出部１４は、全ての列ｊについて、つまり全てのスピーカに対応する各マイクの伝達関数（伝達関数マトリックス（行列Ｑ）における全ての伝達関数）を求めたか否かを判断し（ステップＳ．８８）、全ての列において伝達関数を求めていない場合（ステップＳ．８８においてＮｏの場合）には、ｊに１を加え、μに０を代入した後に（ステップＳ．８９）、処理をステップＳ．８３に移行して上述した処理を繰り返し実行する。
【０１８３】
全ての列において伝達関数を求めた場合（ステップＳ．８８においてＹｅｓの場合）、伝達関数算出部１４は、全ての列の伝達関数の値を求めたものと判断して、伝達関数の計算処理（伝達関数マトリックス（行列Ｑ）の算出処理）を終了する。
【０１８４】
このように、実施の形態２に係る伝達関数算出部１４では、長時間連続して同時に複数のスピーカ９から評価音を出力すると共に、スピーカ９から出力された評価音を連続してマイク１０で同時に測定した場合であっても、時間重みを適用して伝達関数を算出することができる。このため、連続的に測定位置が変更される場合であっても時間変化に応じた適切な伝達関数を求めることが可能となる。
【０１８５】
なお、実施の形態２では、上述したように、長時間連続して同時に複数のスピーカ９から評価音を同時に出力すると共に、スピーカ９から出力された評価音を連続してマイク１０で同時に測定した場合について説明を行ったが、実施の形態２に係る発明は、複数のスピーカ９から同時に評価音を出力する場合に限定されるものではない。例えば、従来の方法を用いて伝達関数を求める場合のように、各スピーカから順番に評価音を出力することにより該当するスピーカとマイクとに対応する伝達関数を算出する場合であっても、各評価音の出力を長時間連続して行うことにより、該当するスピーカと各マイクとに対応する伝達関数の算出を、時間重みを適用して求めることが可能となる。
【０１８６】
[実施の形態３]
実施の形態１に示す方法を用いることにより、複数のスピーカから同時に評価音を出力すると共に、複数のマイクから同時に評価音の収録を行うことによって、一度の測定でそれぞれのマイクとスピーカとに対応させた伝達関数を求めることが可能となる。また、実施の形態２に示す方法を用いることにより、時間重みを適用して伝達関数を算出することができるので、連続的に測定位置が変更される場合であっても時間変化に応じた適切な伝達関数を求めることが可能となる。
【０１８７】
しかしながら、音響空間２にノイズが存在する場合、マイク１０で測定された収録音（出力信号、測定データ）Ｙには、評価音（入力信号）Ｘの成分（寄与成分）に加えて、ノイズ成分が含まれることになる。
【０１８８】
例えば、図１３に示すように、測定データＹ（ｚ）を、入力信号Ｘ（ｚ）と、伝達関数Ｈ（ｚ）と、ノイズ成分Ｎ（ｚ）とを用いて表すと、
【数５３】
と表される。
【０１８９】
このように、測定データＹ（ｚ）に、入力信号Ｘ（ｚ）の寄与成分と、ノイズ成分Ｎ（ｚ）とが含まれている場合において、ノイズの影響を低減させるためには、ノイズ信号のパワーに対して入力信号のパワーを相対的に大きくすることにより、測定Ｓ／Ｎ（Signal to Noise ratio）を大きくすることが必要となる。
【０１９０】
ところで、マイク１０を使用して、音響空間２における音響特性を測定する場合には、低域周波数であるほど、レベルの高い環境ノイズが含まれる傾向があり、高い周波数に比例して低域周波数での測定Ｓ／Ｎが劣化する傾向がある。このため、スピーカ９より出力される評価音（入力信号）の低域周波数成分のレベルを大きくすることにより、収録音（測定データ）に含まれる低い周波数成分のＳ／Ｎを改善することが考えられる。
【０１９１】
従って、本実施の形態３に係る伝達関数算出部１４では、入力信号として低域周波数成分のレベルの高い信号を用いることによって、測定データに含まれる低い周波数成分のＳ／Ｎを改善する方法を採用する。
【０１９２】
図１４（ａ）（ｂ）は、入力信号の周波数特性をそれぞれ変えた場合における入力信号の信号パワー値とマイクで集音されるノイズの信号パワーとの対応を示した図である。（ａ）は、入力信号の周波数特性が平坦な場合（つまり、低域周波数成分の信号パワーの補強を行わない場合）を示しており、（ｂ）は、入力信号の周波数特性を変えることにより、低域周波数成分の信号パワーが大きくなるようにした場合を示している。また、（ａ）（ｂ）には、周波数帯域に応じて信号パワーが変化するノイズの状態を対比して示している。
【０１９３】
図１４（ａ）に示す入力信号とノイズとの信号のパワー比をとると、低域周波数帯域において、ノイズのパワー値と入力信号のパワー値とが近い値となるため、Ｓ／Ｎが小さな値になってしまう。一方で、図１４（ｂ）に示す入力信号とノイズとの信号のパワー比をとると、入力信号の信号パワーが周波数の減少にともない増加する状態を示しているため、低域周波数帯域であっても、ノイズのパワー値と入力信号のパワー値との間隔が一定に保たれることになり、Ｓ／Ｎの減少を抑制することが可能となる。
【０１９４】
しかしながら、低域周波数成分の信号パワーが大きな入力信号を用いると、自己相関特性が劣化して、連立一次方程式の係数行列（ＡｔＡ）の性質が悪くなるため、連立方程式を反復解法で解くときの収束性が劣化してしまう。このように、解の収束性が悪くなると、反復解法に必要な収束回数が多く必要になり、有限の反復回数で得られる解の精度が劣化してしまう。そこで、実施の形態３に係る伝達関数算出部１４では、収束性を改善するために、入力信号（評価音）に対して前処理を行うことにより、収束性の劣化を回避して解の精度を保つ方法が用いられている。
【０１９５】
下記の式７５は、低域成分が高いパワーとなる試験信号を用いる場合において、収束性の改善を実現することが可能な数式を示している。実施の形態３の入力信号Ｘには、予めマイク１０に加えられるノイズの周波数特性を考慮して周波数特性が平坦ではない信号を使用しているものとする。このため、
Ｙｔ（ｚ）＝Ｘｔ（ｚ）Ｈｔ（ｚ）
を解く代わりに、この式の両辺にｇ（ｚ）を畳み込む前処理を実施して、
（ｇ（ｚ）Ｙｔ（ｚ））＝（ｇ（ｚ）Ｘｔ（ｚ））Ｈｔ（ｚ） ・・・式７５
を解く。
【０１９６】
ここで、
ｇ（ｚ）Ｘｔ（ｚ）
＝（ｇ（ｚ）Ｘ１（ｚ） ｇ（ｚ）Ｘ２（ｚ） … ｇ（ｚ）ＸＳ（ｚ））
・・・ 式７６
ｇ（ｚ）Ｙｔ（ｚ）
＝（ｇ（ｚ）Ｙ１（ｚ） ｇ（ｚ）Ｙ２（ｚ） … ｇ（ｚ）ＹＭ（ｚ））
・・・ 式７７
である。
【０１９７】
式７５を解くことは、これを代数方程式で表した
ＧＤ＝（ＧＡ）Ｑ ・・・ 式７８
を解いて、
Ｑ＝（（ＧＡ）ｔＧＡ）−１（（ＧＡ）ｔＤ）・・・ 式７９
を得ることと同じである。
【０１９８】
ここで、Ｇはｇ（ｚ）の畳み込みマトリックスを示しており、
【数５４】
で表される。
【０１９９】
式７９に示す式の係数行列（ＧＡ）ｔ（ＧＡ）は、
（ＧＡ）ｔ（ＧＡ）＝Ａｔ（ＧｔＧ）Ａ
が単位行列に近づくようなＧを与えることで、解の収束性を改善することができる。例えば、前処理に用いるフィルタとして、入力信号Ｘの周波数特性とほぼ逆の周波数特性を持つ最小位相フィルタを利用することができる。
【０２００】
一般的に、周波数特性がフラットな信号の相関特性は、周波数特性がフラットでない信号より鋭い相関ピークを持ち、相関特性がよい。また最小位相フィルタを用いることにより前処理フィルタの段数を必要最小限にし、前処理に必要な計算コストの低減を図ることが可能となる。
【０２０１】
以上、本発明に係る音場測定システムについて、図面を用いて詳細に説明を行ったが、本発明に係る音場測定システムは、上述したものに限定されるわけではない。当業者であれば、特許請求の範囲に記載された範疇内において、各種の変更例または修正例に想到しうることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
【０２０２】
なお、上述した音場測定システムを用いて、多数のスピーカから一度に評価音を出力し、多数のマイクにおいて同時に評価音の収録を行うことによって、各スピーカから各マイクに対応する伝達関数を一度に算出する方法が、従来のように各スピーカより順番に評価音を出力し、多数のマイクで評価音の収録を行うことにより伝達関数を算出する場合に比べて、測定時間の短縮化を図ることが可能であることを、具体的に説明する。
【０２０３】
スピーカ（音源）の数をＳとし、伝達関数を測定するための応答時間がＬで打ち切られると仮定する。従来のように、各スピーカより順番に評価音を出力し、多数のマイクにおいて同時に評価音の収録を行うことにより伝達関数を算出する場合において、１つのスピーカ（音源）より出力される評価音をマイクで収録して測定を行う時間は、概ね
α×Ｌ
程度と見積もることができる。
【０２０４】
ここでαは測定条件（積分処理に要する時間や応答時間等）を考慮した定数であるが、測定データの積分時間や平均回数を大きくしてノイズ環境下での測定Ｓ／Ｎを確保するために、αは１より大きな値とする。全てのスピーカから出力される評価音を測定するためには、
α×Ｌ×Ｓ
程度の時間が必要とされる。
【０２０５】
一方で、本実施の形態１〜実施の形態３に示した音場測定システムを用いる場合における測定時間を求めるためには、
１）ノイズ環境下において測定Ｓ／Ｎを確保するために必要な測定時間
２）方程式の係数行列が逆行列を持つために原理上必要となる入力信号の長さ
を考慮する必要がある。
【０２０６】
まず１つ目の「ノイズ環境下での測定Ｓ／Ｎを確保するために必要な測定時間」について考える。本実施の形態に示した測定方式は、他のスピーカ（音源）から再生される音が測定ノイズとして加算されることがない方式である。このため、測定Ｓ／Ｎを確保するために必要な測定時間は、本実施の形態に示したように複数のスピーカから同時に評価音を出力する方式であっても、従来の方式により１つのスピーカで出力される音を測定する方式であっても、一度の測定で必要とされる測定時間は同等である。従って、測定時間は音源の数に依存せず、
α×Ｌ
と見積もることができる。
【０２０７】
次に２つ目の「方程式の係数行列が逆行列を持つために原理上必要となる入力信号の長さ」は
β×Ｌ×Ｓ
と見積もることができる。
なお、βの値は、入力信号（評価音）の長さとして演算処理上必要な値を示しており、応答時間Ｌに比例した値となる。このβの値は、測定環境のノイズの大きさには無関係であり、β＝２程度で十分であると判断できる。
【０２０８】
以上２つの条件を考慮すると、本実施の形態に示した測定方法を用いる際に必要とされる測定時間は、
ＭＡＸ（α×Ｌ，β×Ｌ×Ｓ）
（：（α×Ｌ）と（β×Ｌ×Ｓ）とを比較して大きな値を有効とする）
と見積もることができる。
【０２０９】
従って、従来方式を用いて測定する時間と本提案方式を用いて測定する時間との比は以下のようになる。
【数５５】
ノイズ環境下においては、α＞βとなるため、本提案方式の測定時間の方が従来方式の測定時間より短くなる。
【符号の説明】
【０２１０】
１ …音場測定システム
２ …音響空間
３ …スピーカ部
４ …音声出力装置
５ …コンピュータ
７ …マイクアレイ
８ …音声入力装置
９ …スピーカ
１０ …マイク
１２ …評価音生成部（評価音生成手段）
１３ …収録音記録部
１４ …伝達関数算出部（伝達関数算出手段）
１５ …伝達関数記録部
１６ …制御部
【特許請求の範囲】
【請求項１】
特定の音響空間に設定される複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録することによって、各スピーカから各マイクに対応する伝達関数を算出するための音場測定システムであって、
前記スピーカより出力するための評価音の生成を行う評価音生成手段と、
前記マイクを用いて収録された収録音と前記評価音生成手段により生成された評価音とに基づいて、各スピーカから各マイクに対応する伝達関数を算出する伝達関数算出手段と
を備え、
前記伝達関数算出手段は、前記収録音と前記評価音とに基づいて、時間領域における反復解法を適用することにより、前記伝達関数の算出処理を実行し、前記反復解法の適用において用いる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させる演算方法を用い、
前記評価音生成手段は、前記評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第１行列と、該第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成する
ことを特徴とする音場測定システム。
【請求項２】
前記伝達関数算出手段は、前記伝達関数の算出処理において、前記伝達関数を求めるための計算式に対して時間窓を構成する対角行列を掛け合わせ、該対角行列の要素を変更することにより、前記複数のスピーカから連続して評価音を出力すると共に前記複数のマイクから連続して収録音の収録を行った場合における任意の時間範囲の伝達関数を算出することを
を特徴とする請求項１に記載の音場測定システム。
【請求項３】
前記評価音生成手段は、前記評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、前記マイクにより収録される収録音の低域周波数帯域におけるＳ／Ｎの改善を図ること
を特徴とする請求項１又は請求項２に記載の音場測定システム。
【請求項４】
前記伝達関数算出手段は、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させて求められた前記相関演算処理または前記畳み込み演算処理の演算結果のうち、前記時間領域のみで前記相関演算処理または前記畳み込み演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、前記時間領域における前記伝達関数の算出処理を行う
ことを特徴とする請求項１乃至請求項３のいずれか１項に記載の音場測定システム。
【請求項５】
前記誤差が生じない範囲は、前記評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定されること
を特徴とする請求項４に記載の音場測定システム。
【請求項１】
特定の音響空間に設定される複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録することによって、各スピーカから各マイクに対応する伝達関数を算出するための音場測定システムであって、
前記スピーカより出力するための評価音の生成を行う評価音生成手段と、
前記マイクを用いて収録された収録音と前記評価音生成手段により生成された評価音とに基づいて、各スピーカから各マイクに対応する伝達関数を算出する伝達関数算出手段と
を備え、
前記伝達関数算出手段は、前記収録音と前記評価音とに基づいて、時間領域における反復解法を適用することにより、前記伝達関数の算出処理を実行し、前記反復解法の適用において用いる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させる演算方法を用い、
前記評価音生成手段は、前記評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第１行列と、該第１行列の転置行列を成す第２行列との積の逆行列が存在し得るように評価音を生成する
ことを特徴とする音場測定システム。
【請求項２】
前記伝達関数算出手段は、前記伝達関数の算出処理において、前記伝達関数を求めるための計算式に対して時間窓を構成する対角行列を掛け合わせ、該対角行列の要素を変更することにより、前記複数のスピーカから連続して評価音を出力すると共に前記複数のマイクから連続して収録音の収録を行った場合における任意の時間範囲の伝達関数を算出することを
を特徴とする請求項１に記載の音場測定システム。
【請求項３】
前記評価音生成手段は、前記評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、前記マイクにより収録される収録音の低域周波数帯域におけるＳ／Ｎの改善を図ること
を特徴とする請求項１又は請求項２に記載の音場測定システム。
【請求項４】
前記伝達関数算出手段は、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させて求められた前記相関演算処理または前記畳み込み演算処理の演算結果のうち、前記時間領域のみで前記相関演算処理または前記畳み込み演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、前記時間領域における前記伝達関数の算出処理を行う
ことを特徴とする請求項１乃至請求項３のいずれか１項に記載の音場測定システム。
【請求項５】
前記誤差が生じない範囲は、前記評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定されること
を特徴とする請求項４に記載の音場測定システム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【公開番号】特開２０１０−１６４８６３（Ｐ２０１０−１６４８６３Ａ）
【公開日】平成２２年７月２９日（２０１０．７．２９）
【国際特許分類】
【出願番号】特願２００９−８３４０（Ｐ２００９−８３４０）
【出願日】平成２１年１月１９日（２００９．１．１９）
【出願人】（０００００１４８７）クラリオン株式会社 (1,722)
【公開日】平成２２年７月２９日（２０１０．７．２９）
【国際特許分類】
【出願日】平成２１年１月１９日（２００９．１．１９）
【出願人】（０００００１４８７）クラリオン株式会社 (1,722)
[ Back to top ]