function y = localreg(data,h,order)

% LOCALREG   Local-regression smoother.
%
%  Y = LOCALREG(DATA,H)
%  Y = LOCALREG(DATA,H,ORDER)
%
%  DATA is the time series to be smoothed. When DATA is a matrix 
%  instead of vector, each column is smoothed separately.
%  The smoothed data, Y, has the same shape as DATA.
%
%  H is a positive integer specifying the bandwidth. 
%  Each local regression is performed within a window of width 2H+1.
%  ORDER is 'linear', 'quadratic', or 'cubic'. When not supplied, 
%  ORDER defaults to 'quadratic'. (Note: The quadratic and cubic 
%  smoothers are identical everywhere except near the edges.)
%
%  Algorithm: A window of width 2H+1 is moved along the time series,
%  centered on each data point in turn. A regression polynomial is
%  fitted within each window and used to smooth the value at the 
%  center. The points near the edges of the time series are predicted 
%  from asymmetric windows of the same width 2H+1. No points are
%  lost near the edges and thus the smoothed series has the same 
%  length as the original series.
%
%  Example:
%    t = [1:100]' ;
%    signal = 5 * sin(t/5) ;
%    noisy_data = signal + 5 * randn(100,1) ;
%    smooth = localreg(noisy_data,10) ;
%    twice = localreg(smooth,10) ;
%    plot(t,signal,'b:',t,noisy_data,'b.',t,smooth,'r',t,twice,'k:') ;
%
%  See also MOVAVG, FILTER.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2000/01/04 18:20 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

%--  Decipher input arguments
if (nargin < 3)
  order = 'quadratic' ;
end

[rows,cols] = size(data) ;
if (rows==1 & cols>1)
  y = localreg(data',h,order)' ;   % row vector --> column vector
  return
end

%--  Prepare predictors for the regressions
winwidth = 2*h + 1 ;
t = [-h : +h]' ;     % zero at the center
switch lower(order)
  case 'linear'
    X = [ones(winwidth,1) t] ;
  case 'quadratic'
    X = [ones(winwidth,1) t t.^2] ;
  case 'cubic'
    X = [ones(winwidth,1) t t.^2 t.^3] ;
end    
left_X = X(1:h,:) ;
right_X = X(h+2:winwidth,:) ;


%--  Prepare a shift-index matrix for the moving windows
if (winwidth > rows)
  error('Window width too large for this number of data points.') ;
end
S = shiftmat(winwidth,rows-winwidth+1) ;

%--  Main loop
y = zeros(rows,cols) ;
for k = 1:cols
  curr_data = data(:,k) ;
  %- Do all the regressions at once
  replicated_data = curr_data(S) ;
  reg_coefs = X \ replicated_data ;
  %- Predict the leftmost H points from the leftmost regression
  left_y = left_X * reg_coefs(:,1) ;
  %- Predict the middle points, each from its own regression
  %  All terms but the intercept vanish at the central point (T=0)
  middle_y = reg_coefs(1,:)' ;
  %- Predict the rightmost H points from the rightmost regression
  right_y = right_X * reg_coefs(:,end) ;
  %- Put all parts together
  y(:,k) = [left_y ; middle_y ; right_y] ;
end  

%-- Return Y
%%%%% End of LOCALREG.M
