% ARMA_COEFF.M
%
% P. Flandrin, Sept. 30, 2003
%
% calcule les coefficients d'un mod�le ARMA
% d'ordres 0, 1, 2, 3 ou 4 � partir du placement 
% de ses p�les et z�ros
%
% input  : - p : ordre de la partie AR du mod�le
%          - q : ordre de la partie MA du mod�le
%          - n : num�ro de la figure
% output : - A : vecteur [1 x p+1] des coeffs AR
%          - B : vecteur [1 x p+1] des coeffs MA

function [A,B] = ARMA_coeff(p,q,n);

N = 500;
theta = linspace(0,pi,N);
x = cos(theta);
y = sin(theta);

A = zeros(1,p+1);
B = zeros(1,q+1);

figure(n)
clf

if p == 0;
	
	A(1) = 1;
	
elseif p == 1

	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
	
	title('pole r�el')
	[X,Y] = ginput(1);
	plot(X,Y,'r.','LineWidth',6)
	
	A(1) = 1;
	A(2) = -X;
	
	hold off
	
elseif p == 2

	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
	
	title('pole complexe ')
	[X,Y] = ginput(1);
	plot(X,Y,'r.','LineWidth',6)

	Z = X +i*Y;
	
	A(1) = 1;
	A(2) = -2*real(Z);
	A(3) = (abs(Z)).^2;
	
	hold off
	
elseif p == 3

	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
	
	title('pole r�el')
	[X,Y] = ginput(1);
	plot(X,Y,'r.','LineWidth',6)
	
	R = X;
	title('pole complexe ')
	[X,Y] = ginput(1);
	plot(X,Y,'r.','LineWidth',6)

	Z = X +i*Y;
	
	A(1) = 1;
	A(2) = -(R+2*real(Z));
	A(3) = 2*R*X+(abs(Z)).^2;
	A(4) = -R*(abs(Z)).^2;

	hold off
	
elseif p == 4

	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
	
	X = [];
	Y = [];
	
	for k = 1:2
		
		title(['pole complexe ',int2str(k),'/2'])
		[X(k),Y(k)] = ginput(1);
		plot(X,Y,'r.','LineWidth',6)
		
	end

	Z = X +i*Y;
	
	A(1) = 1;
	A(2) = -2*real(Z(1)+Z(2));
	A(3) = (abs(Z(1))).^2 + (abs(Z(2))).^2 +4*real(Z(1))*real(Z(2));
	A(4) = -2*((abs(Z(1))).^2*real(Z(2))+(abs(Z(2))).^2*real(Z(1)));
	A(5) = (abs(Z(1)*Z(2))).^2;
	
	hold off
	
end

hold on

if q == 0;
	
	B(1) = 1;
	
elseif q == 1

	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
	
	title('zero reel')
	[X,Y] = ginput(1);
	plot(X,Y,'bo','LineWidth',2)
	
	B(1) = 1;
	B(2) = -X;
	
	hold off
	
elseif q == 2

	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
	
	title('zero complexe ')
	[X,Y] = ginput(1);
	plot(X,Y,'bo','LineWidth',2)

	Z = X +i*Y;
	
	B(1) = 1;
	B(2) = -2*real(Z);
	B(3) = (abs(Z)).^2;
	
	hold off
	
elseif q == 3

	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
	
	title('zero reel')
	[X,Y] = ginput(1);
	plot(X,Y,'bo','LineWidth',2)
	
	R = X;
	title('zero complexe ')
	[X,Y] = ginput(1);
	plot(X,Y,'bo','LineWidth',2)

	Z = X +i*Y;
	
	B(1) = 1;
	B(2) = -(R+2*real(Z));
	B(3) = 2*R*X+(abs(Z)).^2;
	B(4) = -R*(abs(Z)).^2;

	hold off
	
elseif q == 4

	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
	
	X = [];
	Y = [];
	
	for k = 1:2
		
		title(['zero complexe ',int2str(k),'/2'])
		[X(k),Y(k)] = ginput(1);
		plot(X,Y,'bo','LineWidth',2)

		
	end

	Z = X +i*Y;
	
	B(1) = 1;
	B(2) = -2*real(Z(1)+Z(2));
	B(3) = (abs(Z(1))).^2 + (abs(Z(2))).^2 +4*real(Z(1))*real(Z(2));
	B(4) = -2*((abs(Z(1))).^2*real(Z(2))+(abs(Z(2))).^2*real(Z(1)));
	B(5) = (abs(Z(1)*Z(2))).^2;
	
end

hold off