2006年09月11日

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) | 日記 | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
aetheogamous unrevered internalization atmospherology untied adobe tridiapason osmoscope
<a href= http://www.sla.purdue.edu/WAAW/Palmquist/index.htm >Women in Photography Archive</a>
http://www.dreammaker-remodel.com/

きのこな夜
Posted by Iris Morrow at 2007年12月17日 16:51

オタ系の子いいねーwwwww
ちょっとコスプレしてくれって言ったらネコ耳メ イ ドになってくれたしなぁwww
本人もノリノリでしまいには喘ぎながら「にゃーん」とか言ってたしなwww(*゚∀゚)・∵. ブハッ!!
あーでもこれかなり八マるねwwwちなみ次は巫女のカッコさせる予定wwwww

http://EnvI.StrowcruE.net/i1pr5sy/
Posted by ネ コ 耳 モード突入っ!! at 2009年08月18日 21:21

夏も終わりだから寂しいコが多いのかなー?めっさメールくるんだけどwww

・・・で、昨日気づいたんだが8月で会った人数が合計30人を越えてたよ!!Σ( ̄△ ̄;)
1日1人ペースどころじゃねぇなwwwつーかヤりすぎだろ俺wwwww
いや〜ちょっとコレ31日までに1 0 0 マ ソいくんじゃねぇの!?www(σ゚∀゚)σワクワク

http://OmeTroRo.com/Puru/0xal1cj/
Posted by 受信フォルダ整理めんどくせぇ! at 2009年08月26日 01:55

ちょwww奥さんwwwww1人だって言ってたのに3人でした♪とかちょっと待ってよwww
気がついたらプェラ・ア@ルなめ・ベロッチューの三連発にずっとあえぎっぱなしだわwwwww
家に帰ったら1キロ以上体重減ってたしなwww(ヽ゚Д゚/)ゲソゲソー

ま〜3人合わせて20マソくれたから大満足だけどwwwww

http://soiE.ikiSugI.com/a1bg8xl/
Posted by 主婦のトライアングルw at 2009年09月02日 03:38

いや〜うんまそーなマヌ-コだったから1時間くらいずっと夢中でペロッペロしちゃったよwww
ヨダレ垂らしてすんごいアヘ顔になってたけど、オレ全然やめなかったしねwww

おかげでホントは報 酬 3 万だったんだけど6 万にグレードアップ♪♪♪
オレのナメ好きもたまには役に立つんだなwwwww
(●´艸`)ムフフ

http://Erin%2eYumenokuni%2eNet/k3dppit/
Posted by ついつい・・・ at 2009年09月05日 07:26
コメントを書く
お名前: [必須入力]

メールアドレス:

ホームページアドレス:

コメント: [必須入力]

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


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

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

Physical Review

Recent stories in Physical Review Focus

×

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