(*:Title: SplineWavelet *) (*:Author: Susumu Sakakibara *) (*:Copyright: Copyright Susumu Sakakibara *) (*:Package Version: 1.01 1.0 Created Dec. 7, 1993 1.01 Modified Dec. 16, 1993 *) (*:Mathematica Version: 2.2 *) (*:Summary: This package provides functions for analyzing equally spaced data using cubic cardinal B-spline wavelets. *) (*:Keywords: Wavelet analysis, Cardinal B-spline wavelets *) (*:Discussion: This package should also be useful for generic wavelets other than cubic spline wavelets. *) (*:Source: C. K. Chui, An Introduction to Wavelets, Academic Press, 1992 *) (*:Warning: All the functions that deal with lists assumes the list of the form {minimumOrder_Integer, dataSequence_List}. *) BeginPackage["SplineWavelet`"] ScalingSample::usage = "ScalingSample[data] gives the list of {x, y[x]} pairs, where y[x] is the scaling function with data taken as the expansion coefficients, evaluated at the integer knots x = n. ScalingSample[data, j] gives the list of {x, y[2^j x]} pairs, evaluated at the dyadic knots x = 2^(-j) n. ScalingSample[data, {j, k}] gives the list of {x, y[2^j x]} pairs, evaluated at the dyadic points x = 2^(-j - k) n. k must be a nonnegative integer." ScalingInterpolate::usage = "ScalingInterpolate[data] gives the list of values y[x] of the scaling function with data taken as the expansion coefficients, evaluated at the integer knots x = n. ScalingInterpolate[data, {k}] gives the same list, but evaluated at the dyadic points x = 2^(-k) n." DisplayScaling::usage = "DisplayScaling[data] generates a plot of the scaling function y[x] with data taken as the expansion coefficients, using samples at the integer knots x = n. DisplayScaling[data, j], j a negative integer, generates a plot of the scaling function y[2^j x], using samples at integer points x = n." FineDisplayScaling::usage = "FineDisplayScaling[data] generates a plot of the scaling function y[x] with data taken as the expansion coefficients, using the samples at dyadic points x = 2^(-4) n, with n being integer. FineDisplayScaling[data, j] generates a plot of the scaling function y[2^j x], using the samples at dyadic points x = 2^(-j - 4) n." FinerScalingCoefficient::usage = "FinerScalingCoefficient[coef] generates scaling function coefficients at one level higher resolution by using the interpolatory graphical display algorithm. FinerScalingCoefficient[coef, n] generates scaling function coefficients at n level higher resolution." WaveletSample::usage = "WaveletSample[data] gives the list of {x, y[x]} pairs, where y[x] is the wavelet with data taken as the expansion coefficients, evaluated at half-integer knots x = 2^(-1) n, with n being integer. WaveletSample[data, j] gives the list of {x, y[2^j x]} pairs, evaluated at the dyadic knots x = 2^(-j - 1) n. WaveletSample[data, {j, k}] gives the list of {x, y[2^j x]} pairs, evaluated at the dyadic points x = 2^(-j - k) n." WaveletInterpolate::usage = "WaveletInterpolate[data] gives the list of values y[x] of the wavelet with data taken as the expansion coefficients, evaluated at the half-integer knots x = 2^(-1) n, with n being integer. WaveletInterpolate[data, {k}] gives the same list, but evaluated at the dyadic points x = 2^(-k) n. k must be a positive integer." DisplayWavelet::usage = "DisplayWavelet[data] generates a plot of the wavelet y[x] with data taken as the expansion coefficients, using samples at the half-integer knots x = 2^(-1) n, with n being integer. DisplayWavelet[data, j], j a negative integer, generates a plot of the wavelet y[2^j x], using samples at the integer points x = n." FineDisplayWavelet::usage = "FineDisplayWavelet[data] generates a plot of the wavelet y[x] with data taken as the expansion coefficients, using samples at dyadic points x = 2^(-4) n, with n being integer. FineDisplayWavelet[data, j] generates a plot of the wavelet y[2^j x], using the samples at the dyadic points x = 2^(-j - 3) n." FindInterpolation::usage = "FindInterpolation[data] gives the scaling function expansion coefficients which interpolates the data, using the fundamental spline." DecomposeToScaling::usage = "DecomposeToScaling[data] gives the scaling function expansion coefficients of one less resolution level by decomposing the scaling function expansion coefficients data." DecomposeToWavelet::usage = "DecomposeToWavelet[data] gives the wavelet expansion coefficients of one less resolution level by decomposing the scaling function expansion coefficients data." PlotData::usage = "PlotData[data] is similar to ListPlot, but takes into account where the data begins." Begin["`Private`"] interpolationOrder = 5; (* The fundamental spline is truncated by this order, inclusive. *) inverseEFOrder = 9; (* The inverse of the Euler-Frobenius polynomial is truncated by this order, inclusive. *) (* Definitions of several convolution operators *) ListConvolve[a_List, b_List] := Module[{tList}, tList = Partition[ Flatten[ Join[ Table[0, {Length[Last[a]] - 1}], Last[b], Table[0, {Length[Last[a]] - 1}] ] ], Length[Last[a]], 1 ]; {First[a] + First[b], (Reverse[Last[a]] . #)& /@ tList} ] UpsampleConvolve[a_List, b_List] := Module[{tList}, tList = Partition[ Flatten[ Join[ Table[0, {Length[Last[a]] - 2}], (# /. # -> {0, #})& /@ Last[b], Table[0, {Length[Last[a]] - 1}] ] ], Length[Last[a]], 1 ]; {First[a] + 2 First[b], (Reverse[Last[a]] . #)& /@ tList} ] DecomposeConvolve[a_List, b_List] := Module[{tList, tLowest, tHighest}, tLowest = First[a] + First[b]; tHighest = Length[Last[a]] + Length[Last[b]] + tLowest - 2; tList = Partition[ Flatten[ Join[ Table[0, {Length[Last[a]] - 1 - Mod[tLowest, 2]}], Last[b], Table[0, {Length[Last[a]] - 1 - Mod[tHighest, 2]}] ] ], Length[Last[a]], 2 ]; {Sign[tLowest] Floor[Abs[tLowest]/2], (Reverse[Last[a]] . #)& /@ tList} ] (* The recursion relation to generate cardinal B-spline at knots. *) bSp[2, 1] := 1; bSp[2, k_Integer] := 0 /; k =!= 1; bSp[m_Integer, k_Integer] := k/(m - 1) bSp[m - 1, k] + (m - k)/(m - 1) bSp[m - 1, k - 1]; (* Defininion of Euler-Frobenius polynomial *) EulerFrobenius[n_Integer, z_] := Expand[n! Sum[bSp[n + 1, j + 1] z^j, {j, 0, n - 1}]] (* Coefficients used for interpolation by the fundamental spline *) `splAtKnots = {0, Table[bSp[4, n], {n, 0, 4}]}; (* Two-scale relation coefficients, one for scaling function, and the other for wavelets. *) `pCoef = {0, Table[Binomial[4, k]/8, {k, 0, 4}]}; `qCoef = Module[{a, b}, a = {0, Table[bSp[8, n + 1], {n, 0, 6}]}; b = ListConvolve[pCoef, a]; b /. b[[2]] -> b[[2]] {1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1} ]; (* Functions to produce functions in the scaling function space, and functions in the wavelet space, using the interpolatory graphical display algorithm. *) ScalingInterpolate[fmtList_List] := ScalingInterpolate[fmtList, {0}] ScalingInterpolate[fmtList_List, {k_Integer}] := Last[ ListConvolve[ splAtKnots, Nest[UpsampleConvolve[pCoef, #]&, fmtList, k] ] ] ScalingSample[fmtList_List, j_Integer:0] := ScalingSample[fmtList, {j, 0}]; ScalingSample[fmtList_List, {j_, k_Integer:0}] := Module[{ip, dpl}, ip = Nest[UpsampleConvolve[pCoef, #]&, fmtList, k]; dpl = ListConvolve[splAtKnots, ip]; Transpose[{ Range[dpl[[1]], Length[dpl[[2]]] - 1 + dpl[[1]]]/2^(k + j), dpl[[2]]} ] ] FinerScalingCoefficient[coef_List, k_Integer:1] := Nest[UpsampleConvolve[pCoef, #]&, coef, k]; FineDisplayScaling[fmtList_List, j_Integer:0, opt___] := ListPlot[ScalingSample[fmtList, {j, 4}], opt, PlotJoined -> True] DisplayScaling[fmtList_List, j_Integer:0, opt___] := ListPlot[ScalingSample[fmtList, {j, -j}], opt, PlotJoined -> True] WaveletInterpolate[fmtList_List] := WaveletInterpolate[fmtList, {1}] WaveletInterpolate[fmtList_List, {k_Integer}] := Module[{ip, dpl, tmpList}, tmpList = UpsampleConvolve[qCoef, fmtList]; ip = Nest[UpsampleConvolve[pCoef, #]&, tmpList, k - 1]; dpl = ListConvolve[splAtKnots, ip]; Last[dpl] ] WaveletSample[fmtList_List, j_Integer:0] := WaveletSample[fmtList, {j, 1}]; WaveletSample[fmtList_List, {j_Integer, k_Integer:1}] := Module[{ip, dpl, tmpList}, tmpList = UpsampleConvolve[qCoef, fmtList]; ip = Nest[UpsampleConvolve[pCoef, #]&, tmpList, k - 1]; dpl = ListConvolve[splAtKnots, ip]; Transpose[{ Range[dpl[[1]], Length[dpl[[2]]] - 1 + dpl[[1]]]/2^(k + j), dpl[[2]]} ] ] FineDisplayWavelet[fmtList_List, j_Integer:0, opt___] := ListPlot[WaveletSample[fmtList, {j, 4}], opt, PlotJoined -> True] DisplayWavelet[fmtList_List, j_Integer:0, opt___] := Module[{jj}, jj = If[j == 0, 1, -j]; ListPlot[WaveletSample[fmtList, {j, jj}], opt, PlotJoined -> True] ] (* Coefficients defining the fundamental spline. interpolationOrder is given at the beginning. *) `halfSeq = Sqrt[3] Table[(Sqrt[3] - 2)^k, {k, 0, interpolationOrder}] //N; InterpolationSeq = {-interpolationOrder - 2, Join[Reverse[Drop[halfSeq, 1]], halfSeq]}; (* The function that generates interpolation of the data sequence. *) FindInterpolation[data_] := ListConvolve[InterpolationSeq, data]; (* Below are parameters needed to decompose a function into the scaling function space and wavelet space of lower resolution level. *) `alpha = -0.5352804307964381635; `beta = -0.1225546151923266906; `gamma = -0.009148694809608275136; {`a, `b, `c}; `substRule = {a -> alpha, b -> beta, c -> gamma}; `coefA = (a^3*b*c)/ ((-1 + a^2)*(a - b)*(-1 + a*b)*(a - c)*(-1 + a*c)); `coefB = -((a*b^3*c)/ ((a - b)*(-1 + a*b)*(-1 + b^2)*(b - c)*(-1 + b*c))); `coefC = (a*b*c^3)/ ((a - c)*(b - c)*(-1 + a*c)*(-1 + b*c)*(-1 + c^2)); `tmpG = Expand[(1 + z)^4/16 EulerFrobenius[7, z]]; coefListG = {-7, Join[{tmpG /. z -> 0}, Table[Coefficient[tmpG, z^k], {k, 10}] ]}; `tmpH = Expand[-7! (1 - z)^4/16]; coefListH = {-7, Join[{tmpH /. z -> 0}, Table[Coefficient[tmpH, z^k], {k, 4}] ]}; `defineTruncation[n_] := Module[{halfTable}, halfTable = N[ coefA a^Range[0, n] + coefB b^Range[0, n] + coefC c^Range[0, n] /. substRule ]; {-n, Join[Reverse[Drop[halfTable, 1]], halfTable]} ] (* inverseEFOrder is given at the beginning *) `trList = defineTruncation[inverseEFOrder]; `gList = UpsampleConvolve[coefListG, trList]; `hList = UpsampleConvolve[coefListH, trList]; DecomposeToScaling[scalingCoef_] := DecomposeConvolve[gList, scalingCoef]; DecomposeToWavelet[scalingCoef_] := DecomposeConvolve[hList, scalingCoef]; (* Plotting function especially designed for the list in the present format. *) PlotData[data_, opt___] := ListPlot[Transpose[{ Range[data[[1]], Length[data[[2]]] - 1 + data[[1]]], data[[2]]}], opt, PlotStyle -> {PointSize[0.016]}] End[] (* `Private` *) Protect[ ScalingInterpolate, DisplayScaling, FineDisplayScaling, WaveletInterpolate, DisplayWavelet, FineDisplayWavelet, ScalingSample, WaveletSample, FindInterpolation, DecomposeToScaling, DecomposeToWavelet, FinerScalingCoefficient, PlotData ] EndPackage[] (* `SplineWavelet` *)