% ARMA_CHANGE.M
%
% P. Flandrin, Sept. 25, 2003
%
% visualise un changement de mod�le ARMA d'ordres
% 0, 1, 2, 3 ou 4 par calcul de l'erreur de pr�diction
%
% calls  : - ARMA_coeff.m

N = 500; % nombre de points du processus

p1 = 4; % ordre de la partie AR du mod�le 1 (<= 4)
q1 = 0; % ordre de la partie MA du mod�le 1 (<= 4)

p2 = 2; % ordre de la partie AR du mod�le 2 (<= 4)
q2 = 0; % ordre de la partie MA du mod�le 2 (<= 4)

n = 1; % num�ro de la 1�re figure de placement des p�les

% ----- mod�le 1

[A1,B1] = ARMA_coeff(p1,q1,n);

b1 = randn(1,N);
x1 = filter(B1,A1,b1);

% ----- mod�le 2

[A2,B2] = ARMA_coeff(p2,q2,n+1);

b2 = randn(1,N);
x2 = filter(B2,A2,b2);
x2 = x2/sqrt(mean(x2.^2))*sqrt(mean(x1.^2));

figure(n+2)
clf

r = fix(N/4 + N/2*rand);
x = x1;
x(r+1:N) = x2(r+1:N);

subplot(313)
plot(x)
axis([1 N 1.1*min(x) 1.1*max(x)])
Ax = axis;

pause

hold on
plot([r r],[min(x) max(x)],'r')
hold off

pause

subplot(311)
plot(1:r,x1(1:r))
axis(Ax)
hold on
plot(r+1:N,x1(r+1:N),'--')
plot([r r],[min(x) max(x)],'r')
hold off

subplot(312)
plot(1:r,x2(1:r),'--')
axis(Ax)
hold on
plot(r+1:N,x2(r+1:N))
plot([r r],[min(x) max(x)],'r')
hold off

pause

xp = filter(A1,B1,x);
ep = xp.^2;

figure(n+3)

subplot(211)
plot(x)
hold on
plot([r r],[min(x) max(x)],'r')
title('signal')
hold off

subplot(212)

d = 2;

for k = 1:d:N
	
	plot(1:d:k,ep(1:d:k))
	hold on
	plot([1 N],[3.8 3.8],'r--')
	plot([1 N],[6.6 6.6],'r:')
	axis([1 N 0 1.1*max(ep(1:d:N))])
	title('puissance de l erreur de prediction -- seuil de fausse alarme � 95 % (--) et 99 % (..)')
	pause(.01)
	hold off
	
end