(* :Title: Smooth *)

(* :Context: "Smooth`" *)

(* :Author:
	P.H. Lundow
	Bug reports to phl@kth.se
*)
(* :Summary:
	Functions for smoothing equidistant data with a Savitsky-Golay filter.
*)
(* :History:
	020507 Created.
	030417 Smooth extended to lists of y-values.
	030425 Added SmoothAll.
	030619 Added MinimumDistance and EquidistantData.
	030923 Small improvement in SGKernel.
	070618 Clean-up.
*)
(* :Mathematica Version: 4.0 *)
(* :Keywords: *) 
(* :Limitations: *) 
(* :Discussion:
	There is room for improvement here. Returned list of smoothed
	data is shorter than original list. Find a better way to
	produce many kernels that takes care of the margins.
*)


BeginPackage["Smooth`"]

SGKernel::usage =
	"SGKernel[left, right, degree, derivative] returns the Savitsky-Golay
	kernel used by the function Smooth when convoluting the data.
	The kernel has length left+right+1. Left is the number of
	leftward data points and right is the number of rightward
	points. Degree refers to the degree of the polynomial and
	derivative is the order of the desired derivative.\nRef:
	Numerical Recipes, chapter 14.8."

Smooth::usage = 
       	"Smooth[list, window, degree, derivative] returns the smoothed form
	of equally spaced and sorted data. This is done by fitting polynomials
	of a given degree to a moving window of the list. Argument list is
	either a list of pairs {{x1,y1},{x2,y2},...} or a list {y1,y2,...}
	where the x-values are taken to be 1,2,3... If a derivative of the
	data is desired then give the order of the derivative, default is 0.
	\nExample:\n a=Table[{x,Sin[x]+Random[Real,{-1,1}],{x,0,2*Pi,0.001}];
	\n b=Smooth[a,15,2];
	\n c=Smooth[a,15,2,1];
	\n This fits 2:nd degree polynomials to moving windows of width 15 to
	the data. List b is the smoothed data and list c is the derivative of
	the data."

SmoothAll::usage =
	"SmoothAll[{vector,matrix}, window, degree, derivative] works like
	Smooth except that vector is a list of x-values and matrix is a list
	of lists of y-values at the corresponding x-values. The result is
	returned on the same form.\nExample:\n xold={1,2,3,4,5};
	\n yold={{1,3,5,4,4},{2,3,3,2,1},{3,4,6,4,3}};
	\n {xnew,ynew}=SmoothAll[{xold,yold},2,1,0];"

MinimumDistance::usage = 
	"MinimumDistance[data] returns the minimum distance between
	two x-values in the data which is a list of lists. The
	x-values are the first positions of each sublist."

EquidistantData::usage =
	"EquidistantData[data] extracts the largest chunk of
	equidistant data from a list of lists.\nExample:\n
	EquidistantData[{{0,1},{2,8},{3,7},{4,9},{6,3}}] returns the
	list\n {{2,8},{3,7},{4,9}}."

Begin["`Private`"]

definitions = {MinimumDistance, EquidistantData, SGKernel, Smooth, SmoothAll}

Scan[Unprotect, definitions]


MinimumDistance[data:{__List}] :=
	Module[{x},
		x = Map[First, data];
		Min[Drop[x, 1] - Drop[x, -1]]
	]


EquidistantData[data:{__List}] :=
	Module[{min, len, pos, tmp},
		min = MinimumDistance[data];
		tmp = Split[data, (#1[[1]] + min == #2[[1]])&];
		len = Map[Length, tmp];
		pos = Flatten[Position[len, Max[len], Infinity, 1]][[1]];
		tmp[[pos]]
	]


SGKernel[left_?NonNegative, right_?NonNegative, degree_?NonNegative, 
	derivative_?NonNegative] :=
	Module[{i, j, k, l, matrix, vector},
		matrix = Table[  (* matrix is symmetric *)
			l = i + j;
			If[l == 0,
				left + right + 1,
			(*Else*)
				Sum[k^l, {k, -left, right}]
			],
			{i, 0, degree},
			{j, 0, degree}
		];
		vector = LinearSolve[
			matrix, 
			MapAt[1&, Table[0, {degree+1}], derivative+1]
		];
		(* vector = Inverse[matrix][[derivative + 1]]; *)
		Table[
			vector.Table[If[i == 0, 1, k^i], {i, 0, degree}],
			{k, -left, right}
		]
       	] /; derivative <= degree <= left+right
	

Smooth[list_?MatrixQ, window_, degree_, derivative_:0]:=
	Module[{kernel, list1, list2, margin, space},
		margin = Floor[window/2];
		list1 = Take[
			Map[First, list], 
			{margin + 1, Length[list] - margin}
		];
		list2 = Map[Last, list];
		kernel = SGKernel[margin, margin, degree, derivative];
		list2 = ListCorrelate[kernel, list2];
		(* Data _should_ be equally spaced, but... *)
		space = Min[Drop[list1, 1] - Drop[list1, -1]];
		list2 = list2*(derivative!/space^derivative);
		Transpose[{list1, list2}]
	] /; derivative <= degree <= 2*Floor[window/2] && $VersionNumber >= 4.0

Smooth[list_?VectorQ, window_, degree_, derivative_:0]:=
	Module[{pairs},
		pairs = MapThread[List, {Range[Length[list]], list}];
		Map[Last, Smooth[pairs, window, degree, derivative]]
	]


SmoothAll[{x_?VectorQ, y_?MatrixQ}, window_, degree_, derivative_:0]:=
	Module[{old, new, tmp},
		tmp = Map[(
			old = Transpose[{x, #}];
			new = Smooth[old, window, degree, derivative];
			Map[Last, new])&,
			y
		];
		{Map[First, new], tmp}
	]


Scan[Protect, definitions]


End[]

EndPackage[]