Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

複素数データを使用したスペクトル推定 - Marple のテスト ケース

この例では、時系列データでスペクトル推定を実行する方法を説明します。ここでは、Marple のテスト ケースを使用します (The complex data in L. Marple: S.L. Marple, Jr, Digital Spectral Analysis with Applications, Prentice-Hall, Englewood Cliffs, NJ 1987)。

テスト データ

それではまず、テスト データを読み込んでみましょう。

load marple

System Identification Toolbox™ のルーチンの大部分では複素数データを使用できます。しかし、ここではプロットのためデータの実数部と虚数部を分けて調べることにします。

まず、データの概要を見ます。

subplot(211),plot(real(marple)),title('Real part of data.')
subplot(212),plot(imag(marple)),title('Imaginary part of data.')

予備解析ステップとして、データのピリオドグラムを確認します。

per = etfe(marple);
w = per.Frequency;
clf
h = spectrumplot(per,w);
opt = getoptions(h);
opt.FreqScale = 'linear';
opt.FreqUnits = 'Hz';
setoptions(h,opt)

データ レコードが 64 サンプルのみで、ピリオドグラムが 128 周波数に算出されているため、明らかに狭周波数ウィンドウからの振動が見られます。このため、ピリオドグラムに何らかの円滑化を適用します (周波数分解能 1/32 Hz に対応)。

sp = etfe(marple,32);
spectrumplot(per,sp,w);

ここで、ブラックマン-テューキー法でスペクトル推定を行います。

ssm = spa(marple); % Function spa performs spectral estimation
spectrumplot(sp,'b',ssm,'g',w,opt);    
legend({'Smoothed periodogram','Blackman-Tukey estimate'});

既定のウィンドウ長では、この少量のデータに非常に狭いラグ ウィンドウが示されます。以下のように、大きなラグ ウィンドウを選択できます。

ss20 = spa(marple,20);
spectrumplot(sp,'b',ss20,'g',w,opt);
legend({'Smoothed periodogram','Blackman-Tukey estimate'});

自己回帰 (AR) モデルの推定

パラメトリック 5 次 AR モデルは、以下のように算出されます。

t5 = ar(marple,5);

ピリオドグラム推定と比較します。

spectrumplot(sp,'b',t5,'g',w,opt); 
legend({'Smoothed periodogram','5th order AR estimate'});

AR コマンドは、スペクトル推定で 20 の異なる方法に対応します。前述の方法は、Marple の著書で '修正共分散推定' と呼ばれています。

その他のよく知られた方法では、以下のように取得します。

tb5 = ar(marple,5,'burg');      % Burg's method
ty5 = ar(marple,5,'yw');        % The Yule-Walker method
spectrumplot(t5,tb5,ty5,w,opt);
legend({'Modified covariance','Burg','Yule-Walker'})

操作変数法による AR モデルの推定

また、AR モデル化は、操作変数法 (Instrumental Variable Approach) でも実行できます。この場合、関数 ivar を使用します。

ti = ivar(marple,4); 
spectrumplot(t5,ti,w,opt);
legend({'Modified covariance','Instrumental Variable'})

スペクトルの自己回帰移動平均 (ARMA) モデル

さらに、System Identification Toolbox は、スペクトルの ARMA モデル化にも対応しています。

ta44 = armax(marple,[4 4]); % 4 AR-parameters and 4 MA-parameters
spectrumplot(t5,ta44,w,opt);
legend({'Modified covariance','ARMA'})