BeginPackage["Spline`", "Algebra`Trigonometry`"];

(*
   Package for the Battle-LeMarie spline wavelets and associated functions --
   Usage Examples are in Spline.ma
   
   Author:  Jack K. Cohen, Colorado School of Mines, 1993
   	    (jkc@keller.mines.colorado.edu)

   Reference: Daubechies, Ten Lectures, SIAM, 1992
   	See pages 146-148.
   
   Convention:  My Fourier transforms are done with the kernel
   	Exp[-2 Pi I x xi]).  Thus, for example, the orthonormalization
	sum add to  1 instead of to Daubechies' 1/2Pi.  And my half-period
	point is 1/2 instead of her Pi, etc.

   Caveat:  I suspect that I am missing a major efficiency, so if you
   	see a better way, my feelings won't be hurt.  (jkc, 10/93)

   Mystery:  There are some commented off versions of the modules that
   	DON'T work for reasons I only dimly understand.  Another place
	for the Mathematica guru to show me how.   (jkc, 10/93)
*)

(* Declaration of public function names in Spline wavelet package *)

SplineNormSquare::usage = "SplineNormSquare[n, xi] -- the orthonormalization sum, Sum[Abs[phi^(xi + n)]^2, {n , -Infinity, Infinity}] for the Battle-Lemarie spline wavelets."

SplinePhiHat::usage = "SplinePhiHat[n, xi] -- the Fourier Transform of the B-spline (these are NOT orthonormal--see SplinePhiHatSharp).  The argument n is the order of the spline to be used in defining the Battle-Lemarie wavelet.  The argument xi is the Fourier variable.  Note: The transform is done with the kernel Exp[-2 Pi I x xi], thus PhiHat has support on [-1, 1]."

SplinePhiHatSharp::usage = "SplinePhiHatSharp[n, xi] -- the Fourier Transform of the Battle-Lemarie wavelet.  The argument n is the order of the spline to be used in defining the wavelet.  The argument xi is the Fourier variable.  Note: The transform is done with the kernel Exp[-2 Pi I x xi], thus PhiHatSharp has support on [-1, 1]."

SplinePhi::usage = "SplinePhi[nu, x, (L)] -- the B-spline (these are NOT orthonormal--see SplinePhiSharp).  The argument n is the order of the spline to be used in defining to be used in defining the Battle-Lemarie wavelet.  The argument x is the independent variable.  SplinePhi is computed by using a Fourier Cosine Series.  The series is constructed for the interval [-L, L], using the argument L (default 10) to define the `effective' support of the scaling function.  Note: SplinePhi is continued periodically outside [-L, L]."

SplinePhiSharp::usage = "SplinePhiSharp[nu, x, (L)] -- the Battle-Lemarie scaling function.  The argument n is the order of the spline to be used in defining to be used in defining the wavelet.  The argument x is the independent variable.  SplinePhiSharp is computed by using a Fourier Cosine Series.  The series is constructed for the interval [-L, L], using the argument L (default 10) to define the `effective' support of the scaling function.  Note: SplinePhiSharp is continued periodically outside [-L, L]."

SplinePsiHat::usage = "SplinePsiHat[n, xi] -- the Fourier Transform of the Battle-Lemarie mother wavelet.  The argument n is the  order of the spline to be used in defining the wavelet.  The argument xi is the Fourier variable.    SplinePsiHat is a complex valued function with the phase factor chosen as Exp[- I pi xi].    Note: The transform is done with the kernel Exp[-2 Pi I x xi]."

SplinePsi::usage = "SplinePsi[n, x, (L)] -- the Battle-Lemarie wavelet.  The argument n is the order of the spline to be used in defining the wavelet.  The argument x is the independent variable.   SplinePsi is computed by using a Fourier Cosine Series.  The series is constructed for the interval [-L, L], using the argument L (default 10) to define the `effective' support of the scaling function.  Note: SplinePsi is continued periodically outside [-L, L]."

SplineM0::usage = "SplineM0[n, xi] -- the Battle-Lemarie characteristic function (low pass or scaling function filter).  The argument n is the  order of the spline to be used in defining the wavelet.  The argument xi is the Fourier variable.   SplineM0 is 1-periodic."

SplineM1::usage = "SplineM1[n, xi] -- the Battle-Lemarie high pass or wavelet filter.  The argument n is the  order of the spline to be used in defining the wavelet.  The argument xi is the Fourier variable.   SplineM1 is 1-periodic."


Begin["`private`"];

SplineNormSquare[n_Integer?Positive, xi_] :=
Module[
        {m = 2n},
        Together@Expand@TrigReduce[
                Sin[Pi xi]^(m+2) /((m+1)! Pi^m) *
                        D[Csc[Pi xi]^2, {xi, m}]
        ]
]

SplinePhiHat[n_Integer?Positive, xi_] :=
Module[
        {phihat},

        phihat =
        If[ xi == 0,
        (* then *) 1,
        (* else *) (Sin[Pi xi]/(Pi xi))^(n+1)
        ];
        
        If[ OddQ[n],
        (* then *) phihat,
        (* else *) Exp[-I Pi xi] phihat
        ]
]

SplinePhiHatSharp[n_Integer?Positive, xi_] :=
	 SplinePhiHat[n, xi] / Sqrt[SplineNormSquare[n, xi]]

SplinePhi[n_Integer?Positive, x_, L_:10] :=
Module[
        {nmax = Floor[2 (40 Exp[-n/2]) L], xi, absphihat, sum},
        
        absphihat[xi_] := Abs@SplinePhiHat[n, xi];
        sum = 1/2 ;     (* first term *) (* phihat[0] = 1 *)
        
        If[ OddQ[n],
        (* Then *) 
                sum += Sum[             (* add the rest of the terms *)
                        Cos[x Pi k/L] absphihat[k/(2L)],
                {k, 1, nmax}],
        (* Else *) 
                sum += Sum[             (* add the rest of the terms *)
                        Cos[(x - 1/2) Pi k/L] absphihat[k/(2L)],
                {k, 1, nmax}]
        ];

        sum / L
]

SplinePhiSharp[n_Integer?Positive, x_, L_:10] :=
Module[
        {nmax = Floor[2 (30/n) L], xi, absphihat, sum},
        
        absphihat[xi_] = Abs@SplinePhiHatSharp[n, xi];
        sum = 1/2 ;     (* first term *) (* phihat[0] = 1 *)
        
        If[ OddQ[n],
        (* Then *) 
                sum += Sum[             (* add the rest of the terms *)
                        Cos[x Pi k/L] absphihat[k/(2L)],
                {k, 1, nmax}],
        (* Else *) 
                sum += Sum[             (* add the rest of the terms *)
                        Cos[(x - 1/2) Pi k/L] absphihat[k/(2L)],
                {k, 1, nmax}]
        ];

        sum / L
]

SplineM0[n_Integer?Positive, x_] :=
Module[
        {phihat, xi},

        phihat[xi_] = SplinePhiHatSharp[n, xi];
        
	If[ IntegerQ[x], (* Watch out for the zeroes of Sin[Pi x] *)
                (* Then *) 1,
                (* Else *) phihat[2 x]/phihat[x] (* 1-periodic by inspection *)
        ]
]

(* This doesn't work
SplineM1[n_Integer?Positive, xi_] = Exp[- I Pi xi]  SplineM0[xi + 1/2]
*)

SplineM1[n_Integer?Positive, x_] :=
Module[
	{phihat, xi, y = x + 1/2},

    phihat[xi_] = SplinePhiHatSharp[n, xi];
        
	If[ IntegerQ[y],
        (* Then *) Exp[- I Pi y],
        (* Else *) Exp[- I Pi y] phihat[2 y]/phihat[y]
    ]
]

(* This doesn't work
SplinePsiHat[n_Integer?Positive, x_] :=
Module[
        {m0, phihat, xi},

        phihat[xi_] := SplinePhiHatSharp[n, xi];
          m0[xi_] := SplineM0[xi];
	        
        Exp[- I Pi x] m0[(x + 1)/2] phihat[x/2]
]
*)

SplinePsiHat[n_Integer?Positive, x_] :=
Module[
	{m0, phihat, xi},

	phihat[xi_] = SplinePhiHatSharp[n, xi];
	
	m0[xi_] :=
	If[ IntegerQ[xi],
		(* Then *) 1,
		(* Else *) phihat[2 xi]/phihat[xi]
	];
	
	Exp[- I Pi x] m0[(x + 1)/2] phihat[x/2]
]

SplinePsi[n_Integer?Positive, x_, L_:10] :=
Module[
        {nmax = Floor[2 (15/n) L], abspsihat, xi, sum},
        
        abspsihat[xi_] := Abs@SplinePsiHat[n, xi];
        sum = 0;        (* first term *)  (* psihat[0] = 0 *)
        
        sum += Sum[             (* add the rest of the terms *)
                Cos[(x - 1/2) Pi k/L] abspsihat[k/(2L)],
        {k, 1, nmax}];

        sum / L
]
                                         
End[];

Protect[
        SplineNormSquare,
        SplinePhiHat,
        SplinePhiHatSharp,
        SplinePhi,
        SplinePhiSharp,
        SplinePsi,
        SplinePsiHat,
        SplineM0,
        SplineM1
]

EndPackage[];