phc

//この例は'levin'マクロ(対応するヘルプおよび例を参照)による
//結果と比較する際に有用です
//
//サイン波(1Hzおよび2Hz)にガウスノイズを加算することにより
//定義したプロセスを考えます(これが観測量です);
//シミュレーションされるプロセスのサンプリング周波数は10Hzです.

t=0:.1:100;rand('normal');
y=sin(2*%pi*t)+sin(2*%pi*2*t);y=y+rand(y);plot(t,y)

//yの共分散

nlag=128;
c=corr(y,nlag);

//共分散系列からハンケル行列を計算
// (より大きなnおよびmを指定することにより共分散からより多くの
// 情報を得ることもできます; 試して結果を比較してみてください!

n=20;m=20;
h=hank(n,m,c);

//マルコフ表現(mh,mf,mg)を計算
//ここでは,状態量の次元(4)を選択 :
//次元の推定の方がむしろ困難な問題です!
//nsを変えてみてください!
//(観測量の次元はここでは1にします)

ns=4;
[mh,mf,mg]=phc(h,1,ns);

//mfのスペクトルが観測量プロセスyの周波数スペクトルを
//有することを確認
//(この例では,yは10Hz(T=0.1s)で標本化されたことを思い出してください.
//このため,記録から元の周波数を回復し,周波数サンプリングにより正しく
//スケーリングする必要があります)

s=spec(mf);s=log(s);
s=s/2/%pi/.1;

//ここで,スペクトル推定値を得ました
imag(s),