きのこな夜（研究室）

Texas InstrumentsのFilterProを利用すると簡単にできます。実際に回路をつくってトラブったときはQucsに回路図を書き写してシミュレートしてやれば、どこでおかしくなったのか一発で分かります。

これで、高次のフィルタがあっという間にできるぞ♪

ある振動する系からの信号のスペクトルが分かっていて、そのピークの位置から周波数を求めたいというときがあります。でも、普通は時系列は離散化されていているので、スペクトルの周波数の値も飛び飛びで（スペクトル幅があるので）、ある程度以上の精度で決定することは望めません。

ところが、SN比の良い信号であると、かなりの正確に周波数を求めることができます。たとえばLabviewの中にはBuneman Frequency estimatorというのが入っていて、これを使うと、ある種の補完を用いてスペクトルのピークの位置を推定してくれます。

ではBuneman法による周波数推定について考えて見ましょう。

続きを読む

ツールにあるやつを使って設計。適当に回路定数を調整して、実際に作れそうなものに。前後にバッファアンプをかませると簡単にできますね。すばらしい。明日作ろう。

Rを使って統計処理をするときに、たくさん似たような処理を繰り返すのでいちいちコマンドを書くのは面倒ということがあります。こういうときはFOR文を使うとかなりコードを短縮できます。
コツは変数やファイル名に連番をつけることです。あとは、

for(i in 番号) eval(parse(text=paste("
hoge_",i,"egoh
",sep="")))

のように使います。ここでiに連番の番号が入ります。

たとえば、ファイル名がA1.dat〜A100.datというファイルのデータをread.tableで読み込んでa1〜a100という変数に代入したいときは

for(i in 1:5) eval(parse(text=paste("
a",i," <- read.table('a",i,".dat')
",sep="")))

のように使います。「eval」と「parse」は、他にも便利な応用に使えます。たとえば、a1〜a100の2カラム目のデータの二乗平均を取りたいときは

m_sq <- eval(parse(text=paste("a",1:100,"＄V2 ^2",sep="",collapse="+")))/100)

とすればOKです。(＄は半角に直します。)
ヘルプファイルを見ていろいろ試してみてください。

参照
RプログラミングTips大全-(Rjpwiki)

Ubuntuサーバーと格闘(1) 〜インストールの続きです。

■がいっぱい出たのは、日本語のフォントを表示できるようになっていないからです。つーことでjbftermを使って表示してやりましょう。
以下の手順はほとんどココを参照しました。

まずやっとくこと
ディスプレイの説明書を引っ張り出してきて解像度と色数を調べといてください。
ちなみに拙者の環境では1024x768 16bitでやんした。

/boot/grub/menu.lstの編集
また、viさんを使って編集です。viでファイルを編集したいときは

sudo vi ファイル名

でやんす。

# kopt=root=/dev/hdd1 ro

とか言う行があるはずなのでこの部分を

# kopt=root=/dev/hdd1 ro vga=ほにゃらら

に変更します。
「ほにゃらら」に入る部分はディスプレイによりけり。ココにある表を参考にして置き換えてください。拙者の環境では1024x768 16bitなので「0x317 (791)」(xはアルファベットのエックス)で上の行は

# kopt=root=/dev/hdd1 ro vga=0x317

とすればOKです。

編集が終わったら

sudo update-grub
sudo reboot

最後の1行は再起動。

起動後
おそらく高解像度になっているので文字が小さいはずです。わーい。
必要そうなものを片っ端からインストールします。

sudo apt-get install xfonts-base (いらないかも)
sudo apt-get install jfbterm

■がいっぱい現れてやりにくいですが、なんとかインストールが終わったら

modprobe vga16fb
modprobe fbcon
jfbterm -q

と打ち込こみます。新しい画面が表示されて細いフォントの字が表示されるはずです。日本語が表示されるかテストしてみましょう。

man apt-get

（いままで使ってきたapt-getのマニュアル）

古いパソコンに、Ubuntuサーバー版をインストールしてみた。Web上にほとんど記録がなかったので作業記録を書いておきます。ちなみに私、LINUXに関しては超ド素人ですので、質問は受け付けられません。自己責任で・・・
ほとんどここを参照しています。
まずは・・・

ダウンロード

本家サイト⇒ダウンロード⇒Ubuntu 6.06 LTSのなかからJapan⇒Server install CDのなかからPC (Intel x86) server install CD（ま。ここは何にインストールするかによって変わるでしょうが・・・）という手順でダウンロードファイルまでたどり着けます。
が勝手に直リンク(^_^;)
x86系ならここから直接いけらー。
http://ftp.ecc.u-tokyo.ac.jp/UBUNTU-CDS/6.06/ubuntu-6.06.1-server-i386.iso (432MB)
これをCD-Rに焼きます。(イメージファイルとして焼いてください)

インストール
・インストールするPCをCD起動できるように設定したあとネットワークに接続してCD-Rを入れて起動。
(分からない人はここは無視して、CD-Rを入れてから再起動してみてWindowsが起動したらCD起動できないということですので、ここを参照して設定してください。)

Install LAMP server

を選択。

・言語は日本語を選択。あとは適当に質問に答えて、はいはいはい
(ホントはちゃんと設定しないといけないんだろうけどね)
途中、管理者のユーザー名とパスワードを聞かれるので、あらかじめ考えておいてください。


パッケージ管理の設定
とやらをします。

% sudo vi /etc/apt/sources.list

と入力してsources.list を編集。最初にsudoをするとパスワードを聞いてくるのでパスワードを入力。

# dev ...
# dev-src ...

で始まる行のコメントアウトをはずして

dev ...
dev-src ...

としする。保存。viの使い方が分からない人はここを見て要領をつかんだらコレをみながら試行錯誤してください。

とりあえずapt-getとやらをしてみる。

sudo apt-get update

んん■がいっぱい出力されてる・・・

2次のburg法で周波数推定の続きです。

2つの係数が与えられたときの周波数が分かったので、Rでこれを求める関数を作ってみましょう。後で窓関数を入れる事が出来るようにしておきます。

frq.burg <- function(x){
len <- length(x)
av_av <- mean(x)
x <- x - av_av
si <- 0.0
bo <- 0.0
for(k in 1:(len-2)){
# Normal Burg
W <- 1/(len-2)
# main
si <- si + (x[k] + x[k+2]) * x[k+1] *W
bo <- bo + x[k+1]*x[k+1] * W
}
a1 <- -si/(bo *2)
si <- 0.0
bo <- 0.0

for(k in 1:(len-2)){
# Normal Burg
W <- 1/(len-2)
# main
si <- si+(x[k]+a1*x[k+1])*(x[k+2]+a1*x[k+1]) *W
bo <- bo+((x[k]+a1*x[k+1])*(x[k]+a1*x[k+1])+(x[k+2]+a1*x[k+1])*(x[k+2]+a1*x[k+1]))*W
}
a2 <- - 2*si/bo
a1 <- (1 + a2 )* a1
peak <- (pi-acos(a1/(4*a2)+a1/4))/(2*pi)
return(peak)
}

ちょっと無駄がありますが..

使い方ですが、たとえば時系列 ts があったときに

freq <- ar.otburg(ts)*frequency(ts)

のようにして使います。

さて3次の場合を求めることができますか?
（4次まではなんとか出来ますが、やらないほうが身のためです。）

2次のarモデルのパラメータが与えられたとき、スペクトルは

`f = (P_m dt)/|1+sum_(m=1)^M a_(m) exp(i 2 pi m f)|`

で与えられます。つまり

`|1+a_(1) exp(i 2 pi f)+a_(2) exp(i 4 pi f)|`

を最小にするところが、スペクトルの最大を与えることに成ります。
MAXIMAでその点を求める計算をしてみましょう。(手計算でも出来るけど。)

1+a1*exp(%i*2*%pi*f)+a2*exp(%i*4*%pi*f);
abs(%)^2;
trigexpand(%)\$
trigsimp(%)\$
ratsubst((cos(2*%pi*f)+1)/2, cos(%pi*f)^2, %);
rat(%)\$
ratsubst(X, cos(2*%pi*f), %)\$
diff(%,X)\$
solve(%=0,X)[1]\$
ratsubst(cos(2*%pi*f),X,%);
solve(%,f)[1];

いつものように＄の前に￥が見えたら消してくださいね。(^_^;)
最初から`f`で微分しても良かったのですが、それでは例題として面白くないので、一回`X`を`cos(2 pi f)`に置き換え、その最小になる`X`を求めた上で（結局微分を使ってしまった^_^;）`X`から`f`を求めています。これを適当なファイルに保存してMaximaで解いてもらいましょう。

結果

((%i1) batch("*****.mac")\$
batching #p*****.mac
(%i2) a2*exp(%i*4*%pi*f)+a1*exp(%i*2*%pi*f)+1
(%o2) a2*%e^(4*%i*%pi*f)+a1*%e^(2*%i*%pi*f)+1
(%i3) abs(%)^2
(%o3) (a2*sin(4*%pi*f)+a1*sin(2*%pi*f))^2+(a2*cos(4*%pi*f)+a1*cos(2*%pi*f)+1)^2
(%i4) trigexpand(%)
(%i5) trigsimp(%)
(%i6) ratsubst((1+cos(2*%pi*f))/2,cos(%pi*f)^2,%)
(%o6) a2*(4*cos(2*%pi*f)^2+2*a1*cos(2*%pi*f)-2)+2*a1*cos(2*%pi*f)+a2^2+a1^2+1
(%i7) rat(%)
(%i8) ratsubst(X,cos(2*%pi*f),%)
(%i9) diff(%,X)
(%i10) solve(%=0,X)[1]
(%i11) ratsubst(cos(2*%pi*f),X,%)
(%o11) cos(2*%pi*f)=-(a1*a2+a1)/(4*a2)
(%i12) solve(%,f)[1]
\`solve' is using arc-trig functions to get a solution.
Some solutions will be lost.
(%o12) f=(%pi-acos(a1/(4*a2)+a1/4))/(2*%pi)

`f=(pi-acos(a_1/(4a_2)+a_1/4))/(2 pi)`

がスペクトルの最大を与えるということが分かります。

2次のarモデルに当てはめ
arsuiteiは僕の作った関数(もっとちゃんとした名前は無いのかという突っ込みは・・・)

> arsuitei(ts01M[1:300])
-1.826461 0.997897
> ar.burg(ts01M,order=2)

Call:
ar.burg.default(x = ts01M, order.max = 2)

Coefficients:
1 2
1.8264 -0.9977

Order selected 2 sigma^2 estimated as 0.003509

あらら？符号が逆。普通のarでやると

ar(ts01M,order=2)

Call:
ar(x = ts01M, order.max = 2)

Coefficients:
1 2
1.8196 -0.9914

Order selected 2 sigma^2 estimated as 0.01301

やっぱりこっちが正しいのかな(^_^;)
ちなみにDN Swingler - IEEE Trans. Acoustics, Speech & Signal, 1980を見てみると

`a_(22) = - X/Y ~= +1.`

になると書いてある。うぅぅぅぅぅぅぅん?

どうも符号のつけ方に2つの流儀があるようですね(^_^;)
日本でよく知られている日野本(日野, 「スペクトル解析」, 朝倉, 1977)に従うと、僕のやった計算と同じになります。

太陽系惑星「８個案」で調整…国際天文学連合総会 : 科学 : YOMIURI ONLINE（読売新聞）

"太陽系の惑星を９個から１２個に増やすとする惑星の定義案を示していた国際天文学連合（本部・パリ）は２３日、冥王（めいおう）星を惑星の地位から格下げし、８個にする修正案で最終調整に入った。
（略）

〈１〉自分の重力で球形になったもの
〈２〉その軌道領域で主要な天体であること――などを惑星とする新定義案で最終調整に乗り出した。

　新定義によると、軌道が海王星と一部重なり、ほかの惑星よりはるかに小さい冥王星は、惑星の座からすべり落ちることになる。「セレス」「カロン」も惑星の資格はなくなる。"

http://www.yomiuri.co.jp/science/news/20060823it04.htm

惑星の定義は、科学的な根拠より、輝かしい物理の発展に照らして行われるべきだと思います。もともと冥王星は、海王星や天王星の軌道のずれを説明するために第9惑星があるのではないかと考えられ、一生懸命探した結果見つかったもの。
だから惑星なんです。

結局、質量の推定量から、海王星や天王星の軌道のずれを説明することはできないことが分かっているんですが、どんなに探しても「第9の惑星」に該当するものは冥王星以外に無かった。

惑星としての冥王星は、私たちが、宇宙についてそして物理学そのものについて如何に理解していないかを鋭く指摘する存在として、大きな意味を持っていると思います。

IMPORTANT FOCUS NEWS
The Menace of Underground Ice
Nobel Focus: Discovery of Quasicrystals
The Purest X-Ray Beam
Cold Atoms Reveal Their Crystalline Na…
Tsunami Puzzle Explained
Rod-Shaped Nucleus
Sweet Music to your Nerves
Pioneer Anomaly is Fading
Ultrasound Cleans with a Twist
Random Roads Less Travelled

Nonlocal Memory Effects in the Dynamic…
Sub-Heisenberg Estimation Strategies A…
Efficiency of Isothermal Molecular Mac…
Lensing of 21-cm Fluctuations by Primo…
Helical Superconducting Black Holes
Search for Production of Resonant Stat…
Validity of Factorization of the High-…
New Class of Three-Body States
Plasmoid Ejection and Secondary Curren…
Dynamical Aspects of Strong Pinning of…

Unusual angular momentum transfer in e…
Soliton magnetization dynamics in spin…
Cavity-enhanced multiplexed comb spect…
Wentzel-Kramers-Brillouin analysis of …
Description and measurement of observa…
Adaptive Bayesian quantum tomography
Bose-Hubbard model for universal quant…

Entanglement dynamics in a non-Markovi…
Magnetic vortices on closely packed sp…
Spectroscopy of low-frequency noise an…
Anisotropic compression in the high-pr…
Neutron diffraction study of spin and …