%  arma.m
%
% -----------------------------------
%   AR model から ARMA model へ
% ----------------------------------
%
%  ある ARMA過程 を想定して
%
Fs = 1000;
N = 500;
%a = [1 -1.3136  1.4401  -1.0919  0.83527];　「計数・測定」p190
%b = [1  0.13137 0.023543 0.10775 0.03516];
a = [1 -1.31 1.44 -1.10 0.84];
b = [1  0.13 0.02  0.11 0.04];
h = impz(b,a,N);            %  そのインパルス応答 　　    
subplot(521);plot(h)
legend('given impulse response')
%
f = freqz(b,a,N);           %  その周波数応答　　
subplot(522);semilogy(abs(f));
axis([0 500 1 100]);legend('given frequency response')
%
x = filter(b,a,1*randn(1024,1));            %  サンプル時系列生成
subplot(523);plot(x);
axis([0 1000 -20 20]);legend('sample')
[pxx,fr] = psd(x,256,Fs,hanning(256),128);  %  FFT法によるスペクトル推定
subplot(524);semilogy(fr,sqrt(pxx));
axis([0 500 1 100]);legend('FFT spectrum')
%
%  AR過程（次数無限大）として同定
%
ai = aryule(x,50);            % Yule-Walker法  ／次数(50)
%
hai=impz(1,ai,N);   % インパルス応答
subplot(525);plot(impz(1,ai,N))
legend('AR model')
f = freqz(1,ai,N);
subplot(526);semilogy(abs(f));
legend('AR model');axis([0 500 1 100])
%
%  ARモデル　-＞　ARMAモデル
%
%          ARMAモデル　　（prony法）
%
[bci,aci] = prony(hai,4,4)     % 次数(4,4)
subplot(527);plot(impz(bci,aci,N));
legend('ARMA (prony)')
f = freqz(bci,aci,N);
subplot(528);semilogy(abs(f));axis([0 500 1 100])
%
%          ARMAモデル　　（Steiglitz-McBride法）
%
[bs,as] = stmcb(hai,4,4,10)    % 次数(4,4); 10 iterations　
subplot(529);plot(impz(bs,as,N));
legend('ARMA (Steiglitz-McBride)')
f = freqz(bs,as,N);
subplot(5,2,10);semilogy(abs(f));axis([0 500 1 100])
%--------------------------------------------------------
%                                2000.6.13  by K.Tsukada
%                                2002.6.24  revised 

