BeginPackage["WaveletTools`", {"Wavelets`","Graphics`Legend`"}]
(* WaveletTools implements functions to plot scaling and wavelet functions
   in the time and frequency domains. Some numerical wavelet utility
   functions are also implemented.
   Note: the plotting routines make use of the cascade algorithm
   implemented in the Mathematica wavelets package Wavelets, if you
   do not have the package you can still use the other routines
   without it.
   by Bernard Cena
   Version 0.60 - 29/03/2000
*)

plotSW::usage =
  "plotSW[g, n, pr] or plotSW[g, f, n, pr] plots the scaling and wavelet functions derived from filter g (and f) using n number of recursive refinement steps. The plotting range (pr) can also be specified as per PlotRange option in Plot."

printSW::usage =
  "printSW[g, n, s, pr] or printSW[g, f, s, pr] prints and EPS plot of the scaling and wavelet functions derived from filter g (and f) using n number of recursive refinement steps into file with prefix s. The plotting range (pr) can also be specified as per PlotRange option in Plot."

plotTD::usage =
  "plotTD[g, n, pr] or plotTD[g, f, n, pr] plots orthogonal or biorthogonal scaling and wavelet analysis and reconstruction functions based on filters g (and f) using n number of recursive refinement steps. The plotting range (pr) can also be specified as per PlotRange option in Plot."

printTD::usage =
  "printTD[g, n, s, pr] or printTD[g, f, n, s, pr] prints an EPS plot to file with prefix s of an orthogonal or biorthogonal scaling and wavelet analysis and reconstruction functions based on filters g (and f) using n number of recursive refinement steps. The plotting range (pr) can also be specified as per PlotRange option in Plot."

plotFD::usage =
  "plotFD[g, n, pr] or printFD[g, f, n, pr] plots the frequency responses of low- and high-pass filters based on g (and f) and the freqency responses of the scaling and wavelet functions using n-th order refinement. The plotting range (pr) can also be specified as per PlotRange option in Plot."

printFD::usage =
  "printFD[g, n, s, pr] or printFD[g, f, n, s, pr] prints an EPS plot to file with prefix s of the frequency responses of low- and high-pass filters based on g (and f) and the freqency responses of the scaling and wavelet function using n-th order refinement. The plotting range (pr) can also be specified as per PlotRange option in Plot."
 
trimfilt::usage =
  "trimfilt[g] trims leading and trailing zero coefficients from filter g."

zfilt::usage =
  "zfilt[g] returns a centered z-transform of filter coefficients g."

plotRoots::usage =
  "plotRoots[g] plots the roots of the z-transform of the filter g."

calcSP::usage =
  "calcSP[g] calculates approximation order and smoothness (finite energy derivatives) for filter g using the procedure described in Strang and Nguyen, 'Wavelets and Filter Banks'. Uses numerical eigenanalysis routine for speed and rationalizes eigenvalues with accuracy of 10^-5. Returns {s,p}."

calcAppPow::usage =
  "calcAppPow[g, d] calculates the approximation power constant as derived by Michael Unser et al, where g is the filter coefficient list. An optional parameter d can be set to the number of zeros at z=-1 of the filter to avoid numerical problems for long filters with approximated coefficients. Discretisation parameters can be changed in the code."

calcSFMom::usage =
  "calcSFMom[g, n, sh] calculates the n-th moment of the scaling function generated by filter g centered at x=sh (default sh=0)."

visRecFilter::usage =
  "visRecFilter[n] returns the visualisation biorthogonal reconstruction filter of order n using a closed form formula. The even order filters are computed by the closed form formula for n-1 order biorthogonal coiflet as derived by Jun Tian in 'Vanishing Moments and Wavelet Approximation'."

visAnaFilter::usage =
  "visAnaFilter[n] returns the visualisation biorthogonal analysis filter to match of order n using a closed form formula. The even order filters are computed by the closed form formula for n-1 order biorthogonal coiflet as derived by Jun Tian in 'Vanishing Moments and Wavelet Approximation'."

Begin["WaveletTools`Private`"]

plotSW[g_List, refine_Integer, pltrng_:All] := Module[
	{ng = 2 g/((Sqrt[2] Plus@@g)),
	lg = (-Floor[((Length[g] - 1))/2]), filt, phi1, psi1},
	filt = {{ng, lg}, {ng, lg}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True]; 
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	Show[{ListPlot[phi1, PlotRange -> pltrng, PlotJoined -> True,
		DisplayFunction -> Identity],
		ListPlot[psi1, PlotRange -> pltrng, PlotJoined -> True,
        PlotStyle -> {Dashing[{0.01}]}, DisplayFunction -> Identity]},
        Frame->True, DisplayFunction -> $DisplayFunction]
];

printSW[g_List, refine_Integer, prefix_String:"swplot", pltrng_:All] := Module[
	{ng = 2 g/((Sqrt[2] Plus@@g)),
	lg = (-Floor[((Length[g] - 1))/2]), filt, phi1, psi1},
	filt = {{ng, lg}, {ng, lg}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True]; 
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	pl=Show[{ListPlot[phi1, PlotRange -> pltrng, PlotJoined -> True,
		DisplayFunction -> Identity],
		ListPlot[psi1, PlotRange -> pltrng, PlotJoined -> True,
        PlotStyle -> {Dashing[{0.01}]}, DisplayFunction -> Identity]},
        Frame->True];
	Display[prefix <> "_sw.eps", pl, "EPS"];
];

plotSW[g_List, f_List, refine_Integer, pltrng_:All] := Module[
	{ng = 2 g/((Sqrt[2] Plus@@g)), lg = (-Floor[((Length[g] - 1))/2]),
	nf = 2 f/((Sqrt[2] Plus@@f)), lf = (-Floor[((Length[f] - 1))/2]),
	filt, phi1, psi1, phi2, psi2},
	filt = {{ng, lg}, {nf, lf}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True]; 
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	phi2 = ScalingFunction[Reverse[filt], refine,
		InverseTransform -> True]; 
	psi2 = Wavelet[Reverse[filt], refine, InverseTransform -> True];
	Show[{ListPlot[phi1, PlotRange -> pltrng, PlotJoined -> True,
		DisplayFunction -> Identity],
		ListPlot[psi1, PlotRange -> pltrng, PlotJoined -> True,
        PlotStyle -> {Dashing[{0.01}]}, DisplayFunction -> Identity]},
        Frame->True, DisplayFunction -> $DisplayFunction];
	Show[{ListPlot[phi2, PlotRange -> pltrng, PlotJoined -> True,
		DisplayFunction -> Identity],
		ListPlot[psi2, PlotRange -> pltrng, PlotJoined -> True,
        PlotStyle -> {Dashing[{0.01}]}, DisplayFunction -> Identity]},
        Frame->True, DisplayFunction -> $DisplayFunction];
];

printSW[g_List, f_List, refine_Integer, prefix_String:"swplot", pltrng_:All] := Module[
	{ng = 2 g/((Sqrt[2] Plus@@g)), lg = (-Floor[((Length[g] - 1))/2]),
	nf = 2 f/((Sqrt[2] Plus@@f)), lf = (-Floor[((Length[f] - 1))/2]),
	filt, phi1, psi1, phi2, psi2, pl1, pl2},
	filt = {{ng, lg}, {nf, lf}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True]; 
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	phi2 = ScalingFunction[Reverse[filt], refine,
		InverseTransform -> True]; 
	psi2 = Wavelet[Reverse[filt], refine, InverseTransform -> True];
	pl1=Show[{ListPlot[phi1, PlotRange -> pltrng, PlotJoined -> True,
		DisplayFunction -> Identity],
		ListPlot[psi1, PlotRange -> pltrng, PlotJoined -> True,
        PlotStyle -> {Dashing[{0.01}]}, DisplayFunction -> Identity]},
        Frame->True];
	pl2=Show[{ListPlot[phi2, PlotRange -> pltrng, PlotJoined -> True,
		DisplayFunction -> Identity],
		ListPlot[psi2, PlotRange -> pltrng, PlotJoined -> True,
        PlotStyle -> {Dashing[{0.01}]}, DisplayFunction -> Identity]},
        Frame->True];
	Display[prefix <> "_asw.eps", pl1, "EPS"];
	Display[prefix <> "_ssw.eps", pl1, "EPS"];
];

plotTD[g_List, refine_Integer, pltrng_:All] := Module[
	{ng = 2  g/((Sqrt[2] Plus@@g)), lg = (-Floor[((Length[g] - 1))/2]),
	filt, phi1, psi1, pls1, plw1, mn, mx, pr},
	filt = {{ng, lg}, {ng, lg}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True]; 
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	mn = Min[phi1[[1,1]], psi1[[1,1]]];
	mx = Max[phi1[[-1,1]], psi1[[-1,1]]];
	pr = pltrng;
	If[pr == Null, pr = {{mn, mx}, All}];
	pls1 = ListPlot[phi1, PlotJoined -> True, PlotRange -> pr, 
		DisplayFunction -> Identity, 
		PlotLabel -> "Scaling Function"];
	plw1 = ListPlot[psi1, PlotJoined -> True, PlotRange -> pr, 
		DisplayFunction -> Identity, 
		PlotLabel -> "Wavelet Function"];
	Show[GraphicsArray[{pls1, plw1}]]
];

printTD[g_List, refine_Integer, prefix_String:"tdplot", pltrng_:All] := Module[
	{ng = 2  g/((Sqrt[2] Plus@@g)), lg = (-Floor[((Length[g] - 1))/2]),
	filt, phi1, psi1, pls1, plw1, mn, mx, pr},
	filt = {{ng, lg}, {ng, lg}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True];
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	mn = Min[phi1[[1,1]], psi1[[1,1]]];
	mx = Max[phi1[[-1,1]], psi1[[-1,1]]];
	pr = pltrng;
	If[pr == Null, pr = {{mn, mx}, All}];
	pls1 = ListPlot[phi1, PlotJoined -> True, PlotRange -> pr, 
		DisplayFunction -> Identity];
	plw1 = ListPlot[psi1, PlotJoined -> True, PlotRange -> pr, 
		DisplayFunction -> Identity];
(* print all plots *)
	Display[prefix <> "_sf.eps", pls1, "EPS"];
	Display[prefix <> "_wf.eps", plw1, "EPS"];
];

plotTD[g_List, f_List, refine_Integer, pltrng_:Null] := Module[
	{ng = 2 g/((Sqrt[2] Plus@@g)), lg = (-Floor[((Length[g] - 1))/2]),
	nf = 2 f/((Sqrt[2] Plus@@f)), lf = (-Floor[((Length[f] - 1))/2]),
	filt, phi1, psi1, phi2, psi2, pls1, plw1, pls2, plw2, mn, mx, pr},
	filt = {{ng, lg}, {nf, lf}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True]; 
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	phi2 = ScalingFunction[Reverse[filt], refine,
		InverseTransform -> True]; 
	psi2 = Wavelet[Reverse[filt], refine, InverseTransform -> True];
	mn = Min[phi1[[1,1]], psi1[[1,1]], phi2[[1,1]], psi2[[1,1]]];
	mx = Max[phi1[[-1,1]], psi1[[-1,1]], phi2[[-1,1]], psi2[[-1,1]]];
	pr = pltrng;
	If[pr == Null, pr = {{mn, mx}, All}];
	pls1 = ListPlot[phi1, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity,
		PlotLabel -> "Analysis Scaling Function"];
	plw1 = ListPlot[psi1, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity,
		PlotLabel -> "Analysis Wavelet Function"];
	pls2 = ListPlot[phi2, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity,
		PlotLabel -> "Synthesis Scaling Function"];
	plw2 = ListPlot[psi2, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity,
		PlotLabel -> "Synthesis Wavelet Function"];
	Show[GraphicsArray[{{pls1, pls2}, {plw1, plw2}}]]
];

printTD[g_List, f_List, refine_Integer, prefix_String:"tdplot", pltrng_:Null] := Module[
	{ng = 2 g/((Sqrt[2] Plus@@g)), lg = (-Floor[((Length[g] - 1))/2]),
	nf = 2 f/((Sqrt[2] Plus@@f)), lf = (-Floor[((Length[f] - 1))/2]),
	filt, phi1, psi1, phi2, psi2, pls1, plw1, pls2, plw2, mn, mx, pr},
	filt = {{ng, lg}, {nf, lf}}; 
	phi1 = ScalingFunction[filt, refine, InverseTransform -> True]; 
	psi1 = Wavelet[filt, refine, InverseTransform -> True];
	phi2 = ScalingFunction[Reverse[filt], refine,
		InverseTransform -> True]; 
	psi2 = Wavelet[Reverse[filt], refine, InverseTransform -> True];
	mn = Min[phi1[[1,1]], psi1[[1,1]], phi2[[1,1]], psi2[[1,1]]];
	mx = Max[phi1[[-1,1]], psi1[[-1,1]], phi2[[-1,1]], psi2[[-1,1]]];
	pr = pltrng;
	If[pr == Null, pr = {{mn, mx}, All}];
	pls1 = ListPlot[phi1, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity];
	plw1 = ListPlot[psi1, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity];
	pls2 = ListPlot[phi2, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity];
	plw2 = ListPlot[psi2, PlotRange -> pr, PlotJoined -> True,
		DisplayFunction -> Identity];
(* print all plots *)
	Display[prefix <> "_asf.eps", pls1, "EPS"];
	Display[prefix <> "_awf.eps", plw1, "EPS"];
	Display[prefix <> "_ssf.eps", pls2, "EPS"];
	Display[prefix <> "_swf.eps", plw2, "EPS"];
];

plotFD[g_List, refine_Integer, plrng_:All] := Module[
	{ng = g/Plus@@g,gl=Length[g]-1,gt,glf,ghf,sf,wf,pga,pgp,psa,pwa},
(* time domain definition *)
	gt[x_/;x<0||x>gl]:=0;
	gt[x_]:=ng[[x+1]];
(* low-pass frequency response *)
	glf[o_]:=glf[o]=Sum[gt[k] Exp[I k o \[Pi]],{k,0,gl}]//N;
(* high-pass frequency response *)
	ghf[o_]:=ghf[o]=Sum[(-1)^k gt[gl-k] Exp[I k o \[Pi]],{k,0,gl}]//N;
(* plot filter low/high-pass frequency amplitude response *)
	pga=Plot[{Abs[glf[x]], Abs[ghf[x]]}, {x, 0, 1},
    PlotStyle -> {{GrayLevel[0]}, {Dashing[{0.01}]}},
    Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, {1}},
    PlotLabel -> "Amplitude responses", DisplayFunction -> Identity];
(* plot filter low/high-pass frequency phase response *)
	pgp=Plot[Arg[glf[x]], {x, 0, 1},
    Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, None},
    PlotLabel -> "Phase response", DisplayFunction -> Identity];
(* scaling and wavelet function frequency response *)
	sf[o_]:=sf[o]=Product[glf[o/2^j],{j,1,refine}]//N;
	wf[o_]:=wf[o]=ghf[o/2] sf[o/2]//N;
(* plot scaling and wavelet function frequency amplitude response *)
	psa=ListPlot[Table[{x, Abs[sf[x]]}, {x, 0, 20, 1/2^3}],
    PlotJoined -> True, PlotRange -> plrng, DisplayFunction -> Identity,
    Ticks -> {Table[i,{i, 0, 20, 2}], {1}},
    PlotLabel -> "Scaling function amplitude"];
	pwa=ListPlot[Table[{x, Abs[wf[x]]}, {x, 0, 20, 1/2^3}],
    PlotJoined -> True, PlotRange -> plrng, DisplayFunction -> Identity,
    Ticks -> {Table[i,{i, 0, 20, 2}], {1}},
    PlotLabel -> "Wavelet function amplitude"];
(* display all plots *)
    Show[GraphicsArray[{{pga,pgp},{psa,pwa}}]]
];

printFD[g_List, refine_Integer, prefix_String:"fdplot", plrng_:All] := Module[
	{ng = g/Plus@@g,gl=Length[g]-1,gt,glf,ghf,sf,wf,pga,pgp,psa,pwa},
(* time domain definition *)
	gt[x_/;x<0||x>gl]:=0;
	gt[x_]:=ng[[x+1]];
(* low-pass frequency response *)
	glf[o_]:=gf[o]=Sum[gt[k] Exp[I k o \[Pi]],{k,0,gl}]//N;
(* high-pass frequency response *)
	ghf[o_]:=ghf[o]=Sum[(-1)^k gt[gl-k] Exp[I k o \[Pi]],{k,0,gl}]//N;
(* plot filter low/high-pass frequency amplitude response *)
	pga=Plot[{Abs[glf[x]], Abs[ghf[x]]}, {x, 0, 1},
    PlotStyle -> {{GrayLevel[0]}, {Dashing[{0.01}]}},
    Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, {1}},
    DisplayFunction -> Identity];
(* plot filter low/high-pass frequency phase response *)
	pgp=Plot[Arg[glf[x]], {x, 0, 1},
    Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, None},
    DisplayFunction -> Identity];
(* scaling and wavelet function frequency response *)
	sf[o_]:=sf[o]=Product[glf[o/2^j],{j,1,refine}]//N;
	wf[o_]:=wf[o]=ghf[o/2] sf[o/2]//N;
(* plot scaling and wavelet function frequency amplitude response *)
	psa=ListPlot[Table[{x, Abs[sf[x]]}, {x, 0, 20, 1/2^3}],
    PlotJoined -> True, PlotRange -> plrng, DisplayFunction -> Identity,
    Ticks -> {Table[i,{i, 0, 20, 2}], {1}}];
	pwa=ListPlot[Table[{x, Abs[wf[x]]}, {x, 0, 20, 1/2^3}],
    PlotJoined -> True, PlotRange -> plrng, DisplayFunction -> Identity,
    Ticks -> {Table[i,{i, 0, 20, 2}], {1}}];
(* print all plots *)
	Display[prefix <> "_far.eps", pga, "EPS"];
	Display[prefix <> "_fpr.eps", pgp, "EPS"];
	Display[prefix <> "_sar.eps", psa, "EPS"];
	Display[prefix <> "_war.eps", pwa, "EPS"];
];

plotFD[g_List, f_List, refine_Integer, plrng_:All] := Module[
	{ng = g/Plus@@g, gl=Length[g]-1, nf = f/Plus@@f, fl=Length[f]-1,
	gt, ft, glf, ghf, flf, fhf, sf1 ,wf1, sf2, wf2, pga1 ,pga2,
	pgp1, pgp2, psa1, psa2, pwa1, pwa2},
(* time domain definitions *)
	gt[x_/;x<0||x>gl]:=0;
	gt[x_]:=ng[[x+1]];
	ft[x_/;x<0||x>fl]:=0;
	ft[x_]:=nf[[x+1]];
(* low-pass frequency responses *)
	glf[o_]:=glf[o]=Sum[gt[k] Exp[I k o \[Pi]],{k,0,gl}]//N;
	flf[o_]:=flf[o]=Sum[ft[k] Exp[I k o \[Pi]],{k,0,fl}]//N;
(* high-pass frequency response *)
	ghf[o_]:=ghf[o]=Sum[(-1)^k gt[gl-k] Exp[I k o \[Pi]],{k,0,gl}]//N;
	fhf[o_]:=fhf[o]=Sum[(-1)^k ft[fl-k] Exp[I k o \[Pi]],{k,0,fl}]//N;
(* plot filter low/high-pass frequency amplitude response *)
	pga1 = Plot[{Abs[glf[x]], Abs[fhf[x]]}, {x, 0, 1},
	    PlotStyle -> {{GrayLevel[0]}, {Dashing[{0.01}]}},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, {1}},
    	PlotLabel -> "Analysis Amplitude Responses",
    	DisplayFunction -> Identity];
	pga2 = Plot[{Abs[flf[x]], Abs[ghf[x]]}, {x, 0, 1},
	    PlotStyle -> {{GrayLevel[0]}, {Dashing[{0.01}]}},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, {1}},
    	PlotLabel -> "Synthesis Amplitude Responses",
    	DisplayFunction -> Identity];
(* plot filter low/high-pass frequency phase response *)
	pgp1 = Plot[Arg[glf[x]], {x, 0, 1},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, None},
    	PlotLabel -> "Analysis Phase Response",
    	DisplayFunction -> Identity];
	pgp2 = Plot[Arg[flf[x]], {x, 0, 1},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, None},
    	PlotLabel -> "Synthesis Phase Response",
    	DisplayFunction -> Identity];
(* scaling and wavelet function frequency response *)
	sf1[o_]:=sf1[o]=Product[glf[o/2^j],{j,1,refine}]//N;
	wf1[o_]:=wf1[o]=fhf[o/2] sf1[o/2]//N;
	sf2[o_]:=sf2[o]=Product[flf[o/2^j],{j,1,refine}]//N;
	wf2[o_]:=wf2[o]=ghf[o/2] sf2[o/2]//N;
(* plot scaling and wavelet function frequency amplitude response *)
	psa1 = ListPlot[Table[{x, Abs[sf1[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}},
    	PlotLabel -> "Analysis scaling function amplitude"];
	pwa1 = ListPlot[Table[{x, Abs[wf1[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}},
    	PlotLabel -> "Analysis wavelet function amplitude"];
	psa2 = ListPlot[Table[{x, Abs[sf2[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}},
    	PlotLabel -> "Synthesis scaling function amplitude"];
	pwa2 = ListPlot[Table[{x, Abs[wf2[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}},
    	PlotLabel -> "Synthesis wavelet function amplitude"];
(* display all plots *)
    Show[GraphicsArray[{{pga1,pgp1},{pga2, pgp2},
    	{psa1,pwa1},{psa2,pwa2}}]]
];

printFD[g_List, f_List, refine_Integer, prefix_String:"fdplot", plrng_:All] := Module[
	{ng = g/Plus@@g, gl=Length[g]-1, nf = f/Plus@@f, fl=Length[f]-1,
	gt, ft, glf, ghf, flf, fhf, sf1 ,wf1, sf2, wf2, pga1 ,pga2,
	pgp1, pgp2, psa1, psa2, pwa1, pwa2},
(* time domain definitions *)
	gt[x_/;x<0||x>gl]:=0;
	gt[x_]:=ng[[x+1]];
	ft[x_/;x<0||x>fl]:=0;
	ft[x_]:=nf[[x+1]];
(* low-pass frequency responses *)
	glf[o_]:=glf[o]=Sum[gt[k] Exp[I k o \[Pi]],{k,0,gl}]//N;
	flf[o_]:=flf[o]=Sum[ft[k] Exp[I k o \[Pi]],{k,0,fl}]//N;
(* high-pass frequency response *)
	ghf[o_]:=ghf[o]=Sum[(-1)^k gt[gl-k] Exp[I k o \[Pi]],{k,0,gl}]//N;
	fhf[o_]:=fhf[o]=Sum[(-1)^k ft[fl-k] Exp[I k o \[Pi]],{k,0,fl}]//N;
(* plot filter low/high-pass frequency amplitude response *)
	pga1 = Plot[{Abs[glf[x]], Abs[fhf[x]]}, {x, 0, 1},
	    PlotStyle -> {{GrayLevel[0]}, {Dashing[{0.01}]}},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, {1}},
    	DisplayFunction -> Identity];
	pga2 = Plot[{Abs[flf[x]], Abs[ghf[x]]}, {x, 0, 1},
	    PlotStyle -> {{GrayLevel[0]}, {Dashing[{0.01}]}},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, {1}},
    	DisplayFunction -> Identity];
(* plot filter low/high-pass frequency phase response *)
	pgp1 = Plot[Arg[glf[x]], {x, 0, 1},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, None},
    	DisplayFunction -> Identity];
	pgp2 = Plot[Arg[flf[x]], {x, 0, 1},
    	Ticks -> {{{0,"0"},{1/2,"0.5"},{1,"1"}}, None},
    	DisplayFunction -> Identity];
(* scaling and wavelet function frequency response *)
	sf1[o_]:=sf1[o]=Product[glf[o/2^j],{j,1,refine}]//N;
	wf1[o_]:=wf1[o]=fhf[o/2] sf1[o/2]//N;
	sf2[o_]:=sf2[o]=Product[flf[o/2^j],{j,1,refine}]//N;
	wf2[o_]:=wf2[o]=ghf[o/2] sf2[o/2]//N;
(* plot scaling and wavelet function frequency amplitude response *)
	psa1 = ListPlot[Table[{x, Abs[sf1[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}}];
	pwa1 = ListPlot[Table[{x, Abs[wf1[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}}];
	psa2 = ListPlot[Table[{x, Abs[sf2[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}}];
	pwa2 = ListPlot[Table[{x, Abs[wf2[x]]}, {x, 0, 20, 1/2^3}],
    	PlotJoined -> True, PlotRange -> plrng,
    	DisplayFunction -> Identity,
    	Ticks -> {Table[i,{i, 0, 20, 2}], {1}}];
(* print all plots *)
	Display[prefix <> "_afar.eps", pga1, "EPS"];
	Display[prefix <> "_sfar.eps", pga2, "EPS"];
	Display[prefix <> "_afpr.eps", pgp1, "EPS"];
	Display[prefix <> "_sfpr.eps", pgp2, "EPS"];
	Display[prefix <> "_asar.eps", psa1, "EPS"];
	Display[prefix <> "_ssar.eps", psa2, "EPS"];
	Display[prefix <> "_awar.eps", pwa1, "EPS"];
	Display[prefix <> "_swar.eps", pwa2, "EPS"];
];

trimfilt[g_List] := Module[
	(* removes leading and trailing zero coefficients from filter g. *)
	{f=g},
	While[First[f]==0, f=Drop[f,1]];
	While[Last[f]==0, f=Drop[f,-1]];
	f
];

zfilt[g_List] := Module[
	(* returns a centered z-transform of filter. *)
	{gl, hgl, even},
	gl = Length[g]; hgl = Floor[gl/2]; even = If[EvenQ[gl],1,0];
	Plus @@ MapThread[# Global`z^(#2) &,{g, Range[-hgl+even, hgl]}]
];

plotRoots[g_List] :=
	(* plots roots of the z-transform of the filter g. *)
	Show[ Graphics[Circle[{0,0},1]],
		Graphics[{PointSize[0.02], Point[{Re[#],Im[#]}]&
		/@ Last/@First/@((g//trimfilt//zfilt)==0//NSolve)}],
		AspectRatio -> Automatic, Axes -> True
	]

calcSP[g_List] := Module[
	(* returns {p,s} for LPF given eigenvalues of T
	  the autocorrelation matrix of filter g. *)
	{ng = g / Plus@@g,k,pos,ls,maxk,gl,f,t,tm,l,v},
	gl = Length[g]-1;
	f[x_/;x<0||x>gl]:=0;
	f[x_]:=ng[[x+1]];
	fm=Table[2f[2 i-j],{i,0,gl-1},{j,0,gl-1}];
	{l,v}=Eigensystem[N[fm,100]];
	ls=Rationalize[Abs[l],10^(-8)];
	If[Last[Sort[ls]]!=1,Return[{0,0}]];
	t[x_/;x<0||x>2 gl]:=0;
	t[x_]:=Sum[f[k] f[k-x+gl],{k,0,gl}];
	tm=Table[2t[2 i-j+1],{i,0,2 gl-2},{j,0,2 gl-2}];
	{l,v}=Eigensystem[N[tm,100]];
	k=0;ls=Rationalize[Abs[l],10^(-8)];maxk=Length[l];
	While[k<maxk,
		pos=Position[ls,1/2^k];
		If[Length[pos]>0,
			ls=Delete[ls,pos[[1,1]]],
			k--;Break[]];
			k++;
		];
	If[Length[ls]>0,
		{(k+1)/2,N[-Log[4,Last[Sort[ls]]]]},
		{(k+1)/2,N[(k+1)/2-1/2]}
	]
];

calcAppPow[g_List] := Module[
	(* Calculates the approximation power constant as derived by
	   Michael Unser et al. *)
	{ng = g / Plus@@g, scaleres = 6, calcres = 100, len = Length[g],
		even = If[EvenQ[Length[g]], 1, 0], hflen = Floor[Length[g]/2],
		h, dn, dh, g2, f, nf, hf, sf, ss},
	h[z_] := Plus @@ MapThread[# z^(#2) &,
		{ng, Range[-hflen, hflen-even]}];
	dn = 1;
	dh = Chop[D[h[z^(-1)], {z, dn}] /. z -> -1];
	While[dh == 0,
		dn = dn+1;
		dh = Chop[D[h[z^(-1)], {z, dn}] /. z -> -1];
	];
	g2 = CoefficientList[ Expand[ PolynomialQuotient[
		Expand[ z^hflen * h[z] ], ((1+z)/2)^dn, z]] , z];
	f[z_] := Plus@@MapThread[# z^(#2) &, {g2, Range[0, Length[g2] - 1]}];
	nf = Length[ng] - 1;
	hf[w_] := hf[w] = Sum[ng[[k+1]] * E^(I*k*w), {k, 0, nf}];
	sf[w_] := sf[w] = Product[hf[w/2^j], {j, 1, scaleres}];
	ss[k_ /; k==0] := 0;
	ss[k_ /; k<0] := ss[-k];
	ss[k_ /; OddQ[k]] := dn!/2^(2*dn)*f[-1]*sf[k*Pi];
	ss[k_ /; EvenQ[k]] := 1/2^dn*ss[k/2];
        N[Sum[Abs[N[ss[k]]]^2, {k, -calcres, calcres}]^(1/2)]
];

calcAppPow[g_List, dn_Integer] := Module[
	(* Calculates the approximation power constant as derived by
	   Michael Unser et al. Force to perform calculations with dn zeros
	   at z=-1 to avoid numerical errors for long filters. *)
	{ng = g / Plus@@g, scaleres = 6, calcres = 100, len = Length[g],
		even = If[EvenQ[Length[g]], 1, 0], hflen = Floor[Length[g]/2],
		h, dh, g2, f, nf, hf, sf, ss},
	h[z_] := Plus @@ MapThread[# z^(#2) &,
		{ng, Range[-hflen, hflen-even]}];
	g2 = CoefficientList[ Expand[ PolynomialQuotient[
		Expand[ z^hflen * h[z] ], ((1+z)/2)^dn, z]] , z];
	f[z_] := Plus@@MapThread[# z^(#2) &, {g2, Range[0, Length[g2] - 1]}];
	nf = Length[ng] - 1;
	hf[w_] := hf[w] = Sum[ng[[k+1]] * E^(I*k*w), {k, 0, nf}];
	sf[w_] := sf[w] = Product[hf[w/2^j], {j, 1, scaleres}];
	ss[k_ /; k==0] := 0;
	ss[k_ /; k<0] := ss[-k];
	ss[k_ /; OddQ[k]] := dn!/2^(2*dn)*f[-1]*sf[k*Pi];
	ss[k_ /; EvenQ[k]] := 1/2^dn*ss[k/2];
        N[Sum[Abs[N[ss[k]]]^2, {k, -calcres, calcres}]^(1/2)]
];

calcSFMom[g_List, n_Integer, sh_:0] := Module[
(* calculates n-th scaling function moment generated by filter g centered at x=sh *)
	{ng = g / Plus@@g, len = Length[g], hlen,
	mom, m, even = If[EvenQ[Length[g]], 1, 0]},
	hlen = IntegerPart[len/2];
	m[0] := 1;
	m[i_] := Plus@@MapThread[#^i #2&, {Range[-hlen+even-sh,hlen-sh], ng}];
	mom[0] := m[0];
	mom[p_] := 1/(2^p-1) Sum[Binomial[p,i] m[i] mom[p-i], {i, p}];
	mom[n]
];

(* closed form formula for odd degree biorthogonal coiflet reconstruction filters - from Tian. *)
bcwR[kin_, nin_] := Module[
	{k = (kin - 1)/2, n = (nin + 1)/2, c}, 
    If[EvenQ[nin],
    	Message["Even degree not implemented!"]; Return[$Failed]
    ]; 
    If[EvenQ[kin],
    	If[kin == 0, c = 1, c = 0], 
    	c = ((-1)^k*(2*n - 1))/(2^(4*n - 3)*(2*k + 1))*
         Binomial[2*n - 2, n - 1]*Binomial[2*n - 1, n + k];
    ];
    c
]; 

(* closed for formula for matching analysis filter of the same approximation order - from Tian. *)
bcwA[kin_, nin_] := Module[
	{c},
	If[EvenQ[nin],
		Message["Even degree not implemented!"]; Return[$Failed]
	];
	If[OddQ[kin],
		c = bcwR[kin,nin],
		c = If[kin==0, 2, 0] -
			Sum[bcwR[2*l+1,nin]*bcwR[2*l+1-kin,nin],{l,-nin,nin}]
	];
	c
];

coifRCoefs[n_ /; n>0] := Table[bcwR[k, n-1], {k, -n+1, n-1}];

coifACoefs[n_ /; n>0] := Table[bcwA[k, n-1], {k, -2*n+2, 2*n-2}];

visRecFilter[n_ /; n>0] := Module[
	{c, h},
	If[EvenQ[n],
		c = coifRCoefs[n],
		c = coifRCoefs[n+1];
		h[z_] := Plus @@ MapThread[# z^(#2) &, {c, Range[-n, n]}];
		c = CoefficientList[ Expand[ PolynomialQuotient[
			Expand[ z^n * h[z] ], (1+z)/2, z]] , z];
	];
	c
];

visAnaFilter[n_ /; n>0] := Module[
	{},
	If[EvenQ[n],
		coifACoefs[n],
		Message["Odd degrees not implemented yet!"]; Return[$Failed]
	]
];

End[]
EndPackage[]