PPF Piece-wise Polynomial Fit Y_FIT=PPF(X,Y,X0,N,M); fits Y=f(X) to a N-th order polynomial fit with X0 knots or control points. X is a vector and Y can be a vector or a matrix (colums are observations), X0 can be a scalar (number of control points) or a vector (control points), N is the order of the polynomial fit, and M is the order of continuity imposed (one for continuity of first-derivatives, two for continuity in the second-derivatives, etc.). N defaults to 3, M defaults to N-1, and X0 defaults to X [Y_FIT,PARAMS]=PPF(...) also returns the parameters of the fit Y_NEW=PPF(X_NEW,PARAMS) estimates new data using previous fit % example 1: Fitting sequential points in a curve x=2:9; y=rand(size(x)); [y_fit,params]=ppf(x,y); x_new=linspace(1,10,1000); y_new=ppf(x_new,params); plot(x,y,'bo',x_new,y_new,'r-'); set(gca,'xtick',params.control); grid on; % example 2: Approximating data points to a polynomial curve x=sort(rand(1,200)); y=x.*(sin(6*pi*x)+.5*randn(size(x))); [y_fit,params]=ppf(x,y,10); % 10 control points plot(x,y,'bo',[x;x],[y;y_fit],'b:',x,y_fit,'r.-'); set(gca,'xtick',params.control); grid on;
0001 function [y_fit,params]=ppf(x,y,control,order,continuity,weight,type); 0002 % PPF Piece-wise Polynomial Fit 0003 % 0004 % Y_FIT=PPF(X,Y,X0,N,M); 0005 % fits Y=f(X) to a N-th order polynomial fit with X0 knots or control points. 0006 % X is a vector and Y can be a vector or a matrix (colums are observations), 0007 % X0 can be a scalar (number of control points) or a vector (control points), 0008 % N is the order of the polynomial fit, and M is the order of continuity imposed 0009 % (one for continuity of first-derivatives, two for continuity in the second-derivatives, etc.). 0010 % N defaults to 3, M defaults to N-1, and X0 defaults to X 0011 % 0012 % [Y_FIT,PARAMS]=PPF(...) also returns the parameters of the fit 0013 % Y_NEW=PPF(X_NEW,PARAMS) estimates new data using previous fit 0014 % 0015 % % example 1: Fitting sequential points in a curve 0016 % x=2:9; y=rand(size(x)); 0017 % [y_fit,params]=ppf(x,y); 0018 % x_new=linspace(1,10,1000); 0019 % y_new=ppf(x_new,params); 0020 % plot(x,y,'bo',x_new,y_new,'r-'); set(gca,'xtick',params.control); grid on; 0021 % 0022 % % example 2: Approximating data points to a polynomial curve 0023 % x=sort(rand(1,200)); y=x.*(sin(6*pi*x)+.5*randn(size(x))); 0024 % [y_fit,params]=ppf(x,y,10); % 10 control points 0025 % plot(x,y,'bo',[x;x],[y;y_fit],'b:',x,y_fit,'r.-'); set(gca,'xtick',params.control); grid on; 0026 % 0027 0028 % alfnie@bu.edu 0029 % 05/06 0030 0031 x=x(:)'; sX=size(x); 0032 if isstruct(y), 0033 params=y; 0034 ncontrol=length(params.control); 0035 if ~params.type, 0036 x=ppf_normdata(sX(2),params.control,ncontrol,x); 0037 X=ppf_createX(sX(2),0:ncontrol-1,params.order,params.continuity,params.weight,ncontrol,x); 0038 else, 0039 X=ppf_createX(sX(2),params.control,params.order,params.continuity,params.weight,ncontrol,x); 0040 end 0041 y_fit=params.b'*X'; 0042 else, 0043 if sum(size(y)>1)<=1, y=y(:)'; end 0044 if nargin<3 || isempty(control), control=unique(x); end 0045 if nargin<4 || isempty(order), order=3; end 0046 if nargin<5 || isempty(continuity), continuity=order-1; end 0047 if nargin<6 || isempty(weight), weight=.01; end % when underdetermined, (weight<1) weights more highest order terms 0048 if nargin<7 || isempty(type), type=1; end % (type=0) equalizes x data prior to polynomial fitting 0049 sY=size(y); 0050 0051 if prod(size(control))==1, 0052 sx=sort(x); 0053 ncontrol=control; 0054 control=sx(round(linspace(1,sX(2),ncontrol))); 0055 else, ncontrol=length(control); end 0056 if ~type, 0057 x=ppf_normdata(sX(2),control,ncontrol,x); 0058 X=ppf_createX(sX(2),0:ncontrol-1,order,continuity,weight,ncontrol,x); 0059 else, 0060 X=ppf_createX(sX(2),control,order,continuity,weight,ncontrol,x); 0061 end 0062 b=pinv(1e-10*eye(size(X,2))+X'*X)*X'*y'; 0063 y_fit=b'*X'; 0064 if nargout>1, 0065 params.order=order; 0066 params.continuity=continuity; 0067 params.control=control; 0068 params.b=b; 0069 params.weight=weight; 0070 params.type=type; 0071 %params.X=X; 0072 end 0073 end 0074 0075 function X=ppf_createX(sX,control,order,continuity,weight,ncontrol,x); 0076 X=ones(sX,order+(order-continuity)*(ncontrol-1)); for n2=1:order-1, X(:,1+n2)=(x').^n2; end; 0077 for n2=1:order-continuity, 0078 t1=(max(0,x'-control(1)).^(eps+order-n2+1)); 0079 for n1=1:ncontrol-1, 0080 t2=(max(0,x'-control(n1+1)).^(eps+order-n2+1)); 0081 X(:,order+(order-continuity)*(n1-1)+n2)=(weight)*(t1-t2); 0082 t1=t2; 0083 end 0084 end 0085 0086 function x=ppf_normdata(sX,control,ncontrol,x) 0087 r=sum(repmat(control',[1,sX])<repmat(x,[ncontrol,1]),1); 0088 idx=find(~r); x(idx)=0; 0089 idx=find(r==ncontrol); x(idx)=ncontrol-1; 0090 idx=find(r>0 & r<ncontrol); x(idx)=r(idx)-1+(x(idx)-control(r(idx)))./(control(r(idx)+1)-control(r(idx)));