逆二乗特性解析装置及び逆二乗特性評価方法
【課題】測定した音圧値から適切な減衰曲線の推定を行うことが可能な逆二乗特性解析装置を提供する。また、その減衰曲線の推定を少ない演算量によって行うこと、および、適切な推定結果を安定的に得ることが可能な逆二乗特性解析装置を提供する。
【解決手段】逆二乗特性解析装置１１は、試験音場内の複数の位置における音圧の測定値とそれぞれの測定点の位置を入力データとして、音圧と位置の関係を反比例で近似した減衰曲線を推定する減衰曲線推定手段１２と、推定した減衰曲線と音圧の測定値との差を音圧レベルの尺度において算出する差分算出手段１３と、算出した差に基づいて試験音場の特性を評価する特性評価手段１４を備え、さらに、減衰曲線推定手段１２は、音圧の測定値と減衰曲線との差を音圧レベルの尺度において最小化するように減衰曲線の推定を行う。
【解決手段】逆二乗特性解析装置１１は、試験音場内の複数の位置における音圧の測定値とそれぞれの測定点の位置を入力データとして、音圧と位置の関係を反比例で近似した減衰曲線を推定する減衰曲線推定手段１２と、推定した減衰曲線と音圧の測定値との差を音圧レベルの尺度において算出する差分算出手段１３と、算出した差に基づいて試験音場の特性を評価する特性評価手段１４を備え、さらに、減衰曲線推定手段１２は、音圧の測定値と減衰曲線との差を音圧レベルの尺度において最小化するように減衰曲線の推定を行う。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、半無響室および無響室の自由空間特性を評価する解析装置及び方法に関する。
【背景技術】
【０００２】
半無響室および無響室は、自由空間（音の反射のない空間）を模した試験室であり、音の測定環境として広く利用されている。半無響室および無響室が自由空間にどの程度近い特性を有しているかという点は、その試験室における測定精度に影響を及ぼすため、その特性の適切な評価方法の確立は極めて重要な技術課題である。半無響室および無響室の自由空間としての特性を評価する方法としては、ISO 3745（非特許文献１）の附属書Aに記された逆二乗則検定が一般的に利用されている。
【０００３】
＜逆二乗則検定の概要＞
以下では、逆二乗則検定の概要について説明する。
自由空間に置かれた無指向性の点音源から放射される音の強さ（音響インテンシティ）は、音源からの距離の二乗に反比例して減衰する（逆二乗則）。また、このときの音の強さ（音響インテンシティ）は音圧の二乗に比例するため、音圧は音源からの距離に反比例して減衰する。
【０００４】
従って、無響室内に無指向性の点音源を設置して、様々な測定点において音圧を測定した場合、その無響室が完全な自由空間であるならば、各測定点で測定される音圧は音源からの距離に反比例したものとなる。しかしながら、壁面からの反射が生じる場合には、測定される音圧は、重畳される反射波の度合いに応じて反比例曲線から乖離するようになる。そこで、その乖離量を算出し、その大きさによって自由空間としての特性評価を行うことが可能となる。
【０００５】
半無響室については、自由空間に床からの反射が重畳された特性（反射面上の自由空間）が理想特性となるため、その室内に無指向性の点音源を設置したとしても、床からの反射波との干渉により、音圧は必ずしも距離に反比例して減衰しない。しかし、音源を床面上に設置した場合には、床面上方へ半球状に放射される直接波と床面から鏡面反射され上方に半球状に広がる反射波が完全に重なり合い、両波面の干渉による減衰特性の変化が生じないため、自由空間と同様に扱うことが可能となる。すなわち、半無響室が反射面上の自由空間として完全な特性を持つならば、音源を床面上に置いた場合には、各測定点で測定される音圧は音源からの距離に反比例したものとなる。しかしながら、床面以外の壁面からの反射が生じる場合には、測定される音圧は、重畳される床面以外の反射波の度合いに応じて反比例曲線から乖離するようになる。つまり、半無響室の場合も、音源を床面上に置いて測定を行うことによって、無響室と同様に、反比例曲線からの音圧値の乖離量による特性評価が可能となる。
【０００６】
ISO 3745には、この基本原理に基づいて、点音源を設置した試験室内において測定した音圧データから理想的な減衰曲線を推定し、その減衰曲線と測定値との乖離量によって、その試験室の自由空間としての特性評価を行う方法が記載されている。
以上が、逆二乗則検定の概要である。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開2006-350948号公報
【非特許文献】
【０００８】
【非特許文献１】国際規格 ISO 3745:2003
【非特許文献２】茨木俊秀・福島雅夫、最適化の手法（共立出版）ISBN4-320-02664-0
【発明の概要】
【発明が解決しようとする課題】
【０００９】
ISO 3745に記載された方法に基づいて解析を行う従来の逆二乗特性解析装置では、測定した音圧値から適切な減衰曲線が推定されず、理想的な減衰特性からの測定値の乖離量が正確に算出されないことがあった。その場合、半無響室および無響室が良好な特性を有していても、検定結果は不適合となることがあり、合理的な検定結果が得られないという問題があった。
【００１０】
本発明は、このような課題を解決するためになされたものであり、測定した音圧値から適切な減衰曲線の推定を行うことが可能な逆二乗特性解析装置を提供することを目的とする。またさらに、その減衰曲線の推定を少ない計算量によって行うこと、および、適切な推定結果を安定的に得ることが可能な逆二乗特性解析装置を提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記の目的を達成するために、本発明の第一の実施態様による逆二乗特性解析装置は、試験音場内の複数の測定点における音圧の測定値とそれぞれの測定点の位置とを入力データとして、音圧と位置の関係を反比例で近似した減衰曲線を推定する減衰曲線推定手段と、前記推定された減衰曲線と前記音圧の測定値との差を音圧レベルの尺度において算出する差分算出手段と、前記算出した差に基づいて前記試験音場の自由空間特性を評価する特性評価手段を備える逆二乗特性解析装置であって、前記減衰曲線推定手段は、前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように前記減衰曲線の推定を行うことを特徴とする。
【００１２】
本発明の第二の実施態様による逆二乗特性解析装置は、試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、前記pi、ri (i=1〜N)を入力データとして、式(1)で表される減衰曲線の未知パラメータa、r0を推定する減衰曲線推定手段と、前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出する差分算出手段と、前記算出された差が所定の許容値に収まるか否かによって前記試験音場の自由空間特性を評価する特性評価手段を備える逆二乗特性解析装置であって、前記減衰曲線推定手段は、式(3)を最小にするように前記未知パラメータa、r0の推定を行うことを特徴とする。
【数１】
【数２】
【数３】
【００１３】
本発明の第三の実施態様による逆二乗特性解析装置は、上記第二の実施態様による逆二乗特性解析装置であって、前記減衰曲線推定手段は、前記未知パラメータr0に初期値を設定し、式(4)、式(5)及び式(6)に従って前記未知パラメータr0を更新する処理を前記未知パラメータr0が収束するまで繰返し行って前記未知パラメータr0の収束解を算出し、前記算出した未知パラメータr0の収束解を用いて、式(7)に従って前記未知パラメータaを算出することを特徴とする。
【数４】
【数５】
【数６】
【数７】
【００１４】
本発明の第四の実施態様による逆二乗特性解析装置は、上記第二の実施態様による逆二乗特性解析装置であって、前記減衰曲線推定手段は、式(8)、式(9)、式(10)及び式(11)に従って、前記未知パラメータa、r0の推定を行うことを特徴とする。
【数８】
【数９】
【数１０】
【数１１】
【００１５】
本発明の第五の実施態様による試験音場内の逆二乗特性の評価方法は、試験音場内の複数の測定点における音圧の測定値とそれぞれの測定点の位置とを入力データとして、音圧と位置の関係を反比例で近似した減衰曲線を推定するステップと、前記推定された減衰曲線と前記音圧の測定値との差を音圧レベルの尺度において算出するステップと、前記算出した差に基づいて前記試験音場の自由空間特性を評価するステップからなり、前記減衰曲線を推定するステップでは、前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように前記減衰曲線の推定を行うことを特徴とする。
【００１６】
本発明の第六の実施態様による試験音場内の逆二乗特性の評価方法は、試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、前記pi、ri (i=1〜N)を入力データとして、式(1)で表される減衰曲線の未知パラメータa、r0を推定するステップと、前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出するステップと、前記算出された差が所定の許容値に収まるか否かによって前記試験音場の自由空間特性を評価するステップからなり、前記減衰曲線を推定するステップは、式(3)を最小にするように前記未知パラメータa、r0の推定を行うことを特徴とする。
【数１２】
【数１３】
【数１４】
【発明の効果】
【００１７】
本発明の逆二乗特性解析装置によれば、減衰曲線推定手段において音圧の測定値と減衰曲線との差を最小化するように減衰曲線の推定を行うときの尺度が、差分算出手段および特性評価手段において音圧の測定値と減衰曲線の差を評価するときの尺度と一致するため、適切な減衰曲線を推定することが可能となる。その結果として、測定値の減衰曲線からの乖離量が精度よく算出されるため、合理的な検定結果を得ることが可能となる。
【図面の簡単な説明】
【００１８】
【図１】半無響室における逆二乗則検定の様子を説明する模式図である。
【図２】前提技術の逆二乗特性解析装置を説明するブロック構成図である。
【図３】前提技術の逆二乗特性解析装置による解析結果の一例を示すグラフである。
【図４】実施例１の逆二乗特性解析装置を説明するブロック構成図である。
【図５】実施例１の逆二乗特性解析装置の減衰曲線推定部において、未知パラメータr0を算出する処理手順を説明する処理フロー図である。
【図６】実施例１の逆二乗特性解析装置による解析結果の一例を示すグラフである。
【図７】前提技術の逆二乗特性解析装置による解析結果の別の一例を示すグラフである。
【図８】実施例１の逆二乗特性解析装置による解析結果の別の一例を示すグラフである。
【図９】実施例２の逆二乗特性解析装置を説明するブロック構成図である。
【図１０】実施例２の逆二乗特性解析装置による解析結果の一例を示すグラフである。
【図１１】実施例２の逆二乗特性解析装置による解析結果の別の一例を示すグラフである。
【発明を実施するための形態】
【００１９】
まず、以下では、本発明の実施の形態に係る前提技術を説明する。ここでは、ISO 3745に記載された方法を例とした逆二乗則検定の実施手順、および、その中で用いられる逆二乗特性解析装置の構成と動作について説明する。
【００２０】
＜逆二乗則検定の実施手順＞
以下では、ISO 3745に従って逆二乗則検定を行う手順について説明する。なお、ISO 3745は半無響室と無響室を検定対象としているが、ここでは半無響室の検定を例に説明を行う。
【００２１】
図１は検定対象となる半無響室の模式図である。図１を参照して、逆二乗則検定では、試験音場である半無響室の床面の中央に発音体である音源（スピーカ）を設置し、ここから試験音場内（半無響室内）へ音の放射を行う。さらに、この音源を端点とする直線経路上にいくつかの測定点を設定し、それぞれの測定点において音圧を測定する。
【００２２】
直線経路は音源から半無響室上部のいずれかの角に向かうものとし、測定点については、直線経路上の音源から0.5mの位置を第1番目の測定点として、そこから順に直線経路に沿って0.1 m毎に合計N個の測定点が設定されるものとする。直線経路および測定点の位置については、ISO 3745にその決め方が定められており、上記の設定はこれに従うものである。また、このときの第i番目(i=1〜N)の測定点の音源からの距離をri (i=1〜N)、測定された音圧値をpi (i=1〜N)とする。ここで、距離ri の単位はm（メートル）、音圧piの単位はPa（パスカル）である。但し、ISO 3745では所定の周波数帯域（1/1もしくは1/3オクターブ帯域）毎に検定を行うことが定められており、上記のpiは検定対象とする周波数帯域の音圧値を表すものとする。
【００２３】
なお、上記の測定に用いる収音システム（マイクロホンなど）に求められる仕様や、検定を行う各帯域の周波数範囲、帯域分割フィルタの特性などは、いずれもISO 3745に定められており、ここでもその定めに従って測定を行う。また、音源（スピーカ）についても無指向性で点音源に近い特性が要求されるが、ISO 3745に具体的な要求仕様が定められており、それに従うものとする。このような測定は、規格に従い適切な測定装置を使用することで一般的に実施が可能である。
【００２４】
ところで、音圧は、圧力の単位であるパスカル（単位記号Pa）や、20μPaを基準として式(2)に従う対数尺度で示したデシベル（単位記号dB）など、異なる単位で表されることがある。これらは相互に変換することが可能であり、いずれの単位で表されていても本質的な違いはないが、説明の便宜のため、本明細書では、パスカルで表される量を音圧（量記号p）、式(2)に従いデシベルで表される量を音圧レベル（量記号Lp）と区別して呼ぶこととする。また、20μPaを基準音圧と定義し、記号p0で表すこととする。さらに、式(12)に従い、基準音圧20μPaで正規化された音圧の逆数を逆数音圧（量記号q）と定義する。
【数１５】
【数１６】
【００２５】
上記に従ってN個の測定点で音圧を測定した後、各測定点の音圧値pi(i=1〜N)と音源からの距離ri (i=1〜N)を後述する逆二乗特性解析装置に入力する。逆二乗特性解析装置は、入力された測定値の距離減衰特性が反比例曲線からどの程度乖離しているかを解析し、その乖離量に応じて、適合または不適合の結果を出力する。出力結果が適合である場合には、検定を行った周波数帯域において、音源から第N番目の測定点の位置までの範囲で良好な逆二乗特性が成立しており、その範囲は自由音場とみなされる。出力結果が不適合である場合は、検定を行った周波数帯域において、音源から第N番目の測定点の位置までの範囲は自由音場とみなされない。
以上の手順によって、ISO 3745に従う半無響室の逆二乗則検定が行われる。
【００２６】
＜前提技術に係る逆二乗特性解析装置の構成＞
以下では、前提技術に係る逆二乗特性解析装置の構成を説明する。この逆二乗特性解析装置は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして受け、その距離減衰特性が反比例曲線から乖離する量を算出し、その乖離量に基づいて自由空間としての適合／不適合の解析結果を出力するものである。なお、その解析方法は、ISO 3745に記載された方法に従うものである。
【００２７】
図２は、前提技術に係る逆二乗特性解析装置を説明するブロック構成図である。図２を参照して、逆二乗特性解析装置１は、測定値から反比例曲線を推定する減衰曲線推定部２と、推定された反比例曲線と測定値との乖離量を算出する差分算出部３と、算出された乖離量から半無響室の自由空間特性を判定する特性評価部４を備える。
【００２８】
減衰曲線推定部１は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、式(1)でモデル化された反比例曲線の推定を行う。式(1)は、音圧pが距離rに反比例する関係を示した式であるが、未知のパラメータa、r0を含んだ形でモデル化されている。減衰曲線推定部１では、式(13)、式(14)に従って未知パラメータa、r0を算出し、反比例曲線の推定を完了する。但し、qi(i=1〜N)は第i番目の測定点における逆数音圧であり、式(12)に従ってpiから算出される。なお、式(1)のモデル式、および、式(13)、式(14)の未知パラメータa、r0の推定式は、ISO 3745の記載に基づくものである。
【数１７】
【数１８】
【数１９】
【００２９】
差分算出部３は、推定された未知パラメータa、r0、および、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、各測定点における音圧レベルの測定値Lpi (i=1〜N)、音圧レベルの推定値Lp(ri) (i=1〜N)、および、それらの差分δi (i=1〜N)を算出する。ここで、Lpi およびLp(ri)の算出は、それぞれ、式(15)、式(16)に従って行われ、δiの算出は、式(17)に従って行われる。なお、式(15)は、式(2)に基づいて得られるpiとLpiの関係式であり、式(16)は、式(1)と式(2)に基づいて得られるa、r0、riとLp(ri)の関係式である。
【数２０】
【数２１】
【数２２】
【００３０】
特性評価部４は、差分δi (i=1〜N)を入力データとして、これらのδiに対して式(18)の関係が成立するか否かの判定を行い、全ての測定点、すなわち、i=1〜Nの全てのiにおいて式(18)が成立する場合には適合、そうでない場合には不適合の結果を出力する。なお、δlimは許容限界値であり、式(18)は差分δi の絶対値がこの許容限界値以内に収まることを示す関係式である。δlim の値は、ISO 3745では検定対象とする周波数帯域毎に定められており、ここでもそれに従う。
【数２３】
【００３１】
＜前提技術に係る逆二乗特性解析装置の動作＞
以下では、測定点数が26個（N=26）、各測定点の音圧値pi、距離riが、表１の第３列と第２列に示されたデータである場合を例に、前提技術に係る逆二乗特性解析装置で行われる解析動作を説明する。なお、距離riはri=0.5+0.1(i-1) m の関係があるので、第26番目の測定点は音源から3.0mの位置となり、検定範囲は音源から直線経路に沿って3.0 mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【表１】
【００３２】
逆二乗特性解析装置１において、各測定点における音圧値pi、音源からの距離riが入力されると、減衰曲線推定部２では、式(13)、式(14)に従って未知パラメータa、r0の算出が行われる。このとき、表１の第３列と第２列に示されたpi、riの入力データに対しては、a=1037.32、r0=-0.0105871が算出される。
【００３３】
未知パラメータa、r0が推定されると、差分算出部３では、式(15)と式(16)に従い、音圧レベルの測定値Lpiおよび推定値Lp(ri)が算出され、さらに、式(17)に従い、その差分δiが算出される。このときのLpi、Lp(ri)、δiは、表１の第５列、第６列、第７列にそれぞれ記された値となる。なお、図３は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表１の第５列の値）を丸い点で、推定曲線（式(19)にa=1037.32、r0=-0.0105871を代入した曲線）を実線で示してある。
【数２４】
【００３４】
差分δiが算出されると、特性評価部４では、式(18)に従い、それらの値が許容限界値以内に収まるか否かが判定される。検定対象としている中心周波数500Hzの1/3オクターブ帯域では、許容値δlimは2.5dBと定められており、差分δiの算出結果である表１の第７列の値に対しては全て式(18)が成立するため、適合の解析結果が出力される。この解析結果は、中心周波数が500Hzの1/3オクターブ帯域において、直線経路に沿って音源から3.0 mの位置までの範囲が自由音場とみなされることを示している。
【００３５】
以上によって、前提技術に係る逆二乗特性解析装置の構成、および、各測定点の音圧値pi (i=1〜N)と距離ri (i=1〜N)が入力されてから、自由空間としての適合／不適合を示す結果が出力されるまでの解析動作が説明された。なお、前提技術に係る逆二乗特性解析装置で用いられている解析方法はISO 3745に記載された方法に従うものである。
【００３６】
ここで、式(1)のモデル式の未知パラメータa、r0について説明する。aは音源からの音響出力の大きさを示すパラメータであり、減衰曲線が測定時の音量に依存するために設けられているものである。また、r0は音源の音響中心を補正するためのパラメータであり、実際の音響中心とriの算出時に原点として設定される音響中心とが必ずしも一致しないので、その補正のために設けられているものである。これらのパラメータは、一般的に、事前に知ることが難しいため、式(1)のように未知のパラメータとしてモデル化し、減衰曲線推定部において、測定値からこれらのパラメータを推定して反比例曲線を決定する構成となっている。なお、ここで重要な点は、推定されるべき曲線は、測定環境が完全な自由音場であったとした場合の理想的な減衰曲線であるということである。その理由は、測定値との乖離量を算出する際、そのような曲線を基準としなければ、自由空間としての評価を正しく行うことができないからである。しかし、そのような曲線を厳密に求めることはできないので、各測定点で測定された音圧値の近傍を通る曲線が適切な推定曲線であると仮定し、その曲線は式(13)、式(14)に従って求められるという考えに基づいて、前提技術に係る逆二乗特性解析装置は構成されている。
【００３７】
しかしながら、後述するように、前提技術に係る逆二乗特性解析装置において推定される減衰曲線は、各測定点で測定された音圧値の近傍を通る曲線とはならないことがあり、その結果、正しい検定結果が得られないことがある。このような前提技術に係る課題に対処するため、前提技術のさらなる改良として、本発明に係る実施の形態を以下に説明する。
【実施例１】
【００３８】
以下では、本発明の逆二乗特性解析装置の実施例を説明する。なお、逆二乗則検定を行う手順については、前提技術の逆二乗特性解析装置に代えて、本発明の逆二乗特性解析装置を使用すること以外、前提技術で説明した手順と全て同じであるので説明を省略し、逆二乗特性解析装置の構成と解析動作のみを説明する。
【００３９】
＜本実施例の逆二乗特性解析装置の構成＞
以下では、逆二乗則検定において使用される本実施例の逆二乗特性解析装置の構成を説明する。この逆二乗特性解析装置は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして受け、その距離減衰特性が反比例曲線から乖離する量を算出し、その乖離量に基づいて自由空間としての適合／不適合の解析結果を出力するものである。
【００４０】
図４は、本発明の逆二乗特性解析装置を説明するブロック構成図である。図４を参照して、逆二乗特性解析装置１１は、測定値から反比例曲線を推定する減衰曲線推定部１２と、推定された反比例曲線と測定値との乖離量を算出する差分算出部１３と、算出された乖離量から半無響室の自由空間特性を判定する特性評価部１４を備える。
【００４１】
減衰曲線推定部１１は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、式(1)でモデル化された反比例曲線の推定を行う。式(1)は、音圧pが距離rに反比例する関係を示した式であるが、未知のパラメータa、r0を含んだ形でモデル化されている。
【数２５】
【００４２】
未知パラメータr0は、図５の処理フローで示される繰返し演算によって算出される。図５を参照して、処理ステップＳ１は、r0に初期値を設定するステップであり、ここでは、r0=0が設定されるものとする。処理ステップＳ２は、r0の更新を行うステップであり、現時点でのr0をr0,oldとして、式(4)に従ってr0,oldから新しく更新されたr0を算出する。但し、式(4)の中の∂E/∂r0、∂2E/∂r02は、それぞれ、式(5)、式(6)に従い算出されるものである。処理ステップＳ３は、この繰返し演算が収束したか否かを判定するステップである。処理ステップＳ２において算出したr0とr0,oldの関係が、式(20)の関係を満たす場合には、演算が収束したと判断し、その時点でのr0を最終結果として算出を終える。式(20)を満たさない場合には、再びステップＳ２に戻り、r0を更新する演算を繰り返す。なお、εは演算精度などを考慮して定められる充分小さな正定数であるが、ここではε=10-12であるとする。
【数２６】
【数２７】
【数２８】
【数２９】
【００４３】
このような繰返し演算は、DSP（Digital Signal Processor）やRAM（Random Access Memory）によって構成される一般的な演算装置によって実施することが可能であり、減衰曲線推定部１２はこのような演算装置を備えている。
【００４４】
繰返し演算によって未知パラメータr0が算出されると、式(7)に従って未知パラメータaが算出され、式(1)でモデル化された減衰曲線の推定が完了する。
【数３０】
【００４５】
なお、式(1)のモデル式はISO 3745に基づくものであるが、図５の処理フロー、および、式(4)、式(5)、式(6)、式(20)、式(7)に基づく未知パラメータa、r0の推定方法は、本実施例の逆二乗特性解析装置の特徴となる部分である。
【００４６】
差分算出部１３は、推定された未知パラメータa、r0および、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、各測定点における音圧レベルの測定値Lpi (i=1〜N)、音圧レベルの推定値Lp(ri) (i=1〜N)、および、それらの差分δi (i=1〜N)を算出する。ここで、Lpi およびLp(ri)の算出は、それぞれ、式(15)、式(16)に従って行われ、δiの算出は、式(17)に従って行われる。なお、式(15)は、式(2)に基づいて得られるpiとLpiの関係式であり、式(16)は、式(1)と式(2)に基づいて得られるa、r0、riとp(ri)の関係式である。
【数３１】
【数３２】
【数３３】
【００４７】
特性評価部１４は、差分δi (i=1〜N)を入力データとして、これらのδiに対して式(18)の関係が成立するか否かの判定を行い、全ての測定点、すなわち、i=1〜Nの全てのiにおいて式(18)が成立する場合には適合、そうでない場合には不適合の結果を出力する。なお、δlimは許容限界値であり、式(18)は差分δi の絶対値がこの許容限界値以内に収まることを示す関係式である。δlim の値は、ISO 3745では検定対象とする周波数帯域毎に定められており、ここでもそれ従う。
【数３４】
【００４８】
＜本実施例の逆二乗特性解析装置の動作＞
以下では、前提技術で例示したのと同様に、測定点数が26個（N=26）、各測定点の音圧値pi、距離riが表２の第３列と第２列に示されたデータである場合を例に、本実施例の逆二乗特性解析装置で行われる解析動作を説明する。なお、距離riはri=0.5+0.1(i-1) m の関係があるので、第26番目の測定点は音源から3.0 mの位置となり、検定範囲は音源から直線経路に沿って3.0 mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【表２】
【００４９】
逆二乗特性解析装置１１において、各測定点における音圧値pi、音源からの距離riが入力されると、減衰曲線推定部１２では、図５に示された処理フロー、および、式(4)、式(5)、式(6)、式(20)に従って未知パラメータr0の算出が行われ、さらに、式(7)に従って未知パラメータaの算出が行われる。このとき、表２の第３列と第２列に示されたpi、riの入力データに対しては、a=1029.45、r0=0.0043552が算出される。
【００５０】
未知パラメータa、r0が推定されると、差分算出部１３では、式(15)と式(16)に従い、音圧レベルの測定値Lpiおよび推定値Lp(ri)が算出され、さらに、式(17)に従い、その差分δiが算出される。このときのLpi、Lp(ri)、δiは、表２の第５列、第６列、第７列にそれぞれ記された値となる。なお、図６は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表２の第５列の値）を丸い点で、推定曲線（式(19)にa=1029.45、r0=0.0043552を代入した曲線）を実線で示してある。
【数３５】
【００５１】
差分δiが算出されると、特性評価部１４では、式(18)に従い、それらの値が許容限界値以内に収まるか否かが判定される。検定対象としている中心周波数500Hzの1/3オクターブ帯域では、許容値δlimは2.5dBと定められており、差分δiの算出結果である表２の第７列の値に対しては全て式(18)が成立するため、適合の解析結果が出力される。この解析結果は、中心周波数が500Hzの1/3オクターブ帯域において、直線経路に沿って音源から3.0 mの位置までの範囲が自由音場とみなされることを示している。
【００５２】
以上によって、本実施例の逆二乗特性解析装置の構成、および、各測定点の音圧値pi (i=1〜N)と距離ri (i=1〜N)が入力されてから、自由空間としての適合／不適合を示す結果が出力されるまでの解析動作が説明された。
【００５３】
＜本実施例の逆二乗特性解析装置の特徴と効果＞
以下では、本実施例の逆二乗特性解析装置の特徴と効果について詳しく説明する。
最初に、前提技術において未知パラメータa、r0の推定に用いられる式(13)、式(14)の技術的意味を明らかにする。
【００５４】
各測定点における逆数音圧の測定値qiと推定値q(ri)の差の二乗和をEとすると、Eは式(21)のように表される。推定値q(ri)は、式(12)と式(1)に基づいて式(22)のように表されるので、これを式(21)に代入すると、Eは式(23)のように式変形される。
【数３６】
【数３７】
【数３８】
【００５５】
ここで、式(23)のEを最小にする未知パラメータa、r0は、Eをa、r0でそれぞれ微分して0とおいた連立方程式、すなわち、式(24)を解いて求められる。
【数３９】
【００５６】
ところが、式(24)はa、r0について複雑な方程式となり、容易に解くことができない。そこで、式(25)の変数変換を用いて式(23)を式(26)のように式変形し、式(26)のEを最小化するA、Bを求めた後、式(25)の関係を用いてA、Bからa、r0を求める。
【数４０】
【数４１】
【００５７】
式(26)のEを最小化するA、Bを求めるため、式(26)をA、Bで微分すると式(27)が得られ、これらを0とおいてA、Bについて整理すると式(28)の連立方程式が得られる。式(28)はA、Bの線形方程式として容易に式(29)のように式変形することができ、さらにこれを式(25)の関係式を用いてa、r0について解くと式(30)、式(31)となって、これらが式(21)のEを最小化するa、r0の算出式となる。
【数４２】
【数４３】
【数４４】
【数４５】
【数４６】
【００５８】
式(30)、式(31)は、前提技術においてa、r0の推定で用いた式(13)、式(14)に一致しており、以上の検討によって、前提技術の逆二乗特性解析装置において推定される減衰曲線は、式(21)のE、すなわち、逆数音圧の差分の二乗和を最小化したものであることが明らかになった。
【００５９】
ところで、差分算出部および特性評価部では、測定値と推定値の差分を音圧レベルにおいて算出するので、前提技術の逆二乗特性解析装置で推定される減衰曲線の評価尺度である逆数音圧とは一致していない。以下では、このことが解析結果に及ぼす影響について説明する。
【００６０】
逆数音圧qにおける微小な変化をΔq、音圧レベルLpにおける微小な変化をΔLpとすると、両者の間には近似的に式(32)の関係が成立する。式(2)、式(12)を用いて、音圧レベルLpと逆数音圧qの関係は式(33)と表されるので、Lpのqによる微分dLp/dqは式(34)となり、これを式(32)に代入するとΔqとΔLpの間には式(35)の関係が成立する。
【数４７】
【数４８】
【数４９】
【数５０】
【００６１】
式(35)は、ΔqとΔLpの間に1/pに比例する重み係数が掛かることを示している。つまり、逆数音圧における測定値と推定値の差分を、全ての測定点において一様に最小化した場合、音圧レベルの尺度においては、測定点毎に1/pに比例した重みをかけて最小化したこととなる。1/pは音圧が小さいほど大きな値となるので、前提技術の逆二乗特性解析装置における減衰曲線の推定方法は、音圧レベルで見た場合、音圧が小さい測定点における差分を重点的に小さくするような推定をしていることに相当する。一般的に、音源に近い測定点ほど音圧は大きいので、このような推定がなされた場合、特に音源付近において測定値と大きく乖離した減衰曲線が推定される可能性がある。このような各測定点の近傍を一様には通らない減衰曲線は、測定環境が完全な自由音場であった場合の理想的な減衰曲線を推定しているとは言えず、それを基準として算出された乖離量も自由空間としての特性を正しく反映させたものとはならない。
【００６２】
以上の検討によって、前提技術の逆二乗特性解析装置では、減衰曲線を推定する際に最小化する差分の尺度と、測定値と推定値の乖離量を評価する時の尺度の違いが考慮されておらず、この点が誤った検定結果をもたらす原因となり得ることが示された。
【００６３】
上記の検討結果を考慮すると、減衰曲線推定部においては、逆数音圧における差の二乗和ではなく、音圧レベルにおける差の二乗和を最小にするよう減衰曲線の推定を行うことが有効であると考えられる。これは言い換えると、音圧レベルにおける測定値Lpiと推定値Lp(ri)の差の二乗和、すなわち、式(3)のEを最小化する未知パラメータa、r0を算出することである。
【数５１】
【００６４】
式(2)、式(1)に基づいて得られる式(37)の関係を用いて、式(3)を未知パラメータa、r0を含む形に式変形すると、式(36)が得られる。
【数５２】
【数５３】
【００６５】
ここで、式(36)のEを最小にする未知パラメータa、r0は、Eをa、r0、もしくは、これらに適当な変数変換を施した変数でそれぞれ微分して0とおいた連立方程式を解いて求められるが、これらは簡単な形の方程式にはならず、容易に解くことができない。
【００６６】
このような問題は、非線形の最小二乗問題として知られており、一般的な解法としては、「最適化の手法」（非特許文献２）などに開示されたニュートン法を適用することが考えられる。ニュートン法を直接この問題に適用すると、a、r0に適当な初期値を設定し、式(38)に従ってa、r0を更新し、収束した時点でのa、r0を解とするという手順となる。但し、式(38)において、aold、r0,oldは更新前のパラメータ、a、r0は更新後のパラメータである。
【数５４】
【００６７】
しかしながら、式(36)を式(38)に代入して展開すると、この更新式は非常に複雑なものとなるため、計算に大きな演算量を要するという問題が生じる。またさらに、式(38)では、パラメータa、r0を同時に更新することになるが、一般的に複数のパラメータを繰返し演算によって同時に推定する場合は収束性が悪くなる傾向があり、適切な収束解が得られにくいという問題がある。
【００６８】
そこで、簡便な計算によって安定した解を得るために、最小化する評価式の特徴を利用し、その評価式に適した解法を導出することが考えられる。例えば、特開2006-350948号公報（特許文献１）には、主に有理多項式を対象にしたいくつかの解法が開示されている。しかしながら、これらの有理多項式に特化した解法は、式(3)のEの評価式を最小化する問題には直接適用することができない。そこで、式(3)のEを最小化する問題に対しては、この評価式に適した解法を新たに導出する必要がある。
【００６９】
本実施例の逆二乗特性解析装置は、上記のような技術的検討を経て、式(3)のEを最小化する減衰曲線の推定に適した構成を示したものである。以下では、この点について説明する。
【００７０】
式(3)のEを最小化する未知パラメータa、r0の算出方法を導出するため、式(3)から導出された式(36)を、さらに式(39)のように変形する。但し、ここでは式(40)の変数変換を用いている。式(39)のEをA、r0でそれぞれ微分すると、式(41)、式(42)が得られる。式(41)を0とおいて、これをAについて解くと式(43)となり、これを式(42)に代入すると、Aを消去することができて、式(44)が得られる。また、式(40)と式(43)を用いると、aについての算出式である式(45)が得られる。
【数５５】
【数５６】
【数５７】
【数５８】
【数５９】
【数６０】
【数６１】
【００７１】
以上の式変形を行うことにより、式(44)に含まれる未知パラメータはr0のみとなるので、r0に関しては式(44)を利用して１変数の最小二乗問題として求めることが可能となる。また、aに関しては、r0が求められた後、式(45)にr0を代入して直接求めることができる。
【００７２】
そこで、式(44)を基にしたニュートン法をr0の導出に適用すると、r0の更新式は式(46)となり、これに従って繰返し演算を行うことでr0が算出される。但し、∂2E/∂r02は、式(44)をさらにr0で微分することで求められ、式(47)で表される。また、r0の初期値としては、r0が音響中心の補正量であることを考えれば、r0は0に近い値となることが予想されるので、0を設定するのが妥当である。
【数６２】
【数６３】
【００７３】
ここで、r0を導出するための繰返し演算で用いられる式(46)、式(44)、式(47)、および、aを導出するための式(45)は、逆二乗特性解析装置１１の減衰曲線推定部１２において減衰曲線の推定に用いられる式(4)、式(5)、式(6)、および、式(7)と一致していることがわかる。
【００７４】
つまり、本実施例の逆二乗特性解析装置における減衰曲線の推定方法は、上記の技術的検討を経て導出されたものであり、このような構成にすることにより、減衰曲線を適切に推定することが可能となる。また、繰返し演算にて推定されるパラメータが１つのみの構成となっているため、比較的簡便な演算で処理を行うことが可能であり、適切な収束解が安定的に得られやすいという効果が得られる。
【００７５】
なお、減衰曲線を推定する際に最小化する差分の尺度と、測定値と推定値の乖離量を評価する時の尺度を一致させることの効果を実証するため、以下では、模擬的に作成した入力データを用いて、前提技術の逆二乗特性解析装置と本実施例の逆二乗特性解析装置の解析結果を比較する。
【００７６】
ここで使用する入力データは、測定点数がN=46で、各測定点の音源からの距離riはri=0.5+0.1(i-1) mに従って算出したものであり、音圧値piは式(48)に従って算出される音圧レベルLpiを音圧に変換したものである。なお、式(48)のα(r)は微小な誤差を示しており、|α(r)|の最大値は、r=4.5 mにおいて約1.52 dBである。すなわち、式(48)は、音源から0.5m〜2.5mの範囲では、式(1)においてa=1000、r0=0としたモデル式に対応する理想的な減衰曲線を示しており、2.5m〜5.0mの範囲では、その理想的な減衰曲線に微小な誤差が付加された曲線を示している。
【００７７】
なお、第46番目の測定点は音源から5.0mの位置となるので、検定範囲は音源から直線経路に沿って5.0mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【数６４】
【００７８】
表３は、前提技術の逆二乗特性解析装置の解析結果を示す表である。表３を参照して、第２列と第３列が、入力データとなる各測定点の音源からの距離ri、および、音圧の測定値piである。これらの入力データから前提技術に係る逆二乗特性解析装置において推定される未知パラメータa、r0は、a=888.79、r0=0.208608である。表３の第５列、第６列は、各測定点における音圧レベルの測定値Lpiと推定値Lp(ri)を示しており、第７列はそれらの差分δiを示している。
【表３】
【００７９】
図７は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表３の第５列の値）を丸い点で、推定曲線（式(19)にa=888.79、r0=0.208608を代入した曲線）を実線で示してある。
【数６５】
【００８０】
表３の第７列と図７を参照すると、前提技術の逆二乗特性解析装置の解析結果では、音源に近い測定点において測定値と推定値の乖離が大きくなっており、|δi|の最大値はr=0.5mの位置で約3.67dBとなっている。検定対象としている中心周波数500Hzの1/3オクターブ帯域の許容限界値δlimは2.5dBであるので、検定結果は不適合となる。
【００８１】
しかしながら、解析に用いた入力データは、音源から遠い範囲にのみ誤差が付加されており、音源に近い範囲では理想的な減衰曲線に従うものである。また、理想的な減衰曲線からの乖離量も最大1.52dB程度であり、本来、検定結果は適合となるべきものである。従って、前提技術の逆二乗特性解析装置では、適切な減衰曲線が推定されず、その結果として、誤った検定結果が出力される可能性があることが示された。
【００８２】
これに対し、表４は、本実施例の逆二乗特性解析装置の解析結果を示す表である。表４を参照して、第２列と第３列が、入力データとなる各測定点の音源からの距離ri、および、音圧の測定値piである。これらの入力データから本実施例の逆二乗特性解析装置において推定される未知パラメータa、r0は、a=960.208、r0=0.0389398である。表４の第５列、第６列は、各測定点における音圧レベルの測定値Lpiと推定値Lp(ri)を示しており、第７列はそれらの差分δiを示している。
【表４】
【００８３】
図８は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表４の第５列の値）を丸い点で、推定曲線（式(19)にa=960.208、r0=0.0389398を代入した曲線）を実線で示してある。
【００８４】
表４の第７列と図８を参照すると、本実施例の逆二乗特性解析装置の解析結果では、全ての測定点において測定値と推定値の乖離が一様に小さくなっており、|δi|の最大値はr=4.5mの位置で約1.24dBとなっている。検定対象としている中心周波数500Hzの1/3オクターブ帯域の許容限界値δlimは2.5dBであるので、検定結果は適合となる。
【００８５】
この結果は、解析に用いた入力データが、r=4.5mの位置で最大1.52dBの誤差を付加したものであることとよく一致しており、本実施例の逆二乗特性解析装置では、このような入力データに対しても、適切な減衰曲線が推定されて、正しい検定結果が出力されることが示された。
【００８６】
以上のことから、本実施例によれば、式(1)にてモデル化され、式(3)の評価式を最小にする減衰曲線の推定を、繰返し演算によって推定されるパラメータ数が少なく、かつ、簡素化された推定式にて行うため、少ない演算量にて計算が可能となり、さらに、適切な推定結果を安定的に得ることが可能となる。
【実施例２】
【００８７】
以下では、本発明の逆二乗特性解析装置の別の実施例を説明する。なお、逆二乗則検定を行う手順については、前提技術の逆二乗特性解析装置に代えて、本発明の逆二乗特性解析装置を使用すること以外、前提技術で説明した手順と全て同じであるので説明を省略し、逆二乗特性解析装置の構成と解析動作のみを説明する。
【００８８】
＜本実施例の逆二乗特性解析装置の構成＞
以下では、逆二乗則検定において使用される本実施例の逆二乗特性解析装置の構成を説明する。この逆二乗特性解析装置は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして受け、その距離減衰特性が反比例曲線から乖離する量を算出し、その乖離量に基づいて自由空間としての適合／不適合の解析結果を出力するものである。
【００８９】
図９は、本発明の逆二乗特性解析装置を説明するブロック構成図である。図９を参照して、逆二乗特性解析装置２１は、測定値から反比例曲線を推定する減衰曲線推定部２２と、推定された反比例曲線と測定値との乖離量を算出する差分算出部２３と、算出された乖離量から半無響室の自由空間特性を判定する特性評価部２４を備える。
【００９０】
減衰曲線推定部２２は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、式(1)でモデル化された反比例曲線の推定を行う。式(1)は、音圧pが距離rに反比例する関係を示した式であるが、未知のパラメータa、r0を含んだ形でモデル化されている。減衰曲線推定部２２では、式(8)、式(9)に従って未知パラメータa、r0を算出し、反比例曲線の推定を完了する。但し、wi (i=1〜N)は式(10)に従ってpiから算出される量である。また、qi(i=1〜N)は第i番目の測定点における逆数音圧であり、式(11)に従ってpiから算出される。なお、式(1)のモデル式はISO 3745に基づくものであるが、式(8)、式(9)、式(10)、式(11)に基づく未知パラメータa、r0の推定方法は、本実施例の逆二乗特性解析装置の特徴となる部分である。
【数６６】
【数６７】
【数６８】
【数６９】
【数７０】
【００９１】
差分算出部２３は、推定されたパラメータa、r0、および、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、各測定点における音圧レベルの測定値Lpi (i=1〜N)、音圧レベルの推定値Lp(ri) (i=1〜N)、および、それらの差分δi (i=1〜N)を算出する。ここで、Lpi およびLp(ri)の算出は、それぞれ、式(15)、式(16)に従って行われ、δiの算出は、式(17)に従って行われる。なお、式(15)は、式(2)に基づいて得られるpiとLpiの関であり、式(16)は、式(1)と式(2)に基づいて得られるa、r0、riとp(ri)の関係式である。
【数７１】
【数７２】
【数７３】
【００９２】
特性評価部２４は、差分δi (i=1〜N)を入力データとして、これらのδiに対して式(18)の関係が成立するか否かの判定を行い、全ての測定点、すなわち、i=1〜Nの全てのiにおいて式(18)が成立する場合には適合、そうでない場合には不適合の結果を出力する。なお、δlimは許容限界値であり、式(18)は差分δi の絶対値がこの許容限界値以内に収まることを示す関係式である。δlim の値は、ISO 3745では検定対象とする周波数帯域毎に定められており、ここでもそれ従う。
【数７４】
【００９３】
＜本実施例の逆二乗特性解析装置の動作＞
以下では、前提技術で例示したのと同様に、測定点数が26個（N=26）、各測定点の音圧値pi、距離riが表５の第３列と第２列に示されたデータである場合を例に、本実施例の逆二乗特性解析装置で行われる解析動作を説明する。なお、距離riはri=0.5+0.1(i-1) m の関係があるので、第26番目の測定点は音源から3.0mの位置となり、検定範囲は音源から直線経路に沿って3.0mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【表５】
【００９４】
逆二乗特性解析装置２１において、各測定点における音圧値pi、音源からの距離riが入力されると、減衰曲線推定部２２では、式(8)、式(9)、式(10)、式(11)に従って未知パラメータa、r0の算出が行われる。このとき、表５の第３列と第２列に示されたpi、riの入力データに対しては、a=1039.16、r0=-0.001743が算出される。
【００９５】
未知パラメータa、r0が推定されると、差分算出部２３では、式(15)と式(16)に従い、音圧レベルの測定値Lpiおよび推定値Lp(ri)が算出され、さらに、式(17)に従い、その差分δiが算出される。このときのLpi、Lp(ri)、δiは、表５の第５列、第６列、第７列にそれぞれ記された値となる。なお、図１０は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表５の第５列の値）を丸い点で、推定曲線（式(19)にa=1039.16、r0=-0.001743を代入した曲線）を実線で示してある。
【数７５】
【００９６】
差分δiが算出されると、特性評価部２４では、式(18)に従い、それらの値が許容限界値以内に収まるか否かが判定される。検定対象としている中心周波数500Hzの1/3オクターブ帯域では、許容値δlimは2.5dBと定められており、差分δiの算出結果である表５の第７列の値に対しては全て式(18)が成立するため、適合の解析結果が出力される。この解析結果は、中心周波数が500Hzの1/3オクターブ帯域において、直線経路に沿って音源から3.0mの位置までの範囲が自由音場とみなされることを示している。
【００９７】
以上によって、本実施例の逆二乗特性解析装置の構成、および、各測定点の音圧値pi (i=1〜N)と距離ri (i=1〜N)が入力されてから、自由空間としての適合／不適合を示す結果が出力されるまでの解析動作が説明された。
【００９８】
＜本実施例の逆二乗特性解析装置の特徴と効果＞
以下では、本実施例の逆二乗特性解析装置の特徴と効果について詳しく説明する。
実施例１において明らかにしたように、前提技術の逆二乗特性解析装置において推定される減衰曲線は、式(21)のE、すなわち、逆数音圧の差分の二乗和を最小化するものであり、測定値と推定値の乖離量を評価する尺度である音圧レベルとは一致していないため、適切な検定結果を得ることができない。
【００９９】
この点を考慮すると、減衰曲線推定部において、逆数音圧における差の二乗和ではなく、音圧レベルにおける差の二乗和、すなわち、式(3)のEを最小にするよう減衰曲線の推定を行うことが有効であると考えられるが、簡便な計算によって安定した解が得られるように、そのような推定を行うことは困難である。
【数７６】
【０１００】
本実施例の逆二乗特性解析装置は、上記のような技術的検討を経て、式(3)のEを最小化する減衰曲線の推定に適した別の構成を示したものである。以下では、この点について説明する。
【０１０１】
式(3)のEを最小化する未知パラメータa、r0の算出方法を導出するため、逆数音圧における測定値と推定値の差と、音圧レベルにおける測定値と推定値の差の関係について考えると、測定値と推定値の差が小さい場合には近似的に式(49)の関係が成り立つ。但し、音圧レベルLpの逆数音圧qによる微分dLp/dqには式(34)の関係式を用い、さらにその式の中のpを測定値piで近似している。この式(49)を式(3)に代入すると式(50)が得られる。但し、wiは式(51)に従い、piから算出される量である。
【数７７】
【数７８】
【数７９】
【数８０】
【０１０２】
推定値q(ri)は、式(12)、式(1)に基づいて式(22)のように表されるので、これを式(50)に代入すると、Eは式(52)のように式変形される。
【数８１】
【数８２】
【０１０３】
ここで、式(52)のEを最小にする未知パラメータa、r0は、Eをa、r0でそれぞれ微分して0とおいた連立方程式を解いて求められるが、この連立方程式は複雑で容易に解くことができない。そこで、式(25)の変数変換を用いて式(52)を式(53)のように式変形し、式(53)のEを最小化するA、Bを求めた後、式(25)の関係を用いてA、Bからa、r0を求める。
【数８３】
【数８４】
【０１０４】
式(53)のEを最小化するA、Bを求めるため、式(53)をA、Bでそれぞれ微分すると式(54)が得られ、これらを0とおいてA、Bについて整理すると式(55)の連立方程式が得られる。式(55)はA、Bの線形方程式として容易に式(56)のように式変形することができ、さらにこれを式(25)の関係式を用いてa、r0について解くと、式(57)、式(58)が得られる。これらはa、r0を直接算出できる式となっているので、以上の式変形を行うことによって、式(3)のEを最小化するa、r0の算出が可能となる。
【数８５】
【数８６】
【数８７】
【数８８】
【数８９】
【０１０５】
ここで、a、r0を導出するための式(57)、式(58)、式(51)は、逆二乗特性解析装置２１の減衰曲線推定部２２において減衰曲線の推定に用いられる式(8)、式(9)、式(10)と一致していることがわかる。
【０１０６】
つまり、本実施例の逆二乗特性解析装置における減衰曲線の推定方法は、上記の技術的検討を経て導出されたものであり、このような構成にすることにより、減衰曲線を適切に推定することが可能になる。また、繰返し演算を行う必要がない構成となっているため、簡便な演算で処理を行うことが可能であり、解が確実に得られるという効果がある。
【０１０７】
なお、減衰曲線を推定する際に最小化する差分の尺度と、測定値と推定値の乖離量を評価する時の尺度を一致させることの効果を実証するため、以下では、模擬的に作成した入力データを用いて、前提技術の逆二乗特性解析装置と本実施例の逆二乗特性解析装置の解析結果を比較する。
【０１０８】
ここで使用する入力データは、実施例１と同様に、測定点数がN=46で、各測定点の音源からの距離riはri=0.5+0.1(i-1) mに従って算出したものであり、音圧値piは式(48)に従って算出される音圧レベルLpiを音圧に変換したものである。なお、式(48)のα(r)は微小な誤差を示しており、|α(r)|の最大値は、r=4.5 mにおいて約1.52 dBである。すなわち、式(48)は、音源から0.5m〜2.5mの範囲では、式(1)においてa=1000、r0=0としたモデル式に対応する理想的な減衰曲線を示しており、2.5m〜5.0mの範囲では、その減衰曲線に微小な誤差が付加された曲線を示している。
【０１０９】
なお、第46番目の測定点は音源から5.0mの位置となるので、検定範囲は音源から直線経路に沿って5.0mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【数９０】
【０１１０】
実施例１において説明したように、この入力データに対する前提技術の逆二乗特性解析装置の解析結果は、音源に近い測定点において、測定値と推定値の乖離が大きくなり、|δi|の最大値はr=0.5mの位置で約3.67dBとなって、検定結果は不適合となった。つまり、前提技術に係る逆二乗特性解析装置では、適切な減衰曲線が推定されず、その結果として、誤った検定結果が出力される可能性があることが示された。
【０１１１】
これに対し、表６は、本実施例の逆二乗特性解析装置の解析結果を示す表である。表６を参照して、第２列と第３列が、入力データとなる各測定点の音源からの距離ri、および、音圧の測定値piである。これらの入力データから本実施例の逆二乗特性解析装置において推定される未知パラメータa、r0は、a=969.883、r0=0.0309326である。表６の第５列、第６列は、各測定点における音圧レベルの測定値Lpiと推定値Lp(ri)を示しており、第７列はそれらの差分δiを示している。
【表６】
【０１１２】
図１１は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表６の第５列の値）を丸い点で、推定曲線（式(19)にa=969.883、r0=0.0309326を代入した曲線）を実線で示してある。
【数９１】
【０１１３】
表６の第７列と図１１を参照すると、本実施例の逆二乗特性解析装置の解析結果では、全ての測定点において測定値と推定値の乖離が一様に小さくなっており、|δi|の最大値はr=4.5mの位置で約1.32dBとなっている。検定対象としている中心周波数500Hzの1/3オクターブ帯域の許容限界値δlimは2.5dBであるので、検定結果は適合となる。
【０１１４】
この結果は、解析に用いた入力データが、r=4.5mの位置で最大1.52dBの誤差を付加したものであることとよく一致しており、本実施例の逆二乗特性解析装置では、このような入力データに対しても、適切な減衰曲線が推定されて、正しい検定結果が出力されることが示された。
【０１１５】
以上のことから、本実施例によっても、式(1)にてモデル化され、式(3)の評価式を最小にする減衰曲線の推定を、簡素化された推定式にて行うため、少ない演算量にて計算が可能となり、さらに、適切な推定結果を安定的に得ることが可能となる。
【符号の説明】
【０１１６】
１、１１、２１…逆二乗特性解析装置
２、１２、２２…減衰曲線推定部
３、１３、２３…差分算出部
４、１４、２４…特性評価部
Ｓ１…処理ステップ（パラメータr0への初期値設定）
Ｓ２…処理ステップ（パラメータr0の更新）
Ｓ３…処理ステップ（収束判定）
【技術分野】
【０００１】
本発明は、半無響室および無響室の自由空間特性を評価する解析装置及び方法に関する。
【背景技術】
【０００２】
半無響室および無響室は、自由空間（音の反射のない空間）を模した試験室であり、音の測定環境として広く利用されている。半無響室および無響室が自由空間にどの程度近い特性を有しているかという点は、その試験室における測定精度に影響を及ぼすため、その特性の適切な評価方法の確立は極めて重要な技術課題である。半無響室および無響室の自由空間としての特性を評価する方法としては、ISO 3745（非特許文献１）の附属書Aに記された逆二乗則検定が一般的に利用されている。
【０００３】
＜逆二乗則検定の概要＞
以下では、逆二乗則検定の概要について説明する。
自由空間に置かれた無指向性の点音源から放射される音の強さ（音響インテンシティ）は、音源からの距離の二乗に反比例して減衰する（逆二乗則）。また、このときの音の強さ（音響インテンシティ）は音圧の二乗に比例するため、音圧は音源からの距離に反比例して減衰する。
【０００４】
従って、無響室内に無指向性の点音源を設置して、様々な測定点において音圧を測定した場合、その無響室が完全な自由空間であるならば、各測定点で測定される音圧は音源からの距離に反比例したものとなる。しかしながら、壁面からの反射が生じる場合には、測定される音圧は、重畳される反射波の度合いに応じて反比例曲線から乖離するようになる。そこで、その乖離量を算出し、その大きさによって自由空間としての特性評価を行うことが可能となる。
【０００５】
半無響室については、自由空間に床からの反射が重畳された特性（反射面上の自由空間）が理想特性となるため、その室内に無指向性の点音源を設置したとしても、床からの反射波との干渉により、音圧は必ずしも距離に反比例して減衰しない。しかし、音源を床面上に設置した場合には、床面上方へ半球状に放射される直接波と床面から鏡面反射され上方に半球状に広がる反射波が完全に重なり合い、両波面の干渉による減衰特性の変化が生じないため、自由空間と同様に扱うことが可能となる。すなわち、半無響室が反射面上の自由空間として完全な特性を持つならば、音源を床面上に置いた場合には、各測定点で測定される音圧は音源からの距離に反比例したものとなる。しかしながら、床面以外の壁面からの反射が生じる場合には、測定される音圧は、重畳される床面以外の反射波の度合いに応じて反比例曲線から乖離するようになる。つまり、半無響室の場合も、音源を床面上に置いて測定を行うことによって、無響室と同様に、反比例曲線からの音圧値の乖離量による特性評価が可能となる。
【０００６】
ISO 3745には、この基本原理に基づいて、点音源を設置した試験室内において測定した音圧データから理想的な減衰曲線を推定し、その減衰曲線と測定値との乖離量によって、その試験室の自由空間としての特性評価を行う方法が記載されている。
以上が、逆二乗則検定の概要である。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開2006-350948号公報
【非特許文献】
【０００８】
【非特許文献１】国際規格 ISO 3745:2003
【非特許文献２】茨木俊秀・福島雅夫、最適化の手法（共立出版）ISBN4-320-02664-0
【発明の概要】
【発明が解決しようとする課題】
【０００９】
ISO 3745に記載された方法に基づいて解析を行う従来の逆二乗特性解析装置では、測定した音圧値から適切な減衰曲線が推定されず、理想的な減衰特性からの測定値の乖離量が正確に算出されないことがあった。その場合、半無響室および無響室が良好な特性を有していても、検定結果は不適合となることがあり、合理的な検定結果が得られないという問題があった。
【００１０】
本発明は、このような課題を解決するためになされたものであり、測定した音圧値から適切な減衰曲線の推定を行うことが可能な逆二乗特性解析装置を提供することを目的とする。またさらに、その減衰曲線の推定を少ない計算量によって行うこと、および、適切な推定結果を安定的に得ることが可能な逆二乗特性解析装置を提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記の目的を達成するために、本発明の第一の実施態様による逆二乗特性解析装置は、試験音場内の複数の測定点における音圧の測定値とそれぞれの測定点の位置とを入力データとして、音圧と位置の関係を反比例で近似した減衰曲線を推定する減衰曲線推定手段と、前記推定された減衰曲線と前記音圧の測定値との差を音圧レベルの尺度において算出する差分算出手段と、前記算出した差に基づいて前記試験音場の自由空間特性を評価する特性評価手段を備える逆二乗特性解析装置であって、前記減衰曲線推定手段は、前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように前記減衰曲線の推定を行うことを特徴とする。
【００１２】
本発明の第二の実施態様による逆二乗特性解析装置は、試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、前記pi、ri (i=1〜N)を入力データとして、式(1)で表される減衰曲線の未知パラメータa、r0を推定する減衰曲線推定手段と、前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出する差分算出手段と、前記算出された差が所定の許容値に収まるか否かによって前記試験音場の自由空間特性を評価する特性評価手段を備える逆二乗特性解析装置であって、前記減衰曲線推定手段は、式(3)を最小にするように前記未知パラメータa、r0の推定を行うことを特徴とする。
【数１】
【数２】
【数３】
【００１３】
本発明の第三の実施態様による逆二乗特性解析装置は、上記第二の実施態様による逆二乗特性解析装置であって、前記減衰曲線推定手段は、前記未知パラメータr0に初期値を設定し、式(4)、式(5)及び式(6)に従って前記未知パラメータr0を更新する処理を前記未知パラメータr0が収束するまで繰返し行って前記未知パラメータr0の収束解を算出し、前記算出した未知パラメータr0の収束解を用いて、式(7)に従って前記未知パラメータaを算出することを特徴とする。
【数４】
【数５】
【数６】
【数７】
【００１４】
本発明の第四の実施態様による逆二乗特性解析装置は、上記第二の実施態様による逆二乗特性解析装置であって、前記減衰曲線推定手段は、式(8)、式(9)、式(10)及び式(11)に従って、前記未知パラメータa、r0の推定を行うことを特徴とする。
【数８】
【数９】
【数１０】
【数１１】
【００１５】
本発明の第五の実施態様による試験音場内の逆二乗特性の評価方法は、試験音場内の複数の測定点における音圧の測定値とそれぞれの測定点の位置とを入力データとして、音圧と位置の関係を反比例で近似した減衰曲線を推定するステップと、前記推定された減衰曲線と前記音圧の測定値との差を音圧レベルの尺度において算出するステップと、前記算出した差に基づいて前記試験音場の自由空間特性を評価するステップからなり、前記減衰曲線を推定するステップでは、前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように前記減衰曲線の推定を行うことを特徴とする。
【００１６】
本発明の第六の実施態様による試験音場内の逆二乗特性の評価方法は、試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、前記pi、ri (i=1〜N)を入力データとして、式(1)で表される減衰曲線の未知パラメータa、r0を推定するステップと、前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出するステップと、前記算出された差が所定の許容値に収まるか否かによって前記試験音場の自由空間特性を評価するステップからなり、前記減衰曲線を推定するステップは、式(3)を最小にするように前記未知パラメータa、r0の推定を行うことを特徴とする。
【数１２】
【数１３】
【数１４】
【発明の効果】
【００１７】
本発明の逆二乗特性解析装置によれば、減衰曲線推定手段において音圧の測定値と減衰曲線との差を最小化するように減衰曲線の推定を行うときの尺度が、差分算出手段および特性評価手段において音圧の測定値と減衰曲線の差を評価するときの尺度と一致するため、適切な減衰曲線を推定することが可能となる。その結果として、測定値の減衰曲線からの乖離量が精度よく算出されるため、合理的な検定結果を得ることが可能となる。
【図面の簡単な説明】
【００１８】
【図１】半無響室における逆二乗則検定の様子を説明する模式図である。
【図２】前提技術の逆二乗特性解析装置を説明するブロック構成図である。
【図３】前提技術の逆二乗特性解析装置による解析結果の一例を示すグラフである。
【図４】実施例１の逆二乗特性解析装置を説明するブロック構成図である。
【図５】実施例１の逆二乗特性解析装置の減衰曲線推定部において、未知パラメータr0を算出する処理手順を説明する処理フロー図である。
【図６】実施例１の逆二乗特性解析装置による解析結果の一例を示すグラフである。
【図７】前提技術の逆二乗特性解析装置による解析結果の別の一例を示すグラフである。
【図８】実施例１の逆二乗特性解析装置による解析結果の別の一例を示すグラフである。
【図９】実施例２の逆二乗特性解析装置を説明するブロック構成図である。
【図１０】実施例２の逆二乗特性解析装置による解析結果の一例を示すグラフである。
【図１１】実施例２の逆二乗特性解析装置による解析結果の別の一例を示すグラフである。
【発明を実施するための形態】
【００１９】
まず、以下では、本発明の実施の形態に係る前提技術を説明する。ここでは、ISO 3745に記載された方法を例とした逆二乗則検定の実施手順、および、その中で用いられる逆二乗特性解析装置の構成と動作について説明する。
【００２０】
＜逆二乗則検定の実施手順＞
以下では、ISO 3745に従って逆二乗則検定を行う手順について説明する。なお、ISO 3745は半無響室と無響室を検定対象としているが、ここでは半無響室の検定を例に説明を行う。
【００２１】
図１は検定対象となる半無響室の模式図である。図１を参照して、逆二乗則検定では、試験音場である半無響室の床面の中央に発音体である音源（スピーカ）を設置し、ここから試験音場内（半無響室内）へ音の放射を行う。さらに、この音源を端点とする直線経路上にいくつかの測定点を設定し、それぞれの測定点において音圧を測定する。
【００２２】
直線経路は音源から半無響室上部のいずれかの角に向かうものとし、測定点については、直線経路上の音源から0.5mの位置を第1番目の測定点として、そこから順に直線経路に沿って0.1 m毎に合計N個の測定点が設定されるものとする。直線経路および測定点の位置については、ISO 3745にその決め方が定められており、上記の設定はこれに従うものである。また、このときの第i番目(i=1〜N)の測定点の音源からの距離をri (i=1〜N)、測定された音圧値をpi (i=1〜N)とする。ここで、距離ri の単位はm（メートル）、音圧piの単位はPa（パスカル）である。但し、ISO 3745では所定の周波数帯域（1/1もしくは1/3オクターブ帯域）毎に検定を行うことが定められており、上記のpiは検定対象とする周波数帯域の音圧値を表すものとする。
【００２３】
なお、上記の測定に用いる収音システム（マイクロホンなど）に求められる仕様や、検定を行う各帯域の周波数範囲、帯域分割フィルタの特性などは、いずれもISO 3745に定められており、ここでもその定めに従って測定を行う。また、音源（スピーカ）についても無指向性で点音源に近い特性が要求されるが、ISO 3745に具体的な要求仕様が定められており、それに従うものとする。このような測定は、規格に従い適切な測定装置を使用することで一般的に実施が可能である。
【００２４】
ところで、音圧は、圧力の単位であるパスカル（単位記号Pa）や、20μPaを基準として式(2)に従う対数尺度で示したデシベル（単位記号dB）など、異なる単位で表されることがある。これらは相互に変換することが可能であり、いずれの単位で表されていても本質的な違いはないが、説明の便宜のため、本明細書では、パスカルで表される量を音圧（量記号p）、式(2)に従いデシベルで表される量を音圧レベル（量記号Lp）と区別して呼ぶこととする。また、20μPaを基準音圧と定義し、記号p0で表すこととする。さらに、式(12)に従い、基準音圧20μPaで正規化された音圧の逆数を逆数音圧（量記号q）と定義する。
【数１５】
【数１６】
【００２５】
上記に従ってN個の測定点で音圧を測定した後、各測定点の音圧値pi(i=1〜N)と音源からの距離ri (i=1〜N)を後述する逆二乗特性解析装置に入力する。逆二乗特性解析装置は、入力された測定値の距離減衰特性が反比例曲線からどの程度乖離しているかを解析し、その乖離量に応じて、適合または不適合の結果を出力する。出力結果が適合である場合には、検定を行った周波数帯域において、音源から第N番目の測定点の位置までの範囲で良好な逆二乗特性が成立しており、その範囲は自由音場とみなされる。出力結果が不適合である場合は、検定を行った周波数帯域において、音源から第N番目の測定点の位置までの範囲は自由音場とみなされない。
以上の手順によって、ISO 3745に従う半無響室の逆二乗則検定が行われる。
【００２６】
＜前提技術に係る逆二乗特性解析装置の構成＞
以下では、前提技術に係る逆二乗特性解析装置の構成を説明する。この逆二乗特性解析装置は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして受け、その距離減衰特性が反比例曲線から乖離する量を算出し、その乖離量に基づいて自由空間としての適合／不適合の解析結果を出力するものである。なお、その解析方法は、ISO 3745に記載された方法に従うものである。
【００２７】
図２は、前提技術に係る逆二乗特性解析装置を説明するブロック構成図である。図２を参照して、逆二乗特性解析装置１は、測定値から反比例曲線を推定する減衰曲線推定部２と、推定された反比例曲線と測定値との乖離量を算出する差分算出部３と、算出された乖離量から半無響室の自由空間特性を判定する特性評価部４を備える。
【００２８】
減衰曲線推定部１は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、式(1)でモデル化された反比例曲線の推定を行う。式(1)は、音圧pが距離rに反比例する関係を示した式であるが、未知のパラメータa、r0を含んだ形でモデル化されている。減衰曲線推定部１では、式(13)、式(14)に従って未知パラメータa、r0を算出し、反比例曲線の推定を完了する。但し、qi(i=1〜N)は第i番目の測定点における逆数音圧であり、式(12)に従ってpiから算出される。なお、式(1)のモデル式、および、式(13)、式(14)の未知パラメータa、r0の推定式は、ISO 3745の記載に基づくものである。
【数１７】
【数１８】
【数１９】
【００２９】
差分算出部３は、推定された未知パラメータa、r0、および、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、各測定点における音圧レベルの測定値Lpi (i=1〜N)、音圧レベルの推定値Lp(ri) (i=1〜N)、および、それらの差分δi (i=1〜N)を算出する。ここで、Lpi およびLp(ri)の算出は、それぞれ、式(15)、式(16)に従って行われ、δiの算出は、式(17)に従って行われる。なお、式(15)は、式(2)に基づいて得られるpiとLpiの関係式であり、式(16)は、式(1)と式(2)に基づいて得られるa、r0、riとLp(ri)の関係式である。
【数２０】
【数２１】
【数２２】
【００３０】
特性評価部４は、差分δi (i=1〜N)を入力データとして、これらのδiに対して式(18)の関係が成立するか否かの判定を行い、全ての測定点、すなわち、i=1〜Nの全てのiにおいて式(18)が成立する場合には適合、そうでない場合には不適合の結果を出力する。なお、δlimは許容限界値であり、式(18)は差分δi の絶対値がこの許容限界値以内に収まることを示す関係式である。δlim の値は、ISO 3745では検定対象とする周波数帯域毎に定められており、ここでもそれに従う。
【数２３】
【００３１】
＜前提技術に係る逆二乗特性解析装置の動作＞
以下では、測定点数が26個（N=26）、各測定点の音圧値pi、距離riが、表１の第３列と第２列に示されたデータである場合を例に、前提技術に係る逆二乗特性解析装置で行われる解析動作を説明する。なお、距離riはri=0.5+0.1(i-1) m の関係があるので、第26番目の測定点は音源から3.0mの位置となり、検定範囲は音源から直線経路に沿って3.0 mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【表１】
【００３２】
逆二乗特性解析装置１において、各測定点における音圧値pi、音源からの距離riが入力されると、減衰曲線推定部２では、式(13)、式(14)に従って未知パラメータa、r0の算出が行われる。このとき、表１の第３列と第２列に示されたpi、riの入力データに対しては、a=1037.32、r0=-0.0105871が算出される。
【００３３】
未知パラメータa、r0が推定されると、差分算出部３では、式(15)と式(16)に従い、音圧レベルの測定値Lpiおよび推定値Lp(ri)が算出され、さらに、式(17)に従い、その差分δiが算出される。このときのLpi、Lp(ri)、δiは、表１の第５列、第６列、第７列にそれぞれ記された値となる。なお、図３は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表１の第５列の値）を丸い点で、推定曲線（式(19)にa=1037.32、r0=-0.0105871を代入した曲線）を実線で示してある。
【数２４】
【００３４】
差分δiが算出されると、特性評価部４では、式(18)に従い、それらの値が許容限界値以内に収まるか否かが判定される。検定対象としている中心周波数500Hzの1/3オクターブ帯域では、許容値δlimは2.5dBと定められており、差分δiの算出結果である表１の第７列の値に対しては全て式(18)が成立するため、適合の解析結果が出力される。この解析結果は、中心周波数が500Hzの1/3オクターブ帯域において、直線経路に沿って音源から3.0 mの位置までの範囲が自由音場とみなされることを示している。
【００３５】
以上によって、前提技術に係る逆二乗特性解析装置の構成、および、各測定点の音圧値pi (i=1〜N)と距離ri (i=1〜N)が入力されてから、自由空間としての適合／不適合を示す結果が出力されるまでの解析動作が説明された。なお、前提技術に係る逆二乗特性解析装置で用いられている解析方法はISO 3745に記載された方法に従うものである。
【００３６】
ここで、式(1)のモデル式の未知パラメータa、r0について説明する。aは音源からの音響出力の大きさを示すパラメータであり、減衰曲線が測定時の音量に依存するために設けられているものである。また、r0は音源の音響中心を補正するためのパラメータであり、実際の音響中心とriの算出時に原点として設定される音響中心とが必ずしも一致しないので、その補正のために設けられているものである。これらのパラメータは、一般的に、事前に知ることが難しいため、式(1)のように未知のパラメータとしてモデル化し、減衰曲線推定部において、測定値からこれらのパラメータを推定して反比例曲線を決定する構成となっている。なお、ここで重要な点は、推定されるべき曲線は、測定環境が完全な自由音場であったとした場合の理想的な減衰曲線であるということである。その理由は、測定値との乖離量を算出する際、そのような曲線を基準としなければ、自由空間としての評価を正しく行うことができないからである。しかし、そのような曲線を厳密に求めることはできないので、各測定点で測定された音圧値の近傍を通る曲線が適切な推定曲線であると仮定し、その曲線は式(13)、式(14)に従って求められるという考えに基づいて、前提技術に係る逆二乗特性解析装置は構成されている。
【００３７】
しかしながら、後述するように、前提技術に係る逆二乗特性解析装置において推定される減衰曲線は、各測定点で測定された音圧値の近傍を通る曲線とはならないことがあり、その結果、正しい検定結果が得られないことがある。このような前提技術に係る課題に対処するため、前提技術のさらなる改良として、本発明に係る実施の形態を以下に説明する。
【実施例１】
【００３８】
以下では、本発明の逆二乗特性解析装置の実施例を説明する。なお、逆二乗則検定を行う手順については、前提技術の逆二乗特性解析装置に代えて、本発明の逆二乗特性解析装置を使用すること以外、前提技術で説明した手順と全て同じであるので説明を省略し、逆二乗特性解析装置の構成と解析動作のみを説明する。
【００３９】
＜本実施例の逆二乗特性解析装置の構成＞
以下では、逆二乗則検定において使用される本実施例の逆二乗特性解析装置の構成を説明する。この逆二乗特性解析装置は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして受け、その距離減衰特性が反比例曲線から乖離する量を算出し、その乖離量に基づいて自由空間としての適合／不適合の解析結果を出力するものである。
【００４０】
図４は、本発明の逆二乗特性解析装置を説明するブロック構成図である。図４を参照して、逆二乗特性解析装置１１は、測定値から反比例曲線を推定する減衰曲線推定部１２と、推定された反比例曲線と測定値との乖離量を算出する差分算出部１３と、算出された乖離量から半無響室の自由空間特性を判定する特性評価部１４を備える。
【００４１】
減衰曲線推定部１１は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、式(1)でモデル化された反比例曲線の推定を行う。式(1)は、音圧pが距離rに反比例する関係を示した式であるが、未知のパラメータa、r0を含んだ形でモデル化されている。
【数２５】
【００４２】
未知パラメータr0は、図５の処理フローで示される繰返し演算によって算出される。図５を参照して、処理ステップＳ１は、r0に初期値を設定するステップであり、ここでは、r0=0が設定されるものとする。処理ステップＳ２は、r0の更新を行うステップであり、現時点でのr0をr0,oldとして、式(4)に従ってr0,oldから新しく更新されたr0を算出する。但し、式(4)の中の∂E/∂r0、∂2E/∂r02は、それぞれ、式(5)、式(6)に従い算出されるものである。処理ステップＳ３は、この繰返し演算が収束したか否かを判定するステップである。処理ステップＳ２において算出したr0とr0,oldの関係が、式(20)の関係を満たす場合には、演算が収束したと判断し、その時点でのr0を最終結果として算出を終える。式(20)を満たさない場合には、再びステップＳ２に戻り、r0を更新する演算を繰り返す。なお、εは演算精度などを考慮して定められる充分小さな正定数であるが、ここではε=10-12であるとする。
【数２６】
【数２７】
【数２８】
【数２９】
【００４３】
このような繰返し演算は、DSP（Digital Signal Processor）やRAM（Random Access Memory）によって構成される一般的な演算装置によって実施することが可能であり、減衰曲線推定部１２はこのような演算装置を備えている。
【００４４】
繰返し演算によって未知パラメータr0が算出されると、式(7)に従って未知パラメータaが算出され、式(1)でモデル化された減衰曲線の推定が完了する。
【数３０】
【００４５】
なお、式(1)のモデル式はISO 3745に基づくものであるが、図５の処理フロー、および、式(4)、式(5)、式(6)、式(20)、式(7)に基づく未知パラメータa、r0の推定方法は、本実施例の逆二乗特性解析装置の特徴となる部分である。
【００４６】
差分算出部１３は、推定された未知パラメータa、r0および、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、各測定点における音圧レベルの測定値Lpi (i=1〜N)、音圧レベルの推定値Lp(ri) (i=1〜N)、および、それらの差分δi (i=1〜N)を算出する。ここで、Lpi およびLp(ri)の算出は、それぞれ、式(15)、式(16)に従って行われ、δiの算出は、式(17)に従って行われる。なお、式(15)は、式(2)に基づいて得られるpiとLpiの関係式であり、式(16)は、式(1)と式(2)に基づいて得られるa、r0、riとp(ri)の関係式である。
【数３１】
【数３２】
【数３３】
【００４７】
特性評価部１４は、差分δi (i=1〜N)を入力データとして、これらのδiに対して式(18)の関係が成立するか否かの判定を行い、全ての測定点、すなわち、i=1〜Nの全てのiにおいて式(18)が成立する場合には適合、そうでない場合には不適合の結果を出力する。なお、δlimは許容限界値であり、式(18)は差分δi の絶対値がこの許容限界値以内に収まることを示す関係式である。δlim の値は、ISO 3745では検定対象とする周波数帯域毎に定められており、ここでもそれ従う。
【数３４】
【００４８】
＜本実施例の逆二乗特性解析装置の動作＞
以下では、前提技術で例示したのと同様に、測定点数が26個（N=26）、各測定点の音圧値pi、距離riが表２の第３列と第２列に示されたデータである場合を例に、本実施例の逆二乗特性解析装置で行われる解析動作を説明する。なお、距離riはri=0.5+0.1(i-1) m の関係があるので、第26番目の測定点は音源から3.0 mの位置となり、検定範囲は音源から直線経路に沿って3.0 mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【表２】
【００４９】
逆二乗特性解析装置１１において、各測定点における音圧値pi、音源からの距離riが入力されると、減衰曲線推定部１２では、図５に示された処理フロー、および、式(4)、式(5)、式(6)、式(20)に従って未知パラメータr0の算出が行われ、さらに、式(7)に従って未知パラメータaの算出が行われる。このとき、表２の第３列と第２列に示されたpi、riの入力データに対しては、a=1029.45、r0=0.0043552が算出される。
【００５０】
未知パラメータa、r0が推定されると、差分算出部１３では、式(15)と式(16)に従い、音圧レベルの測定値Lpiおよび推定値Lp(ri)が算出され、さらに、式(17)に従い、その差分δiが算出される。このときのLpi、Lp(ri)、δiは、表２の第５列、第６列、第７列にそれぞれ記された値となる。なお、図６は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表２の第５列の値）を丸い点で、推定曲線（式(19)にa=1029.45、r0=0.0043552を代入した曲線）を実線で示してある。
【数３５】
【００５１】
差分δiが算出されると、特性評価部１４では、式(18)に従い、それらの値が許容限界値以内に収まるか否かが判定される。検定対象としている中心周波数500Hzの1/3オクターブ帯域では、許容値δlimは2.5dBと定められており、差分δiの算出結果である表２の第７列の値に対しては全て式(18)が成立するため、適合の解析結果が出力される。この解析結果は、中心周波数が500Hzの1/3オクターブ帯域において、直線経路に沿って音源から3.0 mの位置までの範囲が自由音場とみなされることを示している。
【００５２】
以上によって、本実施例の逆二乗特性解析装置の構成、および、各測定点の音圧値pi (i=1〜N)と距離ri (i=1〜N)が入力されてから、自由空間としての適合／不適合を示す結果が出力されるまでの解析動作が説明された。
【００５３】
＜本実施例の逆二乗特性解析装置の特徴と効果＞
以下では、本実施例の逆二乗特性解析装置の特徴と効果について詳しく説明する。
最初に、前提技術において未知パラメータa、r0の推定に用いられる式(13)、式(14)の技術的意味を明らかにする。
【００５４】
各測定点における逆数音圧の測定値qiと推定値q(ri)の差の二乗和をEとすると、Eは式(21)のように表される。推定値q(ri)は、式(12)と式(1)に基づいて式(22)のように表されるので、これを式(21)に代入すると、Eは式(23)のように式変形される。
【数３６】
【数３７】
【数３８】
【００５５】
ここで、式(23)のEを最小にする未知パラメータa、r0は、Eをa、r0でそれぞれ微分して0とおいた連立方程式、すなわち、式(24)を解いて求められる。
【数３９】
【００５６】
ところが、式(24)はa、r0について複雑な方程式となり、容易に解くことができない。そこで、式(25)の変数変換を用いて式(23)を式(26)のように式変形し、式(26)のEを最小化するA、Bを求めた後、式(25)の関係を用いてA、Bからa、r0を求める。
【数４０】
【数４１】
【００５７】
式(26)のEを最小化するA、Bを求めるため、式(26)をA、Bで微分すると式(27)が得られ、これらを0とおいてA、Bについて整理すると式(28)の連立方程式が得られる。式(28)はA、Bの線形方程式として容易に式(29)のように式変形することができ、さらにこれを式(25)の関係式を用いてa、r0について解くと式(30)、式(31)となって、これらが式(21)のEを最小化するa、r0の算出式となる。
【数４２】
【数４３】
【数４４】
【数４５】
【数４６】
【００５８】
式(30)、式(31)は、前提技術においてa、r0の推定で用いた式(13)、式(14)に一致しており、以上の検討によって、前提技術の逆二乗特性解析装置において推定される減衰曲線は、式(21)のE、すなわち、逆数音圧の差分の二乗和を最小化したものであることが明らかになった。
【００５９】
ところで、差分算出部および特性評価部では、測定値と推定値の差分を音圧レベルにおいて算出するので、前提技術の逆二乗特性解析装置で推定される減衰曲線の評価尺度である逆数音圧とは一致していない。以下では、このことが解析結果に及ぼす影響について説明する。
【００６０】
逆数音圧qにおける微小な変化をΔq、音圧レベルLpにおける微小な変化をΔLpとすると、両者の間には近似的に式(32)の関係が成立する。式(2)、式(12)を用いて、音圧レベルLpと逆数音圧qの関係は式(33)と表されるので、Lpのqによる微分dLp/dqは式(34)となり、これを式(32)に代入するとΔqとΔLpの間には式(35)の関係が成立する。
【数４７】
【数４８】
【数４９】
【数５０】
【００６１】
式(35)は、ΔqとΔLpの間に1/pに比例する重み係数が掛かることを示している。つまり、逆数音圧における測定値と推定値の差分を、全ての測定点において一様に最小化した場合、音圧レベルの尺度においては、測定点毎に1/pに比例した重みをかけて最小化したこととなる。1/pは音圧が小さいほど大きな値となるので、前提技術の逆二乗特性解析装置における減衰曲線の推定方法は、音圧レベルで見た場合、音圧が小さい測定点における差分を重点的に小さくするような推定をしていることに相当する。一般的に、音源に近い測定点ほど音圧は大きいので、このような推定がなされた場合、特に音源付近において測定値と大きく乖離した減衰曲線が推定される可能性がある。このような各測定点の近傍を一様には通らない減衰曲線は、測定環境が完全な自由音場であった場合の理想的な減衰曲線を推定しているとは言えず、それを基準として算出された乖離量も自由空間としての特性を正しく反映させたものとはならない。
【００６２】
以上の検討によって、前提技術の逆二乗特性解析装置では、減衰曲線を推定する際に最小化する差分の尺度と、測定値と推定値の乖離量を評価する時の尺度の違いが考慮されておらず、この点が誤った検定結果をもたらす原因となり得ることが示された。
【００６３】
上記の検討結果を考慮すると、減衰曲線推定部においては、逆数音圧における差の二乗和ではなく、音圧レベルにおける差の二乗和を最小にするよう減衰曲線の推定を行うことが有効であると考えられる。これは言い換えると、音圧レベルにおける測定値Lpiと推定値Lp(ri)の差の二乗和、すなわち、式(3)のEを最小化する未知パラメータa、r0を算出することである。
【数５１】
【００６４】
式(2)、式(1)に基づいて得られる式(37)の関係を用いて、式(3)を未知パラメータa、r0を含む形に式変形すると、式(36)が得られる。
【数５２】
【数５３】
【００６５】
ここで、式(36)のEを最小にする未知パラメータa、r0は、Eをa、r0、もしくは、これらに適当な変数変換を施した変数でそれぞれ微分して0とおいた連立方程式を解いて求められるが、これらは簡単な形の方程式にはならず、容易に解くことができない。
【００６６】
このような問題は、非線形の最小二乗問題として知られており、一般的な解法としては、「最適化の手法」（非特許文献２）などに開示されたニュートン法を適用することが考えられる。ニュートン法を直接この問題に適用すると、a、r0に適当な初期値を設定し、式(38)に従ってa、r0を更新し、収束した時点でのa、r0を解とするという手順となる。但し、式(38)において、aold、r0,oldは更新前のパラメータ、a、r0は更新後のパラメータである。
【数５４】
【００６７】
しかしながら、式(36)を式(38)に代入して展開すると、この更新式は非常に複雑なものとなるため、計算に大きな演算量を要するという問題が生じる。またさらに、式(38)では、パラメータa、r0を同時に更新することになるが、一般的に複数のパラメータを繰返し演算によって同時に推定する場合は収束性が悪くなる傾向があり、適切な収束解が得られにくいという問題がある。
【００６８】
そこで、簡便な計算によって安定した解を得るために、最小化する評価式の特徴を利用し、その評価式に適した解法を導出することが考えられる。例えば、特開2006-350948号公報（特許文献１）には、主に有理多項式を対象にしたいくつかの解法が開示されている。しかしながら、これらの有理多項式に特化した解法は、式(3)のEの評価式を最小化する問題には直接適用することができない。そこで、式(3)のEを最小化する問題に対しては、この評価式に適した解法を新たに導出する必要がある。
【００６９】
本実施例の逆二乗特性解析装置は、上記のような技術的検討を経て、式(3)のEを最小化する減衰曲線の推定に適した構成を示したものである。以下では、この点について説明する。
【００７０】
式(3)のEを最小化する未知パラメータa、r0の算出方法を導出するため、式(3)から導出された式(36)を、さらに式(39)のように変形する。但し、ここでは式(40)の変数変換を用いている。式(39)のEをA、r0でそれぞれ微分すると、式(41)、式(42)が得られる。式(41)を0とおいて、これをAについて解くと式(43)となり、これを式(42)に代入すると、Aを消去することができて、式(44)が得られる。また、式(40)と式(43)を用いると、aについての算出式である式(45)が得られる。
【数５５】
【数５６】
【数５７】
【数５８】
【数５９】
【数６０】
【数６１】
【００７１】
以上の式変形を行うことにより、式(44)に含まれる未知パラメータはr0のみとなるので、r0に関しては式(44)を利用して１変数の最小二乗問題として求めることが可能となる。また、aに関しては、r0が求められた後、式(45)にr0を代入して直接求めることができる。
【００７２】
そこで、式(44)を基にしたニュートン法をr0の導出に適用すると、r0の更新式は式(46)となり、これに従って繰返し演算を行うことでr0が算出される。但し、∂2E/∂r02は、式(44)をさらにr0で微分することで求められ、式(47)で表される。また、r0の初期値としては、r0が音響中心の補正量であることを考えれば、r0は0に近い値となることが予想されるので、0を設定するのが妥当である。
【数６２】
【数６３】
【００７３】
ここで、r0を導出するための繰返し演算で用いられる式(46)、式(44)、式(47)、および、aを導出するための式(45)は、逆二乗特性解析装置１１の減衰曲線推定部１２において減衰曲線の推定に用いられる式(4)、式(5)、式(6)、および、式(7)と一致していることがわかる。
【００７４】
つまり、本実施例の逆二乗特性解析装置における減衰曲線の推定方法は、上記の技術的検討を経て導出されたものであり、このような構成にすることにより、減衰曲線を適切に推定することが可能となる。また、繰返し演算にて推定されるパラメータが１つのみの構成となっているため、比較的簡便な演算で処理を行うことが可能であり、適切な収束解が安定的に得られやすいという効果が得られる。
【００７５】
なお、減衰曲線を推定する際に最小化する差分の尺度と、測定値と推定値の乖離量を評価する時の尺度を一致させることの効果を実証するため、以下では、模擬的に作成した入力データを用いて、前提技術の逆二乗特性解析装置と本実施例の逆二乗特性解析装置の解析結果を比較する。
【００７６】
ここで使用する入力データは、測定点数がN=46で、各測定点の音源からの距離riはri=0.5+0.1(i-1) mに従って算出したものであり、音圧値piは式(48)に従って算出される音圧レベルLpiを音圧に変換したものである。なお、式(48)のα(r)は微小な誤差を示しており、|α(r)|の最大値は、r=4.5 mにおいて約1.52 dBである。すなわち、式(48)は、音源から0.5m〜2.5mの範囲では、式(1)においてa=1000、r0=0としたモデル式に対応する理想的な減衰曲線を示しており、2.5m〜5.0mの範囲では、その理想的な減衰曲線に微小な誤差が付加された曲線を示している。
【００７７】
なお、第46番目の測定点は音源から5.0mの位置となるので、検定範囲は音源から直線経路に沿って5.0mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【数６４】
【００７８】
表３は、前提技術の逆二乗特性解析装置の解析結果を示す表である。表３を参照して、第２列と第３列が、入力データとなる各測定点の音源からの距離ri、および、音圧の測定値piである。これらの入力データから前提技術に係る逆二乗特性解析装置において推定される未知パラメータa、r0は、a=888.79、r0=0.208608である。表３の第５列、第６列は、各測定点における音圧レベルの測定値Lpiと推定値Lp(ri)を示しており、第７列はそれらの差分δiを示している。
【表３】
【００７９】
図７は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表３の第５列の値）を丸い点で、推定曲線（式(19)にa=888.79、r0=0.208608を代入した曲線）を実線で示してある。
【数６５】
【００８０】
表３の第７列と図７を参照すると、前提技術の逆二乗特性解析装置の解析結果では、音源に近い測定点において測定値と推定値の乖離が大きくなっており、|δi|の最大値はr=0.5mの位置で約3.67dBとなっている。検定対象としている中心周波数500Hzの1/3オクターブ帯域の許容限界値δlimは2.5dBであるので、検定結果は不適合となる。
【００８１】
しかしながら、解析に用いた入力データは、音源から遠い範囲にのみ誤差が付加されており、音源に近い範囲では理想的な減衰曲線に従うものである。また、理想的な減衰曲線からの乖離量も最大1.52dB程度であり、本来、検定結果は適合となるべきものである。従って、前提技術の逆二乗特性解析装置では、適切な減衰曲線が推定されず、その結果として、誤った検定結果が出力される可能性があることが示された。
【００８２】
これに対し、表４は、本実施例の逆二乗特性解析装置の解析結果を示す表である。表４を参照して、第２列と第３列が、入力データとなる各測定点の音源からの距離ri、および、音圧の測定値piである。これらの入力データから本実施例の逆二乗特性解析装置において推定される未知パラメータa、r0は、a=960.208、r0=0.0389398である。表４の第５列、第６列は、各測定点における音圧レベルの測定値Lpiと推定値Lp(ri)を示しており、第７列はそれらの差分δiを示している。
【表４】
【００８３】
図８は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表４の第５列の値）を丸い点で、推定曲線（式(19)にa=960.208、r0=0.0389398を代入した曲線）を実線で示してある。
【００８４】
表４の第７列と図８を参照すると、本実施例の逆二乗特性解析装置の解析結果では、全ての測定点において測定値と推定値の乖離が一様に小さくなっており、|δi|の最大値はr=4.5mの位置で約1.24dBとなっている。検定対象としている中心周波数500Hzの1/3オクターブ帯域の許容限界値δlimは2.5dBであるので、検定結果は適合となる。
【００８５】
この結果は、解析に用いた入力データが、r=4.5mの位置で最大1.52dBの誤差を付加したものであることとよく一致しており、本実施例の逆二乗特性解析装置では、このような入力データに対しても、適切な減衰曲線が推定されて、正しい検定結果が出力されることが示された。
【００８６】
以上のことから、本実施例によれば、式(1)にてモデル化され、式(3)の評価式を最小にする減衰曲線の推定を、繰返し演算によって推定されるパラメータ数が少なく、かつ、簡素化された推定式にて行うため、少ない演算量にて計算が可能となり、さらに、適切な推定結果を安定的に得ることが可能となる。
【実施例２】
【００８７】
以下では、本発明の逆二乗特性解析装置の別の実施例を説明する。なお、逆二乗則検定を行う手順については、前提技術の逆二乗特性解析装置に代えて、本発明の逆二乗特性解析装置を使用すること以外、前提技術で説明した手順と全て同じであるので説明を省略し、逆二乗特性解析装置の構成と解析動作のみを説明する。
【００８８】
＜本実施例の逆二乗特性解析装置の構成＞
以下では、逆二乗則検定において使用される本実施例の逆二乗特性解析装置の構成を説明する。この逆二乗特性解析装置は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして受け、その距離減衰特性が反比例曲線から乖離する量を算出し、その乖離量に基づいて自由空間としての適合／不適合の解析結果を出力するものである。
【００８９】
図９は、本発明の逆二乗特性解析装置を説明するブロック構成図である。図９を参照して、逆二乗特性解析装置２１は、測定値から反比例曲線を推定する減衰曲線推定部２２と、推定された反比例曲線と測定値との乖離量を算出する差分算出部２３と、算出された乖離量から半無響室の自由空間特性を判定する特性評価部２４を備える。
【００９０】
減衰曲線推定部２２は、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、式(1)でモデル化された反比例曲線の推定を行う。式(1)は、音圧pが距離rに反比例する関係を示した式であるが、未知のパラメータa、r0を含んだ形でモデル化されている。減衰曲線推定部２２では、式(8)、式(9)に従って未知パラメータa、r0を算出し、反比例曲線の推定を完了する。但し、wi (i=1〜N)は式(10)に従ってpiから算出される量である。また、qi(i=1〜N)は第i番目の測定点における逆数音圧であり、式(11)に従ってpiから算出される。なお、式(1)のモデル式はISO 3745に基づくものであるが、式(8)、式(9)、式(10)、式(11)に基づく未知パラメータa、r0の推定方法は、本実施例の逆二乗特性解析装置の特徴となる部分である。
【数６６】
【数６７】
【数６８】
【数６９】
【数７０】
【００９１】
差分算出部２３は、推定されたパラメータa、r0、および、各測定点における音圧値pi (i=1〜N)と音源からの距離ri (i=1〜N)を入力データとして、各測定点における音圧レベルの測定値Lpi (i=1〜N)、音圧レベルの推定値Lp(ri) (i=1〜N)、および、それらの差分δi (i=1〜N)を算出する。ここで、Lpi およびLp(ri)の算出は、それぞれ、式(15)、式(16)に従って行われ、δiの算出は、式(17)に従って行われる。なお、式(15)は、式(2)に基づいて得られるpiとLpiの関であり、式(16)は、式(1)と式(2)に基づいて得られるa、r0、riとp(ri)の関係式である。
【数７１】
【数７２】
【数７３】
【００９２】
特性評価部２４は、差分δi (i=1〜N)を入力データとして、これらのδiに対して式(18)の関係が成立するか否かの判定を行い、全ての測定点、すなわち、i=1〜Nの全てのiにおいて式(18)が成立する場合には適合、そうでない場合には不適合の結果を出力する。なお、δlimは許容限界値であり、式(18)は差分δi の絶対値がこの許容限界値以内に収まることを示す関係式である。δlim の値は、ISO 3745では検定対象とする周波数帯域毎に定められており、ここでもそれ従う。
【数７４】
【００９３】
＜本実施例の逆二乗特性解析装置の動作＞
以下では、前提技術で例示したのと同様に、測定点数が26個（N=26）、各測定点の音圧値pi、距離riが表５の第３列と第２列に示されたデータである場合を例に、本実施例の逆二乗特性解析装置で行われる解析動作を説明する。なお、距離riはri=0.5+0.1(i-1) m の関係があるので、第26番目の測定点は音源から3.0mの位置となり、検定範囲は音源から直線経路に沿って3.0mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【表５】
【００９４】
逆二乗特性解析装置２１において、各測定点における音圧値pi、音源からの距離riが入力されると、減衰曲線推定部２２では、式(8)、式(9)、式(10)、式(11)に従って未知パラメータa、r0の算出が行われる。このとき、表５の第３列と第２列に示されたpi、riの入力データに対しては、a=1039.16、r0=-0.001743が算出される。
【００９５】
未知パラメータa、r0が推定されると、差分算出部２３では、式(15)と式(16)に従い、音圧レベルの測定値Lpiおよび推定値Lp(ri)が算出され、さらに、式(17)に従い、その差分δiが算出される。このときのLpi、Lp(ri)、δiは、表５の第５列、第６列、第７列にそれぞれ記された値となる。なお、図１０は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表５の第５列の値）を丸い点で、推定曲線（式(19)にa=1039.16、r0=-0.001743を代入した曲線）を実線で示してある。
【数７５】
【００９６】
差分δiが算出されると、特性評価部２４では、式(18)に従い、それらの値が許容限界値以内に収まるか否かが判定される。検定対象としている中心周波数500Hzの1/3オクターブ帯域では、許容値δlimは2.5dBと定められており、差分δiの算出結果である表５の第７列の値に対しては全て式(18)が成立するため、適合の解析結果が出力される。この解析結果は、中心周波数が500Hzの1/3オクターブ帯域において、直線経路に沿って音源から3.0mの位置までの範囲が自由音場とみなされることを示している。
【００９７】
以上によって、本実施例の逆二乗特性解析装置の構成、および、各測定点の音圧値pi (i=1〜N)と距離ri (i=1〜N)が入力されてから、自由空間としての適合／不適合を示す結果が出力されるまでの解析動作が説明された。
【００９８】
＜本実施例の逆二乗特性解析装置の特徴と効果＞
以下では、本実施例の逆二乗特性解析装置の特徴と効果について詳しく説明する。
実施例１において明らかにしたように、前提技術の逆二乗特性解析装置において推定される減衰曲線は、式(21)のE、すなわち、逆数音圧の差分の二乗和を最小化するものであり、測定値と推定値の乖離量を評価する尺度である音圧レベルとは一致していないため、適切な検定結果を得ることができない。
【００９９】
この点を考慮すると、減衰曲線推定部において、逆数音圧における差の二乗和ではなく、音圧レベルにおける差の二乗和、すなわち、式(3)のEを最小にするよう減衰曲線の推定を行うことが有効であると考えられるが、簡便な計算によって安定した解が得られるように、そのような推定を行うことは困難である。
【数７６】
【０１００】
本実施例の逆二乗特性解析装置は、上記のような技術的検討を経て、式(3)のEを最小化する減衰曲線の推定に適した別の構成を示したものである。以下では、この点について説明する。
【０１０１】
式(3)のEを最小化する未知パラメータa、r0の算出方法を導出するため、逆数音圧における測定値と推定値の差と、音圧レベルにおける測定値と推定値の差の関係について考えると、測定値と推定値の差が小さい場合には近似的に式(49)の関係が成り立つ。但し、音圧レベルLpの逆数音圧qによる微分dLp/dqには式(34)の関係式を用い、さらにその式の中のpを測定値piで近似している。この式(49)を式(3)に代入すると式(50)が得られる。但し、wiは式(51)に従い、piから算出される量である。
【数７７】
【数７８】
【数７９】
【数８０】
【０１０２】
推定値q(ri)は、式(12)、式(1)に基づいて式(22)のように表されるので、これを式(50)に代入すると、Eは式(52)のように式変形される。
【数８１】
【数８２】
【０１０３】
ここで、式(52)のEを最小にする未知パラメータa、r0は、Eをa、r0でそれぞれ微分して0とおいた連立方程式を解いて求められるが、この連立方程式は複雑で容易に解くことができない。そこで、式(25)の変数変換を用いて式(52)を式(53)のように式変形し、式(53)のEを最小化するA、Bを求めた後、式(25)の関係を用いてA、Bからa、r0を求める。
【数８３】
【数８４】
【０１０４】
式(53)のEを最小化するA、Bを求めるため、式(53)をA、Bでそれぞれ微分すると式(54)が得られ、これらを0とおいてA、Bについて整理すると式(55)の連立方程式が得られる。式(55)はA、Bの線形方程式として容易に式(56)のように式変形することができ、さらにこれを式(25)の関係式を用いてa、r0について解くと、式(57)、式(58)が得られる。これらはa、r0を直接算出できる式となっているので、以上の式変形を行うことによって、式(3)のEを最小化するa、r0の算出が可能となる。
【数８５】
【数８６】
【数８７】
【数８８】
【数８９】
【０１０５】
ここで、a、r0を導出するための式(57)、式(58)、式(51)は、逆二乗特性解析装置２１の減衰曲線推定部２２において減衰曲線の推定に用いられる式(8)、式(9)、式(10)と一致していることがわかる。
【０１０６】
つまり、本実施例の逆二乗特性解析装置における減衰曲線の推定方法は、上記の技術的検討を経て導出されたものであり、このような構成にすることにより、減衰曲線を適切に推定することが可能になる。また、繰返し演算を行う必要がない構成となっているため、簡便な演算で処理を行うことが可能であり、解が確実に得られるという効果がある。
【０１０７】
なお、減衰曲線を推定する際に最小化する差分の尺度と、測定値と推定値の乖離量を評価する時の尺度を一致させることの効果を実証するため、以下では、模擬的に作成した入力データを用いて、前提技術の逆二乗特性解析装置と本実施例の逆二乗特性解析装置の解析結果を比較する。
【０１０８】
ここで使用する入力データは、実施例１と同様に、測定点数がN=46で、各測定点の音源からの距離riはri=0.5+0.1(i-1) mに従って算出したものであり、音圧値piは式(48)に従って算出される音圧レベルLpiを音圧に変換したものである。なお、式(48)のα(r)は微小な誤差を示しており、|α(r)|の最大値は、r=4.5 mにおいて約1.52 dBである。すなわち、式(48)は、音源から0.5m〜2.5mの範囲では、式(1)においてa=1000、r0=0としたモデル式に対応する理想的な減衰曲線を示しており、2.5m〜5.0mの範囲では、その減衰曲線に微小な誤差が付加された曲線を示している。
【０１０９】
なお、第46番目の測定点は音源から5.0mの位置となるので、検定範囲は音源から直線経路に沿って5.0mの位置までとなる。また、piは中心周波数が500Hzの1/3オクターブ帯域の音圧値であるものとし、この周波数帯域を検定対象とする。この帯域の許容限界値δlimは、ISO 3745では2.5dBと定められている。
【数９０】
【０１１０】
実施例１において説明したように、この入力データに対する前提技術の逆二乗特性解析装置の解析結果は、音源に近い測定点において、測定値と推定値の乖離が大きくなり、|δi|の最大値はr=0.5mの位置で約3.67dBとなって、検定結果は不適合となった。つまり、前提技術に係る逆二乗特性解析装置では、適切な減衰曲線が推定されず、その結果として、誤った検定結果が出力される可能性があることが示された。
【０１１１】
これに対し、表６は、本実施例の逆二乗特性解析装置の解析結果を示す表である。表６を参照して、第２列と第３列が、入力データとなる各測定点の音源からの距離ri、および、音圧の測定値piである。これらの入力データから本実施例の逆二乗特性解析装置において推定される未知パラメータa、r0は、a=969.883、r0=0.0309326である。表６の第５列、第６列は、各測定点における音圧レベルの測定値Lpiと推定値Lp(ri)を示しており、第７列はそれらの差分δiを示している。
【表６】
【０１１２】
図１１は、この例における音圧レベルの測定値と推定曲線をグラフ化した図である。横軸が距離r、縦軸が音圧レベルLpで、各測定点における測定値（表６の第５列の値）を丸い点で、推定曲線（式(19)にa=969.883、r0=0.0309326を代入した曲線）を実線で示してある。
【数９１】
【０１１３】
表６の第７列と図１１を参照すると、本実施例の逆二乗特性解析装置の解析結果では、全ての測定点において測定値と推定値の乖離が一様に小さくなっており、|δi|の最大値はr=4.5mの位置で約1.32dBとなっている。検定対象としている中心周波数500Hzの1/3オクターブ帯域の許容限界値δlimは2.5dBであるので、検定結果は適合となる。
【０１１４】
この結果は、解析に用いた入力データが、r=4.5mの位置で最大1.52dBの誤差を付加したものであることとよく一致しており、本実施例の逆二乗特性解析装置では、このような入力データに対しても、適切な減衰曲線が推定されて、正しい検定結果が出力されることが示された。
【０１１５】
以上のことから、本実施例によっても、式(1)にてモデル化され、式(3)の評価式を最小にする減衰曲線の推定を、簡素化された推定式にて行うため、少ない演算量にて計算が可能となり、さらに、適切な推定結果を安定的に得ることが可能となる。
【符号の説明】
【０１１６】
１、１１、２１…逆二乗特性解析装置
２、１２、２２…減衰曲線推定部
３、１３、２３…差分算出部
４、１４、２４…特性評価部
Ｓ１…処理ステップ（パラメータr0への初期値設定）
Ｓ２…処理ステップ（パラメータr0の更新）
Ｓ３…処理ステップ（収束判定）
【特許請求の範囲】
【請求項１】
試験音場内の複数の測定点における音圧の測定値と
それぞれの測定点の位置とを入力データとして、
音圧と位置の関係を反比例で近似した減衰曲線を推定する減衰曲線推定手段と、
前記推定された減衰曲線と前記音圧の測定値との差を
音圧レベルの尺度において算出する差分算出手段と、
前記算出した差に基づいて前記試験音場の自由空間特性を評価する特性評価手段を
備える逆二乗特性解析装置であって、
前記減衰曲線推定手段は、
前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように
前記減衰曲線の推定を行うことを特徴とする
逆二乗特性解析装置。
【請求項２】
試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、
第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、
発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、
前記pi、ri (i=1〜N)を入力データとして、
式(1)で表される減衰曲線の未知パラメータa、r0を推定する減衰曲線推定手段と、
前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、
式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、
前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出する差分算出手段と、
前記算出された差が所定の許容値に収まるか否かによって
前記試験音場の自由空間特性を評価する特性評価手段を
備える逆二乗特性解析装置であって、
前記減衰曲線推定手段は、式(3)を最小にするように
前記未知パラメータa、r0の推定を行うことを特徴とする
逆二乗特性解析装置。
【数１】
【数２】
【数３】
【請求項３】
請求項２に記載の逆二乗特性解析装置であって、
前記減衰曲線推定手段は、
前記未知パラメータr0に初期値を設定し、
式(4)、式(5)及び式(6)に従って前記未知パラメータr0を更新する処理を
前記未知パラメータr0が収束するまで繰返し行って前記未知パラメータr0の収束解を算出し、
前記算出した未知パラメータr0の収束解を用いて、
式(7)に従って前記未知パラメータaを算出することを特徴とする
逆二乗特性解析装置。
【数４】
【数５】
【数６】
【数７】
【請求項４】
請求項２に記載の逆二乗特性解析装置であって、
前記減衰曲線推定手段は、
式(8)、式(9)、式(10)及び式(11)に従って、前記未知パラメータa、r0の推定を行うことを特徴とする
逆二乗特性解析装置。
【数８】
【数９】
【数１０】
【数１１】
【請求項５】
試験音場内の複数の測定点における音圧の測定値と
それぞれの測定点の位置とを入力データとして、
音圧と位置の関係を反比例で近似した減衰曲線を推定するステップと、
前記推定された減衰曲線と前記音圧の測定値との差を
音圧レベルの尺度において算出するステップと、
前記算出した差に基づいて前記試験音場の自由空間特性を評価するステップからなる
試験音場内の逆二乗特性の評価方法であって、
前記減衰曲線を推定するステップでは、
前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように
前記減衰曲線の推定を行うことを特徴とする方法。
【請求項６】
試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、
第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、
発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、
前記pi、ri (i=1〜N)を入力データとして、
式(1)で表される減衰曲線の未知パラメータa、r0を推定するステップと、
前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、
式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、
前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出するステップと、
前記算出された差が所定の許容値に収まるか否かによって
前記試験音場の自由空間特性を評価するステップからなる
試験音場内の逆二乗特性の評価方法であって、
前記減衰曲線を推定するステップは、式(3)を最小にするように
前記未知パラメータa、r0の推定を行うことを特徴とする方法。
【数１２】
【数１３】
【数１４】
【請求項１】
試験音場内の複数の測定点における音圧の測定値と
それぞれの測定点の位置とを入力データとして、
音圧と位置の関係を反比例で近似した減衰曲線を推定する減衰曲線推定手段と、
前記推定された減衰曲線と前記音圧の測定値との差を
音圧レベルの尺度において算出する差分算出手段と、
前記算出した差に基づいて前記試験音場の自由空間特性を評価する特性評価手段を
備える逆二乗特性解析装置であって、
前記減衰曲線推定手段は、
前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように
前記減衰曲線の推定を行うことを特徴とする
逆二乗特性解析装置。
【請求項２】
試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、
第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、
発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、
前記pi、ri (i=1〜N)を入力データとして、
式(1)で表される減衰曲線の未知パラメータa、r0を推定する減衰曲線推定手段と、
前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、
式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、
前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出する差分算出手段と、
前記算出された差が所定の許容値に収まるか否かによって
前記試験音場の自由空間特性を評価する特性評価手段を
備える逆二乗特性解析装置であって、
前記減衰曲線推定手段は、式(3)を最小にするように
前記未知パラメータa、r0の推定を行うことを特徴とする
逆二乗特性解析装置。
【数１】
【数２】
【数３】
【請求項３】
請求項２に記載の逆二乗特性解析装置であって、
前記減衰曲線推定手段は、
前記未知パラメータr0に初期値を設定し、
式(4)、式(5)及び式(6)に従って前記未知パラメータr0を更新する処理を
前記未知パラメータr0が収束するまで繰返し行って前記未知パラメータr0の収束解を算出し、
前記算出した未知パラメータr0の収束解を用いて、
式(7)に従って前記未知パラメータaを算出することを特徴とする
逆二乗特性解析装置。
【数４】
【数５】
【数６】
【数７】
【請求項４】
請求項２に記載の逆二乗特性解析装置であって、
前記減衰曲線推定手段は、
式(8)、式(9)、式(10)及び式(11)に従って、前記未知パラメータa、r0の推定を行うことを特徴とする
逆二乗特性解析装置。
【数８】
【数９】
【数１０】
【数１１】
【請求項５】
試験音場内の複数の測定点における音圧の測定値と
それぞれの測定点の位置とを入力データとして、
音圧と位置の関係を反比例で近似した減衰曲線を推定するステップと、
前記推定された減衰曲線と前記音圧の測定値との差を
音圧レベルの尺度において算出するステップと、
前記算出した差に基づいて前記試験音場の自由空間特性を評価するステップからなる
試験音場内の逆二乗特性の評価方法であって、
前記減衰曲線を推定するステップでは、
前記音圧の測定値と前記減衰曲線との差を音圧レベルの尺度において最小化するように
前記減衰曲線の推定を行うことを特徴とする方法。
【請求項６】
試験音場内に設置された発音体を端点とする直線経路上のN個の測定点において、
第i番目 (i=1〜N)の測定点におけるパスカルで表される音圧の測定値をpi、
発音体からのメートルで表される距離をriとし、さらに、基準音圧を20μPaとして、これをp0としたとき、
前記pi、ri (i=1〜N)を入力データとして、
式(1)で表される減衰曲線の未知パラメータa、r0を推定するステップと、
前記推定された未知パラメータa、r0と前記pi、ri (いずれもi=1〜N)を入力データとして、
式(2)に従う音圧レベルの測定値Lpi (i=1〜N)と式(1)及び式(2)から求められる音圧レベルの推定値Lp(ri) (i=1〜N)を算出し、
前記Lpi及びLp(ri) (いずれもi=1〜N)を入力データとして、両者の差を算出するステップと、
前記算出された差が所定の許容値に収まるか否かによって
前記試験音場の自由空間特性を評価するステップからなる
試験音場内の逆二乗特性の評価方法であって、
前記減衰曲線を推定するステップは、式(3)を最小にするように
前記未知パラメータa、r0の推定を行うことを特徴とする方法。
【数１２】
【数１３】
【数１４】
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２０１３−８３５１５（Ｐ２０１３−８３５１５Ａ）
【公開日】平成２５年５月９日（２０１３．５．９）
【国際特許分類】
【出願番号】特願２０１１−２２２８６６（Ｐ２０１１−２２２８６６）
【出願日】平成２３年１０月７日（２０１１．１０．７）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【Ｆターム（参考）】
【公開日】平成２５年５月９日（２０１３．５．９）
【国際特許分類】
【出願日】平成２３年１０月７日（２０１１．１０．７）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【Ｆターム（参考）】
[ Back to top ]