2014-06-26

Firefox OS for Raspberry Pi

スマートフォンなる多機能携帯電話の囹圄も及ばぬ不自由さに嫌気が差した昨今、戯れに持ちたる Rapberry Pi をスマフォと PC の中間みたいなものに仕立て上げられないかと思い立っては種々様々に画策し、FirefoxOS を動かせるらしいと聞いては b2g-17.0a1.linuxgl-gnueabi-armhf_v6 なる謎のバイナリを実行するも不安定過ぎ（ブラウザのアドレスバーを触ると死ぬ）、仕様が無いので自らコンパイルするに至りました。

FirefoxOS をコンパイルする

そこら辺に転がってた HP Pavilion dv6（CPU: Pentium 2.00 GHz, Memory: 2 GB）に入れた Ubuntu 14.04LTS 上で Firefox OS for Raspberry Pi に従い

# python_curses は見付からなかったけど問題ないみたい
$ sudo apt-get install git mercurial diffstat chrpath
# ビルドで要求されるもの。texinfo -> autoinfo, libsdl1.2-dev -> sdl-config
$ sudo apt-get install g++ gawk texinfo autoconf2.13 libsdl1.2-dev
# Ubuntu のバグ (http://lostquery.com/questions/728/importerror-no-module-named-_sysconfigdata_nd)
$ cd /usr/lib/python2.7
$ sudo ln -s plat-x86_64-linux-gnu/_sysconfigdata_nd.py .

$ cd ~/src
$ git clone git://yoctoproject.org/poky
$ cd ~/src/poky
$ git clone git://git.yoctoproject.org/meta-raspberrypi; git clone git://git.openembedded.org/meta-openembedded; git clone git://github.com/imphil/meta-b2g.git
$ . ./oe-init-build-env rpi-build
# ここで conf/bblayers.conf と conf/local.conf を編集
$ cd ~/src/poky
$ . ./oe-init-build-env rpi-build
$ bitbake -v rpi-b2g-image

$ sudo umount /dev/sdl*
$ cd ~/src/poky/rpi-build/tmp/deploy/images/raspberrypi
$ sudo dd if=rpi-b2g-image-raspberrypi.rpi-sdimg of=/dev/sdl
$ sync

という感じでした。 4 並行でやりましたが、途中で "ImportError: No module named _sysconfigdata_nd" 等とエラーが出た（上の Ubuntu のバグ）為やり直したのを含めて 7.2 時間という所。

実機で動かす

さてこれを Raspbery Pi の Model B に挿して起動し、ネットに繋がらないのでモデムを再起動したりしつつ

$ B2G_HOMESCREEN=http://www.mozilla.org b2g &

で見事ブラウザが起動しました。 Redirect がどうのと意味不明の供述をしますが一応アクセス出来ているみたいです。しかし全くキー入力を受け付けず操作不能なので電源コードを抜いて強制終了しました。再起動し、何故か & を付けてバックグラウンドで実行してるのを辞めると

$ B2G_HOMESCREEN=http://www.google.com b2g

終了できないのは改善しましたが、依然 Tab、Alt（何故かこれでボタンを押す）、Backspace やマウス（見えないがクリックは出来る）しか入力を受け付けない様子でした。せめて文字を入力できないと検索とか全然できないんですが、下からにゅっとキーボードが出て来たりしないんですかね？

JavaScript のイベントを調べる

ちゃんとキーボードの入力を受け取っているのか不安になってきたので、JavaScript の keydown とかで

keyCode	key
8	Backspace
9	Tab
14	[
32	Alt
48	0
...
57	9
59	l
61	^
65	Control
66	v
67	x
68	s
69	w
70	d
71	f
72	g
73	u
74	h
75	j
76	k
77	n
78	b
79	i
80	o
81	p
82	e
83	a
84	r
85	y
86	c
88	z
89	t
90	]
113	F1
...
121	F9
173	-
188	m
191	.

名状し難き惨状。 vxswdfguhjknbiopearyczt という配列は一体何なのでしょうか。 Google 先生に訊くと更に謎な PASTEBIN 一件だけを返してきます。また表にない , @ p 等は反応が無く（Enter も！）、Printscreen については押すと b2g が終了します。謎仕様。

一方マウスでは mousemove/mousedown/mouseup がちゃんと発生して、mousemove は clientX/clientY の座標も返してくれるので操作性はある程度どうにかなりそうですが、mousedown 時の button 値は常に 0 が返る様です。

跋

まぁ総じて、動くっちゃ動くという感じなんですが、どうにかなりませんかね、という。誰かに真面に操作できる様にして下さい。

それから先般 3.2インチ液晶モジュール（TP付き） - aitendo@shopping を注文したので、続編があるかも知れません。

2014-06-03

LyX でアレがやりたい集

「TeX を捨てよ、LyX で書こう - 跳慮跋考」の補遺。順次更新、リクエスト募集。
モジュールの追加は文書(D)→設定(S)...→モジュールで行う。

定理環境とかの中に箇条書きを入れたい

「定理」モジュールを使って普通に
定理 1.
1.
2.
←ここへ続きを書こうとすると定理 2. ってなって困る。
これは途中の箇条書きが定理 1. の外に出ているのが原因。
箇条書き部分を選択して（1. が選べなくても気にしない）、編集(E)→リストの階層を下げる(I) をするといい感じに。

箇条書きの番号を変えたい

「調節可能な箇条書き(enumitem)」モジュールを追加する。
箇条書き(連番) を作成して右クリックすると「箇条書き(連番)のオプション」がある。
これを選ぶと 1. の右に箇条書き(連番)のオプション（←ここにボックス）が出来るので、ボックスのとこに label= と入力。
ここで Ctrl+L して TeX コードを書くボックスを作る。これをやらないと下の書式が解釈されない。
そこに \roman* 等と書けば、出力では番号の書式が変更される。
書式の指定ついては「enumerate 環境の箇条書きを、括弧付きにしたり英語にしたり - joker8phoenix's diary」で言う

\renewcommand{\labelenumi}{(\roman{enumi})}

とかの \roman{enumi} を \roman* に替える感じ。つまり * が適宜 {enumi} とか {enumii} とかに置換されるらしい。
セクション番号を入れたりしたい場合は「LaTeXしよう！ - 番号の変更」とかを参照されたし。

分からない・上手く出来ないっぽい事

「上手く出来ない」ってのはプリアンブルとかを弄る必要があるって事。

定理 1. とかの後ろで改行する

2014-05-31

TeX を捨てよ、LyX で書こう

「本物の物書きは TeX を使う」とは二十世紀末の著名な格言ですが（要出典）、実際問題これが中々に面倒な訳です。
TeX は随分と原始的なので、見た侭を編集できない。
少し書いてはコンパイル、どこの $ が閉じてないんだか分からない。
まぁ凡そ人間が直接書くべきではないのです。
ではどうするか。
その一つの解が GUI を備えた TeX/LaTeX 文書作成ソフト LyX です。「捨てよ」って言った割に内部で文書生成するのは TeX ですが気にしない。
使い方も分からない道具について WYSIWYM がどうのと思想を語っても面白くないので、ここでは兎に角も楽して文書を作る為の方法を解説します。

初期設定

インストールはまあいいとして、日本語で文書を作る準備です。
一先ず適当に新規（Ctrl+N）して、メニューバーから文書(D)→設定(S)...（以下では「文書の設定」と呼称）を開きます。
最初に〔文書クラス〕から文書クラス(C) を選択しましょう。適当に数式交じりの文書とかレポートとかを書くならば「日本語Article (jsarticle)」で良いと思います。
次に〔言語〕から言語(L) に「日本語」、文字コードでその他(E) に「日本語(pLaTeX) (SJIS)」*1、言語パッケージ(K) に「なし」を設定して、［文書の既定値として保存］したら［OK］。

基本操作

以上の設定をした上で何か打って Ctrl+R または文書(D)→[DVI]を表示(V) すると、ちゃんとプレビューが出来ると思います。
タイトル・著者とか章や節などの要素は、ツールバーの〔標準〕ってなっているプルダウンメニューから挿入できます。
それにファイル(F)→読み込み(I)、書き出し(E) の各種で入出力すれば、最低限文書作成の用途には使えるでしょう。

数式

数式は Ctrl+M で行内、Alt-M D で別行立てにて挿入されます。尚別行立て数式の後ろに文字を打つと自動的に改行されます。このとき Enter によって改行すると、そこから別の段落が始まると解釈されるので注意しましょう。
挿入(I)→数式(H) から align 環境とか gather 環境とかも勿論使えます。
また数式内では Alt-M N により数式番号の有り無しが切り替わります。
ツール(T)→設定(S)... の〔編集〕→捷径*2から他にも色々とショートカットがあるっぽい事が分かりますが、個人的にはこういう「見て分からない UI」って好きじゃないですねえ。
挿入(I)→数式(H)→マクロにより数式内で使うマクロが作成でき、これは入力するとちゃんと解釈してくれます。\R := \mathrm{R} とすると、それ以降で \R と打てば ℝ が表示されると。

定理環境

節毎に定義とか定理とか問とかの番号を振る奴です。理数系ならば必須と言えるでしょう。
「文書の設定」で〔モジュール〕を開き、選択可能のリストから「定理」*3を選んで［追加(D)］します。［OK］して〔標準〕のとこを開くと、下の方に定理とかが追加されています。
「定理(節毎連番)」とかも追加すると使えるんですが、もっと細かく*4設定できないものですかね。

図

初めに挿入(I)→フロート(A)→図でキャプションとかのスペースを作ります。それで「図1:」とかの前に移動して、挿入(I)→画像(G)... からファイルやオプションを設定し画像を入れます。eps とかいう謎形式じゃなくても読み込めるので優秀です。
編集(E)→段落設定(P)... から配置の中央揃え(E) を選ぶと図を中央に置けます。
因みに dviout で画像がモノクロなのは Option→Setup Parameters...→Graphic の GIF を BMP(full color) とかにすると直ります。

相互参照

式(1) とか図2 というのを本文に直接書くのは保守性とか文書構造を分かってない人間のやる事です。
例えば図や番号付きの数式内で挿入(I)→ラベル(L)... を選ぶと、数式に eq:Maxwellsou とか名前を付けられます。こうしておくと、文書中にて挿入(I)→相互参照(R)... から既存のラベルを選択して、出力ではそのラベルが付いたものの番号を表示する様に出来ます。数式の場合は形式(F) で括弧付きを選ぶとよいでしょう。

参考文献

巨人の肩に立つって奴ですね。
BibTeX の文献データベースファイル（JabRef とかで管理すると便利）があれば、それを導入して簡単に参考文献の一覧が作成できます。
先ず「文書の設定」から〔書誌情報〕の引用様式を NatBib にします。右の NatBib様式(S) で「連番」を選ぶと [1] とか、「著者-年」を選ぶと [Sondhi, 1974] とかが選べます。
挿入(I)→一覧/目次(I)→BibTeX書誌情報... から［追加(A)］→［一覧(B)...］で .bib ファイルを指定し［追加(A)］します。様式(Y) に「plainnat」を指定し［OK］を押下。これで本文中に挿入(I)→文献引用(C)... から引用を置く事が出来ます。

その他

他にも色々と機能がありますが、ここまでの操作で結構見掛ける事になったと思うので適当に弄ってみて下さい。
ツールバーのアイコンにも「数式を挿入」「ラベルを挿入」等があるので、人によってはそっちの方が便利かもしれません。
ヘルプの用途別説明書には xy-pic で可換図式を書く話もあるのですが、面倒だし結局邪悪なコマンドを憶える必要があるみたいですね。そんな事よりも TypeMath を使いましょう！

跋

LyX って日本語の解説が無い訳ではないんですが、散逸していて一々探さないとなんですよね。ヘルプの文書は充実してるんですけど、長いし。
ですのでライトな TeX ユーザーが気軽に乗り換えられる程度の必要最低限な解説があれば、より良い（この場合 platex と格闘する日本全国の人々が救われる様な）世界に貢献できるだろうと思い書いてみました。
アレが書いてないなんてあり得ない！　とかあったら twitter の方とかで。対応出来るか分かりませんが……。
TeX をゴリゴリ書く旧時代ワークスから一人でも多く解放される事を願って。

*1:ここで別のエンコーディングを選ぶと文字化けするのだが、どうしたものか。

*2:short cut ってまぁそういう事なんだけど何とも cool ですね。

*3:「定理」は theorem.sty で「定理(AMS)」は amsthm.sty だろうが違いはよく知らない。theorem.sty でも名前付き定理は使えるはずだが、「定理名付き定理」には「定理(AMS)」が必要と言っている。

*4:\theorembodyfont{\normalfont} で斜体になるのを防ぐとか、\theoremstyle くらいは定理環境に設定したいところですが、プリアンブルを弄っても上手くいかない気配。如何すべし。

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%以上を言うんだとか。