説明

音場の再構成

音場を再構成する方法を開示する。この方法は、一組の測定位置で測定した第1の音響量の測定値を受け取ること、および標的位置に対する第2の音響量を平面波の重ね合せによって計算することを含む。この方法は、それぞれの関数の一組の補間表現を記憶することを含み、関数はそれぞれ、入力パラメータが2つ以下の関数であり、計算することが、一組の相関関数のそれぞれを計算することを含み、それぞれの相関関数が、前記測定位置のうちの第1の測定位置における平面波と第2の位置における平面波との相関関係を、前記一組の補間表現から得た値の1次結合として指示する。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は音場(acoustic field)の再構成(reconstruction)に関する。
【背景技術】
【0002】
近距離場音響ホログラフィ(NAH:Near−field Acoustical Holography)は、音波放射の3D視覚化、および音源(sound source)に近いある表面を横切る測定に基づく雑音源の精確な位置同定のための非常に有用なツールである。エバネッセント波(evanescent wave)成分をも再構成することができる近距離場音響ホログラフィの能力は、非常に高い空間分解能を保証する。
【0003】
知られているある近距離場音響ホログラフィ法は、分離可能な座標系内のある水準面(level surface)を横切る規則格子(regular grid)測定に基づき、この測定は、空間離散フーリエ変換(DFT:Discrete Fourier Transform)によってNAH計算を実行することを可能にする。例えば、E.G.Williams、J.D.MaynardおよびE.J.Skudrzyk、「Sound source reconstruction using a microphone array」、J.Acoust.Soc.Am.68、340〜344、(1980)を参照されたい。DFTを使用するため処理は非常に高速だが、DFTを使用することの副作用には、測定エリアが音圧の高いエリアを完全にカバーしない場合の重大な空間ウィンドウイング効果(spatial windowing effect)が含まれる。場合によっては、測定エリアに関するこの要件を満たすことができず、多くの場合には、必要なサイズが極端に大きくなる。
【0004】
空間ウィンドウイング効果を低減させ、同時に、複雑さおよび計算要求量の増大と引替えにDFT空間処理をなおも維持する一連の技法が提案されている。例えば、J.Hald、「Reduction of spatial windowing effects in acoustical holography」、Proceedings of Inter−Noise、1994を参照されたい。一般に、最初に反復的手順を使用して測定音圧を測定エリアの外側へ外挿し、続いて、その拡張されたデータ・ウィンドウに、DFTに基づくホログラフィ法を適用する。
【0005】
空間DFTの使用を回避し、必要な測定エリアを低減させることを追求する他の方法が提案されている。
【0006】
そのような1つの方法が、球面波関数に関する音場の局所モデルを使用するヘルムホルツ方程式最小2乗(HELS:Helmholtz’Equation Least Squares)法である。例えば、米国特許第6,615,143号、Z.WangおよびS.F.Wu、「Helmholtz equation−least−squares method for reconstructing the acoustic pressure field」、J.Acoust.Soc.Am、102(4)、2020〜2032(1997)、またはS.F.Wu、「On reconstruction of acoustic fields using the Helmholtz equation−least−squares method」、J.Acoust.Soc.Am、107、2511〜2522(2000)を参照されたい。しかしながら、共通の原点を有する球面波関数だけを使用して音場を表現するため、音源面(source surface)が球面でなく、音源面の中心が同じ原点にない場合には、音源面における音場の再構成に誤差が導入される。この先行技術の方法の他の欠点は、この推定問題が劣決定(under−determined)でない必要があることであり、このことは、測定位置の数が、場の表現(field representation)において使用される球面波関数の数に等しいか、またはそれよりも多くなければならないことを意味する。これにより、十分に正確なモデルを得るために多数の測定位置が必要となることがある。第3の欠点は、モデル領域においてより強く減衰する関数が、同じ領域においてより低い振幅にスケーリング(scaling)されるような方式では、波動関数のスケーリングが適用されないことである。このようなスケーリングがなされないため、チホノフ正則化(Tikhonov regularization)のような伝統的な正則化法が適切に機能しない。その代わりに、上記の先行技術の方法は、正則化されていない最小2乗解と組み合わせた球面波展開の最適な打切り(truncation)のための計算コストが高い反復的な検索を使用する。
【0007】
以前に提案された別の方法が、R.SteinerおよびJ.Hald、「Near−field Acoustical Holography without the errors and limitations caused by the use of spatial DFT」、Intern.J.Acoust.Vib.6、83〜89(2001)において提唱されている統計的最適化近距離場音響ホログラフィ(SONAH:Statistically Optimized Near−field Acoustical Holography)法である。この先行技術の方法は平面波関数に基づき、全ての伝搬波および重み付けした一組のエバネッセント波が最適な平均正確さで射影されるような方法で定義された伝達行列(transfer matrix)を使用することによって、測定表面の近くのマッピング表面上の音響量を計算する。この伝達行列は、1)全ての測定位置対間の自己および相互相関の自己相関行列、ならびに2)それぞれの測定点から全ての再構成位置への相互相関の相互相関行列から得られる。この自己および相互相関は、平面伝搬波および平面エバネッセント波の領域(domain)にある。
【0008】
特定の平面測定ジオメトリ(geometry)へのこの枠組みの適用が、J.Hald、「Patch near−field acoustical holography using a new statistically optimal method」、Proceedings of Inter−Noise 2003、およびJ.Hald、「Patch holography in cabin environments using a two−layer handheld array and an extended SONAH algorithm」、Proceedings of Euronoise 2006に記載されている。
【0009】
これらの先行技術の文献は、測定平面に垂直な圧力および粒子速度成分の計算を開示している。しかしながら、これらの先行技術の文献に開示されている相関式の数値評価は計算効率が悪いため、この方法は計算コストが高い。
【先行技術文献】
【特許文献】
【0010】
【特許文献1】米国特許第6,615,143号
【非特許文献】
【0011】
【非特許文献1】E.G.Williams、J.D.MaynardおよびE.J.Skudrzyk、「Sound source reconstruction using a microphone array」、J.Acoust.Soc.Am.68、340〜344、(1980)
【非特許文献2】J.Hald、「Reduction of spatial windowing effects in acoustical holography」、Proceedings of Inter−Noise、1994
【非特許文献3】Z.WangおよびS.F.Wu、「Helmholtz equation−least−squares method for reconstructing the acoustic pressure field」、J.Acoust.Soc.Am、102(4)、2020〜2032(1997)
【非特許文献4】S.F.Wu、「On reconstruction of acoustic fields using the Helmholtz equation−least−squares method」、J.Acoust.Soc.Am、107、2511〜2522(2000)
【非特許文献5】R.SteinerおよびJ.Hald、「Near−field Acoustical Holography without the errors and limitations caused by the use of spatial DFT」、Intern.J.Acoust.Vib.6、83〜89(2001)
【非特許文献6】J.Hald、「Patch near−field acoustical holography using a new statistically optimal method」、Proceedings of Inter−Noise 2003
【非特許文献7】J.Hald、「Patch holography in cabin environments using a two−layer handheld array and an extended SONAH algorithm」、Proceedings of Euronoise 2006
【非特許文献8】J.GomesおよびP.C.Hansen、「A study on regularization parameter choice in Near−field Acoustical Holography」、Proceedings of Acoustics’08、2008
【非特許文献9】J.Hald、「Array designs optimized for both low−frequency NAH and high−frequency beamforming」、Proceedings of Inter−Noise 2004
【発明の概要】
【課題を解決するための手段】
【0012】
本明細書では、少なくとも1つの音源が発生させた音場を再構成する方法を開示する。この方法の実施形態は、
− 一組の測定位置のそれぞれで測定した少なくとも1つの第1の音響量の一組の測定値を受け取ること、
− 前記測定位置のうちの第1の測定位置における少なくとも一組の平面伝搬波および平面エバネッセント波と第2の位置における少なくとも一組の平面波との相関関係をそれぞれが指示する一組の相関関数のそれぞれを計算すること、
− 少なくとも計算した相関関数と一組の測定値とから、標的位置における再構成された音場を指示する第2の音響量を計算すること
を含み、さらに、
相関関数に対するそれぞれの寄与(contribution)を補間する一組の補間関数(interpolation function)の表現を記憶することを含み、それぞれの寄与が、入力パラメータが2つ以下の関数であり、一組のそれぞれの相関関数を計算することが、記憶した表現から得た値の1次結合としてそれぞれの相関関数を計算することを含む。
【0013】
本発明の発明者は、本明細書に記載した方法を、一般的なジオメトリおよびいくつかの異なる音響量の計算に、一様でかつ計算効率がよい方式で適用することができることを理解した。多様な測定ジオメトリおよび音響量に対する相関関数を、入力値が2つ以下のそれぞれの補間関数によって少なくとも局所的に近似することができる1次元および/または2次元関数の1次結合として表現することができることも分かった。具体的には、補間関数が2次元以下の関数であるため、計算効率がよく、それにも関わらず計算が正確な補間が達成され、補間関数を記憶するのに必要な記憶容量がわずかですむ。したがって、数値計算の基礎をなす諸式が、補間と組み合わされた近似関数またはテーブルを使用することによる効率的な計算によく適した形を有する。本明細書に記載した方法の実施形態は、先行技術の方法の約10倍の速度で実行することができる。
【0014】
音源は、音響放射を放出し、または反射する任意の物体とすることができる。物体の形状は任意であり、音は、あらゆる種類の音、例えば雑音、可聴音、不可聴音(超音波、超低周波音など)など、あるいはこれらの音の組合せとすることができる。音源(1つまたは複数)を源のない領域(source−free region)から分離する音源面によって、1つまたは複数の音源を含む領域を定義することができる。音源面は一般に、音源(1つまたは複数)を含む物理的物体の表面とすることができる。
【0015】
この方法の実施形態は、一組の平面伝搬波関数および平面エバネッセント波関数をそれぞれが表現する1つまたは複数の仮想源平面を定義することを含み、これらの平面波はそれぞれ波数を有し、仮想源平面に対してある方向に伝搬しかつ/または減衰する(evanescent)。一般に、仮想源平面に関して、源から、測定位置および再構成位置を含む領域(この領域を測定領域とも呼ぶ)内へ向かう法線を定義することができる。対応する平面伝搬波は一般に、垂直方向から垂直方向に対して90度までの全ての伝搬方向に伝搬することができる。エバネッセント波は、この垂直方向において、すなわち測定領域内へ向かって指数関数的に減衰する。エバネッセント波は、異なる減衰率および異なる接線伝搬方向を有する連続スペクトルを形成することができる。
【0016】
例えば、仮想源平面は、音源を含む物体の表面の一部分に対して、すなわち音源面に対して定義することができる。例えば、音源面から短い距離を置いて離隔した、物体の内部に位置する表面として、仮想源平面を定義することができる。
【0017】
測定位置における音響量は、適当な音響測定装置、例えばマイクロホン、水中マイクロホン、圧力勾配変換器、粒子速度変換器などによって、またはこれらの装置の組合せによって測定することができる。ある実施形態では、この測定が、音響測定装置のアレイ(array)、例えば規則格子または不規則格子(irregular grid)、例えば2次元または3次元格子として配置された一組の音響測定装置によって実行される。測定する第1の音響量は、音圧、音圧勾配、粒子速度および/または他の同種の音響量とすることができる。
【0018】
再構成された音場は、音圧、音の強さ、粒子速度などの第2の音響量の空間分布を指示する適当な任意のデータ・セットによって表現することができる。第1の音響量と第2の音響量は、同じ音響量でもまたは異なる音響量でもよい。このデータ・セットを、任意の表面、例えば解析対象の音/雑音放出物体の表面またはその物体の近くの表面の実際の表面ジオメトリ上に、マップ、例えば第2の音響量のマップとして直接に表現することができる。
【0019】
硬い表面上の粒子速度は、その表面自体の実際の振動に密接に対応するため、構造モデルとの相関に対してそれらの結果を直接に使用することができる。ある実施形態では、音場が、共形マップ(conformal map)によって表現される。ある表面、一般に音源面または音源面の近くの表面上において音場パラメータを再構成し、マップすることができる。あるいは、またはこれに加えて、再構成された音場の特性の他の表現、例えば波動関数表現の3D妥当性領域(3D region of validity)内の単一位置におけるスペクトルを生成することもできる。
【0020】
本明細書に記載した方法の実施形態では、計算点における圧力または粒子速度の再構成が、1)全ての測定位置対間の自己および相互相関の自己相関行列の計算、ならびに2)それぞれの測定点から標的位置への相互相関の相互相関ベクトルの計算を含む。この自己相関行列は、全ての標的位置に対して再使用することができる。自己および相互相関は、場表現において使用される平面伝搬波関数および平面エバネッセント波関数の領域内で定義される。したがって、第2の位置は、標的位置と測定位置の中から選択することができる。具体的には、この計算は、それぞれの測定点間の一組の平面波の相関関数の自己および相互相関行列の計算、ならびにそれぞれの測定点と標的点との間の一組の平面波の相関関数の相互相関ベクトルの計算を含むことができる。
【0021】
それぞれの寄与を補間する補間関数という用語は、それらの寄与を局所的に近似する任意の関数または任意の一組の関数を指すことが意図されている。補間関数の表現が、それぞれの関数の予め計算された局所近似をそれぞれが含む一組のルックアップ・テーブルを含むときには、特に効率的な計算が提供される。それぞれのルックアップ・テーブルは、2つの入力パラメータによって索引が付けられた関数値の2次元配列を指示することができ、それぞれの関数値は、それらの2つの入力値によって定義される入力点における2次元配列の値を指示することができる。これらのそれぞれの関数はそれぞれ、2つの入力パラメータによってパラメータ表示された積分を含む2次元積分関数とすることができる。したがって、テーブル化された関数値間の補間、例えば3次補間(cubic interpolation)によって、別の点における関数値を効率的に得ることができる。このルックアップ・テーブルは低次元のテーブル、すなわち2以下の次元を有するテーブルであるため、この補間法は正確でかつ計算効率がよく、このルックアップ・テーブルは記憶領域をわずかしか必要とせず、少ない計算資源を用いて作成することができる。
【0022】
2つの入力パラメータのうちの第1のパラメータは、波数を指示する係数(factor)によってスケーリングされた、第1の測定位置と仮想源平面に平行な平面内の第2の位置の間の射影された距離を指示することができる。2つの入力パラメータのうちの第2のパラメータは、波数を指示する係数によってスケーリングされた、仮想源平面から第1の測定位置までの距離と仮想源平面から第2の位置までの距離の1次結合(一般に和または差)を指示することができる。ある実施形態では、このパラメータ表示が、相関関数を、予め計算された一組の2次元関数によって表現することを可能にする。第2の入力パラメータは、仮想源平面に対する法線(通常はz軸)上へ射影され、波数を指示する係数によってスケーリングされた、第1の測定位置と第2の位置の間の射影された距離を指示することができる。
【0023】
ある実施形態では、それぞれの平面波波動関数が、所定の重み付けスキーム(weighting scheme)に基づくそれぞれの重み係数(weighting factor)によってスケーリング/重み付けされる。例えば、この所定の重み付けスキームは、一様重み付け(uniform weighting)と方向性重み付け(directional weithting)の中から選択することができる。一様重み付けスキームでは、それぞれの平面波が同じ定数係数によってスケーリングされる。平面波の方向性重み付けでは、それぞれの平面波が、平面波の伝搬方向に依存したそれぞれの重み係数によって、例えば、源平面に平行なより大きな成分を有する伝搬方向を有する平面波により高い重み係数を提供することによってスケーリングされる。このような波は、仮想源平面の表面法線(surface normal)に対して平行でないため、離軸波(off−axis wave)とも呼ばれる。方向性重み付けは、再構成法の正確さを向上させることが分かっている。
【0024】
ある実施形態によれば、所与の重み付けスキームに関して、さまざまな測定ジオメトリおよび再構成ジオメトリに対する音圧と粒子速度の両方を、6つの2次元関数の一組のテーブル表現に基づいて計算することができることが分かった。それによって補間された関数値を計算しなければならない関数パラメータだけがジオメトリとともに変化する。
【0025】
ある実施形態では、仮想源平面が、測定領域とも呼ぶ、測定位置および標的位置を含む第1の空間領域と、源領域とも呼ぶ、第1の空間領域とは異なる(第1の空間領域とは反対の)第2の空間領域とを定義するように選択され、第2の空間領域は一般に音源の大部分を含む。すなわち、標的位置および測定位置は測定領域に位置し、音場は、測定領域において再構成される。仮想源平面は一般に、源平面の内側に短い距離を置いて定義されるが、この再構成法の目的上、仮想源平面によって定義される測定領域は源を含まないと仮定する。
【0026】
状況によっては、単一の仮想源平面によって源のない空間領域を定義することができない、またはそうすることが実際的でないことがある。したがって、ある実施形態では、この方法が、それぞれの平面波がそれぞれの仮想源平面に対してある方向に伝搬し、減衰する2組以上の平面波に基づく。例えば、それぞれの仮想源平面対間に源のない領域が位置する間隔を置いて配置された1つまたは複数の平行仮想源平面対として、仮想源平面を定義することができる。本発明の発明者は、本発明の方法およびテーブル化された同じ一組の補間関数を、この状況にも適用することができ、それによって効率的で柔軟な再構成法を提供することができることを理解した。
【0027】
この一組の測定位置を、1つまたは複数の測定平面内、例えば単一の平面または平行な2つ以上の平面内に配置することができる。それぞれの平面において、測定位置は、規則格子状にまたは不規則パターンとして、あるいは他の適当な方式で配置することができる。さらに、本明細書に記載した方法を、非平面測定ジオメトリ、すなわち測定位置が、1つまたは複数の平行平面内に配置されず、例えば曲面上に配置される配置に対して適用することもできる。
【0028】
ある実施形態では、全ての仮想源平面に平行なそれぞれの平面内に測定点および計算点を配置すると、効率をさらに増大させることができる。このような構成では、平面波のそれぞれの周波数に対して、相関関数に対する寄与の入力パラメータのうちの1つの入力パラメータの値が少数だけ生じ、したがって、それぞれのテーブルが前記入力パラメータのそれぞれの値に対応する1Dルックアップ・テーブルの作成を作成することが魅力的になる。したがって、ある種のジオメトリに関しては、この方法を1次元補間関数に基づく方法とすることさえでき、したがってこの方法を特に高い計算効率で実行することができる。
【0029】
第2の音響量が、仮想源平面に平行な方向の成分を有する粒子速度であるときには、曲面の音源面を有する物体の改良された解析が容易になる。具体的には、この実施形態は、(例えば音源面と共形の)曲面の再構成表面を横切って分布する標的位置における改良された音場の再構成を可能にする。再構成法のこの実施形態は、仮想源平面に平行な粒子速度の成分を提供するため、再構成表面に垂直な粒子速度(および/または強さ)の成分を、曲面の再構成表面を横切って決定することができる。
【0030】
本明細書に開示した方法の実施形態を、標的点における音場パラメータの線形予測を一組の測定量に基づいて実行するものとみなすことができ、予測係数は、最小2乗法の意味において、全ての平面伝搬波および重み付けされた一組のエバネッセント波に対する最良推定値を提供するように決定される。最小2乗フィット(least−squares fit)は、優決定(over−determined)の一組の1次方程式を近似的に解く一組の複素係数を決定する数値的プロセスとみなすことができる。
【0031】
次いで、計算点において、これらの平面波波動関数の値が、平均して、どれくらいよく予測されるのかを知ることが望ましい。したがって、ある実施形態では、再構成法が、計算された音響量の推定誤差を指示する量を計算することを含む。以前に説明した自己および相互相関を使用して、すなわち計算資源をほとんど追加せずに、誤差インジケータ(error indicator)を計算することができる。SONAH法の予測の正確さを推定する上記の方法は、所与の測定ジオメトリを有するSONAH予測の妥当性領域を視覚化し、決定するのに非常に有用であることが立証されている。
【0032】
上述の方法および後述する方法の諸特徴は、ソフトウェアまたはファームウェア内において少なくとも部分的に実現することができ、コンピュータ実行可能命令などのプログラム・コード手段を実行することによって、データ処理装置または他の処理手段上で実施することができることに留意されたい。この説明および以下の説明では、用語「処理手段」が、上記の機能を実行するように適切に適合された任意の回路および/または装置を含む。具体的には、用語「処理手段」は、汎用または専用プログラマブル・マイクロプロセッサ、ディジタル信号プロセッサ(DSP)、特定用途向けIC(ASIC)、プログラマブル論理アレイ(PLA)、フィールド・プログラマブル・ゲート・アレイ(FPGA)、専用電子回路など、あるいはこれらの組合せを含む。
【0033】
本発明の実施形態は、上述の方法および後述する方法、ならびに最初に述べた方法に関して説明した利益および利点のうちの1つまたは複数の利益および利点をそれぞれがもたらし、最初に述べた方法に関して説明した実施形態および/または従属請求項に開示された実施形態に対応する1つまたは複数の実施形態をそれぞれが有するシステム、装置および製品手段を含む、さまざまな方法で実現することができる。
【0034】
具体的には、音場を再構成する処理装置の実施形態は、一組の測定位置のそれぞれで測定した少なくとも1つの第1の音響量の一組の測定値を受け取るインタフェースと、
− 前記測定位置のうちの第1の測定位置における少なくとも一組の平面波と第2の位置における少なくとも一組の平面波との相関関係をそれぞれが指示する一組の相関関数のそれぞれを計算し、
− 少なくとも計算した相関関数と一組の測定値とから、標的位置における再構成された音場を指示する第2の音響量を計算する
ように構成された処理ユニットとを備え、
処理ユニットは、相関関数に対するそれぞれの寄与を補間する一組の補間関数の表現を記憶する記憶媒体を備え、それぞれの寄与は、入力パラメータが2つ以下の関数であり、処理ユニットは、一組のそれぞれの相関関数を、記憶した表現から得た値の1次結合として計算するように適合されている。
【0035】
記憶された表現は、再構成法を実行するコンピュータ・プログラムのソフトウェア・コンパイル時間に(例えばコンピュータ・プログラムの一部としてまたは別個の1つもしくは複数のファイルとして)生成することができ、あるいは、コンパイルの結果生成されたプログラム・コードによって実行時に生成することができ、あるいは、これらを組み合わせて生成することができることが理解される。例えば、この表現は、コンピュータ・プログラムまたはインスレーション・プログラムによって生成し、そのコンピュータ・プログラムの最初の実行/インストレーションの間にコンピュータ上に記憶することができる。
【0036】
音場を再構成するシステムは、上に開示し、後に開示する装置と、一組の測定位置で第1の音響量を測定する一組の変換器であり、測定した値を装置へ送るように、装置に通信接続で接続可能な一組の変換器とを備える。
【0037】
コンピュータ・プログラムは、プログラム・コード手段を含むことができ、このプログラム・コード手段は、このプログラム・コード手段がデータ処理システム上で実行されたときに、データ処理システムが、上に開示し、後に開示する方法のステップを実行するように適合される。このコンピュータ・プログラムは、記憶手段上に記憶することができ、またはデータ信号として具体化することができる。記憶手段は、RAM、ROM、EPROM、EEPROM、フラッシュ・メモリ、磁気または光記憶装置(CD ROM、DVD、ハード・ディスクおよび/または他の同種の装置など)など、データを記憶する適当な回路または装置を備えることができる。
【0038】
上記の態様およびその他の態様は、図面を参照して以下で説明する実施形態によって明白となる。
【図面の簡単な説明】
【0039】
【図1】音場を再構成するシステムの概略ブロック図である。
【図2】再構成された音場を計算する方法の流れ図である。
【図3】仮想源平面が1つだけの状況、2つの平行な仮想源平面が源のない均質な領域を限定している状況、および箱形の均質な源のない領域の6つ全ての面に仮想源平面がある状況の測定ジオメトリを示す図である。
【図4】(k,k)領域および伝搬角領域における平面波波動関数の密度を示す図である。
【図5】アレイ平面内の音圧予測の固有相対誤差水準E(r)を、2つの不規則平面アレイについて示す図である。
【図6】単層および2重層8×8要素アレイのxz対称平面における相対誤差水準を示す図である。左列:単層アレイの固有誤差水準E(r)。中央列:単層アレイによって推定される単極子圧力の相対誤差。右列:2重層アレイの固有誤差水準E(r)。
【図7a】(a)8×8要素平面アレイのアレイ平面に垂直な対称軸上の固有相対誤差水準Eを示す図である。
【図7b】8×8要素平面アレイのアレイ平面に垂直な対称軸上の振幅利得Aを示す図である。
【図8】8×8要素平面アレイのアレイ平面に垂直な対称軸上の固有相対誤差水準E(r)および振幅利得A(r)を、低周波数および高周波数について示す図である。これらの曲線は、再構成法で使用する異なるダイナミック・レンジDを表す。
【図9】8×8要素平面アレイのアレイ平面に垂直な対称軸上の固有相対誤差水準E(r)および振幅利得A(r)を、低周波数および高周波数について示す図である。これらの曲線は、再構成法で使用する異なる仮想源平面距離d=−zを表し、エバネッセント波の実際の内容を定義する仮想源距離は一定の6cmに維持される。
【図10】2重層8×8アレイの中心平面内の単一の単極子点源から3Dインテンシティベクトルを計算する一実施形態において、2つの面だけからではなく6つの全ての面から平面波波動関数を追加することによって得られる正確さの向上を示す図である。点源はアレイ・エリアから6cm外側にある。
【発明を実施するための形態】
【0040】
全図面を通じ、可能な場合には常に、等しい参照符号は、等しい要素、特徴または成分、あるいは対応する要素、特徴または成分を指す。
【0041】
図1は、音場を再構成するシステムの概略ブロック図を示す。このシステムは、一組の音響受信器(aoustic receiver)108と、音響受信器に接続された解析ユニット(analysing unit)103とを備える。
【0042】
以下では、音響受信器108を変換器(transducer)とも呼ぶ。のみならず、音響受信器は、マイクロホン、水中マイクロホン、あるいは音圧、音圧勾配、粒子速度または他の線形量などの音響特性を測定する他の適当な装置とすることができることが理解される。図1の実施形態では、変換器108が、変換器108が規則格子、例えば1、2または3次元格子として配置された変換器のアレイ102として実現される。規則格子として配置すると効率的な数値実現が可能だが、他のジオメトリを使用することもできる。ある実施形態では、単一の変換器を使用して、異なる点における測定を逐次的に実施することもできる。このような逐次測定は例えば、音源が静止している状況および1つまたは複数の基準信号に対する位相を測定することができる状況において使用することができる。変換器の数およびアレイのジオメトリ、例えば変換器間の間隔は、解析する物体のサイズおよび幾何学的複雑さ、関心の周波数範囲ならびに所望の空間分解能に応じて選択することができる。
【0043】
変換器アレイ102は、変換器108が、測定した信号を、例えば有線または無線信号接続を介して解析ユニットへ送ることができるような態様で、解析ユニット103に接続される。
【0044】
解析ユニット103は、測定された信号を変換器アレイ102から受け取り、これを処理するインタフェース回路104と、インタフェース回路104とデータ通信する処理ユニット105と、記憶媒体112と、処理ユニット105とデータ通信する出力ユニット106とを含む。図1には単一のユニットとして示されているが、解析ユニット103を、2つの別個の装置、例えば収集フロントエンド装置とコンピュータとに物理的に分割し、または3つ以上の装置に分割することもできることが理解される。同様に、解析ユニットの異なる下位ブロックに関して記載した機能を、代替のまたは追加の機能ユニット/モジュールまたはハードウェア・ユニット/モジュールに分割することができることも理解される。
【0045】
インタフェース回路は、変換器108からの出力信号を受け取り、受け取った信号を、処理ユニット105による後続の解析のために処理するのに適した信号処理回路を含む。インタフェース回路104は、以下の構成要素のうちの1つまたは複数の構成要素を備えることができる:受け取った信号を増幅する1つまたは複数の前置増幅器、受け取った信号をディジタル信号に変換する1つまたは複数のアナログ−ディジタル(A/D)変換器、1つまたは複数のフィルタ、例えば帯域幅フィルタ、および/あるいは他の同種の構成要素。インタフェース回路は例えば、出力データとして、変換器ごとの振幅および位相を、周波数の関数として提供することができる。ある実施形態では、このインタフェース・ユニットが同時データ収集を実行し、次いで、処理ユニット105が、周波数領域へのデータ変換を含む残りの全ての処理を、一般にFFTを使用して実行することができる。
【0046】
処理ユニット105は、適当にプログラムされたマイクロプロセッサ、コンピュータの中央処理ユニット、あるいはインタフェース・ユニット104から受け取った信号を処理する他の適当な装置、例えばASIC、DSPおよび/または他の同種の装置とすることができる。この処理ユニットは、インタフェース回路104を介して受け取った変換器信号を処理して、再構成された音場を本明細書に記載されたとおりに計算するように適合される。
【0047】
記憶媒体112は、RAM、ROM、EPROM、EEPROM、フラッシュ・メモリ、磁気または光記憶装置(CD ROM、DVD、ハード・ディスクおよび/または他の同種の装置など)など、一組の補間関数の表現を指示するデータを記憶する適当な回路または装置を備えることができる。図1には、処理ユニットとは別個だが、処理ユニットに通信接続された記憶媒体が示されている。しかしながら、記憶媒体112は、処理ユニット105の一部として、例えば内部記憶装置として具体化することもできる。
【0048】
出力ユニット106は、表示装置、あるいは再構成された音場の視覚表現を提供する他の適当な装置または回路、例えば印刷された表現を提供するプリンタおよび/またはプリンタ・インタフェースを備えることができる。あるいは、またはこれに加えて、出力ユニット106は、RAM、ROM、EPROM、EEPROM、フラッシュ・メモリ、磁気または光記憶装置(CD ROM、DVD、ハード・ディスクなど)、有線または無線データ通信インタフェース、例えばコンピュータまたは電気通信ネットワーク(LAN、ワイド・エリア・ネットワーク、インターネットなど)とのインタフェース、および/あるいは他の同種の回路または装置など、再構成された音場を指示するデータを伝達し、かつ/または記憶する適当な回路または装置を備えることができる。
【0049】
解析ユニット103は、適当にプログラムされたコンピュータ、例えば適当な信号収集ボードまたは回路を含むPCとして実現することができる。
【0050】
このシステムはさらに、解析ユニット103に接続された適当な位置検出装置110であって、変換器108の位置、例えば基準座標系に対する変換器108の位置を検出する位置検出装置110を備えることができる。この位置検出装置は、変換器アレイを所定の位置および向きに取り付けることができ、かつ/あるいはアレイの位置および向きを手動で測定し、または例えば適当なセンサによって自動的に測定することができる、試験台、フレームまたは他の構造を含むことができる。あるいは、変換器アレイ102が、位置検出装置または位置検出装置の少なくとも一部分を含んでもよい。例えば、変換器アレイが、磁場位置検出システムのコイル、ジャイロスコープ、光学的位置検出システムおよび/または他の同種のシステムを含むことができる。
【0051】
操作時には、音響放射を放出する音源111であり、音場を再構成する対象である音源111を含む物体100の表面の近くに、変換器アレイ102を配置する。変換器の数、アレイのジオメトリ、例えば変換器間の間隔、およびアレイと音源面との間の距離は、解析する物体のサイズおよび幾何学的複雑さ、関心の周波数範囲ならびに所望の空間分解能に応じて選択することができる。例えば、自動車エンジンの音場の再構成に対しては、要素間間隔が3cmの8×8行列のマイクロホン・アレイを使用することができる。しかしながら、他のタイプのアレイおよび他のサイズのアレイを使用することもできることが理解される。
【0052】
例えば位置検出装置110によってアレイ102の位置が決定され、解析ユニット103に供給される。アレイ102の変換器108が、決められた位置において音圧または他の適当な音響量を測定し、その結果得られた変換器信号が解析ユニット103へ送られる。
【0053】
例えば、変換器アレイを、位置検出装置が組み込まれた手持ち式アレイとすることができ、こうすると、物体の周囲に分布した異なるアクセス可能位置で測定することができ、試験に対する源の特別な準備にかかる時間が短縮される。例えば、試験台上の音源、例えば圧縮機、ポンプまたは他の機械をマップするときには、面全体に障害なくアクセスできるようにするために試験構成を変更するのに必要な時間、例えば燃料管路の配置、電子回路の配線などを変更するのに必要な時間が、排除されることはないにしても、短縮される。他の典型的な用途は、自動車の室内での使用であり、そこでは、3Dアレイ格子を使用して、全ての方向の源を識別することができる。例えば(例えば8×8×2センサを含む)2層アレイを使用することができる。
【0054】
解析ユニット103は、この入力データから、再構成された音場を計算し、再構成された音場の表現を記憶し、かつ/または出力する。
【0055】
次に、図2を参照し、図1を続けて参照して、再構成された音場を生成する方法の一実施形態を説明する。
【0056】
本明細書に記載した再構成法を例示するため、最初に、いわゆるSONAH法の概要を示す。
【0057】
この説明の目的上、複素時間調和音圧(complex time−harmonic sound pressure)p(r)は、均質流体によって占められた源のない領域Ω内の一組の位置r、i=1,2,...,Iにおいて測定されたものであると仮定する。測定した音圧に基づいてΩ内における音場を再構成する問題について考える。これを達成するため、最初に、同次波動方程式(homogeneous wave equation)を満たし、Ω内に存在しうる全ての音圧場を十分に高い正確さ(accuracy)で表現することができる、離散した一組の基本波動関数(elementary wave funciton)Ψ、n=1,...,Nを導入する。一般に、これらの基本波動関数は、平面波、円筒波または球面波とすることができ、あるいは、Ωの外側に分布した単極子点源または双極子点源などの基本源からの音場とすることができる。
【0058】
いわゆるHELS法では、再構成の最初のステップが、音場の表現中のそれぞれの基本波動関数の振幅および位相を決定するステップである。このステップは、基本波動モデルが、測定された圧力をできるだけ正確に表現するようにすることによって実行される。
【0059】
【数1】

は複素展開係数(complex expansion coefficient)である。I≧Nであり、測定点が適切に選択された場合には、(1)の最小2乗解によって、展開係数を見つけることができる。そうでない場合、この問題は劣決定であり、(1)の無限の解の中で、最小の2ノルムを有する解を採用することができる。この解を、ヘルムホルツ方程式最小ノルム(HELN:Helmholtz’Equation Least Norm)解と呼ぶ。
【0060】
行列公式(matrix formulation)を得るため、測定位置における波動関数値の行列Bを定義し、
【0061】
【数2】

測定された圧力および展開係数を用いてベクトルpおよびaをそれぞれ定義する。
【0062】
【数3】

これらの定義を使用して、式(1)を次の形に書き直すことができる。
Ba=p (4)
I>Nにおいて、最小2乗(HELS)解は、
a=(BB)−1
であり、I<Nにおいて、最小ノルム(HELN)解は、
a=B(BB−1
であり、記号Hは、エルミート転置(Hermitian transpose)を表す。実際には、これらの解に、正則化を適用しなければならない。他の正則化スキームを使用することもできるが、この説明の目的上、チホノフ正則化だけを考え、その場合、HELSおよびHELN式はそれぞれ、
a=(BB+εI)−1p (5)
および
a=B(BB+εI)−1p (6)
となる。εは正の正則化パラメータ、Iは適当な次元の単位対角2次行列である。実際に、正則化されたこれらの2つの式(5)と(6)は、IおよびNの値とは無関係に全く同じ解を与えることが分かる。しかしながら、それでも、I<Nに関しては、行列BBの次元の方が行列BBの次元よりも小さいため、HELN式(6)を使用した方が計算上有利である。
【0063】
展開係数が分かれば、波動関数の寄与の和を求めることによって、Ω内の任意の点rにおける音圧を計算することができる。
【0064】
【数4】

Tは行列またはベクトルの転置を表す。次いで、オイラーの方程式を適用することによって、音場の表現(7)から、χ方向の粒子速度を推定することができる。
【0065】
【数5】

上式で、ωは角周波数、ρは媒質の密度であり、陰時間変動(implicit time variation)は、ejωtである。続いて使用する他の記号は、媒質中における音の伝搬速度cおよび波数k≡ω/cである。
【0066】
平面波は、波数領域において連続スペクトルを構成し、そのため、上記の公式におけるこれらの使用は、波数領域におけるサンプリングによって得られた波動関数の部分集合を含まなければならない。本明細書に記載した方法の実施形態は、例えばJ.Hald、「Patch near−field acoustical holography using a new statistically optimal method」、Proceedings of Inter−Noise 2003およびJ.Hald、「Patch holography in cabin environments using a two−layer handheld array and an extended SONAH algorithm」、Proceedings of Euronoise 2006に記載されているSONAH法を利用した2重の方法で場推定問題を解くことによって、このサンプリングの問題を解決する。
【0067】
SONAH法によれば、Ω内の点rにおける未知の音圧p(r)は、既知の音圧値の1次結合として得られる。
p(r)=pc(r) (9)
上式で、ベクトルc(r)は複素推定重み(complex estimation weight)を含む。次いで、予測式(9)が、最小2乗法の意味において、全ての展開関数Ψに対する最良平均予測を提供するようにすることによって、ベクトルcを得る。式(2)および(7)を参照すると、これは、以下の一組の1次方程式に対する最小2乗解としてc(r)を見つけなければならないことを意味する。
c(r)=α(r) (10)
この式は、HELS式(4)が優決定であるときは常に劣決定であり、HELS式(4)が劣決定であるときは常に優決定である。この説明の目的上、チホノフ正則化を適用した解を考える。このことは、(10)に対する最小2乗解と最小ノルム解の両方を、以下の最小2乗式に関して表現することができることを意味する。
c(r)=(AA+εI)−1α(r) (11)
上式で、行列Aは行列Bの転置である。
A≡B (12)
定義(7)および(12)から、続いて、要素AAおよびAα(r)が、それぞれ下式(13)および(14)によって与えられる。「」は複素共役化(complex conjugation)を表す。
【0068】
【数6】

行列AAは、基本波動関数(展開関数)の領域内の測定点間の相互相関の行列と見ることができ、ベクトルAα(r)は、測定点と推定位置rとの間の相互相関を含む(例えば、R.SteinerおよびJ.Hald、「Near−field Acoustical Holography without the errors and limitations caused by the use of spatial DFT」、Intern.J.Acoust.Vib.6、83〜89(2001)を参照されたい)。さらに、AAが、負でない固有値を有するエルミート対称行列であり、その次元が、展開関数の数とは無関係に、測定点の数Iに等しいことに気づくこともある。
【0069】
上記の予測公式をSONAH予測公式と呼び、この公式はHELS/HELN式と同じ解を提供する。実際、HELN式(7)および(6)から、このSONAH式は導かれる。
【0070】
【数7】

ここでは、任意の非特異2次行列(non−singular quadratic matrix)について、逆(inverse)の転置は転置の逆に等しいことを利用した。SONAHは展開係数aを明示的には計算しないため、一組の展開係数に対して重み付け/フィルタリングを実行することはできない。このようなフィルタリングは、重みベクトルc(r)の計算に組み込まれている。これを実行するツールは正則化だが、これを適切に機能させるためには、行列AAを定義する基本波動関数の適当なスケーリングが必要である。一般に、より大きな正則化パラメータε値を選択することによって、音源面における再構成において最大の増幅を必要とするまだ抑制されていない基本エバネッセント波を抑制することが望ましい。このことが起こるためには、これらの基本波が、チホノフ正則化のカットオフ・レベルよりもまだ大きいAAの最も小さな特異値(singular value)に対応していなければならない。基本波動関数のスケーリングを通して、再構成においてより大きな増幅を必要とする波動関数がマイクロホン位置を横切ってより小さな振幅を有することを保証することによって、このような振舞を支持することができる。
【0071】
この実施形態では、基本波が、平面伝搬波およびエバネッセント波であり、最初に、図3aに示すようなフリー・フィールド条件(free−field condition)、すなわち、z=z、z<0の源平面107の後ろの1つまたは複数の源が、Ωで示す源のない均質な半空間z≧zに放射する条件において、音場を再構成する場合を考える。この場合には、
【0072】
【数8】

の形の平面波波動関数によって放射音場を表現することができる。上式で、Fは振幅スケーリング関数、z=z、z≦Zは、平面波波動関数を定義するいわゆる仮想源平面101であり、例えばエバネッセント成分はこの仮想源平面に垂直な方向に減衰する。仮想源平面は一般に、実際の音源面から内側に短い距離を置いて位置するように定義される。源平面107は、依然として源のない領域にある仮想源平面に平行な最も近い平面である。したがって源平面は、実際の音源面にちょうど触れる(すなわち実際の音源面に接する)平行な平面である。単一の仮想源平面だけがあるフリー・フィールド条件の場合には、xy平面が一般に仮想源平面に平行である。式(16)では、r≡(x,y,z)およびk≡(k,k,k)であり、kおよびkは実波数、kは、下式で定義される複素数である。
【0073】
【数9】

【0074】
したがって、完全な一組の展開関数は、2つの波数変数kおよびkにおいて連続スペクトルを構成する。仮想源平面z=zは実源平面z=zと一致することがあるが、仮想源平面は一般に実源平面の後ろにある。すなわちz<zである。仮想源平面の座標上の上付き添字「+」は、z軸の正方向に伝搬している波が表現されていることを示す。
【0075】
式(16)では、スケーリング関数F(k)が、仮想源平面z=z上の全ての波動関数の振幅を定義することに留意されたい。
【0076】
上述の離散HELSまたは離散SONAH公式を使用するため、最初に、2D波数領域(k,k)における規則的サンプリング(regular sampling)によって離散部分集合を選択する。
【0077】
【数10】

サンプリング間隔を選択するときには、波数領域における規則的な離散表現が、ちょうど空間DFT処理に基づく伝統的なNAHの場合のようにラップアラウンド誤差(wrap−around error)につながる空間的に周期的な場を表現することを知っておくべきである。サンプリング間隔はしたがって、2πを最小許容周期長(smallest acceptable period length)で割ったものよりも小さくすべきである。同時に、サンプリング・エリアは、全ての伝搬波、および測定領域においてかなりの振幅を有する全てのエバネッセント波をカバーすべきである。これによって一般に、展開関数の数Nが測定位置の数Iよりも大きくなる。したがって、等価(equivalent)のHELS問題は劣決定となり、解(6)は支配的に最小ノルム解となる。
【0078】
ある実施形態では、F(k)=1の場合、波動関数の単位重み付け(unit weighting)を使用することができる。この場合、式(16)によれば、全ての平面波波動関数が、仮想源平面z=zにおける振幅に等しい振幅を有する。したがって、HELS展開係数ベクトルaは単純に、その平面における音場表現の波数スペクトルであり、行列Bを、その平面波スペクトルからマイクロホン位置における音圧値への伝達関数の行列とみなすことができる。Bの要素は、式(16)にマイクロホン座標を挿入することよって与えられる。伝達行列Bによって、エバネッセント波成分(Φ、|k|>k)の減衰が、波数とともに指数関数的に増大することに言及することは興味深い。したがって、全体的な傾向は、これらの波数成分は、行列Bのますます小さな特異値に対応するというものである。この特性は、波数スペクトルaに対する逆解(inverse solution)におけるチホノフ正則化によって、エバネッセント波成分のフィルタリングを実行することができることを保証する。伝搬波のある組合せ、例えば仮想源平面内の測定位置から遠いエリアからの音波放射を表現する組合せも小さな特異値に対応する。このような組合せも、正則化によって抑制される。
【0079】
HELS問題の上述の劣決定性のため、解は、支配的に最小ノルム解となる。正則化によって、解ベクトルaのノルムのさらなる低減が導入される。したがって、この解法は、波数スペクトルにおけるエネルギーを最小化し、結果的に、マイクロホンにおける圧力がうまく表現されているという制約条件下で、仮想源平面内および平行な平面内の圧力分布におけるエネルギーも最小化する。これにより、測定位置の周囲の体積内において良好な音場表現が達成され、測定位置から離れるにつれて増大する過小推定(under−estimation)が見られる。SONAHとHELSは同じ解(15)を提供するため、後に説明する結果によってさらに示すように、SONAHに対する結果も同じである。
【0080】
図4は、k=0におけるk波数軸に沿った式(18)の規則的サンプリングを示す。放射円内において、すなわち
【数11】

に関して、(k,k)における特定の標本は、ベクトル(k,k,k)の方向に伝搬している平面波に対応する。kは式(17)によって与えられる。図4では、示された半径kの半円上に標本位置kを射影することによって、この方向を見つけることができる。(k,k)平面における標本点の密度が一定のとき、平面波の角密度はz方向の近くで最も高く、軸角度(off−axis angle)θが90°に向かって増大すると低下する。したがって、一定の重み関数Fを使用すると、軸方向の近くで伝搬している平面波は、SONAHの予測重みに対する最小2乗解(11)において、軸から遠く外れた波よりも高い重みを受け取る。これは、離軸波により高い重みFを与えることによって補償することができる。(k,k)平面上の無限に小さな面積の弓形(segment)を考えると、それは、1/cos(θ)=k/|k|倍だけ大きな方向半球(directional hemisphere)
【数12】

上のエリアを表す。したがって、
|F(k)|k/|k
が一定になるように、例えば
【数13】

となるようにFを選択することによって、一定の方向パワー・スペクトル密度を得ることができる。上式で、Fは定数である。式(13)、(14)、(16)および(18)は、|F|によってスケーリングされた波数寄与の全体にわたって和を求めることにより、行列AAおよびAα(r)の要素に対する寄与が得られることを示している。このことが、振幅密度に関してではなくパワー密度に関して等化(equalization)を選択することが有利であることがある理由である。放射円の外側にも適用される式(19)の「全方向重み付け(omnidirectional weighting)」は、平均して、一定の重み付けF(k)=1よりもわずかによい正確さを提供することが分かっている。したがって、ある実施形態では、全方向重み付けが使用され、これはさらに、この説明の全体を通じてこの方法を示すために使用される。
【0081】
波数領域における離散表現によって生じるラップアラウンド問題を完全に除去するため、次に、(k,k)平面におけるサンプリング間隔をゼロにする。このことは、式(13)および(14)における求和が(k,k)平面にわたって積分になることを意味する。式(11)は、正則化係数がそれに応じて調整される場合、AAおよびAα(r)の要素に適用される共通のスケーリング係数が、SONAHの結果を変化させないことを示している。したがって、極限積分式(limiting integral expression)を、以下のようにスケーリングすることを選択することができる。
【0082】
【数14】

行列要素の極限値(20)および(21)は、無限の次元を有する行列Aおよびベクトルα(r)に対応する。とにかくもAおよびα(r)は、後続のいくつかの微分において使用され、その場合、これらを、極限N→∞の場合と理解すべきである。
【0083】
式(16)、(17)および(19)を適用することによって、式(20)および(21)の2次元積分を1次元積分にすることができる。最初に、波数領域において極座標(K,Ψ)への変換を実行する(k,k)=(Kcos(Ψ),Ksin(Ψ))。これにより、極角(polar angle)Ψにおける積分を解析的に解くことができる。その結果、式(20)および(21)の積分を、以下の2つの2変数関数の形で表現することができる。
【0084】
【数15】

上式で、Jは、0次のベッセル関数(Bessel function)である。関数Iは、積分式において置換K=ksin(θ)を実施した波数領域における放射円の内側からの寄与を表し、関数IIは、積分式において置換
【数16】

を実施した放射円の外側からの寄与を表す。これらの変数置換は、(22)および(23)の中の積分が、ガウス型の数値積分により適合するように導入したものである。関数IおよびIIを使用して、行列要素に対する式(20)および(21)を、以下のように再公式化することができる。
【0085】
【数17】

ここで、r≡(x,y,z)およびr≡(x,y,z)は2つのマイクロホン位置であり、圧力推定位置はr≡(x,y,z)であり、
【数18】

は、xy平面上(仮想源平面)上におけるマイクロホン位置iとjの射影間の距離、
【数19】

は、xy平面上におけるマイクロホン位置iと圧力計算位置の射影間の距離である。こうすると、式(9)および(11)を使用して、位置rにおける圧力を推定することができる。式(11)の正則化パラメータεを決定するさまざまな方法が、J.GomesおよびP.C.Hansen、「A study on regularization parameter choice in Near−field Acoustical Holography」、Proceedings of Acoustics’08、2008において検討され、SONAHに関しては、L曲線法(L−curve method)が全体的に最も安定した結果を与えることが分かった。これよりも初期の非相関(ncorrelated)測定誤差に関する検討では、一般化交差検証(GCV:Generalized Cross Validation)を指向していた。
【0086】
粒子速度の推定値を得るため、式(9)および(11)によって与えられるSONAH圧力推定値にオイラーの方程式(8)を適用する。その結果として、
【数20】

を得、上式で
【数21】

である。
【0087】
βχ(r)に対する陽式(explicit formula)は、式(27)に式(25)を挿入し、続いて関数IおよびIIに関する式(22)および(23)を適用することによって得ることができる。z成分に関しては、
【数22】

が得られ、xおよびy成分に関しては、
【数23】

が得られる。新たな関数I、II、IxyおよびIIxyは、IおよびIIとともに下表1に記載されている。表1は、基本波動関数の一定重み付け(F(k)=1)を実施する一実施形態においてこれらの関数が有する形も含んでいる。
【0088】
【表1】

【0089】
上表において、JおよびJはそれぞれ0次および1次のベッセル関数である。本発明の発明者は、積分IおよびIxyの実数部に対する解析的式だけを見つけることができたが、これらの積分は一般に、例えば適当な数値積分法、例えばガウス型アルゴリズムを使用した数値積分によって計算することができる。上記の積分のうちの3つの積分における無限大まで及ぶ積分の代わりに、指数関数が十分に低いレベルに達する点においてそれらの積分の上限を打ち切ることができ、または測定格子内の頻度および間隔に基づいて上限を選択することができる。
【0090】
上述のとおり、変数ρおよびζの離散値の2Dテーブルを(例えば数値積分によって)作成し、記憶することによって、多数の積分をかなり高速に実施することができる。特定のρおよびζ値に対する積分を計算するため、次いで、このテーブルの補間を実行しなければならない。あるいは、またはこれに加えて、異なる領域をカバーする2D近似関数を構築してもよい。
【0091】
積分I、IxyおよびIは特異性(singularity)を持たず、そのため、これらの積分は、直接にテーブル形態に入れ、補間することができる。両方のパラメータが0の場合、すなわちこのパラメータ・セットの原点(ρ,ζ)=(0,0)にある場合、積分II、IIxyおよびIIは特異性を有する。補間スキームがこの原点の近くで適正に機能するようにするためには、テーブル化されたこれらの関数の中の特異性を抽出すべきである。
【0092】
波動関数の方向性重み付けの場合を考えると、非常に小さなρおよびζ値に関しては以下のようになる。
【0093】
【数24】

【0094】
これらの特異性は、補間の間に以下のように除去することができる。IIの場合には、テーブルに記憶する一組のパラメータの組合せに対する数値積分によってこの積分を解き、それらの値をテーブルに記憶する前に、それらの値に
【数25】

を掛ける。このようにして、スケーリングされた積分に対して補間を実行する。補間の後、その結果を
【数26】

で割る。残りの積分も同様に処理することができる。単位重み付けを用いて、特異性の同様の除去を、それらをテーブル化する前、補間の間に実行する。
【0095】
全ての仮想源平面と平行なそれぞれの平面内に測定点および計算点を配置すると、効率をさらに向上させることができることに言及することも興味深い。このような構成では、それぞれの周波数に対して第2のパラメータζの値が少数だけ生じ、したがって、それぞれのテーブルがそれぞれのζ値に対応する1Dテーブルを作成することが魅力的になる。例えば、単一の仮想源平面、単一の平行な測定平面および測定平面とは異なる単一の平行な計算平面を用いて、2組の1Dテーブルを作成することができる。これは、実行中に、それぞれの周波数に対して効率的に実行することができる。同様に、1つ仮想源平面、2つの平行な測定平面および別個の平行な計算平面を用いて、全ての周波数に対して4組の1Dテーブルを作成することができる。例えば、これらの1Dテーブルは、予め計算した2Dテーブルの補間によって有利に作成することができる。
【0096】
一般に、上述のとおり、測定位置において音圧が測定されているとき、標的位置rの圧力は、
p(r)=p(AA+εI)−1α(r)=(Aα(r))((AA)+εI)−1
として得られる。上式で、pは、所与の位置における測定/既知音圧データのベクトル、εは正則化パラメータである。同様の式
【数27】

からχ方向の粒子速度を得ることができる。上記の公式は、チホノフ型の正則化を使用しているが、同様の正則化スキーム、例えばいわゆる強化チホノフ正則化(Enhanced Tikhonov regularization)を使用した別の公式をしようすることもできる。
【0097】
次に、図2を参照し、図1を続けて参照して、上述の方法論に基づく再構成法の一実施形態を説明する。
この方法の最初の初期設定ステップS1では、適当な座標系および仮想源平面101を定義する。
【0098】
一般に、試験準備の間に、調査対象の物体の実際の音源面に沿って位置するように測定位置を選択することができる。例えば、フリー・フィールド条件および平面変換器アレイを仮定すると、実際の音源面の一部分にほぼ平行な測定平面を選択することができる。同様に、3Dアレイを使用する非フリー・フィールド条件の場合には、アレイがかなり扁平であると仮定し、実際の音源面の一部分と平行に近くなるように、アレイの中央平面(median plane)を配置することができる。
【0099】
ステップS1において、この方法は、例えば上述の適当な位置検出システムによって、測定点の位置に関する情報、すなわち一組の位置r、i=1,2,...,Iの座標を受け取ることができる。この方法は次いで、その測定ジオメトリに対する座標系を選択することができる。例えば、平面変換器アレイに関しては、xy平面を、測定平面内に位置するように(すなわち測定平面がz=0平面に位置するように)定義し、原点が平面アレイの中心に位置し、正のz軸が、音源面から離れる方向、すなわち源のない領域内へ向かう方向を指すと定義することができる。同様に、3Dアレイに関しては、z=0平面を、アレイの中央平面となるように定義することができる。しかしながら、別の方法で座標系を選択することもできることが理解される。
【0100】
初期設定ステップS1の間に、この方法はさらに、選択した座標系に対して、例えばxy平面に平行な仮想源平面101を定義することができる。この仮想源平面の位置決めは、自動的に実行し、または少なくとも部分的にユーザ入力に基づいて実行することができる。
【0101】
用途によっては、実際の源ジオメトリが分かっていない、または少なくとも正確には分かっていないことがある。このような状況では、アレイからの所定の距離のところにあるアレイに平行な平面上において音源を再構成することができる。この場合には、アレイと平行な仮想源平面との間の距離を、この方法への入力パラメータ、例えばユーザ定義の入力パラメータとして指定することができる。
【0102】
他の状況では、実際の源ジオメトリが既知であり、それをインポート(import)し、整列させることができることがあり、または、アレイ位置を測定するために使用する位置測定システムを用いて、実際の源ジオメトリを測定することもできる。アレイ位置と物体ジオメトリの両方が既知である場合には、処理ソフトウェアが、これら2つの間の距離を自動的に決定することができる。これに基づいて、この方法は、アレイが平面であると仮定して、例えば実際の音源面の内側の指定された距離のところに、アレイに平行な仮想源平面を自動的に位置決めすることができる。アレイが平面でない場合には、中央平面を定義することができる。
【0103】
一般に、仮想源平面は、物体内の物体の表面から短い距離のところに定義することができ、選択する距離は、測定ジオメトリ、例えば隣接する測定点間の典型的な距離または平均距離に依存することがある。例えば、仮想源平面は、音源面から源物体内へ、すなわち測定点の位置および音場を再構成する標的位置(1つまたは複数)とは反対側の実際の音源面の側に、最も近い隣接測定点間の平均距離よりも大きい距離(例えば約1倍から約3倍)だけ離れるように、選択することができる。最適な正確さは、実際の音源面から後ろに測定格子間隔の約1.5倍の距離を置いたところ、すなわち物理的な音源面よりも測定格子から離れたところに仮想源平面を配置することによって、得ることができることが分かっている。低周波数では、これよりもわずかに大きな距離がしばしばより良好な結果を与えることが分かっている。この再構成法の目的上、仮想源平面の片側は、源がなく均質であるとして取り扱われ、そのため、その仮想源平面に関連付けられた一組の平面伝搬波およびエバネッセント波によって、音場を表現することができる。
【0104】
ステップS1の初期設定はさらに、適用する重み付けスキーム、適用する正則化スキーム、使用する測定システムのタイプおよび特性、例えば変換器の数、音場を再構成する標的位置の座標、計算する1つまたは複数の音響特性、周波数範囲、関数テーブルの補間する部分の初期設定(他の部分には、この方法を実行するソフトウェアを供給することができる)など、1つまたは複数のパラメータの設定を含むことができる。例えば、解析ユニット103の適当なユーザ・インタフェースを介して、これらのパラメータおよび追加のまたは代替のパラメータをユーザが定義することができる。
【0105】
ステップS2において、この方法は、一組の位置rにおいて測定したそれぞれの複素時間調和音圧p(r)を受け取る。rは、上記の座標系に対するそれらの位置、例えば位置検出デバイス110から受け取った入力から決定される変換器位置の座標ベクトルを表す。この説明の目的上、位置rは、空気、水などの均質流体によって占められている源のない領域Ωに位置すると仮定する。
【0106】
上で論じたフリー・フィールド状況では、全ての源が平面z=0の一方の側にあり、原則として、例えばこの平面内での測定に基づいて音場を完全に再構成することができるが、本明細書に記載した方法を他の任意の3D測定格子に関して適用することもできる。さらに、後述するとおり、この方法を、源平面が複数ある状況に適用することもできる。
【0107】
複素時間調和音圧p(r)は、受け取った変換器信号からさまざまな方法で決定することができる。例えば、全ての変換器を用いた1回の同時記録を実施し、続いて、全てのFFT線(FTT line)に対する複素(時間調和)音場を提供する全ての信号のFFTを実施することができる。あるいは、基準信号から全てのアレイ変換器への相互スペクトル(cross spectrum)および基準信号の自己スペクトル(auto−spectrum)の測定を実行することもできる。相互スペクトルを自己スペクトルの平方根で割ったものは、適用したスペクトルの全ての周波数線に対する複素(時間調和)音場を提供する。図2の方法は、測定音圧に基づくΩにおける音場の再構成を実行する。代替実施形態は、音圧勾配、粒子速度またはこれらの組合せなどの1つまたは複数の代替のまたは追加の音響量を測定することができる。
【0108】
ステップS3において、この方法は、式(24)を使用して、すなわち
【数28】

および
【数29】

でそれぞれ評価した2次元関数IとIIの1次結合として、自己相関行列AAを計算する。そのために、この方法は、ρおよびζに対するそれぞれの値を、それぞれIおよびIIの予め計算した関数値のルックアップ・テーブル112への索引として使用する。このルックアップ・テーブルを、解析ユニットの記憶領域または他の適当な記憶媒体に記憶することができる。このようにしてこの方法は、ρおよびζに最も近い記憶された関数値を得ることができ、ρおよびζの2次元1次補間、または他の適当な補間スキーム、例えば3次補間を実行することができる。
【0109】
あるいは、この方法は、その代わりに、転置行列BBまたは他の適当な変換行列を計算することができることが理解される。
【0110】
ステップS4において、この方法は、ベクトルqを、
q=((AA)+εI)−1
またはその転置
=p(AA+εI)−1
として計算する。
【0111】
ここで、εは、上述のとおり、正則化パラメータである。この正則化パラメータは、指定された信号対雑音比(SNR)から計算することができ、または雑音が多い(noisy)測定データに基づいて、一般化交差検証(GCV)、L曲線法または同種のパラメータ推定手順を使用して、自動的に計算することができる。SNRが指定されているときには、行列BB内の対角要素の平均値または最大値に10^(SNR/10)を掛けた積として、このパラメータεを得ることができる。フリー・フィールド条件の場合、単一の仮想源平面だけが使用されるときには全ての対角要素が等しく、そのため、平均値または最大値の代わりに単一の対角要素の値を使用することができる。
【0112】
上記の行列演算は、行列を逆行列にするなど行列方程式を解く適当な任意の数値技法によって実行することができる。
【0113】
次いで、この方法は、計算されたq値から音圧を計算することができる。音圧の代わりに、または音圧に加えて、この方法は、別の音響量、例えば上述の粒子速度を計算することができることが理解される。圧力および粒子速度の計算値から音の強さを得ることができる。
【0114】
したがって、この方法は、音場を再構成する一組の標的点rのそれぞれの標的点における音圧および/または粒子速度を計算することができる。標的点の2つの例が図1の109に示されている。多くの状況において、調査対象の物体の音源面で音場を再構成することが重要である。例えば、この方法は、規則格子の一組の格子点のそれぞれの格子点で音圧を計算することができる。この説明の目的上、再構成された音場値を計算する点を計算点と呼ぶ。
【0115】
この方法は、全ての計算点を循環し、所望の音響量を計算することができる。具体的には、ステップS5において、この方法は、現在の計算位置rで所望の音響量を計算する。
【0116】
例えば、この方法は、ベクトル
φp(r)≡Aα(r)
および/または
Ψχ(r)≡Aβχ(r)
を計算することができ、音圧を
p(r)=qφ(r)
として計算し、かつ/または粒子速度を
χ(r)=qΨχ(r)
として計算する。
【0117】
この方法は、式(25)および(28)を使用して、それぞれ相関ベクトルAα(r)およびAβχ(r)を、上記ステップS3におけるAAの計算に関して説明したのと同様の方法で、すなわちそれぞれ入力パラメータρおよびζのそれぞれの値で評価した2次元関数I、IIおよびI、IIまたはIxy、IIxyの1次結合として計算する。上記のとおり、関数を評価する入力パラメータの値は、2D関数の予め計算された関数値を得るためのそれぞれのルックアップ・テーブル112への索引として使用され、それぞれテーブル化された値I、IIおよびI、IIまたはIxy、IIxy間の3次補間または他のタイプの補間を実行することによって使用される。
【0118】
ステップS6において、この方法は、未処理の計算点が残っているかどうかを判定する。残っている場合には、この方法はステップS5へ戻り、残っていない場合にはステップS7へ進み、ステップS7において、計算した一組の音圧および/または粒子速度を、例えば図形表現として、一連の数値として、および/または他の適当な形式で出力する。あるいは、またはこれに加えて、音圧および粒子速度から音の強さが容易に得られる。
【0119】
以上では、図3aに示す状況のフリー・フィールド条件に関する一実施形態を説明した。しかしながら、本明細書に記載した方法の実施形態は、別の状況においても、例えば源および反射物体が半空間だけに限定されないときにも使用することができる。以下では、非フリー・フィールド条件に対する平面波公式を提供する、したがって式(16)に定義した平面波波動関数が測定領域内の音場に対する完全な基礎をもはや構成しない状況を考慮する、上記の実施形態に対する変更を説明する。
【0120】
最初に、位置zおよびz、z≦z≦zにある平行な2平面101aと101bの間の均質な源のない媒質中で測定を実施する図3bに示すケースを考える。z≦zに関する源を有する音場の部分は、式(16)の基礎関数によって表現することができる。さらにz≧zに関する源を有する音場成分を表現するため、結合された組を与える対応する一組の平面波波動関数を追加する。
【0121】
【数30】

上式では、式(17)のkの定義が依然として使用され、両方の組に対して全方向重み付けが適用される。
【0122】
【数31】

したがって、仮想源平面の追加は、SONAH(またはHELS)法において一組の平面波波動関数を追加することに対応する。式(13)、(14)、(20)および(21)は、行列AAおよびAα(r)の要素に2つの組からの寄与を追加することによって、結合された一組の基礎関数(30)および(31)を考慮したSONAH推定手順を得ることができることを示している。これを実行した結果は次のとおりである。
【0123】
【数32】

粒子速度を推定するに必要なベクトルAβχ(r)は、フリー・フィールド条件の場合とちょうど同じように導出することができ、結果は以下のとおりである。
【0124】
【数33】

上記の公式は測定位置の選択を制限しないが、正確な予測を得るのに平面アレイでは不十分なことがある。このことは、測定平面内において対称な像である源は、それぞれのマイクロホンにおいて等しい圧力を発生させるため、互いに区別することができないことによって、例証することができる。一般に、全く同じ平行な2つの長方形アレイからなる2重層アレイを使用することができる。アレイの両側に源があるときの波に対する重みは、定数係数F0,z+およびF0,z−によって調整することができる。一般に、これらのパラメータの値はともに1として選択されるが、アレイの片側に有意な源が存在しない場合には、これらのパラメータの一方をゼロに設定し、それによって一方の側だけに源が分布するという情報をこの公式に提供することによって、より良好な推定を達成することができる。
【0125】
無限平面セクション(例えばスラブ(slab)またはスライス(slice))の両側に源があるケースから、箱形体積の6つの全ての面に源があるケースへの拡張を、図3cの断面図に示す。この状況では単に、箱のそれぞれの側に属する平面波スペクトルからの行列要素への寄与を追加する。対応する式を下表2に示す。
【0126】
【表2】

【0127】
それぞれの面は半球からの全ての平面波を等しい重みで含むため、2面にだけ源がある場合の式は、全ての方向からの平面伝搬波を考慮することが認められることがある。しかし、これらの2つの面が、全てのタイプのエバネッセント波を含むというわけではない。xおよびy方向に減衰する波は失われる。
【0128】
下記の実施例は、波がxy平面において減衰することを考慮することによって達成可能な正確さの向上を例示するシミュレートした測定を示す。
【0129】
したがって、表1に示した2D関数と同じ2D関数を使用して、図2の方法を、測定位置の複数の面に源があるケースに適用することができる。自己および相互相関行列の要素を構成する相関関数は依然として、これらの2D関数の関数値の1次結合として得ることができる。
【0130】
固有推定誤差
以下では、計算点において、本明細書に開示した再構成法の平面波波動関数の値が、平均して、どれくらいよく予測されるのかの少なくとも1つの推定を決定する方法の一実施形態を説明する。
【0131】
最初に、フリー・フィールド条件を考えると、式(16)によって、一組の平面伝搬波関数およびエバネッセント波関数が与えられる。ベクトルα(r)は、式(7)の中の定義に従って、計算点rにおけるこれらの全ての基本波の真の値を含む。式(18)、(12)および(2)は、行列Aの1つの行が、測定位置を横切る波動関数のうちの1つの波動関数の値を含むことを示し、したがって、式(9)および(11)から、下式のように定義されるベクトル
【数34】

【0132】
【数35】

は、計算点における波動関数の推定値を含むことになる。このベクトルの2ノルムは、予測の平均振幅の測度(measure)である。
【0133】
【数36】

平均過小推定(average underestimation)は、予測ノルム
【数37】

が真のノルム‖α(r)‖よりも小さいことが特徴である。真のノルム‖α(r)‖を見つけるため、式(18)、(12)および(2)から、ちょうどα(r)が計算点における全ての関数値を含むように、行列Aの1つの列が、マイクロホン位置における全ての波動関数の真の値を含むことが分かる。したがって、単にマイクロホン位置rを計算位置rに置き換えることによって、行列AAの対角要素[AA]ijの式(24)を使用して、ααの値を計算することができる。
【0134】
【数38】

ここで、定義(22)および(23)を使用して、上記の2つの積分に対する解析的式を容易に検証することができる。固有平均振幅利得(inherent average amplitude gain)は下式のように定義される。
【0135】
【数39】

推定誤差ベクトル
【数40】

の2ノルムは、式(37)および(38)を適用することによって以下のように計算することができ、
【0136】
【数41】

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

単純に、行列AAおよびAα(r)の要素に対する拡張式、例えば2つの面に源がある場合の式(33)および(34)を適用することにより、以上の説明を、非フリー・フィールド条件に拡張することができる。この場合には、ちょうど式(24)を式(33)に拡張したときのように、ααに対する式(39)を次のように拡張しなければならない。
【0138】
【数43】

特定の源構成に対して認められる相対予測誤差は、以下の2つの理由から、上記の固有誤差とは異なる。
【0139】
1)その特定の音場の平面波振幅スペクトルは、この方法が最適化される平面波振幅スペクトルとは異なる。フリー・フィールド条件に関しては、この最適振幅スペクトルが、仮定した仮想源平面z=z上で|F(k)|の形を有する。
2)特定の測定音場に関しては、予測誤差が、平面波成分全体にわたるRMSとしてではなく、全体の音圧に関して計算される。
【0140】
これらの2つの差異のそれぞれの効果については後に示す。
【0141】
実施例:予測誤差の数値的検討
以下では、本明細書に記載した方法の予測の妥当性領域に関する情報の概要を、いくつかの典型的なアレイ・ジオメトリについて、周波数の関数として説明する。
【0142】
測定誤差は考慮せず、正則化パラメータの自動決定法(GCV、L曲線、...)についても調査しない。この実施例の目的上、式(11)および(26)で使用する正則化パラメータεは、指定されたダイナミック・レンジD(信号対雑音比)から、下式を使用して決定される(例えばR.SteinerおよびJ.Hald、「Near−field Acoustical Holography without the errors and limitations caused by the use of spatial DFT」、Intern.J.Acoust.Vib.6、83〜89(2001)を参照されたい)。
ε=[AA]ij10−D/10 (44)
上式で、[AA]ijは、行列AAの対角要素である。この関係は、平面アレイおよびフリー・フィールド条件に対して導出したものであり、その場合、全ての対角要素は等しい。ここでは、一般に、対角要素の平均値を使用する。Dの値は一般に15から30dBの範囲にある。別の値が指定されない限り、以下ではD=30dBを使用する。
【0143】
実施例A:測定平面内の正確な圧力表現の領域
最初に、2つの異なる不規則平面アレイ設計のアレイ平面z=0における総称(generic)固有誤差分布を調べる。これらの2つの異なる不規則平面アレイ設計の直径はともに0.75mであり、これらの設計はともに、少なくとも20kHzまでのビーム形成(beamforming)用途に対して最適化されている。一方の設計は、ホイール内のスポークとして配置された全く同じ13個の6要素線アレイからなるスポーク・ホイール・アレイ(Spoke Wheel array)である。図5を参照されたい。もう一方の設計は、全く同じ7つの扇形からなり、それぞれの扇形に12個のマイクロホンが不規則だが一様に分布したセクタ・ホイール・アレイ(Sector Wheel Array)(例えばJ.Hald、「Array designs optimized for both low−frequency NAH and high−frequency beamforming」、Proceedings of Inter−Noise 2004を参照されたい)である。SONAHの適用に関して良好な成果を得るため、セクタ・ホイール・アレイの一様分布は、アレイ・ジオメトリ最適化の間中、維持した。これらのアレイは平面であるため、フリー・フィールド条件を仮定し、単一の仮想源平面だけを使用する。平均要素間隔は、要素あたりの平均アレイ面積を計算することによって定義し、両アレイ・ジオメトリについて、平均要素間隔は約7cmである。規則的な長方形アレイに関して、この間隔は、約2kHzまでの周波数範囲をサポートすると考えられ、その周波数まで空間エイリアシング(spatial aliasing)を回避するため、アレイ格子間隔よりも小さくないスタンドオフ距離(standoff distance)を選択すべきである。仮想源平面は、z=−7cmにある実際の音源面から、マイクロホン間隔およそ1つ分後ろに配置すべきであるため、アレイ平面から14cm、すなわちz=−14cmのところの仮想源平面を選択する。
【0144】
図5の輪郭プロットは、アレイ平面内の音圧予測の固有相対誤差水準E(r)を示し、マイクロホン位置は丸によって示されている。白色のエリアは、固有誤差水準E(r)が−20dB未満であることを示し、灰色の輪郭は、黒で表現された−6dBよりも大きな誤差水準への2dB刻みでの正確さの低下を示す。低周波数では、両方アレイが、アレイ・エリア全体にわたって良好な音圧表現を提供することができるが、1kHzでは既に、スポーク・ホイール・アレイのマイクロホン分布の隙間が大きすぎて、それらの隙間において良好な表現を提供することができない。セクタ・ホイール・アレイは、2kHzまでずっと、アレイ・エリア全体にわたって許容できる圧力表現を提供するが、スポーク・ホイール・アレイ・ジオメトリに関してはそうでないことが明らかである。しかし、スポーク・ホイール・アレイはビーム形成用途に対して非常に有効であり、大きなサイズでは、その支持構造の構築がより容易である。
【0145】
実施例B:測定平面の外側の正確な圧力表現の領域
次に、アレイ平面に対して垂直な平面内における正確な音圧予測の領域を調べる。ビーム形成用に設計した不規則アレイの代わりに、xおよびy方向の要素間隔が3cmで、0≦x≦21cm、0≦y≦21cmのエリアをカバーする規則的な8×8要素の単層(8×8×1)および2重層(8×8×1)アレイを考える。これらの両タイプのアレイはxy平面z=0にマイクロホン層を有し、2層アレイは、z=3cmに全く同じ追加の層を有する。アレイの対称平面であるxz平面y=10.5cmにおける相対音圧予測誤差に注目する。ちょうど不規則平面アレイと同様に、単層規則アレイに関しては、アレイから、アレイ要素間隔の2倍に等しい距離にある仮想源平面、すなわちz=−6cmにある仮想源平面を選択する。2重層アレイに関しては、前(front)仮想源平面に対しては同じ距離、後(rear)仮想源平面に対しては、要素間隔の19倍の距離、すなわちz=−6cm、z=60cmを選択する。
【0146】
図6の左列は、(単一の仮想源平面を有する)単層アレイの固有誤差水準E(r)を示し、右列は、(2つの仮想源平面を有する)2重層アレイの固有誤差水準E(r)を示す。中央の列は、特定の音場、すなわち(x,y,z)=(7.5,7.5,−6.0)cmにおける単極子点源の音場に対する、単層アレイによる音圧予測の相対誤差を示す。このマッピング・エリアは、x軸に沿って両方の方向に、アレイ・エリアを過ぎて6cm延び、z方向には、z=−3cmから、正のz方向にアレイ要素間隔の4倍の距離まで延びる。明らかに、特定の音場予測の誤差分布は、総称固有誤差よりも不規則であり、一般に総称固有誤差よりも低い。あるエリアでは、固有誤差水準によって予測される誤差よりもかなり小さい誤差が達成される。これは、異なる平面波成分からの誤差寄与が互いに打ち消しあうためと考えることができる。固有誤差水準はこのような相殺効果を回避し、本発明の方法の推定の正確さのはるかに一般的な全体像を提供する。2重層アレイは、単層アレイに比べて、正確な表現の領域をかなり拡張するように見えるが、追加の仮想源平面を追加することによって、音場の自由度も増大した。低い予測誤差を有する領域のz方向の広がりがいくぶん制限される理由は、SONAHが、xy方向の源分布については何も仮定していないためである。良好な遠距離場推定を提供するため、源エリアのサイズについてのなんらかの仮定をこの方法に組み込むことができる。これは例えば、平面伝搬波およびエバネッセント波の代わりに単極子源の限定された分布に基づくSONAH公式を使用することによって達成することができる。
【0147】
実施例C:パラメータの変更に対する予測の正確さの感度
パラメータの変更に対する感度の調査におけるデータ量を限定するため、アレイに垂直な中心に沿った予測誤差だけに注目する。実施例Bの単層8×8×1アレイだけを考える。図7aは、図6の左列に示した条件と同じ条件に対する固有誤差水準を示し、そのため、図6の左列の輪郭プロットの中央を通る一組の水平スライスだけしか示されていない。しかし、図7aにはいくつかの周波数、例えば間隔3cmの規則アレイ格子のより高い5kHzの限界周波数が追加されている。明らかに、5kHzでは、測定平面でさえも、固有誤差レベルの急速な増大が始まっているが、依然として許容できる予測を達成することができる。図7bは、対応する振幅利得A(r)の曲線を示す。5kHzでさえも、アレイ平面からアレイ格子間隔1つ分の範囲内では、依然として振幅は正確に予測されるが、その間隔の外側では、次第に増大する過小推定が観察される。
【0148】
図7では、ダイナミック・レンジDは30dBに固定した式(44)から正則化パラメータを決定した。次のステップは、ダイナミック・レンジDを変化させたときに予測誤差がどれくらい変化するのかを見ることである。図8の左列は、Dの変化に対する固有誤差水準Eの感度を、低周波数(500Hz)および高周波数(4kHz)について示す。周波数ごとに、10dBから40dBの範囲の対応するDの値が記された一組の曲線が示されている。一般に、測定平面よりも源領域から離れた位置での圧力の予測は感度が低く、一方、源領域により近い位置での予測は、D値の増大(より小さな正則化パラメータ値に対応する)とともに向上し、特に低周波数でそうである。より小さなD値はエバネッセント波の再構成を制限し、エバネッセント波は、低周波数においておよび源領域の近くで音場のより大きな部分を構成することから、このことを予想することができる。図8の右列は、SONAH推定における対応する固有振幅利得Aの曲線を示す。500Hzでは、源領域の直ぐ近くを除いて、Dの変化に対する振幅予測の感度はあまり高くないが、4kHzでは、D値が約20dBよりも小さくなったときに、全ての位置で(源領域から離れた位置でも)、次第に増大する過小推定が観察される。このことは、高周波数では平面伝搬波が平面波スペクトルのより大きい部分を構成し、したがってSONAH自己相関行列AAの特異値スペクトルのより大きい部分を構成することによって説明されるであろう。Dがある点よりも小さくなる(正則化パラメータが増大する)と、伝搬波成分の部分さえも、正則化によってかなり減衰する。
【0149】
以前に定義した固有振幅利得Aおよび誤差水準Eは、推定位置における誤差を、全ての平面伝搬波およびエバネッセント波にわたって平均したものを表し、エバネッセント波成分には、SONAH推定アルゴリズムの導出に使用するのと同じ重みが与えられる。どちらの場合も、波動関数は、同じ仮想源平面(1つまたは複数)(フリー・フィールド条件ではz=z)内においてスケーリングされる。次に調べるべき課題は、SONAH計算において使用する仮想源距離が、測定用途においてエバネッセント波の実際の内容を表していない場合に、誤差水準がどれくらい変化するかということである。フリー・フィールド条件および単一の仮想源平面の使用だけを考える。この調査を実施するため、式(38)および(41)において、相関行列AAおよびAαが推定誤差を平均する平面波場を表し、推定重みベクトルc(r)がSONAH推定を表すことに注目する。したがって、1)(38)および(41)でのAAおよびAαの計算ならびに2)SONAH推定重みベクトルc(r)の計算に対して異なる仮想源距離を使用することによって、1)測定におけるエバネッセント波の実際の内容と2)SONAH計算におけるエバネッセント波の重みとの間の不整合の程度の変化に対するSONAH推定の感度を調べることができる。図9に示す結果に関しては、実際の音場に対して、すなわち式(38)および(41)で使用するAAおよびAαの計算において、以前に使用した仮想源平面z=z=−6cmを維持し、SONAH処理に対しては、変化する仮想源平面距離dを使用した。大きなd値は、SONAH推定におけるエバネッセント波に対する重みを小さくする。図9の左列は、dの変化に対する固有誤差水準Eの感度を、低周波数(500Hz)および高周波数(4kHz)について示す。周波数ごとに、3cmから12cmの範囲の対応するdの値が記された一組の曲線が示されている。右列は、対応する振幅利得の曲線を含む。全ての場合に、正則化パラメータは、ダイナミック・レンジDを30dBに固定することによって指定される。
【0150】
500Hzでは、過度に小さな仮想源距離、すなわち再構成におけるエバネッセント波に対する過度に大きな重みに対する誤差水準の感度はいくぶん低いが、過度に大きな距離は、最終的に、測定平面よりも源領域から離れた位置において誤差水準のわずかな過大推定(A>1)をもたらす。これは、波数漏れ(wave−number leakage)によって、小さなエバネッセント波成分が伝搬成分として処理され、それによってそのエバネッセント波成分がより離れた場に寄与することによって説明される。測定点の数がいくぶん少ないため、伝搬波成分とエバネッセント波成分との間の区別はあまり鮮明ではない。(重みベクトルc(r)によって表現される)SONAH処理が、エバネッセント波に非常に小さな重みを与えると、エバネッセント波成分の伝搬成分への漏れは増大する。dが大きいときにも起こるz<0に対する誤差の増大に注目すると、説明は異なってくる。SONAH計算行列においてエバネッセント波には過度に小さな重みが与えられるため、選択された正則化パラメータは、これらの波成分の一部が十分に再構成されることを妨げる。4kHzにおいて、エバネッセント波は一般に、音場のはるかに小さな部分を構成し、そのため、SONAH計算においてエバネッセント波により大きなまたはより小さな重みを与えることに対する感度はより小さい。最大の効果は、再構成において、実際の平面波スペクトルが示唆するものよりも大きな重みをエバネッセント波に与える過度に小さな仮想源距離を使用するときの振幅利得に対する効果である。源領域の近くでは過大推定が観察され、より遠くでは、通常よりも大きな過小推定が見られる。これは、500Hzでの過大推定に対して示した説明と同じ説明によって説明され、ここでは、漏れが反対方向に進む。一部の小さな伝搬波成分がエバネッセント波として処理され、したがってそれらは源エリアから離れると減衰し、源エリアの近くで増幅される。
【0151】
実施例D:アレイ平面内の源−6つの仮想源平面の使用
図10は、2つの面だけからではなく6つの全ての面から平面波波動関数を追加し、それによってxy平面内の方向に沿って減衰するエバネッセント波成分の適当なモデル化を可能にすることによって得られる正確さの向上の一例を示す。実施例Bに関して説明したものと同じ8×8×2 2重層アレイ・ジオメトリを用いたシミュレートされた測定を実行する。(x,y,z)=(−6.0,7.5,1.5)cm、すなわち2つのアレイ層間の中心平面内のアレイ・エリアから6cm外側に、単一の単極子点源を配置する。明らかに、この単極子は、xy平面内の方向に沿って減衰するエバネッセント波型を生成する。SONAHを使用して、前アレイ平面z=0cm内の測定位置における強さベクトル(intensity vector)を、128個の測定圧力値に基づいて予測する。次いで、下式を使用して相対平均誤差水準を計算する。
【0152】
【数44】

上式で、
【数45】

は、計算点番号iにおける真の強さベクトル、Iは、同じベクトルのSONAH推定値であり、求和はi全体についてである。ダイナミック・レンジ30dB、D=30dBを使用してSONAH計算を実行し、全ての方向にアレイから6cm離れた以下の仮想源平面を選択する。
=9cm、z=−6cm、x=27cm、x=−6cm、y=27cm、y=−6cm
【0153】
したがって、単極子は、これらの6つの仮想源平面によって定義された箱の境界上に位置する。図10を見ると、低周波数における強さ推定は、xy平面に平行な仮想源平面だけ存在するときよりも、全ての6つの仮想源平面が存在するときの方がかなり正確であることが分かる。しかしながら、周波数が増大するにつれてこの向上は消滅する。これは、高周波数ではエバネッセント波が音場のより小さな部分を構成するためである。
【0154】
いくつかの実施形態を詳細に説明し、示したが、本発明はそれらの実施形態に限定されるもではなく、下記の特許請求の範囲に定義された主題の範囲に含まれる他の方法で、本発明を具体化することもできる。例えば、本明細書に開示した方法の実施形態は主として、測定音響量が音圧である実施形態に関して説明した。しかしながら、他の実施形態では別の音響量を測定することができることが理解される。例えば、音圧の代わりに、測定位置rにおけるχ方向の粒子速度uχ(r)を測定することもできる。このような一実施形態では、展開係数ベクトルaについて解く1次方程式系が用意されているとき、式(1)の代わりに、位置rに対して以下の式を使用することができる。
【0155】
【数46】

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

これは、行列AおよびBにおいて、関数Ψ(r)が、Γχ,n(r)、n=1,2,...,Nに置き換えられ、式(3)のベクトルpにおいて、圧力p(r)が測定粒子速度uχ(r)に置き換えられたことを意味する。次いで、これに従って、自己相関行列のAAの含まれる要素を変更しなければならない。
【0157】
本明細書に記載した方法および装置を使用して、例えば機械、モータ、エンジン、自動車などの車両および/または他の同種の振動物体の音響特性を解析するときに、これらの振動物体などのさまざまな音源/雑音源の音場を再構成することができる。
【0158】
本明細書に記載した方法の実施形態は、異なるいくつかの要素を含むハードウェアによって実現することができ、かつ/または適当にプログラムされたマイクロプロセッサによって少なくとも部分的に実現することができる。
【0159】
いくつかの手段を列挙した装置請求項では、これらの手段のいくつかを、同じ1つのハードウェア要素、構成要素またはアイテムによって具体化することができる。単に、ある種の手段が相異なる従属請求項に記載されていること、またはある種の手段が異なる実施形態に記載されていることは、これらの手段の組合せを有利に使用することができないことを意味しない。
【0160】
本明細書で使用されるとき、用語「含む(comprises/comprising)」は、明示された特徴、要素、ステップまたは構成要素の存在を指定すると解釈されるが、この用語は、1つまたは複数の他の特徴、要素、ステップ、構成要素またはこれらの集団の存在または追加を排除しないことを強調しておく。

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

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate

【図7a】
image rotate

【図7b】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate


【公表番号】特表2011−527422(P2011−527422A)
【公表日】平成23年10月27日(2011.10.27)
【国際特許分類】
【出願番号】特願2011−517085(P2011−517085)
【出願日】平成21年6月26日(2009.6.26)
【国際出願番号】PCT/EP2009/058039
【国際公開番号】WO2010/003837
【国際公開日】平成22年1月14日(2010.1.14)
【出願人】(504349973)ブリュエル アンド ケアー サウンド アンド ヴァイブレーション メジャーメント エー/エス (6)
【Fターム(参考)】