Home > fcn > ppf.m

ppf

PURPOSE ^

PPF Piece-wise Polynomial Fit

SYNOPSIS ^

function [y_fit,params]=ppf(x,y,control,order,continuity,weight,type);

DESCRIPTION ^

 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;

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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)));

Generated on Tue 27-Mar-2007 12:06:24 by m2html © 2003