知財求人 - 知財ポータルサイト「IP Force」
特開2023-62922モーダルパラメータの同定方法、モーダルパラメータの同定装置およびプログラム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023062922
(43)【公開日】2023-05-09
(54)【発明の名称】モーダルパラメータの同定方法、モーダルパラメータの同定装置およびプログラム
(51)【国際特許分類】
G10D 3/02 20060101AFI20230427BHJP
G10D 13/22 20200101ALI20230427BHJP
G10C 3/06 20060101ALI20230427BHJP
【ＦＩ】
G10D3/02
G10D13/22
G10C3/06 100
【審査請求】未請求
【請求項の数】4
【出願形態】ＯＬ
(21)【出願番号】P 2021173110
(22)【出願日】2021-10-22
(71)【出願人】
【識別番号】000004075
【氏名又は名称】ヤマハ株式会社
(74)【代理人】
【識別番号】110003177
【氏名又は名称】弁理士法人旺知国際特許事務所
(72)【発明者】
【氏名】富永 英嗣
【テーマコード（参考）】
5D002
【Ｆターム（参考）】
5D002CC01
(57)【要約】
【課題】「空気中に存在する場合の構造物」または「構造物が付随する空気」の時間領域シミュレーションを高速に行う。
【解決手段】コンピューター（１０）が、空気中に存在する場合の構造物の周波数応答、または、構造物が付随する場合の前記空気の周波数応答を算出する第１過程と、コンピューター（１０）が、周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第２過程と、コンピューター（１０）が、周波数応答をリファレンスとして、構造物における加振点の複素固有モード、または、空気における加振点の複素固有モードを同定する第３過程と、を含む。
【選択図】図３
【解決手段】コンピューター（１０）が、空気中に存在する場合の構造物の周波数応答、または、構造物が付随する場合の前記空気の周波数応答を算出する第１過程と、コンピューター（１０）が、周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第２過程と、コンピューター（１０）が、周波数応答をリファレンスとして、構造物における加振点の複素固有モード、または、空気における加振点の複素固有モードを同定する第３過程と、を含む。
【選択図】図３
【特許請求の範囲】
【請求項1】
コンピューターが、空気中に存在する場合の構造物の物理モデル、または、構造物が付随する場合の空気の物理モデルにおけるモーダルパラメータを同定するモーダルパラメータ同定方法であって、
前記コンピューターが、前記空気中に存在する場合の前記構造物の周波数応答、または、前記構造物が付随する場合の前記空気の周波数応答を算出する第１過程と、
前記コンピューターが、前記周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第２過程と、
前記コンピューターが、前記周波数応答をリファレンスとして、前記構造物における加振点の複素固有モード、または、前記空気における加振点の複素固有モードを同定する第３過程と、
を含む物理モデルにおけるモーダルパラメータの同定方法。
前記コンピューターが、前記空気中に存在する場合の前記構造物の周波数応答、または、前記構造物が付随する場合の前記空気の周波数応答を算出する第１過程と、
前記コンピューターが、前記周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第２過程と、
前記コンピューターが、前記周波数応答をリファレンスとして、前記構造物における加振点の複素固有モード、または、前記空気における加振点の複素固有モードを同定する第３過程と、
を含む物理モデルにおけるモーダルパラメータの同定方法。
【請求項2】
前記第３過程の後に、
前記コンピューターが、前記周波数応答をリファレンスとして、前記構造物における前記加振点とは異なる応答点の複素固有モード、または、前記空気における前記加振点とは異なる応答点の複素固有モードを同定する第４過程を、
を含む請求項１に記載の物理モデルにおけるモーダルパラメータの同定方法。
前記コンピューターが、前記周波数応答をリファレンスとして、前記構造物における前記加振点とは異なる応答点の複素固有モード、または、前記空気における前記加振点とは異なる応答点の複素固有モードを同定する第４過程を、
を含む請求項１に記載の物理モデルにおけるモーダルパラメータの同定方法。
【請求項3】
コンピューターに、空気中に存在する場合の構造物の物理モデル、または、構造物が付随する場合の空気の物理モデルにおけるモーダルパラメータを同定させるプログラムであって、
前記コンピューターに、
前記空気中に存在する場合の前記構造物の周波数応答、または、前記構造物が付随する場合の前記空気の周波数応答を算出する算出部と、
前記周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第１同定部、および、
前記周波数応答をリファレンスとして、前記構造物における加振点の複素固有モード、または、前記空気における加振点の複素固有モードを同定する第２同定部、
として機能させるプログラム。
前記コンピューターに、
前記空気中に存在する場合の前記構造物の周波数応答、または、前記構造物が付随する場合の前記空気の周波数応答を算出する算出部と、
前記周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第１同定部、および、
前記周波数応答をリファレンスとして、前記構造物における加振点の複素固有モード、または、前記空気における加振点の複素固有モードを同定する第２同定部、
として機能させるプログラム。
【請求項4】
空気中に存在する場合の構造物の物理モデル、または、構造物が付随する場合の空気の物理モデルにおけるモーダルパラメータを同定するモーダルパラメータ同定装置であって、
前記空気中に存在する場合の前記構造物の周波数応答、または、前記構造物が付随する場合の前記空気の周波数応答を算出する算出部と、
前記周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第１同定部と、
前記周波数応答をリファレンスとして、前記構造物における加振点の複素固有モード、または、前記空気における加振点の複素固有モードを同定する第２同定部と、
を含むモーダルパラメータの同定装置。
前記空気中に存在する場合の前記構造物の周波数応答、または、前記構造物が付随する場合の前記空気の周波数応答を算出する算出部と、
前記周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第１同定部と、
前記周波数応答をリファレンスとして、前記構造物における加振点の複素固有モード、または、前記空気における加振点の複素固有モードを同定する第２同定部と、
を含むモーダルパラメータの同定装置。
【発明の詳細な説明】
【技術分野】
【0001】
本開示は、モーダルパラメータの同定方法、モーダルパラメータの同定装置およびプログラムに関する。
【背景技術】
【0002】
構造物の物理モデルにおける動特性を求める手法としては、次のような技術が挙げられる。例えば、非特許文献１には、１点加振多点参照（SIMO: Single-Input Multi-Output）方式で求めた構造物のインパルス応答をリファレンスとして用いることで、一般粘性減衰系モデルにおける複素固有値を同定する技術が記載されている。また、非特許文献２には、多点加振多点参照（MIMO: Multi-Input Multi-Output）方式で求めた構造物の周波数応答をリファレンスとして用いることで、一般粘性減衰系モデルにおける複素固有値と複素固有モードとを同定する技術が記載されている。
【先行技術文献】
【非特許文献】
【0003】
【非特許文献1】Kerem Ege, Xavier Boutillon, Bertrand David. High-resolution modal analysis. Journal of Sound and Vibration, Elsevier, 2009, 325 (4-5), p.852-869.
【非特許文献2】日本音響学会 編、「音・振動のモード解析と制御」コロナ社、1996年、p.110-114.
【発明の概要】
【発明が解決しようとする課題】
【0004】
しかしながら、非特許文献１に記載された技術では、１点加振多点参照方式であるために、複素固有モードの高精度な同定が困難である。非特許文献２に記載された技術では、多点加振多点参照方式であるために、非特許文献１に記載された技術を用いた場合に比べ複素固有モードの同定精度は高いものの、複素固有値の高精度な同定が困難であり、また、複素固有値と複素固有モードとを同時に同定するので、大規模な構造物への適用が困難である。
【課題を解決するための手段】
【0005】
本開示の一態様に係るモーダルパラメータの同定方法は、コンピューターが、空気中に存在する場合の構造物の物理モデル、または、構造物が付随する場合の空気の物理モデルにおけるモーダルパラメータを同定するモーダルパラメータ同定方法であって、前記コンピューターが、前記空気中に存在する場合の前記構造物の周波数応答、または、前記構造物が付随する場合の前記空気の周波数応答を算出する第１過程と、前記コンピューターが、前記周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第２過程と、前記コンピューターが、前記周波数応答をリファレンスとして、前記構造物における加振点の複素固有モード、または、前記空気における加振点の複素固有モードを同定する第３過程と、を含む。
【図面の簡単な説明】
【0006】
【図1】物理モデルを説明するための図である。
【図2】モーダルパラメータの同定処理を実行するコンピューターを示す図である。
【図3】コンピューターにおける機能的な構成を示すブロック図である。
【図4】モーダルパラメータの同定処理を実行するコンピューターを示す図である。
【図5】モーダルパラメータの同定処理を実行するコンピューターを示す図である。
【図6】参照ＦＲとモード合成によるＦＲとがほぼ一致していることを示す図である。
【図7】参照ＦＲとモード合成によるＦＲとがほぼ一致していることを示す図である。
【図8】参照ＦＲとモード合成によるＦＲとがほぼ一致していることを示す図である。
【図9】参照ＦＲとモード合成によるＩＲとがほぼ一致していることを示す図である。
【図10】参照ＦＲとモード合成によるＩＲとがほぼ一致していることを示す図である。
【図11】参照ＦＲとモード合成によるＩＲとがほぼ一致していることを示す図である。
【図12】真空中の周波数応答と空気中の周波数応答とを示す図である。
【発明を実施するための形態】
【0007】
以下、本開示の実施形態に係るモーダルパラメータの同定方法について図面を参照して説明する。
なお、各図において、各部の寸法および縮尺は、実際のものと適宜に異ならせてある。また、以下に述べる実施の形態は、好適な具体例であるから、技術的に好ましい種々の限定が付されているが、本発明の範囲は、以下の説明において特に本発明を限定する旨の記載がない限り、これらの形態に限られるものではない。
なお、各図において、各部の寸法および縮尺は、実際のものと適宜に異ならせてある。また、以下に述べる実施の形態は、好適な具体例であるから、技術的に好ましい種々の限定が付されているが、本発明の範囲は、以下の説明において特に本発明を限定する旨の記載がない限り、これらの形態に限られるものではない。
【0008】
まず、線形を仮定することが可能な構造物の時間領域シミュレーションを実行する場合、構造物モデルの変位を固有モードの重ね合わせとして表現する「モード合成法」は、コストパフォーマンスの最も高い解析手法として知られている。真空中に置かれた構造物モデルであれば、当該構造物のモーダルパラメータ、即ち、実固有値および実固有モードを高い精度で算出する手法が既に確立している。有限要素法（Finite Element Method : FEM）を用いた固有値解析は、このための最も強力な手法として知られている。
【0009】
一方、空気中、詳細には、ある程度以上の大きさの室内、無音響室または屋外に置かれた大規模な構造物モデルの高精度なモーダルパラメータの算出は、極めて困難であると言わざるを得ない。この理由は、空気中に置かれた構造物モデルに対して、空気と構造との連成を考慮した上で有限要素法を用いて固有モードを同定することは、モデルが大きくなるにつれ、あるいは、解析周波数が高くなるにつれ、計算の規模が著しく大きくなり、その結果、計算コストが非現実的に高くなってしまうためである。音響構造連成問題に対する有限要素法にかわるコストパフォーマンスの高い数値解析手法として高速多重極境界要素法 (Fast Multipole Boundary Element Method : FMBEM）が知られているが、この手法で求まるものは周波数応答であって、固有モードではない。要するに、大規模な「空気中に置かれた構造物」の固有モードを現実的な計算コストで高精度に算出する方法はこれまで提案されてこなかった。
【0010】
図１２は、構造物の一例として、ある楽器の響板について、真空中の周波数応答と、空気中の周波数応答とを計算した結果を示す図である。この例が示す通り、構造物によっては、真空中の周波数応答と、空気中の周波数応答の間に無視できない差がある場合がある。このような差が生じる主な理由は、空気中では、真空中と比較して、空気の質量効果によって固有振動数が低くなることに加え、振動のエネルギーが音響放射のために消費されて減衰が大きくなるためである。
なお、図１２において、当該構造物の真空中における周波数応答は、当該構造物の有限要素固有値解析の結果として得られる固有モードの重ね合わせとして算出した応答である。空気中における周波数応答は、上記の真空中構造物の有限要素固有値解析の結果を用いながら、構造物が無響室に置かれている場合を仮定して、高速多重極境界要素法による音響構造連成解析を実施することで算出した応答である。
なお、図１２において、当該構造物の真空中における周波数応答は、当該構造物の有限要素固有値解析の結果として得られる固有モードの重ね合わせとして算出した応答である。空気中における周波数応答は、上記の真空中構造物の有限要素固有値解析の結果を用いながら、構造物が無響室に置かれている場合を仮定して、高速多重極境界要素法による音響構造連成解析を実施することで算出した応答である。
【0011】
さて、構造物の固有モードを実験的に同定する手法は、非特許文献１あるいは非特許文献２に記載されている方法を含め、これまでに多くの手法が提案されており、実験モード同定ソフトとして製品化もされている。
しかしながら、実験では高品質の周波数応答を大量に取得することがそもそも困難なために、実験的モード同定を目的として開発されてきた従来の手法では、せいぜい数個から数十個程度の固有モードを同定することしか想定されておらず、必然的に大規模モデルには対応できない。
しかしながら、実験では高品質の周波数応答を大量に取得することがそもそも困難なために、実験的モード同定を目的として開発されてきた従来の手法では、せいぜい数個から数十個程度の固有モードを同定することしか想定されておらず、必然的に大規模モデルには対応できない。
【0012】
図１は、空気中に置かれた構造物、あるいは、構造物が付随する空気に対する入力と応答を示す概念図である。
【0013】
本実施形態では、第１に、空気と構造物との連成系を特に線形連成系Ｌｃと表記し、当該線形連成系Ｌｃを一般粘性減衰系モデルで表現することにした。
本体と空気とから成る連成系、即ち、音響構造連成系を一般粘性減衰系モデルで表現できることの理論的根拠は、構造変位と空気速度ポテンシャルを「時間および空間に依存する変数」とした場合の有限要素法による音響構造連成運動方程式が、形式的に実対称の質量行列、減衰行列、剛性行列によって書き表せることにある。この３つの実対称行列から構成される運動方程式の形は、（比例粘性減衰ではない）一般的な粘性減衰を有する構造物のモデル、即ち一般粘性減衰系モデルの有限要素法による運動方程式の形と数学的に同じである。
なお、音響構造連成系を一般粘性減衰系モデルで表現することの利点は、このモデルの固有モードが広義の直交性を有することにある。このとき、モード合成法を用いた高速な時間領域シミュレーションが可能になるのである。
本体と空気とから成る連成系、即ち、音響構造連成系を一般粘性減衰系モデルで表現できることの理論的根拠は、構造変位と空気速度ポテンシャルを「時間および空間に依存する変数」とした場合の有限要素法による音響構造連成運動方程式が、形式的に実対称の質量行列、減衰行列、剛性行列によって書き表せることにある。この３つの実対称行列から構成される運動方程式の形は、（比例粘性減衰ではない）一般的な粘性減衰を有する構造物のモデル、即ち一般粘性減衰系モデルの有限要素法による運動方程式の形と数学的に同じである。
なお、音響構造連成系を一般粘性減衰系モデルで表現することの利点は、このモデルの固有モードが広義の直交性を有することにある。このとき、モード合成法を用いた高速な時間領域シミュレーションが可能になるのである。
【0014】
一般粘性減衰系モデルにおけるインパルス応答（Impulse Response）は、次式（１）のように示される。
また、一般粘性減衰系モデルにおける周波数応答（Frequency Response）は、次式（２）のように示される。
【数1】
【数2】
【0015】
式（１）または式（２）において、ｊは虚数単位であり、＊は複素共役であり、ｒはモード次数である。ｉは加振点のインデクッスであり、ｌは応答点のインデックスである。
λｒは複素固有値であり、λｒ＝－σｒ＋ｊωｄｒで示される。
σｒはモード減衰率であり、ωｄｒは減衰固有角振動数である。
λｒは複素固有値であり、λｒ＝－σｒ＋ｊωｄｒで示される。
σｒはモード減衰率であり、ωｄｒは減衰固有角振動数である。
【0016】
Ｒｒｉｌはモード留数であり、次式（３）で示される。
なお、式（３）において、
ｘｒｉは、加振点ｉにおける複素固有モードの実部であり、
ｙｒｉは、加振点ｉにおける複素固有モードの虚部である。
【数3】
ｘｒｉは、加振点ｉにおける複素固有モードの実部であり、
ｙｒｉは、加振点ｉにおける複素固有モードの虚部である。
【0017】
線形連成系Ｌｃを一般粘性減衰系モデルで表現する場合に重要となる点は、当該一般粘性減衰系モデルのモーダルパラメータをいかに精度よく決定することができるか、という点にある。
このため、本実施形態では、第２に「高速多重極境界要素法を用いて多点加振多点参照（MIMO）方式で算出した線形連成系の周波数応答」をリファレンスとして、一般粘性減衰系モデルにフィットするようにモーダルパラメータである複素固有値と複素固有モードを同定することにした。
このため、本実施形態では、第２に「高速多重極境界要素法を用いて多点加振多点参照（MIMO）方式で算出した線形連成系の周波数応答」をリファレンスとして、一般粘性減衰系モデルにフィットするようにモーダルパラメータである複素固有値と複素固有モードを同定することにした。
【0018】
図２は、本実施形態に係るモーダルパラメータの同定方法を実行するコンピューター１０のブロック図である。本実施形態に係るモーダルパラメータの同定方法は、コンピューターとソフトウェアとの協働により実現される。
コンピューター１０は、制御装置１１と記憶装置１２とを含む。
コンピューター１０は、制御装置１１と記憶装置１２とを含む。
【0019】
制御装置１１は、例えばＣＰＵ（Central Processing Unit）等の単数または複数の処理回路で構成され、コンピューター１０の各要素を統括的に制御する。なお、制御装置１１は、ＣＰＵのほか、ＤＳＰ（Digital Signal Processor）やＡＳＩＣ（Application Specific Integrated Circuit）等の回路によって構成されてもよい。
【0020】
記憶装置１２は、例えば磁気記録媒体または半導体記録媒体等の公知の記録媒体で構成された単数または複数のメモリーであり、制御装置１１が実行するプログラムと制御装置１１が使用する各種のデータとを記憶する。なお、複数種の記録媒体の組合せにより記憶装置１２を構成してもよい。また、コンピューター１０に対して着脱可能な可搬型の記録媒体、または、コンピューター１０が通信網を介して通信可能な外部記録媒体（例えばオンラインストレージ）を、記憶装置１２として利用してもよい。
【0021】
図３は、制御装置１１の機能的な構成を例示するブロック図である。制御装置１１は、記憶装置１２に記憶されたプログラムを実行することで、モーダルパラメータを同定するための複数の要素（処理制御部１１０、前処理部１１１、算出部１１２、第１同定部１１５、第２同定部１１６および第３同定部１１７）として機能する。
なお、制御装置１１における機能の一部を他の装置、例えばネットワークを介して接続されたサーバ装置に負担させる構成としてもよい。
なお、制御装置１１における機能の一部を他の装置、例えばネットワークを介して接続されたサーバ装置に負担させる構成としてもよい。
【0022】
次に、コンピューター１０においてモーダルパラメータを同定する手順について説明する。
図４は、モーダルパラメータを同定する手順を例示するフローチャートである。
図４は、モーダルパラメータを同定する手順を例示するフローチャートである。
【0023】
まず、コンピューター１０における処理制御部１１０は、前処理部１１１に対して次のステップＳ１１からステップＳ１４までの処理を実行させて、線形連成系Ｌｃを表す一般粘性減衰系モデルのモーダルパラメータの同定に必要な準備をする。
具体的には、前処理部１１１は、構造物を構成する材料についての弾性係数およびヒステリシス減衰係数の三次元直交異方性を考慮した上で、当該構造物の三次元有限要素モデルを作成し、当該三次元有限要素モデルの質量行列、減衰行列および剛性行列を算出する（ステップＳ１１）。
前処理部１１１は、当該三次元有限要素モデルに対し、接着、ボルト締めなどの製造工程に伴う接触非線形と幾何学的非線形とを考慮して、非線形静解析を実施し、ステップＳ１１で算出した剛性行列を修正する（ステップＳ１２）。なお、ステップＳ１２は、処理の簡略化のためにスキップしてもよい。
具体的には、前処理部１１１は、構造物を構成する材料についての弾性係数およびヒステリシス減衰係数の三次元直交異方性を考慮した上で、当該構造物の三次元有限要素モデルを作成し、当該三次元有限要素モデルの質量行列、減衰行列および剛性行列を算出する（ステップＳ１１）。
前処理部１１１は、当該三次元有限要素モデルに対し、接着、ボルト締めなどの製造工程に伴う接触非線形と幾何学的非線形とを考慮して、非線形静解析を実施し、ステップＳ１１で算出した剛性行列を修正する（ステップＳ１２）。なお、ステップＳ１２は、処理の簡略化のためにスキップしてもよい。
【0024】
前処理部１１１は、ステップＳ１１で算出した、または、ステップＳ１２で修正した、質量行列および剛性行列をインプットとして、当該三次元有限要素モデルが真空中に存在する場合の、即ち空気が存在しない場合の、当該三次元有限要素モデルの実固有値および実固有モードを固有値解析により算出する（ステップＳ１３）。
前処理部１１１は、ステップＳ１1で算出した減衰行列とステップＳ１３で算出した実固有モードとをインプットとして、当該構造物の三次元有限要素モデルが真空中に存在する場合のモード減衰行列を算出する（ステップＳ１４）。
前処理部１１１は、ステップＳ１1で算出した減衰行列とステップＳ１３で算出した実固有モードとをインプットとして、当該構造物の三次元有限要素モデルが真空中に存在する場合のモード減衰行列を算出する（ステップＳ１４）。
【0025】
次に、処理制御部１１０は、算出部１１２、第１同定部１１５、第２同定部１１６および第３同定部１１７に、次の処理を実行させて、一般粘性減衰系モデルのモーダルパラメータ（複素固有値および複素固有モード）を同定する（ステップＳ１５）。
【0026】
図５は、図４におけるステップＳ１５の手順を詳細に例示するフローチャートである。
算出部１１２は、線形連成系Ｌｃの三次元モデルを解析して、当該線形連成系Ｌｃの周波数応答を、例えば高速多重境界法（Fast Multipole Boundary Element Method）または有限要素法を用いて算出する（ステップＳ１５１）。この解析で算出される周波数応答は、線形連成系Ｌｃにおける加振点ｉでの力に対して、応答点ｌの変位の応答である。
算出部１１２は、線形連成系Ｌｃの三次元モデルを解析して、当該線形連成系Ｌｃの周波数応答を、例えば高速多重境界法（Fast Multipole Boundary Element Method）または有限要素法を用いて算出する（ステップＳ１５１）。この解析で算出される周波数応答は、線形連成系Ｌｃにおける加振点ｉでの力に対して、応答点ｌの変位の応答である。
【0027】
ここでは、加振点ｉの総数をｍ個、応答点ｌの総数をｎ個とする。ある一つの加振位置において加振方向が例えばｘ方向、ｙ方向およびｚ方向のように複数ある場合には、それぞれを１個の加振点として区別することにする。応答点ｌについても同様である。応答点ｌの集合は、加振点ｉの集合を含んでいるものとする。ｍとｎとは、ｎ≧ｍ≧２を満たす整数である。
【0028】
第１同定部１１５は、算出した周波数応答の逆高速フーリエ変換（Inverse Fast Fourier Transform）をリファレンスとして、線形連成系における複素固有値λｒ（＝－σｒ＋ｊωｄｒ）を例えばESPRIT（Estimation of Signal Parameters via Rotational Invariance Techniques）により同定する（ステップＳ１５２）。
なお、複素固有値の同定に際し、例えば上記ｍ点における自己インパルス応答が参照される。また、複素固有値の同定についてはProny法を用いてもよいが、ESPRITに比べ同定精度が劣る傾向がある。
なお、複素固有値の同定に際し、例えば上記ｍ点における自己インパルス応答が参照される。また、複素固有値の同定についてはProny法を用いてもよいが、ESPRITに比べ同定精度が劣る傾向がある。
【0029】
第２同定部１１６は、算出した周波数応答をリファレンスとして、線形連成系Ｌｃの加振点ｉにおける複素固有モードを、非線形最適化法を用いて同定する（ステップＳ１５３）。非線形最適化法の例としては、ＢＦＧＳアルゴリズム（Broyden Fletcher Goldfarb Shanno algorithm)を用いた、直線探索法を伴う準ニュートン法、または、信頼領域法を伴う修正ガウスニュートン法（Levenberg-Marquardt法とも呼ばれる）が挙げられる。大規模な問題に対しては、前者の方法の方がより強力である。加振点ｉにおける複素固有モードの実部および虚部は、具体的には、後述する式（５）の右辺第１項におけるｘｒｉおよびｙｒｉで表される。なお、加振点ｉにおける複素固有モードの同定に際し、例えば ｍ・（ｍ＋１）／２点における自己および相互周波数応答が参照される。
【0030】
第３同定部１１７は、周波数応答をリファレンスとして、線形連成系Ｌｃにおいて加振点ｉ以外の応答点ｌにおける複素固有モードを、線形最適化法を用いて同定する（ステップＳ１５４）。線形最適化法の例としては、ＱＲ分解法、または、特異値分解法が挙げられる。応答点ｌにおける複素固有モードの実部および虚部は、具体的には、後述する式（５）の右辺第１項におけるｘｒｌおよびｙｒｌで表される。なお、応答点ｌにおける複素固有モードの同定に際し、例えば ｍ・（ｎ－ｍ）点における相互周波数応答が参照される。
【0031】
なお、線形連成系Ｌｃの規模がそれほど大きくない場合、あるいは、固有モードを出力したい点がそれほど多くない場合は、ステップＳ１５１における周波数応答解析の計算コストがそれほど大きくならないため、ｎ＝ｍとしてもよい。この場合、ステップＳ１５４はスキップされる。
【0032】
ステップＳ１５１～Ｓ１５４によって一般粘性減衰系モデルのモーダルパラメータ（複素固有値および複素固有モード）が同定される。
【0033】
同定について、次の式（４）および式（５）を参照して説明する。
【数4】
【数5】
【0034】
同定とは、式（５）において左辺を、右辺第１項の 「モデルを用いて再合成した周波数応答」から右辺第２項の「リファレンスの周波数応答」を減じた残差、として表した場合に、式（４）の右辺で示されるように、当該残差の二乗和を最小とする値（実部、虚部）を求めることをいう。
【0035】
次に、本実施形態において一般粘性減衰系モデルにおいて同定されたモーダルパラメータが、どの程度までリファレンスに近づくか、について説明する。
【0036】
図６乃至図８は、異なる３点（point １～３）において、リファレンスとした周波数応答（実線）に対して、モード合成による周波数応答（○印）を示す図である。図６乃至図８では、周波数に対する位相および振幅の特性が左欄に、周波数に対する実部および虚部の特性が右欄に、それぞれ示される。
図６乃至図８に示されるように、リファレンスとした周波数応答に対して、モード合成による周波数応答が、ほぼ一致しており、高い精度で同定が成立していることが判る。
図６乃至図８に示されるように、リファレンスとした周波数応答に対して、モード合成による周波数応答が、ほぼ一致しており、高い精度で同定が成立していることが判る。
【0037】
図９乃至図１１は、上記３点（point １～３）において、リファレンスとした周波数応答のＩＦＦＴとしてのインパルス応答（実線）に対して、モード合成によるインパルス応答（○印）を示す図である。
図９乃至図１１に示されるように、リファレンスとしたインパルス応答に対して、モード合成によるインパルス応答が、ほぼ一致しており、高い精度で同定が成立していることが判る。
図９乃至図１１に示されるように、リファレンスとしたインパルス応答に対して、モード合成によるインパルス応答が、ほぼ一致しており、高い精度で同定が成立していることが判る。
【0038】
なお、実施形態では、「空気中に存在する場合の構造物」の物理モデルにおけるモーダルパラメータを同定したが、当該「構造物が付随する場合の空気」の物理モデルにおけるモーダルパラメータについても同様に同定することができる。
前者において、ステップＳ１５１の「空気中構造物の周波数応答解析」を実施するところを、後者においては、「構造物が付随する場合の空気の周波数応答解析」を実施する、という違いがあるのみである。ただし、前者において算出すべき周波数応答は「構造変位のフーリエ変換」／「構造に与える加振力のフーリエ変換」であるのに対し、後者において算出すべき周波数応答は「空気の圧力のフーリエ変換」／「空気の流量のフーリエ変換」となる。
前者において、ステップＳ１５１の「空気中構造物の周波数応答解析」を実施するところを、後者においては、「構造物が付随する場合の空気の周波数応答解析」を実施する、という違いがあるのみである。ただし、前者において算出すべき周波数応答は「構造変位のフーリエ変換」／「構造に与える加振力のフーリエ変換」であるのに対し、後者において算出すべき周波数応答は「空気の圧力のフーリエ変換」／「空気の流量のフーリエ変換」となる。
【0039】
本実施形態に係るモーダルパラメータの同定方法によれば、複素固有値の同定と複素固有モードの同定とが分離して実行されるので、加振点および応答点を多数として、「空気中に存在する大規模な構造物」の物理モデル、あるいは、「大規模な構造物が付随する空気」の物理モデルへの適用が容易になる。このため、従来手法では困難だった「空気中に置かれた大規模な構造物」、あるいは、「大規模な構造物が付随する空気」の動特性を、シミュレーションによって効率よく予測することが可能となる。
【0040】
このように、本実施形態によれば、振動や音響の特性がその価値を左右するような製品、あるいは、構造物と空気との連成がその性能を左右するような製品について、当該製品の試作前に、物理モデルを用いたシミュレーションによって動特性を予測することでき、原因と結果との因果関係が明確化されるので、当該製品の開発効率を大きく向上させることができる。
上記製品として例えばアコースティック楽器が一例として挙げられる。このようなアコースティック楽器の本体（響板、管体など）を製作するにあたり、当該本体の物理モデルの設計を変更し、変更した物理モデルを用いて楽音信号を合成することにより、試作前の楽器の楽音をシミュレーションにより聴くことができる。したがって、例えば新規な楽器を効率良く開発することができる。
上記製品として例えばアコースティック楽器が一例として挙げられる。このようなアコースティック楽器の本体（響板、管体など）を製作するにあたり、当該本体の物理モデルの設計を変更し、変更した物理モデルを用いて楽音信号を合成することにより、試作前の楽器の楽音をシミュレーションにより聴くことができる。したがって、例えば新規な楽器を効率良く開発することができる。
【0041】
＜変形例＞
以上に例示した実施形態は多様に変形され得る。実施形態に適用され得る具体的な変形の態様を以下に例示する。以下の例示から任意に選択された２以上の態様を、相互に矛盾しない範囲で併合してもよい。
以上に例示した実施形態は多様に変形され得る。実施形態に適用され得る具体的な変形の態様を以下に例示する。以下の例示から任意に選択された２以上の態様を、相互に矛盾しない範囲で併合してもよい。
【0042】
楽音信号を合成するコンピューター１０の機能は、上述した通り、制御装置１１を構成する単数または複数のプロセッサと、記憶装置１２に記憶されたプログラムとの協働により実現される。
このプログラムは、コンピューター１０が読取可能な記録媒体に格納された形態で提供されてコンピューターにインストールされ得る。記録媒体は、例えば非一過性（non-transitory）の記録媒体であり、ＣＤ－ＲＯＭ等の光学式記録媒体（光ディスク）が好例であるが、半導体記録媒体または磁気記録媒体等の公知の任意の形式の記録媒体も包含される。なお、非一過性の記録媒体とは、一過性の伝搬信号（transitory, propagating signal）を除く任意の記録媒体を含み、揮発性の記録媒体も除外されない。また、配信装置が通信網を介してプログラムを配信する構成では、当該配信装置においてプログラムを記憶する記憶装置１２が、前述の非一過性の記録媒体に相当する。
このプログラムは、コンピューター１０が読取可能な記録媒体に格納された形態で提供されてコンピューターにインストールされ得る。記録媒体は、例えば非一過性（non-transitory）の記録媒体であり、ＣＤ－ＲＯＭ等の光学式記録媒体（光ディスク）が好例であるが、半導体記録媒体または磁気記録媒体等の公知の任意の形式の記録媒体も包含される。なお、非一過性の記録媒体とは、一過性の伝搬信号（transitory, propagating signal）を除く任意の記録媒体を含み、揮発性の記録媒体も除外されない。また、配信装置が通信網を介してプログラムを配信する構成では、当該配信装置においてプログラムを記憶する記憶装置１２が、前述の非一過性の記録媒体に相当する。
【0043】
＜付記＞
以上の記載から、例えば以下のように本発明の好適な態様が把握される。なお、各態様の理解を容易にするために、以下では、図面の符号を便宜的に括弧書で併記するが、本発明を図示の態様に限定する趣旨ではない。
以上の記載から、例えば以下のように本発明の好適な態様が把握される。なお、各態様の理解を容易にするために、以下では、図面の符号を便宜的に括弧書で併記するが、本発明を図示の態様に限定する趣旨ではない。
【0044】
本開示のひとつの態様（態様１）に係る楽音合成方法は、コンピューター（１０）が、「空気中に存在する場合の構造物」の物理モデル、または、「構造物が付随する場合の空気」の物理モデルにおけるモーダルパラメータを同定するモーダルパラメータ同定方法であって、コンピューター（１０）が、空気中に存在する場合の構造物の周波数応答、または、構造物が付随する場合の空気の周波数応答を算出する第１過程と、コンピューター（１０）が、周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第２過程と、コンピューター（１０）が、周波数応答をリファレンスとして、構造物における加振点の複素固有モード、または、空気における加振点の複素固有モードを同定する第３過程と、を含む。
この態様１によれば、複素固有値の同定と複素固有モードの同定とを同時ではなく、分離して実行するので、加振点および応答点を多数として、空気中に存在する大規模な構造物の物理モデル、あるいは、大規模な構造物が付随する空気の物理モデルへの適用が容易となる。
この態様１によれば、複素固有値の同定と複素固有モードの同定とを同時ではなく、分離して実行するので、加振点および応答点を多数として、空気中に存在する大規模な構造物の物理モデル、あるいは、大規模な構造物が付随する空気の物理モデルへの適用が容易となる。
【0045】
態様１の具体例（態様２）において、第３過程の後に、コンピューター（１０）が、周波数応答をリファレンスとして、構造物における加振点とは異なる応答点の複素固有モード、または、空気における加振点とは異なる応答点の複素固有モードを同定する第４過程を、含む。
この態様２によれば、加振点を除いた応答点の数を多くして、より大規模な物理モデルへの適用が可能になる。
この態様２によれば、加振点を除いた応答点の数を多くして、より大規模な物理モデルへの適用が可能になる。
【0046】
また、態様１に係るモーダルパラメータの同定方法は、態様３に係るプログラムとして把握することが可能である。すなわち、態様３４に係るプログラムは、コンピュータ（１０）に、「空気中に存在する場合の構造物」の物理モデル、または、「構造物が付随する場合の空気」の物理モデルにおけるモーダルパラメータを同定させるプログラムであって、コンピュータ（１０）に、「空気中に存在する場合の構造物」の周波数応答、または、「構造物が付随する場合の空気」の周波数応答を算出する算出部（１１２）、周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第１同定部（１１５）、および、周波数応答をリファレンスとして、構造物における加振点の複素固有モード、または、空気における加振点の複素固有モードを同定する第２同定部（１１６）、として機能させる。
【0047】
また、態様１に係るモーダルパラメータの同定方法は、態様４に係るモーダルパラメータの同定装置として把握することが可能である。すなわち、態様４に係るモーダルパラメータの同定装置は、「空気中に存在する場合の構造物」の物理モデル、または、「構造物が付随する場合の空気」の物理モデルにおけるモーダルパラメータを同定するモーダルパラメータ同定装置であって、「空気中に存在する場合の構造物」の周波数応答、または、「構造物が付随する場合の空気」の周波数応答を算出する算出部（１１２）と、周波数応答の逆フーリエ変換をリファレンスとして複素固有値を同定する第１同定部（１１５）と、周波数応答をリファレンスとして、構造物における加振点の複素固有モード、または、空気における加振点の複素固有モードを同定する第２同定部（１１６）と、を含む。
【符号の説明】
【0048】
１０…コンピューター、１１…制御装置、１２…記憶装置、１３…演奏情報出力装置、１４…放音装置、１１０…処理制御部、１１１…前処理部、１１２…算出部、１１５…第１同定部、１１６…第２同定部、１１７…第３同定部、Ｌｃ…線形連成系。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】