2006年08月22日

Data-Adaptive Weighted Burg法で周波数推定

昨日の続きです。きれいな正弦波が入ってきたときは線形予測フィルタの`a_(11)`から`a_(22)`を解析的に計算できるので、予測フィルタの次数は1で十分。
ということで、前よりも簡単に周波数推定ができます。
Burg法には位相の問題があったのですがData-Adaptive Weighted Burg法を用いると解決できます。

http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1165046

を参考にしてまず1次の線形予測フィルタの係数を計算します。

Rscript
spdab2 <- function(x){
len <- length(x)
S10 <- 0.0
S11 <- 0.0
for(k in 1:(len-2)){
W <- x[k+1]*x[k+1]
# Sij <- x[k+i] * x[k+j] + x[k+2-i] *x[k+2-j]
S10 <- S10 + x[k+1] * x[k] + x[k+1] * x[k+2]
S11 <- S11 + x[k+1] * x[k+1] +x[k+1] * x[k+1]
}
a11 <- -S10/S11
return(a11)
}


ここで、昨日考えたような波が入ってくると、
`a_(11)= - cos theta `

となるので、あとは時系列データfsに対して
fn <- acos(-spdab2(ts))/(2*pi) * frequency(ts)

とするだけ。おしまい。

ちょっと遅いです。高速化はまたいつか。
続きを読む
posted by metola at 17:16| Comment(0) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする

2006年08月21日

Burgアルゴリズムで周波数推定

DN Swingler - IEEE Trans. Acoustics, Speech & Signal, 1980
を見てみると、時系列データ
`x(k)=cos(theta k + phi)`
に対して2次のPrediction Error Filter (PEF) のパラメータを計算すると
`2a_(11) = -2 cos (theta - delta)`
となる。ここで`delta`の・はメンドいから説明しない。上記論文を読んで。


というわけで、適当に周波数推定ができるか試してみた。


Rにはar.burgという関数がある。burg法を使って時系列データにarモデルを当てはめてくれるという代物。`theta = 2 pi f \/ N_s`(`N_s`は1秒あたりのサンプリング数)になることを考慮に入れつつ(符号は気にしない♪←コラ)こいつを使って周波数推定をしてみる。

Freq <- acos(ar.burg(ts, order.max=2)\$ar[1]/2)/(2*pi) * frequency(ts)


・・・なんかうまくいってます。($の前の¥は消してくださいね^_^;)
ちょーてきとーなのでご利用は計画的に。
posted by metola at 12:48| Comment(0) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする

2006年08月18日

spec.ar

Rのspec.ar関数。なにかで修正しなきゃいけんかったんじゃけど、なぜ修正が必要だったか忘れた。

spec.ar <- function (x, n.freq, order = NULL, plot = TRUE, na.action = na.fail,
method = "yule-walker", ...)
{
if (!is.list(x)) {
series <- deparse(substitute(x))
x <- na.action(as.ts(x))
xfreq <- frequency(x)
nser <- NCOL(x)
x <- ar(x, is.null(order), order, na.action = na.action,
method = method)
}
else {
cn <- match(c("ar", "var.pred", "order"), names(x))
if (any(is.na(cn)))
stop("'x' must be a time series or an ar() fit")
series <- x\$series
xfreq <- x\$frequency
if (is.array(x\$ar))
nser <- dim(x\$ar)[2]
else nser <- 1
}
order <- x\$order
if (missing(n.freq))
n.freq <- 500
freq <- seq(0, 0.5, length = n.freq)
if (nser == 1) {
coh <- phase <- NULL
if (order >= 1) {
cs <- outer(freq, 1:order, function(x, y) cos(2 *
pi * x * y)) %*% x\$ar
sn <- outer(freq, 1:order, function(x, y) sin(2 *
pi * x * y)) %*% x\$ar
spec <- x\$var.pred/(xfreq * ((1 - cs)^2 + sn^2))
}
else spec <- rep(x\$var.pred/(xfreq), length(freq))
}
else .NotYetImplemented()
freq<-freq*xfreq
spg.out <- list(freq = freq, spec = spec, coh = coh, phase = phase,
n.used = nrow(x), series = series, method = paste("AR (",
order, ") spectrum ", sep = ""))
class(spg.out) <- "spec"
if (plot) {
plot(spg.out, ci = 0, ...)
return(invisible(spg.out))
}
else return(spg.out)
}



$の前に¥が見えたら消してくださいね(^_^;)
posted by metola at 14:35| Comment(0) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする

2006年08月10日

C# GSL 覚書

CsharpからのGSLの利用
http://www.tmps.org/index.php?Csharp%A4%AB%A4%E9%A4%CEGSL%A4%CE%CD%F8%CD%D1

うぅぅぅぅん
ただFFTがしたいだけのためにここまでするのはめんどくさそうだ(^_^;)



関係ないけど面白そう。でも英語がやたら読みにくいのはなぜ?・・・というか開発(^_^;)
http://numerator.sourceforge.net/
posted by metola at 21:31| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年08月09日

覚書

講義ノートが面白い
http://www.tac.tsukuba.ac.jp/~yaoki/
posted by metola at 23:41| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年08月03日

LaTeXMathMLって・・・

久しぶりにASCIIMathMLのサイトに行ったら子分のLaTeXMathMLというのができていた。

TeXで書いた原稿をそのままぺたっと貼り付けて投稿したら数式入りのエントリーの出来上がりという優れもの。ここも導入しようか知らん。
ただASCIIMathMLですらかなりマニアックだからなぁ。

実は、携帯電話で見たときの事を考えて、ここではTeX形式で式を書かず、アクサングラーブではさんだところを数式とみなして処理してます。べた書き段階での可読性としてはASCIIMathMLはすごくいい。TeXを使ったことのある人ならわかるとおもうけど、TeXで書き始めると\やらわけわからんコマンドが数式の中に入って、慣れるまではすらすらと読めないんですよねぇ。

たとえば分数も
TeXだったら
\frac{1}{\sqrt{2}}

と書くところをASCIIMathMLなら
1/sqrt(2)

でOKという感じ。ちなみに上の式を数式にすると
`1/sqrt(2)`

です。ASCIIMathMLの方が読みやすくないですか?

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

2006年08月01日

検証:テポドンは何キロ飛んだ?

高校物理の問題でといてみよー♪


1.ミサイルは角度`theta`傾いて発射
2.その角度を保ったまま等加速度運動をして発射から40秒後、高度1.5kmに達した。またそのときの飛行距離は、発射台から見て`x_1` kmであった。
3.その直後、ミサイルは制御不能に陥り爆発した。

このとき以下の問いに答えなさい。
(1) 40秒後の`x`(横)方向の速度`v_x`と`y`(高さ)方向の速度`v_y`を km/s の単位で求めなさい。
(2) 爆発後、落下するまでにかかる時間`t'`を求めなさい。
(3) 最終的にミサイルの残骸が落下する位置は`x_1`の何倍になるか計算しなさい。
(4) `theta`が60度、45度、30度のときの`x_2`を計算しなさい。

ただし、ミサイルの運動に伴う空気抵抗は考えないものとし、重力加速度`g=`10/1000 km/s2とする。

問題設定は下記の図を参照してください。



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

2006年06月30日

Google Scholarが…

Google Scholarが…
google.png

日本語対応になっている!
posted by metola at 13:40| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年06月29日

気になる言語

GNU Sather
http://www.gnu.org/software/sather/
http://gnu.hands.com/brave-gnu-world/issue-18.ja.html
http://www.geocities.jp/shido_takafumi/sather/index.html
↑ではwindowsはないと書いてあるけど実際にはあるみたい。
ftp://ftp.icsi.berkeley.edu/pub/sather/Sather-1.1-Win32.ZIP
posted by metola at 22:02| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年06月23日

TiddlyWikiで知的財産権法の条文集

TiddlyWikiを使っていて思うんですが、これで条文集を作ったら便利だろうなぁと思いました。んで、実際にやってみたのがこれ。
IPtiddly.zip
解凍して、index.htmlをご覧ください。マシンのスペックによっては止まるかも。あとTiddlyWikiはjavascriptの裏技のバケモノのようなものなのでブラウザはFireFoxをお勧めします。

収録してある条文は
本 特許法
本 実用新案法
本 意匠法
本 商標法

条文は一部RONの六法全書 on LINEからコピペさせていただきました。(^_^;)

ほかにもTiddlyWiki周りでお世話になったサイトがたくさんあるのですが、
http://hsj.jp/junknews/archives/tiddlywiki_susume.html
http://flow.dip.jp/mt/archives/u/twmemo.html
の他たくさんありすぎて覚えていないので割愛させてください。

特許法を最初に作ったので、タグふりをまじめにやってて、実用新案法は最後にやったので、慣れてきてきれいにできているというのは無視の方向で・・・

変なところがたくさんあると思いますが、あとはご自身で編集してくださいませ。(サポート放棄w)
あ、でもコメントにここが変!と書いてくださると個人的に助かります。
posted by metola at 14:15| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

Physical Review

Recent stories in Physical Review Focus

広告


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

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

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


×

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