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年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 | このブログの読者になる | 更新情報をチェックする

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月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 | このブログの読者になる | 更新情報をチェックする

Physical Review

Recent stories in Physical Review Focus

広告


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

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

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


×

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