2006年12月25日

アクティブ・フィルタの設計は・・・

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

これで、高次のフィルタがあっという間にできるぞ♪
posted by metola at 22:54| Comment(21) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年12月22日

Buneman周波数推定

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

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

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

続きを読む
タグ:MAXIMA
posted by metola at 15:11| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年12月20日

QucsでLPFを作ってみた。

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

タグ:Qucs
posted by metola at 23:53| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年12月18日

Rで連番処理

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)
posted by metola at 11:45| Comment(0) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする

2006年10月09日

Ubuntuサーバーと格闘(2) 〜日本語をなんとか

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のマニュアル)
posted by metola at 15:44| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

Ubuntuサーバーと格闘(1) 〜インストール

古いパソコンに、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


んん■がいっぱい出力されてる・・・がく〜(落胆した顔)
posted by metola at 14:39| Comment(0) | TrackBack(0) | 管理 | このブログの読者になる | 更新情報をチェックする

2006年09月11日

2次のburg法で周波数推定(2)

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次まではなんとか出来ますが、やらないほうが身のためです。)
posted by metola at 13:28| Comment(17) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする

2次のburg法で周波数推定

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)`

がスペクトルの最大を与えるということが分かります。
タグ:MAXIMA
posted by metola at 13:02| Comment(5) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年09月06日

Rのar.burgが変

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)に従うと、僕のやった計算と同じになります。
posted by metola at 17:47| Comment(0) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする

2006年08月24日

惑星の数は?

太陽系惑星「8個案」で調整…国際天文学連合総会 : 科学 : YOMIURI ONLINE(読売新聞)
"太陽系の惑星を9個から12個に増やすとする惑星の定義案を示していた国際天文学連合(本部・パリ)は23日、冥王(めいおう)星を惑星の地位から格下げし、8個にする修正案で最終調整に入った。
(略)

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

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

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

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

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

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

posted by metola at 11:18| Comment(0) | TrackBack(1) | 日記 | このブログの読者になる | 更新情報をチェックする

Physical Review

Recent stories in Physical Review Focus

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。