(* :Name: NumericalMath`TrigFit` *)

(* :Title: Curve Fitting over Trigonometric Series *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
This package uses the function Fourier to fit trigonometric
series to data that is assumed to be periodic.  The list of
data points is suppossed to cover exactly one period.
*)

(* :Context: NumericalMath`TrigFit` *)

(* :Package Version: 2.0 *)

(* :Copyright: Copyright 1990-1999, Wolfram Research, Inc.
*)

(* :History:
	Originally written by Jerry B. Keiper, March 1994.
*)

(* :Keywords: Fit, least squared error, trigonometric approximation *)

(* :Source:
	Exploring Mathematics with Mathematica, by Theodore W. Gray and
	Jerry Glynn.  Addison-Wesley, 1991
*)

(* :Warnings: None. *)

(* :Mathematica Version: 2.0 *)

(* :Limitations:
*)


BeginPackage["NumericalMath`TrigFit`"]

TrigFit::usage =
"TrigFit[data, n, x] gives the least squares trigonometric fit to data much
as Fit[data, {1, Cos[x], Sin[x], Cos[2x], Sin[2x, ... Cos[n x], Sin[n x]}, x]
would do, except that the calculation is done using Fourier[ ].  The data must
be a list of numbers sampled at equally spaced points over one period
starting at 0.  TrigFit[data, n, {x, L}] can be used if the fundamental
period is not 2Pi.  TrigFit[data, n, {x, x0, x1}] can be used for
data starting at x0 and ending at the sample point immediately before x1."

Begin["NumericalMath`TrigFit`Private`"]

TrigFit[data_List, n_Integer, {x_, x0_:0, x1_}] :=
    Module[{answer = SinAndCos[data, Min[n, Floor[(Length[data]-1)/2]],
							2Pi(x-x0)/(x1-x0)]},
	answer /; answer =!= $Failed
	]

TrigFit[data_List, n_Integer, x_] :=
    Module[{answer = SinAndCos[data, Min[n, Floor[(Length[data]-1)/2]], x]},
	answer /; answer =!= $Failed
	]

SinAndCos[data_, n_, x_] :=
    Module[{fcslist, clist, slist, c0, len, k, rflag},
	len = Length[data];
	fcslist = If[Precision[data] == Infinity, N[data], data];
	rflag = (Union[Im[fcslist], SameTest -> SameQ] === {0});
	clist = Drop[fcslist, 1];
	slist = Reverse[clist];
	fcslist = Prepend[(clist+slist)/2 + I (clist-slist)/2, fcslist[[1]]];
	fcslist = Fourier[fcslist/Sqrt[len]];
	If[!VectorQ[fcslist, NumberQ], Return[$Failed]];
	c0 = First[fcslist];
	clist = Drop[fcslist, 1];
	slist = Reverse[clist];
	{clist, slist} = {Take[clist + slist, n], Take[clist - slist, n]};
	If[rflag, c0 = Re[c0]; clist = Re[clist]; slist = Re[slist]];
	c0 + Table[Cos[k x], {k, n}] . clist - Table[Sin[k x], {k, n}] . slist
	];

End[]  (* NumericalMath`TrigFit`Private` *)

EndPackage[] (* NumericalMath`TrigFit` *)

(* :Tests:
data = Table[1 + 2Sin[x] + 3Cos[2x] + 4 Sin[3x], {x, 0, 2Pi-2Pi/7, 2Pi/7}];
TrigFit[data, 0, x]
TrigFit[data, 1, x]
TrigFit[data, 2, x]
TrigFit[data, 3, x]
TrigFit[data, 4, x]

TrigFit[data, 0, {x, L}]
TrigFit[data, 1, {x, L}]
TrigFit[data, 2, {x, L}]
TrigFit[data, 3, {x, L}]

TrigFit[data, 0, {x, x0, x1}]
TrigFit[data, 1, {x, x0, x1}]
TrigFit[data, 2, {x, x0, x1}]
TrigFit[data, 3, {x, x0, x1}]

TrigFit[data, 3, (x-x0)2Pi/L]
TrigFit[N[data, 40], 4, x]

data = Table[1 + 2Sin[x] + 3Cos[2x] + 4 Sin[3x], {x, 0, 2Pi-2Pi/8, 2Pi/8}];
TrigFit[data, 0, x]
TrigFit[data, 1, x]
TrigFit[data, 2, x]
TrigFit[data, 3, x]
TrigFit[data, 4, x]

TrigFit[data, 0, {x, L}]
TrigFit[data, 1, {x, L}]
TrigFit[data, 2, {x, L}]
TrigFit[data, 3, {x, L}]

TrigFit[data, 0, {x, x0, x1}]
TrigFit[data, 1, {x, x0, x1}]
TrigFit[data, 2, {x, x0, x1}]
TrigFit[data, 3, {x, x0, x1}]

TrigFit[data, 3, (x-x0)2Pi/L]
TrigFit[N[data, 40], 4, x]
*)

