function [phi,del]=peryule(x,T,p,norm)
%
%		[phi,del]=peryule(x,T,p,norm)
%
%	x = univariate time series
%	T = period of PC structure
%  p = order of autoregression; assumed constant over time
%  norm = normalization flag; if norm > 0, normalize by sample periodic std's
%  phi = T x p matrix of estimated periodic PAR coefficients
%  del = T x 1 vector of noise weights
%
% NOTES: it is assumed that periodic means have been removed
%			it is required that p <= length(x)
%			the basic idea appears in Pagano, but this is mainly from Vecchia
%			but with our own notation;
%
%	the Periodic Yule-Walker equations are solved for each t in the period
%  NOTE : an error occurs when p>T due to code around NOTE(A); indexes are wrong
%  so do not use for p > T
nx=length(x);
if norm
	stdx=zeros(T,1);
    for t=1:T
        index=[t:T:nx];				% this will get them all, including remainder
        stdx(t)=std(x(index));
        x(index)=x(index)/stdx(t);	% mean not removed as assumed zero
   end
end


r=zeros(p,1);								% col vector
phi=zeros(T,p);
del=zeros(T,1);						% col vector
for t=1:T
   baseind=[t:-1:t-p];					% get the base index set for this t
   if t-p <= 0								% smallest one is neg
      baseind=baseind+T;				% NOTE(A) put in next period assuming one T will do it
   end										% all must be increased by T together
   inum=floor((nx-baseind)/T)+1;		% number of periods in nx when you start at baseind
   inummin=min(inum);
   yt=x([baseind(1):T:nx])';			% sample x starting at t; make it col
   yt=yt(1:inummin);						% chop it off if too long
   y=zeros(inummin,p);
   for s=1:p
      isamp=[baseind(s+1):T:nx];		% baseind(s)+1 is t-1 
      isamp=isamp(1:inummin);			% chop if too long
      y(:,s)=x(isamp)';					% sample x at isamp; make a col
      r(s)=yt'*y(:,s)/(inummin-1);	% get the sth element of r(p,1)	
   end
   R=cov(y);
   phi(t,:)=(inv(R)*r)';				% make it a row
   del(t)=var(yt)-phi(t,:)*r;	        % recall r is a col
end