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

Physical Review

Recent stories in Physical Review Focus

広告


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

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

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


×

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