% MEAN_CHANGE.M
%
% P. Flandrin, Sept. 25, 2003
%
% d�tecte un saut de moyenne par somme cumul�e
%
% calls : - AR_coeff.m

N = 500; % nombre de points du processus
H = .3; % amplitude du saut de moyenne % .3
saut = zeros(1,N);
h = H; % amplitude de saut mini � d�tecter
thresh = 15;

% ----- mod�le

p = 4; % ordre du mod�le
n = 1;

figure(n)
clf

A1 = AR_coeff(p,n);

b1 = randn(1,N+5*p);
xx = filter(1,A1,b1);
x1 = xx(5*p+1:N+5*p);
x1 = x1/sqrt(mean(x1.^2));

figure(n+1)

clf

r = fix(N/4 + N/2*rand);
x = x1;
x(r+1:N) = x1(r+1:N)+H;
saut(r+1:N) = ones(1,N-r);

subplot(313)
plot(x)
title(['signal AR(',int2str(p),') avec saut de moyenne'])

pause

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

subplot(311)
plot(x1)
title(['signal AR(',int2str(p),') sans saut de moyenne'])
Ax = axis;

subplot(312)
plot(saut,'r')
title('saut')
axis(Ax)

pause

y = x-h/2;
ys = cumsum(y);

T = [];

for k = 1:N
	
	T(k) = ys(k) - min(ys(1:k));
	
end

ind = find(T>=thresh);
ysmin = min(ys(1:ind(1)));
re = find(ys == ysmin);

figure(n+3)

subplot(311)
plot(x)
hold on
plot([r r],[min(x) max(x)],'r')
title(['signal AR(',int2str(p),') avec saut de moyenne'])
hold off

d = 4;
m = (max(ys)-min(ys))/20;

for k = 1:d:N
	
	subplot(312)

	plot(1:d:k,ys(1:d:k))
	hold on
	axis([1 N min(ys)-m max(ys)+m])
	title('somme cumulee')
	pause(.1)
	hold off
	
	subplot(313)

	plot(1:d:k,T(1:d:k))
	hold on
	plot([1 N],[thresh thresh],'r--')
	title('statistique de decision')
	axis([1 N 0 1.1*max(T)])

	if k >= r
		
		plot([r r],[0 1.1*max(T)],'r--')
		
	end
		
	if k >= ind(1)
		
		plot([re re],[0 1.1*max(T)],'r')
		
	end
	
	pause(.01)
	hold off
	
end