2014-02-14

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

音声合成

以前の記事に書いた 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 の線形結合で書ける事、即ち
$\begin{pmatrix} P_{\mathrm{in}} \\ U_{\mathrm{in}} \end{pmatrix} = \begin{pmatrix} \phi(0,s) & \psi(0,s) \\ -\frac{A}{\rho s} \phi'(0,s) & -\frac{A}{\rho s} \psi'(0,s) \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix}$ …(5)
$\begin{pmatrix} P_{\mathrm{out}} \\ U_{\mathrm{out}} \end{pmatrix}=\begin{pmatrix}\phi(L,s)&\psi(L,s) \\ -\frac{A}{\rho s}\phi'(L,s)& -\frac{A}{\rho s}\psi'(L,s) \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix}$ …(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) より
$\begin{pmatrix} P_{\mathrm{out}} \\ U_{\mathrm{out}} \end{pmatrix} = \begin{pmatrix} \cosh(\sigma L) & \sinh(\sigma L) \\ -\frac{A \sigma}{\rho s} \sinh(\sigma L) & -\frac{A \sigma}{\rho s} \cosh(\sigma L) \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & -\frac{A \sigma}{\rho s} \end{pmatrix}^{-1} \begin{pmatrix} P_{\mathrm{in}} \\ U_{\mathrm{in}} \end{pmatrix}$
即ち、入出力の関係式を
$\begin{pmatrix} P_{\mathrm{out}} \\ U_{\mathrm{out}} \end{pmatrix} = K \begin{pmatrix} P_{\mathrm{in}} \\ U_{\mathrm{in}} \end{pmatrix}$ …(7)
と置くと
$K = \begin{pmatrix} \cosh(\sigma L) & -\frac{\rho s}{A \sigma} \sinh(\sigma L) \\ -\frac{A \sigma}{\rho s} \sinh(\sigma L) & \cosh(\sigma L) \end{pmatrix}$
となります。K は縦続行列（chain matrix）や ABCD 行列*4と呼ばれる様です。
尚
$\det K = \{ \cosh(\sigma L) \}^2 - \{ \sinh(\sigma L) \}^2 = 1$ …(8)
が成り立っています。

声道モデル

前節により円筒に入出力する音圧及び体積速度の関係が分かったので、声道を輪切りにしたものそれぞれを円筒で近似する縦続音響管モデルにおいて (7) の関係が各円筒で成り立っている事、つまり
$\begin{pmatrix} P_{n+1,\mathrm{out}} \\ U_{n+1,\mathrm{out}} \end{pmatrix} = K_n \begin{pmatrix} P_{n,\mathrm{in}} \\ U_{n,\mathrm{in}} \end{pmatrix} \quad (n = 1, 2, \ldots, N-1)$ …(9)
が分かります。一方で
$\begin{pmatrix} P_{n+1,\mathrm{in}} \\ U_{n+1,\mathrm{in}} \end{pmatrix} = \begin{pmatrix} P_{n,\mathrm{out}} \\ U_{n,\mathrm{out}} \end{pmatrix}$
ですから、(9) を纏めると
$K_{\mathrm{tract}} = K_{N-1} K_{N-2} \cdots K_1$
が声道全体の縦続行列になります。
実際には K は s の、つまり周波数の函数なので、求めたい点毎に行列積の計算が必要になりますが、これにより声道の音響特性を調べる事が出来ます。また音響管の接続と看做せればよいので、スピーカーなどの特性を調べる場合にも適用可能でしょう。
例えば音圧の伝達函数は、
$K = \begin{pmatrix} k_{11} & k_{12} \\ k_{21} & k_{22} \end{pmatrix}$
と置くと
$\begin{pmatrix} 1 \\ U_{\mathrm{out}} / P_{\mathrm{out}} \end{pmatrix} = \begin{pmatrix} k_{11} & k_{12} \\ k_{21} & k_{22} \end{pmatrix} \begin{pmatrix} P_{\mathrm{in}} / P_{\mathrm{out}} \\ U_{\mathrm{in}} / P_{\mathrm{out}} \end{pmatrix}$
$\begin{pmatrix} P_{\mathrm{in}} / P_{\mathrm{out}} \\ U_{\mathrm{in}} / P_{\mathrm{out}} \end{pmatrix} = \frac{1}{\det K} \begin{pmatrix} k_{22} & -k_{12} \\ -k_{21} & k_{11} \end{pmatrix} \begin{pmatrix} 1 \\ U_{\mathrm{out}} / P_{\mathrm{out}} \end{pmatrix}$
より口での放射インピーダンス $P_{\mathrm{out}} / U_{\mathrm{out}}$ を $Z_{\mathrm{out}}$ と置くと、(8) より $\det K = 1$ なので
$\frac{P_{\mathrm{out}}}{P_{\mathrm{in}}} = \frac{1}{ k_{22} - k_{12} / Z_{\mathrm{out}} }$
と求められます。

シミュレータ

触って動かせる声道シミュレータを作ってみたので、宜しければこちらをどうぞ。
JavaScript ライブラリの jQuery と Highcharts を使っていて、Chrome、Firefox、Opera、IE の最新版で動作する事を確認しています。
円筒を表す四角の上端をドラッグする事でその断面積を操作できます。
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]

2013-11-25

まどか☆マギカ新編感想、或いは暁美ほむらと私

アニメ

（この文章は『劇場版魔法少女まどか☆マギカ [新編] 叛逆の物語』の感想でありネタバレを含みます。あと色々なものに酔って書いてるのでかなり読みにくいかも知れません）

私は暁美ほむらが好きだ。
何故と言われても困るが、あの第十話を観てからというもの「その、か、かっこいいな、なんて」「あなたは、どんな願い事をして魔法少女になったの？」「昨日、助けてくれたこと……絶対忘れたりしないもん！」等々の科白を聞くにつけてはほむらの心中が察せられて胸の締め付けられる思いをしてしまうのだ。
そうして私は『叛逆の物語』を観た。唖然であった。「愛よ」――知ってる、よく知っている。だが一体何をしたというのか、何がしたかったというのか。正直、救済前の奇妙に噛み合わない会話からずっとほむらが何を想い願って行動したのか全くの意味不明であった。
前半では悠々と我々を歓迎し、剰え「これが見たかったのでしょう？」「楽しいでしょう？」と問いかけておいて「でもそうじゃない」「それで分かったつもりか」と叩き落とす。今思い返してみるとこれは作り手からの、或いは暁美ほむら自身からの観る者への叛逆であったのかもしれない。
それでも彼女があの、一見理想の塊が如き世界で鬱々としてる事だけは痛切に伝わってくる。何を想うか分からなくとも、いや分からないからこそその感じだけが心の底に蟠り続けるのだ。
だから何度も観た。絶望の淵から希望を手繰る彼女の声が聞こえるまで。やがてほむらの行動だけは何となく腑に落ちるようになった。何故自死しなければならなかったかと言えば、改変後のほむら唯一の拠り所だったまどかの祈りを踏み躙ったからだろう、何故救済されてはならなかったかと言えば、まどかの祈りの重過ぎる代償を知ってしまったからだろうと。鹿目まどかを人間に戻すというのは、そういう壁に阻まれ道を閉ざされ追い詰められたほむらが「まどかを救う」為に見出したたった一つの冴えたやり方だったのではないかと。
だがそうすると、己の言わば「神堕とし」をほむらは何故「欲望」と呼んだのか、何故そんな（負の方向であるにせよ）積極的な言葉で表現したのか、私には彼女の言動がどうもしっくり来なかった。「鹿目まどかのままでいればいい」という想いが理想の押しつけだったとして、それに暁美ほむらが気づいたとしても、「それでも、あなたが幸せな世界を望むから」と言い切ったほむらが自らの苦心の選択肢を「欲望」と切り捨てるとは思えないのだ。
独立不羈を旗印に掲げる私ではあるが今回ばかりは如何ともし難く、色々と人の感想を読み漁り意見を求めた。その中で最も重要な示唆を与えたのは救済前の会話についてのもので、ほむらは「もう一度、あなたに会いたい」は「人間の」まどかと会いたいという事だ、と述べていた。成程そう考えれば彼女の行動にちゃんと筋が通ってくるのではあるが、果たしてそれが暁美ほむらの心からの願いなのか。「もう誰にも頼らない」と心に決めてからずっと距離を置き続け、ワルプルギスの夜を越えれば見滝原を去ると仄めかし、宇宙改変後はいつしか終わりの時に再開する事だけを頼りとしてきたあの暁美ほむらが何故人間の鹿目まどかを追い求めるのか。前編・後編を観返してみても今一つその実感は湧いてこなかった。はてさて進退窮まった、演繹が駄目なら発見的手法に頼らねばならない。

さて私は人の心中を推し量るのが上手くない。所謂心の理論が未熟なのかも知れないが、兎に角そんな私でも呆けて生きてきた中で多少は覚えた法則がある。それは、女の子は一般に「貸し借り」の天秤の傾きにかなり敏感であるという事だ。男の子はよく漫画で「貸し」だ何だと言うが、それと同じくらいに女の子は「貸し」を作らないで必ず「お返し」をする、というのが私にとっては大変リアルな事柄なのである。
そんな事を思い起こしつつまたほむらの歩んだ道をなぞってみると、どうやら彼女の最初の祈りというのは「鹿目さんを助ける」というよりは寧ろ鹿目さんに助けらた分ちゃんとお返し出来るような、そんな本当の友達、「対等な友達」になりたいという想いがあったのだろうと思えてくる。読み返してみればノベライズでのまどかは全く同じような悩みをずっと抱えているし、その観点からすると第十話において描かれたループの三周目は、魔法少女としてまどかとほむらが共闘して強大な敵を辛くも撃退し共に力尽きるという、ある種理想的な終局を迎えている。それはもう、暁美ほむらからすれば「何もかも滅茶苦茶に」してもいい程に。だが鹿目まどかはそうではなかった。彼女は最期まで「魔法少女としてみんなを守る」というアイデンティティを失いたくなかったのだ。故に、皮肉な事に「対等な友達」だからこそ、まどかはほむらに「助けて」と願ってしまう。友殺しの重き枷までつけて。斯くしてまどかを手に掛けたお下げに眼鏡の暁美ほむらは心の奥底へ封じ込められ、ただ約束を果たし罪を雪がんとする冷淡な暁美ほむらが生まれた。
だが契約を阻めば阻むほど、まどかは魔法少女の残酷な真実を知ってゆく。そして最後には、自らの存在自体を犠牲にして暁美ほむらを、世界の全ての魔法少女を救いたいと祈るのだ。それは魔法少女みんなの祈りを尊んだ結果かもしれないけれど、暁美ほむらだけはそこにはいなくて。「今までのほむらちゃん」の無駄にしていないのは行動だけで、ほむらの願いは結局何も叶えられる事など無かった。
神に等しくして、対等からは程遠く、感じることさえ出来ない。事ここに至っては、ほむらの願いはどうあっても叶えられなくなってしまったが故、彼女はまどかの最後の祈りを胸に抱いて生きるしかなくなった。「最高の友達」とは裏腹に、結局何もしてあげられなかった、罪を償えなかったという後悔を残して。

こうして漸く新編へと戻ってくる。
前半において綿密に描かれる夢の世界、あれは最初の暁美ほむらが願った世界だ。まどかと共に、みんなと共に楽しく戦い楽しく勝つ。
しかし、長い永い時を魔法少女の残酷な運命と戦ってきたほむらにとって、その世界は揺らぐ水面の鏡に映るが如く不安な違和感に満ちる世界だった。鹿目まどかの祈りを受け継いだ事だけが唯一の拠り所であった彼女にとって、それを無下にする世界には安住など出来よう筈もなかったのだ。
だがしかし、改変後のほむらの意志は夢の世界を巡ると共に瓦解してゆく。鹿目まどかの本当の気持ち、そして自分の魔女化。唯一残った拠り所さえも自らの浅ましさ、弱さが故に踏み躙ってしまったと、それはまどかが何よりも大切な全てを捧げてまで願ったものなのだと、そう気づいた彼女の絶望は如何に深かっただろうか。それでも「諦めないで」と言われた時、一体何を思ったか。
草原の二つの椅子は、きっと二人の理想の関係なのだろう。お下げのほむらはただ一緒に、髪を解いたほむらは救いたいと。しかしその手は届くことなく、高く舞い上がるあなたを貶めるだけ……。
夢の裡に自分の一番初めの祈りを、まどかの言葉に結局果たせなかった約束を思い起こさせられた暁美ほむらは、折り重なった時の中で極限を超えて濃縮された想いのままに、まどかを人間へと貶めたのだ。それは「対等の友達になりたい」「残酷な運命から救いたい」という二つの願いを叶える究極の選択肢である筈だった。
だが叶えてみるとどうしたことか、結局ほむらに残ったのは「まどかの祈りを踏み躙った」という事実だけではないか。だからこそ彼女は絶望したのだろう。私には受け継げなかったとリボンを返して、こんな結末を呼んだ自らの最初の祈りを「欲望」と罵って。それでも、まどかの祈りの大切さを知って尚、暁美ほむらは「あなたの幸せな世界を望むから」と宣言する。全てを手に入れたからにはもう失う事を恐れるしかないというのに、捨て切れない想い。後戻りはしないと新たな罪だけを抱えて生きるその姿は余りに痛々しい。

斯くして暁美ほむら視点では大したバッドエンドを迎えた訳であるが、鹿目まどかに眼を向けてみると彼女は相変わらずほむらの祈りの大切さ、或いは欲望というものを知らないままなのだ。だからきっと、ほむらの「愛」に勝つ事は出来ない。現状を打ち砕く手立ては恐らく、ほむらがまどかを落としでもしない限り見えてこないだろう。ほむらが、の部分は議論百出かも知れないが、兎角私は割と本気だ。愛を知らずして愛には勝てない。故に叫ぼう、まどほむに栄光あれ！

2013-10-13

Reading "Topographic Independent Component Analysis"

信号処理

序

感覚信号の効率的な(即ちスパースな)符号化に関して「視覚野・聴覚野地図の同一適応アルゴリズムによる解釈」に興味を持ったので、論文で使用されているトポグラフィック独立成分分析(TICA)の論文「Topographic Independent Component Analysis」を読もうとしたものの正直よく解らなかった為、Hyvarinen 先生の著「詳解独立成分分析―信号解析の新しい世界」で勉強してからリベンジした次第。以下に論文の要諦を纏めてみる。

いい加減 blog というメディアに適合していない気がするが……。

モデル

共通の分散

通常の ICA では信号源たる確率変数 s_i が互いに独立であることを仮定するが、現実には全く独立である場合と言うのは少なくて ICA を行っても幾許かの従属性が「残存」してしまうので、それを使って独立成分の二次元構造を決定しようというのが TICA の趣旨である。
特にこの論文ではエネルギー、即ち s_i² の相関に着目し、その相関が強いほど近い位置にあるとする。これは近い位置の成分ほどより共起しやすい(が値自体一方から一方を予測しがたい)と考えている事となる。
例えば変数 σ により成分 s_i, s_j (i ≠ j) に共通な分散が与えられるとし、z_i, z_j を白色*1な確率変数として
$\begin{cases} s_i = z_i \sigma \\ s_j = z_j \sigma \end{cases}$
とおく。ここで σ や z_i, z_j は平均 0 で互いに独立とする。
この時、s_i と s_j の値は互いに無相関だがエネルギーはそうではない。実際、独立性と白色性の仮定から
$E\{s_is_j\} = E\{ z_i \sigma \cdot z_j \sigma \} = E\{ z_i z_j \} E\{ \sigma^2 \} = 0$
$E\{ s_i^2 s_j^2 \} - E\{ s_i^2 \} E\{ s_j^2 \} = E\{ z_i^2 z_j^2 \} E\{ \sigma^4 \} - E\{z_i^2\} E\{z_j^2\} E\{ \sigma^2 \}^2$
　 $= E\{ \sigma^4 \} - E\{ \sigma^2 \}^2$
となり、最後の式は σ² の分散なので非負である。

成分の空間構造

s_i の分散 σ_i² を確率変数とし、各 s_i はその分散が与えられた下で互いに独立とする。
σ_i² の従属性は近傍函数 h(i,j) により記述する。通常近傍函数はある種の距離について単調減少する函数と定義され、対称性 h(i,j) = h(j,i) や、任意の i について h(i,i) = const. が要求される。
例えば、i 番目の成分の座標が (x_i, x_j) であるとすると
$h(i,j) = \begin{cases} 1 & (|x_i - x_j| \le 1\ \textrm{or}\ |y_i - y_j| \le 1) \\ 0 & (\textrm{otherwise}) \end{cases}$
とすれば自身と上下左右及び斜めの9近傍が取れる。
これにより、非線形函数 φ を*2用いて
$\sigma_i = \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right)$
で分散を定義する。(正直この φ が何を表すのかよく分からないのだが……)
以上の結果として TICA のモデルは、
$s_i = \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) z_i$
と書くことができる。

尤度函数

導出

まずトポグラフィックな成分 s と「高次」の独立成分 u の結合密度分布
$p(\textbf{s}, \textbf{u}) = \prod_i p(s_i | u_i) p(u_i)$
　 $= \prod_i \frac{1}{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } p^{s|u}_i \left( \frac{s_i}{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } \right) \cdot p^u_i(u_i)$
から(p^s|u_i は p^z_i なので確率密度の変換をしている)、s の周辺密度分布が
$p_s(\textbf{s}) = \int \prod_i p^{s|u}_i\left( \frac{s_i}{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } \right) \frac{ p^u_i(u_i) }{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } d\textbf{u}$
で得られる。W = (w₁, ..., w_n)^T として、 ICA の基本モデル s = Wx より
$p_x\left( \textbf{x}(t) \right) = p_s\left( W \textbf{x}(t) \right) |\det W|$
となるから、尤度函数は
$L(W) = \prod_{t=1}^{T} p_x\left( \textbf{x}(t) \right)$
　 $= \prod_{t=1}^{T} \int \prod_i p^{s|u}_i\left( \frac{ \textbf{w}_i^T \textbf{x}(t) }{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } \right) \frac{ p^u_i(u_i) }{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right)} |\det W| d\textbf{u}$
となる。

近似

上の尤度函数をそのまま用いるのは計算量が多く賢明でないので、扱い易い近似を求める。
まず p^s|u_i や p^u_i は全ての i について等しいとし、p_s|u = p^s|u_i はガウス分布とする。これは条件付きでない s の分布が優ガウス的*3である事を意味する。また非線形函数 φ を
$\phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) = \left( \sum_{k=1}^{n} h(i,k) u_k \right)^{-1/2}$
で定める。*4
これらの近似より s の周辺密度分布は
$p_s(\textbf{s}) = \int \frac{1}{(2\pi)^{n/2}} \exp\left[ -\frac{1}{2} \sum_i s_i^2 \left( \sum_k h(i,k) u_k \right) \right] \prod_i p_u(u_i) \sqrt{ \sum_k h(i,k) u_k } d\textbf{u}$
となり、和の順序を交換すると
$p_s(\textbf{s}) = \int \frac{1}{(2\pi)^{n/2}} \exp\left[ -\frac{1}{2} \sum_k u_k \left( \sum_i h(i,k) s_i^2 \right) \right] \prod_i p_u(u_i) \sqrt{ \sum_k h(i,k) u_k } d\textbf{u}$
と書ける。
ここで更に
$\sqrt{ \sum_k h(i,k) u_k } \approx \sqrt{ h(i,i) u_i }$
なる近似によって(実際にはこれは下界になっている)
$\widetilde{p}_s(\textbf{s}) = \int \frac{1}{(2\pi)^{n/2}} \exp\left[ -\frac{1}{2} \sum_k u_k \left( \sum_i h(i,k) s_i^2 \right) \right] \prod_i p_u(u_i) \sqrt{ h(i,i) u_i } d\textbf{u}$
となる。h(i,i) は定数であったので h₀ とおき、この積分を各 u_k について分解すると
$\widetilde{p}_s(\textbf{s}) = \prod_{k=1}^{n} \int \frac{1}{\sqrt{2\pi}} \exp\left[ -\frac{1}{2} u_k \left( \sum_i h(i,k) s_i^2 \right) \right] p_u(u_k) \sqrt{ h_0 u_k } du_k$
だが、u_k についての対称性から
$G(y) = \log \int \frac{1}{\sqrt{2\pi}} \exp\left( -\frac{1}{2} uy \right) p_u(u) \sqrt{ h_0 u } \,du$
という函数を定義すれば
$\widetilde{p}_s(\textbf{s}) = \prod_{k=1}^{n} \exp\left[ G\left( \sum_i h(i, k) s_i^2 \right) \right]$
と簡潔な形になる。

すると、対数尤度函数は
$\log \widetilde{L}(W) = \sum_{t=1}^{T} \sum_{j=1}^{n} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x}(t))^2 \right) + T \log|\det W|$
により近似できる。*5
ICA に於いて独立成分の確率分布(の対数)にあたる函数 G は、概形さえ合っていれば正しく収束する事が知られており、著者は指数分布に関連した函数として
$G_1(y) = -\alpha_1 \sqrt{y + \epsilon} + \beta_1$
を挙げている。ε は数値的な安定性の為に加えられる定数である。係数は s_i の分散が 1 であるという条件を満たす様に決められる定数だが、以下で W を直交行列とする制約を課す為にこれらは対数尤度を線形変換するだけなので、α₁=1, β₁=0 などと適当に決めても尤度の最大点には影響しない。

学習則

さて尤度が得られたので、勾配法によりこれを最大化しよう。まず両辺を T で割り
$\frac{1}{T} \log \widetilde{L}(W) = E\left\{ \sum_{j=1}^{n} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) \right\} + \log|\det W|$
と書き直す。(右辺第1項は標本平均なので、厳密にはこう書くべきでない。)すると w_k に関する勾配は
$\frac{\partial}{\partial w_{kl}} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) = G'\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) \cdot 2 h(k, j) (\textbf{w}_k^T \textbf{x}) x_l$
となるから、ベクトル勾配に直すと
$\frac{\partial}{\partial \textbf{w}_k} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) = 2 G'\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) h(k, j) \left( \textbf{w}_k^T \textbf{x}(t) \right) \textbf{x}$
となる。
観測値は z = Vx により白色化されているとし、W は直交行列とする。この時 det W は定数なので、結局
$\frac{\partial}{\partial \textbf{w}_k} \frac{1}{T} \log \widetilde{L}(W) = 2 E\left\{ r_k \left( \textbf{w}_k^T \textbf{z} \right) \textbf{z} \right\}$
となる。ここで
$r_k = \sum_{j=1}^{n} h(k, j) g\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{z})^2 \right)$
であり、g は G の導関数である。また W は仮定を満たす為に
$W \leftarrow (WW^T)^{-1/2} W$
で毎回正規直交化されなければならない。

毎回標本平均を取るのは実用的でない場面が多いので、確率的勾配降下法に基づくオンライン学習則として
$\Delta \textbf{w}_k \propto r_k \left( \textbf{w}_k^T \textbf{z} \right) \textbf{z}$
を用いる。これは単に期待値計算を外しただけだが、学習率などに適当な条件を課せばちゃんと収束する事が証明されている。

跋

以上により取り敢えず学習則が導出できたので、一旦筆を擱くとする。やっと原理が解ってきたので愈々実装と行きたいところ。
ただ φ はどうも納得がいかない……。

*1:共分散行列が単位行列に等しい。

*2:不幸にも環境によっては \phi ではなく \varphi が表示されているかもしれないが大目に見て欲しい。

*3:正の尖度を持つ。ガウス分布(正規分布)よりも平均近くの密度が高い。

*4:この函数が単調「減少」である事はどうにも不思議なのだが、私の「u_i が大きいほど活性が大きい」という認識が間違っているのだろうか？

*5:尤度が近傍のエネルギーの和のみの函数でよいというのは Kohonen の不変特徴部分空間の「部分空間の射影のノルムが不変特徴を定める」という話と関連があるらしい。

2013-08-17

『ソ・ラ・ノ・ヲ・ト』第四話より、人は微妙に組成が違うレンズの音を聴き分けられるか

アニメ

レンズの固有振動数ってどんなもんかな？　ってお話。
データも根拠も怪しいものですが、参考程度に。

基本振動だけ分かればいいので、レンズは円筒が連なって円柱を成しているものと看做します。
すると円筒のずれに対する復元力を求めるには、ガラスの剛性率(横弾性率、剪断弾性率とも)が必要。
山善特殊硝子製作所にデータがあるけど所々抜けているので、ガラスが等方的と仮定してヤング率とポアソン比から剛性率を求める*1。
$G = \frac{ E }{ 2 (1 + \gamma) }$
G：剛性率、E：ヤング率、γ：ポアソン比。

すると、長さ l で断面積 A の物体が ⊿x だけずれた時の剪断力 F は
$F = \frac{ GA \, \mathit{\Delta} x }{ l }$
となり、円盤の変形 x を半径 r の関数とし、l = ⊿r として円筒に適応すると
$\rho A \mathit{\Delta}r \frac{d^2x}{dt^2} = \frac{ GA }{ \mathit{\Delta}r } \{ x(r + \mathit{\Delta}r) - x(r) \} - \frac{ GA }{ l } \{ x(r) - x(r - \mathit{\Delta}r) \}$
$\rho \frac{d^2x}{dt^2} = G \frac{ x(r + \mathit{\Delta}r) + x(r - \mathit{\Delta}r) - 2 x(r) }{ \mathit{\Delta}r^2 }$
ここで右辺は二階差分であり、⊿r→0 の極限において二回微分と一致する。
$\frac{ \rho }{G} \, \frac{d^2x}{dt^2} = \frac{d^2x}{dr^2}$
即ち一次元の波動方程式と同等であるから、波の速さが v = √(G/ρ) と求まり、固定端条件での基本振動を考えるとその波長は 4R なので、v/(4R) で固有振動数が求まる。筈。

で、最後にレンズの半径が問題となる。劇中の描画からは直接は推定し難そうなので、統計値に頼る。
平成24年度の学校保健統計調査によると 15 歳の女の子の平均身長は 157.2 cm*2、身長別ボディサイズ平均なるサイトのデータから線形内挿すると手長は 17.56 cmとなる。大体その半分がレンズの直径と考えると、半径は 4.4 cmと推定される。
まぁ周波数比を考えるんなら半径は打ち消されちゃうけど。

ネットで組成比が分かったガラスについてデータを示す*3。

名称	パイレックス7740	テンパックス	パイコール7913	石英ガラス
密度 (g/cm³)	2.23	2.20	2.18	2.20
ヤング率 (10⁶ kg/cm²)	0.61	0.63	0.68	0.74
ポアソン比	0.20	0.22	0.19	0.16
屈折率 (589.3nm)	1.474	1.472	1.458	1.459
剛性率 (10⁶ kg/cm²)	0.254	0.258	0.286	0.319
波速 (m/s)	3342	3391	3584	3769
レンズの固有振動数 (Hz)	760	771	815	857
SiO₂	81.2	81	96.4	100
B₂O₃	11.8	13	3.0	0
Al₂O₃	2.4	2	0.5	0
Na₂O/K₂O	4.5	4	-	0

さて、ここで成分の左の2つに注目すると、これらは成分が似通っており屈折率も近い一方で固有振動数には 11 Hz、比にして 25 セントだから半音の 1/4 もの差があることが分かる。
音の環境心理学の図 2.6 を見ると、この音域での弁別閾は 1～2 Hzであるから、これは十分に聴き分けることができると考えられる。

*1:実はパイレックスとかこの計算値と表の値が結構喰い違ったりするんだけど……。データの斉一性を重視して全部計算値を採用しておきます。

*2:日本人のデータだが

*3:それぞれの組成比は平岡特殊硝子製作株式會社、Corning(pdf)より。石英ガラスは純度99.99%以上を言うんだとか。

2013-07-07

声帯の数値シミュレーションによる音声生成

プログラミング音声合成

初音ミク嬢を筆頭として音声サンプルの編集による音声合成は随分栄華を誇っていますが、力学モデルのシミュレーションによる音声合成はどうも人気が無い様ですね。私が世情に疎いだけかも知れませんけれど。
実用程度*1ならそうした経験的帰納的な手法も使えるでしょうが、完璧を求めるならば、整然とした理論モデルによる演繹的アプローチが必要不可欠だと私は思っています。そしてそれは声以外でも。

モデル

で、何をしたかというと、石坂とFlanaganによる声帯の二質量モデル*2を実装してみました。簡単に言うと、声帯を発条で振動する二つの壁の組み合わせとしてモデル化したものです。
もっと発展した研究の論文*3を参考にしたので、元の論文よりはL_sを無視してたり声道をもっと細かくしたりしています。それと声道のパラメタの式が納得いかなかったので修正したところも。
あと電気的等価回路については、これ(pdf)を見るとちょっと分かるかもですね。
声道断面積関数については色々ネットでデータを探して、MRIによる計測データ*4があったのでこれを使いました。

実装

んで微分方程式をプログラムに解かせる訳です。最初は連立一階常微分方程式を表すラムダ式とパラメタによって解を求めていくクラスを作ったんですが、声道の各要素に対応するラムダ式を簡単に書けず、Boost.Preprocessor などという黒魔法にも手を出してみましたがどうもきもちよくないので結局 GSL のを使いました。そしたら全ての変数を纏めて扱う ℝⁿ→ℝⁿ 函数みたいな設計になってて、なんかもうあっはいって感じですね。

兎に角そうして実装したのがこのコード(Gist)です。パラメタまで色々ハードコードしてるのは気にしないで下さい。C++ よく分かんないので適当に const つけまくったりしてますが意味あるんでしょうかね。なんか速くなったりしませんか？
論文に"The time derivative of the mouth volume velocity (i.e., through the radiation load) is good approximation to the radiated sound pressure."(p.1248)とあるので、dydt の最後の要素の値である dU_O/dt が最終的な音声波形出力(放射音圧)と看做せます。
途中に書いてありますが、y の中身は { H1: 質量1の変位, V1: 質量1の速度, H2: 質量2の変位, V2: 質量2の変位, Ug: 声門内で平均した体積速度, U1, ..., Un: 声道の各要素での体積速度, P1, ..., Pn: 声道の各要素での圧力, Uo: 口での体積速度 } となっています。U1,V1,U2,V2,...にした方がいい気もしますが面倒な割に御利益が無い気もしますね。
本当はパラメタとか最初に全部出力してますが省略。

さて実装で何が困ったって、普通に計算させると発散するんですよね。変位が小さくなると体積速度が上がる、体積速度が上がると変位速度が上がる、という感じの循環が起きてるってことだと思うんですが。
でどうしたかって、まぁ現実に照らし合わせれば数μmの隙間を物凄い勢いで風が通る訳がないので、変位に下限値(コード中の Hmin)を設定してやりました。何ともびみょーですがまず動かなくては仕様がありません。
論文に全然書いてないからもしかしたら私の実装した式が間違ってるかも知れないんですが……。

出力

g++ ならば最低限 -std=gnu++0x -lgsl -lgslcblas -lm あたりのコンパイラオプションつければ大丈夫だと思います。あとは最適化とか。
声道に/a/を選んで(上のコードは/o/ですが)、引数(初期値)を 0.018 0 0.018 0 0 にして実行すると以下の波形を得ます。上から二つの声門面積、声門内の平均体積速度、放射音圧です。
面積が負の方に行き過ぎだったり(質量が剛体として"衝突"するモデルではないので負になるの自体はいいのですが)体積速度がギザギザし過ぎだったりする気もしますが大体は再現できていますね。
f:id:quinoh:20130707011626p:plain

それで波形が出来たのはいいのですが、テキストファイルだから .wav か何かにしないと聴けませんので、適当な変換プログラム(Gist)を書きました。本当に適当なので多少サンプリングレートが狂いますがまぁ聴いても分からないでしょう。最大値くらいは自動で計算すべきかも。
斯くして漸く出来たのがこれ(TwitSound)。結構「あ」っぽいと思いますがどうでしょう？　「あいうえお」全部繋げてみるとこんな感じ(TwitSound)になって「い」「え」とかはもっと不自然だったりするので、この辺りは声道面積の値次第なのかどうなのかという感じですね。スペクトルとかを観察すると「不自然さ」が定量できるのでしょうか。

跋

取り敢えずこれで音声を生成する為のフレームワーク的なものは出来た訳ですので、これをどう自然な声に近づけるかというところが次の問題になります。
筋肉の緊張度を入力パラメータとして声道断面積を制御する事で波形を自然な声に制限したり発声を学習したり逆に音声認識の為の符号化で使えないか……など色々あるのですが、なかなかどれも大変そう。まぁまず声がどんなものかを知らねばなりませんかね。
それから計算速度が遅いのも困りものです。どうにか1秒を1秒以下で計算して欲しいところですが、これより時間刻みを大きくしても発散したりしてどうしたものか。
それとここまで頑張ったところでアレですが、「本質を簡潔に」モデル化するという観点からすると声道をバラバラにしてシミュレートしているのはちょっと泥臭過ぎる嫌いもあります。そうした考えを貫くならば、声道を音響フィルタと看做す「音響フィルタ理論」などを使った方が良いのかも知れません。複雑なモデルほど実装も操作も難しいので。

*1:「程度」なんて言うと色々な人に殺されそうですね

*2:Ishizaka, K.; Flanagan, J. L. "Synthesis of Voiced Sounds From a Two-Mass Model of the Vocal Cords". Bell System Tech. J. 1972, Vol. 51, No. 6, p. 1233-1268. [pdf]

*3:古賀博之; 中川匡弘. "カオス音声生成モデル". 電子情報通信学会技術研究報告. NLP, 非線形問題. 1998, Vol. 98, No. 343, p. 25-32. [CiNii]

*4: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]

2013-07-05

Cross compiling a program with PortAudio for Windows by MinGW on Linux

プログラミング音声解析

最近は VirtualBox 上の ubuntu で開発することで Windows の呪縛から解き放たれていたんですが(まあ Windows がホストOSではありますけど)、どうもマイク入力が VirtualBox ではサポートされていない様子。USB マイクなら大丈夫とか聞くんですが。
開発環境は他のプログラムとも合わせて Linux にしておきたいので、プログラムを Windows 向けにクロスコンパイルしようと思い立ちました。
環境や使用した(された)バージョンは以下です：

Microsoft Windows 8 [Version 6.2.9200]
VirtualBox 4.2.12
ubuntu 12.04 LTS
Microsoft Visual Studio Professional 2012 [Version 11.0.50727.1]
PortAudio V19 (2011-11-21)
MingGW 4.2.1

PortAudio のコンパイル

公式のチュートリアル(http://portaudio.com/docs/v19-doxydocs/compile_windows.html)に従えば大体大丈夫でしょう。
preprocessor definitions の設定は以下のようにしました：

PA_USE_ASIO=1
PA_USE_DS=1
PA_USE_WMME=1
PA_USE_WASAP=1
PA_USE_WDMKS=1
PA_USE_SKELETON=1
PA_WDMKS_NO_KSGUID_LIB

最後のは「エラー 38 error LNK1104: ファイル 'ksguid.lib' を開くことができません。」なるエラーを回避する為です(http://music.columbia.edu/pipermail/portaudio/2011-August/012848.html)。
これで運が良ければ Win32/Release や、構成マネージャを弄れば x64/Release にライブラリファイルが生成されているでしょう。

テストコードのコンパイル

まず MinGW を入れます。

$ sudo apt-get install mingw32

そして上で生成された portaudio_x86.lib とかを portaudio/test なり何なりに持ってきて

$ i586-mingw32msvc-gcc patest_read_record.c -I../include portaudio_x86.lib

生成されたプログラムを Windows 上で実行すれば、ちゃんと録音したりしてくれる筈です。

2013-05-04

共感能力に欠く人間の弁明

戯言

〈この物語は事実を基にしたフィクションです。〉
的な注意書きを見かける事がある。『月光の夏』だとか『電車男』だとか。
というのは恐らく、ノンフィクションと銘打ちながら脚色があったりすると「騙された！」だのと騒ぐ人が居るからなのだろうけれど、では「フィクションに対する感動」と「ノンフィクションに対する感動」は果たして別物なのだろうか、と私は思う。

「物語に於ける死」とは「現実に於ける死」よりも軽いのだろうか？　両者は異質でありその境は侵されざるべきものなのであろうか？　私にとってその答は全き否である。
仮令物語の中であろうと人が傷付き悲しめばその分心は暗く沈む。時には本を読んで一週間くらい鬱々としたり。
寧ろ逆に、「現実」の情報の方こそ無感動に受け取ってしまう。地球の裏側での大災害。へぇそうか。知らない土地の知らない人々が苦しむ姿、それにどれ程の実感がある？　国内だろうと変わりはしない。報道の映像はまるで映画の一場面の様だ。それは「本当にあった事」なのか？
私にとって「知らない土地」なんていうのは物語の舞台と何ら変わりはない。この目で見、この手で触れる事だけが存在証明となる。だがそれは、作家の精神活動と何ら違わないのではないだろうか。彼等は彼等の中の世界を見、触れて物語と成すのだ。決して創造主としてではなく、一人の語り部として。そして物語が共有された時、その世界は「実在」せずして何だと謂うのか。
私にとっての「現実」とは、つまりその程度のものなのだ。この手が届く小さな世界。あなたはそこにいますか？

跳慮跋考

興味も思考も行先不明