function [y,yi_or_p] = orthofit(x,y,n,xi)

%ORTHOFIT Fit polynomial to data.
%   YS = ORTHOFIT(X,Y,N) smooths/fits data Y(X) in a least-squares sense
%   using a polynomial of degree N and returns the smoothed data YS.
%
%   [YS,YI] = ORTHOFIT(X,Y,N,XI) also returns the values YI of the fitting
%   polynomial at the points of a different XI array.
%
%   YI = ORTHOFIT(X,Y,N,XI) returns only the values YI of the fitting
%   polynomial at the points of the XI array.
%
%   [YS,P] = ORTHOFIT(X,Y,N) returns the polynomial coefficients P for use
%   with POLYVAL.
%
%   Notes:
%   -----
%   ORTHOFIT smooths/fits data using a polynomial of degree N developed in
%   a sequence of orthogonal polynomials. ORTHOFIT is more appropriate than
%   POLYFIT for polynomial fitting and smoothing since this method does not
%   involve any matrix linear system but a simple recursive procedure.
%   Degrees much higher than 30 could be used with orthogonal polynomials,
%   whereas badly conditioned matrices may appear with a classical
%   polynomial fitting of degree typically higher than 10.
%
%   To avoid using unnecessarily high degrees, you may let the function
%   POLYDEG choose it for you. POLYDEG finds an optimal polynomial degree
%   according to the Akaike's information criterion (available <a
%   href="matlab:web('http://www.biomecardio.com/matlab/polydeg.html')">here</a>).
%
%   Example:
%   -------
%   x = linspace(0,10,300);
%   y = sin(x.^3/100).^2 + 0.05*randn(size(x));
%   ys = orthofit(x,y,25);
%   plot(x,y,'.',x,ys,'k')
%   % try POLYFIT for comparison...
%
%   % Automatic degree determination with <a
%   href="matlab:web('http://www.biomecardio.com/matlab/polydeg.html')">POLYDEG</a>
%   n = polydeg(x,y);
%   ys = orthofit(x,y,n);
%   plot(x,y,'.',x,ys,'k')
%
%   Reference: Méthodes de calcul numérique 2. JP Nougier. Hermes Science
%   Publications, 2001. Section 4.7 pp 116-121
%
%   Damien Garcia, 09/2007, revised 01/2010
%
%   See also POLYDEG, POLYFIT, POLYVAL.


%% Check input arguments
error(nargchk(3,4,nargin));

if ~isequal(size(x),size(y))
    error('MATLAB:orthofit:XYSizeMismatch',...
          'X and Y must have the same size.')
elseif ~isscalar(n) || n~=round(abs(n))
    error('MATLAB:orthofit:WrongNValue',...
        'N must be a positive integer.')
elseif ~isfloat(x) || ~isfloat(y)
    error('MATLAB:orthofit:NonFloatingXY',...
    'X and Y must be floating point arrays.')
elseif nargin==4 && ~isfloat(xi)
    error('MATLAB:orthofit:NonFloatingXI',...
    'XI must be a floating point array.')
end

%% Particular case: n=0
if n==0
    p = mean(y(:));
    y = ones(size(y))*p;
    if nargin<4
        yi_or_p = p;
    elseif nargout==2
        yi_or_p = ones(size(xi))*p;
    else
        y = ones(size(xi))*p;
    end
    return
end

%% Reshape
x = x(:);
siz0 = size(y);
y = y(:);

%% Coefficients of the orthogonal polynomials
p = zeros(3,n+1);
p(2,2) = mean(x);

N = length(x);
PL = ones(N,n+1);
PL(:,2) = x-p(2,2);

for i = 3:n+1
    p(2,i) = sum(x.*PL(:,i-1).^2)/sum(PL(:,i-1).^2);
    p(3,i) = sum(x.*PL(:,i-2).*PL(:,i-1))/sum(PL(:,i-2).^2);
    PL(:,i) = (x-p(2,i)).*PL(:,i-1)-p(3,i)*PL(:,i-2);
end
p(1,:) = sum(bsxfun(@times,PL,y))./sum(PL.^2);

%% ys = smoothed y
if nargin<4 || (nargin==4 && nargout==2)
    y = sum(bsxfun(@times,PL,p(1,:)),2);
    y = reshape(y,siz0);
end

%% Coefficients of the polynomial in its final form
if nargin<4
    yi = zeros(n+1);
    yi(1,n+1) = 1;
    yi(2,n:n+1) = [1 -p(2,2)];
    for i = 3:n+1
        yi(i,:) = [yi(i-1,2:end) 0]-p(2,i).*yi(i-1,:)-p(3,i)*yi(i-2,:);
    end
    p = sum(bsxfun(@times,p(1,:).',yi),1);
    yi_or_p = p;
    return
end

%% yi values for a given xi
siz0 = size(xi);
xi = xi(:);

N = length(xi);
n = size(p,2);
yi = ones(N,n);

yi(:,2) = xi-p(2,2);
for i = 3:n
    yi(:,i) = (xi-p(2,i)).*yi(:,i-1) - bsxfun(@times,p(3,i),yi(:,i-2));
end
yi = sum(bsxfun(@times,p(1,:),yi),2);
yi = reshape(yi,siz0);

if nargout==2
    yi_or_p = yi;
else
    y = yi;
end