JP5332296B2 - 楽音合成装置およびプログラム - Google Patents

Description

本発明は、リードの振動に応じて発音する管楽器の楽音を合成する技術に関する。

楽器の発音の原理を模擬することで楽音を合成する物理モデル方式の楽音合成装置（物理モデル音源）が従来から提案されている。例えば、全体が自由に移動する剛体の空気弁としてリードをモデル化したうえでクラリネットの挙動を模擬する技術（非特許文献１）や、ひとつの端部が固定された長板状の振動体（片持ち梁）でリードをモデル化したうえでクラリネットの挙動を模擬する技術（非特許文献２）が提案されている。
R.T.Schumacher, "Ab Initio Calculations of the Oscillations of a Clarinet", ACUSTICA, 1981, Volume 48 No.2, p.75- p.85 S.D.Sommerfeldt, W.J.Strong, "Simulation of a player-clarinet system", Acoustical Society of America, 1988, 83(5), p.1908- p.1918

ところで、実際の管楽器のリードは演奏者の唇に連成して複雑に挙動するのに対し、非特許文献１や非特許文献２の技術は、リードに対する単純な外力の作用を模擬しているに過ぎない。したがって、リードの挙動は忠実に再現されず、実際の管楽器の楽音に充分に近似する楽音を合成することは困難である。以上の事情を背景として、本発明は、演奏者の唇の作用を反映した忠実な楽音を合成することを目的としている。

以上の課題を解決するために、本発明に係る楽音合成装置は、吹奏時に唇に接触するリードの振動に応じて発音する管楽器の楽音を合成する装置であって、外力が唇に作用した平衡時におけるリードの挙動を表わす第１運動方程式（例えば図４における式Ａ1）と平衡時における唇の挙動を表わす第２運動方程式（例えば図４における式Ａ2）とを解くことで平衡時における唇の変位およびリードの変位を算定する第１演算手段と、第１演算手段による算定の結果を唇の変位およびリードの変位の初期値として、唇とリードとが相互に連成的に振動する連成振動を表現する運動方程式（例えば図４における式Ｂ）を解くことで、リードの変位を算定する第２演算手段と、第２演算手段が算定した変位に基づいて楽音を合成する楽音合成手段（例えば図４の第３演算部３１３および第４演算部３１４と図１の管体模擬部３３および伝達模擬部３５）とを具備する。

以上の構成においては、唇とリードとの連成振動の運動方程式に基づいてリードの変位が算定されるから、唇の挙動を反映しない運動方程式に基づいてリードの変位を算定する構成と比較してリードの挙動が正確に模擬される。したがって、実際の管楽器の楽音を忠実に再現することが可能となる。

本発明の好適な態様においては、唇に作用する外力の強度が変化するたびに、第１演算手段は、変化後の外力の強度に対応した唇の変位を第１運動方程式および第２運動方程式に基づいて算定し、第２演算手段は、第１演算手段が算定した唇の変位を連成振動の運動方程式に代入することでリードの変位を算定する。以上の態様によれば、唇に作用する外力の変化がリードの変位に反映されるから、唇に対する押圧力を変化させる奏法に対応した多彩な楽音を合成することが可能である。

本発明の好適な態様において、第１運動方程式および第２運動方程式は、唇における位置と押圧力の強度とに応じて変化する唇のバネ定数を含む。以上の態様によれば、バネ定数が押圧力の強度と位置に応じて変化するという現実の唇の特性が忠実に模擬されるから、管楽器の楽音を正確に合成することが可能となる。

本発明の好適な態様において、第１運動方程式は、リードの位置に応じて変化する曲げ剛性を含む。以上の態様においては、曲げ剛性（リードの断面二次モーメントとリードのヤング率との積）がリードの位置に応じて変化するという現実のリードの特性が忠実に模擬されるから、断面の形状が変化しない単純な長板状の振動体でリードを模擬する場合と比較して、管楽器の楽音を正確に合成することが可能となる。

本発明の好適な態様において、第２演算手段は、リードの変位を所定の範囲内に制限する。以上の態様においては、連成振動の運動方程式に基づいて算定されたリードの変位が所定の範囲内に制限されるから、現実のリードの変位の範囲を逸脱した位置までリードが変位するといった状況の模擬が防止される。したがって、現実の管楽器の楽音を正確に再現することが可能である。なお、リードの変位を制限する範囲としては、例えば、唇の底面からマウスピースの表面（リードとの対向面）までの範囲が好適である。

本発明の好適な態様において、連成振動の運動方程式は、唇における位置に応じて変化する当該唇の内部抵抗と、リードにおける位置に応じて変化する当該リードの内部抵抗との少なくとも一方を含む。以上の態様においては、唇およびリードの内部抵抗が位置に応じて変化する様子が模擬されるから、唇やリードの内部抵抗が固定値に設定された構成と比較して、実際の管楽器の楽音を忠実に再現することが可能となる。

なお、唇やリードの変形が小さい場合（より具体的には弾性限界内にある場合）、各々に作用する押圧力の強度が内部抵抗に与える影響は無視できる。しかし、唇やリードの変形が大きい場合（例えば弾性限界外にある場合）、唇やリードの内部抵抗は、各々における位置に加えて押圧力の強度によっても変化する。そこで、本発明の好適な態様において、連成振動の運動方程式は、唇における位置および当該唇に作用する押圧力の強度に応じて変化する当該唇の内部抵抗と、リードにおける位置および当該リードに作用する押圧力の強度に応じて変化する当該リードの内部抵抗との少なくとも一方を含む。以上の態様においては、唇およびリードの内部抵抗が位置や押圧力の強度に応じて変化する様子が模擬されるから、唇やリードの内部抵抗が位置のみに依存する数値や固定値に設定された構成と比較して、実際の管楽器の楽音を忠実に再現することが可能となる。

本発明に係る楽音合成装置は、各処理に専用されるＤＳＰ（Digital Signal Processor）などのハードウェア（電子回路）によって実現されるほか、ＣＰＵ（Central Processing Unit）などの汎用の演算処理装置とプログラムとの協働によっても実現される。本発明に係るプログラムは、吹奏時に唇に接触するリードの振動に応じて発音する管楽器の楽音を合成するためのプログラムであって、外力が唇に作用した平衡時におけるリードの挙動を表わす第１運動方程式と平衡時における唇の挙動を表わす第２運動方程式とを解くことで平衡時における唇の変位およびリードの変位を算定する第１演算処理と、第１演算処理による算定の結果を唇の変位およびリードの変位の初期値として、唇とリードとが相互に連成的に振動する連成振動を表現する運動方程式を解くことで、リードの変位を算定する第２演算処理と、第２演算処理で算定した変位に基づいて楽音を合成する楽音合成処理とをコンピュータに実行させる。以上のプログラムによっても、本発明に係る楽音合成装置と同様の作用および効果が奏される。なお、本発明のプログラムは、コンピュータが読取可能な記録媒体に格納された形態で利用者に提供されてコンピュータにインストールされるほか、通信網を介した配信の形態で提供されてコンピュータにインストールされる。

＜Ａ：第１実施形態＞
図１は、本発明の第１実施形態に係る楽音合成装置の構成を示すブロック図である。本形態の楽音合成装置１００は、サックスやクラリネットといったシングルリードの管楽器の発音の原理を演算によって模擬（シミュレート）することで楽音を合成する。図１に示すように、楽音合成装置１００は、演算処理装置１０と記憶装置４２と入力装置４４と放音装置４６とで構成されるコンピュータシステムで実現される。

演算処理装置（例えばＣＰＵ（Central Processing Unit））１０は、記憶装置４２に格納されたプログラムを実行することで、管楽器の楽音の時間波形（音圧の時間的な変動）を表す楽音データを生成および出力する。記憶装置４２は、演算処理装置１０が実行するプログラムや演算処理装置１０が使用するデータを記憶する。磁気記憶装置や半導体記憶装置など公知の記憶媒体が記憶装置４２として任意に採用される。

入力装置４４は、利用者が操作する複数の操作子で構成される。利用者は、楽音の合成に使用される様々なパラメータを入力装置４４から演算処理装置１０に指示することが可能である。キーボードまたはやマウスなどの入力機器や、管楽器の演奏に関する情報を入力するための楽器型の入力機器（例えばMIDI（Musical Instrument Digital Interface）コントローラ）が入力装置４４として採用される。

放音装置４６は、演算処理装置１０が出力する楽音データに応じた音波を放射する。なお、実際には楽音データをアナログの楽音信号に変換するＤ/Ａ変換器や楽音信号を増幅して放音装置４６に出力する増幅器が設置されるが、図１においてはＤ/Ａ変換器や増幅器の図示が省略されている。

図１の演算処理装置１０は、設定部１２および合成部１４として機能する。なお、演算処理装置１０の各機能が複数の集積回路で分散的に実現される構成も採用される。また、演算処理装置１０の機能の一部は、楽音の合成に専用される電子回路（ＤＳＰ）によって実現されてもよい。

設定部１２は、楽音の合成に必要なパラメータを設定する。合成部１４は、設定部１２が設定したパラメータに基づいて楽音データを生成する。合成部１４はリード模擬部３１と管体模擬部３３と伝達模擬部３５とを含む。リード模擬部３１は、演奏者の唇とリードとの連成振動を模擬する。管体模擬部３３は、管楽器のうちマウスピースからベルまでの管状の部分（すなわちリード以外の部分である。以下「管体部」という）の作用を模擬する。伝達模擬部３５は、ベルや各トーンホールからの放射音に対する伝達特性の付与を模擬する。

図２は、リード模擬部３１が模擬する管楽器のリードの近傍を示す概念図である。リードＭRは、ひとつの端部がマウスピースＭPに固定された長板状の振動体である。図２に示すように、リードＭRの先端部における幅方向の中心を原点としてＸ軸とＹ軸とＺ軸とを想定する。Ｚ軸はリードＭRの幅方向に延在する。Ｘ軸は、外力が作用しない状態におけるリードＭRの上面（マウスピースＭPとの対向面）内においてＺ軸に直交する。また、Ｙ軸は、Ｘ軸およびＺ軸に直交する（リードＭRに対して上下の方向に延在する）。

図３は、管楽器の吹奏時に演奏者の唇ＭLがリードＭRに接触する様子をＺ方向からみた模式図である。図３に示すように、リード模擬部３１は、管楽器の吹奏時に演奏者が唇ＭLを歯ＭTでリードＭRに押付けた状態を模擬する。唇ＭLは、リードＭRのうちＸ方向における位置ｘlip1（リードＭRの先端側）から位置ｘlip2（リードＭRの根元側）までの区間に接触する。また、演奏者の歯ＭTは、唇ＭLの下面のうちＸ方向における位置ｘteeth1（リードＭRの先端側）から位置ｘteeth2（リードＭRの根元側）までの区間に接触して押圧力ｆlip(x)を均等に作用させる。

図４は、リード模擬部３１の機能を示すブロック図である。図４の左側には、設定部１２が設定して記憶装置４２に格納したパラメータが羅列されている。各パラメータの意味を以下に説明する。

まず、リードＭRに関連するパラメータ（Stiff(x)，ｂreed(x)，Ａ(x)，μreed(x)，ρreed）を説明する。Stiff(x)は、Ｘ方向の位置ｘにおけるリードＭRの曲げ剛性［N・m²］である。すなわち、曲げ剛性Stiff(x)は、リードＭRのヤング率Ｅreed［Pa］と位置ｘにおけるリードＭRの断面二次モーメントＩ(x)［m⁴］との乗算値に相当する。図２に示すように、ｂreed(x)は、位置ｘにおけるリードＭRの横幅（Ｚ方向の寸法）［m］であり、Ａ(x)は、位置ｘにおけるリードＭRの断面積（位置ｘを通過するＹＺ平面内の面積）［m²］である。本形態ではリードＭRの断面の形状がＸ方向の位置ｘに応じて変化する。したがって、曲げ剛性Stiff(x)の計算に使用される断面二次モーメントＩ(x)とリードＭRの横幅ｂreed(x)と断面積Ａ(x)とは位置ｘの関数となる。また、図２のμreed(x)は、リードＭRの内部抵抗の分布［(kg/sec)/m］であり、ρreedはリードＭRの密度［kg/m³］である。

次に、唇ＭLに関連するパラメータ（ｋlip(x)，ｄlip(x)，μlip(x)，ｍlip(x)）を説明する。ｋlip(x)は、唇ＭLのＸ方向におけるバネ定数の分布［N/m²］（例えばＸ方向の単位長あたりのバネ定数）である。ｄlip(x)は、外力が作用しない場合の位置ｘにおける唇ＭLのＹ方向の寸法（厚さ）［m］である。μlip(x)は、位置ｘにおける唇ＭLの内部抵抗の分布［(kg/sec)/m］である。ｍlip(x)は、唇ＭLのＸ方向における質量の分布［kg/m］（例えばＸ方向の単位長あたりの質量）である。バネ定数の分布ｋlip(x)と厚さｄlip(x)と内部抵抗の分布μlip(x)と質量の分布ｍlip(x)とは位置ｘに応じて変化する。

図４のＰは演奏者の口腔内の圧力［Pa］であり、ρairは常温（例えば２５℃）における空気の密度［kg/m³］である。Ｈ(x)は、図２に示すように、マウスピースＭPのうちリードＭRに対向する表面のＹ方向における位置（以下「フェーシング位置」という）である。リードＭRのＹ方向の変位ｙ(x,t)がフェーシング位置Ｈ(x)に到達するとリードＭRの上面はマウスピースＭPに接触するから、フェーシング位置Ｈ(x)はリードＭRの変位の限界値（下限値）に相当する。また、Ｚcは、マウスピースＭPのうち管体とみなせる部分の始点（リードＭRの根元）における空気の流動に対する特性インピーダンスである。

図４に示すように、リード模擬部３１は、第１演算部３１１と第２演算部３１２と第３演算部３１３と第４演算部３１４とで構成される。第１演算部３１１は、唇ＭLのうちＸ方向における位置ｘfに押圧力ｆlip(xf)を静的に作用させて平衡した場合のリードＭRの変位ｙ0(xf)と唇ＭLの底面の変位ｙb(xf)とを算定する。第２演算部３１２は、第１演算部３１１が算定した変位ｙ0(xf)および変位ｙb(xf)をリードＭRおよび唇ＭLの底面の変位の初期値（ｔ＝０における数値）として唇ＭLとリードＭRとの連成振動の運動方程式を解くことで、Ｘ方向におけるリードＭRの各位置ｘでの時刻ｔにおけるＹ方向の変位ｙ(x,t)を算定する。第３演算部３１３および第４演算部３１４は、リードＭRから管体部（マウスピースＭP側）に出力する音波の圧力ＰOUTをリードＭRの変位ｙ(x,t)に基づいて算定する。リード模擬部３１による処理の詳細を以下に説明する。

図３に示すように、演奏者の唇ＭLの位置ｘf（ｘteeth1≦ｘf≦ｘteeth2）に歯ＭTから押圧力ｆlip(xf)を作用させて平衡した状態を想定する。押圧力ｆlip(xf)の作用でリードＭRがＹ方向に距離ｄ1だけ変形するとともに唇ＭLがＹ方向に距離ｄ2だけ変形したとすると、リードＭRから唇ＭLに作用する弾性力Ｒ1と唇ＭLからリードＭRに作用する弾性力Ｒ2とは、リードＭRの曲げ剛性Stiff(xf)と唇ＭLのバネ定数ｋlip(xf)とを含む以下の各式で表現される。なお、実際には唇ＭLの上面はリードＭRの下面に接触するが、図３においては唇ＭLの上面がリードＭRの上面に位置するように単純化されている。

リードＭRと唇ＭLとの接触点（位置ｘf）における力の釣合いから、
Ｒ1−Ｒ2＝０
が成立し、唇ＭLと歯ＭTとの接触点（位置ｘf）における力の釣合いから、
ｆlip(xf)＝Ｒ2
が成立する。また、リードＭRの変形と変位との関係から、
ｄ1＝ｙ0(xf)
が成立し、唇ＭLの変形と変位との関係から、
ｄ2＝｛ｙb(xf)−ｄlip(xf)｝−ｙ0(xf)
が成立する。

以上の各式から以下の運動方程式Ａ1および運動方程式Ａ2が導出される。

図４の第１演算部３１１は、設定部１２が設定した曲げ剛性Stiff(xf)と押圧力ｆlip(xf)とバネ定数ｋlip(xf)と厚さｄlip(xf)とを代入して運動方程式Ａ1と運動方程式Ａ2との連立方程式を解くことで唇ＭLの底面の変位ｙb(xf)とリードＭRの変位ｙ0(xf)とを算定する。さらに詳述すると、第１演算部３１１は、差分方程式化やGaussの消去法などを利用して運動方程式Ａ1からリードＭRの変位ｙ0(xf)を算定し、当該変位ｙ0(xf)を運動方程式Ａ2に代入することで唇ＭLの変位ｙb(xf)を算定する。なお、運動方程式Ａ1の解法については後述する。

演奏者が管楽器を吹奏することで唇ＭLとリードＭRとが連成的に振動したときの動特性は以下の運動方程式Ｂで表現される。

第２演算部３１２は、第１演算部３１１が算定した変位ｙ0(xf)を運動方程式ＢにおけるリードＭRの変位ｙ(x,t)の初期値に設定するとともに第１演算部３１１が算定した変位ｙb(xf)を運動方程式Ｂにおける唇ＭLの変位ｙb(x)に代入して運動方程式Ｂを解くことで、リードＭRの変位ｙ(x,t)を算定する。運動方程式Ｂの右辺はＸ方向におけるリードＭRの位置ｘに作用する外部力ｆex(x)に相当する。第２演算部３１２は、第１に、設定部１２が設定した各パラメータ（ｂreed(x)，Ｐ，ｋlip(x)，ｄlip(x)）と第４演算部３１４が算定した圧力ｐ(t)とを運動方程式Ｂの右辺に代入するとともに、第１演算部３１１が算定した変位ｙ0(xf)と変位ｙb(xf)とを運動方程式Ｂの右辺における変位ｙ(x,t)および変位ｙb(x)の初期値として代入することで外部力ｆex(x)を算定する。圧力ｐ(t)は、リードＭRとマウスピースＭPとの間隙の空間のうちリードＭRの先端の近傍（以下「リード直上部」という）における圧力を意味する。なお、第４演算部３１４による圧力ｐ(t)の算定については後述する。

第２に、第２演算部３１２は、設定部１２が設定した各パラメータ（ｍlip(x)，Ａ(x)，μlip(x)，μreed(x)，Stiff(x)，ρreed）を運動方程式Ｂの左辺に代入するとともに先に演算した外部力ｆex(x)を運動方程式Ｂの右辺に設定することでリードＭRの変位ｙ(x,t)を算定する。運動方程式Ｂを解く具体的な方法を以下に例示する。

運動方程式Ｂの左辺における第２項は以下のように変形される。

したがって、運動方程式Ｂは以下の式Ｂ1に変形される。

次に、時間ｔを整数ｉと所定値Δｔとの乗算値として離散化（ｔ＝ｉ・Δｔ）したうえで、時間微分を以下の差分に置換する。

また、図５に示すように、相互に等しい間隔Δｘをあけて分布するようにＸ方向における位置ｘを離散化する。すなわち、位置ｘを整数ｎと所定値Δｘとの乗算値として離散化（ｘ＝ｎ・Δｘ）したうえで、位置微分を以下の差分に置換する。

なお、以上におけるｙ(n,i)は、ｙ(n・Δx，i・Δt)を略記した記号である。

したがって、式Ｂ1は以下の式Ｂ2のように差分方程式化される。

ただし、式Ｂ2においては各項が以下のように置換されている。

また、式Ｂ2の各文字に付加された記号(n，i)は(n・Δx，i・Δt)の略記である。

次に、式Ｂ2における左辺の第２項目から第４項目までに１／２を乗算した方程式と、式Ｂ2のｉを(ｉ＋１)に置換したうえで左辺の第２項目から第４項目までに１／２を乗算した方程式とを加算することで、式Ｂ2を近似的に表現する式Ｂ3が導出される。

式Ｂ3の各項を変数ｙの種類毎に整理して変形すると以下の式Ｂ4が導出される。

ただし、式Ｂ4においては各項が以下のように置換されている。

図５に示すように、リードＭRが位置ＮにてマウスピースＭPに固定されるとすれば、ｙ(N,i)やｙ(N+1,i)は任意の時点ｉにてゼロとなる。また、図５に示すように、外力が作用しないリードＭRの先端（ｎ＝０）においては加速度（∂²y(0,i)／∂x²）およびせん断力（∂³y(0,i)／∂x³）がゼロとなるから、以下の式Ｂ4_1および式Ｂ4_2が成立する。

さらに、式Ｂ4_1と式Ｂ4_2を加算することで以下の式Ｂ4_3が導出され、式Ｂ4_3の３倍から式Ｂ4_2を減算することで以下の式Ｂ4_4が導出される。
０・ｙ(0,i)＋ｙ(1,i)−２ｙ(2,i)＋ｙ(3,i)＝０ ……Ｂ4_3
ｙ(0,i)＋０・ｙ(1,i)−３ｙ(2,i)＋２ｙ(3,i)＝０ ……Ｂ4_4

また、式Ｂ4のｎに２を代入すると以下の式Ｂ4_5が導出される。

ｎ＝３〜Ｎ−１を同様に式Ｂ4に代入して導出される式と前述の式Ｂ4_3および式Ｂ4_4とから以下の式Ｂ5が導出される。

式Ｂ5の解法としてはGaussの消去法が好適である。なお、式Ｂ4_1および式Ｂ4_2から式Ｂ4_3および式Ｂ4_4を導出することで式Ｂ5の左上部の２行２列は対角行列となるから、Gaussの消去法における演算量が削減されるという利点がある。

第２演算部３１２は、第１演算部３１１による算定の結果（ｙ0(xf)，ｙb(xf)）を変位ｙ(x,y)および変位ｙb(x)の初期値として式Ｂ5を解くことでリードＭRの変位ｙ(x,t)を算定する。さらに詳述すると、第１に、第２演算部３１２は、式Ｂ5の右辺のうち現在の変位に相当する変数ｙ(0,i)〜ｙ(N-1,i)および式B5の右辺のうち過去の変位に相当する変数ｙ(2,i-1)〜ｙ(N-1,i-1)の双方に第１演算部３１１が算定したｙ0(0)〜ｙ0(N-1)およびｙ0(2)〜ｙ0(N-1)を代入するとともに式Ｂ5のｙb(2)〜ｙb(N-1)に第１演算部３１１が算定したｙb(xf)を代入して式Ｂ5を解くことで、式Ｂ5の左辺における未来の変位に相当する変数ｙ(0,i+1)〜ｙ(N-1,i+1)を算定する。第２に、第２演算部３１２は、時間をΔｔだけ進めるために、現在の変位に相当する変数ｙ(2,i)〜ｙ(N-1,i)を式Ｂ5の右辺のうち過去の変位に相当する変数ｙ(2,i-1)〜ｙ(N-1,i-1)に代入するとともに、直前に算定した未来の変位に相当する変数ｙ(0,i+1)〜ｙ(N-1,i+1)を式Ｂ5の右辺うち現在の変位に相当する変数ｙ(0,i)〜ｙ(N-1,i)に代入して式Ｂ5を解くことで、式Ｂ5の左辺のうち未来の変位に相当する変数ｙ(0,i+1)〜ｙ(N-1,i+1)を算定する。以上のように時点ｉでの変位ｙ(0,i)〜ｙ(N-1,i)と時点(i-1)での変位ｙ(2,i-1)〜ｙ(N-1,i-1)を代入したうえで式Ｂ5を解いて時点(i+1)での変位ｙ(0,i+1)〜ｙ(N-1,i+1)を算定するという演算を反復することで、第２演算部３１２は、リードＭRの各位置ｘにおける変位ｙ(x,t)の経時的な変化を算定する。

また、設定部１２の設定する押圧力ｆlip(x)が変化するたびに、第１演算部３１１は、変化後の押圧力ｆlip(x)を運動方程式Ａ1および運動方程式Ａ2における押圧力ｆlip(xf)に代入することで新たなｙ0(xf)およびｙb(xf)を算定する。第２演算部３１２は、第１演算部３１１が変位ｙb(xf)を算定するたびに、式Ｂ5のｙb(2)〜ｙb(N-1)に代入する数値を新たな変位ｙb(xf)に更新する。以上の構成によれば、押圧力ｆlip(x)を任意に変化させる演奏法を忠実に再現した楽音を合成することが可能である。一方、第２演算部３１２は、押圧力ｆlip(x)の変化時に第１演算部３１１が新たな変位ｙ0(xf)を算定しても、式Ｂ5における変位ｙ(0,i)〜ｙ(N-1,i)については、第１演算部３１１が算定した変位ｙ0(xf)を反映させない。以上の構成によれば、変位ｙ(x,t)の不連続な変化が回避されるから、聴感上において自然な楽音を生成することが可能である。

図４に示すように、第２演算部３１２は、リードＭRの変位ｙ(x、t)を所定の範囲内に制限する範囲制限部３２を含む。範囲制限部３２は、式Ｂ5から算定されたリードＭRの変位ｙ(x,t)を、第１演算部３１１が算定した唇ＭLの変位ｙb(xf)（唇ＭLのうち歯ＭTが接触する底面の位置）から、設定部１２が設定したフェーシング位置Ｈ(x)までの範囲に制限する。すなわち、範囲制限部３２は、変位ｙ(x,t)が変位ｙb(xf)を超える（Ｙ軸を下向きに正として変位ｙb(xf)を上回る）場合には変位ｙ(x,t)を変位ｙb(xf)に変更するとともに、変位ｙ(x,t)がフェーシング位置Ｈ(x)を超える（下回る）場合には変位ｙ(x,t)をフェーシング位置Ｈ(x)に変更する。以上の構成によれば、リードＭRが唇ＭLの底面よりも下やマウスピースＭPよりも上に位置するといった不条理な状況が模擬されることは回避される。なお、以上においては変位ｙ(x,t)の限界値（上限値）を唇ＭLの底面の変位ｙb(x)としたが、唇ＭLには実際には厚さがあるから、唇ＭLの厚さに相当する所定値（唇ＭLの厚さの最小値に相当する固定値や、押圧力ｆlip(x)に応じて変化する唇ＭLの厚さの最小値に相当する可変値）だけ変位ｙb(x)からフェーシング位置Ｈ(x)に近い位置を変位ｙ(x,t)の限界値（上限値）としてもよい。

なお、第１演算部３１１による変位ｙ0(x)の算定（運動方程式Ａ1の解法）には、以下に概説するように、第２演算部３１２による変位ｙ(x,t)と同様の方法が利用される。まず、式Ａ1は、運動方程式Ｂ1から式Ｂ2の変形と同様の方法で以下の式Ａ1_1のように差分方程式化される。

式Ａ1_1の各項を変数ｙの種類毎に整理して変形すると以下の式Ａ1_2が導出される。

ただし、式Ａ1_2においては各項が以下のように置換されている。

式Ｂ4から式Ｂ5の変形と同様の方法で、式Ａ1_2は以下の式Ａ1_3に変形される。

第１演算部３１１は、Gaussの消去法などの解法を利用して式Ａ1_3を解くことで変位ｙ0(x)（式Ａ1_3におけるｙ(0)〜ｙ(N-1)）を算定する。以上が運動方程式Ａ1の解法の具体例である。

図４の第３演算部３１３は、設定部１２が設定した各パラメータ（Ｈ(x)，ρair，ｂreed(x)，Ｚc）と第２演算部３１２が算定した変位ｙ(x,t)とに基づいてリード直上部における体積流速ｆ(t)を算定する。本形態の第３演算部３１３は、リードＭRの上面と下面との圧力差に起因して発生する体積流速Ｕ(t)と、リードＭRの各部が変位（ｙ(x,t)）することで発生する体積流速ｕ(t)との差分値をリード直上部の体積流速ｆ(t)として算定する（ｆ(t)＝Ｕ(t)−ｕ(t)）。

体積流速ｕ(t)は以下の式Ｃ1で表現される。なお、式Ｃ1におけるｌeffは、リードＭRの先端から支点までの距離（リードＭRの有効長）である。

第３演算部３１３は、設定部１２が設定したリードＭRの横幅ｂreed(x)と第２演算部３１２が算定した変位ｙ(x,t)の時間微分（すなわちリードＭRの速度）とを式Ｃ1に代入してSimpson法などの数値積分を実行することで体積流速ｕ(t)を算定する。

また、体積流速Ｕ(t)は以下の手順で算定される。まず、第３演算部３１３は、リードＭRの先端におけるマウスピースＭPとリードＭRとの間隔ξ(t)［m］を算定する。間隔ξ(t)は、第２演算部３１２が算定したリードＭRの変位ｙ(x,t)のうちリードＭRの先端（ｘ＝０）における変位ｙ(0,t)とリードＭRの先端（ｘ＝０）におけるフェーシング位置Ｈ(0)との差分値（ξ(t)＝ｙ(0,t)−Ｈ(0)）として算定される。

次いで、第３演算部３１３は、リードＭRの先端におけるマウスピースＭPとリードＭRとの間隙を通過する空気の有効質量Ｍ(t)［kg］を算定する。有効質量Ｍ(t)は以下の式Ｃ2で表現される。

式Ｃ2のＲ(t)は、リードＭRの先端における横幅ｂreed(0)と間隔ξ(t)との相対比（Ｒ(t)＝ｂreed(0)／ξ(t)）である。第３演算部３１３は、設定部１２が設定したリードＭRの横幅ｂreed(0)および空気の密度ρairと相対比Ｒ(t)とを式Ｃ2に代入することで有効質量Ｍ(t)を算定する。

有効質量Ｍ(t)と体積流速Ｕ(t)とについては以下の式Ｃ3が成立する（例えば非特許文献１参照）。第３演算部３１３は、式Ｃ3を解くことで体積流速Ｕ(t)を算定する。

式Ｃ3のＡは、所定の係数（例えばＡ＝0.0797）である。式Ｃ3を利用した体積流速Ｕ(t)の算定には例えば以下の方法が採用される。

式Ｃ3は、後述する式Ｄ1および式Ｄ2を利用して以下の式Ｃ4に変形される。

式Ｃ4における微分を後退差分で離散化すると以下の式Ｃ5が導出される。第３演算部３１３は、非線形方程式の数値解法（例えばNewton Raphson法）を利用することで式Ｃ5から体積流速Ｕ(t)を算定する。

第３演算部３１３は、図４に示すように、以上の手順で算定した体積流速Ｕ(t)と体積流速ｕ(t)との差分値を体積流速ｆ(t)として算定する。

図４の第４演算部３１４は、出射波圧力ＰOUT(t)とリード直上部の音圧ｐ(t)とを算定する。出射波圧力ＰOUT(t)は、リードＭRから管体部内に進行する音波（以下「出射波」という）の圧力である。管体部内を進行した音波の一部（以下「反射波」という）は管楽器の開放端（ベル）にて反射することで逆方向に進行してマウスピースＭP内に到達する。したがって、出射波圧力ＰOUT(t)は、体積流速ｆ(t)によって発生する圧力と管体内からリードＭR側に進行する反射波の圧力（以下「反射波圧力」という）ＰIN(t)との加算に相当する。なお、反射波圧力ＰIN(t)は管体模擬部３３が算定する。

体積流速ｆ(t)による圧力は体積流速ｆ(t)と特性インピーダンスＺcとの乗算値であるから、出射波圧力ＰOUT(t)は以下の式Ｄ1で表現される。
ＰOUT(t)＝Ｚc・ｆ(t)＋ＰIN(t) ……Ｄ1
第４演算部３１４は、設定部１２が設定した特性インピーダンスＺcと第３演算部３１３が算定した体積流速ｆ(t)と管体模擬部３３が算定した反射波圧力ＰIN(t)とを式Ｄ1に代入することで出射波圧力ＰOUT(t)を算定する。

また、リード直上部には出射波圧力ＰOUT(t)と反射波圧力ＰIN(t)とが作用するから、リード直上部の圧力ｐ(t)は以下の式Ｄ2で表現される。
ｐ(t)＝ＰOUT(t)＋ＰIN(t) ……Ｄ2
第４演算部３１４は、式Ｄ1に基づいて算定した反射波圧力ＰOUT(t)と管体模擬部３３が算定した反射波圧力ＰIN(t)とを式Ｄ2に代入することで圧力ｐ(t)を算定する。第４演算部３１４が算定した圧力ｐ(t)は、第２演算部３１２による外部力ｆex(x)の算定（式Ｂ）や第３演算部３１３による体積流速Ｕ(t)の演算（式Ｃ3）にフィードバックされる。

次に、管体模擬部３３の機能について説明する。図６に示すように、管楽器の実際の管体部（マウスピースからベルまでの部分）は、ｋ個の管状の単位部Ｕ（Ｕ[1]〜Ｕ[k]）を直列に連結した構造体で近似される（ｋは自然数）。各単位部Ｕの直径および全長（すなわち管体部の形状）は可変に設定される。管体模擬部３３は、図６の構造体を模擬する物理モデル（以下「管体モデル」という）を利用して、管体部の内側における音波の挙動を実現する。

図７は、管体模擬部３３が使用する管体モデルの構成を示すブロック図である。図７に示すように、管体モデルは、経路ｒ1上に単位部Ｕ毎に配置された遅延素子ＤA（ＤA[1]〜ＤA[k]）と、経路ｒ2上に単位部Ｕ毎に配置された遅延素子ＤB（ＤB[1]〜ＤB[k]）と、相隣接する遅延素子ＤA間および相隣接する遅延素子ＤB間に配置された接続部Ｊ（Ｊ[1]〜Ｊ[k-1]）と、(ｋ-1)個の接続部Ｊ[1]〜Ｊ[k-1]のうち管楽器のトーンホールに対応した位置の接続部Ｊに連結されたホール部ＴH（ＴH[1]〜ＴH[k-1]）と、管楽器のベルに対応するベル部ＢLとを含んで構成される。経路ｒ1は、管体部の内側でマウスピースＭPからベル側に進行する出射波の挙動（出射波圧力ＰOUT(k,t)）を模擬し、経路ｒ2は、管体部の内側でベルからマウスピースＭP側に進行する反射波の挙動（反射波圧力ＰIN(k,t)）を模擬する。

第ｉ段目（ｉ＝１〜ｋ）の遅延素子ＤA[i]は、前段から供給される出射波圧力ＰOUT(i,t)を遅延量ｄA[i]だけ遅延させる素子（例えば遅延量ｄA[i]に応じて段数が変化するシフトレジスタ）である。リード模擬部３１（第４演算部３１４）が算定した出射波圧力ＰOUT(t)は初期値ＰOUT(1,t)として第１段目の遅延素子ＤA[1]に供給され、各段の遅延素子ＤA[1]〜ＤA[k]にて順次に遅延が付与されたうえでベル部ＢLに到達する。すなわち、遅延素子ＤA[i]は第ｉ番目の単位部Ｕ[i]における出射波圧力ＰOUT(i,t)の伝播遅延を模擬する。

ベル部ＢLは、管楽器のベルからの音波の放射とベルの先端における音波の反射とを模擬する。図８に示すように、ベル部ＢLは、フィルタ部６２と乗算部６４とを含む。経路ｒ1の第ｋ段目（最終段）の遅延素子ＤA[k]から出力された出射波圧力ＰOUT(k,t)がベル部ＢLに入力される。フィルタ部６２は、ローパスフィルタ部６２１と減算部６２２とで構成される。ローパスフィルタ部６２１は、遅延素子ＤA[k]の出力に相当する出射波圧力ＰOUT(k,t)の時間波形において遮断周波数ｆCBを上回る成分を抑制する。ローパスフィルタ６２１内の乗算器による乗算値ＣBは、ＣB＝２π・ｆCB・Δｔを満たす。減算部６２２は、遅延素子ＤA[k]が出力する出射波圧力ＰOUT(k,t)からローパスフィルタ部６２１の出力を減算することで放射音圧ＰB(t)を算定する。すなわち、減算部６２２は、出射波圧力ＰOUT(k,t)のうち遮断周波数ｆCBを下回る成分を抑制するハイパスフィルタとして機能する。放射音圧ＰB(t)は、ベルから放射される音波の圧力に相当する。

乗算部６４は、管楽器のベルの内外の境界における音波の反射を模擬する。すなわち、乗算部６４は、ローパスフィルタ部６２１からの出力に係数ｒBを乗算することで反射波圧力ＰIN(k,t)を算定して経路ｒ2（図７の遅延素子ＤB[k]）に出力する。反射時に音波の位相が逆転するとともに損失が発生するから、係数ｒBは、例えば絶対値が１よりも小さい負数に設定される。

図７の遅延素子ＤB[i]は、遅延素子ＤA[i]と同様に、前段（ベル部ＢL側）から入力される反射波圧力ＰIN(i,t)を遅延量ｄB[i]だけ遅延させる。すなわち、遅延素子ＤB[i]は、単位部Ｕ[i]における反射波圧力ＰIN(i,t)の伝播遅延を模擬する。ベル部ＢLが算定した反射波圧力ＰIN(k,t)は、遅延素子ＤB[k]〜ＤB[1]にて順次に遅延が付与され、第１段目の遅延素子ＤB[1]から出力された反射波圧力ＰIN(1,t)が反射波圧力ＰIN(t)としてリード模擬部３１（第４演算部３１４）の演算に利用される。

接続部（ジャンクション）Ｊは、管体部の内径の変化に起因した出射波の拡散やエネルギの損失を模擬する。接続部Ｊには、図９の部分(A)の２ポート型と図９の部分(B)の３ポート型とが利用される。２ポート型の接続部Ｊ[i]は、経路ｒ1から入力される出射波圧力ＰOUT(i,t)に係数αiを乗算する乗算部７１と、経路ｒ2から入力される反射波圧力ＰIN(i+1,t)に係数βiを乗算する乗算部７２と、乗算部７１の出力（αi・ＰOUT(i,t)）と乗算部７２の出力（βi・ＰIN(i+1,t)）とを加算する加算部７３と、加算部７３からの出力と出射波圧力ＰOUT(i,t)との差分を新たな反射波出力ＰIN(i,t)として経路ｒ2に出力する減算部７４と、加算部７３からの出力と反射波圧力ＰIN(i+1,t)との差分を新たな出射波圧力ＰOUT(i+1,t)として経路ｒ1に出力する減算部７５とで構成される。以上の２ポート型の接続部Ｊ[i]はトーンホール部ＴHが連結されない部分（例えば図７における接続部Ｊ[1]や接続部Ｊ[2]）に適用される。

一方、図９の部分(B)における３ポート型の接続部Ｊ[i]は、トーンホール部ＴHが連結される部分（例えば図７における接続部Ｊ[3]や接続部Ｊ[4]に採用される。３ポート型の接続部Ｊ[i]は、２ポート型の接続部Ｊ[i]の要素に加えて、加算部７３の出力と第ｉ段目のトーンホール部ＴH[i]が出力する音圧Ｒi(t)との差分を音圧Ｑi(t)として当該トーンホール部ＴH[i]に出力する減算部７６と、音圧Ｒi(t)に係数γiを乗算して加算部７３に出力する乗算部７７とを含む。

トーンホール部ＴH[i]は、第ｉ番目のトーンホールからの音波の放射と当該トーンホールにおける音波の反射とを模擬する。図１０に示すように、トーンホール部ＴH[i]は、図８のベル部ＢLと同様に、遅延素子ＤE1と遅延素子ＤE2とフィルタ部６６と乗算部６８とで構成される。遅延素子ＤE1は、３ポート型の接続部Ｊ[i]から供給される音圧Ｑi(t)を遅延量ｄE1だけ遅延させる。フィルタ部６６は、遅延後の音圧Ｑi(t)のうち遮断周波数ｆCTHを上回る成分を抑制するローパスフィルタ部６６１と、音圧Ｑi(t)からローパスフィルタ部６６１の出力を減算することで放射音圧ＰHi(t)を算定する減算部６６２（ハイパスフィルタ）とで構成される。ローパスフィルタ６６１内の乗算器による乗算値ＣTHは、ＣTH＝２π・ｆCTH・Δｔを満たす。放射音圧ＰHi(t)は、第ｉ番目のトーンホールから外部に放射される音波の圧力に相当する。一方、第ｉ番目のトーンホールの閉塞時に音波の位相が反転しない様子や、トーンホールの開放時に音波の損失や位相の反転が発生する様子が模擬されるように、乗算部６８は、ローパスフィルタ部６６１からの出力に係数ｒHi（例えば絶対値が１を下回る正数や負数）を乗算することで音圧Ｒi(t)を算定する。すなわち、乗算部６８は、トーンホールの内外の境界における音波の反射を模擬する。音圧Ｒi(t)は、遅延素子ＤE2にて遅延量ｄE2だけ遅延されたうえで３ポート側の接続部Ｊ[i]（乗算部７７）に出力される。以上が管体模擬部３３の機能である。

図１の伝達模擬部３５は、管楽器のベルや各トーンホールからの放射音に対する伝達特性の付与を模擬する。図１１に示すように、伝達模擬部３５は、ベルに対応する乗算部３５１と、単位部Ｕ[1]〜Ｕ[k]に対応するｋ個の乗算部３５３と、乗算部３５１の出力とｋ個の乗算部３５３の出力とを加算することで受聴音圧Ｐmix(t)を算定する加算部３５５とで構成される。乗算部３５１は、ベル部ＢLが算定した放射音圧ＰB(t)に係数ＭBを乗算する。第ｉ番目の乗算部３５３は、トーンホール部ＴH[i]が算定した放射音圧ＰHi(t)に係数ＭHiを乗算する。係数ＭHiは、第ｉ番目のトーンホールが閉塞されている場合や第ｉ番目のトーンホールが管楽器に存在しない場合にはゼロに設定されるとともに第ｉ番目のトーンホールが開放されている場合にはゼロを上回る所定値（例えば１）に設定される。したがって、加算部３５５が算定する受聴音圧Ｐmix(t)は、ベルからの放射音と演奏者が開放しているトーンホールからの放射音とを混合した音波（受聴音）の音圧となる。受聴音圧Ｐmix(t)が楽音データとして演算処理装置１０から放音装置４６に出力される。

次に、設定部１２について説明する。図１に示すように、設定部１２は、特性パラメータ変換部２１と形状パラメータ変換部２３とで構成される。特性パラメータ変換部２１は、リードＭRや唇ＭLの物性に関する様々なパラメータを楽音の合成に必要なパラメータに変換する。また、形状パラメータ変換部２３は、管楽器の形状や寸法に関する様々なパラメータを楽音の合成に必要なパラメータに変換する。

図１２は、特性パラメータ変換部２１の具体的な機能を示すブロック図である。利用者は、入力装置４４を適宜に操作することで、図１２の左側に羅列された様々なパラメータを演算処理装置１０に指示する。利用者が指示するパラメータは、空気に関する物性値（ｃair，ρair）と、唇ＭLに関する物性値（ρlip，Ｅlip，tanδlip）と、特定の唇のサンプル（以下「唇サンプル」という）に関する寸法（ｂlip_sample）と、リードＭRに関する物性値（ρreed，Ｅreed，tanδreed）と、特定のリードのサンプル（以下「リードサンプル」という）の寸法（ｂreed_sample，ｌreed_sample，ｄreed_sample）と、息圧Ｐ0および音高ｆnとを含む。

ｃairは空気中における音速［m/sec］であり、ρairは空気の密度［kg/m³］である。息圧Ｐ0は、吹奏時の利用者の口内における空気の圧力である。また、音高ｆnは、演算処理装置１０が合成すべき楽音の高低を示す数値である。音高ｆnを適宜に変化させることで楽曲の演奏音を合成することが可能である。

唇ＭLに関する物性値は、唇ＭLの密度ρlip［kg/m³］と唇ＭLのヤング率Ｅlip［Pa］と唇ＭLの損失係数tanδlipとを含む。唇サンプルの特性値は、横幅（Ｚ方向における寸法）ｂlip_sample［m］を含む。唇サンプルは、実際の人間の唇と同等の物性を持つ材質で形状のみを簡素な立体（本形態では直方体）に単純化した構造体である。したがって、横幅ｂlip_sampleはＸ方向の位置に依存しない固定値である。なお、唇ＭLや唇サンプルに関する各物性値や寸法を利用者が個別に入力する構成のほか、以上の各パラメータ（ρlip，Ｅlip，tanδlip，ｂlip_sample）の数値を複数種の唇ＭLの各々について事前に記憶装置４２に格納しておき、利用者が入力装置４４で選択した種類の唇ＭLに関する各パラメータを特性パラメータ変換部２１が記憶装置４２から取得する構成も好適に採用される。

リードＭRに関する物性値は、リードＭRの密度ρreed［kg/m³］とリードＭRのヤング率Ｅreed［Pa］とリードＭRの損失係数tanδreedとを含む。リードサンプルの特性値は、横幅（Ｚ方向における寸法）ｂreed_sample［m］と長さ（Ｘ方向における寸法）ｌreed_sample［m］と厚さ（Ｙ方向における寸法）ｄreed_sample［m］とを含む。リードサンプルは、実際のリードと同等の物性を持つ材質で形状のみを簡素な立体（本形態では直方体）に単純化した構造体である。したがって、リードサンプルに関する特性値（ｂreed_sample，ｌreed_sample，ｄreed_sample）は固定値である。なお、リードＭRやリードサンプルに関する各物性値や寸法を利用者が個別に入力する構成のほか、以上の各パラメータ（ρreed，Ｅreed，tanδreed，ｂreed_sample，ｌreed_sample，ｄreed_sample）の数値を複数種のリードＭRの各々について事前に記憶装置４２に格納しておき、利用者が入力装置４４で選択した種類のリードＭRに関する各パラメータを特性パラメータ変換部２１が記憶装置４２から取得する構成も好適である。

管楽器のマウスピースＭPにおける特性インピーダンスＺcは以下の式(a1)で表現される。なお、下式のＳinは、マウスピースＭPのうち管体とみなせる部分の始点における面積である。
Ｚc＝（ρair・ｃair）／Ｓin
＝（ρair・ｃair）／｛π・（φin／２）²｝ ……(a1)
図１２に示すように、特性パラメータ変換部２１は、音速ｃairと密度ρairと直径φinとについて式(a1)の演算を実行することで特性インピーダンスＺcを算定する。なお、φinは、リードＭRの根元（マウスピースＭPに固定された部分）におけるマウスピースＭPの内径［m］である。例えば管体モデルにおける第１番目の単位部Ｕ[1]の内径φ1が式(a1)の直径φinとして適用される。

また、唇ＭLのバネ定数の分布ｋlip(x)［N/m²］は以下の式(a2)で表現される。
ｋlip(x)＝{（Ｅlip・ｂlip(x)・ｌlip(x)）／ｄlip(x)}／ｌlip(x)
＝Ｅlip・ｂlip(x)／ｄlip(x) ……(a2)
特性パラメータ変換部２１は、図１２に示すように、唇ＭLの物性値や寸法（Ｅlip，ｂlip(x)，ｄlip(x)）について式(a2)の演算を実行することで唇ＭLのバネ定数の分布ｋlip(x)［N/m²］を算定する。式(a2)において、Ｘ方向の位置ｘにおける唇ＭLの横幅ｂlip(x)および厚さｄlip(x)は音高ｆnから特定される（詳細は後述する）。

唇ＭLの内部抵抗の分布μlip(x)は以下の式(a3)で表現される。なお、式(a3)におけるｍlip_sampleは唇サンプルの質量［kg］であり、ｌlip_sampleは唇サンプルのＸ方向における長さ［m］であり、ｋlip_sampleは唇サンプルのバネ定数［N/m］である。

特性パラメータ変換部２１は、図１２に示すように、唇ＭLの物性値（ρlip，Ｅlip，tanδlip）と唇サンプルの寸法（ｂlip_sample）とについて式(a3)の演算を実行することで唇ＭLの内部抵抗の分布μlip(x)を算定する。なお、本形態においては、単純な直方体の唇サンプルを対象とした式(a3)の算定値で内部抵抗の分布μlip(x)を代表しているので、内部抵抗の分布μlip(x)は位置ｘに依存しない固定値となる。

一方、リードＭRの内部抵抗の分布μreed(x)は以下の式(a4)で表現される。ｍreed_sampleはリードサンプルの質量［kg］であり、Ｉreed_sampleはリードサンプルの断面二次モーメント［m⁴］であり、ｋreed_sampleはリードサンプルのバネ定数［N/m］である。

特性パラメータ変換部２１は、図１２に示すように、リードＭRの物性値（ρreed，Ｅreed，tanδreed）とリードサンプルの寸法（ｂreed_sample，ｄreed_sample，ｌreed_sample）とについて式(a4)の演算を実行することでリードＭRの内部抵抗の分布μreed(x)を算定する。なお、本形態においては、単純な直方体のリードサンプルを対象とした式(a4)の算定値で内部抵抗の分布μreed(x)を代表しているので、内部抵抗の分布μreed(x)は位置ｘに依存しない固定値となる。

また、図１２に示すように、特性パラメータ変換部２１は、アンブシュア（吹奏時の唇ＭLの状態）に関する複数のパラメータ（ｂlip(x)，ｄlip(x)，ｘteeth1，ｘteeth2，ｘlip1，ｘlip2，Ｆlip(x)）と、息圧Ｐ0を補正するための係数ｐmulと、管楽器の運指に関する複数のパラメータ（ｒH1〜ｒHk，ｒB，ＭH1〜ＭHk，ＭB）とをキースケール処理（図１２の記号“KSC”）で音高ｆnに基づいて特定する。キースケール処理は、音高ｆnがとり得る各数値とパラメータの各数値とが対応づけられたテーブルから、実際に指示された音高ｆnに対応する数値を各パラメータについて特定する処理である。

アンブシュアに関する複数のパラメータは、唇ＭLの横幅（Ｚ方向の寸法）ｂlip(x)［m］と、外力が作用しないときの唇ＭLの厚さ（Ｙ方向の寸法）ｄlip(x)［m］と、演奏者の歯ＭTが唇ＭLを押圧する力Ｆlip(x)［N］と、リードＭRに対する演奏者の唇ＭLや歯ＭTの位置に関するパラメータ（ｘlip1，ｘlip2，ｘteeth1，ｘteeth2）とを含む。

特性パラメータ変換部２１は、音高ｆnに対応する横幅ｂlip(x)および厚さｄlip(x)をキースケール処理で特定するとともに、横幅ｂlip(x)と厚さｄlip(x)との乗算値に唇ＭLの密度ρlipを乗算することで唇ＭLの質量の分布ｍlip(x)［kg/m］を算定する。また、横幅ｂlip(x)および厚さｄlip(x)は、前述した唇ＭLのバネ定数ｋlip(x)の算定にも適用される。

図５に示したようにＸ方向の各位置ｘを離散化するために、特性パラメータ変換部２１は、唇ＭLの位置（ｘlip1，ｘlip2）を間隔Δｘで除算した数値を離散後の位置（ｎlip1，ｎlip2）として算定し、歯ＭTの位置（ｘteeth1，ｘteeth2）を間隔Δｘで除算した数値を離散後の位置（ｎteeth1，ｎteeth2）として算定する。さらに、特性パラメータ変換部２１は、位置ｘteeth1と位置ｘteeth2との差分値をＸ方向における歯ＭTの長さｌteethとして算定し、位置ｘlip1と位置ｘlip2との差分値をＸ方向における唇ＭLの長さｌlipとして算定する。そして、特性パラメータ変換部２１は、押圧力Ｆlip(x)を歯ＭTの長さｌteethで除算することで、歯ＭTから唇ＭLの単位長あたりに作用する押圧力ｆlip(x)［N/m］を算定する（ｆlip(x)＝Ｆlip(x)／ｌteeth）。

特性パラメータ変換部２１は、音高ｆnに対応する係数ｐmulをキースケール処理で特定するとともに息圧Ｐ0と係数ｐmulとを乗算することで演奏者の口内の圧力Ｐを算定する。係数ｐmulは、音高ｆnに応じて変化する係数である。実際の管楽器においては管楽器を発音させるための演奏者の息圧の範囲が楽音の音高に応じて相違する（例えば、高音の演奏時の息圧のほうが低音の演奏時と比較して息圧の範囲が広い）という傾向がある。本形態においては、息圧Ｐ0に乗算される係数ｐmulが音高ｆnに応じた可変値であるから、息圧Ｐ0を音高ｆnとは無関係に選定した場合であっても、管楽器の以上の特性が忠実に模擬されるという利点がある。

また、特性パラメータ変換部２１は、管体模擬部３３のトーンホール部ＴH[1]〜ＴH[k]やベル部ＢLにて使用される係数ｒH1〜ｒHkおよび係数ｒBと、伝達模擬部３５にて使用される係数ＭH1〜ＭHkおよび係数ＭBとをキースケール処理で特定する。例えば、係数ＭHiは、音高ｆnの演奏時に第ｉ番目のトーンホールが閉塞される場合にはゼロに設定されるとともに第ｉ番目のトーンホールが開放される場合にはゼロを上回る所定値（例えば１）に設定される。同様に、係数ｒHiは、音高ｆnの演奏時に第ｉ番目のトーンホールが閉塞および開放の何れの状態に操作されるかに応じて別個の数値に設定される。

次に、図１３は、形状パラメータ変換部２３の具体的な機能を示すブロック図である。図１３に示すように、形状パラメータ変換部２３にはリードＭRや管体部の形状や寸法に関する様々なパラメータが供給される。形状パラメータ変換部２３に指示されるパラメータは、管体部を区分した各単位部Ｕ[i]の形状のパラメータ（Ｌi，φi，ｔi，ψi）と、リードＭRの厚みｙd(x,z)およびＺ方向における左端部および右端部の位置（ｚleft(x)，ｚright(x)）と、断面二次モーメントＩ(x)の基準となる軸線のＹ方向における位置ｙc(x)とを含む。

第ｉ番目の単位部Ｕ[i]の形状については、図６に示すように、単位部Ｕ[i]の長さＬiおよび内径φiとトーンホールの深さｔiおよび内径ψiとが指定される。形状パラメータ変換部２３は、第１に、接続部Ｊ[i]に関する係数（２ポート型については係数αiおよびβi、３ポート型については係数αi，βiおよびγi）を以上の各係数から特定する。第２に、形状パラメータ変換部２３は、各単位部Ｕ[i]の長さＬiから遅延素子ＤA[i]の遅延量ｄA[i]および遅延素子ＤB[i]の遅延量ｄB[i]を特定する。なお、以上のパラメータに加えて、ベル部ＢLにおける遮断周波数ｆCBとトーンホール部ＴH[i]における遮断周波数ｆCTHや遅延量（ｄE1，ｄE2）とを形状パラメータ変換部２３が可変に設定してもよい。

第３に、形状パラメータ変換部２３は、以下の式(b1)にリードＭRの各端部の位置（ｚleft(x)，ｚleft(x)）を代入することでリードＭRの横幅ｂreed(x)を算定する。
ｂreed(x)＝ｚright(x)−ｚleft(x) ……(b1)
第４に、形状パラメータ変換部２３は、以下の式(b2)に示すように、リードＭRの左端部の位置ｚleft(x)から右端部の位置ｚright(x)までの区間にわたってリードＭRの厚さｙd(x,z)を積分することで位置ｘにおけるリードＭRの断面積Ａ(x)を算定する。

第５に、形状パラメータ変換部２３は、位置ｙc(x)の軸線に関する断面二次モーメントＩ(x)を以下の式(b3)の演算によって算定する。式(b3)におけるdＡは面積分を意味する。

以上に説明したように、本形態においては、リードＭRと唇ＭLとの連成振動を表現する運動方程式Ｂに基づいてリードＭRの変位ｙ(x,t)が算定されるから、全体が自由に移動する剛体の空気弁としてリードＭRをモデル化する非特許文献１や長板状の振動体でリードＭRをモデル化する非特許文献２の技術と比較してリードＭRの挙動が忠実に模擬される。したがって、実際の管楽器の楽音に近い特性の楽音を高い精度で合成することが可能である。しかも、運動方程式Ｂにおける唇ＭLの変位ｙb(x)は、唇ＭLからリードＭRに作用する押圧力ｆlip(x)が変化するたびに、変化後の押圧力ｆlip(x)から運動方程式Ａ1および運動方程式Ａ2に基づいて算定された結果に更新されるから、押圧力ｆlip(x)を変化させる奏法を忠実に再現することができる。一方、押圧力ｆlip(x)が変化した場合であっても運動方程式ＢにおけるリードＭRの変位ｙ(x,t)は維持されるから、変位ｙ(x,t)の不連続な変化に起因した楽音の違和感は有効に抑制される。

＜Ｂ：第２実施形態＞
次に、本発明の第２実施形態について説明する。第１実施形態においては唇ＭLのバネ定数ｋlip(x)が歯ＭTからの押圧力ｆlip(x)に依存しない構成を例示したが、本形態においては押圧力ｆlip(x)に依存するバネ定数ｋlip(x,flip(x))を利用する。なお、以下の各形態において作用や機能が第１実施形態と同等である要素については、以上と同じ符号を付して各々の詳細な説明を適宜に省略する。

唇ＭLのバネ定数ｋlip(x,flip(x))と押圧力ｆlip(x)との関係は実測によって特定される。図１４は、バネ定数ｋlip(x,flip(x))の実測の方法を説明するための概念図である。図１４に示すように、作業台８０に載置された試験片８２の表面を加圧体８４で押圧する。試験片８２は、唇ＭLと同等の弾性特性をもつ弾性体である。演奏者の歯ＭTによる唇ＭLの押圧と同様に、加圧体８４は、試験片８２の表面の一部のみを押圧する。押圧力ｆlip(x)の強度を変化させながら試験片８２の変形量を測定してバネ定数ｋlip(x,flip(x))を算定する作業を、加圧体８４による押圧の位置ｘを変化させた複数の場合について反復する。以上の試験によって、押圧力ｆlip(x)とバネ定数ｋlip(x,flip(x))との関係が位置ｘ毎に実測される。

図１５は、試験片８２の特定の位置ｘを加圧体８４で押圧した場合に観察された押圧力ｆlip(x)とバネ定数ｋlip(x,flip(x))との関係を示すグラフである。図１５に示すように、試験片８２のバネ定数ｋlip(x,flip(x))は押圧力ｆlip(x)の強度に応じて変化する。すなわち、押圧力ｆlip(x)の強度が上昇するほどバネ定数ｋlip(x,flip(x)は増加する。

以上の測定が完了すると、押圧力ｆlip(x)とバネ定数ｋlip(x,flip(x))との関係を近似する関数（例えばスプライン関数）が、複数の位置ｘの各々について特定される。さらに、物性値や寸法が相違する複数の試験片８２について以上の作業が反復されることで、押圧力ｆlip(x)の作用する位置ｘと押圧力ｆlip(x)の強度とバネ定数ｋlip(x,flip(x))との関係を規定する関数（以下「弾性関数」という）が複数種の唇ＭLの各々について特定される。各弾性関数は楽音合成装置１００の記憶装置４２に格納される。

利用者は、入力装置４４を適宜に操作することで複数種の唇ＭLの何れかを選択する。図１の特性パラメータ変換部２１は、利用者が選択した唇ＭLに対応する弾性関数を記憶装置４２から取得し、キースケール処理で算定した押圧力ｆlip(x)を弾性関数に代入することでバネ定数ｋlip(x,flip(x))を算定する。特性パラメータ変換部２１が算定したバネ定数ｋlip(x,flip(x))は、リード模擬部３１（第１演算部３１１および第２演算部３１２）による演算に使用される。

以上に説明したように、本形態においては、押圧力ｆlip(x)が作用する位置ｘに加えて押圧力ｆlip(x)の強度によってもバネ定数ｋlip(x,flip(x))が変化する。すなわち、演奏時に歯から唇に作用する押圧力（ｆlip(x)）の強度や唇に対する歯の位置（ｘ）に応じて楽音が変化するという実際の管楽器と同様の作用が忠実に再現される。したがって、様々な奏法に対応した多様な楽音を忠実に合成できるという利点がある。

なお、以上においては試験片８２に対して部分的に押圧力ｆlip(x)を作用させたが、試験片８２の上面の全域に対して均等に押圧力ｆlip(x)を作用させてバネ定数ｋlip(x,flip(x))を測定する方法も採用される。以上の方法を採用した場合には、押圧力ｆlip(x)に応じて可変であるが位置ｘに依存しないバネ定数ｋlip(x,flip(x))が弾性関数によって規定される。したがって、歯から唇に作用する押圧力に応じて楽音が変化するという作用を再現することが可能である。

＜Ｃ：第３実施形態＞
唇サンプルやリードサンプルを使用した第１実施形態においては、唇ＭLの内部抵抗μlip(x)やリードの内部抵抗μreed(x)は位置ｘに依存しない固定値となる。本発明の第３実施形態においては、内部抵抗μlip(x)や内部抵抗μreed(x)を位置ｘに応じて変化させる。

唇ＭLの内部抵抗μlip(x)を定義する式(a3)のうち唇サンプルの横幅ｂlip_sampleを、位置ｘに応じた横幅ｂlip(x)に置換すると、以下の式(a3-1)が導出される。

同様に、リードＭRの内部抵抗μreed(x)については、位置ｘに応じて変化するリードＭRの断面積Ａ(x)およびバネ定数ｋreed(x)を変数とする以下の式(a4-1)が導出される。

図１６は、本形態に係る特性パラメータ変換部２１のブロック図である。図１６に示すように、特性パラメータ変換部２１は、唇ＭLの物性値や寸法（tanδlip，ｂlip(x)，ρlip，Ｅlip）について式(a3-1)の演算を実行することで位置ｘに応じた内部抵抗μlip(x)を算定する。式(a3-1)の横幅ｂlip(x)は、第１実施形態と同様に、音高ｆnからキースケール処理で算定される。

また、特性パラメータ変換部２１は、図１６に示すように、リードＭRの物性値（tanδreed，ρreed，Ａ(x)，ｋreed(x)）について式(a4-1)の演算を実行することで、位置ｘに応じた内部抵抗μreed(x)を算定する。形状パラメータ変換部２３が式(b2)の演算で算定した断面積Ａ(x)が式(a4-1)の演算に使用される。式(a4-1)におけるリードＭRのバネ定数ｋreed(x)［N/m］としては、例えば記憶装置４２に格納された数値や入力装置４４から指示された数値が使用される。

以上の手順で算定された内部抵抗μlip(x)および内部抵抗μreed(x)が、第２演算部３１２による運動方程式Ｂの演算で使用される。本形態においては、唇ＭLの内部抵抗μlip(x)およびリードＭRの内部抵抗μreed(x)が位置ｘに応じて変化するから、両者が固定値に設定された構成（例えば第１実施形態）と比較して、実際の管楽器の楽音を忠実に再現することが可能である。

＜Ｄ：第４実施形態＞
唇ＭLやリードＭRの変形が小さい場合（すなわち弾性限界内で変形する場合）、内部抵抗μlip(x)および内部抵抗μreed(x)を位置ｘのみに依存させた第３実施形態でも管楽器の楽音を忠実に再現することが可能である。しかし、唇ＭLやリードＭRの変形が大きい場合（すなわち変形が弾性限界外に到達した場合）、唇ＭLの内部抵抗μlip(x,flip(x))は、位置ｘに加えて唇ＭLに対する押圧力ｆlip(x)にも依存し、リードＭRの内部抵抗μreed(x,freed(x))は、位置ｘに加えてリードＭRに対する押圧力ｆreed(x)にも依存する。

図１７は、リードＭRに作用する押圧力ｆreed(x)とリードＭRの変位（変形量）との関係を示すグラフである。図１７に示すように、押圧力ｆreed(x)が所定値ｆTHを上回ると（すなわち弾性限界に到達すると）、リードＭRの変位は、押圧力ｆreed(x)の強度に応じて非線形に変化するようになる。すなわち、押圧力ｆreed(x)の強度が上昇するほどバネ定数ｋlip(x,flip(x))は減少する（変形し易くなる）。なお、唇ＭLからリードＭRに作用する押圧力ｆreed(x)は、リードＭRから唇ＭLに作用する押圧力ｆlip(x)と同等であるから、以下の説明では押圧力ｆreed(x)を便宜的に押圧力ｆlip(x)として表記する。

唇ＭLの内部抵抗μlip(x,flip(x))は以下の式(a3-2)で定義される。式(a3-2)のバネ定数ｋlip(x,flip(x))は押圧力ｆlip(x)の関数であるから、内部抵抗μlip(x,flip(x))は、位置ｘと押圧力ｆlip(x)とに応じて変化する。同様に、リードＭRの内部抵抗μreed(x,flip(x))は、以下の式(a4-2)で定義されるように、位置ｘと押圧力ｆlip(x)（バネ定数ｋreed(x,flip(x))）とに応じて変化する。

図１８は、本形態に係る特性パラメータ変換部２１のブロック図である。図１８に示すように、特性パラメータ変換部２１は、２種類のテーブル（Ｔlip，Ｔreed）を保持する。テーブルＴlipは、押圧力ｆlip(x)と唇ＭLのバネ定数ｋlip(x,flip(x))とを対応させる。テーブルＴreedは、押圧力ｆlip(x)（＝ｆreed(x)）とリードＭRのバネ定数ｋreed(x,flip(x))とを対応させる。テーブルＴlipやテーブルＴreedの内容は、例えば実際の唇やリードに押圧力を作用させた実験の結果に応じて設定される。特性パラメータ変換部２１は、キースケール処理で算定した押圧力Ｆlip(x)を歯ＭTの長さｌteethで除算することで算定した単位長あたりの押圧力ｆlip(x)に対応するバネ定数ｋlip(x,flip(x))をテーブルＴlipから検索し、押圧力ｆlip(x)に対応するバネ定数ｋreed(x,flip(x))をテーブルＴreedから検索する。

特性パラメータ変換部２１は、テーブルＴlipから検索したバネ定数ｋlip(x,flip(x))と唇ＭLの物性値（ｍlip(x)，tanδlip）とについて式(a3-2)の演算を実行することで、位置ｘと押圧力ｆlip(x)とに応じた内部抵抗μlip(x,flip(x))を算定する。式(a3-2)における質量の分布ｍlip(x)は、第１実施形態と同様に、横幅ｂlip(x)と密度ρlipとの乗算値である。また、特性パラメータ変換部２１は、テーブルＴreedから検索したバネ定数ｋreed(x,flip(x))とリードＭRの物性値や寸法（tanδreed，ρreed，Ａ(x)）とについて式(a4-2)の演算を実行することで、位置ｘと押圧力ｆlip(x)とに応じた内部抵抗μreed(x,flip(x))を算定する。

以上の手順で算定された内部抵抗μlip(x,flip(x))および内部抵抗μreed(x,flip(x))が、第２演算部３１２による運動方程式Ｂの演算で使用される。本形態においては、内部抵抗μlip(x,flip(x))および内部抵抗μreed(x,flip(x))が位置ｘおよび押圧力ｆlip(x)の強度に応じて変化するから、各内部抵抗が固定値に設定された構成（第１実施形態）や位置ｘのみに依存する構成（例えば第３実施形態）と比較して、実際の管楽器の楽音を忠実に再現することが可能である。なお、以上においては唇ＭLやリードＭRの変形が弾性限界外である場合を想定したが、弾性限界内の変形しか想定しない場合にも、図１８の構成を採用することは可能である。

＜Ｅ：変形例＞
以上の形態には以下に例示するような様々な変形を加えることができる。なお、以下の例示から２以上の態様を任意に選択して組合わせてもよい。

（１）変形例１
以上の各形態においては、利用者が入力したパラメータを特性パラメータ変換部２１や形状パラメータ変換部２３が楽音の合成に必要なパラメータに変換する構成を例示したが、合成部１４による演算に使用される各種のパラメータを利用者が直接的に入力する構成も採用される。例えば、図１２においてはアンブシュアや運指に関するパラメータをキースケール処理で算定する構成を例示したが、アンブシュアや運指に関するパラメータを利用者が入力装置４４から直接的に演算処理装置１０に指示する構成も好適である。

（２）変形例２
以上の各形態においてはリードＭRのヤング率Ｅreedと断面二次モーメントＩ(x)との乗算値を曲げ剛性Stiff(x)として算定したが、実測の結果から曲げ剛性Stiff(x)を特定する構成も好適である。例えば、リードＭRを模した試験片の各位置ｘに押圧力を作用させたうえで試験片の変位を測定した結果から曲げ剛性Stiff(x)を算定し、位置ｘと曲げ剛性Stiff(x)との関係を近似する関数（以下「剛性関数」という）を作成する。物性値や寸法が相違する複数種のリードＭRについて剛性関数が以上の手順で順次に生成されて記憶装置４２に格納される。演算処理装置１０のリード模擬部３１（第１演算部３１１や第２演算部３１２）は、複数種のリードＭRのうち何れか（例えば利用者が選択したリードＭR）に対応する剛性関数を記憶装置４２から取得して演算に使用する。以上の構成によっても第１実施形態や第２実施形態と同様の効果が奏される。

（３）変形例３
第２演算部３１２の算定した変位ｙ(x,t)から楽音を合成する方法は任意である。例えば、トーンホールやベルの内外の境界における音波の損失の模擬を省略した構成も採用される。

本発明の実施形態に係る楽音合成装置のブロック図である。 リード模擬部が模擬する管楽器のリードの近傍を示す概念図である。 管楽器の吹奏時における唇とリードとの接触の模式図である。 リード模擬部のブロック図である。 Ｘ方向の位置の離散化について説明する概念図である。 管楽器の管体部の模式図である。 管体モデルのブロック図である。 管体モデルにおけるベル部のブロック図である。 管体モデルにおける接続部のブロック図である。 管体モデルにおけるトーンホール部のブロック図である。 伝達模擬部のブロック図である。 特性パラメータ変換部のブロック図である。 形状パラメータ変換部のブロック図である。 唇のバネ定数を実測する構成の模式図である。 唇（試験片）に作用する押圧力とバネ定数との関係を示すグラフである。 第３実施形態における特性パラメータ変換部のブロック図である。 リードに作用する押圧力と変形量との関係を示すグラフである。 第４実施形態における特性パラメータ変換部のブロック図である。

符号の説明

１００……楽音合成装置、１０……演算処理装置、１２……設定部、１４……合成部、２１……特性パラメータ変換部、２３……形状パラメータ変換部、３１……リード模擬部、３１１……第１演算部、３１２……第２演算部、３１３……第３演算部、３１４……第４演算部、３２……範囲制限部、３３……管体模擬部、３５……伝達模擬部、４２……記憶装置、４４……入力装置、４６……放音装置。

Claims (8)

1. 吹奏時に唇に接触するリードの振動に応じて発音する管楽器の楽音を合成する装置であって、
   外力が唇に作用した平衡時におけるリードの挙動を表わす第１運動方程式と前記平衡時における唇の挙動を表わす第２運動方程式とを解くことで前記平衡時における前記唇の変位および前記リードの変位を算定する第１演算手段と、
   前記第１演算手段による算定の結果を前記唇の変位および前記リードの変位の初期値として、前記唇と前記リードとが相互に連成的に振動する連成振動を表現する運動方程式を解くことで、前記リードの変位を算定する第２演算手段と、
   前記第２演算手段が算定した変位に基づいて楽音を合成する楽音合成手段と
   を具備する楽音合成装置。
2. 前記唇に作用する外力の強度が変化するたびに、
   前記第１演算手段は、変化後の外力の強度に対応した前記唇の変位を前記第１運動方程式および前記第２運動方程式に基づいて算定し、
   前記第２演算手段は、前記第１演算手段が算定した前記唇の変位を前記連成振動の運動方程式に代入することで前記リードの変位を算定する
   請求項１の楽音合成装置。
3. 前記第１運動方程式および前記第２運動方程式は、前記唇における位置と押圧力の強度とに応じて変化する前記唇のバネ定数を含む
   請求項１または請求項２の楽音合成装置。
4. 前記第１運動方程式は、前記リードの位置に応じて変化する曲げ剛性を含む
   請求項１から請求項３の何れかの楽音合成装置。
5. 前記第２演算手段は、前記リードの変位を所定の範囲内に制限する
   請求項１から請求項４の何れかの楽音合成装置。
6. 前記連成振動の運動方程式は、前記唇における位置に応じて変化する当該唇の内部抵抗と、前記リードにおける位置に応じて変化する当該リードの内部抵抗との少なくとも一方を含む
   請求項１から請求項５の何れかの楽音合成装置。
7. 前記連成振動の運動方程式は、前記唇における位置および当該唇に作用する押圧力の強度に応じて変化する当該唇の内部抵抗と、前記リードにおける位置および当該リードに作用する押圧力の強度に応じて変化する当該リードの内部抵抗との少なくとも一方を含む
   請求項１から請求項５の何れかの楽音合成装置。
8. 吹奏時に唇に接触するリードの振動に応じて発音する管楽器の楽音を合成するためのプログラムであって、
   外力が唇に作用した平衡時におけるリードの挙動を表わす第１運動方程式と前記平衡時における唇の挙動を表わす第２運動方程式とを解くことで前記平衡時における前記唇の変位および前記リードの変位を算定する第１演算処理と、
   前記第１演算処理による算定の結果を前記唇の変位および前記リードの変位の初期値として、前記唇と前記リードとが相互に連成的に振動する連成振動を表現する運動方程式を解くことで、前記リードの変位を算定する第２演算処理と、
   前記第２演算処理で算定した変位に基づいて楽音を合成する楽音合成処理と
   をコンピュータに実行させるプログラム。

Images