cspline1slow.m

模式识别工具包

%	cspline1slow	- Cubic smoothing spline 1-D (unaccelerated)
%
%	[fx]	= cspline1slow(x,y,p)
%
%	_____OUTPUT_____________________________________________________________
%	fx	smoothed y corresponding to f(x)		(col vectors)
%
%	_____INPUT______________________________________________________________
%	x	independent scalars				(col vector)
%	y	dependent vector				(col vectors)
%		(independently smoothed in each row/dim)
%	p	smoothing parameter				(scalar)
%		[0,1] p=1 implies point interpolation
%
%	_____NOTES______________________________________________________________
%	Assumptions
%	1 {x,y} are sorted according to {x}
%	2 {x} are unique (no repetitions!)
%		
%	To do:
%
%	implement short-cut/speed-ups
%		Qt*Q	(figure out!!!!)
%
%	implement weighted smoothing via matrix W
%	implement end point conditions (start and end equal)
%
%	_____SEE ALSO___________________________________________________________
%
%	(C) 1998.09.02 Kui-yu Chang
%	http://lans.ece.utexas.edu/~kuiyu

%	This program is free software; you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation; either version 2 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program; if not, write to the Free Software
%	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
%	or check
%			http://www.gnu.org/


function	[fx]	= cspline1slow(x,y,p)

[n d]	= size(y);

%-----	Find duplicate x and weigh accordingly
%w	= [1 1 10 1 1];		% more weight => more important, less variation
%w	= w.*w;
%w	= ones(n,1)

%-----	Common terms
h	= diff(x);		% 1 ... n-1
hi	= 1./h;

%-----	form Qt	(n-2,n)		optimized 1998 Sep 01
dleft			= zeros(n-2,n);
dleft(1:n-2,1:n-2)	= diag(hi(1:n-2));
dright			= zeros(n-2,n);
dright(1:n-2,3:n)	= diag(hi(2:n-1));
dcent			= zeros(n-2,n);
dcent(:,2:n-1)		= diag(-hi(2:n-1)-hi(1:n-2));
Qt	= dleft+dright+dcent;
Q	= Qt';

%-----	form R	(n-2,n-2)	optimized 1998 Sep 01
dleft	= diag([h(2:n-2)],-1);
dcent	= diag(2*(h(2:n-1)+h(1:n-2)));
dright	= diag([h(2:n-2)],1);
R	= dleft+dcent+dright;

%-----	form K
%W	= diag(w);
K	= Q*inv(R)*Q';
fx	= (eye(size(K))+p*K)\y;
return