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)

とするだけ。おしまい。

ちょっと遅いです。高速化はまたいつか。
・・・・Data-AdaptiveとかいいながらWの値を掛けるの忘れてました。でもなぜかこっちのほうが周波数がよく決定できる。うぅぅぅぅぅぅん?
posted by metola at 17:16| Comment(0) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
コメントを書く
お名前: [必須入力]

メールアドレス:

ホームページアドレス:

コメント: [必須入力]

認証コード: [必須入力]


※画像の中の文字を半角で入力してください。
この記事へのトラックバックURL
http://blog.seesaa.jp/tb/22668625

この記事へのトラックバック

Physical Review

Recent stories in Physical Review Focus

×

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