跳慮跋考

興味も思考も行先不明

縦続行列による声道の音響特性の計算

以前の記事に書いた Ishizaka & Flanagan のモデルでは声道を輪切りにして声帯と一緒くたに微分方程式を解いていた訳ですが、これだとかなり巨大な方程式を解く事になって計算コストが馬鹿にならない、もっと言えばリアルタイムで音声合成をやろうとした場合に現実的でないのです。そこで Sondhi & Shroeter が考案したのが、声道部分は時間領域で微分方程式を解くんじゃなく周波数領域で特性を計算して、Fourier 逆変換で時間領域に戻す方法です。

Webster のホーン方程式

まず一般に断面積の変化する管を音がどう伝わるかを考えます。音は空気の圧力変化ですから流体力学によって記述されます。
音波は声道内の三次元空間を伝達する訳ですが、波長が声道の直径に比して十分に大きいとすれば、平面波と看做して一次元に落とし込む事ができます。

管の中央軸に沿ってX座標を取り、2つの底面を x=0, L とします。また管の断面積、音圧と体積速度(流速×断面積)をそれぞれ A(x), p(x,t), u(x,t) と置きます。
運動方程式として圧力差により体積速度が引き起こされる式
 \frac{\partial p}{\partial x} = -\frac{\rho}{A} \frac{\partial u}{\partial t} …(1)
と、連続の式として体積速度差により密度変化が生まれる式
 \frac{\partial u}{\partial x} = -\frac{A}{\rho} \frac{\partial \delta}{\partial t}
が基礎方程式となります。ここで δ は音圧による空気密度の変化量です。
音の伝播は断熱過程なので、空気を理想気体とすると気圧 p'(音圧は気圧の変化量ですから、これは平衡状態の大気圧に p を足したものです)と比熱比 γ について  p'/\rho^{\gamma} は一定、よって  \rho^{-\gamma} dp' - \gamma p' \rho^{-\gamma-1} d\rho = 0 です。すると  \partial p' / \partial\rho = \gamma p' / \rho が成り立つので、 \Delta p' = p = (\gamma p' / \rho) \Delta \rho = (\gamma p' / \rho) \delta となります。一方音速について  c^2 = \gamma p' / \rho ですから、連続の式は
 \frac{\partial u}{\partial x} = -\frac{A}{\rho c^2} \frac{\partial p}{\partial t} …(2)
と書けます。これらの式より u を消去すると
 \frac{\partial}{\partial x} A \frac{\partial p}{\partial x} = \frac{A}{c^2} \frac{\partial^2 p}{\partial t^2} …(3)
が得られます。この式を Webster のホーン方程式と呼びます。

P, U をそれぞれ p, u の Laplace 変換とすると、時間微分が s 倍に置き換わる事から (1), (2), (3) に対応して
 \frac{\partial P}{\partial x} = -\frac{\rho s}{A} U …(1')
 \frac{\partial U}{\partial x} = -\frac{As}{\rho c^2} P …(2')
 \frac{\partial}{\partial x} A \frac{\partial P}{\partial x} = \frac{s^2}{c^2} A P …(3')
の式が得られます。

損失のある管

管の壁による粘性摩擦を考慮すると運動方程式
 \frac{\partial p}{\partial x} = -\frac{\rho s}{A} U - R(x,s) U
の形で、また熱伝導と壁面の変形の効果を考慮すると連続の式は
 \rho \frac{dU}{dx} = -\frac{As}{c^{2}} P - Y(x,s) P
の形で上手く近似できるらしいです。この時、ホーン方程式は
 \frac{d}{dx} \frac{A}{\rho s + AR} \frac{dP}{dx} = \left( \frac{As}{\rho c^2} + Y \right) P …(4)
となります。ここまでは本題じゃないので手短に書いていますが、詳しくは Sondhi & Shroeter 自身による解説(pdf)など*1を参照して下さい。

縦続行列

さて、(4) は二階線形常微分方程式なので、その一般解は2つの独立な解(基本解)φ(x,s), ψ(x,s) の線形結合で書く事が出来ます。即ち
 P = a\phi(x,s) + b\psi(x,s)
です。この時 (1') 式より
 U = -\frac{A}{\rho s} \left( a \frac{\partial \phi}{\partial x} + b \frac{\partial \psi}{\partial x} \right)
が成り立ちます。
ここで x=0, L を代入し管の両端での音圧及び体積速度を考えると、それが a, b の線形結合で書ける事、即ち
…(5)
…(6)
が成り立っている事が分かります( P_{\mathrm{in}}(s) = P(0,s) という感じです)。すると、ここから a, b を削除すれば管の入出力の関係を表す式が導かれます。つまり解の正体は分からなくても、一点での値から別の一点での値が求められる訳です(分からないというか知る必要が無いという感じですが)。

とは言え基本解が分からないとと具体的な値を求められません。断面積 A が x によって変化する場合の解を見つけるのは困難ですが、これが定数であったとするとホーン方程式は
 \frac{d^2P}{dx^2} = \left( \frac{\rho s}{A} + R \right) \left( \frac{As}{\rho c^2} + Y \right) P
 = \frac{1}{c^2} \left( s + \frac{A}{\rho}R \right) \left( s + \frac{\rho c^2}{A} Y \right) P
となるので*2 \sigma^2 = \frac{1}{c^2} \left( s + \frac{A}{\rho}R \right) \left( s + \frac{\rho c^2}{A} Y \right) とすれば*3 cosh(σx), sinh(σx) が基本解として取れます。
consh の方を φ とすると、 (5), (6) より

即ち、入出力の関係式を
…(7)
と置くと

となります。K は縦続行列(chain matrix)や ABCD 行列*4と呼ばれる様です。

…(8)
が成り立っています。

声道モデル

前節により円筒に入出力する音圧及び体積速度の関係が分かったので、声道を輪切りにしたものそれぞれを円筒で近似する縦続音響管モデルにおいて (7) の関係が各円筒で成り立っている事、つまり
…(9)
が分かります。一方で

ですから、(9) を纏めると

が声道全体の縦続行列になります。
実際には K は s の、つまり周波数の函数なので、求めたい点毎に行列積の計算が必要になりますが、これにより声道の音響特性を調べる事が出来ます。また音響管の接続と看做せればよいので、スピーカーなどの特性を調べる場合にも適用可能でしょう。
例えば音圧の伝達函数は、

と置くと


より口での放射インピーダンス  P_{\mathrm{out}} / U_{\mathrm{out}}  Z_{\mathrm{out}} と置くと、(8) より  \det K = 1 なので

と求められます。

シミュレータ

触って動かせる声道シミュレータを作ってみたので、宜しければこちらをどうぞ。
JavaScript ライブラリの jQuery と Highcharts を使っていて、ChromeFirefoxOperaIE の最新版で動作する事を確認しています。
円筒を表す四角の上端をドラッグする事でその断面積を操作できます。
f:id:quinoh:20140214211313p:plain
f:id:quinoh:20140214211335p:plain
管 2 本だけの場合とより詳しいデータの場合を比較すると、案外に近い特性を与えている事など分かるかと思います。 2 管モデルの形状は Macquarie univ. のサイトを、複雑な方は以前と同じく MRI による計測データ*5参考にしました。

*1:他に日本語では「音響管共鳴の可視化と基本的な声道モデルの構成手順」がネット上で読めます。というかこれらを読めばこの記事は別に要りませんが

*2:R(x,s) は管の壁面による効果なので、断面積が x に依存しないならば R も依存しないと考えてよい

*3:右辺は複素数であるが、σはこの式を満たしさえすればよいので、根号の枝は自由に取れる

*4:各成分を ABCD で表した為だと思われる

*5:Story, Brad H.; Titze, Ingo R. "Vocal tract area functions from magnetic resonance imaging". J. Acoust. Soc. Am. 1996, Vol. 100, No. 1, p. 537-554. [ResearchGate]