2006年12月22日

Buneman周波数推定

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

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

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

目的
DFTのデータから信号の周波数`beta`を推定する。

理論
`m`点DFTの計算結果から、信号の周波数`beta`が、周波数`f_{b}`から`f_{b+1}`の間にあることが分かっています。`f_{b}`と`f_{b+1}`におけるフーリエ変換の値`F_{b}`と`F_{b+1}`の絶対値はそれぞれ解析的に次のように求めることが出来ます。

`|F_{b}| = sin( pi * ( beta - b ) ) / sin ( pi * ( beta - b )/m)`
`|F_{b+1}| = sin( pi * (beta - b - 1))/ sin ( pi * ( beta - b - 1)/m)`

`delta = beta - b `と置き、`|F_{b}|`/`|F_{b+1}|`を計算すると
`|F_{b}|/|F_{b+1}| =sin( pi delta ) * sin( (pi ( delta - 1 ) )/ m ) / ( sin ( pi ( delta - 1 ) ) * sin( (pi delta) / m))`

これを`delta`について解けばよいわけです。
ちょっとこれは面倒なので、MAXIMAを使いながらやってみましょう。

MAXIMAで計算

(%i2) F_b:sin(%pi*(%beta-b))/sin(%pi*(%beta-b)/m)
(%i3) F_c:sin(%pi*(-1-b+%beta))/sin(%pi*(-1-b+%beta)/m)
(%i4) FF = ratsubst(%delta,%beta-b,F_b/F_c)
(%o4) FF = sin(%pi*%delta)*sin(%pi*(%delta-1)/m)
/(sin(%pi*(%delta-1))*sin(%pi*%delta/m))
(%i5) expand(%)
(%i6) trigexpand(%)
(%i7) expand(%)
(%i8) trigreduce(%)
(%o8) FF = sin(%pi/m)*cot(%pi*%delta/m)-cos(%pi/m)
(%i9) cos(%pi/m)+%
(%i10) 1/%
(%i11) trigreduce(%)
(%o11) 1/(FF+cos(%pi/m)) = csc(%pi/m)*tan(%pi*%delta/m)
(%i12) solve(%,%delta)[1]
`solve' is using arc-trig functions to get a solution.
Some solutions will be lost.
(%o12) %delta = m*atan(1/(csc(%pi/m)*FF+cos(%pi/m)*csc(%pi/m)))/%pi
(%i13) trigsimp(%)
(%o13) %delta = m*atan(sin(%pi/m)/(FF+cos(%pi/m)))/%pi
(%i14) ratsubst(%beta-b,%delta,%)
(%o14) %beta-b = m*atan(sin(%pi/m)/(FF+cos(%pi/m)))/%pi
(%i15) solve(%,%beta)[1]
(%o15) %beta = (m*atan(sin(%pi/m)/(FF+cos(%pi/m)))+%pi*b)/%pi
(%i16) expand(%)
(%o16) %beta = m*atan(sin(%pi/m)/(FF+cos(%pi/m)))/%pi+b

最後の(%o16)を見れば分かるように、周波数`beta`は
`beta = m/pi*tan^(-1)(sin( pi/m)/(|F_{b}|/|F_{b+1}|+cos( pi/m)))+b`

によって計算することができます。

参考文献
DH Bailey, PN Swarztrauber, "The Fractional Fourier Transform and Applications", SIAM Review, 33(3) 189-404 (1991)
ラベル:MAXIMA
posted by metola at 15:11| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
コメントを書く
お名前: [必須入力]

メールアドレス:

ホームページアドレス:

コメント: [必須入力]

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


※画像の中の文字を半角で入力してください。

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

Physical Review

Recent stories in Physical Review Focus

×

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