function P = obsv(A,C,p)
% P = obsv(A,C,p)
% Form the observability matrix P = [ C ; CA ; CA^2 ; CA^3 ; ... ; C A^(p-1) ]
% for the discrete time system, x(k+1) = A*x(k); y(k) = C*x(k) 

 [m,n] = size(C);                               % number of outputs and states
 if nargin < 3, p = n; else p = max(p,n); end   % number of blocks in P
 P = zeros(m*p,n);				% dimension of P
 P(1:m,:) = C;                                  % first block is output matrix
 for k=1:p-1
	CAk = P((k-1)*m+1:(k-1)*m+m,:);
	P(k*m+1:k*m+m,:) = CAk * A;
 end
% ---------------------------------------------------------------------- OBSV
