function [y,E,r,s,bE] = polyfilt(d,m,x,t,tt,alpha,fig)
%               FILTERING USING STEPWISE POLYNOMIAL REGRESSION
% for the data model x(t) = y(t) + e with random disturbances e: N(0,sigma^2*I)
% -----------------------------------------------------------------------------
% CALLING:
% 1) SMOOTHING PROCESS : [y  ,E ,r,s,bE] = polyfilt(d,m,x ,t,tt,alpha,fig)
% 2) SMOOTHING OPERATOR: [FIR,E0]        = polyfilt(d,m,[],t,tt,alpha)
% INPUT:
% d  ... degree of the moving polynomial, d>=0 (integer)
% m  ... length of the filter impulse response (FIR-length), m>d (integer)
% x  ... is an n x 1 vector of data observed at (time) t (optional: default=[])
%        x==[] => FIR and E0 are matrices defining the smoothing operator
% t  ... is an n x 1 vector of time mesh points for observations x
%        (optional: default=[]=uniform mesh (1:n)')
%        n>=m ... filter length cannot exceed data size
%        x==[] & t==[] => n=m (t are time coordinates of filter coefficients)
% 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:
% 1) SMOOTHING PROCESS:
% y ...  is an N x 1 vector of smoothed (filtered) data
% E ...  [Es,Ep] is an N x 2 matrix, where Ep>Es and
%        [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
% 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 ... =[b1,E1,b2,E2,...] is an (d+1) x 2*(n-m+1) matrix, where
%         b# are coefficient estimates of #-th regression polynomial
%         p#(t)=b#(1)*(t-ts)^d+...+b#(d)*(t-ts)+b#(d+1) where
%         ts is the center of the #-th regression interval (window), and
%         [b#-E#,b#+E#] are their confidence intervals of (100-alpha) %
%
% 2) SMOOTHING OPERATOR:
% FIR... N x n matrix defining the smoothing operator y = FIR*x
%        FIR(k,:) ... FIR for computing y(k) from x, k=1,2,...,N
% E0...  [E0s,E0p] is an N x 2 matrix determining the confidence matrix E
%        explained above by E=sigma.*E0 where sigma is an estimate of the
%        standard deviation of random disturbances,
%        usually sigma=norm(x-y)/sqrt(n-d-1) where y are evaluated at
%        tt=t to obtain the residuals r=x-y.

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

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

p = d + 1; q = fix(m/2);
if m<p, error('Filter length is not greater than the degree ''d'''); end

n = max(length(t),length(x)); n = n+m*(n==0);
if n<m, error('Filter length ''m'' exceeds the data size'); end

equispaced = isempty(t);
if equispaced, t = (1:n).';
else, t = reshape(t,length(t),1);
end
[t,I] = sort(t);
if ~xempty, x = x(I); end
clear I;

if isempty(tt), tt= t;
else, tt = reshape(tt,length(tt),1); end
N = length(tt);
tt = sort(tt);

convfilter = 0;
if N==n
   convfilter = (all(t==tt) & equispaced);
end

if ~xempty & length(x)~=length(t)
   error('The lengths of observed data ''x'' and time vector ''t'' differ ');
end

D(:,p) = ones(m,1);
if m==n                                 % Standard polynomial regression
   t0 = (t(1)+t(m))./2;
   C(:,p) = ones(N,1);
   for k=d:-1:1
       D(:,k) = (t -t0).*D(:,k+1);
       C(:,k) = (tt-t0).*C(:,k+1);
   end
   [bE,rE,sE,PbE,CbE] = linregr(x,D,C,[],alpha);
   if xempty
      y = C*bE;
      E = [CbE,PbE(1:N)];
   else
      y = CbE(:,1);
      E = [CbE(:,2),PbE(1:N,2)];
      r = rE(:,1);
      s = sE(1);
   end
else                                    % Stepwise polynomial regression
   if ~xempty & convfilter              % -using fast convolution filtering
      [FIR,E] = polyfilt(d,m,[],[],[],alpha);
      E = [E(1:q,:);ones(n-m+1,1)*E(q+1,:);E(q+2:m,:)];
      y = convfilt(FIR,x);
      r = x - y;
   else                                 % -using direct approach
      y = zeros(N,1+xempty*(n-1)); r = zeros(n,1); E = zeros(N,2);
      if nargout>4, bE = zeros(p,2*(n-m+1)); end
      a = 1;
      for s=m:n
          sm = s - m; sm1 = sm + 1;
          ts = t(sm1:s);
          tC = (s<n)*(ts(q+1)+ts(m-q+1))./2 + (s==n)*tt(N);
          t0 = (ts(1)+ts(m))./2;
          ts = ts - t0; xs = x(sm1:(~xempty)*s);
          b  = a + sum(tt(a:N)<=tC) - 1;
          if b>=a | (~xempty  & nargout>1)
             sD = 1+(m<s)*q : m-(s<n)*(m-q-1); % Indices of residue segment
             sC = a:b;   % compute index range in tt for evaluation at step s
             a  = b+1;
             C  = [zeros(length(sC),p-1),ones(length(sC),1)];
             for k=d:-1:1
                 D(:,k) = ts.*D(:,k+1);
                 if length(sC)>0, C(:,k) = (tt(sC)-t0).*C(:,k+1); end
             end
             [b,rE,sE,PbE,CbE] = linregr(xs,D,C,[],alpha);
             if xempty
                y(sC,sm1:s) = C*b;
                E(sC,:)     = [CbE(:,1),PbE(1:length(sC),1)];
             else
                r(sD+sm) = rE(sD,1);
                if length(sC)>0
                   y(sC)  = CbE(:,1);
                   E(sC,:)=[CbE(:,2),PbE(1:length(sC),2)]./sE(1);
                end
                if nargout>4
                   bE(:,sm+sm1:sm1+sm1) = b;
                end
             end
          end
      end
   end
   if ~xempty, s = norm(r)./sqrt(n-p); E = E*s; end
end
if ~xempty & ~isempty(fig)    % plot regression curves and confidence bands
   plot(t,x,'go',tt,y,'w',...
        tt,y-E(:,1),'y--',tt,y+E(:,1),'y--',tt,y-E(:,2),'y:',tt,y+E(:,2),'y:');
end
