(*  :Title:	Signal Processing Primitive Objects  *)

(*  :Authors:	Brian Evans, James McClellan  *)

(*
    :Summary:   Provide functions which are common in DSP but
		not standard in Mathematica.  Functions like
		Step and Impulse, as well as operators like
		Upsample and Shift.
 *)

(*  :Context:	SignalProcessing`Support`SigProc`  *)

(*  :PackageVersion:  2.6	*)

(*
    :Copyright:	Copyright 1989-1991 by Brian L. Evans
		Georgia Tech Research Corporation

	Permission to use, copy, modify, and distribute this software
	and its documentation for any purpose and without fee is
	hereby granted, provided that the above copyright notice
	appear in all copies and that both that copyright notice and
	this permission notice appear in supporting documentation,
	and that the name of the Georgia Tech Research Corporation,
	Georgia Tech, or Georgia Institute of Technology not be used
	in advertising or publicity pertaining to distribution of the
	software without specific, written prior permission.  Georgia
	Tech makes no representations about the suitability of this
	software for any purpose.  It is provided "as is" without
	express or implied warranty.
 *)

(*  :History:	signal processing, operators, functions, symbolic processing *)

(*  :Keywords:	*)

(*
    :Source:	     Myers, C.  {Signal Representation for Symbolic and
		Numeric Processing}.  MIT Ph.D. Thesis.  1986.
		     Covell, M.  {An Algorithm Design Environment}.
		MIT Ph. D. Thesis, 1989.
		     Bamberger, R. {The Directional Filter Bank:  A Multirate
		Filter Bank for the Directional Decompostion of Images},
		Georgia Tech Ph. D. Thesis, 1990.
 *)

(*  :Warning:	*)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*
    :Discussion:	This package introduces the functions (signals) and
		  operators (systems) common in signal processing but missing
		  in Mathematica.
			Thanks to John Mitchell at the University of Chicago
		  for writing patches for the Summation operator and the Delta
		  function.
*)

(*  :Functions: Aliasby			operator
		AliasSinc		alias
		ASinc			alias
		CConvolve		operator
		CircularShift		operator
		Convolve		operator
		Continuous		option value
		CPulse			function
		CStep			function
		Delta			function
		DeltaPlot		function
		DFT			operator
		DTFT			operator
		Discrete		option value
		Difference		operator
		Dirichlet		function
		DiscreteGraphics	routine
		Domain			option
		Downsample		operator
		DummyVariables		function
		Extent			function
		FT			operator
		GetDeltaFunctions	function
		Impulse			function
		Interleave		operator
		IntervalsToFunction	routine
		InvDFT			operator
		InvFT			operator
		InvL			operator
		InvZ			operator
		L			operator
		LineImpulse		function
		OperatorInOperExpr	routine
		OperatorName		routine
		ParametersInOperExpr	routine
		PolyphaseDownsample	operator
		PolyphaseUpsample	operator
		Periodic		operator
		Pulse			function
		RationalGCD		routine
		Rev			operator
		RewriteSummations	routine
		ScaleAxis		operator
		ScalingFactor		routine
		SequenceToFunction	routine
		SequencePlot		routine
		Shift			operator
		SignalPlot		function
		SignalsInOperExpr	routine
		Sinc			function
		SPSimplify		function
		Step			function
		Summation		operator
		TheFunction		routine
		ToContinuous		routine
		ToDiscrete		routine
		Unit			function
		Upsample		operator
		UpsampledFunction	routine
		UpsampleFactor		routine
		UpsampleSequence	routine
		Z			operator

	This package provides signal primitives (functions), system
	primitives (operators), and supporting routines.   Each new
	object is placed into the SPsignals, SPsystems, or SPfunctions
	list accordingly (at end of file).
	
	Each signal primitive (function) is specified as follows:

	(1) default arguments (if any),
	(2) definition (so that it behaves as a function),
	(3) multidimensional separability (if any),
	(4) automatic simplification rules,
	(5) manual simplification rules (augment Simplify),
	(6) output form rules (including TeX form rules), and
	(7) derivative and integration rules (if any)

	Each system primitive (operator) is specified as follows:

	(1) multidimensional separability (if any),
	(2) manual simplification rules (augment Simplify),
	(3) definition (so that it can be converted into a function form),
	(4) output form rules (including TeX form rules), and
	(5) specify which parameter is (are) the variable(s)
 *)


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ]


(*  B E G I N     P A C K A G E  *)

BeginPackage[ "SignalProcessing`Support`SigProc`",
	      "SignalProcessing`Support`SupCode`" ]


(*  U S A G E     I N F O R M A T I O N  *)

Aliasby::usage =
	"Aliasby[m,w][x] represents the aliasing of x, a continuous-frequency \
	signal, by replicating it every (2 Pi / m) in w and dividing the \
	result by m."

CConvolve::usage =
	"CConvolve[t][x1, x2, ...] represents the continuous convolution in t \
	of functions x1, x2, .... \
	CConvolve[{t1, t2, ...}][x1, x2, ...] represents the multidimensional \
	continuous convolution in t1, t2, ... of functions x1, x2, ..."

CircularShift::usage =
	"CircularShift[n0, N, n] represents the circular shifting of a \
	sequence of length N in index variable n by n0. \
	That is, the sequence is treated as a ring buffer of length N \
	so each sample at index n is shifted to the right by
	(n + n0) mod N places. \
	See also Shift and RotateRight."

Convolve::usage =
	"Convolve[n][x1, x2, ...] represents the discrete convolution in n of \
	functions x1, x2, .... \
	Convolve[{n1, n2, ...}][x1, x2, ...] represents the multidimensional \
	discrete convolution in n1, n2, ... of functions x1, x2, ..."

Continuous::usage =
	"Continuous is a possible value for a signal's domain."

CPulse::usage = 
	"CPulse[l,t] defines a pulse which begins at t=0 and ends at t = l. \
	The CPulse has value 1 within the range (0,l), \
	0 outside this range, and 1/2 at the points t=0 and t=l. \
	A continuous-pulse center at the origin is written as \
	CPulse[l, t + l/2] or Shift[-l/2, t][CPulse[l, t]]."

CStep::usage =
	"CStep[t], a.k.a. Unit[-1][t], is the unit step function \
	which is 1 for t > 0, 0 for t < 0, and 1/2 at t = 0. \
	It is commonly used for continuous expressions t. \
	See also Step and Unit."

Delta::usage =
	"Delta[expr] is the Dirac delta function. \
	The area under this functions is 1 but it only has value \
	at the origin. \
	That is, Integrate[ Delta[t] g[t], {t, t1, t2} ] is g[0] \
	if t1 <= 0 <= t2, 0 otherwise. \
	It differs from the Kronecker delta function Impulse[t]."

DeltaIntegrate::usage =
	"DeltaIntegrate[fun, {x, xlower, xupper}] and DeltaIntegrate[fun, x] \
	are special hooks for the built-in Integrate routine."

DeltaPlot::usage =
	"DeltaPlot is a object that can generate a plot of only Dirac delta \
	functions or a combination of Dirac delta functions and another plot. \
	In the first case, either four or six arguments are passed: \
	deltafuns, t, tmin, tmax, ymin, ymax.\
	In the second case, either five or seven arguments are passed: \
	deltafuns, t, tmin, tmax, plot, ymin, ymax. \
	In both cases, ymin and ymax are the optional arguments. \
	DeltaPlot returns a list of plots suitable for sending to Show. \
	Note that 1-D Delta functions are plotted as arrows."

DFT::usage =
	"DFT[L, n, k] is the forward discrete Fourier transform operator. \
	Here, L is the DFT length(s), n is the time-domain variable(s), and \
	k is the frequency-domain variable(s). \
	Multidimensional DFT's are supported. \
	Applying TheFunction to the DFT operator will invoke \
	the DFT rule base (object DFTransform) if loaded."

DTFT::usage =
	"DTFT[n, w] is the forward discrete-time Fourier transform operator. \
	Here, n is the time-domain variable(s), and w is the frequency-domain \
	variable(s). \
	Multidimensional DTFT's are supported. \
	Applying TheFunction to the DTFT operator will invoke the DTFT \
	rule base (object DTFTransform) if loaded."

Difference::usage =
	"Difference[k,n] represents the k-th backward difference with respect \
	to the variable n. \
	The first-order difference of f equals f[n] - f[n-1], \
	the second-order difference is f[n] - 2 f[n-1] + f[n-2], and so on. \
	Applying this operator to the sequence x[n] would be written as \
	Difference[k,n][ x[n] ]."

Dirichlet::usage =
	"Dirichlet[N, w] is the aliased sinc function. \
	which is defined as Sin[N w / 2] / ( N Sin [w / 2] ). \
	Aliases for Dirichlet are ASinc, AliasSinc, and AliasedSinc."

Discrete::usage =
	"Discrete is a possible value for a signal's domain."

DiscreteGraphics::usage =
	"DiscreteGraphics[{{x1,y1}, {x2,y2}, ... {xN,yN}}] will plot \
	the points as lines drawn from the x-axis to the y-value. \
	It takes two optional arguments:  a line thickness and \
	a point thickness. \
	Each is a real value from 0 to 1 indicating a width that \
	is a fraction of the size of the entire graph."

Domain::usage =
	"Domain is an option for many signal processing functions. \
	It is either Discrete or Continuous."

Downsample::usage =
	"Downsample[m,n] represents the downsampling operator. \
	Downsampling resamples a function at the points n = 0, m, 2m, 3m, ... \
	In N dimensions, n is a list of N indices and m is an N x N matrix. \
	Applying TheFunction to f[n] which is downsampled in one dimension, \
	(written as Downsample[m, n][ f[n] ]) will produce f[m n]."

DummyVariables::usage =
	"DummyVariables[dimension, variable] generates dimension dummy \
	variables with variable as the root. \
	For example, n can be generated by DummyVariables[1, n] or \
	DummyVariables[0, n]. \
	As another example, DummyVariables[3, t] returns {t1, t2, t3}. \
	In this cases, the generated variables have the same context \
	as variable. \
	DummyVariables[dimension, variable, context] can \
	be used to specify another context/package."

Extent1D::usage =
	"Extent1D[expression, var] returns an interval \
	{lower-limit, upper-limit} outside of which the function is zero. \
	This routine assumes that var is a real variable."

FT::usage =
	"FT[t, w] is the forward continuous Fourier transform operator. \
	Here, t is the time-domain variable(s) and w is the \
	frequency-domain variable(s). \
	Applying TheFunction to this operator will invoke the continuous-time \
	Fourier transform rule base (object CTFTransform) if loaded."

GetDeltaFunctions::usage =
	"GetDeltaFunctions[function, variable(s)] will extract \
	all delta functions from the argument function.  See also DeltaPlot."

Impulse::usage =
	"Impulse[n] is the unit Kronecker Delta function. \
	At n = 0, the function evaluates to 1. \
	Elsewhere, it evaluates to 0."

Interleave::usage =
	"Interleave[n][x0, x1, ...] interleaves samples of signals \
	x0, x1, ..., which are functions of n. \
	This is only applicable to discrete signals."

IntervalsToFunction::usage =
	"IntervalsToFunction[interlist, n, step, pulse] \
	will translate the list of functions and intervals (interlist) \
	into a signal expression as a function of n. \
	Each element in the interlist is a list of three elements: \
	<function_of_n>, <n_minus>, and <n_plus>. \
	The function <function_of_n> will be defined for \
	<n_minus> <= n <= <n_plus>. \
	Note that <n_minus> can be an integer or -Infinity and \
	that <n_plus> can be an integer or Infinity. \
	The default values are n = Global`n, \
	step = Step, and pulse = Pulse. \
	This function also works for when n represents a continuous variable \
	(step = CStep and pulse = CPulse)."

InvDFT::usage =
	"InvDFT[L, k, n] is the inverse discrete Fourier transform operator. \
	Here, L is the DFT length(s), n is the time-domain variable(s), and \
	k is the frequency-domain variable(s). \
	Multidimensional inverse DFT's are supported. \
	Applying TheFunction to the InvDFT operator will invoke the \
	inverse DFT rule base (object InvDFTransform) if loaded."

InvDTFT::usage =
	"InvDTFT[w, n] is the inverse discrete-time Fourier transform \
	operator. \
	Here, n is the time-domain variable(s), and w is the \
	frequency-domain variable(s). \
	Multidimensional inverse DTFT's are supported. \
	Applying TheFunction to the InvDTFT operator will invoke the \
	inverse DTFT rule base InvDTFTransform if loaded."

InvFT::usage =
	"InvFT[w, t] is the inverse discrete Fourier transform operator. \
	Here, t is the time-domain variable(s) and w is the \
	frequency-domain variable(s). \
	Multidimensional Fourier transforms are supported. \
	Applying TheFunction to this operator will invoke the inverse \
	continuous-time Fourier transform rule base InvCTFTransform if loaded."

InvL::usage =
	"InvL[s, t] is the inverse Laplace transform operator. \
	Applying TheFunction to the InvL operator will invoke the \
	bilateral inverse Laplace rule base InvLaPlace if loaded."

InvZ::usage =
	"InvZ[z, n] is the inverse z-transform operator. \
	Applying TheFunction to the InvZ operator will invoke the inverse \
	z-transform rule base InvZTransform if loaded."

L::usage =
	"L[t, s] is the forward Laplace transform operator. \
	Applying TheFunction to the L operator will invoke the bilateral \
	Laplace rule base if loaded."

LineImpulse::usage =
	"LineImpulse[nlist, coefflist] represents a line impulse where \
	nlist is a list of variables like {n1,n2,n3} and coefflist is \
	a corresponding list of coefficients like {1,2,3}. \
	For the previous lists, the line impulse is a set of impulse \
	along the line n1 = 2 n2 = 3 n3."

OperatorInOperExpr::usage =
	"OperatorInOperExpr[ operator_expression ] returns the full \
	operator in operator_expression. \
	For example, Shift[2,n] would be returned from \
	OperatorInOperExpr[ Shift[2,n][x[n]] ]."

OperatorName::usage =
	"OperatorName[ operator_expression ] returns the head of the \
	expression. \
	For example, OperatorName[ Shift[2,n][x[n]] ] will return Shift."

ParametersInOperExpr::usage =
	"ParametersInOperExpr[ operator_expression ] returns the parameters \
	of the operator in operator_expression. \
	For example, (2,n) would be returned from \
	ParametersInOperExpr[ Shift[2,n][x[n]] ]."

Periodic::usage =
	"Periodic[k,n][f] indicates that f is a function of n which is \
	periodic with period of k."

PolyphaseDownsample::usage =
	"PolyphaseDownsample[m, h, n][x] is equivalent to downsampling \
	x by m (w/r to n) and then filtering (w/r to n) by a polyphase \
	form of the FIR filter h. \
	Hence, the filtering is carried out at the lower sampling rate."

PolyphaseUpsample::usage =
	"PolyphaseUpsample[l, h, n][x] is equivalent to filtering (w/r to n) \
	x by a polyphase form of the FIR filter h and then upsampling the \
	result by l (w/r to n). \
	Hence, the filtering is carried out at the lower sampling rate."

Pulse::usage =
	"Pulse[len] and Pulse[len, n] define a pulse which begins at \
	n = 0 and ends at n = len - 1. \
	A discrete pulse centered at the origin is written as \
	Pulse[len, n + (len - 1)/2] or Shift[-(len - 1)/2, n][Pulse[len, n]]. \
	Sometime, a pulse will be automatically simplified; \
	e.g., Pulse[0, n] is zero and Pulse[1, n] is really Impulse[n]."

RationalGCD::usage =
	"RationalGCD[e1, e2, ...] returns the greatest common divisor of the \
	numerators of elements ei divided by the greatest common divisor of \
	the denominators of elements ei. \
	Every element must be a polynomial or a number. \
	Each real-valued and complex-valued element is rationalized \
	before GCD's are computed."

Rev::usage =
	"Rev[n][x] represents the operation of reversing the direction of x, \
	with respect to the variables(s) n."

RewriteSummations::usage =
	"RewriteSummations[f, t, tlower, tupper] rewrites Summation \
	and Periodic operators that appear in f so that the domain \
	of the new expression includes t in [tlower, tupper]. \
	It relies on the Extent1D function to compute the domain \
	of the arguments to the Periodic and Summation operators."

ScaleAxis::usage =
	"ScaleAxis[l,w][x] represents the compression by a factor of l (an \
	integer) of the continuous w axis of x, which is usually a \
	continuous-frequency signal."

ScalingFactor::usage =
	"ScalingFactor[expr, variable] returns the greatest common divisor \
	of all coefficients of powers of z in expr. \
	The powers of z that are functions of z are not considered. \
	For example, ScalingFactor[ a z + a^2 z^2, z ] returns a. \
	This function enables the implementation of the similarity \
	property for transforms."

SequencePlot::usage =
	"SequencePlot[f, {n, start, end}] plots real-valued 1-D sequences \
	and has the same options as Show. \
	You will have to apply Re or Im to f to plot the real or \
	imaginary values of the sequence, respectively. \
	SequencePlot[f, {n1, start1, end1}, {n2, start2, end2}] plots \
	2-D sequences and supports the same options as Plot3D."

SequenceToFunction::usage =
	"SequenceToFunction[seq], SequenceToFunction[seq, n], and \
	SequenceToFunction[seq, n, noffset] return the sequence seq, \
	{x[noffset], ..., x[noffset + N - 1]}, in symbolic form as a sum of \
	impulse functions. \
	For example, SequenceToFunction[{1, 2, 3}] -> \
	Impulse[n] + 2 Impulse[n - 1] + 3 Impulse[n - 2]. \
	For multidimensional sequences, the argument n must be a list of \
	variables of appropriate length. \
	For example, SequenceToFunction[{{1,2},{3,4}}, {n1, n2}] -> \
	Impulse[n1,n2] + 2 Impulse[n1,n2-1] + 3 Impulse[n1-1,n2] + \
	4 Impulse[n1-1, n2-2]."

Shift::usage =
	"Shift[m,v][x] represents the shifting of the expression x by m \
	samples to the right in the v direction, i.e. x[v] --> x[v - m]."

SignalPlot::usage =
	"SignalPlot[f, {t, start, end}] will plot f(t) as an one-dimensional, \
	continuous-time function. \
	It will show the real part as solid lines, \
	and the imaginary part as dashed lines. \
	Delta functions are plotted as upward pointing arrows. \
	SignalPlot[f, {t1, start1, end1}, {t2, start2, end2}] treats f \
	as a function of two variables t1 and t2. \
	SignalPlot supports the same options as Plot for 1-D signals \
	(functions) and Plot3D for 2-D signals (functions)."

SignalsInOperExpr::usage =
	"SignalsInOperExpr[ operator_expression ] returns the signal(s) \
	of the operator in operator_expression. \
	For example, { x[n] } would be returned from \
	SignalsInOperExpr[ Shift[2,n][x[n]] ]. \
	Another version SignalsInOperExpr[ operator_expression, operator ] \
	behaves as the first version if operator is the head of \
	operator_expression; otherwise, operator_expression is returned."

Sinc::usage =
	"Sinc[t] is the unit sinc function which is the ratio of \
	Sin[t] over t.  Note that Sinc[0] = 1."

SPSimplify::usage =
	"SPSimplify[e] will simplify expression e according to Mathematica's \
	simplification rules and the rules in SPSimplificationRules. \
	At each iteration, Simplify is called and SPSimplificationRules \
	is applied to the result. \
	Unlike Simplify, SPSimplify supports a Dialogue option. \
	You can specify which symbols to treat as real variables by using \
	the Variables option."

Step::usage =
	"Step[n] is the unit step function which is 1 for n >= 0 \
	and 0 otherwise. \
	It is commonly used for discrete n. \
	This function differs from CStep only at the origin."

Summation::usage =
	"Summation[i][f], Summation[i, iend][f], and \
	Summation[i, istart, iend, inc][f] represents the summation of f \
	with respect to variable i:  i = (istart) to (iend) step (inc). \
	The Summation operator is an abstraction of Mathematica's Sum object. \
	Applying the object TheFunction to a Summation operator will invoke \
	the Sum object if istart, iend, and inc are numbers."

TheFunction::usage =
	"TheFunction[object] returns the mathematical function embedded \
	in object. \
	For list and transform objects, this function returns \
	the first element after the head of the object; otherwise, the \
	contents of the slot TheFunction are returned. \
	For signal processing operations, \
	the equivalent mathematical operation is returned."

ToContinuous::usage =
	"ToContinuous[expr] replaces any discrete operator or signal \
	in expr with its continuous equivalent."

ToDiscrete::usage =
	"ToDiscrete[expr] replaces any continuous operator or signal \
	in expr with its discrete equivalent."

Unit::usage =
	"Unit[n][t] is the nth order continuous unit function in t. \
	Unit[0][t] is the Dirac delta function (see Delta). \
	Unit[-1][t] is the continuous step function (see CStep)."

Upsample::usage =
	"Upsample[k,n][f] represents the upsampling operation on f, \
	which is either an expression of n or a sequence (list) of values. \
	Upsampling inserts k-1 zeroes between every two samples. \
	In N dimensions, n is a list of N indices and k is an N x N matrix."

UpsampledFunction::usage =
	"UpsampledFunction[m, fill, n, upf] represents the upsampling of f. \
	In one-dimension, the upsampled form of f, upf, \
	is obtained by evaluating n -> n / m."

UpsampleFactor::usage =
	"UpsampleFactor[f, n] returns the upsampling index which can be \
	a negative or positive number. \
	If this routine returns 1 or -1, then no upsampling has occurred; \
	if it returns 0, then f is constant with respect to n."

UpsampleSequence::usage =
	"UpsampleSequence[seq,m,fill] inserts m-1 fill values between \
	every two elements of seq."

Z::usage =
	"Z[n, z] is the forward z-transform operator. \
	Applying TheFunction to the Z operator will invoke the z-transform \
	rule base ZTransform if loaded."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin["`Private`"]


(*  C O N S T A N T S  *)

axesdefaultvalue = If [ TrueQ[$VersionNumber >= 2.0], True, Automatic]


(*  M E S S A G E S  *)

CircularShift::badlength =
	"The sequence length parameter of the CircularShift \
	operator must be non-zero."

DeltaPlot::unresolved = "These symbols were assigned a value of 1: ``."

Downsample::convert =
    "Converting downsampling vector into a diagonal matrix: ``."
Downsample::badmD =
    "Downsampling factor has different dimensions than downsampled variables."

Pulse::notsym = "Last argument to Pulse is not a symbol."

Shift::badmD =
    "Shift vector has different dimensions than the specified variables."

SequencePlot::nolist =
    "It does not make sense to plot a list of 2-D sequences."

SignalPlot::nolist =
    "It does not make sense to plot a list of 2-D functions."

Summation::extent =
    "Could not compute which part of the summation should be plotted."

Upsample::convert =
    "Converting upsampling vector into a diagonal matrix: ``."
Upsample::badmD =
    "Upsampling factor has different dimensions than upsampled variables."



(*  L O C A L     F U N C T I O N S  *)

(*  breakup  *)
breakup[a_, z_] := a /; FreeQ[a, z]
breakup[a_Times, z_] := Select[a, breakupq[z]]
breakup[x_, z_] := 1

breakupq[z_][a_] := FreeQ[a, z]

(*  gcdfilter --  rationalize floating point numbers		*)
gcdfilter[x_Complex] := ToCollection[ gcdfilter[Re[x]], gcdfilter[Im[x]] ]
gcdfilter[x_Real] := Rationalize[x]
gcdfilter[x_] := x

(*  informatargs and informat  *)
SetAttributes[informatargs, {Listable}]
informatargs[m_, n_Symbol] := ToString[ StringForm["``  in  ``", m, n] ]
informatargs[m_, n_String] := ToString[ StringForm["``  in  ``", m, n] ]
informatargs[m_, n_Pattern] := ToString[ StringForm["``  in  ``", m, n] ]

informat[op_, m_, n_] :=
Format[ SequenceForm[op, Subscript[ ColumnForm[ToList[informatargs[m, n]]] ]] ]

(* intdivide --  kludge for Quotient primitive which has a problem  *)
(*               under Mma 1.2 when the second argument is negative *)
If [ TrueQ[ $VersionNumber >= 2.0 ],
     intdivide[n_, m_] := Quotient[n, m],
     intdivide[n_, m_] := Sign[m] Quotient[n, m Sign[m]] ]

(*  inTeXformat  *)
inTeXformat[op_, m_, n_Symbol] := Subscripted[ op[m, n] ]
inTeXformat[op_, m_, n_List] := Subscripted[ op[TeXfix[m], TeXfix[n]] ]

(*  TeXfix  *)
TeXfix[a_?VectorQ] := MatrixForm[ Transpose[{a}] ]
TeXfix[a_] := DelimitedMatrix[ a ]

(*  DelimitedMatrix  *)
DelimitedMatrix/: Format[ DelimitedMatrix[a_], TeXForm ] :=
	StringForm["\\left[ `` \\right]", MatrixForm[a] ]


(*  D E T E R M I N A T I O N     O F      V A R I A B L E S	*)

DummyVariables[i_, v_] :=
	If [ i < 2,
	     v,
	     Table[GenerateSymbol[v, index], {index, 1, i}] ] /;
	IntegerQ[i] && ( i >= 0 )

DummyVariables[i_, v_, context_] :=
	If [ i < 2,
	     GenerateSymbol[v, " ", context],
	     Table[GenerateSymbol[v, index, context], {index, 1, i}] ] /;
	IntegerQ[i] && ( i >= 0 )


(*  F U N C T I O N S  *)

(*  Aliasby  *)
Aliasby/: Aliasby[1, w_][x_] := x

Aliasby/: Format[ Aliasby[m_, w_] ] := informat[ Aliasby, m, w ]
Aliasby/: Format[ Aliasby[m_, w_], TeXForm ] := inTeXformat[ Aliasby, m, w ]

Aliasby/: TheFunction[ Aliasby[m_, w_][f_] ] :=
		1/m Periodic[2 Pi / m, w][ TheFunction[f] ]

Aliasby/: GetOperatorVariables[ Aliasby[m_, w_] ] := w

(*  CConvolve --  continuous convolution  *)
CConvolve/: CConvolve[{t_Symbol}] := CConvolve[t]
CConvolve/: CConvolve[t1_, trest__] := CConvolve[{t1, trest}]

CConvolve/: CConvolve[t_Symbol] [x_, c_. Delta[t_ + m_.]] :=
	(c x) /. t -> (-m) /; FreeQ[m, t]
CConvolve/: CConvolve[t_Symbol] [c_. Delta[t_ + m_.], x_] :=
	(c x) /. t -> (-m) /; FreeQ[m, t]

CConvolve/: Extent1D[ CConvolve[t_][x__], t_ ] :=
	Apply[ Plus, Map[ Extent1D[t], {x} ] ]

CConvolve/: CConvolve[tlist_List] [x_, c_. Delta[t_ + m_.]] :=
	CConvolve[ Complement[tlist, {t}] ][x /. t -> (-m), c] /;
	MemberQ[tlist, t] && FreeQ[m, t]
CConvolve/: CConvolve[tlist_List] [c_. Delta[t_ + m_.], x_] :=
	CConvolve[ Complement[tlist, {t}] ][c, x /. t -> (-m)] /;
	MemberQ[tlist, t] && FreeQ[m, t]

CConvolve/: CConvolve[t_][x_, z_?ZeroQ] := z
CConvolve/: CConvolve[t_][z_?ZeroQ, x_] := z

CConvolve/: Format[CConvolve[t_]] :=
		Format[ Subscripted[CConvolve[t]] ]
CConvolve/: Format[CConvolve[t_][x1_, x2_], TeXForm] :=
		StringForm[ "`` \\star_{``} ``", Format[x1], t, Format[x2] ]
CConvolve/: Format[CConvolve[t_][x1_, x2_, rest__], TeXForm] :=
		Block [	{format, i, length, operstr, xlist},
			xlist = {x1, x2, rest};
			length = Length[xlist];
			operstr = ToString[StringForm["\\star_{``}", t]];
			str = StringJoin[" ", operstr, " ``"];
			format = Apply[ StringJoin,
					Table[ str, {i, 1, length} ] ];
			format = StringJoin["`` ", operstr, " ``", format];
			Apply[StringForm, {format, Map[Format, xlist]}] ]

(*  CircularShift  *)
CircularShift/: CircularShift[n0_, N_] := CircularShift[n0, N, Global`n]

CircularShift/: CircularShift[n0_, 0, n_][x_] :=
	MyMessage[CircularShift::badlength, x]

CircularShift/: Extent1D[ CircularShift[n0_, N_, n_][f_], n_ ] := { 0, N-1 }

CircularShift/: TheFunction[ CircularShift[n0_Integer, N_Integer, n_][f_List] ] :=
	RotateRight[ Take[f, N], n0 ]

CircularShift/: TheFunction[ CircularShift[n0_, N_, n_][f_] ] :=
	TheFunction[f Pulse[N, n]] /. ( n -> Mod[n + n0, N] )

CircularShift/: N[ CircularShift[n0_, N_, n_][f_] ] :=
	N[ TheFunction[CircularShift[n0, N, n][f]] /. ( n -> Mod[n + n0, N] ) ]

CircularShift/: Format[ CircularShift[n0_, N_, n_] ] :=
	informat[CircularShift, Mod[n + n0, N], n]
CircularShift/: Format[ CircularShift[m_,n_], TeXForm ] :=
	inTeXformat[CircularShift, Mod[n + n0, N], n]

(*  CPulse  *)
CPulse[l_] := CPulse[l, Global`t]			(* default args     *)

CPulse/: CPulse[l_, t_] := 1	 /; N[0 < t < l]	(* definition	    *)
CPulse/: CPulse[l_, t_] := 0	 /; N[(t < 0) || (t > l)]
CPulse/: CPulse[l_, t_] := 1/2   /; N[(t == 0) || (t == l)]
CPulse/: CPulse[Infinity, t_] := CStep[t]
CPulse/: TheFunction[CPulse[l_, t_]] := CStep[t] - CStep[t - l]
CPulse/: Extent1D[ CPulse[l_, a_. t_ + b_.], t_ ] :=
	Which [ Positive[a], {-b/a, (l-b)/a},
		Negative[a], {(l-b)/a, -b/a},
		True,	     Extent1DOrder[-b/a, (l-b)/a] ] /;
	FreeQ[{a,b,l}, t]

CPulse/: Format[ CPulse[l_, n_] ] :=			(* output forms     *)
		Format[ CPulse Subscript[l] [n] ]

CPulse/: Format[ CPulse, TeXForm ] := "\\sqcap"	
CPulse/: Format[ CPulse[l_, n_], TeXForm ] :=
		StringForm[ "\\sqcap_{``} (``)", l, n ]

Derivative[0,1][CPulse] :=				(* derivative rule *)
	( Delta[#2] + Delta[#2 - #1] )&

Unprotect[Integrate]					(* integration rules *)
Integrate/: Integrate[f_. CPulse[L_, a_. x_ + b_.] + c_.,
		vars1___, x_Symbol, vars2___] :=
	1/a CPulse[a x + b] ( Integrate[f, vars1, x, vars2] /. x -> a x + b ) +
	  Integrate[c, vars1, x, vars2] /;
	FreeQ[{a, b, L}, x]
Integrate/: Integrate[expr_. CPulse[L_, a_. x_ + b_.] + c_.,
			vars1___, {x_Symbol, x1_, x2_}, vars2___] :=
	Sign[x2 - x1] *
	cpulseintersect[ Min[x1, x2], Max[x1, x2],
			 Min[-b/a, (L - b)/a], Max[-b/a, (L - b)/a] ] *
	Integrate[expr,
		  vars1,
		  {x, Max[Min[x1,x2], Min[-b/a, (L - b)/a]],
		      Min[Max[x1,x2], Max[-b/a, (L - b)/a]]},
	          vars2] +
	Integrate[c, vars1, {x, x1, x2}, vars2] /;
	FreeQ[{a, b, L}, x]
Protect[Integrate]

cpulseintersect[x1_, x2_, x3_, x4_] := Step[x4 - x1] Step[x2 - x3]

(* CStep  *)
CStep/: CStep[Infinity] := 1				(* definition 	    *)
CStep/: CStep[t_] := 1		/; N[ t > 0 ]
CStep/: CStep[t_] := 0		/; N[ t < 0 ]
CStep/: CStep[0] := 1/2
CStep/: CStep[0.] := 0.5
CStep/: CStep[-Infinity] := 0

CStep/: Extent1D[ CStep[a_. t_ + b_.], t_ ] :=
	Which [ Positive[a], {-b/a, Infinity},
		Negative[a], {-Infinity, -b/a},
		True,	     Extent1DOrder[-b/a, Sign[a] Infinity] ] /;
	FreeQ[{a,b,l}, t]

CStep/: CStep[t_, args__] := CStep[t] CStep[args]	(* multidimensional *)
CStep/: CStep[List[args__]] := CStep[args]		(* separability     *)

CStep/: CStep[n_Symbol + m_.] CStep[n_Symbol  + k_.] :=	(* automatic	    *)
	CStep[n + Min[m,k]] /;				(* simplification   *)
	NumberQ[m] && NumberQ[k]			(* rules	    *)
CStep/: CStep[n_Symbol + m_.] CStep[-n_Symbol + k_.] :=
	CPulse[k + m, n + m] /;
	N[ -m <= k ]
CStep/: CStep[n_Symbol + m_.] CStep[-n_Symbol + k_.] := 0 /;
	N[ -m > k ]

CStep/: Simplify[CStep[a_ n_Symbol + b_.]] :=		(* man. simp. rules *)
	CStep[n + b/a] /;
	MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ] 
CStep/: Simplify[CStep[- a_ n_Symbol + b_.]] :=
	CStep[-n + b/a] /;
	MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ]

CStep/: Format[ CStep, TeXForm ] := "u_{-1}"		(* output form     *)

Derivative[1][CStep] := Delta				(* derivative rule *)

Unprotect[Integrate]					(* integration rules *)
Integrate/: Integrate[f_. CStep[a_. x_ + b_.] + c_.,
		vars1___, x_Symbol, vars2___] :=
	1/a CStep[a x + b] ( Integrate[f, vars1, x, vars2] /. x -> a x + b ) +
	  Integrate[c, vars1, x, vars2] /;
	FreeQ[{a, b}, x]
Integrate/: Integrate[expr_. CStep[a_. x_ + b_.] + c_.,
			vars1___, {x_Symbol, x1_, x2_}, vars2___] :=
	Sign[x2 - x1] *
	cstepintersect[x1, x2, a, b] *
	Integrate[expr,
		  vars1,
		  {x, csteplbound[x1, x2, a, b], cstepubound[x1, x2, a, b]},
	          vars2] +
	Integrate[c, vars1, {x, x1, x2}, vars2] /;
	FreeQ[{a, b}, x]
Protect[Integrate]

cstepintersect[x1_, x2_, a_, b_] :=
	If [ Positive[a],
	     Step[Max[x1, x2] + b/a],
             Step[-b/a - Min[x1,x2]] ]

csadj[x_] := x //. SignalProcessing`Support`SupCode`Private`MinMaxRules

csteplbound[x1_, x2_, a_, b_] := csadj[Max[Min[x1, x2], -b/a]] /; Positive[a]
csteplbound[x1_, x2_, a_, b_] := csadj[Min[x1, x2]] /; Negative[a]
cstepubound[x1_, x2_, a_, b_] := csadj[Max[x1, x2]] /; Positive[a]
cstepubound[x1_, x2_, a_, b_] := csadj[Min[Max[x1, x2], -b/a]] /; Negative[a]

Format[csteplbound[x1_, x2_, a_, b_]] :=
	SequenceForm["LeftEndpointOf[", {x1,x2}, " and step at ", -b/a, "]"]
Format[cstepubound[x1_, x2_, a_, b_]] :=
	SequenceForm["RightEndpointOf[", {x1,x2}, " and step at ", -b/a, "]"]

(*  Convolve --  discrete convolution  *)
Convolve/: Convolve[{n_Symbol}] := Convolve[n]
Convolve/: Convolve[n1_, nrest__] := Convolve[{n1, nrest}]

Convolve/: Convolve[n_] [x_, c_. Impulse[n_ + m_.]] :=
	(c x) /. n -> (-m) /; FreeQ[m, t]
Convolve/: Convolve[n_] [c_. Impulse[n_ + m_.], x_] :=
	(c x) /. n -> (-m) /; FreeQ[m, t]

Convolve/: Extent1D[ Convolve[n_][x__], n_ ] :=
	Apply[ Plus, Map[ Extent1D[n], {x} ] ]

Convolve/: Convolve[nlist_List] [x_, c_. Impulse[n_ + m_.]] :=
	Convolve[ Complement[nlist, {n}] ][x /. n -> (-m), c] /;
	MemberQ[nlist, n] && FreeQ[m, n]
Convolve/: Convolve[nlist_List] [c_. Impulse[n_ + m_.], x_] :=
	Convolve[ Complement[nlist, {n}] ][c, x /. n -> (-m)] /;
	MemberQ[nlist, n] && FreeQ[m, n]

Convolve/: Convolve[n_][x_, z_?ZeroQ] := z
Convolve/: Convolve[n_][z_?ZeroQ, x_] := z

	(*  Formatting the convolution operator *)

Convolve/: Format[Convolve[n_]] :=
		Format[ Subscripted[Convolve[n]] ]
Convolve/: Format[Convolve[n_][x1_, x2_], TeXForm] :=
		StringForm[ "`` \\star_{``} ``", Format[x1], n, Format[x2] ]
Convolve/: Format[Convolve[n_][x1_, x2_, rest__], TeXForm] :=
		Block [	{format, i, length, operstr, xlist},
			xlist = {x1, x2, rest};
			length = Length[xlist];
			operstr = ToString[StringForm["\\star_{``}", n]];
			str = StringJoin[" ", operstr, " ``"];
			format = Apply[ StringJoin,
					Table[ str, {i, 1, length} ] ];
			format = StringJoin["`` ", operstr, " ``", format];
			Apply[StringForm, {format, Map[Format, xlist]}] ]

(*  Delta  *)
Delta/: Delta[Infinity] := 0				(* definition	     *)
Delta/: Delta[t_] := 0		/; N[ t != 0 ]
Delta/: Delta[-Infinity] := 0

Delta/: Delta[a_. x_Symbol + b_.] Delta[c_. x_ + d_.] := 0 /; ( b c != a d )
Delta/: Delta[a_. x_Symbol + b_.] Delta[a_. x_ + d_.] := 0 /; ( b != d )

Delta/: Extent1D[ Delta[a_. t_ + b_.], t_ ] := { -b/a, -b/a } /;
	FreeQ[{a,b}, t]

Delta/: Delta[t_, args__] := Delta[t] Delta[args]	(* multidimensional  *)
Delta/: Delta[List[args__]] := Delta[args]		(* separability	     *)

Delta/: Format[ Delta, TeXForm ] := "\\delta"		(* output form       *)

Derivative[1][Delta] := Unit[1]				(* derivative rule   *)

	(* The following Delta rules were written by John Mitchell *)
	(* and modified by Brian Evans.				   *)
getivars[ x_Symbol ] := x
getivars[ {x_Symbol, xl_, xu_} ] := x
Unprotect[Integrate]					(* integration rules *)
Integrate/: Integrate[a_. + expr_. Delta[term_], vars_] :=
	Integrate[a, vars] + DeltaIntegrate[expr Delta[term], vars] /;
	! FreeQ[ term, getivars[vars] ]
Integrate/: Integrate[a_. + expr_. Delta[term_], vars1___, vars_, vars2___] :=
	Integrate[a, vars1, vars, vars2] + 
	  Integrate[ DeltaIntegrate[expr Delta[term], vars], vars1, vars2 ] /;
	! FreeQ[ term, getivars[vars] ]
Protect[Integrate]

(*  DeltaIntegrate  *)
DeltaIntegrate[expr_. Delta[a_. x_ + b_.], x_Symbol] :=
	Limit[expr, x -> -b/a] CStep[x + b/a] / Abs[a] /;
	FreeQ[{a,b}, x] && FreeQ[expr, Delta]
DeltaIntegrate[expr_. Delta[a_. x_ + b_.], {x_Symbol, x1_, x1_}] := 0
DeltaIntegrate[expr_. Delta[a_ x_ + b_.], {x_Symbol, x1_, x2_}] :=
	If [ N[x1 > x2],
	     -1 DeltaIntegrate[expr Delta[a x + b], {x, x2, x1}],
	     Integrate[expr Delta[x + b/a] / Abs[a], {x, x1, x2}],
	     Integrate[expr Delta[x + b/a] / Abs[a], {x, x1, x2}] ] /;
	FreeQ[{a,b}, x] && FreeQ[expr, Delta]
DeltaIntegrate[expr_. Delta[x_ + c_.], {x_Symbol, x1_, x2_}] :=
	If [ N[x1 > x2],
	     -1 DeltaIntegrate[expr Delta[x + c], {x, x2, x1}],
	     Limit[expr, x -> -c] CStep[-c - x1] CStep[c + x2],
	     Limit[expr, x -> -c] CStep[-c - x1] CStep[c + x2] ] /;
	FreeQ[c, x] && FreeQ[expr, Delta]
		(* CStep[-c - x1] CStep[c + x2] <--> x1 <= -c <= x2 *)
DeltaIntegrate[Delta[a_. x_ + b_.] Delta[c_.x_ + d_.],
		      {x_Symbol, x1_, x2_}] :=
	If [ N[x1 > x2],
	     -1 DeltaIntegrate[Delta[a x + b] Delta[c x + d], {x, x2, x1}],
	     (1/Abs[a c]) Delta[-b/a + d/c] *
		CStep[-(b/a) - x1] CStep[(b/a) + x2],
	     (1/Abs[a c]) Delta[-b/a + d/c] *
		CStep[-(b/a) - x1] CStep[(b/a) + x2] ] /;
	FreeQ[{a,b,c,d}, x]
DeltaIntegrate[Delta[a_. x_ + b_.]^2, {x_Symbol, x1_, x2_}] :=
	If [ N[x1 > x2],
	     -1 DeltaIntegrate[Delta[a x + b]^2, {x, x2, x1}],
	     (1/a^2) Delta[0] CStep[-(b/a) - x1] CStep[(b/a) + x2],
	     (1/a^2) Delta[0] CStep[-(b/a) - x1] CStep[(b/a) + x2] ] /;
	FreeQ[{a,b}, x]

(*  DeltaPlot--  plot Delta functions as special case of continuous plotting  *)
dpcleanup[deltafuns_, t_, tmin_, tmax_] :=
	Block [	{delfuns, varlist = {}},
		delfuns = deltafuns /.
			    ( d_ Summation[index_, istart_, iend_, inc_][c_] :>
			      Summation[index, istart, iend, inc][c d] );
		delfuns = delfuns /.
		{ Summation[index_, istart_, iend_, inc_]
			   [c_. Delta[b_. t + a_.]] :>
  		  Block [ {aa = a, adjfun, bb = b, cc = c, defaults, f1, f2,
			   i, indexfunoft, minindex, maxindex, rules, varlist},
			  varlist = Select[GetVariables[{a, b, c}],
				   	   ((! SameQ[#1, t]) &&
					    (! SameQ[#1, index]))&];
			  If [ ! EmptyQ[varlist],
			       defaults = Map[(#1 -> 1)&, varlist];
			       aa = a /. defaults;
			       bb = b /. defaults;
			       cc = c /. defaults ];
			  rules = Solve[ bb t + aa == 0, index ];
			  indexfunoft = Replace[index, First[rules] ];
			  f1 = N[ indexfunoft /. t -> tmin ];
			  f2 = N[ indexfunoft /. t -> tmax ];
			  minindex = Max[istart, inc Floor[Min[f1, f2] / inc]];
			  maxindex = Min[iend, inc Ceiling[Max[f1, f2] / inc]];
			  adjfun = ((cc /. t -> - aa / bb) Delta[bb t + aa]) /.
				   index -> i;
			  Apply[Plus,
				Table[adjfun, {i, minindex, maxindex, inc}]]] };

		varlist = Sort[ varlist ~Join~
				Select[ GetVariables[delfuns],
					(! SameQ[#1, t])& ] ];
		If [ ! EmptyQ[varlist],
		     Message[ DeltaPlot::unresolved, varlist ];
		     delfuns = delfuns /. Map[(#1 -> 1)&, varlist] ];

		delfuns ]

getlocation[ c_. Delta[a_. t_ + b_.], t_ ] := N[ - b / a ]

DeltaPlot[deltafuns_, t_, start_, end_] :=
	DeltaPlot[deltafuns, t, start, end, 0, 1]

DeltaPlot[deltafuns_, t_, start_, end_, plot_] := plot /;
	FreeQ[deltafuns, t]

    (* Use PlotRange as function to find range of plot [see MJ v1 #4 p47] *)
DeltaPlot[deltafuns_, t_, start_, end_, plot_] :=
	Block [	{xmax, xmin, ylist, ymax, ymin},
		{{xmin, xmax}, {ymin, ymax}} = PlotRange[plot];
		DeltaPlot[deltafuns, t, start, end, plot, ymin, ymax] ]

DeltaPlot[deltafuns_, t_, start_, end_, min_, max_] := NullPlot /;
	FreeQ[deltafuns, t]

	(* every version of DeltaPlot calls this one *)
DeltaPlot[deltafuns_, t_, start_, end_, min_, max_] :=
	Block [	{bottom, delfuns, height, length,
		 locations, nend, nstart, points, width},
		nstart = N[start];
		nend = N[end];
		delfuns = Expand[ dpcleanup[deltafuns, t, nstart, nend] ];
		height = N[max - min];
		width = nend - nstart;
		points = If [ SameQ[Head[delfuns], Plus],
			      Map[getlocation[#1, t]&, Apply[List, delfuns]],
			      { getlocation[delfuns, t] } ];
		locations = Select[ points, N[LessEqual[nstart, #1, nend]]& ];
		If [ min <= 0,
		     bottom = 0; length = max,
		     bottom = min; length = height ];
		length = Min[height, max];
		If [ SameQ[locations, {}],
		     NullPlot,
		     Map [ Arrow2D[{#1, bottom}, width, height, length]&,
			   locations ] ] ]

DeltaPlot[deltafuns_, t_, start_, end_, plot_, min_, max_] :=
	Block [	{deltaplot, ymax = max, ymin = min},
		If [ ymax == ymin,
		     Which [ ymin > 0,  ymin = -ymax,
			     ymin == 0, ymax = 1,
			     ymin < 0,  ymax = -ymin ] ];
		deltaplot = DeltaPlot[deltafuns, t, start, end, ymin, ymax];
		If [ SameQ[deltaplot, NullPlot],
		     plot,
		     Prepend[ deltaplot, plot ] ] ]

(*  Difference  *)
Difference/: TheFunction[Difference[k_, n_][f_]] :=
	Block [	{jj, newf},
		newf = TheFunction[f] /. n ~ReplaceWith~ (n - jj);
		Apply[ Plus,
		       Table[(-1)^jj Binomial[k,jj] newf, {jj, 0, k}] ] ] /;
	IntegerQ[k] && k > 0

Difference/: Extent1D[ Difference[k_, n_][f_], n_ ] := Extent1D[f, n]

Difference/: N[Difference[1, n_][f_]] := f - ( f /. n ~ReplaceWith~ (n - 1) )
Difference/: N[Difference[k_,n_][f_]] :=
	Apply[Plus,
	      Table[(-1)^jj Binomial[k,jj] f /. n ~ReplaceWith~ (n - jj),
		    {jj, 0, k}] ] /;
	IntegerQ[k] && k > 0

Difference/: Format[ Difference[k_, n_] ] := informat[Difference, k, n]

Difference/: Format[ Difference, TeXForm ] := "\\Delta"
Difference/: Format[ Difference[k_,n_], TeXForm ] :=
		inTeXformat[ Difference, k, n ]

Difference/: GetOperatorVariables[ Difference[k_, n_] ] := n

(*  Dirichlet  *)
Dirichlet/: Dirichlet[len_, 0] := len			(* definition	*)
Dirichlet/: Dirichlet[len_, k_?EvenQ Pi] := len (-1)^(len - 1)
Dirichlet/: Dirichlet[len_Integer, w_?NumberQ] :=
		Sin[len w / 2] / ( len Sin[w / 2] )

Dirichlet/: TheFunction[Dirichlet[len_, w_]] :=
		Sin[len w / 2] / ( len Sin[w / 2] )

Dirichlet/: N[ Dirichlet[len_, w_] ] := N [ TheFunction[Dirichlet[len, w]] ]

Dirichlet/: Format[ Dirichlet[l_, n_] ] :=		(* output forms     *)
		Format[ Dirichlet Subscript[l] [n] ]

(*  DiscreteGraphics  *)
oneline[xycoord_] :=
	ToCollection[ Line[{{xycoord[[1]], 0}, xycoord}], Point[xycoord] ]
DiscreteGraphics[coordlist_List, lthick_:0.008, pthick_:0.024 ] :=
	Graphics[ Join[ { Thickness[lthick], PointSize[pthick] },
			N [ Map[oneline, coordlist] ] ] ]

(*  Downsample [Myers, 219] [Bamberger, ch. 2]  *)

	(* Error checking of downsample parameters  *)

Downsample/: Downsample[m_List, n_List][f_] :=
	MyMessage[ Downsample::convert,
		   Downsample[DiagonalMatrix[m], n][f],
		   m ] /;
	VectorQ[m] && SameQ[Length[m], Length[n]] && ! PatternQ

Downsample/: Downsample[m_List, n_List][f_] :=
	MyMessage[Downsample::badmD, f] /;
	VectorQ[m] && ! SameQ[Length[m], Length[n]]

Downsample/: Downsample[m_List, n_List][f_] :=
	MyMessage[Downsample::badmD, f] /;
	MatrixQ[m] && ! Apply[SameQ, Append[Dimensions[m], Length[n]]]

Downsample/: Downsample[m_List, n_Symbol][f_] :=
	MyMessage[Downsample::badmD, f]

		(*  degenerate cases  *)

Downsample/: Downsample[1, n_Symbol][c_] := c
Downsample/: Downsample[-1, n_Symbol][c_] := Rev[n][c]

	(* convert upsample into a functional form  *)

Downsample/: TheFunction[Downsample[m_, n_Symbol][f_]] :=
	TheFunction[f] /. ( n -> m n )
Downsample/: TheFunction[Downsample[m_, n_List][f_]] :=
	TheFunction[f] /. ( n ~ReplaceWith~ (m . n) ) /;
	MatrixQ[m] && VectorQ[n]

Downsample/: N [ Downsample[m_, n_][f_] ] :=
	N [ TheFunction[Downsample[m,n][f]] ]

	(*  extent of downsample operation  *)

Downsample/: Extent1D[ Downsample[m_, n_][f_], n_ ] :=
	intdivide[ Extent1D[f, n], m ]

	(* formatting of the downsample operator *)

Downsample/: Format[ Downsample[m_, n_List] ] :=
	Block [	{bar},
		bar = TableForm[ Apply[ Table,
					{ "  |  ", Dimensions[n] } ] ];
		Format[ SequenceForm[ Downsample,
				      Subscript[ TableForm[m] ],
				      Subscript[ bar ],
				      Subscript[ TableForm[n] ] ] ] ]
Downsample/: Format[ Downsample[m_, n_] ] :=
	Format[ Subscripted[Downsample[m,n]] ]

Downsample/: Format[ Downsample, TeXForm ] := "\\downarrow"
Downsample/: Format[ Downsample[m_, n_], TeXForm ] :=
	inTeXformat[ Downsample, m, n ]

Downsample/: GetOperatorVariables[ Downsample[m_, n_] ] := n

(*  DFT  [Myers, 221]  *)
DFT/: DFT[L_, n_] := DFT[L, n, DummyVariables[Length[n], Global`k]]
DFT/: DFT[L_, n_, k_] [ InvDFT[L_, ki_, n_] [f_] ] := 
	If [ SameQ[k, ki], f, f /. ReplaceWith[ki, k] ] /;
	Length[k] == Length[ki]

DFT/: DFT^-1[a__] := InvDFT[a]

DFT/: Format[ DFT, TeXForm ] := "{\\cal DFT}"
DFT/: Format[ DFT[L_, n_, k_], TeXForm ] :=
	StringForm[ "{\\cal DFT}_{``, L = ``}", n, L]

DFT/: GetOperatorVariables[ DFT[L_, n_, k_] ] := n

DFT/: Extent1D[ DFT[L_, n_, k_][f_], n_ ] := {0, L}
DFT/: Extent1D[ DFT[L_, n_, k_][f_], k_ ] := {0, L}

(*  DTFT  *)
DTFT/: DTFT[n_] := DTFT[n, DummyVariables[Length[n], Global`w]]
DTFT/: DTFT[n_, w_] [ InvDTFT[wi_, n_] [f_] ] :=
	If [ SameQ[w, wi], f, f /. ReplaceWith[wi, w] ] /;
	Length[w] == Length[wi]

DTFT/: DTFT^-1[a__] := InvDTFT[a]

DTFT/: Format[ DTFT, TeXForm ] := "{\\cal DTFT}"
DTFT/: Format[ DTFT[n_, w_], TeXForm ] := StringForm[ "{\\cal DTFT}_{``}", n ]

(*  Extent1D  *)
Extent1DUnion[ {x1_, y1_}, {x2_, y2_} ] := { Min[x1, x2], Max[y1, y2] }
Extent1DIntersection[ {x1_, y1_}, {x2_, y2_} ] := { Max[x1, x2], Min[y1, y2] }
Extent1DOrder[ x1_, y1_ ] := { Min[x1, y1], Max[x1, y1] }

Extent1D[t_][x_] := Extent1D[x, t]

Extent1D[x_ + y_, n_] := Extent1DUnion[ Extent1D[x, n], Extent1D[y, n] ]
Extent1D[x_ y_, n_] := Extent1DIntersection[ Extent1D[x, n], Extent1D[y, n] ]
Extent1D[0, n_] := {0, 0}
Extent1D[x_, n_] := {-Infinity, Infinity}

(*  FT  *)
FT/: FT[n_] := FT[n, DummyVariables[Length[n], Global`w]]
FT/: FT[n_, w_] [ InvFT[wi_, n_] [f_] ] :=
	If [ SameQ[w, wi], f, f /. ReplaceWith[wi, w] ] /;
	Length[w] == Length[wi]

FT/: FT^-1[a__] := InvFT[a]

FT/: Format[ FT, TeXForm ] := "{\\cal F}"
FT/: Format[ FT[n_, w_], TeXForm ] := StringForm[ "{\\cal F}_{``}", n ]

(*  GetDeltaFunctions  *)
choose[{a_, aterm_}] := If [ SameQ[aterm, 0], a, aterm ]

	(* rules specific for 1-D functions of t  *)
GetDeltaFunctions[ c_. Delta[a_. t_Symbol + b_.], t_Symbol ] :=
	c Delta[a t + b]

	(* rules specific for m-D functions of t1, t2, ...  *)
GetDeltaFunctions[ c_. Delta[a_. t_Symbol + b_.], tlist_List ] :=
	c Delta[a t + b] /;
	MemberQ[tlist, t]

	(* rules for all dimensions  *)
GetDeltaFunctions[ a_ + b_, rest__ ] :=
	GetDeltaFunctions[a, rest] + GetDeltaFunctions[b, rest]
GetDeltaFunctions[ h_[b__], rest__ ] :=
	Block [ {arglist, deltalist},
		arglist = {b};
		deltalist = Map[ GetDeltaFunctions[#1, rest]&, arglist ];
		If [ SameQ[Apply[Plus, deltalist], 0],
		     0,
		     Apply[h, Map[choose, Transpose[{arglist, deltalist}]]] ] ]
GetDeltaFunctions[ x_, rest__ ] := 0


(*  Impulse  *)
Impulse/: Impulse[Infinity] := 0		 	  (* definition       *)
Impulse/: Impulse[n_] := 0	/; N[ n != 0 ]
Impulse/: Impulse[n_] := 1	/; N[ n == 0 ]
Impulse/: Impulse[-Infinity] := 0

Impulse/: Extent1D[ Impulse[a_. n_ + b_.], n_ ] := { -b/a, -b/a } /;
	FreeQ[{a,b}, n]

Impulse/: Impulse[n_, args__] := Impulse[n] Impulse[args] (* multidimensional *)
Impulse/: Impulse[List[args__]] := Impulse[args]	  (* separability     *)

Impulse/: Format[ Impulse[n_], TeXForm ] :=		  (* output form rule *)
	  StringForm["\\delta[``]", n]

(*  Interleave [Myers, 219-220]  *)
Interleave/: Interleave[n_][x_] := x 			  (* definition	      *)
Interleave/: TheFunction[Interleave[n_][x1_, x__]] :=
	Block [	{index = Unique["k"], L, xlist},
		xlist = List[x1, x];
		L = Length[xlist];
		Summation[index, 0, L-1, 1]
			 [ Shift[index, n] [Upsample[L, n]
					     [ Element[xlist, index + 1] ] ] ] ]
Interleave/: N [ Interleave[n_][x1_, x__] ] :=
	N [ TheFunction[Interleave[n][x1,x]] ]

Interleave/: Format[ Interleave[n_] ] :=		(* output form rules *)
		Format[ Subscripted[Interleave[n]] ]
Interleave/: Format[ Interleave[n_], TeXForm ] :=
		Format[ Subscripted[Interleave[n]] ]


(*  IntervalsToFunction  *)
IntervalsToFunction[interlist_, n_:Global`n, step_:Step, pulse_:Pulse] :=
	Block [	{int, len, signal},

		MakeOne[interval_] :=
		  Block [
		    {f, nminus, nplus, pulselen},
		    f = interval[[1]];
		    nminus = interval[[2]];
		    nplus = interval[[3]];
		    pulselen = nplus - nminus + 1;
		    Which [ SameQ[nminus, -Infinity],
			      f step[nplus - n],
			    SameQ[nplus, Infinity],
			      f step[n - nminus],
			    True,
			      f pulse[pulselen, n - nminus] ] ];

		len = Length[interlist];
		signal = MakeOne[ interlist[[1]] ];
		For [ int = 2, int <= len, int++,
		      signal += MakeOne[ interlist[[int]] ] ];

		signal ]

(*  InvDFT  [Myers, 221]  *)
InvDFT/: InvDFT[L_, k_] :=
	InvDFT[L, k, DummyVariables[Length[k], Global`n]]
InvDFT/: InvDFT[L_, k_, ni_] [ DFT[L_, n_, k_] [f_] ] :=
	If [ SameQ[ni, n], f, f /. ReplaceWith[n, ni] ] /;
	Length[ni] == Length[n]

InvDFT/: Format[ InvDFT, TeXForm ] := "{\\cal DFT}^{-1}"
InvDFT/: Format[ InvDFT[L_, k_, n_], TeXForm ] :=
	  StringForm[ "{\\cal DFT}^{-1}_{``, L = ``}", k, L ]

InvDFT/: GetOperatorVariables[ InvDFT[L_, k_, n_] ] := k

InvDFT/: Extent1D[ InvDFT[L_, k_, n_][f_], n_ ] := {0, L}
InvDFT/: Extent1D[ InvDFT[L_, k_, n_][f_], k_ ] := {0, L}

(*  InvDTFT  *)
InvDTFT/: InvDTFT[w_] := InvDTFT[w, DummyVariables[Length[w], Global`n]]
InvDTFT/: InvDTFT[w_, ni_] [ DTFT[n_, w_] [f_] ] :=
	If [ SameQ[ni, n], f, f /. ReplaceWith[n, ni] ] /;
	Length[ni] == Length[n]

InvDTFT/: Format[ InvDTFT, TeXForm ] := "{\\cal DTFT}^{-1}"
InvDTFT/: Format[ InvDTFT[w_, n_], TeXForm ] :=
	  StringForm[ "{\\cal DTFT}^{-1}_{``}", w ]

(*  InvFT  *)
InvFT/: InvFT[w_] := InvFT[w, DummyVariables[Length[w], Global`t]]
InvFT/: InvFT[w_, ti_] [ FT[t_, w_] [f_] ] :=
	If [ SameQ[ti, t], f, f /. ReplaceWith[t, ti] ] /;
	Length[ti] == Length[t]

InvFT/: Format[ InvFT, TeXForm ] := "{\\cal F}^{-1}"
InvFT/: Format[ InvFT[w_, t_], TeXForm ] :=
	  StringForm[ "{\\cal F}^{-1}_{``}", w ]

(*  InvL  *)
InvL/: InvL[s_] := InvL[s, DummyVariables[Length[s], Global`t]]
InvL/: InvL[s_, ti_] [ L[t_, s_] [f_] ] :=
	If [ SameQ[ti, t], f, f /. ReplaceWith[t, ti] ] /;
	Length[t] == Length[ti]

InvL/: Format[ InvL, TeXForm ] := "{\\cal L}^{-1}"
InvL/: Format[ InvL[s_, t_], TeXForm ] := StringForm[ "{\\cal L}^{-1}_{``}", s ]

(*  InvZ  [Myers, 221]  *)
InvZ/: InvZ[zi_] := InvZ[zi, DummyVariables[Length[zi], Global`n]]
InvZ/: InvZ[z_, ni_] [ Z[n_, z_] [f_] ] :=
	If [ SameQ[ni, n], f, f /. ReplaceWith[n, ni] ] /;
	Length[n] == Length[ni]

InvZ/: Format[ InvZ[z_, n_], TeXForm ] := "{\\cal Z}^{-1}"
InvZ/: Format[ InvZ[z_, n_], TeXForm ] := StringForm[ "{\\cal Z}^{-1}_{``}", z ]

(*  L  *)
L/: L[t_] := L[t, DummyVariables[Length[t], Global`s]]
L/: L[t_, s_] [ InvL[si_, t_] [f_] ] :=
	If [ SameQ[s, si], f, f /. ReplaceWith[si, s] ] /;
	Length[s] == Length[si]

L/: L^-1[a__] := InvL[a]

L/: Format[ L, TeXForm ] := "{\\cal L}"
L/: Format[ L[t_, s_], TeXForm ] := StringForm[ "{\\cal L}_{``}", t ]

(*  LineImpulse  *)
Unprotect[samecoefficient];
Unprotect[mergecoefficient];

SetAttributes[samecoefficient, Listable];
samecoefficient[c1_, c2_] :=
	If [ SameQ[c1,c2],
	     ! SameQ[c1, Null],
	     SameQ[c1, Null] || SameQ[c2, Null] ]

SetAttributes[mergecoefficient, Listable];
mergecoefficient[c1_, c2_] := 
	Which [ SameQ[c1, Null], c2,
		SameQ[c2, Null], c1,
		True,            c1 ]

Protect[samecoefficient];
Protect[mergecoefficient];

LineImpulseMerge[nlist_, ncoeffs_, klist_, kcoeffs_] :=
	Block [	{commonvars, kscale, nscale, scaledkcoeffs, scaledncoeffs},
		scaledncoeffs = ncoeffs;
		scaledkcoeffs = kcoeffs;
		commonvars = Intersection[nlist, klist];
		If [ ! EmptyQ[commonvars],
		     normvar = First[commonvars];
		     nscale = AssociateItem[normvar, nlist, ncoeffs];
		     scaledncoeffs = ncoeffs / nscale;
		     kscale = AssociateItem[normvar, klist, kcoeffs];
		     scaledkcoeffs = kcoeffs / kscale ];
		allvars = Union[nlist, klist];
		coeffs1 = AssociateItem[allvars, nlist, scaledncoeffs];
		coeffs2 = AssociateItem[allvars, klist, scaledkcoeffs];
		{ Apply[And, samecoefficient[coeffs1, coeffs2]],
		  LineImpulse[allvars, mergecoefficient[coeffs1, coeffs2]] } ]

LineImpulse/: LineImpulse[n_?AtomQ, c_] := Impulse[c n]

LineImpulse/: LineImpulse[nlist_, ncoeffs_] LineImpulse[klist_, kcoeffs_] :=
	LineImpulseMerge[nlist, ncoeffs, klist, kcoeffs] [[2]] /;
	! EmptyQ[Intersection[nlist,klist]] && LineImpulseMerge[nlist, ncoeffs, klist, kcoeffs] [[1]]

LineImpulse/: TheFunction[ LineImpulse[nlist_, ncoeffs_] ] :=
	Impulse[ Apply[Plus, nlist, ncoeffs] ]

LineImpulse/: N [ LineImpulse[nlist_, ncoeffs_] ] :=
	N [ TheFunction[ LineImpulse[nlist, ncoeffs] ] ]

LineImpulse/: Format[ LineImpulse[nlist_List, ncoeffs_List] ] :=
	Format[Impulse[Apply[Plus, nlist ncoeffs]]] /;
	Length[nlist] == Length[ncoeffs]

(*  OperatorInOperExpr  *)
OperatorInOperExpr[ h_[a__][s__] ] := h[a]
OperatorInOperExpr[ h_[s__] ] := h

(*  OperatorName  *)
OperatorName[ h_[a__][s__] ] := h
OperatorName[ h_[a__] ] := h
OperatorName[ x_ ] := x

(*  ParametersInOperExpr  *)
ParametersInOperExpr[ h_[a__][s__] ] := a
ParametersInOperExpr[ h_[s__] ] := Null

(*  Periodic  *)
Periodic/: Format[ Periodic[k_, n_] ] := informat[Periodic, k, n]
Periodic/: Format[ Periodic[k_, n_], TeXForm ] := inTeXformat[Perodic, k, n]

Periodic/: GetOperatorVariables[ Periodic[k_, n_] ] := n

(*  PolyphaseDownsample  *)
PolyphaseDownsample/: GetOperatorVariables[ PolyphaseDownsample[m_, h_, n_] ] := n

(*  PolyphaseUpsample  *)
PolyphaseUpsample/: GetOperatorVariables[ PolyphaseUpsample[l_, h_, n_] ] := n

(*  Pulse  *)
Pulse[l_] := Pulse[l, Global`n]				(* default args    *)

Pulse[0, n_] := 0 					(* definition      *)
Pulse[1, n_] := Impulse[n]
Pulse[Infinity, n_] := Step[n]
Pulse/: Pulse[l_, n_] := 1 /; N[0 <= n < l]
Pulse/: Pulse[l_, n_] := 0 /; N[( n < 0 ) || ( n >= l )]

Pulse/: TheFunction[Pulse[l_, n_]] := Step[n] - Step[n - l]

Pulse/: Extent1D[ Pulse[l_, a_. n_ + b_.], n_ ] :=
	Which [ Positive[a], {-b/a, (l-b)/a},
		Negative[a], {(l-b)/a, -b/a},
		True,	     Extent1DOrder[-b/a, (l-b)/a] ] /;
	FreeQ[{a,b,l}, n]

Pulse/: Format[ Pulse[l_, n_] ] :=			(* output forms    *)
		Format[ Pulse Subscript[l] [n] ]

Pulse/: Format[ Pulse, TeXForm ] := "\\sqcap"
Pulse/: Format[ Pulse[l_, n_], TeXForm ] :=
		StringForm[ "\\sqcap_{``} [``]", l, n ]

Derivative[0,1][Pulse] :=				(* derivative rule *)
	( Impulse[#2] + Impulse[#2 - #1] )&

(*  RationalGCD  *)
RationalGCD[ ] := 0
RationalGCD[ x_ ] := x
RationalGCD[ arg1_, rest__ ] := rationalGCD[{arg1, rest}]

rationalGCD[ numlist_List ] :=
	Block [	{denomgcd, gcdfun, gcdlist, newlist, numergcd},
		newlist = Flatten[numlist];
		gcdlist = Map[ gcdfilter, newlist ];
		gcdfun = If [ TrueQ[$VersionNumber >= 2.0],
			      PolynomialGCD, GCD ];
		numergcd = Apply[gcdfun, Map[Numerator, gcdlist]];
		denomgcd = Apply[gcdfun, Map[Denominator, gcdlist]];
		numergcd / denomgcd ]

(*  Rev  *)
Rev/: Rev[n_][Rev[n_][x_]] := x
Rev/: Rev[n_][Rev[m_][x_]] :=
	Block [	{remaining},
		remaining = SetExclusion[ToList[n], ToList[m]];
		Which [ EmptyQ[remaining],
			  x,
			SameQ[Length[remaining], 1],
			  Rev[ remaining[[1]] ] [x],
			True,
			  Rev[remaining][x] ] ] /;
	! EmptyQ[Intersection[ToList[n], ToList[m]]]

Rev/: Extent1D[ Rev[n_][f_], n_ ] :=
	Block [	{extent},
		extent = Extent1D[f, n];
		{ -Second[extent], -First[extent] } ]

Rev/: TheFunction[Rev[n_][x_]] := ( TheFunction[x] /. n ~ReplaceWith~ (-n) )
Rev/: N[ Rev[n_][x_] ] := N [ TheFunction[Rev[n][x]] ]

Rev/: Format[ Rev[n_] ] := Format[ Subscripted[Rev[n]] ]
Rev/: Format[ Rev[n_], TeXForm ] := Format[ Subscripted[Rev[n]] ]

(*  RewriteSummations  *)
donotrewrite[f_, t_] := FreeQ[f, Summation] && FreeQ[f, Periodic[k_,t] ]

RewriteSummations[ f_, t_, tl_, tu_ ] := f /; donotrewrite[f, t]

RewriteSummations[ f_, t_, tl_, tu_ ] := 
	Block [	{ii},
		f /. { Periodic[k_, v_][g_] :>		     (* use v_ not t *)
			 rewriteonesummation[ g /. t -> t + ii k,
					      t, tl, tu,
					      ii, -Infinity, Infinity, 1 ] /;
			 SameQ[v, t] && donotrewrite[g, t],
		       Summation[i_, il_, iu_, inc_][g_] :>
			 rewriteonesummation[ g, t, tl, tu, i, il, iu, inc ] /;
			 donotrewrite[g, t] } ]

rewriteonesummation[ g_, t_, tl_, tu_, i_, il_, iu_, inc_ ] :=
	Block [ {ext, il0, ilact, iu0, iuact, lext, lsol, uext, usol, valid},
		ext = Simplify[ Extent1D[g, t] ] //.
			SignalProcessing`Support`SupCode`Private`MinMaxRules;
		valid = ListQ[ext] && SameQ[Length[ext], 2];
		If [ ! valid, ext = { -Infinity, Infinity } ];

		{lext, uext} = ext;
		If [ FreeQ[lext, i] && ! FreeQ[uext, i], lext = uext ];
		If [ FreeQ[uext, i] && ! FreeQ[lext, i], uext = lext ];

		If [ FreeQ[lext, i] && FreeQ[uext, i],
		     Message[ Summation::extent ];
		       Summation[i, il, iu, inc][g],

		     lsol = Solve[ lext == tl, i ];
		     usol = Solve[ uext == tu, i ];
		     If [ SameQ[Head[lsol], Solve] || SameQ[Head[usol], Solve],
			  Message[ Summation::extent ];
		            Summation[i, il, iu, inc][g],

		          il0 = Min[ Map[ GetRoot, lsol ] ];
		          iu0 = Max[ Map[ GetRoot, usol ] ];
		          ilact = Max[il, Min[il0, iu0]];
		          iuact = Min[iu, Max[il0, iu0]];
			  If [ InfinityQ[-Infinity],
		               ilact = ilact - Mod[ilact, inc];
		                 iuact = iuact - Mod[iuact, inc] + inc,
		               ilact = ilact - Mod[ilact - il, inc];
		                 iuact = iuact - Mod[iuact - il, inc] + inc ];
		          Sum[ g, {i, ilact, iuact, inc} ] ] ] ]


(*  ScaleAxis [Myers, 220]  *)
ScaleAxis/: TheFunction[ScaleAxis[l_, w_Symbol][x_]] :=
		TheFunction[x] /. w -> l w
ScaleAxis/: TheFunction[ScaleAxis[{l__},{w__}][x_]] :=
		( TheFunction[x] /. {w} ~ReplaceWith~ ({l} . {w}) )

ScaleAxis/: N[ ScaleAxis[l_, w_][x_] ] := N[ TheFunction[ ScaleAxis[l,w][x] ] ]

ScaleAxis/: Format[ ScaleAxis[l_,w_] ] := informat[ScaleAxis, l, w]
ScaleAxis/: Format[ ScaleAxis[l_,w_], TeXForm ] := inTeXformat[ScaleAxis, l, w]

ScaleAxis/: GetOperatorVariables[ ScaleAxis[l_, w_] ] := w

(*  ScalingFactor --  must be independent of variable z  *)
ScalingFactor[f_, z_] :=
	breakup[ Apply[ RationalGCD, GetAllFactors[f, z] ], z ]

(*  SequencePlot  *)
Options[ SequencePlot ] :=
	{ Axes -> axesdefaultvalue,
	  DisplayFunction :> $DisplayFunction,
	  PlotLabel -> "Sequence Plot",
	  PlotRange -> All }

baseticks = { 1, 2, 2, 5, 5, 5, 5, 5, 10, 10 };
getticks[nstart_, nend_] :=
	Block [	{factor, i, ninc, ns, nterm},
		ninc = Round[(nend - nstart) / 5];
		If [ ninc == 0, ninc = 1 ];
		If [ ninc > 2,
		     factor = 10;
		     While [ ninc > factor, factor *= 10 ];
		     factor /= 10;
		     nterm = intdivide[ninc, factor];
		     ninc = factor * baseticks[[nterm]] ];
		ns = Sign[nstart] ninc intdivide[Abs[nstart], ninc];
		Table[i, {i, ns, nend, ninc}] ]

dgetvalue[f_, n_Symbol, n0_Integer] := GetValue[f, n, n0]
dgetvalue[f_, n_Symbol, n0_Real] := GetValue[f, n, Round[n0]]
dgetvalue[f_, {n1_Symbol, n2_Symbol}, {n01_?NumberQ, n02_?NumberQ}] :=
	GetValue[f, {n1, n2}, {Round[n01], Round[n02]}]

SequencePlot[f_List, {n1_Symbol, st1_, end1_},
		     {n2_Symbol, st2_, end2_}, op___] :=
	MyMessage[ SequencePlot::nolist, NullPlot ]

SequencePlot[f_, {n1_Symbol, st1_, end1_}, {n2_Symbol, st2_, end2_}, op___] :=
	Block [	{funct, ftemp, funs, n1length, n2length, n1str, n2str,
		 n1temp, n2temp, nstart, nend, oplist, plot},

		nstart = Floor[N[{st1, st2}]];
		nend   = Ceiling[N[{end1, end2}]];
		n1length = nend[[1]] - nstart[[1]];
		n2length = nend[[2]] - nstart[[2]];
		n1str = StripPackage[n1];
		n2str = StripPackage[n2];
		oplist = Join [	ToList[op],
				Options[SequencePlot],
				{ PlotPoints -> {n1length + 1, n2length + 1},
				  AxesLabel -> {n1str, n2str, " "},
				  Ticks -> Automatic } ];

		funct = ToDiscrete[f];
		ftemp = N[ TheFunction[ funct ] ];

		(*  Generate three-dimensional coordinates for plot  *)
		Off[Power::infy, Infinity::indt];
		plot = Apply[ Plot3D,
			      Join[ { dgetvalue[ ftemp, {n1, n2},
						 {n1temp, n2temp} ],
				      { n1temp, nstart[[1]], nend[[1]] },
				      { n2temp, nstart[[2]], nend[[2]] } },
				    oplist ] ];
		On[Power::infy, Infinity::indt];

		plot ]


SequencePlot[f_, {n_Symbol, start_, end_}, op___] :=
	Block [	{factor, ftemp, funct, funs, length, linethickness,
		 nend, nstart, nstr, ntemp, oplist, plot,
		 pointthickness, showoptions},

		nstart = Floor[N[start]];
		nend   = Ceiling[N[end]];

		(*  Set up the options  *)
		nticks = getticks[nstart, nend];
		nstr = StripPackage[n];
		oplist = Join [	Flatten[ToList[op]],
				Options[SequencePlot],
			        { AxesLabel -> { nstr, " " },
				  Ticks -> { nticks, Automatic } } ];

		factor = N[ Log[10, nend - nstart + 3 ] ];
		linethickness = 0.008 / factor;
		pointthickness = 0.024 / factor;

		(*  Extract the functional form of the expression  *)
		funct = ToDiscrete[f];
		ftemp = N[ TheFunction[ funct ] ];

		(*  Handle any infinite summations  *)
		ftemp = RewriteSummations[ftemp, n, nstart, nend];

		(*  Generate form for 1-D function  *)
		SetAttributes[seqFun1D, Listable];    (* avoid fun pointers *)
		seqFun1D[g_] :=
			DiscreteGraphics [
				Table [ {ntemp, dgetvalue[g, n, ntemp]},
					{ntemp, nstart, nend}],
				linethickness,
				pointthickness ];

		(*  Plot "time"-domain	*)
		Off[Power::infy, Infinity::indt];
		showoptions = RemoveOptions[oplist, specialPlotOptions];
		plot = Show [ seqFun1D[ftemp], showoptions ];
		On[Power::infy, Infinity::indt];

		plot ]

(*  SequenceToFunction  *)
SequenceToFunction[seq_, n_:Global`n, noffset_:0] :=
	Block [ {},
		(* Variable no. of arguments -- could not use Function *)
		impulsegen[l__] := seq[[l]] Apply[Impulse,
						  n - {l} - noffset + 1];
		Apply[ Plus, Flatten [ Array[impulsegen, Dimensions[seq]] ] ] ]

(*  SignalPlot  *)
specialPlotOptions =
	Complement[ Sort[ Map[First, Options[Plot]] ],
		    Sort[ Map[First, Options[Graphics]] ] ]

sigPlot1D[args__][g_] := SignalPlot[g, args]

Options[ SignalPlot ] :=
	{ Axes -> axesdefaultvalue,
	  DisplayFunction :> $DisplayFunction,
	  PlotLabel -> "Signal Plot",
	  PlotRange -> All,
	  Ticks -> Automatic }

SignalPlot[ f_List, {t1_Symbol, start1_, end1_},
		    {t2_Symbol, start2_, end2_}, op___ ] :=
	MyMessage[ SignalPlot::nolist, NullPlot ]

SignalPlot[ f_,	{t1_Symbol, start1_, end1_},
		{t2_Symbol, start2_, end2_}, op___ ] :=
	Block [	{deltafuns, freqresp, ftemp, funct, oplist,
		 t1str, t2str, t1temp, t2temp},

		(*  Set up for plotting  *)
		t1str = StripPackage[t1];
		t2str = StripPackage[t2];
		oplist = Join[ ToList[op],
			       { AxesLabel -> { t1str, t2str, " " } },
			       Options[ SignalPlot ] ];

		(*  Make function continuous and replace non-sep delta's *)
		funct = ToContinuous[f] /.
		     Delta[a_. t1 + b_. t2 + c_.] :> Impulse[a t1 + b t2 + c];

		(*  Separate delta functions from rest of funct  *)
		deltafuns = GetDeltaFunctions[funct, {t1, t2}];
		ftemp = funct /.
			{ h_[ c_. Delta[a_. t1 + b_.] ] :> 0,
			  h_[ c_. Delta[a_. t1 + b_.] + x_ ] :> h[x],
			  Delta[a_. t1 + b_.] :> 0,
			  h_[ c_. Delta[a_. t2 + b_.] ] :> 0,
			  h_[ c_. Delta[a_. t2 + b_.] + x_ ] :> h[x],
			  Delta[a_. t2 + b_.] :> 0 };

		(*  Convert freq. response into plottable function form  *)

		ftemp = N [ TheFunction[funct] ];

		(*  Plot delta functions separately from rest of function  *)

		Plot3D[ GetValue[ftemp, {t1, t2}, {t1temp, t2temp}],
			{ t1temp, start1, end1 },
			{ t2temp, start2, end2 },
			Release[oplist] ] ]

SignalPlot[f_List, {t_Symbol, start_, end_}, op___] :=
	Block [	{fun, oplist, showoptions, theplots},
		fun = sigPlot1D[{t,start,end}, DisplayFunction->Identity, op];
		theplots = Map[ fun, f ];
		oplist = Join[ToList[op], Options[SignalPlot]];
		showoptions = RemoveOptions[oplist, specialPlotOptions];
		Show[ theplots, ToCollection[showoptions] ] ]

SignalPlot[f_, {t_Symbol, start_, end_}, op___] :=
	Block [	{deltafuns, ftemp, funct, i, imstyle,
		 nend = N[end], nstart = N[start], oplist,
		 plot, showoptions, restyle, tstr, ttemp },

		(*  Set up for plotting  *)

		tstr = StripPackage[t];
		restyle = Thickness[0.006];
		imstyle = Dashing[{0.05,0.05}]; 
		oplist = Join[ ToList[op],
			       { AxesLabel -> { tstr, " " },
				 DisplayFunction -> Identity,
				 PlotStyle -> { restyle, imstyle } },
			       Options[ SignalPlot ] ];

		(*  Separate delta functions from rest of funct  *)

		funct = ToContinuous[f];
		deltafuns = GetDeltaFunctions[funct, t];
		ftemp = funct /.
			{ h_[ c_. Delta[a_. t + b_.] ] :> 0,
			  h_[ c_. Delta[a_. t + b_.] + x_ ] :> h[x],
			  Delta[a_. t + b_.] :> 0 };

		(*  Convert function into plottable functional form  *)

		ftemp = N[ TheFunction[ftemp] ];

		(*  Handle any infinite summations  *)

		ftemp = RewriteSummations[ftemp, t, nstart, nend];
	
		(*  Plot delta functions separately from rest of ftemp  *)

		Off[Power::infy, Infinity::indt];

		(*  1.  Generate the graphics for the entire signal  *)

		If [ ZeroQ[ftemp],
		     plot = If [ ! AtomQ[deltafuns],   (* delta funs found *)
				 DeltaPlot[deltafuns, t, nstart, nend],
				 NullPlot ];
		     plot = { Plot[ 0,
				    { ttemp, nstart, nend},
				    Release[oplist] ],
			      plot },

		     plot = Plot[ { Re[GetValue[ftemp, t, ttemp]],
				    Im[GetValue[ftemp, t, ttemp]] },
				  { ttemp, nstart, nend },
				  Release[oplist] ];
		     If [ ! AtomQ[deltafuns],      (* delta funs found *)
		          plot = DeltaPlot[deltafuns, t, nstart, nend, plot]] ];

		(*  2.  Plot the entire signal  *)

		oplist = Join[ ToList[op], Options[SignalPlot] ];
		showoptions = RemoveOptions[oplist, specialPlotOptions];

		plot = Show[ plot, ToCollection[showoptions] ];

		On[Power::infy, Infinity::indt];

		plot ]


(*  SignalsInOperExpr  *)
SignalsInOperExpr[ h_[s__] ] := s
SignalsInOperExpr[ h_[a__][s__] ] := s

SignalsInOperExpr[ h_[s__], h_ ] := s
SignalsInOperExpr[ h_[s__], op_ ] := h[s]		/; ! SameQ[h, op]

SignalsInOperExpr[ h_[a__][s__], h_ ] := s
SignalsInOperExpr[ h_[a__][s__], op_ ] := h[a][s]	/; ! SameQ[h, op]

(*  Shift [Myers, 217]  *)
Shift/: Shift[m_, n_List] := Shift[ Table[m, Length[n]], n ] /; ! ListQ[m]

Shift/: Shift[m_List, n_List][x_] := MyMessage[ Shift::badmD, x ] /;
	Length[m] != Length[n]

Shift/: Shift[0, n_][x_] := x

Shift/: Extent1D[ Shift[m_, n_][f_], n_ ] := Extent1D[f, n] + m

Shift/:	TheFunction[Shift[k_,n_][x_]] :=
	( TheFunction[x] /. n ~ReplaceWith~ (n - k) )
Shift/: N[ Shift[k_,n_][x_] ] := N [ TheFunction[ Shift[k,n][x] ] ]

Shift/: Format[ Shift[m_,n_] ] := informat[Shift, m, n]
Shift/: Format[ Shift[m_,n_], TeXForm ] := inTeXformat[Shift, m, n]

Shift/: GetOperatorVariables[ Shift[l_, n_] ] := n


(*  Sinc  *)
Sinc/: Sinc[-Infinity] := 0
Sinc/: Sinc[0] := 1
Sinc/: Sinc[0.] := 1.
Sinc/: Sinc[Infinity] := 0

Sinc/: Sinc[x_?NumberQ] := Sin[x] / x

Sinc/: Sinc[x_, y__] := Sinc[x] Sinc[y]

Sinc/: Format[Sinc, TeXForm] := "{\\rm sinc}"

Sinc/: Derivative[1][Sinc] := (( Cos[#1] / #1 ) - ( Sinc[#1] / #1 ))&


(*  SPSimplify  *)
Options[SPSimplify] :=
	Join[{ Dialogue -> False, Variables -> {} }, Options[Simplify]]

SPSimplify[ x_, op___ ] :=
	Block [	{dialflag, dialogue, different, laste, newe,
		 oplist, trig, varrules, varflag, variables},
		oplist = ToList[op] ~Join~ Options[SPSimplify];
		trig = Replace[ Trig, oplist ];
		dialogue = Replace[ Dialogue, oplist ];
		dialflag = InformUserQ[dialogue];
		variables = Replace[ Variables, oplist ];
		varflag = ! EmptyQ[variables];
		If [ varflag,
		     varrules = makesimplifyrules[variables] ];

		newe = x;
		different = True;
		While [ different,
			laste = newe;
			newe = If [ TrueQ[ $VersionNumber >= 2.0 ],
				    Simplify[laste, Trig -> trig],
				    Simplify[laste] ];
			newe = newe /. SPSimplificationRules;
			If [ varflag, newe = newe /. varrules ];
			different = ! SameQ[laste, newe];
			If [ dialflag && different,
			     Print["which simplifies to"]; Print[newe] ] ];
		newe ]

makesimplifyrules[var_Symbol] := makesimplifyrules[var, {}]
makesimplifyrules[var1_Symbol, var2_Symbol, rest___] :=
	makesimplifyrules[{var1, var2, rest}]
makesimplifyrules[varlist_List] :=
	Apply[ Join,
	       Map[ makesimplifyrules[#,
				      Complement[varlist, ToList[#]]]&,
				      varlist] ]
makesimplifyrules[t_, varlist_] :=
	{ Delta[a_ t + b_.] :> Delta[t + b/a] / Abs[a] /;
	    FreeQ[a, t] && MyFreeQ[b, varlist],
	  d_. (y_. + amp1_. Delta[t + t0_.]) :>
	    d y + Limit[d amp1, t -> -t0] Delta[t + t0] /;
	    FreeQ[t0, t] && (! FreeQ[d amp1, t]) &&
	    MyFreeQ[t0, varlist] && ! InfinityQ[ Limit[d amp1, t -> -t0] ],
	  c_. Sign[b_. t] + c_ :> 2 c CStep[b t] /;
	    FreeQ[b, t],
	  CStep[a_ t + b_.] :> CStep[t + b/a] / Abs[a] /;
	    Positive[a] && FreeQ[a, t] && MyFreeQ[b, varlist],
	  CStep[a_ t + b_.] :> CStep[-t - b/a] / Abs[a] /;
	    Negative[a] && FreeQ[a, t] && MyFreeQ[b, varlist] }

(*  Step  *)
Step/: Step[Infinity] := 1				(* definition	    *)
Step/: Step[t_] := 1		/; N[ t > 0 ]
Step/: Step[t_] := 0		/; N[ t < 0 ]
Step/: Step[0] := 1
Step/: Step[0.] := 1
Step/: Step[-Infinity] := 0

Step/: Extent1D[ Step[a_. n_ + b_.], n_ ] :=
	Which [ Positive[a], {-b/a, Infinity},
		Negative[a], {-Infinity, -b/a},
		True,	     Extent1DOrder[-b/a, Sign[a] Infinity] ] /;
	FreeQ[{a,b,l}, n]

Step/: Step[t_, args__] := Step[t] Step[args]		(* multidimensional *)
Step/: Step[List[args__]] := Step[args]			(* separability	    *)

Step/: Step[t_]^a_ := Step[t]				(* automatic	    *)
Step/: Step[n_Symbol + m_.] Step[n_Symbol  + k_.] :=	(* simplification   *)
	Step[n + Min[m,k]] /;				(* rules	    *)
	NumberQ[m] && NumberQ[k]
Step/: Step[n_Symbol + m_.] Step[-n_Symbol + k_.] :=
	Pulse[k + m, n + m] /;
	N[ -m <= k ]
Step/: Step[n_Symbol + m_.] Step[-n_Symbol + k_.] := 0 /;
	N[ -m > k ]

Step/: Simplify[Step[a_ n_Symbol + b_.]] :=		(* man. simp. rules *)
	Step[n + b/a] /;
	MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ]
Step/: Simplify[Step[- a_ n_Symbol + b_.]] :=
	Step[-n + b/a] /;
	MyFreeQ[{a,b}, n] && Implies[ NumberQ[N[a]], N[a > 0] ]

Step/: Format[ Step[x_], TeXForm ] :=			(* output form      *)
	StringForm["u[``]", x]

Derivative[1][Step] := Impulse				(* derivative rule  *)

(*  Summation  *)
Summation/: Summation[i_, n_] := Summation[i, 1, n, 1]
Summation/: Summation[i_, istart_, iend_] := Summation[i, istart, iend, 1]

Summation/: TheFunction[Summation[i_Symbol, istart_, iend_, inc_][f_]] :=
	Block [	{tempsym},
		If [ NumberQ[istart] && NumberQ[iend] && NumberQ[inc],
		     Sum[ TheFunction[f /. i -> tempsym],
			  {tempsym, istart, iend, inc} ],
		     Summation[i, istart, iend, inc][ TheFunction[f] ] ] ]

Summation/: Extent1D[ Summation[i_, istart_, iend_, inc_][f_], n_ ] :=
	Block [	{extent},
		extent = Extent1D[f, n];
		Extent1DOrder[ extent /. i -> istart,
			       extent /. i -> iend ] ] /;
	FreeQ[{i,istart,iend,inc}, n]

Summation/: N[Summation[i_Symbol, istart_?NumberQ, iend_?NumberQ, inc_?NumberQ][f_]] :=
	Sum[ N[f], {i, istart, iend, inc} ]

Summation/: Format[ Summation[i_, istart_, iend_, 1], TeXForm ] :=
	StringForm[ "\\sum_{`` = ``}^{``}", i, istart, iend ]

Summation/: Format[ Summation[i_, istart_, iend_, inc_], TeXForm ] :=
	StringForm["\\sum_{``=``, `` = `` + ``}^{``}",
		   i, istart, i, i, inc, iend]

	(* The following rules were suggested by John Mitchell *)
getdvars[ x_Symbol ] := x
getdvars[ {x_Symbol, k_} ] := x
Unprotect[D]
D/: D[expr_. Summation[i_, istart_, iend_, inc_][f_], x__] :=
	Summation[i, istart, iend, inc][ D[expr f, x] ] /;
	MyFreeQ[ i, Map[getdvars, {x}] ]
Protect[D]
Unprotect[Integrate]
Integrate/: Integrate[ expr_. Summation[i_, istart_, iend_, inc_][f_], x__] :=
	 Summation[i, istart, iend, inc][ Integrate[expr f, x] ] /;
	 MyFreeQ[ i, Map[getdvars, {x}] ]
Protect[Integrate]

(*  TheFunction  *)
TheFunction[x_?AtomQ] := x
TheFunction[Hold[x__]] := Hold[x]
TheFunction[Sign[x_]] := CStep[x] - CStep[-x]	(* Step or CStep will work  *)
TheFunction[h_[a_]] := h [ TheFunction[a] ]
TheFunction[h_[a1_, a__]] := Apply[ h, Map[TheFunction, { a1, a }] ]

(*  ToContinuous  *)
ToContinuous[expr_] := ReplaceRepeated[expr, ToContinuousRules]

ToContinuousRules =
     {
	Impulse[n_]   :> Delta[n],
	Step[n_]      :> CStep[n],
	Pulse[l_, n_] :> CPulse[l, n],
	Convolve[n_]  :> CConvolve[n]
     }

(*  ToDiscrete  *)
ToDiscrete[expr_] := ReplaceRepeated[expr, ToDiscreteRules]

ToDiscreteRules =
     {
	CStep[t_]      :> Step[t],
	Delta[t_]      :> Impulse[t],
	CPulse[l_, t_] :> Pulse[l, t],
	CConvolve[t_]  :> Convolve[t]
     }

(*  Unit *)
Unit/: Unit[0][t_] := Delta[t]				(* definition       *)
Unit/: Unit[-1][t_] := CStep[t]
Unit/: Unit[n_][t_?NumberQ] :=
	 t^-(n + 1) / Factorial[-(n + 1)] CStep[t] /; N[n < 0]

Unit/: Extent1D[ Unit[n_Integer][a_. t_ + b_.], t_ ] := {-b/a, -b/a} /;
	FreeQ[{a,b}, t] && ( n <= 0 )
Unit/: Extent1D[ Unit[n_Integer][a_. t_ + b_.], t_ ] :=
	Extent1DOrder[ -b/a, Sign[a] Infinity ] /;
	FreeQ[{a,b}, t] && ( n > 0 )

Unit/: Unit[n_][List[args]] := Unit[n][args]		(* multidimensional *)
Unit/: Unit[n_][t1_, t__] := Unit[n][t1] Unit[n][t]	(* separability     *)

Unit/: TheFunction[Unit[n_][t__]] :=
	 t^-(n + 1) / Factorial[-(n + 1)] CStep[t] /; N[n < 0]

Unit/: Format[ Unit[n_] ] := Subscripted[ Unit[n] ]	(* output forms     *)
Unit/: Format[ Unit[n_], TeXForm ] := StringForm[ "u_{``}", n ]

Derivative[1][Unit[n_]] := Unit[n+1]			(* derivative rule  *)

Unprotect[Integrate]					(* integration rules *)
Integrate/: Integrate[ expr_. Unit[n_][a_. x_ + b_.], x_Symbol ] :=
	Limit[expr, x -> -b/a] Unit[n-1][a x + b] / a /;
	FreeQ[{a,b,c}, x]
Integrate/: Integrate[ expr_. Unit[n_][a_ x_ + b_.], {x_Symbol, x1_, x2_}] :=
	Integrate[expr Unit[n][x + b/a] / a, x] /;
	FreeQ[{a,b}, x]
Integrate/: Integrate[expr_. Unit[n_][x_ + c_.], {x_Symbol, x1_, x2_}] :=
	Limit[expr, x -> -c] CStep[-c - x1] CStep[c + x2] Unit[n-1][x + c] /;
	FreeQ[c, x]	(* CStep[-c - x1] CStep[c + x2] <--> x1 <= -c <= x2 *)
Protect[Integrate]

(*  Upsample [Myers, 218-219] [Bamberger, ch. 2]  *)

	(* Error checking of upsample parameters  *)

Upsample/: Upsample[m_List, n_List][f_] :=
	MyMessage[Upsample::convert, Upsample[DiagonalMatrix[m], n][f], m] /;
	VectorQ[m] && SameQ[Length[m], Length[n]]

Upsample/: Upsample[m_List, n_List][f_] :=
	MyMessage[Upsample::badmD, f] /;
	VectorQ[m] && ! SameQ[Length[m], Length[n]]

Upsample/: Upsample[m_List, n_List][f_] :=
	MyMessage[Upsample::badmD, f] /;
	MatrixQ[m] && ! Apply[SameQ, Append[Dimensions[m], Length[n]]]

Upsample/: Upsample[m_List, n_Symbol][f_] :=
	MyMessage[Upsample::badmD, f]

		(*  degenerate cases  *)

Upsample/: Upsample[1, n_][f_] := f
Upsample/: Upsample[-1, n_][f_] := Rev[n][f]

	(* convert upsample into a functional form  *)

Upsample/: N[Upsample[m_,n_][f_List]] := UpsampleSequence[f, m, 0]

Upsample/: TheFunction[Upsample[m_, n_Symbol][f_List]] :=
		UpsampleSequence[f, m, 0]
Upsample/: TheFunction[Upsample[m_, n_Symbol][f_]] :=
		UpsampledFunction[m, 0, n, TheFunction[f] /. n -> n / m]

	(* extent of upsampled function  *)

Upsample/: Extent1D[ Upsample[m_, n_][f_], n_ ] := m Extent1D[f, n]

	(* formatting of the upsample operator *)

Upsample/: Format[ Upsample[m_, n_List] ] :=
	Block [	{bar},
		bar = TableForm[ Table[ "  |  ", Release[Dimensions[n]] ] ];
		Format[ SequenceForm[ Upsample,
				      Subscript[ TableForm[m] ],
				      Subscript[ bar ],
				      Subscript[ TableForm[n] ] ] ] ]
Upsample/: Format[ Upsample[m_, n_] ] :=
		Format[ Subscripted[Upsample[m, n]] ]

Upsample/: Format[ Upsample, TeXForm ] := "\\uparrow"
Upsample/: Format[ Upsample[m_, n_], TeXForm ] :=
		inTeXformat[ Upsample, m, n ]

Upsample/: GetOperatorVariables[ Upsample[l_, n_] ] := n

(*  UpsampledFunction  *)
UpsampledFunction[m_Integer, l_, n_Integer, f_] :=
	If [ IntegerQ[n / m], f, l ]
UpsampledFunction[m_?NumberQ, l_, n_?NumberQ, f_] :=
	If [ IntegerQ[Rationalize[n / m]], f, l ]

(*  UpsampleFactor --  must be independent of variable z  *)
UpsampleFactor[f_, z_] :=
	Block [	{ratgcd},
		If [ FreeQ[f, z], Return[0] ];
		ratgcd = Apply[ RationalGCD, GetAllExponents[f, z] ];
		ratgcd = breakup[ratgcd, z];
		If [ IntegerQ[ratgcd],
		     ratgcd,
		     Numerator[ratgcd] ] ]

(*  UpsampleSequence  *)
UpsampleSequence[seq_, m_:2, fill_:0] :=
	Block [	{fillcount, input, len, newseq},

		len = Length[seq];
		newseq = List[ seq[[1]] ];

		(*  Build the new sequence newseq by appending
		    m - 1 fill elements and then by appending
		    a value from the original sequence.  *)

		For [ input = 2, input <= len, input++,
		      For [ fillcount = 1, fillcount < m, fillcount++,
			    AppendTo[newseq, fill] ];
		      AppendTo[newseq, seq[[input]] ] ];

		newseq ]  /;
	IntegerQ[m] && ( m >= 2 )

(*  Z  *)
Z/: Z[n_] := Z[n, DummyVariables[Length[n], Global`z]]
Z/: Z[n_, z_] [ InvZ[zi_, n_] [f_] ] := 
	If [ SameQ[z, zi], f, f /. ReplaceWith[zi, z] ] /;
	Length[z] == Length[zi]

Z/: Z^-1[a__] := InvZ[a]

Z/: Format[ Z, TeXForm ] := "{\\cal Z}"
Z/: Format[ Z[n_, z_], TeXForm ] := StringForm[ "{\\cal Z}_{``}", n ]

(*  E N D     P A C K A G E  *)

End[]
EndPackage[]


(*  A L I A S E S  *)

Unprotect[ { AliasedSinc, AliasSinc, ASinc } ]
Clear[ AliasedSinc, AliasSinc, ASinc ]
AliasedSinc = Dirichlet
AliasedSinc::usage = Dirichlet::usage
AliasSinc = Dirichlet
AliasSinc::usage = Dirichlet::usage
ASinc = Dirichlet
ASinc::usage = Dirichlet::usage
Protect[ { AliasedSinc, AliasSinc, ASinc } ]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ]


(*  H E L P     I N F O R M A T I O N  *)

Block [	{newfuns},
	newfuns = 
	  { DeltaIntegrate,	DeltaPlot,		DiscreteGraphics,
	    DummyVariables,	GetDeltaFunctions,	IntervalsToFunction,
	    OperatorInOperExpr, OperatorName,		ParametersInOperExpr,
	    ScalingFactor,	SequencePlot,		SequenceToFunction,
	    SignalPlot,		SignalsInOperExpr,	"Signum",
	    SPSimplify,		TheFunction,		ToContinuous,
	    ToDiscrete,		UpsampledFunction,	UpsampleFactor,
	    UpsampleSequence };
	Combine[ SPfunctions, newfuns ];
	Apply[ Protect, newfuns ];
	(* Under Mma 2.0, Delta is used by Integrate *)
	If [ TrueQ[ $VersionNumber >= 2.0 ],
	     Unprotect[Delta] ] ]

Block [	{newoperators},
	newoperators =
	  { Aliasby,		CConvolve,		CircularShift,
	    Convolve,		DFT,		  	Difference,
	    Downsample,		DTFT,			FT,
	    Interleave,		InvDFT,			InvDTFT,
	    InvFT,		InvL,			InvZ,
	    L,			PolyphaseDownsample,	PolyphaseUpsample,
	    Periodic,		Rev,		  	RewriteSummations,
	    ScaleAxis,		Shift,		  	Summation,
	    Upsample,		Z };
	Combine[ SPoperators, newoperators ];
	Apply[ Protect, newoperators ] ]

Block [	{newsignals},
	newsignals =
	  { "AliasSinc",	"ASinc",		CPulse,
	    CStep,		Delta,			Dirichlet,
	    Impulse,		LineImpulse,		Pulse,
	    Sinc,		Step,			Unit };
	Combine[ SPsignals, newsignals ];
	Apply[ Protect, newsignals ] ]


(*    When encoded using ":=", these next rules causes Simplify to   *)
(*  carry out incorrect simplifications when an expression contains  *)
(*  an Impulse term and a term containing 1/n because Simplify will  *)
(*  group the expression over a common denominator thereby putting   *)
(*  the n in front of Impulse[n] which whould force it to zero,      *)
(*  thereby eliminating a perfectly valid term.			     *)

Unprotect[ SPSimplificationRules ]
AppendTo[ SPSimplificationRules,
	  ( f_ Impulse[n_Symbol + m_.] :> (f /. n -> -m) Impulse[n + m] /;
	    NumberQ[m] && ! FreeQ[f, n] ) ]
AppendTo[ SPSimplificationRules, ( x_ Sinc[b_. x_] :> Sin[b x] / b ) ]
AppendTo[ SPSimplificationRules, ( Sinc[-x_] :> Sinc[x] ) ]
AppendTo[ SPSimplificationRules, ( Delta[-x_] :> Delta[x] ) ]
AppendTo[ SPSimplificationRules, ( Delta[a_?Positive x_] :> Delta[x] / a ) ]
AppendTo[ SPSimplificationRules, ( Impulse[-x_] :> Impulse[x] ) ]
AppendTo[ SPSimplificationRules, ( CPulse[L_, L_ - x_] :> CPulse[L, x] ) ]

(*  Simplifications when an impulse is the input      *)

AppendTo[ SPSimplificationRules,
	  ( Downsample[L_, n_][a_. Impulse[n_]] :> a Impulse[n] ) ]
AppendTo[ SPSimplificationRules,
	  ( Upsample[L_, n_][a_. Impulse[n_]] :> a Impulse[n] ) ]

(*  Simplifications when an exponential is the input  *)

AppendTo[ SPSimplificationRules,
	  ( Rev[n_Symbol][Exp[ Complex[0, a_] b_. n_ ] ] :>
	    Exp[-I a b n] ) ]
AppendTo[ SPSimplificationRules,
	  ( ScaleAxis[l_,w_Symbol][Exp[Complex[0, a_] b_. w_]] :>
	    Exp[I a b w l] ) ]

(*  Other simplifications:  when these are encoded    *)
(*  using ":=", they can cause rules in some of the   *)
(*  packages to confuse Mathematica to no end.  For   *)
(*  example, hard wiring these rules will cause       *)
(*  the "RewriteRules.m" package to hang.             *)
AppendTo[ SPSimplificationRules,
	  ( Aliasby[m_, w_][c_] :> c /; MyFreeQ[c,w] ) ]
AppendTo[ SPSimplificationRules,
	  ( Difference[k_, n_][c_] :> 0 /; MyFreeQ[c,n] ) ]
AppendTo[ SPSimplificationRules,
	  ( Downsample[m_, n_][c_] :> c /; MyFreeQ[c, n] ) ]
AppendTo[ SPSimplificationRules,
	  ( Periodic[k_,n_][x_] :> x /; MyFreeQ[x, n] ) ]
AppendTo[ SPSimplificationRules,
	  ( Rev[n_][x_] :> x /; MyFreeQ[x, n] ) ]
AppendTo[ SPSimplificationRules,
	  ( ScaleAxis[l_, w_][c_] :> c /; MyFreeQ[c, w] ) ]
AppendTo[ SPSimplificationRules,
	  ( Shift[m_, n_][x_] :> x /; MyFreeQ[x, n] ) ]

Protect[ SPSimplificationRules ]


(*  E N D     M E S S A G E  *)

Print["The knowledge base of signal processing operators has been loaded."]
Null

