JP5220922B2 - 音場の再構成 - Google Patents

Description

本発明は音場（ａｃｏｕｓｔｉｃ ｆｉｅｌｄ）の再構成（ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ）に関する。

近距離場音響ホログラフィ（ＮＡＨ：Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ）は、音波放射の３Ｄ視覚化、および音源（ｓｏｕｎｄ ｓｏｕｒｃｅ）に近いある表面を横切る測定に基づく雑音源の精確な位置同定のための非常に有用なツールである。エバネッセント波（ｅｖａｎｅｓｃｅｎｔ ｗａｖｅ）成分をも再構成することができる近距離場音響ホログラフィの能力は、非常に高い空間分解能を保証する。

知られているある近距離場音響ホログラフィ法は、分離可能な座標系内のある水準面（ｌｅｖｅｌ ｓｕｒｆａｃｅ）を横切る規則格子（ｒｅｇｕｌａｒ ｇｒｉｄ）測定に基づき、この測定は、空間離散フーリエ変換（ＤＦＴ：Ｄｉｓｃｒｅｔｅ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）によってＮＡＨ計算を実行することを可能にする。例えば、Ｅ．Ｇ．Ｗｉｌｌｉａｍｓ、Ｊ．Ｄ．ＭａｙｎａｒｄおよびＥ．Ｊ．Ｓｋｕｄｒｚｙｋ、「Ｓｏｕｎｄ ｓｏｕｒｃｅ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｕｓｉｎｇ ａ ｍｉｃｒｏｐｈｏｎｅ ａｒｒａｙ」、Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．６８、３４０〜３４４、（１９８０）を参照されたい。ＤＦＴを使用するため処理は非常に高速だが、ＤＦＴを使用することの副作用には、測定エリアが音圧の高いエリアを完全にカバーしない場合の重大な空間ウィンドウイング効果（ｓｐａｔｉａｌ ｗｉｎｄｏｗｉｎｇ ｅｆｆｅｃｔ）が含まれる。場合によっては、測定エリアに関するこの要件を満たすことができず、多くの場合には、必要なサイズが極端に大きくなる。

空間ウィンドウイング効果を低減させ、同時に、複雑さおよび計算要求量の増大と引替えにＤＦＴ空間処理をなおも維持する一連の技法が提案されている。例えば、Ｊ．Ｈａｌｄ、「Ｒｅｄｕｃｔｉｏｎ ｏｆ ｓｐａｔｉａｌ ｗｉｎｄｏｗｉｎｇ ｅｆｆｅｃｔｓ ｉｎ ａｃｏｕｓｔｉｃａｌ ｈｏｌｏｇｒａｐｈｙ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｉｎｔｅｒ−Ｎｏｉｓｅ、１９９４を参照されたい。一般に、最初に反復的手順を使用して測定音圧を測定エリアの外側へ外挿し、続いて、その拡張されたデータ・ウィンドウに、ＤＦＴに基づくホログラフィ法を適用する。

空間ＤＦＴの使用を回避し、必要な測定エリアを低減させることを追求する他の方法が提案されている。

そのような１つの方法が、球面波関数に関する音場の局所モデルを使用するヘルムホルツ方程式最小２乗（ＨＥＬＳ：Ｈｅｌｍｈｏｌｔｚ’Ｅｑｕａｔｉｏｎ Ｌｅａｓｔ Ｓｑｕａｒｅｓ）法である。例えば、米国特許第６，６１５，１４３号、Ｚ．ＷａｎｇおよびＳ．Ｆ．Ｗｕ、「Ｈｅｌｍｈｏｌｔｚ ｅｑｕａｔｉｏｎ−ｌｅａｓｔ−ｓｑｕａｒｅｓ ｍｅｔｈｏｄ ｆｏｒ ｒｅｃｏｎｓｔｒｕｃｔｉｎｇ ｔｈｅ ａｃｏｕｓｔｉｃ ｐｒｅｓｓｕｒｅ ｆｉｅｌｄ」、Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ、１０２（４）、２０２０〜２０３２（１９９７）、またはＳ．Ｆ．Ｗｕ、「Ｏｎ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｏｆ ａｃｏｕｓｔｉｃ ｆｉｅｌｄｓ ｕｓｉｎｇ ｔｈｅ Ｈｅｌｍｈｏｌｔｚ ｅｑｕａｔｉｏｎ−ｌｅａｓｔ−ｓｑｕａｒｅｓ ｍｅｔｈｏｄ」、Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ、１０７、２５１１〜２５２２（２０００）を参照されたい。しかしながら、共通の原点を有する球面波関数だけを使用して音場を表現するため、音源面（ｓｏｕｒｃｅ ｓｕｒｆａｃｅ）が球面でなく、音源面の中心が同じ原点にない場合には、音源面における音場の再構成に誤差が導入される。この先行技術の方法の他の欠点は、この推定問題が劣決定（ｕｎｄｅｒ−ｄｅｔｅｒｍｉｎｅｄ）でない必要があることであり、このことは、測定位置の数が、場の表現（ｆｉｅｌｄ ｒｅｐｒｅｓｅｎｔａｔｉｏｎ）において使用される球面波関数の数に等しいか、またはそれよりも多くなければならないことを意味する。これにより、十分に正確なモデルを得るために多数の測定位置が必要となることがある。第３の欠点は、モデル領域においてより強く減衰する関数が、同じ領域においてより低い振幅にスケーリング（ｓｃａｌｉｎｇ）されるような方式では、波動関数のスケーリングが適用されないことである。このようなスケーリングがなされないため、チホノフ正則化（Ｔｉｋｈｏｎｏｖ ｒｅｇｕｌａｒｉｚａｔｉｏｎ）のような伝統的な正則化法が適切に機能しない。その代わりに、上記の先行技術の方法は、正則化されていない最小２乗解と組み合わせた球面波展開の最適な打切り（ｔｒｕｎｃａｔｉｏｎ）のための計算コストが高い反復的な検索を使用する。

以前に提案された別の方法が、Ｒ．ＳｔｅｉｎｅｒおよびＪ．Ｈａｌｄ、「Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ ｗｉｔｈｏｕｔ ｔｈｅ ｅｒｒｏｒｓ ａｎｄ ｌｉｍｉｔａｔｉｏｎｓ ｃａｕｓｅｄ ｂｙ ｔｈｅ ｕｓｅ ｏｆ ｓｐａｔｉａｌ ＤＦＴ」、Ｉｎｔｅｒｎ．Ｊ．Ａｃｏｕｓｔ．Ｖｉｂ．６、８３〜８９（２００１）において提唱されている統計的最適化近距離場音響ホログラフィ（ＳＯＮＡＨ：Ｓｔａｔｉｓｔｉｃａｌｌｙ Ｏｐｔｉｍｉｚｅｄ Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ）法である。この先行技術の方法は平面波関数に基づき、全ての伝搬波および重み付けした一組のエバネッセント波が最適な平均正確さで射影されるような方法で定義された伝達行列（ｔｒａｎｓｆｅｒ ｍａｔｒｉｘ）を使用することによって、測定表面の近くのマッピング表面上の音響量を計算する。この伝達行列は、１）全ての測定位置対間の自己および相互相関の自己相関行列、ならびに２）それぞれの測定点から全ての再構成位置への相互相関の相互相関行列から得られる。この自己および相互相関は、平面伝搬波および平面エバネッセント波の領域（ｄｏｍａｉｎ）にある。

特定の平面測定ジオメトリ（ｇｅｏｍｅｔｒｙ）へのこの枠組みの適用が、Ｊ．Ｈａｌｄ、「Ｐａｔｃｈ ｎｅａｒ−ｆｉｅｌｄ ａｃｏｕｓｔｉｃａｌ ｈｏｌｏｇｒａｐｈｙ ｕｓｉｎｇ ａ ｎｅｗ ｓｔａｔｉｓｔｉｃａｌｌｙ ｏｐｔｉｍａｌ ｍｅｔｈｏｄ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｉｎｔｅｒ−Ｎｏｉｓｅ ２００３、およびＪ．Ｈａｌｄ、「Ｐａｔｃｈ ｈｏｌｏｇｒａｐｈｙ ｉｎ ｃａｂｉｎ ｅｎｖｉｒｏｎｍｅｎｔｓ ｕｓｉｎｇ ａ ｔｗｏ−ｌａｙｅｒ ｈａｎｄｈｅｌｄ ａｒｒａｙ ａｎｄ ａｎ ｅｘｔｅｎｄｅｄ ＳＯＮＡＨ ａｌｇｏｒｉｔｈｍ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｅｕｒｏｎｏｉｓｅ ２００６に記載されている。

これらの先行技術の文献は、測定平面に垂直な圧力および粒子速度成分の計算を開示している。しかしながら、これらの先行技術の文献に開示されている相関式の数値評価は計算効率が悪いため、この方法は計算コストが高い。

米国特許第６，６１５，１４３号

Ｅ．Ｇ．Ｗｉｌｌｉａｍｓ、Ｊ．Ｄ．ＭａｙｎａｒｄおよびＥ．Ｊ．Ｓｋｕｄｒｚｙｋ、「Ｓｏｕｎｄ ｓｏｕｒｃｅ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｕｓｉｎｇ ａ ｍｉｃｒｏｐｈｏｎｅ ａｒｒａｙ」、Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．６８、３４０〜３４４、（１９８０） Ｊ．Ｈａｌｄ、「Ｒｅｄｕｃｔｉｏｎ ｏｆ ｓｐａｔｉａｌ ｗｉｎｄｏｗｉｎｇ ｅｆｆｅｃｔｓ ｉｎ ａｃｏｕｓｔｉｃａｌ ｈｏｌｏｇｒａｐｈｙ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｉｎｔｅｒ−Ｎｏｉｓｅ、１９９４ Ｚ．ＷａｎｇおよびＳ．Ｆ．Ｗｕ、「Ｈｅｌｍｈｏｌｔｚ ｅｑｕａｔｉｏｎ−ｌｅａｓｔ−ｓｑｕａｒｅｓ ｍｅｔｈｏｄ ｆｏｒ ｒｅｃｏｎｓｔｒｕｃｔｉｎｇ ｔｈｅ ａｃｏｕｓｔｉｃ ｐｒｅｓｓｕｒｅ ｆｉｅｌｄ」、Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ、１０２（４）、２０２０〜２０３２（１９９７） Ｓ．Ｆ．Ｗｕ、「Ｏｎ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｏｆ ａｃｏｕｓｔｉｃ ｆｉｅｌｄｓ ｕｓｉｎｇ ｔｈｅ Ｈｅｌｍｈｏｌｔｚ ｅｑｕａｔｉｏｎ−ｌｅａｓｔ−ｓｑｕａｒｅｓ ｍｅｔｈｏｄ」、Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ、１０７、２５１１〜２５２２（２０００） Ｒ．ＳｔｅｉｎｅｒおよびＪ．Ｈａｌｄ、「Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ ｗｉｔｈｏｕｔ ｔｈｅ ｅｒｒｏｒｓ ａｎｄ ｌｉｍｉｔａｔｉｏｎｓ ｃａｕｓｅｄ ｂｙ ｔｈｅ ｕｓｅ ｏｆ ｓｐａｔｉａｌ ＤＦＴ」、Ｉｎｔｅｒｎ．Ｊ．Ａｃｏｕｓｔ．Ｖｉｂ．６、８３〜８９（２００１） Ｊ．Ｈａｌｄ、「Ｐａｔｃｈ ｎｅａｒ−ｆｉｅｌｄ ａｃｏｕｓｔｉｃａｌ ｈｏｌｏｇｒａｐｈｙ ｕｓｉｎｇ ａ ｎｅｗ ｓｔａｔｉｓｔｉｃａｌｌｙ ｏｐｔｉｍａｌ ｍｅｔｈｏｄ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｉｎｔｅｒ−Ｎｏｉｓｅ ２００３ Ｊ．Ｈａｌｄ、「Ｐａｔｃｈ ｈｏｌｏｇｒａｐｈｙ ｉｎ ｃａｂｉｎ ｅｎｖｉｒｏｎｍｅｎｔｓ ｕｓｉｎｇ ａ ｔｗｏ−ｌａｙｅｒ ｈａｎｄｈｅｌｄ ａｒｒａｙ ａｎｄ ａｎ ｅｘｔｅｎｄｅｄ ＳＯＮＡＨ ａｌｇｏｒｉｔｈｍ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｅｕｒｏｎｏｉｓｅ ２００６ Ｊ．ＧｏｍｅｓおよびＰ．Ｃ．Ｈａｎｓｅｎ、「Ａ ｓｔｕｄｙ ｏｎ ｒｅｇｕｌａｒｉｚａｔｉｏｎ ｐａｒａｍｅｔｅｒ ｃｈｏｉｃｅ ｉｎ Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ａｃｏｕｓｔｉｃｓ’０８、２００８ Ｊ．Ｈａｌｄ、「Ａｒｒａｙ ｄｅｓｉｇｎｓ ｏｐｔｉｍｉｚｅｄ ｆｏｒ ｂｏｔｈ ｌｏｗ−ｆｒｅｑｕｅｎｃｙ ＮＡＨ ａｎｄ ｈｉｇｈ−ｆｒｅｑｕｅｎｃｙ ｂｅａｍｆｏｒｍｉｎｇ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｉｎｔｅｒ−Ｎｏｉｓｅ ２００４

本明細書では、少なくとも１つの音源が発生させた音場を再構成する方法を開示する。この方法の実施形態は、
− 一組の測定位置のそれぞれで測定した少なくとも１つの第１の音響量の一組の測定値を受け取ること、
− 前記測定位置のうちの第１の測定位置における少なくとも一組の平面伝搬波および平面エバネッセント波と第２の位置における少なくとも一組の平面波との相関関係をそれぞれが指示する一組の相関関数のそれぞれを計算すること、
− 少なくとも計算した相関関数と一組の測定値とから、標的位置における再構成された音場を指示する第２の音響量を計算すること
を含み、さらに、
相関関数に対するそれぞれの寄与（ｃｏｎｔｒｉｂｕｔｉｏｎ）を補間する一組の補間関数（ｉｎｔｅｒｐｏｌａｔｉｏｎ ｆｕｎｃｔｉｏｎ）の表現を記憶することを含み、それぞれの寄与が、入力パラメータが２つ以下の関数であり、一組のそれぞれの相関関数を計算することが、記憶した表現から得た値の１次結合としてそれぞれの相関関数を計算することを含む。

本発明の発明者は、本明細書に記載した方法を、一般的なジオメトリおよびいくつかの異なる音響量の計算に、一様でかつ計算効率がよい方式で適用することができることを理解した。多様な測定ジオメトリおよび音響量に対する相関関数を、入力値が２つ以下のそれぞれの補間関数によって少なくとも局所的に近似することができる１次元および／または２次元関数の１次結合として表現することができることも分かった。具体的には、補間関数が２次元以下の関数であるため、計算効率がよく、それにも関わらず計算が正確な補間が達成され、補間関数を記憶するのに必要な記憶容量がわずかですむ。したがって、数値計算の基礎をなす諸式が、補間と組み合わされた近似関数またはテーブルを使用することによる効率的な計算によく適した形を有する。本明細書に記載した方法の実施形態は、先行技術の方法の約１０倍の速度で実行することができる。

音源は、音響放射を放出し、または反射する任意の物体とすることができる。物体の形状は任意であり、音は、あらゆる種類の音、例えば雑音、可聴音、不可聴音（超音波、超低周波音など）など、あるいはこれらの音の組合せとすることができる。音源（１つまたは複数）を源のない領域（ｓｏｕｒｃｅ−ｆｒｅｅ ｒｅｇｉｏｎ）から分離する音源面によって、１つまたは複数の音源を含む領域を定義することができる。音源面は一般に、音源（１つまたは複数）を含む物理的物体の表面とすることができる。

この方法の実施形態は、一組の平面伝搬波関数および平面エバネッセント波関数をそれぞれが表現する１つまたは複数の仮想源平面を定義することを含み、これらの平面波はそれぞれ波数を有し、仮想源平面に対してある方向に伝搬しかつ／または減衰する（ｅｖａｎｅｓｃｅｎｔ）。一般に、仮想源平面に関して、源から、測定位置および再構成位置を含む領域（この領域を測定領域とも呼ぶ）内へ向かう法線を定義することができる。対応する平面伝搬波は一般に、垂直方向から垂直方向に対して９０度までの全ての伝搬方向に伝搬することができる。エバネッセント波は、この垂直方向において、すなわち測定領域内へ向かって指数関数的に減衰する。エバネッセント波は、異なる減衰率および異なる接線伝搬方向を有する連続スペクトルを形成することができる。

例えば、仮想源平面は、音源を含む物体の表面の一部分に対して、すなわち音源面に対して定義することができる。例えば、音源面から短い距離を置いて離隔した、物体の内部に位置する表面として、仮想源平面を定義することができる。

測定位置における音響量は、適当な音響測定装置、例えばマイクロホン、水中マイクロホン、圧力勾配変換器、粒子速度変換器などによって、またはこれらの装置の組合せによって測定することができる。ある実施形態では、この測定が、音響測定装置のアレイ（ａｒｒａｙ）、例えば規則格子または不規則格子（ｉｒｒｅｇｕｌａｒ ｇｒｉｄ）、例えば２次元または３次元格子として配置された一組の音響測定装置によって実行される。測定する第１の音響量は、音圧、音圧勾配、粒子速度および／または他の同種の音響量とすることができる。

再構成された音場は、音圧、音の強さ、粒子速度などの第２の音響量の空間分布を指示する適当な任意のデータ・セットによって表現することができる。第１の音響量と第２の音響量は、同じ音響量でもまたは異なる音響量でもよい。このデータ・セットを、任意の表面、例えば解析対象の音／雑音放出物体の表面またはその物体の近くの表面の実際の表面ジオメトリ上に、マップ、例えば第２の音響量のマップとして直接に表現することができる。

硬い表面上の粒子速度は、その表面自体の実際の振動に密接に対応するため、構造モデルとの相関に対してそれらの結果を直接に使用することができる。ある実施形態では、音場が、共形マップ（ｃｏｎｆｏｒｍａｌ ｍａｐ）によって表現される。ある表面、一般に音源面または音源面の近くの表面上において音場パラメータを再構成し、マップすることができる。あるいは、またはこれに加えて、再構成された音場の特性の他の表現、例えば波動関数表現の３Ｄ妥当性領域（３Ｄ ｒｅｇｉｏｎ ｏｆ ｖａｌｉｄｉｔｙ）内の単一位置におけるスペクトルを生成することもできる。

本明細書に記載した方法の実施形態では、計算点における圧力または粒子速度の再構成が、１）全ての測定位置対間の自己および相互相関の自己相関行列の計算、ならびに２）それぞれの測定点から標的位置への相互相関の相互相関ベクトルの計算を含む。この自己相関行列は、全ての標的位置に対して再使用することができる。自己および相互相関は、場表現において使用される平面伝搬波関数および平面エバネッセント波関数の領域内で定義される。したがって、第２の位置は、標的位置と測定位置の中から選択することができる。具体的には、この計算は、それぞれの測定点間の一組の平面波の相関関数の自己および相互相関行列の計算、ならびにそれぞれの測定点と標的点との間の一組の平面波の相関関数の相互相関ベクトルの計算を含むことができる。

それぞれの寄与を補間する補間関数という用語は、それらの寄与を局所的に近似する任意の関数または任意の一組の関数を指すことが意図されている。補間関数の表現が、それぞれの関数の予め計算された局所近似をそれぞれが含む一組のルックアップ・テーブルを含むときには、特に効率的な計算が提供される。それぞれのルックアップ・テーブルは、２つの入力パラメータによって索引が付けられた関数値の２次元配列を指示することができ、それぞれの関数値は、それらの２つの入力値によって定義される入力点における２次元配列の値を指示することができる。これらのそれぞれの関数はそれぞれ、２つの入力パラメータによってパラメータ表示された積分を含む２次元積分関数とすることができる。したがって、テーブル化された関数値間の補間、例えば３次補間（ｃｕｂｉｃ ｉｎｔｅｒｐｏｌａｔｉｏｎ）によって、別の点における関数値を効率的に得ることができる。このルックアップ・テーブルは低次元のテーブル、すなわち２以下の次元を有するテーブルであるため、この補間法は正確でかつ計算効率がよく、このルックアップ・テーブルは記憶領域をわずかしか必要とせず、少ない計算資源を用いて作成することができる。

２つの入力パラメータのうちの第１のパラメータは、波数を指示する係数（ｆａｃｔｏｒ）によってスケーリングされた、第１の測定位置と仮想源平面に平行な平面内の第２の位置の間の射影された距離を指示することができる。２つの入力パラメータのうちの第２のパラメータは、波数を指示する係数によってスケーリングされた、仮想源平面から第１の測定位置までの距離と仮想源平面から第２の位置までの距離の１次結合（一般に和または差）を指示することができる。ある実施形態では、このパラメータ表示が、相関関数を、予め計算された一組の２次元関数によって表現することを可能にする。第２の入力パラメータは、仮想源平面に対する法線（通常はｚ軸）上へ射影され、波数を指示する係数によってスケーリングされた、第１の測定位置と第２の位置の間の射影された距離を指示することができる。

ある実施形態では、それぞれの平面波波動関数が、所定の重み付けスキーム（ｗｅｉｇｈｔｉｎｇ ｓｃｈｅｍｅ）に基づくそれぞれの重み係数（ｗｅｉｇｈｔｉｎｇ ｆａｃｔｏｒ）によってスケーリング／重み付けされる。例えば、この所定の重み付けスキームは、一様重み付け（ｕｎｉｆｏｒｍ ｗｅｉｇｈｔｉｎｇ）と方向性重み付け（ｄｉｒｅｃｔｉｏｎａｌ ｗｅｉｔｈｔｉｎｇ）の中から選択することができる。一様重み付けスキームでは、それぞれの平面波が同じ定数係数によってスケーリングされる。平面波の方向性重み付けでは、それぞれの平面波が、平面波の伝搬方向に依存したそれぞれの重み係数によって、例えば、源平面に平行なより大きな成分を有する伝搬方向を有する平面波により高い重み係数を提供することによってスケーリングされる。このような波は、仮想源平面の表面法線（ｓｕｒｆａｃｅ ｎｏｒｍａｌ）に対して平行でないため、離軸波（ｏｆｆ−ａｘｉｓ ｗａｖｅ）とも呼ばれる。方向性重み付けは、再構成法の正確さを向上させることが分かっている。

ある実施形態によれば、所与の重み付けスキームに関して、さまざまな測定ジオメトリおよび再構成ジオメトリに対する音圧と粒子速度の両方を、６つの２次元関数の一組のテーブル表現に基づいて計算することができることが分かった。それによって補間された関数値を計算しなければならない関数パラメータだけがジオメトリとともに変化する。

ある実施形態では、仮想源平面が、測定領域とも呼ぶ、測定位置および標的位置を含む第１の空間領域と、源領域とも呼ぶ、第１の空間領域とは異なる（第１の空間領域とは反対の）第２の空間領域とを定義するように選択され、第２の空間領域は一般に音源の大部分を含む。すなわち、標的位置および測定位置は測定領域に位置し、音場は、測定領域において再構成される。仮想源平面は一般に、源平面の内側に短い距離を置いて定義されるが、この再構成法の目的上、仮想源平面によって定義される測定領域は源を含まないと仮定する。

状況によっては、単一の仮想源平面によって源のない空間領域を定義することができない、またはそうすることが実際的でないことがある。したがって、ある実施形態では、この方法が、それぞれの平面波がそれぞれの仮想源平面に対してある方向に伝搬し、減衰する２組以上の平面波に基づく。例えば、それぞれの仮想源平面対間に源のない領域が位置する間隔を置いて配置された１つまたは複数の平行仮想源平面対として、仮想源平面を定義することができる。本発明の発明者は、本発明の方法およびテーブル化された同じ一組の補間関数を、この状況にも適用することができ、それによって効率的で柔軟な再構成法を提供することができることを理解した。

この一組の測定位置を、１つまたは複数の測定平面内、例えば単一の平面または平行な２つ以上の平面内に配置することができる。それぞれの平面において、測定位置は、規則格子状にまたは不規則パターンとして、あるいは他の適当な方式で配置することができる。さらに、本明細書に記載した方法を、非平面測定ジオメトリ、すなわち測定位置が、１つまたは複数の平行平面内に配置されず、例えば曲面上に配置される配置に対して適用することもできる。

ある実施形態では、全ての仮想源平面に平行なそれぞれの平面内に測定点および計算点を配置すると、効率をさらに増大させることができる。このような構成では、平面波のそれぞれの周波数に対して、相関関数に対する寄与の入力パラメータのうちの１つの入力パラメータの値が少数だけ生じ、したがって、それぞれのテーブルが前記入力パラメータのそれぞれの値に対応する１Ｄルックアップ・テーブルの作成を作成することが魅力的になる。したがって、ある種のジオメトリに関しては、この方法を１次元補間関数に基づく方法とすることさえでき、したがってこの方法を特に高い計算効率で実行することができる。

第２の音響量が、仮想源平面に平行な方向の成分を有する粒子速度であるときには、曲面の音源面を有する物体の改良された解析が容易になる。具体的には、この実施形態は、（例えば音源面と共形の）曲面の再構成表面を横切って分布する標的位置における改良された音場の再構成を可能にする。再構成法のこの実施形態は、仮想源平面に平行な粒子速度の成分を提供するため、再構成表面に垂直な粒子速度（および／または強さ）の成分を、曲面の再構成表面を横切って決定することができる。

本明細書に開示した方法の実施形態を、標的点における音場パラメータの線形予測を一組の測定量に基づいて実行するものとみなすことができ、予測係数は、最小２乗法の意味において、全ての平面伝搬波および重み付けされた一組のエバネッセント波に対する最良推定値を提供するように決定される。最小２乗フィット（ｌｅａｓｔ−ｓｑｕａｒｅｓ ｆｉｔ）は、優決定（ｏｖｅｒ−ｄｅｔｅｒｍｉｎｅｄ）の一組の１次方程式を近似的に解く一組の複素係数を決定する数値的プロセスとみなすことができる。

次いで、計算点において、これらの平面波波動関数の値が、平均して、どれくらいよく予測されるのかを知ることが望ましい。したがって、ある実施形態では、再構成法が、計算された音響量の推定誤差を指示する量を計算することを含む。以前に説明した自己および相互相関を使用して、すなわち計算資源をほとんど追加せずに、誤差インジケータ（ｅｒｒｏｒ ｉｎｄｉｃａｔｏｒ）を計算することができる。ＳＯＮＡＨ法の予測の正確さを推定する上記の方法は、所与の測定ジオメトリを有するＳＯＮＡＨ予測の妥当性領域を視覚化し、決定するのに非常に有用であることが立証されている。

上述の方法および後述する方法の諸特徴は、ソフトウェアまたはファームウェア内において少なくとも部分的に実現することができ、コンピュータ実行可能命令などのプログラム・コード手段を実行することによって、データ処理装置または他の処理手段上で実施することができることに留意されたい。この説明および以下の説明では、用語「処理手段」が、上記の機能を実行するように適切に適合された任意の回路および／または装置を含む。具体的には、用語「処理手段」は、汎用または専用プログラマブル・マイクロプロセッサ、ディジタル信号プロセッサ（ＤＳＰ）、特定用途向けＩＣ（ＡＳＩＣ）、プログラマブル論理アレイ（ＰＬＡ）、フィールド・プログラマブル・ゲート・アレイ（ＦＰＧＡ）、専用電子回路など、あるいはこれらの組合せを含む。

本発明の実施形態は、上述の方法および後述する方法、ならびに最初に述べた方法に関して説明した利益および利点のうちの１つまたは複数の利益および利点をそれぞれがもたらし、最初に述べた方法に関して説明した実施形態および／または従属請求項に開示された実施形態に対応する１つまたは複数の実施形態をそれぞれが有するシステム、装置および製品手段を含む、さまざまな方法で実現することができる。

具体的には、音場を再構成する処理装置の実施形態は、一組の測定位置のそれぞれで測定した少なくとも１つの第１の音響量の一組の測定値を受け取るインタフェースと、
− 前記測定位置のうちの第１の測定位置における少なくとも一組の平面波と第２の位置における少なくとも一組の平面波との相関関係をそれぞれが指示する一組の相関関数のそれぞれを計算し、
− 少なくとも計算した相関関数と一組の測定値とから、標的位置における再構成された音場を指示する第２の音響量を計算する
ように構成された処理ユニットとを備え、
処理ユニットは、相関関数に対するそれぞれの寄与を補間する一組の補間関数の表現を記憶する記憶媒体を備え、それぞれの寄与は、入力パラメータが２つ以下の関数であり、処理ユニットは、一組のそれぞれの相関関数を、記憶した表現から得た値の１次結合として計算するように適合されている。

記憶された表現は、再構成法を実行するコンピュータ・プログラムのソフトウェア・コンパイル時間に（例えばコンピュータ・プログラムの一部としてまたは別個の１つもしくは複数のファイルとして）生成することができ、あるいは、コンパイルの結果生成されたプログラム・コードによって実行時に生成することができ、あるいは、これらを組み合わせて生成することができることが理解される。例えば、この表現は、コンピュータ・プログラムまたはインスレーション・プログラムによって生成し、そのコンピュータ・プログラムの最初の実行／インストレーションの間にコンピュータ上に記憶することができる。

音場を再構成するシステムは、上に開示し、後に開示する装置と、一組の測定位置で第１の音響量を測定する一組の変換器であり、測定した値を装置へ送るように、装置に通信接続で接続可能な一組の変換器とを備える。

コンピュータ・プログラムは、プログラム・コード手段を含むことができ、このプログラム・コード手段は、このプログラム・コード手段がデータ処理システム上で実行されたときに、データ処理システムが、上に開示し、後に開示する方法のステップを実行するように適合される。このコンピュータ・プログラムは、記憶手段上に記憶することができ、またはデータ信号として具体化することができる。記憶手段は、ＲＡＭ、ＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ、フラッシュ・メモリ、磁気または光記憶装置（ＣＤ ＲＯＭ、ＤＶＤ、ハード・ディスクおよび／または他の同種の装置など）など、データを記憶する適当な回路または装置を備えることができる。

上記の態様およびその他の態様は、図面を参照して以下で説明する実施形態によって明白となる。

音場を再構成するシステムの概略ブロック図である。 再構成された音場を計算する方法の流れ図である。 仮想源平面が１つだけの状況、２つの平行な仮想源平面が源のない均質な領域を限定している状況、および箱形の均質な源のない領域の６つ全ての面に仮想源平面がある状況の測定ジオメトリを示す図である。 （ｋ_ｘ，ｋ_ｙ）領域および伝搬角領域における平面波波動関数の密度を示す図である。 アレイ平面内の音圧予測の固有相対誤差水準Ｅ（ｒ）を、２つの不規則平面アレイについて示す図である。 単層および２重層８×８要素アレイのｘｚ対称平面における相対誤差水準を示す図である。左列：単層アレイの固有誤差水準Ｅ（ｒ）。中央列：単層アレイによって推定される単極子圧力の相対誤差。右列：２重層アレイの固有誤差水準Ｅ（ｒ）。 （ａ）８×８要素平面アレイのアレイ平面に垂直な対称軸上の固有相対誤差水準Ｅを示す図である。 ８×８要素平面アレイのアレイ平面に垂直な対称軸上の振幅利得Ａを示す図である。 ８×８要素平面アレイのアレイ平面に垂直な対称軸上の固有相対誤差水準Ｅ（ｒ）および振幅利得Ａ（ｒ）を、低周波数および高周波数について示す図である。これらの曲線は、再構成法で使用する異なるダイナミック・レンジＤを表す。 ８×８要素平面アレイのアレイ平面に垂直な対称軸上の固有相対誤差水準Ｅ（ｒ）および振幅利得Ａ（ｒ）を、低周波数および高周波数について示す図である。これらの曲線は、再構成法で使用する異なる仮想源平面距離ｄ＝−ｚ^＋を表し、エバネッセント波の実際の内容を定義する仮想源距離は一定の６ｃｍに維持される。 ２重層８×８アレイの中心平面内の単一の単極子点源から３Ｄインテンシティベクトルを計算する一実施形態において、２つの面だけからではなく６つの全ての面から平面波波動関数を追加することによって得られる正確さの向上を示す図である。点源はアレイ・エリアから６ｃｍ外側にある。

全図面を通じ、可能な場合には常に、等しい参照符号は、等しい要素、特徴または成分、あるいは対応する要素、特徴または成分を指す。

図１は、音場を再構成するシステムの概略ブロック図を示す。このシステムは、一組の音響受信器（ａｏｕｓｔｉｃ ｒｅｃｅｉｖｅｒ）１０８と、音響受信器に接続された解析ユニット（ａｎａｌｙｓｉｎｇ ｕｎｉｔ）１０３とを備える。

以下では、音響受信器１０８を変換器（ｔｒａｎｓｄｕｃｅｒ）とも呼ぶ。のみならず、音響受信器は、マイクロホン、水中マイクロホン、あるいは音圧、音圧勾配、粒子速度または他の線形量などの音響特性を測定する他の適当な装置とすることができることが理解される。図１の実施形態では、変換器１０８が、変換器１０８が規則格子、例えば１、２または３次元格子として配置された変換器のアレイ１０２として実現される。規則格子として配置すると効率的な数値実現が可能だが、他のジオメトリを使用することもできる。ある実施形態では、単一の変換器を使用して、異なる点における測定を逐次的に実施することもできる。このような逐次測定は例えば、音源が静止している状況および１つまたは複数の基準信号に対する位相を測定することができる状況において使用することができる。変換器の数およびアレイのジオメトリ、例えば変換器間の間隔は、解析する物体のサイズおよび幾何学的複雑さ、関心の周波数範囲ならびに所望の空間分解能に応じて選択することができる。

変換器アレイ１０２は、変換器１０８が、測定した信号を、例えば有線または無線信号接続を介して解析ユニットへ送ることができるような態様で、解析ユニット１０３に接続される。

解析ユニット１０３は、測定された信号を変換器アレイ１０２から受け取り、これを処理するインタフェース回路１０４と、インタフェース回路１０４とデータ通信する処理ユニット１０５と、記憶媒体１１２と、処理ユニット１０５とデータ通信する出力ユニット１０６とを含む。図１には単一のユニットとして示されているが、解析ユニット１０３を、２つの別個の装置、例えば収集フロントエンド装置とコンピュータとに物理的に分割し、または３つ以上の装置に分割することもできることが理解される。同様に、解析ユニットの異なる下位ブロックに関して記載した機能を、代替のまたは追加の機能ユニット／モジュールまたはハードウェア・ユニット／モジュールに分割することができることも理解される。

インタフェース回路は、変換器１０８からの出力信号を受け取り、受け取った信号を、処理ユニット１０５による後続の解析のために処理するのに適した信号処理回路を含む。インタフェース回路１０４は、以下の構成要素のうちの１つまたは複数の構成要素を備えることができる：受け取った信号を増幅する１つまたは複数の前置増幅器、受け取った信号をディジタル信号に変換する１つまたは複数のアナログ−ディジタル（Ａ／Ｄ）変換器、１つまたは複数のフィルタ、例えば帯域幅フィルタ、および／あるいは他の同種の構成要素。インタフェース回路は例えば、出力データとして、変換器ごとの振幅および位相を、周波数の関数として提供することができる。ある実施形態では、このインタフェース・ユニットが同時データ収集を実行し、次いで、処理ユニット１０５が、周波数領域へのデータ変換を含む残りの全ての処理を、一般にＦＦＴを使用して実行することができる。

処理ユニット１０５は、適当にプログラムされたマイクロプロセッサ、コンピュータの中央処理ユニット、あるいはインタフェース・ユニット１０４から受け取った信号を処理する他の適当な装置、例えばＡＳＩＣ、ＤＳＰおよび／または他の同種の装置とすることができる。この処理ユニットは、インタフェース回路１０４を介して受け取った変換器信号を処理して、再構成された音場を本明細書に記載されたとおりに計算するように適合される。

記憶媒体１１２は、ＲＡＭ、ＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ、フラッシュ・メモリ、磁気または光記憶装置（ＣＤ ＲＯＭ、ＤＶＤ、ハード・ディスクおよび／または他の同種の装置など）など、一組の補間関数の表現を指示するデータを記憶する適当な回路または装置を備えることができる。図１には、処理ユニットとは別個だが、処理ユニットに通信接続された記憶媒体が示されている。しかしながら、記憶媒体１１２は、処理ユニット１０５の一部として、例えば内部記憶装置として具体化することもできる。

出力ユニット１０６は、表示装置、あるいは再構成された音場の視覚表現を提供する他の適当な装置または回路、例えば印刷された表現を提供するプリンタおよび／またはプリンタ・インタフェースを備えることができる。あるいは、またはこれに加えて、出力ユニット１０６は、ＲＡＭ、ＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ、フラッシュ・メモリ、磁気または光記憶装置（ＣＤ ＲＯＭ、ＤＶＤ、ハード・ディスクなど）、有線または無線データ通信インタフェース、例えばコンピュータまたは電気通信ネットワーク（ＬＡＮ、ワイド・エリア・ネットワーク、インターネットなど）とのインタフェース、および／あるいは他の同種の回路または装置など、再構成された音場を指示するデータを伝達し、かつ／または記憶する適当な回路または装置を備えることができる。

解析ユニット１０３は、適当にプログラムされたコンピュータ、例えば適当な信号収集ボードまたは回路を含むＰＣとして実現することができる。

このシステムはさらに、解析ユニット１０３に接続された適当な位置検出装置１１０であって、変換器１０８の位置、例えば基準座標系に対する変換器１０８の位置を検出する位置検出装置１１０を備えることができる。この位置検出装置は、変換器アレイを所定の位置および向きに取り付けることができ、かつ／あるいはアレイの位置および向きを手動で測定し、または例えば適当なセンサによって自動的に測定することができる、試験台、フレームまたは他の構造を含むことができる。あるいは、変換器アレイ１０２が、位置検出装置または位置検出装置の少なくとも一部分を含んでもよい。例えば、変換器アレイが、磁場位置検出システムのコイル、ジャイロスコープ、光学的位置検出システムおよび／または他の同種のシステムを含むことができる。

操作時には、音響放射を放出する音源１１１であり、音場を再構成する対象である音源１１１を含む物体１００の表面の近くに、変換器アレイ１０２を配置する。変換器の数、アレイのジオメトリ、例えば変換器間の間隔、およびアレイと音源面との間の距離は、解析する物体のサイズおよび幾何学的複雑さ、関心の周波数範囲ならびに所望の空間分解能に応じて選択することができる。例えば、自動車エンジンの音場の再構成に対しては、要素間間隔が３ｃｍの８×８行列のマイクロホン・アレイを使用することができる。しかしながら、他のタイプのアレイおよび他のサイズのアレイを使用することもできることが理解される。

例えば位置検出装置１１０によってアレイ１０２の位置が決定され、解析ユニット１０３に供給される。アレイ１０２の変換器１０８が、決められた位置において音圧または他の適当な音響量を測定し、その結果得られた変換器信号が解析ユニット１０３へ送られる。

例えば、変換器アレイを、位置検出装置が組み込まれた手持ち式アレイとすることができ、こうすると、物体の周囲に分布した異なるアクセス可能位置で測定することができ、試験に対する源の特別な準備にかかる時間が短縮される。例えば、試験台上の音源、例えば圧縮機、ポンプまたは他の機械をマップするときには、面全体に障害なくアクセスできるようにするために試験構成を変更するのに必要な時間、例えば燃料管路の配置、電子回路の配線などを変更するのに必要な時間が、排除されることはないにしても、短縮される。他の典型的な用途は、自動車の室内での使用であり、そこでは、３Ｄアレイ格子を使用して、全ての方向の源を識別することができる。例えば（例えば８×８×２センサを含む）２層アレイを使用することができる。

解析ユニット１０３は、この入力データから、再構成された音場を計算し、再構成された音場の表現を記憶し、かつ／または出力する。

次に、図２を参照し、図１を続けて参照して、再構成された音場を生成する方法の一実施形態を説明する。

本明細書に記載した再構成法を例示するため、最初に、いわゆるＳＯＮＡＨ法の概要を示す。

この説明の目的上、複素時間調和音圧（ｃｏｍｐｌｅｘ ｔｉｍｅ−ｈａｒｍｏｎｉｃ ｓｏｕｎｄ ｐｒｅｓｓｕｒｅ）ｐ（ｒ_ｉ）は、均質流体によって占められた源のない領域Ω内の一組の位置ｒ_ｉ、ｉ＝１，２，．．．，Ｉにおいて測定されたものであると仮定する。測定した音圧に基づいてΩ内における音場を再構成する問題について考える。これを達成するため、最初に、同次波動方程式（ｈｏｍｏｇｅｎｅｏｕｓ ｗａｖｅ ｅｑｕａｔｉｏｎ）を満たし、Ω内に存在しうる全ての音圧場を十分に高い正確さ（ａｃｃｕｒａｃｙ）で表現することができる、離散した一組の基本波動関数（ｅｌｅｍｅｎｔａｒｙ ｗａｖｅ ｆｕｎｃｉｔｏｎ）Ψ_ｎ、ｎ＝１，．．．，Ｎを導入する。一般に、これらの基本波動関数は、平面波、円筒波または球面波とすることができ、あるいは、Ωの外側に分布した単極子点源または双極子点源などの基本源からの音場とすることができる。

いわゆるＨＥＬＳ法では、再構成の最初のステップが、音場の表現中のそれぞれの基本波動関数の振幅および位相を決定するステップである。このステップは、基本波動モデルが、測定された圧力をできるだけ正確に表現するようにすることによって実行される。

ａ_ｎは複素展開係数（ｃｏｍｐｌｅｘ ｅｘｐａｎｓｉｏｎ ｃｏｅｆｆｉｃｉｅｎｔ）である。Ｉ≧Ｎであり、測定点が適切に選択された場合には、（１）の最小２乗解によって、展開係数を見つけることができる。そうでない場合、この問題は劣決定であり、（１）の無限の解の中で、最小の２ノルムを有する解を採用することができる。この解を、ヘルムホルツ方程式最小ノルム（ＨＥＬＮ：Ｈｅｌｍｈｏｌｔｚ’Ｅｑｕａｔｉｏｎ Ｌｅａｓｔ Ｎｏｒｍ）解と呼ぶ。

行列公式（ｍａｔｒｉｘ ｆｏｒｍｕｌａｔｉｏｎ）を得るため、測定位置における波動関数値の行列Ｂを定義し、

測定された圧力および展開係数を用いてベクトルｐおよびａをそれぞれ定義する。

これらの定義を使用して、式（１）を次の形に書き直すことができる。
Ｂａ＝ｐ （４）
Ｉ＞Ｎにおいて、最小２乗（ＨＥＬＳ）解は、
ａ＝（Ｂ^ＨＢ）^−１Ｂ^Ｈｐ
であり、Ｉ＜Ｎにおいて、最小ノルム（ＨＥＬＮ）解は、
ａ＝Ｂ^Ｈ（ＢＢ^Ｈ）^−１ｐ
であり、記号Ｈは、エルミート転置（Ｈｅｒｍｉｔｉａｎ ｔｒａｎｓｐｏｓｅ）を表す。実際には、これらの解に、正則化を適用しなければならない。他の正則化スキームを使用することもできるが、この説明の目的上、チホノフ正則化だけを考え、その場合、ＨＥＬＳおよびＨＥＬＮ式はそれぞれ、
ａ＝（Ｂ^ＨＢ＋εＩ）^−１Ｂ^Ｈｐ （５）
および
ａ＝Ｂ^Ｈ（ＢＢ^Ｈ＋εＩ）^−１ｐ （６）
となる。εは正の正則化パラメータ、Ｉは適当な次元の単位対角２次行列である。実際に、正則化されたこれらの２つの式（５）と（６）は、ＩおよびＮの値とは無関係に全く同じ解を与えることが分かる。しかしながら、それでも、Ｉ＜Ｎに関しては、行列ＢＢ^Ｈの次元の方が行列Ｂ^ＨＢの次元よりも小さいため、ＨＥＬＮ式（６）を使用した方が計算上有利である。

展開係数が分かれば、波動関数の寄与の和を求めることによって、Ω内の任意の点ｒにおける音圧を計算することができる。

Ｔは行列またはベクトルの転置を表す。次いで、オイラーの方程式を適用することによって、音場の表現（７）から、χ方向の粒子速度を推定することができる。

上式で、ωは角周波数、ρ_０は媒質の密度であり、陰時間変動（ｉｍｐｌｉｃｉｔ ｔｉｍｅ ｖａｒｉａｔｉｏｎ）は、ｅ^ｊωｔである。続いて使用する他の記号は、媒質中における音の伝搬速度ｃおよび波数ｋ≡ω／ｃである。

平面波は、波数領域において連続スペクトルを構成し、そのため、上記の公式におけるこれらの使用は、波数領域におけるサンプリングによって得られた波動関数の部分集合を含まなければならない。本明細書に記載した方法の実施形態は、例えばＪ．Ｈａｌｄ、「Ｐａｔｃｈ ｎｅａｒ−ｆｉｅｌｄ ａｃｏｕｓｔｉｃａｌ ｈｏｌｏｇｒａｐｈｙ ｕｓｉｎｇ ａ ｎｅｗ ｓｔａｔｉｓｔｉｃａｌｌｙ ｏｐｔｉｍａｌ ｍｅｔｈｏｄ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｉｎｔｅｒ−Ｎｏｉｓｅ ２００３およびＪ．Ｈａｌｄ、「Ｐａｔｃｈ ｈｏｌｏｇｒａｐｈｙ ｉｎ ｃａｂｉｎ ｅｎｖｉｒｏｎｍｅｎｔｓ ｕｓｉｎｇ ａ ｔｗｏ−ｌａｙｅｒ ｈａｎｄｈｅｌｄ ａｒｒａｙ ａｎｄ ａｎ ｅｘｔｅｎｄｅｄ ＳＯＮＡＨ ａｌｇｏｒｉｔｈｍ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｅｕｒｏｎｏｉｓｅ ２００６に記載されているＳＯＮＡＨ法を利用した２重の方法で場推定問題を解くことによって、このサンプリングの問題を解決する。

ＳＯＮＡＨ法によれば、Ω内の点ｒにおける未知の音圧ｐ（ｒ）は、既知の音圧値の１次結合として得られる。
ｐ（ｒ）＝ｐ^Ｔｃ（ｒ） （９）
上式で、ベクトルｃ（ｒ）は複素推定重み（ｃｏｍｐｌｅｘ ｅｓｔｉｍａｔｉｏｎ ｗｅｉｇｈｔ）を含む。次いで、予測式（９）が、最小２乗法の意味において、全ての展開関数Ψ_ｎに対する最良平均予測を提供するようにすることによって、ベクトルｃを得る。式（２）および（７）を参照すると、これは、以下の一組の１次方程式に対する最小２乗解としてｃ（ｒ）を見つけなければならないことを意味する。
Ｂ^Ｔｃ（ｒ）＝α（ｒ） （１０）
この式は、ＨＥＬＳ式（４）が優決定であるときは常に劣決定であり、ＨＥＬＳ式（４）が劣決定であるときは常に優決定である。この説明の目的上、チホノフ正則化を適用した解を考える。このことは、（１０）に対する最小２乗解と最小ノルム解の両方を、以下の最小２乗式に関して表現することができることを意味する。
ｃ（ｒ）＝（Ａ^ＨＡ＋εＩ）^−１Ａ^Ｈα（ｒ） （１１）
上式で、行列Ａは行列Ｂの転置である。
Ａ≡Ｂ^Ｔ （１２）
定義（７）および（１２）から、続いて、要素Ａ^ＨＡおよびＡ^Ｈα（ｒ）が、それぞれ下式（１３）および（１４）によって与えられる。「^＊」は複素共役化（ｃｏｍｐｌｅｘ ｃｏｎｊｕｇａｔｉｏｎ）を表す。

行列Ａ^ＨＡは、基本波動関数（展開関数）の領域内の測定点間の相互相関の行列と見ることができ、ベクトルＡ^Ｈα（ｒ）は、測定点と推定位置ｒとの間の相互相関を含む（例えば、Ｒ．ＳｔｅｉｎｅｒおよびＪ．Ｈａｌｄ、「Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ ｗｉｔｈｏｕｔ ｔｈｅ ｅｒｒｏｒｓ ａｎｄ ｌｉｍｉｔａｔｉｏｎｓ ｃａｕｓｅｄ ｂｙ ｔｈｅ ｕｓｅ ｏｆ ｓｐａｔｉａｌ ＤＦＴ」、Ｉｎｔｅｒｎ．Ｊ．Ａｃｏｕｓｔ．Ｖｉｂ．６、８３〜８９（２００１）を参照されたい）。さらに、Ａ^ＨＡが、負でない固有値を有するエルミート対称行列であり、その次元が、展開関数の数とは無関係に、測定点の数Ｉに等しいことに気づくこともある。

上記の予測公式をＳＯＮＡＨ予測公式と呼び、この公式はＨＥＬＳ／ＨＥＬＮ式と同じ解を提供する。実際、ＨＥＬＮ式（７）および（６）から、このＳＯＮＡＨ式は導かれる。

ここでは、任意の非特異２次行列（ｎｏｎ−ｓｉｎｇｕｌａｒ ｑｕａｄｒａｔｉｃ ｍａｔｒｉｘ）について、逆（ｉｎｖｅｒｓｅ）の転置は転置の逆に等しいことを利用した。ＳＯＮＡＨは展開係数ａを明示的には計算しないため、一組の展開係数に対して重み付け／フィルタリングを実行することはできない。このようなフィルタリングは、重みベクトルｃ（ｒ）の計算に組み込まれている。これを実行するツールは正則化だが、これを適切に機能させるためには、行列Ａ^ＨＡを定義する基本波動関数の適当なスケーリングが必要である。一般に、より大きな正則化パラメータε値を選択することによって、音源面における再構成において最大の増幅を必要とするまだ抑制されていない基本エバネッセント波を抑制することが望ましい。このことが起こるためには、これらの基本波が、チホノフ正則化のカットオフ・レベルよりもまだ大きいＡ^ＨＡの最も小さな特異値（ｓｉｎｇｕｌａｒ ｖａｌｕｅ）に対応していなければならない。基本波動関数のスケーリングを通して、再構成においてより大きな増幅を必要とする波動関数がマイクロホン位置を横切ってより小さな振幅を有することを保証することによって、このような振舞を支持することができる。

この実施形態では、基本波が、平面伝搬波およびエバネッセント波であり、最初に、図３ａに示すようなフリー・フィールド条件（ｆｒｅｅ−ｆｉｅｌｄ ｃｏｎｄｉｔｉｏｎ）、すなわち、ｚ＝ｚ_０、ｚ_０＜０の源平面１０７の後ろの１つまたは複数の源が、Ωで示す源のない均質な半空間ｚ≧ｚ_０に放射する条件において、音場を再構成する場合を考える。この場合には、

の形の平面波波動関数によって放射音場を表現することができる。上式で、Ｆは振幅スケーリング関数、ｚ＝ｚ^＋、ｚ^＋≦Ｚ_０は、平面波波動関数を定義するいわゆる仮想源平面１０１であり、例えばエバネッセント成分はこの仮想源平面に垂直な方向に減衰する。仮想源平面は一般に、実際の音源面から内側に短い距離を置いて位置するように定義される。源平面１０７は、依然として源のない領域にある仮想源平面に平行な最も近い平面である。したがって源平面は、実際の音源面にちょうど触れる（すなわち実際の音源面に接する）平行な平面である。単一の仮想源平面だけがあるフリー・フィールド条件の場合には、ｘｙ平面が一般に仮想源平面に平行である。式（１６）では、ｒ≡（ｘ，ｙ，ｚ）およびｋ≡（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）であり、ｋ_ｘおよびｋ_ｙは実波数、ｋ_ｚは、下式で定義される複素数である。

したがって、完全な一組の展開関数は、２つの波数変数ｋ_ｘおよびｋ_ｙにおいて連続スペクトルを構成する。仮想源平面ｚ＝ｚ^＋は実源平面ｚ＝ｚ_０と一致することがあるが、仮想源平面は一般に実源平面の後ろにある。すなわちｚ^＋＜ｚ_０である。仮想源平面の座標上の上付き添字「＋」は、ｚ軸の正方向に伝搬している波が表現されていることを示す。

式（１６）では、スケーリング関数Ｆ（ｋ_ｚ）が、仮想源平面ｚ＝ｚ^＋上の全ての波動関数の振幅を定義することに留意されたい。

上述の離散ＨＥＬＳまたは離散ＳＯＮＡＨ公式を使用するため、最初に、２Ｄ波数領域（ｋ_ｘ，ｋ_ｙ）における規則的サンプリング（ｒｅｇｕｌａｒ ｓａｍｐｌｉｎｇ）によって離散部分集合を選択する。

サンプリング間隔を選択するときには、波数領域における規則的な離散表現が、ちょうど空間ＤＦＴ処理に基づく伝統的なＮＡＨの場合のようにラップアラウンド誤差（ｗｒａｐ−ａｒｏｕｎｄ ｅｒｒｏｒ）につながる空間的に周期的な場を表現することを知っておくべきである。サンプリング間隔はしたがって、２πを最小許容周期長（ｓｍａｌｌｅｓｔ ａｃｃｅｐｔａｂｌｅ ｐｅｒｉｏｄ ｌｅｎｇｔｈ）で割ったものよりも小さくすべきである。同時に、サンプリング・エリアは、全ての伝搬波、および測定領域においてかなりの振幅を有する全てのエバネッセント波をカバーすべきである。これによって一般に、展開関数の数Ｎが測定位置の数Ｉよりも大きくなる。したがって、等価（ｅｑｕｉｖａｌｅｎｔ）のＨＥＬＳ問題は劣決定となり、解（６）は支配的に最小ノルム解となる。

ある実施形態では、Ｆ（ｋ_ｚ）＝１の場合、波動関数の単位重み付け（ｕｎｉｔ ｗｅｉｇｈｔｉｎｇ）を使用することができる。この場合、式（１６）によれば、全ての平面波波動関数が、仮想源平面ｚ＝ｚ^＋における振幅に等しい振幅を有する。したがって、ＨＥＬＳ展開係数ベクトルａは単純に、その平面における音場表現の波数スペクトルであり、行列Ｂを、その平面波スペクトルからマイクロホン位置における音圧値への伝達関数の行列とみなすことができる。Ｂの要素は、式（１６）にマイクロホン座標を挿入することよって与えられる。伝達行列Ｂによって、エバネッセント波成分（Φ_ｋ、｜ｋ｜＞ｋ）の減衰が、波数とともに指数関数的に増大することに言及することは興味深い。したがって、全体的な傾向は、これらの波数成分は、行列Ｂのますます小さな特異値に対応するというものである。この特性は、波数スペクトルａに対する逆解（ｉｎｖｅｒｓｅ ｓｏｌｕｔｉｏｎ）におけるチホノフ正則化によって、エバネッセント波成分のフィルタリングを実行することができることを保証する。伝搬波のある組合せ、例えば仮想源平面内の測定位置から遠いエリアからの音波放射を表現する組合せも小さな特異値に対応する。このような組合せも、正則化によって抑制される。

ＨＥＬＳ問題の上述の劣決定性のため、解は、支配的に最小ノルム解となる。正則化によって、解ベクトルａのノルムのさらなる低減が導入される。したがって、この解法は、波数スペクトルにおけるエネルギーを最小化し、結果的に、マイクロホンにおける圧力がうまく表現されているという制約条件下で、仮想源平面内および平行な平面内の圧力分布におけるエネルギーも最小化する。これにより、測定位置の周囲の体積内において良好な音場表現が達成され、測定位置から離れるにつれて増大する過小推定（ｕｎｄｅｒ−ｅｓｔｉｍａｔｉｏｎ）が見られる。ＳＯＮＡＨとＨＥＬＳは同じ解（１５）を提供するため、後に説明する結果によってさらに示すように、ＳＯＮＡＨに対する結果も同じである。

図４は、ｋ_ｙ＝０におけるｋ_ｘ波数軸に沿った式（１８）の規則的サンプリングを示す。放射円内において、すなわち
に関して、（ｋ_ｘ，ｋ_ｙ）における特定の標本は、ベクトル（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）の方向に伝搬している平面波に対応する。ｋ_ｚは式（１７）によって与えられる。図４では、示された半径ｋの半円上に標本位置ｋ_ｘを射影することによって、この方向を見つけることができる。（ｋ_ｘ，ｋ_ｙ）平面における標本点の密度が一定のとき、平面波の角密度はｚ方向の近くで最も高く、軸角度（ｏｆｆ−ａｘｉｓ ａｎｇｌｅ）θが９０°に向かって増大すると低下する。したがって、一定の重み関数Ｆを使用すると、軸方向の近くで伝搬している平面波は、ＳＯＮＡＨの予測重みに対する最小２乗解（１１）において、軸から遠く外れた波よりも高い重みを受け取る。これは、離軸波により高い重みＦを与えることによって補償することができる。（ｋ_ｘ，ｋ_ｙ）平面上の無限に小さな面積の弓形（ｓｅｇｍｅｎｔ）を考えると、それは、１／ｃｏｓ（θ）＝ｋ／｜ｋ_ｚ｜倍だけ大きな方向半球（ｄｉｒｅｃｔｉｏｎａｌ ｈｅｍｉｓｐｈｅｒｅ）
上のエリアを表す。したがって、
｜Ｆ（ｋ_ｚ）｜^２ｋ／｜ｋ_ｚ｜
が一定になるように、例えば
となるようにＦを選択することによって、一定の方向パワー・スペクトル密度を得ることができる。上式で、Ｆ_０は定数である。式（１３）、（１４）、（１６）および（１８）は、｜Ｆ｜^２によってスケーリングされた波数寄与の全体にわたって和を求めることにより、行列Ａ^ＨＡおよびＡ^Ｈα（ｒ）の要素に対する寄与が得られることを示している。このことが、振幅密度に関してではなくパワー密度に関して等化（ｅｑｕａｌｉｚａｔｉｏｎ）を選択することが有利であることがある理由である。放射円の外側にも適用される式（１９）の「全方向重み付け（ｏｍｎｉｄｉｒｅｃｔｉｏｎａｌ ｗｅｉｇｈｔｉｎｇ）」は、平均して、一定の重み付けＦ（ｋ_ｚ）＝１よりもわずかによい正確さを提供することが分かっている。したがって、ある実施形態では、全方向重み付けが使用され、これはさらに、この説明の全体を通じてこの方法を示すために使用される。

波数領域における離散表現によって生じるラップアラウンド問題を完全に除去するため、次に、（ｋ_ｘ，ｋ_ｙ）平面におけるサンプリング間隔をゼロにする。このことは、式（１３）および（１４）における求和が（ｋ_ｘ，ｋ_ｙ）平面にわたって積分になることを意味する。式（１１）は、正則化係数がそれに応じて調整される場合、Ａ^ＨＡおよびＡ^Ｈα（ｒ）の要素に適用される共通のスケーリング係数が、ＳＯＮＡＨの結果を変化させないことを示している。したがって、極限積分式（ｌｉｍｉｔｉｎｇ ｉｎｔｅｇｒａｌ ｅｘｐｒｅｓｓｉｏｎ）を、以下のようにスケーリングすることを選択することができる。

行列要素の極限値（２０）および（２１）は、無限の次元を有する行列Ａおよびベクトルα（ｒ）に対応する。とにかくもＡおよびα（ｒ）は、後続のいくつかの微分において使用され、その場合、これらを、極限Ｎ→∞の場合と理解すべきである。

式（１６）、（１７）および（１９）を適用することによって、式（２０）および（２１）の２次元積分を１次元積分にすることができる。最初に、波数領域において極座標（Ｋ，Ψ）への変換を実行する（ｋ_ｘ，ｋ_ｙ）＝（Ｋｃｏｓ（Ψ），Ｋｓｉｎ（Ψ））。これにより、極角（ｐｏｌａｒ ａｎｇｌｅ）Ψにおける積分を解析的に解くことができる。その結果、式（２０）および（２１）の積分を、以下の２つの２変数関数の形で表現することができる。

上式で、Ｊ_０は、０次のベッセル関数（Ｂｅｓｓｅｌ ｆｕｎｃｔｉｏｎ）である。関数Ｉ_ｐは、積分式において置換Ｋ＝ｋｓｉｎ（θ）を実施した波数領域における放射円の内側からの寄与を表し、関数ＩＩ_ｐは、積分式において置換
を実施した放射円の外側からの寄与を表す。これらの変数置換は、（２２）および（２３）の中の積分が、ガウス型の数値積分により適合するように導入したものである。関数Ｉ_ｐおよびＩＩ_ｐを使用して、行列要素に対する式（２０）および（２１）を、以下のように再公式化することができる。

ここで、ｒ_ｉ≡（ｘ_ｉ，ｙ_ｉ，ｚ_ｉ）およびｒ_ｊ≡（ｘ_ｊ，ｙ_ｊ，ｚ_ｊ）は２つのマイクロホン位置であり、圧力推定位置はｒ≡（ｘ，ｙ，ｚ）であり、
は、ｘｙ平面上（仮想源平面）上におけるマイクロホン位置ｉとｊの射影間の距離、
は、ｘｙ平面上におけるマイクロホン位置ｉと圧力計算位置の射影間の距離である。こうすると、式（９）および（１１）を使用して、位置ｒにおける圧力を推定することができる。式（１１）の正則化パラメータεを決定するさまざまな方法が、Ｊ．ＧｏｍｅｓおよびＰ．Ｃ．Ｈａｎｓｅｎ、「Ａ ｓｔｕｄｙ ｏｎ ｒｅｇｕｌａｒｉｚａｔｉｏｎ ｐａｒａｍｅｔｅｒ ｃｈｏｉｃｅ ｉｎ Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ａｃｏｕｓｔｉｃｓ’０８、２００８において検討され、ＳＯＮＡＨに関しては、Ｌ曲線法（Ｌ−ｃｕｒｖｅ ｍｅｔｈｏｄ）が全体的に最も安定した結果を与えることが分かった。これよりも初期の非相関（ｎｃｏｒｒｅｌａｔｅｄ）測定誤差に関する検討では、一般化交差検証（ＧＣＶ：Ｇｅｎｅｒａｌｉｚｅｄ Ｃｒｏｓｓ Ｖａｌｉｄａｔｉｏｎ）を指向していた。

粒子速度の推定値を得るため、式（９）および（１１）によって与えられるＳＯＮＡＨ圧力推定値にオイラーの方程式（８）を適用する。その結果として、
を得、上式で
である。

Ａ^Ｈβ_χ（ｒ）に対する陽式（ｅｘｐｌｉｃｉｔ ｆｏｒｍｕｌａ）は、式（２７）に式（２５）を挿入し、続いて関数Ｉ_ｐおよびＩＩ_ｐに関する式（２２）および（２３）を適用することによって得ることができる。ｚ成分に関しては、
が得られ、ｘおよびｙ成分に関しては、
が得られる。新たな関数Ｉ_ｚ、ＩＩ_ｚ、Ｉ_ｘｙおよびＩＩ_ｘｙは、Ｉ_ｐおよびＩＩ_ｐとともに下表１に記載されている。表１は、基本波動関数の一定重み付け（Ｆ（ｋ_ｚ）＝１）を実施する一実施形態においてこれらの関数が有する形も含んでいる。

上表において、Ｊ_０およびＪ_１はそれぞれ０次および１次のベッセル関数である。本発明の発明者は、積分Ｉ_ｐおよびＩ_ｘｙの実数部に対する解析的式だけを見つけることができたが、これらの積分は一般に、例えば適当な数値積分法、例えばガウス型アルゴリズムを使用した数値積分によって計算することができる。上記の積分のうちの３つの積分における無限大まで及ぶ積分の代わりに、指数関数が十分に低いレベルに達する点においてそれらの積分の上限を打ち切ることができ、または測定格子内の頻度および間隔に基づいて上限を選択することができる。

上述のとおり、変数ρおよびζの離散値の２Ｄテーブルを（例えば数値積分によって）作成し、記憶することによって、多数の積分をかなり高速に実施することができる。特定のρおよびζ値に対する積分を計算するため、次いで、このテーブルの補間を実行しなければならない。あるいは、またはこれに加えて、異なる領域をカバーする２Ｄ近似関数を構築してもよい。

積分Ｉ_ｐ、Ｉ_ｘｙおよびＩ_ｚは特異性（ｓｉｎｇｕｌａｒｉｔｙ）を持たず、そのため、これらの積分は、直接にテーブル形態に入れ、補間することができる。両方のパラメータが０の場合、すなわちこのパラメータ・セットの原点（ρ，ζ）＝（０，０）にある場合、積分ＩＩ_ｐ、ＩＩ_ｘｙおよびＩＩ_ｚは特異性を有する。補間スキームがこの原点の近くで適正に機能するようにするためには、テーブル化されたこれらの関数の中の特異性を抽出すべきである。

波動関数の方向性重み付けの場合を考えると、非常に小さなρおよびζ値に関しては以下のようになる。

これらの特異性は、補間の間に以下のように除去することができる。ＩＩ_ｐの場合には、テーブルに記憶する一組のパラメータの組合せに対する数値積分によってこの積分を解き、それらの値をテーブルに記憶する前に、それらの値に
を掛ける。このようにして、スケーリングされた積分に対して補間を実行する。補間の後、その結果を
で割る。残りの積分も同様に処理することができる。単位重み付けを用いて、特異性の同様の除去を、それらをテーブル化する前、補間の間に実行する。

全ての仮想源平面と平行なそれぞれの平面内に測定点および計算点を配置すると、効率をさらに向上させることができることに言及することも興味深い。このような構成では、それぞれの周波数に対して第２のパラメータζの値が少数だけ生じ、したがって、それぞれのテーブルがそれぞれのζ値に対応する１Ｄテーブルを作成することが魅力的になる。例えば、単一の仮想源平面、単一の平行な測定平面および測定平面とは異なる単一の平行な計算平面を用いて、２組の１Ｄテーブルを作成することができる。これは、実行中に、それぞれの周波数に対して効率的に実行することができる。同様に、１つ仮想源平面、２つの平行な測定平面および別個の平行な計算平面を用いて、全ての周波数に対して４組の１Ｄテーブルを作成することができる。例えば、これらの１Ｄテーブルは、予め計算した２Ｄテーブルの補間によって有利に作成することができる。

一般に、上述のとおり、測定位置において音圧が測定されているとき、標的位置ｒの圧力は、
ｐ（ｒ）＝ｐ^Ｔ（Ａ^ＨＡ＋εＩ）^−１Ａ^Ｈα（ｒ）＝（Ａ^Ｈα（ｒ））^Ｔ（（Ａ^ＨＡ）^Ｔ＋εＩ）^−１ｐ
として得られる。上式で、ｐは、所与の位置における測定／既知音圧データのベクトル、εは正則化パラメータである。同様の式
からχ方向の粒子速度を得ることができる。上記の公式は、チホノフ型の正則化を使用しているが、同様の正則化スキーム、例えばいわゆる強化チホノフ正則化（Ｅｎｈａｎｃｅｄ Ｔｉｋｈｏｎｏｖ ｒｅｇｕｌａｒｉｚａｔｉｏｎ）を使用した別の公式をしようすることもできる。

次に、図２を参照し、図１を続けて参照して、上述の方法論に基づく再構成法の一実施形態を説明する。
この方法の最初の初期設定ステップＳ１では、適当な座標系および仮想源平面１０１を定義する。

一般に、試験準備の間に、調査対象の物体の実際の音源面に沿って位置するように測定位置を選択することができる。例えば、フリー・フィールド条件および平面変換器アレイを仮定すると、実際の音源面の一部分にほぼ平行な測定平面を選択することができる。同様に、３Ｄアレイを使用する非フリー・フィールド条件の場合には、アレイがかなり扁平であると仮定し、実際の音源面の一部分と平行に近くなるように、アレイの中央平面（ｍｅｄｉａｎ ｐｌａｎｅ）を配置することができる。

ステップＳ１において、この方法は、例えば上述の適当な位置検出システムによって、測定点の位置に関する情報、すなわち一組の位置ｒ_ｉ、ｉ＝１，２，．．．，Ｉの座標を受け取ることができる。この方法は次いで、その測定ジオメトリに対する座標系を選択することができる。例えば、平面変換器アレイに関しては、ｘｙ平面を、測定平面内に位置するように（すなわち測定平面がｚ＝０平面に位置するように）定義し、原点が平面アレイの中心に位置し、正のｚ軸が、音源面から離れる方向、すなわち源のない領域内へ向かう方向を指すと定義することができる。同様に、３Ｄアレイに関しては、ｚ＝０平面を、アレイの中央平面となるように定義することができる。しかしながら、別の方法で座標系を選択することもできることが理解される。

初期設定ステップＳ１の間に、この方法はさらに、選択した座標系に対して、例えばｘｙ平面に平行な仮想源平面１０１を定義することができる。この仮想源平面の位置決めは、自動的に実行し、または少なくとも部分的にユーザ入力に基づいて実行することができる。

用途によっては、実際の源ジオメトリが分かっていない、または少なくとも正確には分かっていないことがある。このような状況では、アレイからの所定の距離のところにあるアレイに平行な平面上において音源を再構成することができる。この場合には、アレイと平行な仮想源平面との間の距離を、この方法への入力パラメータ、例えばユーザ定義の入力パラメータとして指定することができる。

他の状況では、実際の源ジオメトリが既知であり、それをインポート（ｉｍｐｏｒｔ）し、整列させることができることがあり、または、アレイ位置を測定するために使用する位置測定システムを用いて、実際の源ジオメトリを測定することもできる。アレイ位置と物体ジオメトリの両方が既知である場合には、処理ソフトウェアが、これら２つの間の距離を自動的に決定することができる。これに基づいて、この方法は、アレイが平面であると仮定して、例えば実際の音源面の内側の指定された距離のところに、アレイに平行な仮想源平面を自動的に位置決めすることができる。アレイが平面でない場合には、中央平面を定義することができる。

一般に、仮想源平面は、物体内の物体の表面から短い距離のところに定義することができ、選択する距離は、測定ジオメトリ、例えば隣接する測定点間の典型的な距離または平均距離に依存することがある。例えば、仮想源平面は、音源面から源物体内へ、すなわち測定点の位置および音場を再構成する標的位置（１つまたは複数）とは反対側の実際の音源面の側に、最も近い隣接測定点間の平均距離よりも大きい距離（例えば約１倍から約３倍）だけ離れるように、選択することができる。最適な正確さは、実際の音源面から後ろに測定格子間隔の約１．５倍の距離を置いたところ、すなわち物理的な音源面よりも測定格子から離れたところに仮想源平面を配置することによって、得ることができることが分かっている。低周波数では、これよりもわずかに大きな距離がしばしばより良好な結果を与えることが分かっている。この再構成法の目的上、仮想源平面の片側は、源がなく均質であるとして取り扱われ、そのため、その仮想源平面に関連付けられた一組の平面伝搬波およびエバネッセント波によって、音場を表現することができる。

ステップＳ１の初期設定はさらに、適用する重み付けスキーム、適用する正則化スキーム、使用する測定システムのタイプおよび特性、例えば変換器の数、音場を再構成する標的位置の座標、計算する１つまたは複数の音響特性、周波数範囲、関数テーブルの補間する部分の初期設定（他の部分には、この方法を実行するソフトウェアを供給することができる）など、１つまたは複数のパラメータの設定を含むことができる。例えば、解析ユニット１０３の適当なユーザ・インタフェースを介して、これらのパラメータおよび追加のまたは代替のパラメータをユーザが定義することができる。

ステップＳ２において、この方法は、一組の位置ｒ_ｉにおいて測定したそれぞれの複素時間調和音圧ｐ（ｒ_ｉ）を受け取る。ｒ_ｉは、上記の座標系に対するそれらの位置、例えば位置検出デバイス１１０から受け取った入力から決定される変換器位置の座標ベクトルを表す。この説明の目的上、位置ｒ_ｉは、空気、水などの均質流体によって占められている源のない領域Ωに位置すると仮定する。

上で論じたフリー・フィールド状況では、全ての源が平面ｚ＝０の一方の側にあり、原則として、例えばこの平面内での測定に基づいて音場を完全に再構成することができるが、本明細書に記載した方法を他の任意の３Ｄ測定格子に関して適用することもできる。さらに、後述するとおり、この方法を、源平面が複数ある状況に適用することもできる。

複素時間調和音圧ｐ（ｒ_ｉ）は、受け取った変換器信号からさまざまな方法で決定することができる。例えば、全ての変換器を用いた１回の同時記録を実施し、続いて、全てのＦＦＴ線（ＦＴＴ ｌｉｎｅ）に対する複素（時間調和）音場を提供する全ての信号のＦＦＴを実施することができる。あるいは、基準信号から全てのアレイ変換器への相互スペクトル（ｃｒｏｓｓ ｓｐｅｃｔｒｕｍ）および基準信号の自己スペクトル（ａｕｔｏ−ｓｐｅｃｔｒｕｍ）の測定を実行することもできる。相互スペクトルを自己スペクトルの平方根で割ったものは、適用したスペクトルの全ての周波数線に対する複素（時間調和）音場を提供する。図２の方法は、測定音圧に基づくΩにおける音場の再構成を実行する。代替実施形態は、音圧勾配、粒子速度またはこれらの組合せなどの１つまたは複数の代替のまたは追加の音響量を測定することができる。

ステップＳ３において、この方法は、式（２４）を使用して、すなわち
および
でそれぞれ評価した２次元関数Ｉ_ｐとＩＩ_ｐの１次結合として、自己相関行列Ａ^ＨＡを計算する。そのために、この方法は、ρおよびζに対するそれぞれの値を、それぞれＩ_ｐおよびＩＩ_ｐの予め計算した関数値のルックアップ・テーブル１１２への索引として使用する。このルックアップ・テーブルを、解析ユニットの記憶領域または他の適当な記憶媒体に記憶することができる。このようにしてこの方法は、ρおよびζに最も近い記憶された関数値を得ることができ、ρおよびζの２次元１次補間、または他の適当な補間スキーム、例えば３次補間を実行することができる。

あるいは、この方法は、その代わりに、転置行列ＢＢ^Ｈまたは他の適当な変換行列を計算することができることが理解される。

ステップＳ４において、この方法は、ベクトルｑを、
ｑ＝（（Ａ^ＨＡ）^Ｔ＋εＩ）^−１ｐ
またはその転置
ｑ^Ｔ＝ｐ^Ｔ（Ａ^ＨＡ＋εＩ）^−１
として計算する。

ここで、εは、上述のとおり、正則化パラメータである。この正則化パラメータは、指定された信号対雑音比（ＳＮＲ）から計算することができ、または雑音が多い（ｎｏｉｓｙ）測定データに基づいて、一般化交差検証（ＧＣＶ）、Ｌ曲線法または同種のパラメータ推定手順を使用して、自動的に計算することができる。ＳＮＲが指定されているときには、行列ＢＢ^Ｈ内の対角要素の平均値または最大値に１０＾（ＳＮＲ／１０）を掛けた積として、このパラメータεを得ることができる。フリー・フィールド条件の場合、単一の仮想源平面だけが使用されるときには全ての対角要素が等しく、そのため、平均値または最大値の代わりに単一の対角要素の値を使用することができる。

上記の行列演算は、行列を逆行列にするなど行列方程式を解く適当な任意の数値技法によって実行することができる。

次いで、この方法は、計算されたｑ値から音圧を計算することができる。音圧の代わりに、または音圧に加えて、この方法は、別の音響量、例えば上述の粒子速度を計算することができることが理解される。圧力および粒子速度の計算値から音の強さを得ることができる。

したがって、この方法は、音場を再構成する一組の標的点ｒのそれぞれの標的点における音圧および／または粒子速度を計算することができる。標的点の２つの例が図１の１０９に示されている。多くの状況において、調査対象の物体の音源面で音場を再構成することが重要である。例えば、この方法は、規則格子の一組の格子点のそれぞれの格子点で音圧を計算することができる。この説明の目的上、再構成された音場値を計算する点を計算点と呼ぶ。

この方法は、全ての計算点を循環し、所望の音響量を計算することができる。具体的には、ステップＳ５において、この方法は、現在の計算位置ｒで所望の音響量を計算する。

例えば、この方法は、ベクトル
φｐ（ｒ）≡Ａ^Ｈα（ｒ）
および／または
Ψ_χ（ｒ）≡Ａ^Ｈβ_χ（ｒ）
を計算することができ、音圧を
ｐ（ｒ）＝ｑ^Ｔφ（ｒ）
として計算し、かつ／または粒子速度を
ｕ_χ（ｒ）＝ｑ^ＴΨ_χ（ｒ）
として計算する。

この方法は、式（２５）および（２８）を使用して、それぞれ相関ベクトルＡ^Ｈα（ｒ）およびＡ^Ｈβ_χ（ｒ）を、上記ステップＳ３におけるＡ^ＨＡの計算に関して説明したのと同様の方法で、すなわちそれぞれ入力パラメータρおよびζのそれぞれの値で評価した２次元関数Ｉ_ｐ、ＩＩ_ｐおよびＩ_ｚ、ＩＩ_ｚまたはＩ_ｘｙ、ＩＩ_ｘｙの１次結合として計算する。上記のとおり、関数を評価する入力パラメータの値は、２Ｄ関数の予め計算された関数値を得るためのそれぞれのルックアップ・テーブル１１２への索引として使用され、それぞれテーブル化された値Ｉ_ｐ、ＩＩ_ｐおよびＩ_ｚ、ＩＩ_ｚまたはＩ_ｘｙ、ＩＩ_ｘｙ間の３次補間または他のタイプの補間を実行することによって使用される。

ステップＳ６において、この方法は、未処理の計算点が残っているかどうかを判定する。残っている場合には、この方法はステップＳ５へ戻り、残っていない場合にはステップＳ７へ進み、ステップＳ７において、計算した一組の音圧および／または粒子速度を、例えば図形表現として、一連の数値として、および／または他の適当な形式で出力する。あるいは、またはこれに加えて、音圧および粒子速度から音の強さが容易に得られる。

以上では、図３ａに示す状況のフリー・フィールド条件に関する一実施形態を説明した。しかしながら、本明細書に記載した方法の実施形態は、別の状況においても、例えば源および反射物体が半空間だけに限定されないときにも使用することができる。以下では、非フリー・フィールド条件に対する平面波公式を提供する、したがって式（１６）に定義した平面波波動関数が測定領域内の音場に対する完全な基礎をもはや構成しない状況を考慮する、上記の実施形態に対する変更を説明する。

最初に、位置ｚ^＋およびｚ⁻、ｚ^＋≦ｚ≦ｚ⁻にある平行な２平面１０１ａと１０１ｂの間の均質な源のない媒質中で測定を実施する図３ｂに示すケースを考える。ｚ≦ｚ^＋に関する源を有する音場の部分は、式（１６）の基礎関数によって表現することができる。さらにｚ≧ｚ⁻に関する源を有する音場成分を表現するため、結合された組を与える対応する一組の平面波波動関数を追加する。

上式では、式（１７）のｋ_ｚの定義が依然として使用され、両方の組に対して全方向重み付けが適用される。

したがって、仮想源平面の追加は、ＳＯＮＡＨ（またはＨＥＬＳ）法において一組の平面波波動関数を追加することに対応する。式（１３）、（１４）、（２０）および（２１）は、行列Ａ^ＨＡおよびＡ^Ｈα（ｒ）の要素に２つの組からの寄与を追加することによって、結合された一組の基礎関数（３０）および（３１）を考慮したＳＯＮＡＨ推定手順を得ることができることを示している。これを実行した結果は次のとおりである。

粒子速度を推定するに必要なベクトルＡ^Ｈβ_χ（ｒ）は、フリー・フィールド条件の場合とちょうど同じように導出することができ、結果は以下のとおりである。

上記の公式は測定位置の選択を制限しないが、正確な予測を得るのに平面アレイでは不十分なことがある。このことは、測定平面内において対称な像である源は、それぞれのマイクロホンにおいて等しい圧力を発生させるため、互いに区別することができないことによって、例証することができる。一般に、全く同じ平行な２つの長方形アレイからなる２重層アレイを使用することができる。アレイの両側に源があるときの波に対する重みは、定数係数Ｆ_０，ｚ＋およびＦ_０，ｚ−によって調整することができる。一般に、これらのパラメータの値はともに１として選択されるが、アレイの片側に有意な源が存在しない場合には、これらのパラメータの一方をゼロに設定し、それによって一方の側だけに源が分布するという情報をこの公式に提供することによって、より良好な推定を達成することができる。

無限平面セクション（例えばスラブ（ｓｌａｂ）またはスライス（ｓｌｉｃｅ））の両側に源があるケースから、箱形体積の６つの全ての面に源があるケースへの拡張を、図３ｃの断面図に示す。この状況では単に、箱のそれぞれの側に属する平面波スペクトルからの行列要素への寄与を追加する。対応する式を下表２に示す。

それぞれの面は半球からの全ての平面波を等しい重みで含むため、２面にだけ源がある場合の式は、全ての方向からの平面伝搬波を考慮することが認められることがある。しかし、これらの２つの面が、全てのタイプのエバネッセント波を含むというわけではない。ｘおよびｙ方向に減衰する波は失われる。

下記の実施例は、波がｘｙ平面において減衰することを考慮することによって達成可能な正確さの向上を例示するシミュレートした測定を示す。

したがって、表１に示した２Ｄ関数と同じ２Ｄ関数を使用して、図２の方法を、測定位置の複数の面に源があるケースに適用することができる。自己および相互相関行列の要素を構成する相関関数は依然として、これらの２Ｄ関数の関数値の１次結合として得ることができる。

固有推定誤差
以下では、計算点において、本明細書に開示した再構成法の平面波波動関数の値が、平均して、どれくらいよく予測されるのかの少なくとも１つの推定を決定する方法の一実施形態を説明する。

最初に、フリー・フィールド条件を考えると、式（１６）によって、一組の平面伝搬波関数およびエバネッセント波関数が与えられる。ベクトルα（ｒ）は、式（７）の中の定義に従って、計算点ｒにおけるこれらの全ての基本波の真の値を含む。式（１８）、（１２）および（２）は、行列Ａの１つの行が、測定位置を横切る波動関数のうちの１つの波動関数の値を含むことを示し、したがって、式（９）および（１１）から、下式のように定義されるベクトル

は、計算点における波動関数の推定値を含むことになる。このベクトルの２ノルムは、予測の平均振幅の測度（ｍｅａｓｕｒｅ）である。

平均過小推定（ａｖｅｒａｇｅ ｕｎｄｅｒｅｓｔｉｍａｔｉｏｎ）は、予測ノルム
が真のノルム‖α（ｒ）‖_２よりも小さいことが特徴である。真のノルム‖α（ｒ）‖_２を見つけるため、式（１８）、（１２）および（２）から、ちょうどα（ｒ）が計算点における全ての関数値を含むように、行列Ａの１つの列が、マイクロホン位置における全ての波動関数の真の値を含むことが分かる。したがって、単にマイクロホン位置ｒ_ｉを計算位置ｒに置き換えることによって、行列Ａ^ＨＡの対角要素［Ａ^ＨＡ］_ｉｊの式（２４）を使用して、α^Ｈαの値を計算することができる。

ここで、定義（２２）および（２３）を使用して、上記の２つの積分に対する解析的式を容易に検証することができる。固有平均振幅利得（ｉｎｈｅｒｅｎｔ ａｖｅｒａｇｅ ａｍｐｌｉｔｕｄｅ ｇａｉｎ）は下式のように定義される。

推定誤差ベクトル
の２ノルムは、式（３７）および（３８）を適用することによって以下のように計算することができ、

最終的に、固有平均推定誤差水準は、以下のように定義される。

単純に、行列Ａ^ＨＡおよびＡ^Ｈα（ｒ）の要素に対する拡張式、例えば２つの面に源がある場合の式（３３）および（３４）を適用することにより、以上の説明を、非フリー・フィールド条件に拡張することができる。この場合には、ちょうど式（２４）を式（３３）に拡張したときのように、α^Ｈαに対する式（３９）を次のように拡張しなければならない。

特定の源構成に対して認められる相対予測誤差は、以下の２つの理由から、上記の固有誤差とは異なる。

１）その特定の音場の平面波振幅スペクトルは、この方法が最適化される平面波振幅スペクトルとは異なる。フリー・フィールド条件に関しては、この最適振幅スペクトルが、仮定した仮想源平面ｚ＝ｚ^＋上で｜Ｆ（ｋ_ｚ）｜の形を有する。
２）特定の測定音場に関しては、予測誤差が、平面波成分全体にわたるＲＭＳとしてではなく、全体の音圧に関して計算される。

これらの２つの差異のそれぞれの効果については後に示す。

実施例：予測誤差の数値的検討
以下では、本明細書に記載した方法の予測の妥当性領域に関する情報の概要を、いくつかの典型的なアレイ・ジオメトリについて、周波数の関数として説明する。

測定誤差は考慮せず、正則化パラメータの自動決定法（ＧＣＶ、Ｌ曲線、．．．）についても調査しない。この実施例の目的上、式（１１）および（２６）で使用する正則化パラメータεは、指定されたダイナミック・レンジＤ（信号対雑音比）から、下式を使用して決定される（例えばＲ．ＳｔｅｉｎｅｒおよびＪ．Ｈａｌｄ、「Ｎｅａｒ−ｆｉｅｌｄ Ａｃｏｕｓｔｉｃａｌ Ｈｏｌｏｇｒａｐｈｙ ｗｉｔｈｏｕｔ ｔｈｅ ｅｒｒｏｒｓ ａｎｄ ｌｉｍｉｔａｔｉｏｎｓ ｃａｕｓｅｄ ｂｙ ｔｈｅ ｕｓｅ ｏｆ ｓｐａｔｉａｌ ＤＦＴ」、Ｉｎｔｅｒｎ．Ｊ．Ａｃｏｕｓｔ．Ｖｉｂ．６、８３〜８９（２００１）を参照されたい）。
ε＝［Ａ^ＨＡ］_ｉｊ１０^{−Ｄ／１０} （４４）
上式で、［Ａ^ＨＡ］_ｉｊは、行列Ａ^ＨＡの対角要素である。この関係は、平面アレイおよびフリー・フィールド条件に対して導出したものであり、その場合、全ての対角要素は等しい。ここでは、一般に、対角要素の平均値を使用する。Ｄの値は一般に１５から３０ｄＢの範囲にある。別の値が指定されない限り、以下ではＤ＝３０ｄＢを使用する。

実施例Ａ：測定平面内の正確な圧力表現の領域
最初に、２つの異なる不規則平面アレイ設計のアレイ平面ｚ＝０における総称（ｇｅｎｅｒｉｃ）固有誤差分布を調べる。これらの２つの異なる不規則平面アレイ設計の直径はともに０．７５ｍであり、これらの設計はともに、少なくとも２０ｋＨｚまでのビーム形成（ｂｅａｍｆｏｒｍｉｎｇ）用途に対して最適化されている。一方の設計は、ホイール内のスポークとして配置された全く同じ１３個の６要素線アレイからなるスポーク・ホイール・アレイ（Ｓｐｏｋｅ Ｗｈｅｅｌ ａｒｒａｙ）である。図５を参照されたい。もう一方の設計は、全く同じ７つの扇形からなり、それぞれの扇形に１２個のマイクロホンが不規則だが一様に分布したセクタ・ホイール・アレイ（Ｓｅｃｔｏｒ Ｗｈｅｅｌ Ａｒｒａｙ）（例えばＪ．Ｈａｌｄ、「Ａｒｒａｙ ｄｅｓｉｇｎｓ ｏｐｔｉｍｉｚｅｄ ｆｏｒ ｂｏｔｈ ｌｏｗ−ｆｒｅｑｕｅｎｃｙ ＮＡＨ ａｎｄ ｈｉｇｈ−ｆｒｅｑｕｅｎｃｙ ｂｅａｍｆｏｒｍｉｎｇ」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ Ｉｎｔｅｒ−Ｎｏｉｓｅ ２００４を参照されたい）である。ＳＯＮＡＨの適用に関して良好な成果を得るため、セクタ・ホイール・アレイの一様分布は、アレイ・ジオメトリ最適化の間中、維持した。これらのアレイは平面であるため、フリー・フィールド条件を仮定し、単一の仮想源平面だけを使用する。平均要素間隔は、要素あたりの平均アレイ面積を計算することによって定義し、両アレイ・ジオメトリについて、平均要素間隔は約７ｃｍである。規則的な長方形アレイに関して、この間隔は、約２ｋＨｚまでの周波数範囲をサポートすると考えられ、その周波数まで空間エイリアシング（ｓｐａｔｉａｌ ａｌｉａｓｉｎｇ）を回避するため、アレイ格子間隔よりも小さくないスタンドオフ距離（ｓｔａｎｄｏｆｆ ｄｉｓｔａｎｃｅ）を選択すべきである。仮想源平面は、ｚ_０＝−７ｃｍにある実際の音源面から、マイクロホン間隔およそ１つ分後ろに配置すべきであるため、アレイ平面から１４ｃｍ、すなわちｚ^＋＝−１４ｃｍのところの仮想源平面を選択する。

図５の輪郭プロットは、アレイ平面内の音圧予測の固有相対誤差水準Ｅ（ｒ）を示し、マイクロホン位置は丸によって示されている。白色のエリアは、固有誤差水準Ｅ（ｒ）が−２０ｄＢ未満であることを示し、灰色の輪郭は、黒で表現された−６ｄＢよりも大きな誤差水準への２ｄＢ刻みでの正確さの低下を示す。低周波数では、両方アレイが、アレイ・エリア全体にわたって良好な音圧表現を提供することができるが、１ｋＨｚでは既に、スポーク・ホイール・アレイのマイクロホン分布の隙間が大きすぎて、それらの隙間において良好な表現を提供することができない。セクタ・ホイール・アレイは、２ｋＨｚまでずっと、アレイ・エリア全体にわたって許容できる圧力表現を提供するが、スポーク・ホイール・アレイ・ジオメトリに関してはそうでないことが明らかである。しかし、スポーク・ホイール・アレイはビーム形成用途に対して非常に有効であり、大きなサイズでは、その支持構造の構築がより容易である。

実施例Ｂ：測定平面の外側の正確な圧力表現の領域
次に、アレイ平面に対して垂直な平面内における正確な音圧予測の領域を調べる。ビーム形成用に設計した不規則アレイの代わりに、ｘおよびｙ方向の要素間隔が３ｃｍで、０≦ｘ≦２１ｃｍ、０≦ｙ≦２１ｃｍのエリアをカバーする規則的な８×８要素の単層（８×８×１）および２重層（８×８×１）アレイを考える。これらの両タイプのアレイはｘｙ平面ｚ＝０にマイクロホン層を有し、２層アレイは、ｚ＝３ｃｍに全く同じ追加の層を有する。アレイの対称平面であるｘｚ平面ｙ＝１０．５ｃｍにおける相対音圧予測誤差に注目する。ちょうど不規則平面アレイと同様に、単層規則アレイに関しては、アレイから、アレイ要素間隔の２倍に等しい距離にある仮想源平面、すなわちｚ^＋＝−６ｃｍにある仮想源平面を選択する。２重層アレイに関しては、前（ｆｒｏｎｔ）仮想源平面に対しては同じ距離、後（ｒｅａｒ）仮想源平面に対しては、要素間隔の１９倍の距離、すなわちｚ^＋＝−６ｃｍ、ｚ⁻＝６０ｃｍを選択する。

図６の左列は、（単一の仮想源平面を有する）単層アレイの固有誤差水準Ｅ（ｒ）を示し、右列は、（２つの仮想源平面を有する）２重層アレイの固有誤差水準Ｅ（ｒ）を示す。中央の列は、特定の音場、すなわち（ｘ，ｙ，ｚ）＝（７．５，７．５，−６．０）ｃｍにおける単極子点源の音場に対する、単層アレイによる音圧予測の相対誤差を示す。このマッピング・エリアは、ｘ軸に沿って両方の方向に、アレイ・エリアを過ぎて６ｃｍ延び、ｚ方向には、ｚ＝−３ｃｍから、正のｚ方向にアレイ要素間隔の４倍の距離まで延びる。明らかに、特定の音場予測の誤差分布は、総称固有誤差よりも不規則であり、一般に総称固有誤差よりも低い。あるエリアでは、固有誤差水準によって予測される誤差よりもかなり小さい誤差が達成される。これは、異なる平面波成分からの誤差寄与が互いに打ち消しあうためと考えることができる。固有誤差水準はこのような相殺効果を回避し、本発明の方法の推定の正確さのはるかに一般的な全体像を提供する。２重層アレイは、単層アレイに比べて、正確な表現の領域をかなり拡張するように見えるが、追加の仮想源平面を追加することによって、音場の自由度も増大した。低い予測誤差を有する領域のｚ方向の広がりがいくぶん制限される理由は、ＳＯＮＡＨが、ｘｙ方向の源分布については何も仮定していないためである。良好な遠距離場推定を提供するため、源エリアのサイズについてのなんらかの仮定をこの方法に組み込むことができる。これは例えば、平面伝搬波およびエバネッセント波の代わりに単極子源の限定された分布に基づくＳＯＮＡＨ公式を使用することによって達成することができる。

実施例Ｃ：パラメータの変更に対する予測の正確さの感度
パラメータの変更に対する感度の調査におけるデータ量を限定するため、アレイに垂直な中心に沿った予測誤差だけに注目する。実施例Ｂの単層８×８×１アレイだけを考える。図７ａは、図６の左列に示した条件と同じ条件に対する固有誤差水準を示し、そのため、図６の左列の輪郭プロットの中央を通る一組の水平スライスだけしか示されていない。しかし、図７ａにはいくつかの周波数、例えば間隔３ｃｍの規則アレイ格子のより高い５ｋＨｚの限界周波数が追加されている。明らかに、５ｋＨｚでは、測定平面でさえも、固有誤差レベルの急速な増大が始まっているが、依然として許容できる予測を達成することができる。図７ｂは、対応する振幅利得Ａ（ｒ）の曲線を示す。５ｋＨｚでさえも、アレイ平面からアレイ格子間隔１つ分の範囲内では、依然として振幅は正確に予測されるが、その間隔の外側では、次第に増大する過小推定が観察される。

図７では、ダイナミック・レンジＤは３０ｄＢに固定した式（４４）から正則化パラメータを決定した。次のステップは、ダイナミック・レンジＤを変化させたときに予測誤差がどれくらい変化するのかを見ることである。図８の左列は、Ｄの変化に対する固有誤差水準Ｅの感度を、低周波数（５００Ｈｚ）および高周波数（４ｋＨｚ）について示す。周波数ごとに、１０ｄＢから４０ｄＢの範囲の対応するＤの値が記された一組の曲線が示されている。一般に、測定平面よりも源領域から離れた位置での圧力の予測は感度が低く、一方、源領域により近い位置での予測は、Ｄ値の増大（より小さな正則化パラメータ値に対応する）とともに向上し、特に低周波数でそうである。より小さなＤ値はエバネッセント波の再構成を制限し、エバネッセント波は、低周波数においておよび源領域の近くで音場のより大きな部分を構成することから、このことを予想することができる。図８の右列は、ＳＯＮＡＨ推定における対応する固有振幅利得Ａの曲線を示す。５００Ｈｚでは、源領域の直ぐ近くを除いて、Ｄの変化に対する振幅予測の感度はあまり高くないが、４ｋＨｚでは、Ｄ値が約２０ｄＢよりも小さくなったときに、全ての位置で（源領域から離れた位置でも）、次第に増大する過小推定が観察される。このことは、高周波数では平面伝搬波が平面波スペクトルのより大きい部分を構成し、したがってＳＯＮＡＨ自己相関行列Ａ^ＨＡの特異値スペクトルのより大きい部分を構成することによって説明されるであろう。Ｄがある点よりも小さくなる（正則化パラメータが増大する）と、伝搬波成分の部分さえも、正則化によってかなり減衰する。

以前に定義した固有振幅利得Ａおよび誤差水準Ｅは、推定位置における誤差を、全ての平面伝搬波およびエバネッセント波にわたって平均したものを表し、エバネッセント波成分には、ＳＯＮＡＨ推定アルゴリズムの導出に使用するのと同じ重みが与えられる。どちらの場合も、波動関数は、同じ仮想源平面（１つまたは複数）（フリー・フィールド条件ではｚ＝ｚ^＋）内においてスケーリングされる。次に調べるべき課題は、ＳＯＮＡＨ計算において使用する仮想源距離が、測定用途においてエバネッセント波の実際の内容を表していない場合に、誤差水準がどれくらい変化するかということである。フリー・フィールド条件および単一の仮想源平面の使用だけを考える。この調査を実施するため、式（３８）および（４１）において、相関行列Ａ^ＨＡおよびＡ^Ｈαが推定誤差を平均する平面波場を表し、推定重みベクトルｃ（ｒ）がＳＯＮＡＨ推定を表すことに注目する。したがって、１）（３８）および（４１）でのＡ^ＨＡおよびＡ^Ｈαの計算ならびに２）ＳＯＮＡＨ推定重みベクトルｃ（ｒ）の計算に対して異なる仮想源距離を使用することによって、１）測定におけるエバネッセント波の実際の内容と２）ＳＯＮＡＨ計算におけるエバネッセント波の重みとの間の不整合の程度の変化に対するＳＯＮＡＨ推定の感度を調べることができる。図９に示す結果に関しては、実際の音場に対して、すなわち式（３８）および（４１）で使用するＡ^ＨＡおよびＡ^Ｈαの計算において、以前に使用した仮想源平面ｚ＝ｚ^＋＝−６ｃｍを維持し、ＳＯＮＡＨ処理に対しては、変化する仮想源平面距離ｄを使用した。大きなｄ値は、ＳＯＮＡＨ推定におけるエバネッセント波に対する重みを小さくする。図９の左列は、ｄの変化に対する固有誤差水準Ｅの感度を、低周波数（５００Ｈｚ）および高周波数（４ｋＨｚ）について示す。周波数ごとに、３ｃｍから１２ｃｍの範囲の対応するｄの値が記された一組の曲線が示されている。右列は、対応する振幅利得の曲線を含む。全ての場合に、正則化パラメータは、ダイナミック・レンジＤを３０ｄＢに固定することによって指定される。

５００Ｈｚでは、過度に小さな仮想源距離、すなわち再構成におけるエバネッセント波に対する過度に大きな重みに対する誤差水準の感度はいくぶん低いが、過度に大きな距離は、最終的に、測定平面よりも源領域から離れた位置において誤差水準のわずかな過大推定（Ａ＞１）をもたらす。これは、波数漏れ（ｗａｖｅ−ｎｕｍｂｅｒ ｌｅａｋａｇｅ）によって、小さなエバネッセント波成分が伝搬成分として処理され、それによってそのエバネッセント波成分がより離れた場に寄与することによって説明される。測定点の数がいくぶん少ないため、伝搬波成分とエバネッセント波成分との間の区別はあまり鮮明ではない。（重みベクトルｃ（ｒ）によって表現される）ＳＯＮＡＨ処理が、エバネッセント波に非常に小さな重みを与えると、エバネッセント波成分の伝搬成分への漏れは増大する。ｄが大きいときにも起こるｚ＜０に対する誤差の増大に注目すると、説明は異なってくる。ＳＯＮＡＨ計算行列においてエバネッセント波には過度に小さな重みが与えられるため、選択された正則化パラメータは、これらの波成分の一部が十分に再構成されることを妨げる。４ｋＨｚにおいて、エバネッセント波は一般に、音場のはるかに小さな部分を構成し、そのため、ＳＯＮＡＨ計算においてエバネッセント波により大きなまたはより小さな重みを与えることに対する感度はより小さい。最大の効果は、再構成において、実際の平面波スペクトルが示唆するものよりも大きな重みをエバネッセント波に与える過度に小さな仮想源距離を使用するときの振幅利得に対する効果である。源領域の近くでは過大推定が観察され、より遠くでは、通常よりも大きな過小推定が見られる。これは、５００Ｈｚでの過大推定に対して示した説明と同じ説明によって説明され、ここでは、漏れが反対方向に進む。一部の小さな伝搬波成分がエバネッセント波として処理され、したがってそれらは源エリアから離れると減衰し、源エリアの近くで増幅される。

実施例Ｄ：アレイ平面内の源−６つの仮想源平面の使用
図１０は、２つの面だけからではなく６つの全ての面から平面波波動関数を追加し、それによってｘｙ平面内の方向に沿って減衰するエバネッセント波成分の適当なモデル化を可能にすることによって得られる正確さの向上の一例を示す。実施例Ｂに関して説明したものと同じ８×８×２ ２重層アレイ・ジオメトリを用いたシミュレートされた測定を実行する。（ｘ，ｙ，ｚ）＝（−６．０，７．５，１．５）ｃｍ、すなわち２つのアレイ層間の中心平面内のアレイ・エリアから６ｃｍ外側に、単一の単極子点源を配置する。明らかに、この単極子は、ｘｙ平面内の方向に沿って減衰するエバネッセント波型を生成する。ＳＯＮＡＨを使用して、前アレイ平面ｚ＝０ｃｍ内の測定位置における強さベクトル（ｉｎｔｅｎｓｉｔｙ ｖｅｃｔｏｒ）を、１２８個の測定圧力値に基づいて予測する。次いで、下式を使用して相対平均誤差水準を計算する。

上式で、
は、計算点番号ｉにおける真の強さベクトル、Ｉ_ｉは、同じベクトルのＳＯＮＡＨ推定値であり、求和はｉ全体についてである。ダイナミック・レンジ３０ｄＢ、Ｄ＝３０ｄＢを使用してＳＯＮＡＨ計算を実行し、全ての方向にアレイから６ｃｍ離れた以下の仮想源平面を選択する。
ｚ⁻＝９ｃｍ、ｚ^＋＝−６ｃｍ、ｘ⁻＝２７ｃｍ、ｘ^＋＝−６ｃｍ、ｙ⁻＝２７ｃｍ、ｙ^＋＝−６ｃｍ

したがって、単極子は、これらの６つの仮想源平面によって定義された箱の境界上に位置する。図１０を見ると、低周波数における強さ推定は、ｘｙ平面に平行な仮想源平面だけ存在するときよりも、全ての６つの仮想源平面が存在するときの方がかなり正確であることが分かる。しかしながら、周波数が増大するにつれてこの向上は消滅する。これは、高周波数ではエバネッセント波が音場のより小さな部分を構成するためである。

いくつかの実施形態を詳細に説明し、示したが、本発明はそれらの実施形態に限定されるもではなく、下記の特許請求の範囲に定義された主題の範囲に含まれる他の方法で、本発明を具体化することもできる。例えば、本明細書に開示した方法の実施形態は主として、測定音響量が音圧である実施形態に関して説明した。しかしながら、他の実施形態では別の音響量を測定することができることが理解される。例えば、音圧の代わりに、測定位置ｒ_ｉにおけるχ方向の粒子速度ｕ_χ（ｒ_ｉ）を測定することもできる。このような一実施形態では、展開係数ベクトルａについて解く１次方程式系が用意されているとき、式（１）の代わりに、位置ｒ_ｉに対して以下の式を使用することができる。

上式には、展開関数の以下の空間導関数が導入されている。

これは、行列ＡおよびＢにおいて、関数Ψ_ｎ（ｒ_ｉ）が、Γ_χ，ｎ（ｒ_ｉ）、ｎ＝１，２，．．．，Ｎに置き換えられ、式（３）のベクトルｐにおいて、圧力ｐ（ｒ_ｉ）が測定粒子速度ｕ_χ（ｒ_ｉ）に置き換えられたことを意味する。次いで、これに従って、自己相関行列のＡ^ＨＡの含まれる要素を変更しなければならない。

本明細書に記載した方法および装置を使用して、例えば機械、モータ、エンジン、自動車などの車両および／または他の同種の振動物体の音響特性を解析するときに、これらの振動物体などのさまざまな音源／雑音源の音場を再構成することができる。

本明細書に記載した方法の実施形態は、異なるいくつかの要素を含むハードウェアによって実現することができ、かつ／または適当にプログラムされたマイクロプロセッサによって少なくとも部分的に実現することができる。

いくつかの手段を列挙した装置請求項では、これらの手段のいくつかを、同じ１つのハードウェア要素、構成要素またはアイテムによって具体化することができる。単に、ある種の手段が相異なる従属請求項に記載されていること、またはある種の手段が異なる実施形態に記載されていることは、これらの手段の組合せを有利に使用することができないことを意味しない。

本明細書で使用されるとき、用語「含む（ｃｏｍｐｒｉｓｅｓ／ｃｏｍｐｒｉｓｉｎｇ）」は、明示された特徴、要素、ステップまたは構成要素の存在を指定すると解釈されるが、この用語は、１つまたは複数の他の特徴、要素、ステップ、構成要素またはこれらの集団の存在または追加を排除しないことを強調しておく。

Claims (30)

1. 少なくとも１つの音源が発生させた音場を再構成する方法であって、
   一組の測定位置のそれぞれで測定した少なくとも１つの第１の音響量の一組の測定値を受け取ること、
   前記測定位置のうちの第１の測定位置における少なくとも一組の平面波と第２の位置における前記少なくとも一組の平面波との相関関係をそれぞれが指示する一組の相関関数のそれぞれを計算すること、
   少なくとも計算した前記相関関数と前記一組の測定値とから、標的位置における再構成された音場を指示する第２の音響量を計算すること
   を含み、さらに、
   前記相関関数に対するそれぞれの寄与を補間する一組の補間関数の表現を記憶することを含み、それぞれの寄与が、入力パラメータが２つ以下の関数であり、一組のそれぞれの相関関数を計算することが、記憶した前記表現から得た値の１次結合としてそれぞれの相関関数を計算することを含む
   方法。
2. 前記一組の平面波が、一組の平面伝搬波および平面エバネッセント波である、請求項１に記載の方法。
3. 前記一組の平面波が、平面伝搬波および平面エバネッセント波の連続スペクトルである、請求項２に記載の方法。
4. 前記第２の位置が、前記標的位置と前記測定位置の中から選択された、請求項１乃至３のいずれか１項に記載の方法。
5. それぞれの補間関数の前記表現がルックアップ・テーブルを含み、それぞれのルックアップ・テーブルが、予め計算されたそれぞれの寄与の値を含む、請求項１乃至４のいずれか１項に記載の方法。
6. それぞれのルックアップ・テーブルが、２つの入力パラメータによって索引が付けられた関数値の２次元配列を指示し、それぞれの関数値が、前記２つの入力値によって定義される入力点における２次元配列の値を指示する、請求項５に記載の方法。
7. 前記少なくとも一組の平面波のうちのそれぞれの平面波が、所定の仮想源平面に対して定義された伝搬方向を有する、請求項１乃至６のいずれか１項に記載の方法。
8. 前記少なくとも１つの音源を含む領域が、前記少なくとも１つの音源を源のない領域から分離する音源面によって定義され、前記仮想源平面が、前記源のない領域の方を向いた第１の側を有し、前記平面波の前記伝搬方向が、仮想源平面から前記源のない領域の方を向いた、請求項７に記載の方法。
9. 前記仮想源平面が、前記音源面の前記少なくとも１つの音源に近い方の側にあるように定義された、請求項７または８に記載の方法。
10. 前記仮想源平面の前記音源面からの距離が、前記一組の測定位置のうちの隣接する最も近い測定位置間の平均距離よりも大きくなるように、前記仮想源平面が定義された、請求項９に記載の方法。
11. 前記相関関数に対する前記それぞれの寄与がそれぞれ、２つの入力パラメータによってパラメータ表示された積分を含む２次元積分関数である、請求項７乃至１０のいずれか１項に記載の方法。
12. 前記２つの入力パラメータのうちの第１のパラメータが、前記仮想源平面内へ射影された、前記第１の測定位置と前記第２の位置の間の射影された距離を指示し、前記射影された距離が、波数を指示する係数によってスケーリングされる、請求項１１に記載の方法。
13. 前記２つの入力パラメータのうちの第２のパラメータが、波数を指示する係数によってスケーリングされた、前記仮想源平面から前記第１の測定位置までの距離と前記仮想源平面から前記第２の位置までのそれぞれの距離の１次結合を指示するパラメータである、請求項１１または１２に記載の方法。
14. それぞれの平面波が、所定の重み付けスキームに基づくそれぞれの重み係数によって重み付けされる、請求項７乃至１３のいずれか１項に記載の方法。
15. 前記所定の重み付けスキームが、一様重み付けと方向性重み付けの中から選択され、平面波の前記方向重み付けが、前記平面波の伝搬方向の関数である、請求項１４に記載の方法。
16. 一組の補間関数の前記表現が、それぞれの重み付けスキームに関する６つのそれぞれの２次元関数を補間する６つの補間関数の表現を含む、請求項１４または１５に記載の方法。
17. 前記一組の測定位置が、１つの測定平面または平行な複数の測定平面内に配置され、前記測定平面が、１つまたは２つの前記仮想源平面と平行である、請求項７乃至１６のいずれか１項に記載の方法。
18. 前記相関関数に対する前記それぞれの寄与がそれぞれ、第１および第２の入力パラメータによってパラメータ表示された２次元関数であり、測定平面が仮想源平面と平行であり、前記標的位置が、前記測定平面と平行な再構成平面として配置され、一組の補間関数の表現を記憶することが、前記第２の入力パラメータのそれぞれの値に対するルックアップ・テーブルを作成することを含み、それぞれのルックアップ・テーブルが、前記第１の入力パラメータによって索引が付けられた前記２次元関数の関数値の１次元配列を指示し、前記配列が、前記第２の入力パラメータの対応値および前記第１の入力パラメータのそれぞれの値における前記２次元関数の関数値を含む、請求項１７に記載の方法。
19. 前記第１の入力パラメータが、相関関数を計算する前記２つの位置の射影間の距離を指示し、前記射影が、前記仮想源平面上の射影であり、前記距離に波数が掛けられる、請求項１８に記載の方法。
20. 前記第２の音響量が、前記仮想源平面に平行な方向の成分を有する粒子速度である、請求項７乃至１９のいずれか１項に記載の方法。
21. それぞれの一組の平面波にそれぞれが対応する２つ以上の仮想源平面を定義することを含み、前記相関関数がそれぞれ、前記測定位置のうちの第１の測定位置における前記一組の平面波と第２の位置における前記一組の平面波との相関関係を指示する、請求項７乃至２０のいずれか１項に記載の方法。
22. 前記一組の相関関数を使用して、計算された第１の音響量の推定誤差を指示する量を計算することを含む、請求項１乃至２１のいずれか１項に記載の方法。
23. 前記測定位置が、不規則平面アレイとして配置された、請求項１乃至２２のいずれか１項に記載の方法。
24. 前記第１の音響量が、音圧、音圧勾配および粒子速度の中から選択される、請求項１乃至２３のいずれか１項に記載の方法。
25. 前記第２の音響量が、音圧、粒子速度および音の強さの中から選択される、請求項１乃至２４のいずれか１項に記載の方法。
26. 少なくとも１つの補間関数の前記表現が、前記それぞれの寄与の予め計算された関数値を含むルックアップ・テーブルを含み、前記寄与が、少なくとも一対の入力値に対する特異性を含む関数であり、記憶した前記表現から得た値の１次結合として相関関数を計算することが、
    前記特異性を補償するために、前記寄与に対する前記それぞれの入力値に依存するそれぞれのスケーリング係数を掛けたいくつかの前記記憶した関数値を得ること、
    前記スケーリングされた関数値間の補間を実行して、スケーリングされた補間値を得、前記補間値を対応するスケーリング係数で割ること
    を含む、請求項１乃至２５のいずれか１項に記載の方法。
27. 音場を再構成する処理装置であって、一組の測定位置のそれぞれで測定した少なくとも第１の音響量の一組の測定値を受け取るインタフェースと、
    前記測定位置のうちの第１の測定位置における少なくとも一組の平面波と第２の位置における前記少なくとも一組の平面波との相関関係をそれぞれが指示する一組の相関関数のそれぞれを計算し、
    少なくとも計算した前記相関関数と前記一組の測定値とから、標的位置における再構成された音場を指示する第２の音響量を計算する
    ように構成された処理ユニットとを備え、
    前記処理ユニットが、前記相関関数に対するそれぞれの寄与を補間する一組の補間関数の表現を記憶する記憶媒体を備え、それぞれの寄与が、入力パラメータが２つ以下の関数であり、前記処理ユニットが、前記一組のそれぞれの相関関数を、記憶した前記表現から得た値の１次結合として計算するように適合された
    処理装置。
28. 音場を再構成するシステムであって、請求項２７に記載の処理装置と、一組の測定位置で第１の音響量を測定する一組の変換器であり、測定した前記第１の音響量を前記処理装置へ送るように、前記装置に通信接続で接続可能な一組の変換器とを備えるシステム。
29. 車両の雑音源が発生させた音場を再構成するために請求項２７に記載の装置を使用すること。
30. プログラム・コード手段を含むコンピュータ・プログラムであって、前記プログラム・コード手段がデータ処理システム上で実行されたときに、前記データ処理システムが、請求項１乃至２６のいずれか１項に記載の方法の前記ステップを実行するように、前記プログラム・コード手段が適合されたコンピュータ・プログラム。