% ARMA_PZ.M
%
% P. Flandrin, 24 sept. 2003
%
% trace la DSP d'un processus ARMA ˆ partir de son
% diagramme p™les-zŽros

subplot(211)
clf
subplot(212)
clf

p = 6; % nombre de p™les (>= 0)
q = 6; % nombres de zŽros (>= 0)

N = 1000;
theta = linspace(0,pi,N);
f = linspace(0,.5,N);
z = exp(i*2*pi*f);

x = cos(theta);
y = sin(theta);

figure(1)

subplot(211)

	plot([-1 1],[0 0])
	hold on
	plot(x,y)
	plot([0 0],[0 0],'+')
	plot(x/2,y/2,':')
	plot(3*x/4,3*y/4,':')
	plot(7*x/8,7*y/8,':')
	plot(15*x/16,15*y/16,':')
	plot(31*x/32,31*y/32,':')
	plot([0 cos(pi/16)],[0 sin(pi/16)],':')
	plot([0 cos(3*pi/16)],[0 sin(3*pi/16)],':')
	plot([0 cos(5*pi/16)],[0 sin(5*pi/16)],':')
	plot([0 cos(7*pi/16)],[0 sin(7*pi/16)],':')
	plot([0 cos(9*pi/16)],[0 sin(9*pi/16)],':')
	plot([0 cos(11*pi/16)],[0 sin(11*pi/16)],':')
	plot([0 cos(13*pi/16)],[0 sin(13*pi/16)],':')
	plot([0 cos(15*pi/16)],[0 sin(15*pi/16)],':')
	plot([0 cos(pi/8)],[0 sin(pi/8)],':')
	plot([0 cos(pi/4)],[0 sin(pi/4)],':')
	plot([0 cos(3*pi/8)],[0 sin(3*pi/8)],':')
	plot([0 cos(pi/2)],[0 sin(pi/2)],':')
	plot([0 cos(5*pi/8)],[0 sin(5*pi/8)],':')
	plot([0 cos(3*pi/4)],[0 sin(3*pi/4)],':')
	plot([0 cos(7*pi/8)],[0 sin(7*pi/8)],':')
	axis([-1.1 1.1 -0.1 1.1])
	set(gca,'XTick',[])
	set(gca,'YTick',[])
	axis equal

S = 100.*ones(1,N);

if p == 0
	
	S = S;
	
else

	Sp = [];

	for k = 1:p
	
		subplot(211)
		title(['p™le ',int2str(k),'/',int2str(p)])
		[X,Y] = ginput(1);
		Z = X + i*Y;
		plot(X,Y,'r.','LineWidth',6)
		hold on
	
		subplot(212)
		S1 = 1./(abs(1-Z./z)).^2;
		Sp = [Sp' S1']';
		plot(f,log(Sp),'r-')
		hold on
		S = S.*S1;
		plot(f,log(S),'k','LineWidth',2)
		xlabel('frequence')
		ylabel('log(DSP')
		title(['ARMA(',int2str(2*k),',',int2str(0),')'])
		hold off
		
	end

end

if q == 0
	
	S = S;
	
	else

		Sq = [];
		
		for k = 1:q
	
		subplot(211)
		title(['zero ',int2str(k),'/',int2str(q)])
		[X,Y] = ginput(1);
		Z = X + i*Y;
		plot(X,Y,'bo','LineWidth',4)
		hold on
	
		subplot(212)
		S1 = (abs(1-Z./z)).^2;
		Sq = [Sq' S1']';
		plot(f,log(Sq),'b-')
		hold on
		
		if p == 0
			
		else
		
			plot(f,log(Sp),'r-')
			
		end
		
		S = S.*S1;
		plot(f,log(S),'k','LineWidth',3)
		xlabel('frequence')
		ylabel('log(DSP')
		title(['ARMA(',int2str(2*p),',',int2str(2*k),')'])
		hold off
	
	end

end

if p == 0 & q == 0

	subplot(212)
	plot(f,log(S),'k','LineWidth',3)
	xlabel('frequence')
	ylabel('log(DSP')
	title(['ARMA(',int2str(p),',',int2str(q),')'])
	
end

subplot(211)
hold off
subplot(212)
hold off