2009-03-02

Multitaper method の直観的な説明・その３

科学研究

前回の続き。

有限長のデータを使う限り、スペクトルの「漏れ」は避けられず、その漏れ方は窓関数の関数形によって決まるということを見た。
これを踏まえ multitaper method (MTM) では、漏れに対する許容値をあらかじめパラメータとして与えた上で、「その許容値の範囲内で漏れが最小になるような窓関数を構成する」ことでスペクトル推定の精度を上げるというアプローチをとる。

データ長 N の時系列データ $X_t \quad (0 \le t \le N-1)$ を考えよう(簡単のためサンプリング間隔は1とする)。これをそのまま離散フーリエ変換すると、周波数解像度は $(2 N)^{-1}$ となる。言い換えれば、周波数領域では $(2 N)^{-1}$ Hz ごとにデータ点が存在する。

スペクトル漏れに対する許容値 W を、この解像度を単位として指定することにする。つまり、漏れの許容値を $W=1/(2N)$ と指定したら、それは周波数領域において、あるデータ点で表わされるスペクトルが隣のデータ点に漏れることは許容する（がそれ以上の漏れは許さない）ということで、 $W=2/(2N)$ ならさらにその隣に漏れてもよいとする、ということ。

さてここで、仮に許容値 W を完璧に実現するような窓関数が存在するとしよう。そのような窓関数のスペクトル分布は、周波数 f が $- W \lt f \lt W$ の範囲にあるときのみパワーを持ち、それ以外では 0 になっているはずだ。有限長の窓関数でそのようなスペクトル分布にできるだけ近い分布を持つものを求めるために、窓関数を N 次元のベクトル $\mathbf{w} = (w_1, w_2, ..., w_N)$ で表わし、この窓関数の周波数範囲 $- W \lt f \lt W$ でのパワー
$\int_{-W}^{W} |U(f)|^2 df \quad\quad \bigl(\quad U(f) = \sum_{t=1}^{N} w_t \exp (-2 \pi i f t)\quad \bigr)$
を最大化することで $\mathbf{w}$ を決めることを考える。

窓関数の規格化を制約条件として Lagrange の未定乗数法を使うと、
$\frac{\partial}{\partial \mathbf{w}} [ \int_{-W}^{W} |U(f)|^2 df - \lambda (|\mathbf{w}|^2 - 1) ] =0$
の解が求める窓関数となる。これは手計算で簡単に
$\sum_{t'=1}^N \frac{\sin [2 \pi W (t-t') ]}{\pi (t-t')} w_{t'} = \lambda w_t$
となることが示せるが、ここで $A_{tt'}=\frac{\sin [2 \pi W (t-t') ]}{\pi (t-t')}$ として行列 A を定義すると、この式は $A\mathbf{w}=\lambda \mathbf{w}$ という行列 A に関する固有値問題の解が、求める窓関数であることを表していることがわかる。

長くなったので続きは次回。

2009-02-19

Ubuntu8.10 で vim7.2 (+ vimgdb patch) をビルド

コンピュータ

Clewn というプロジェクトで vim と gdb を統合した vimgdb というのを作っていると知り、さっそく導入してみた。

デバッグだけに使うなら Clewn というアプリをインストールすればいいだけなのだが、普段使う vim に gdb を統合しようとすると、vim のソースに vimgdb のパッチをあててコンパイルする必要がある。ちょっと引っかかったところがあったので、備忘録としてメモしておく。

Ubuntu8.10 に入ってる vim は version 7.1 だが、どうせだから最新版の version 7.2 を使おう、ということで本家から tarball を落としてきた。さっそく展開して、vimgdb のパッチを当てて*1、とりあえず configure してみた*2。が、"terminal library が見つからない、ncursesとか…" みたいなことを言われて止まってしまう。

仰せのとおりに libncurses5-dev を入れてから仕切り直し。今度は configure、make、make install とエラーなしで進んだが、なぜか gvim ができてない。XもGTKも開発用ライブラリは入れてあるはずなんだが…。configure のログファイルを調べると、どうもXのライブラリが見つからないせいらしい。入ってるはずなんだが、と思いつつさらに詳しく調べると、Intrinsic.h というヘッダが見つからないせいでコケてるようだ。確かにそんなファイルはどこにもない。

どのライブラリに入ってるんだ、とgoogle先生に聞いてみると、すぐ答えてくれた。libxt-dev に入ってるようだ。てことでこれも入れて、今度はつつがなくインストール終了。

vimgdb についてはまたいずれ。

*1:vim の tarball を展開したのと同じディレクトリで vimgdb の tarball も展開し、patch -d vim72 --backup -p0 < vimgdb/vim72.diff とすればOK。

*2:vimgdb を入れる場合は --enable-gdb のオプションを付ける。日本語を使いたい場合はさらに --enable-multibyte --enable-xim --enable-fontset --with-features=big なども付ける。

2009-02-15

Multitaper method の直観的な説明・その２

科学研究

前回の続き。

Multitaper method (MTM) がどうしてうまくいくのかを理解するために、まず、スペクトルの「漏れ」がどうして生じるのかということと、従来の窓関数はどのようにしてこの「漏れ」を目立たなくしているのかということを説明する。

短時間のデータ、たとえば長さ T 秒のデータをフーリエ変換するというのは、無限長のデータに長さ T 秒の箱型の窓関数を掛けてフーリエ変換することに等しい。任意の無限長のデータを $x_t$ 、長さ T の箱型窓関数を $w_t$ 、それぞれのフーリエ変換を $X_f$ 、 $W_f$ とすると、
$w_t \cdot x_t \quad \Longleftrightarrow \quad W_f * X_f$
の関係が成り立つ( $*$ は畳み込みを表し、 $\Leftrightarrow$ は両辺をフーリエ変換・逆フーリエ変換で行き来できるという意味)*1。つまり、データに窓関数をかけてフーリエ変換すると、その周波数成分は、元のデータのフーリエ変換に窓関数のフーリエ変換を畳み込んだものになる。

純粋な正弦波（周波数: $f_0$ )に箱型窓関数を掛けてフーリエ変換することを考えてみる。元のデータの周波数成分は $f_0$ にしかないが、窓関数と掛け合わされると、窓関数の周波数成分が畳み込まれることによって本来の周波数以外の成分が出てきてしまう。これがスペクトルの「漏れ」が生じるメカニズムの直観的な説明。

箱形の窓関数でどういうスペクトル漏れが生じるかは、wikipedia の窓関数の項目を見れば図が載っているのでよく分かる(「矩形窓」のところ)。ここにはその他のいろんな窓関数についても周波数成分が示されていて、スペクトルの漏れ具合が窓関数によっていろいろ異なっているのがよくわかる。が、結局のところ、どんな窓関数を使っても本来の周波数の周りにスペクトルが漏れてしまうこと自体は避けられない。つまり、有限長のデータを用いる限り、スペクトルはどうしても漏れてしまうということだ。

以上を踏まえた上で、MTMではどのようにしてスペクトルの「漏れ」と「ばらつき」に対処するのかについて、次回以降説明していく。

*1:普通は時間領域での畳み込みが周波数領域で積に置き換えられるという形であらわされるが、どっちが時間領域でどっちが周波数領域かというのは計算する人のほうの決めごとに過ぎないので、 $t$ と $f$ を入れ替えてもそのまま成り立つ。

2009-02-12

Multitaper method の直感的な説明・その１

科学研究

名前はよく聞いていたがどんなものかよく知らなかった multitaper method (MTM) について、いろいろ文献を読んでみて、とりあえず「気持ち」の理解ができたようなので忘れないうちにメモ。

まずどういうときにMTMを使うかというと、たとえば時系列データの場合だと

短いデータからパワースペクトルを推定したいとき
データが長くても、パワーの時間変動を知るために(time-frequency analysis)、短時間のセグメントからパワー推定をしたいとき

などに使われる。2.の場合には wavelet が広く使われているが、wavelet とMTMの優劣については今のところよくわからない。だが実用的には wavelet のほうが便利だと思うし、もともとMTMはピークの検出のために開発されたっぽいので、time-frequency analysis には向いてなさそう*1。ということで主に1.の用途を念頭において説明していくことにする。

短いデータに対するパワースペクトル推定で問題となるのは、主に以下の2点。

短いデータをむりやり無限の繰り返し時系列の一単位とみなしてフーリエ変換することに起因する、スペクトルの「漏れ」(spectral leakage)。
たまたまそのデータに乗ってしまったノイズに起因する、スペクトルの「ばらつき」。

1.の「漏れ」というのは、たとえ元データが純粋な正弦波であっても、有限長のデータしかないと、フーリエ変換してもデルタピークにならず、周りの周波数にもパワーがあるように見えてしまうということ。これをなんとかする(ごまかす)ために、普通は窓関数というのを元データに掛けてやってからフーリエ変換するのだが、この窓関数の作りかた＋使い方を工夫することによって、スペクトルの「漏れ」を最小限に抑え、かつ2.の「ばらつき」をも抑えてしまうことができる、というのがMTMのメリット（だと理解した）。

次回に続く。

*1:たとえば、MTMでは基本的にセグメントの幅を決めてしまうと周波数の解像度はどの周波数帯でも一定なので、wavelet のように時間･周波数解像度のトレードオフを周波数ごとに調節することができない。(そのへんは short-time fourier transform と同じ。)

2009-01-25

ダイソーで買えるポスターケース(バズーカ)

研究

今日ダイソーで買い物してたら、なんとポスターを運ぶ例のバズーカを売ってるのを見つけてしまった。

「図面ケース」という品名で、525円で販売している。
事務用品店とかだと、一番安くて1、600円ぐらいからなので、えらく安い。
安いだけあってストラップがちょっと頼りない感じだが、本体部分は意外としっかりしている。
もちろん伸縮可能で、A0サイズに伸ばすこともできる。

ポスター＋ポスターケースは発表が終わるととたんに邪魔になってしまうので、いちいち持って帰ってこないで処分してしまいたいことがよくある。
そういうときはボール紙の筒に入れて運んで、発表が終わると捨ててしまうのだが、海外に行くときなどはボール紙筒だと手がふさがってしまって不便なので、ポスターケースを使うことになる。
で、毎回「じゃまだなあ」と思いながら持って帰ってくるわけだが、これだけ安いやつなら使い捨てしてしまってもいいかもと思えてくる。
ただダイソーは品替わりが激しいから、新しいのを買おうと思ったときにまだ売ってるかどうかは疑問。
どう考えても定番商品にはなりそうもないんだが、ネットで調べてみたら釣り師のひとがサオ入れに使ってるなんてこともあるようで、意外と本来の用途以外にもいろいろ使い道はあるのかも。

2009-01-18

PubMed のRSSフィードを読むのには NewsFox が便利

科学研究論文

生命科学系の研究者なら誰でも知ってる PubMed。
論文データベースなんだが、キーワードを登録しておくと、それを含む新着論文をRSSフィードでお知らせしてくれるという便利なサービスを提供している。
Firefox 使いの僕はもともと Sage でこのフィードを読んでいたんだが、あるとき NewsFox を試してみたところこれがめっぽう便利なので、以降ずっと使い続けている。
何が Sage と違うかというと、Sage だと新着論文のリストがだら〜っと表示されるだけなんだが、NewsFox だと

まずそれぞれのフィードについてどの雑誌に何報の新着論文があるかが表示され、
次に雑誌名をクリックするとその雑誌の新着論文のタイトル・著者等の情報が表示され、
さらに論文のタイトルをクリックするとその論文のアブストラクトが表示される。

これの何が便利かというと、PubMed でキーワード検索すると、かなりの割合で実際に自分の興味のあるものとはほとんど関係ないような論文がヒットしてしまう*1ので、どれが重要そうでどれがそうでないのかをいちいちふるいわけしなくてはいけないのだけれど、Sage のようにヒットした論文をすべていっぺんに表示されてしまうと、そのふるいわけを論文ひとつづつを見てやってやらないとけないので結構面倒くさいのに対して、NewsFox のように段階的に表示されると、それぞれの段階で大雑把にふるいわけができるので非常に楽なのである。
上記のように NewsFox では新着がまず雑誌ごとに分けられているので、この時点でまずまったく自分の分野と関係ない雑誌や、あまりにマイナーな雑誌の論文は無視できる。
次に雑誌名をクリックすると、論文のタイトルと著者が表示されるので、ここでまた、あまり面白そうでないタイトルのものなら無視するとか、面白そうでないタイトルでも著者が大物だったり知っている人だったりしたらチェックするとかで絞り込める。
で、最後にタイトルをクリックしてアブストを表示し、それを読んで最終的にちゃんとチェックすべき論文かどうかが判断できる。
もちろん、ここまでの段階で無視された論文たちは未読の状態で残るので、重要な論文を見落としてないかは、後で時間のあるときにでもチェックしなおすことができる。

Sage や NewsFox 以外にもいくつかフィードリーダを試してみたけど、NewsFox のように PubMed のRSSフィードを雑誌ごとに分けて表示してくれるものは見当たらなかった。
NewsFox 自体はフィードリーダとしてはあまり出来のいいものではないようなので、他にいいものがあったら乗り換えたいんだが、いまだに見つからないので、PubMed のフィードだけ NewsFox で読んで、他のは Sage-too で読むという中途半端なことをしている。
誰か代わりになるようなものをご存知でしたら教えてください。

*1:キーワードを細かく設定すれば絞り込むことができるのだが、そうすると今度は重要な論文を見落としてしまう可能性が高くなる。