function [yE,trE,pE,r,s,bE,cE,sE] = polfs(x,d,Tk,t,tt,alpha,fig)
%         ESTIMATING POLYNOMIAL TREND tr(t) AND PERIODIC COMPONENT p(t)
% using linear regression for the data model
% x(t) = y(t) + e with random disturbances e: N(0,sigma^2*I)
% where y(t)=tr(t)+p(t), tr(t)=polynomial, p(t)=goniometric polynomial
% -----------------------------------------------------------------------------
% CALLING:
% [yE,trE,pE,r,s,bE,cE,sE] = polfs(x,d,Tk,t,tt,alpha,fig)
% INPUT:
% x  ... is an n x 1 vector of data observed at (time) t
% d  ... degree of the polynomial tr(t), d>=0 (integer)
% Tk ... is an 1 x q vector of periods of significant harmonics
%        (see kF or kS returned by 'tperiod' - typically Tk=n./kF or Tk=n./kS)
%        (optional: default=[]=no periodic component)
% t  ... is an n x 1 vector of time mesh points for observations x
%        (optional: default=[]=uniform mesh (1:n)')
% tt ... N x 1 vector of time mesh for smoothed values y
%        (optional: default=[]=t)
% alpha ... is a percentual confidence risk (optional: default=[]=5%)
% fig   ... figure handle for the regression plot (optional: default=[]=no plot)
%
% OUTPUT:
% yE ... =[y,Es,Ep] is an N x 3 matrix where
%        y are smoothed values y(tt)=tr(tt)+p(tt)
%        [y-Es,y+Es] are confidence intervals of (100-alpha) %
%                    for smoothed values y(tt)
%        [y-Ep,y+Ep] are confidence intervals of (100-alpha) %
%                    for predicted values y(tt)+e
% trE... =[tr,Es,Ep] is an N x 3 matrix describing trend component tr(tt)
%         analogously
% pE ... =[p, Es,Ep] is an N x 3 matrix describing periodic component p(tt)
%         analogously
% r  ...  n x 1 vector of residulas x-y where y are evaluated at tt=t,
%         r is an estimate of the random error e.
% s  ...  the estimate of the standard deviation sigma (scalar)
% bE ... =[b,E] is an (d+1) x 2 matrix, where
%         b are coefficient estimates of the regression polynomial
%         tr(t)=b(1)*(t-ts)^d+...+b(d)*(t-ts)+b(d+1) where ts=(t(1)+t(n))/2
%         [b-E,b+E] are their confidence intervals of (100-alpha) %
% cE ... =[c,E] is an q x 2 matrix and
% sE ... =[s,E] is an q x 2 matrix where c and s are coefficient estimates
%         of the the goniometric polynomial
%         p(t) =c(1)*cos(2*pi*k(1)*(t-ts)/n)+...+c(q)*cos(2*pi*k(q)*(t-ts)/n) +
%               s(1)*sin(2*pi*k(1)*(t-ts)/n)+...+s(q)*sin(2*pi*k(q)*(t-ts)/n)
%         [c-E,c+E] and [s-E,s+E] are their confidence intervals of (100-alpha)%

% (C) 1997, V. Vesely, Masaryk University of Brno (Czech Republic)
%

error(nargchk(2,7,nargin));
if nargin<3, Tk= []; end
if nargin<4, t = []; end
if nargin<5, tt= []; end
if nargin<6, alpha=[]; end
if nargin<7, fig  =[]; end

p = d + 1;
q = length(Tk);
kh= reshape(1./Tk,1,q);

n = length(x);
if isempty(t), t=(1:n).'; else t=reshape(t,length(t),1); end
if n~=length(t), error('x and t must be of equal size'); end
if n<p+2*q, error('x is too short: model is overparametrized'); end
if isempty(tt), tt=t; else tt = reshape(tt,length(tt),1); end
N = length(tt);
ts=(t(1)+t(n))./2;

D(:,p) = ones(n,1);
C      = zeros(N,p+2*q);
C(:,p) = ones(N,1);
for k=d:-1:1
    D(:,k) = (t -ts).*D(:,k+1);
    C(:,k) = (tt-ts).*C(:,k+1);
end

if q==0
   [bE,rE,sEE,PbE,CbE] = linregr(x,D,C,[],alpha);
   yE  = [CbE(1:N,:),PbE(1:N,2)];
   trE = yE;
   pE  = zeros(N,1);
   sE  = zeros(0,2);
   cE  = zeros(0,2);
else
   C2 = 2*pi*(t-ts)*kh;
   D  = [D,cos(C2),sin(C2)];
   C2 = 2*pi*(tt-ts)*kh;
   C2 = [zeros(N,p),cos(C2),sin(C2)];
   C1 = C;
   C  = [C1+C2;C1;C2];
   [bE,rE,sEE,PbE,CbE] = linregr(x,D,C,[],alpha);
   yE  = [CbE(1:N,:),PbE(1:N,2)];
   trE = [CbE(N+1:2*N,:)  ,PbE(N+1:2*N,2)];
   pE  = [CbE(2*N+1:3*N,:),PbE(2*N+1:3*N,2)];
   cE  = bE(p+1:p+q,:);
   sE  = bE(p+q+1:p+2*q,:);
   bE  = bE(1:p,:);
end
r   = rE(:,1);
s   = sEE(1);

if ~isempty(fig)    % plot regression curves and confidence bands
   y = yE(:,1);
   plot(t,x,'co',tt,y,'w',...
     tt,y-yE(:,2),'r--',tt,y+yE(:,2),'r--',tt,y-yE(:,3),'g:',tt,y+yE(:,3),'g:');
   title('Smoothed curve');
   pos0 = get(fig,'Position'); pos=pos0;
   if q>0
      pause(2); figure; pos(1:2)=[1.2,0.8].*pos0(1:2); set(gcf,'Position',pos);
      plot(t,x,'co',tt,trE(:,1),'w',...
           tt,trE(:,1)-trE(:,2),'r--',tt,trE(:,1)+trE(:,2),'r--');  % ,...
%           tt,trE(:,1)-trE(:,3),'g:', tt,trE(:,1)+trE(:,3),'g:');
      title('Trend component');
      pause(2); figure; pos(1:2)=[1.4,0.6].*pos0(1:2); set(gcf,'Position',pos);
      plot(tt,pE(:,1),'w',...
           tt,pE(:,1)-pE(:,2),'r--',tt,pE(:,1)+pE(:,2),'r--');      % ,...
%           tt,pE(:,1)-pE(:,3),'g:',tt, pE(:,1)+pE(:,3),'g:');
      title('Periodic component');
   end
   pause(2); figure; pos(1:2)=[1.6,0.4].*pos0(1:2); set(gcf,'Position',pos);
   plot(t,rE(:,1),'y');
   title('Error component');
end
