しげ日記: Multi taper (3)

2006年02月17日

昨日の続き。

なんとなく Matlab で簡単にコードが書けそうな気がしたので、試しに書いてみました。以下は、N=100 point、NW=5 の Multitaper を作るプログラムです。

　　N=100;W=0.05;

　　for i=1:N
　　　for j=1:N
　　　　A=2*W*sinc(2*W*(i-j));
　　　end
　　end

　　[V,D]=eig(A);

　　figure(1)
　　for k=1:5
　　　subplot(5,1,k)
　　　plot(V(:,N-k+1))
　　end

　　figure(2)
　　plot(diag(D))

sinc 関数は matlab 関数で、sinc(t) = sin(pi*t)/(pi*t) です。t = 0 のとき sinc(0)=1 となることに注意。N が大きくなるとループ部が重たくなるので、

　　i=repmat(1:N,N,1);
　　A=2*W*sinc(2*W*(i-i'));

とした方がいいかもです。

それで、固有ベクトルを固有値の大きいものから順に５つならべたのが以下の図です。Mitra 論文の Fig4 と同じ感じですね。これを窓関数として使えばいいのですね。

あと、固有値の分布を表示したのがその次の図です。固有値が小さい順に並んでいます。特徴は、大きいものから 2NW 番目（ここでは 2NW=10）までの固有値は、λ　～ 1 となり、それ以外は λ　～ 0 となる点です。

↑ アソハンがたまたま持ってた。Mitra 論文はこの本をベースにしてるみたい。

投稿者 sfujisawa : 2006年02月17日 20:32

サイン・インを確認しました、 . さん。コメントしてください。 (サイン・アウト)

(いままで、ここでコメントしたとがないときは、コメントを表示する前にこのウェブログのオーナーの承認が必要になることがあります。承認されるまではコメントは表示されません。そのときはしばらく待ってください。)