This is TransSupport.m in view mode; [Download] [Up]
(* :Title: General Platform for Transforms *)
(* :Authors: Brian Evans, James McClellan *)
(*
:Summary: Platform for symbolic transforms
Plots magnitude/phase responses
Generates pole-zero diagrams and root loci
Tools for stability analysis based on transform objects.
*)
(* :Context: SignalProcessing`Support`TransSupport` *)
(* :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: *)
(* :Keywords: pole-zero diagrams, root locus, frequency response, stability *)
(*
:Source: {Multidimensional Digital Signal Processing}, 1984,
by D. Dudgeon and R. Mersereau
{Digital Signal Processing}, 1975, Oppenheim and Schafer
*)
(*
:Warning: The transform objects LTransData and ZTransData must
be defined before this file is loaded (see "ROC.m").
*)
(* :Mathematica Version: 1.2 or 2.0 *)
(* :Limitation: *)
(*
:Discussion: This file has five sections:
A. Routines to aid in the finding of symbolic transforms.
B. Magnitude and phase plots (for 1-D and 2-D signals)
C. Root-locus plots
D. Pole-zero diagrams and root maps (for 1-D and
2-D transforms)
E. Stability analysis which rely on transform objects.
Section B and D display the results of transform objects.
Section E extracts information from transform objects.
*)
(*
:Functions: AddT
ConjT
ConvolveT
DerivativeT
IntegrateT
LineImpulsemDT
MagPhasePlot
MultiDInvTransform
MultiDTransform
MultT
PoleZeroPlot
ROCPlot
RootLocus
ScaleT
ShadedAnnulus
Stable
SubT
SubstituteForT
TransformFixUp
*)
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`TransSupport`",
"SignalProcessing`Support`ROC`",
"SignalProcessing`Support`SigProc`",
"SignalProcessing`Support`SupCode`" ]
(* U S A G E I N F O R M A T I O N *)
AddT::usage =
"AddT[transq, t1, t2] adds the two transforms t1 and t2 together \
and determines the new region of convergence. \
The new transform is returned as a list. \
AddT[transq, t1, t2, lowerlimit] uses lowerlimit as the lower limit \
on region of convergence values when combining ROC's (default is 0)."
ConjT::usage =
"ConjT[transq, X(s), s] implements T{ Conj{x(t)} } --> X*(s*), \
where X(s) is the transform of x(t)."
ConvolveT::usage =
"ConvolveT[transq, convop, t1, t2] implements one dimension of the \
transform of a convolution, Convolution[All, All, t][x1, x2], where \
convop is the convolution operator, t1 is the transform of x1, \
and t2 is the transform of x2."
DerivativeT::usage =
"DerivativeT[transq, expr, s, m] applies the derivative operator \
to the transform expr which is a function of s. \
The derivative is taken m times if and only if expr satisfies transq."
DomainScale::usage =
"DomainScale is an option for MagPhasePlot indicating the \
scaling of the independent variable axis, Linear or Log, \
for the magnitude plot."
IntegrateT::usage =
"IntegrateT[transq, t, variable, lower-limit, upper-limit] integrates \
the transform t (with respect to variable) from lower-limit to \
upper-limit. \
The resulting transform is returned as a list."
Linear::usage =
"Linear is a possible value for axis scaling. See MagPhasePlot."
LineImpulsemDT::usage =
"LineImpulsemDT[transq, expr, s, slist, nleft, fun] applies rules for \
evaluating multidimensional transforms which involve line impulses, \
like the z-transform of f[n1,n2] Impulse[n1 - n2] or the Laplace \
transform of f[t1,t2] Delta[t1 - t2]. \
In either case, the impulse can be removed by setting n1=n2 or t1=t2, \
which yields a one-dimensional function. \
The resulting 1-D transform is then altered by the rule z = z1 z2 \
or s = s1 + s2, respectively. \
Here, the argument fun is Times for the z-transform and Plus for \
the Laplace transform. \
Note that these rules are only applied if expr is a valid transform, \
which occurs if transq[expr] evaluates to True."
MagRangeScale::usage =
"MagRangeScale is an option for MagPhasePlot indicating the \
scaling of the dependent variable axis, Linear or Log, \
for the magnitude plot. \
A setting of Null disables the generation of the magnitude plot."
MagPhasePlot::usage =
"MagPhasePlot[freqresp, {w, wmin, wmax }] plots the magnitude \
and phase response over the specified range of frequencies. \
It returns a list of two elements: the magnitude plot and the \
phase plot. \
The two-dimensional version has the form \
MagPhasePlot[freqresp, {w1, wmin1, wmax1}, {w2, wmin2, wmax2}]. \
The default options are initially biased for continuous-time \
frequency responses of digital signals. \
MagnitudePhasePlot is an alias for MagPhasePlot."
MultiDInvTransform::usage =
"MultiDInvTransform[expression, transvar, timevar, options, transq, \
transform, makeobject, tdefault] finds the multidimensional inverse \
transform of expression by one call to transform per dimension. \
Here, transq, transform, and makeobject are function heads."
MultiDTransform::usage =
"MultiDTransform[makefun, transform, transtest, expr, expr-vars, \
def-trans-var, trans-vars] finds the multidimensional transform \
by calling transform once per dimension. \
If a call to transform produces an incomplete (invalid) transform, \
then the multidimensional transform does not exist and this routine is \
exited. \
The dimension of the transform is determined from \
the number of expr-vars (which equals the number of trans-vars). \
Here, expr is transformed."
MultT::usage =
"MultT[transq, t1, t2] multiplies the transforms t1 and t2 together. \
MultT[transq, t1, t2, lowerlimit] uses lowerlimit as the lower \
limit on the region of convergence (ROC) when combining the ROC's of \
t1 and t2."
PhaseRangeScale::usage =
"PhaseRangeScale is an option for MagPhasePlot indicating the \
scaling of the dependent variable axis, Linear or Log,
for the phase plot. \
A setting of Null disables the generation of the phase plot."
PoleZeroPlot::usage =
"PoleZeroPlot[transform-object] plots the poles and zeros of the \
transform function as well as the region of convergence. \
For the z-transform, this function will also plot the unit circle and \
shade the region of convergence (whenever the ROC does not cover the \
entire z plane). \
To plot a transfer function f, use \
PoleZeroPlot[f, transform-variable, rminus, rplus, zdomainflag]. \
In the two-dimensional case, the transform-variable, rminus, and \
rplus are two-element lists. \
The two-dimensional pole-zero plot has three cases: \
separable transform, symmetric transform, and root map. \
The function returns the poles of the transform."
Radian::usage =
"Radian is a possible value for axis scaling. See MagPhasePlot."
ROCPlot::usage =
"ROCPlot[{rm, rp}] will plot the region of convergence as a \
shaded annular region. \
Possible options for ROCPlot are the same as those for the \
Show (or Plot) command."
RootLocus::usage =
"RootLocus[f(z), z, {freeparam, start, end, step}] plots \
the root locus of f(z) with respect to freeparam over the \
range of start to end evaluated at increments of step. \
RootLocus will generate three plots: the pole zero plot for \
freeparam = start, the root locus, and the pole-zero plot for \
freeparam = end. \
The function returns a list of all three plots. \
RootLocus supports the same options as does Show. \
RootLocusPlot is an alias of RootLocus."
ScaleT::usage =
"ScaleT[transq, t, c] multiplies the transform t by c while leaving \
the region of convergence unaltered. \
The resulting transform is returned as a list."
ShadedAnnulus::usage =
"ShadedAnnulus[rm, rp] and ShadedAnnulus[rm, rp, angulartilt] \
will create a 2-D graphics object that is an annulus \
(rm < radius < rp) filled with uniformly spaced horizontal \
lines rotated by angulartilt degrees."
Stable::usage =
"Stable[transexpr] returns True if the transform transexpr \
represents a stable time-domain signal."
SubT::usage =
"SubT[transq, t1, t2] subtracts the transforms t2 from t1 and \
determines the new region of convergence. \
The resulting transform is returned as a list. \
SubT[transq, t1, t2, lowerlimit] uses lowerlimit as the lower limit \
on region of convergence values when combining ROC's (default is 0)."
SubstituteForT::usage =
"SubstituteForT[transq, t, s, news] substitutes news at every \
occurrence of s in the transform t. \
The resulting transform is returned as a list."
TransformFixUp::usage =
"TransformFixUp[transform, options, time_var, transform_var, \
transform_head, triplet_flag, transform_name, lower, upper] \
attempts to resolve incomplete transforms according to the \
transform pairs provided by the user via the TransformLookup option. \
If the transform is a triplet of information (funct, Rminus, Rplus), \
e.g. the forward z- and Laplace transforms, then triplet_flag is True \
and the lower and upper arguments are used for the second and \
third fields of the triplet."
TransformLookup::usage =
"TransformLookup is an option for the transform rule bases. \
It allows the user to specify a list of transform pairs to augment \
those that already exist in the transform rule bases. \
For example, ZTransform[ Shift[L, n][x[n]], n, z, \
TransformLookup -> { x[n] :> X[z] } ].
This example allows the user to work x[n] without every defining it. \
The multidimensional case is more tricky: \
ZTransform[ x[n1, n2], n, z, TransformLookup -> \
{ x[n1, n2] :> X[z1, n2], X[z1, n2] :> X[z1, z2] } ]. \
For the z- and Laplace transform rule bases, one can also specify \
a region of convergence associated with the transform pairs."
(* E N D U S A G E I N F O R M A T I O N *)
If [ TrueQ[ $VersionNumber >= 2.0 ],
Terms::usage =
"An option for the inverse z- and Laplace transform rule bases, \
InvZTransform and InvLaPlace, which specifies if the inverse \
is to be approximated by a series expansion.";
Protect[ Terms ] ]
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 *)
MagPhasePlot::unresolved =
"These symbols were assigned a value of 1: ``."
MagPhasePlot::nodeltas =
"Delta functions are not shown on the graph"
MultiDInvTransform::notenough =
"Not enough non-transform (time) variables specified."
MultiDTransform::notenough =
"Not enough transform domain variables specified."
PoleZeroPlot::invalidROC =
"The poles `` are inside the specified region of convergence {``, ``}."
PoleZeroPlot::multidimensions =
"Can not produce a pole-zero plot for the multidimensional transform."
PoleZeroPlot::notrational =
"Transform is not a rational polynomial."
PoleZeroPlot::noplot =
"A pole-zero plot cannot be generated."
RootLocus::constant =
"The function `` does not depend on the variable ``."
Transform::badvar =
"The `` variable given to `` is not a symbol or list of symbols: ``"
Transform::incomplete =
"The rule base could not compute the `` of `` with respect to ``."
Transform::integrate =
"Integrating the function `` with respect to ``...."
Transform::novariables =
"No `` variables specified for transform. Possible variables are ``."
Transform::notenough =
"Not enough `` variables specified for the transform."
Transform::twosided =
"The two-sided transform was applied to ``."
TransformLookup::notrule =
"The expression `` passed in the TransformLookup option of the `` \
object is not a rule."
(* A. B A S I C R O U T I N E S F O R T R A N S F O R M S *)
(* AddT *)
AddT[transq_, t1_, t2_, lowerlimit_:0] :=
Transform[ TheFunction[t1] + TheFunction[t2],
FindRMinus[ GetRMinus[t1], GetRMinus[t2], lowerlimit ],
FindRPlus[ GetRPlus[t1], GetRPlus[t2], lowerlimit ] ] /;
transq[t1] && transq[t2]
Format[ AddT[transq_, t1_, t2_, lowerlimit_:0] ] := t1 + t2
(* ConjT *)
ConjT[transq_, t_, s_] :=
Transform[ Conjugate[TheFunction[t]] /. s -> Conjugate[s],
GetRMinus[t], GetRPlus[t] ] /;
transq[t]
Format[ ConjT[transq_, t_, s_] ] :=
SequenceForm[ ToString[StringForm["{ Conj{``} }", t]],
Subscript[s -> news] ]
(* ConvolveT *)
ConvolveT[transq_, convop_, t1_, t2_] :=
Block [ {result},
result = MultT[transq, t1, t2];
Transform[ convop [t1, t2], GetRMinus[t], GetRPlus[t] ] ] /;
transq[t1] && transq[t2]
Format[ ConvolveT[transq_, convop_, t1_, t2_] ] :=
ToString[StringForm["`` ** ``", t1, t2]]
(* DerivativeT *)
DerivativeT[transq_, x_, s_, m_] :=
Transform[ D[TheFunction[x], {s, m}], GetRMinus[x], GetRPlus[x] ] /;
transq[x]
Format[ DerivativeT[transq_, x_, s_, m_] ] :=
SequenceForm[ ColumnForm[{"D" Superscript[m],
" " ~StringJoin~ ToString[s]}],
{ x } ]
(* IntegrateT *)
IntegrateT[transq_, t_, var_, lower_, upper_] :=
Transform[ Integrate[ TheFunction[t], {var, lower, upper} ],
GetRMinus[t], GetRPlus[t] ] /;
transq[t]
Format[ IntegrateT[transq_, t_, var_, l_, u_] ] :=
ToString[StringForm["Integrate[``, {``, ``, ``}]", t, var, l, u]]
(* LineImpulsemDT *)
LineImpulsemDT[transq_, t_, s_, slist_, nleft_, fun_] :=
Block [ {dim, dimsleft, lastrminus, lastrplus, result, rminus, rplus},
rminus = GetRMinus[t];
rplus = GetRPlus[t];
lastrminus = If [ ListQ[rminus], Last[rminus], rminus ];
lastrplus = If [ ListQ[rplus], Last[rplus], rplus ];
dimsleft = Length[nleft];
result = SubstituteForT[transq, t, s, Apply[fun, slist]];
For [ dim = 1, dim <= dimsleft, dim++,
result = ScaleT [ transq, result, SignalProcessing`ROCinfo[nleft[[dim]], lastrminus, lastrplus]] ];
result ] /;
transq[t]
(* MultiDInvTransform *)
MultiDInvTransform[e_, s_, t_, op_, transq_, trans_, makeobject_, tdefault_] :=
Block [ {dim, numdims, result, rminus, roclist,
rplus, sexpr, timedims, tlist},
numdims = Length[s];
timedims = Length[t];
tlist = If [ TrueQ[ numdims > timedims ],
Message[ MultiDInvTransform::notenough ];
DummyVariables[ numdims, tdefault],
Take[t, numdims] ];
If [ transq[e],
rplus = GetRPlus[e];
rminus = GetRMinus[e];
result = TheFunction[e];
For [ dim = 1, dim <= numdims, dim++,
roclist = { result, rminus[[dim]], rplus[[dim]] };
sexpr = makeobject[ roclist, s[[dim]] ];
result = trans[sexpr, s[[dim]], tlist[[dim]], op] ],
result = e;
For [ dim = 1, dim <= numdims, dim++,
result = trans[result, s[[dim]], tlist[[dim]], op] ] ];
result ]
(* MultiDTransform *)
MultiDTransform[makeobjectfun_, transform_, transq_, e_, time_, s_, sl_, op_] :=
Block [ {index, numdims, result, slist },
numdims = Length[time];
result = e;
slist = If [ AtomQ[sl], Table[sl, {index, 1, numdims}], sl ];
If [ Length[slist] < numdims,
Message[MultiDTransform::notenough];
slist = DummyVariables[numdims, s] ];
For [ index = 1, index <= numdims, index++,
result = transform [ op, result, time[[index]],
slist[[index]], time, slist ];
If [ ! TrueQ[transq[result]], Break[] ] ];
makeobjectfun[result] ]
MultiDTransform[makeobjectfun_, transform_, transq_, e_, time_, s_, sl_] :=
Block [ {index, numdims, result, slist },
numdims = Length[time];
result = e;
slist = If [ AtomQ[sl], Table[sl, {index, 1, numdims}], sl ];
If [ Length[slist] < numdims,
Message[MultiDTransform::notenough];
slist = DummyVariables[numdims, s] ];
For [ index = 1, index <= numdims, index++,
result = transform [ result, time[[index]],
slist[[index]], time, slist ];
If [ ! TrueQ[transq[result]], Break[] ] ];
makeobjectfun[result] ]
(* MultT *)
MultT[transq_, t1_, t2_, lowerlimit_:0] :=
Transform[ TheFunction[t1] TheFunction[t2],
FindRMinus[ GetRMinus[t1], GetRMinus[t2], lowerlimit ],
FindRPlus[ GetRPlus[t1], GetRPlus[t2], lowerlimit ] ] /;
transq[t1] && transq[t2]
Format[MultT[transq_, t1_, t2_, lowerlimit_:0]] := t1 t2
(* ScaleT *)
ScaleT[transq_, t_, c_] :=
Transform[ c TheFunction[t], GetRMinus[t], GetRPlus[t] ] /;
transq[t]
ScaleT[transq_, ScaleT[transq_, t_, c_], d_] := ScaleT[transq, t, c d]
Format[ScaleT[transq_, t_, c_]] := c t
(* SubstituteForT *)
SubstituteForT[transq_, t_, s_, news_] :=
Transform[ TheFunction[t] /. s -> news, GetRMinus[t], GetRPlus[t] ] /;
transq[t]
Format[SubstituteForT[transq_, t_, s_, news_]] :=
SequenceForm[ {t}, Subscript[s -> news] ]
(* SubT *)
SubT[transq_, t1_, t2_, lowerlimit_:0] :=
Transform[ TheFunction[t1] - TheFunction[t2],
FindRMinus[ GetRMinus[t1], GetRMinus[t2], lowerlimit ],
FindRPlus[ GetRPlus[t1], GetRPlus[t2], lowerlimit ] ] /;
transq[t1] && transq[t2]
Format[SubT[transq_, t1_, t2_, lowerlimit_:0]] := t1 - t2
(* Transform *)
Transform/: TheFunction[ Transform[fun_, rest___] ] := fun
(* TransformFixUp: resolves functions according to user's transform pairs *)
TransformFixUp[ n_, z_, {}, head_, triplet_, name_, low_, up_ ] := {}
makefixuprules[n_] := { (a_ -> b_) :> ToCollection[{}] /; FreeQ[a, n],
(a_ :> b_) :> ToCollection[{}] /; FreeQ[a, n],
x_ :> x }
TransformFixUp[ n_, z_, options_, head_, triplet_, name_, low_, up_ ] :=
Block [ {fixrules, nolookup, op, rules = {}, validop},
op = Replace[TransformLookup, options];
nolookup = SameQ[op, {}] || SameQ[op, TransformLookup];
If [ ! nolookup,
fixrules = makefixuprules[n];
validop = Map[ Replace[#1, fixrules]&, ToList[op] ];
rules = makepostrules[ validop, n, z, head,
triplet, name, low, up ] ];
rules ]
TransformFixUp[ trans_, n_, z_, options_, head_, triplet_, name_, low_, up_ ] :=
Block [ {rules},
rules = TransformFixUp[ n, z, options, head,
triplet, name, low, up ];
If [ EmptyQ[rules],
trans,
ReplaceRepeated[trans, rules] ] ]
SetAttributes[makepostrules, Listable]
(* rules for transforms returning triplets (forward z- and Laplace) *)
makepostrules[ a_ -> b_Transform, n_, z_, h_, True, transname_, low_, up_ ] :=
( h[ a, n, z, __ ] :> b )
makepostrules[ a_ -> b_, n_, z_, h_, True, transname_, low_, up_ ] :=
( h[ a, n, z, __ ] :> Transform[b, low, up] )
makepostrules[ a_ :> b_Transform, n_, z_, h_, True, transname_, low_, up_ ] :=
( h[ a, n, z, __ ] :> b )
makepostrules[ a_ :> b_, n_, z_, h_, True, transname_, low_, up_ ] :=
( h[ a, n, z, __ ] :> Transform[b, low, up] )
(* rules for transforms not returning triplets *)
makepostrules[ a_ -> b_, n_, z_, h_, False, transname_, low_, up_ ] :=
( h[ a, n, z, __ ] :> b )
makepostrules[ a_ :> b_, n_, z_, h_, False, transname_, low_, up_ ] :=
( h[ a, n, z, __ ] :> b )
(* invalid rule encountered *)
makepostrules[ x_, n_, z_, h_, triplet_, transname_, low_, up_ ] :=
Message[ TransformLookup::notrule, x, transname ]
(* B. F R E Q U E N C Y R E S P O N S E *)
(* getlocation -- Delta functions, a special case of continuous plotting *)
getlocation[ c_. Delta[a_. t_ + b_.], t_ ] := N[ - b / a ]
(* MagPhasePlot *)
MagPhasePlot/: Options[MagPhasePlot] :=
{ DisplayFunction :> $DisplayFunction,
Domain -> Continuous, DomainScale -> Linear,
MagRangeScale -> Linear, PhaseRangeScale -> Degree,
PlotRange -> All }
badoptions[ MagPhasePlot ] =
{ Domain, DomainScale, MagRangeScale, PhaseRangeScale }
(* Supporting local functions *)
fmag[w0_Real, omega_, fresp_, logflag_] :=
Block [ {value},
value = Abs[ GetValue[fresp, omega, w0] ];
If [ logflag, 20 Log[10, value], value ] ]
fmag[w01_Real, omega1_, w02_Real, omega2_, fresp_, logflag_] :=
Block [ {value},
value = Abs[ GetValue[fresp, {omega1, omega2}, {w01, w02}] ];
If [ logflag, 20 Log[10, value], value ] ]
fphase[w0_Real, omega_, fresp_, indegrees_] :=
Block [ {value},
value = GetValue[fresp, omega, w0];
If [ ZeroQ[value], Return[0.] ];
value = N [ Arg[value] ];
Chop[ If [ indegrees, value / Degree, value ] ] ]
fphase[w01_Real, omega1_, w02_Real, omega2_, fresp_, indegrees_] :=
Block [ {value},
value = GetValue[fresp, {omega1, omega2}, {w01, w02}];
If [ ZeroQ[value], Return[0.] ];
value = N [ Arg[value] ];
Chop[ If [ indegrees, value / Degree, value ] ] ]
waxis[w_Real, logflag_] := If [ logflag, Log[10, w], w ]
loglabel[w_] := ToString[StringForm["````", Subscripted[Log[10]], w]]
(* Supporting rule bases for MagPhase Plot *)
(* Interprets signal processing expressions so that they can be plotted *)
ScaleAxisRules = {
ScaleAxis[sc_, w_][x_] :>
Block [ {i, result},
result = x /. w -> sc w;
For [ i = 1, i < sc, i++,
result += x /. w -> (sc w + i 2 Pi) ];
result ]
}
PeriodicRules = {
Periodic[k_, w_Symbol][freqresp_] :>
Block [ {fresp},
fresp = TheFunction[freqresp /. ScaleAxisRules];
fresp + ( fresp /. w -> w + k ) +
( fresp /. w -> w - k ) ],
Periodic[{k1_, k2_}, {w1_Symbol, w2_Symbol}][freqresp_] :>
Block [ {dim1, dim2, fresp, rules},
fresp = TheFunction[freqresp /. ScaleAxisRules];
rules = { w1 -> w1 + dim1 k1,
w2 -> w2 + dim2 k2 };
Apply[Plus,
Table[fresp /. rules,
{dim1, -1, 1}, {dim2, -1, 1}]] ]
}
(* Magnitude and phase plots for 2-D signals (must be defined before 1-D) *)
MagPhasePlot[freqresp_, {w1_Symbol, wmin1_, wmax1_},
{w2_Symbol, wmin2_, wmax2_}, op___] :=
Block [ {arglist, bogus, degrees, dim1, dim2, fresp, ftemp, logrange,
magfun, magplot = NullPlot, omega1, omega2, omitplot, options,
phasefun, phaseplot = NullPlot, plotlabel, plotoptions,
plotrange, result = 0, rules, w1str, w2str},
(* Check for invalid arguments *)
If [ N[wmin1 > wmax1],
Message[Plot::limits, {w1, wmin1, wmax1}]; Return[Null] ];
If [ N[wmin2 > wmax2],
Message[Plot::limits, {w2, wmin2, wmax2}]; Return[Null] ];
(* Set up for plotting *)
Off[Power::infy, Infinity::indt];
w1str = StripPackage[w1];
w2str = StripPackage[w2];
options = ToList[op] ~Join~ Options[MagPhasePlot];
plotoptions = RemoveOptions[options, badoptions[MagPhasePlot]];
(* Adjust frequency response so that plot shows periodicity *)
(* -- need to remove Aliasby operators *)
fresp = freqresp /. PeriodicRules;
(* Extract delta functions from the freq. resp. *)
deltafuns = GetDeltaFunctions[funct, {w1, w2}];
If [ ! SameQ[deltafuns, 0],
Message[ MagPhasePlot::nodeltas ];
fresp = fresp /.
{ h_[ c_. Delta[a_. w1 + b_.] ] :> 0,
h_[ c_. Delta[a_. w1 + b_.] + x_ ] :> h[x],
Delta[a_. w1 + b_.] :> 0,
h_[ c_. Delta[a_. w2 + b_.] ] :> 0,
h_[ c_. Delta[a_. w2 + b_.] + x_ ] :> h[x],
Delta[a_. w2 + b_.] :> 0 } ];
(* Convert expression to plottable/functional form *)
(* and replace constants like Pi with numerical value *)
fresp = N [ TheFunction[fresp] ];
(* Find valueless symbols other than frequency variables *)
(* The first parameter in Summation is a local symbol *)
ftemp = fresp /.
( Summation[n_, l_, u_, inc_][expr_] :>
bogus[l, u, inc, GetVariables[expr /. n -> l]] );
varlist = Select[GetVariables[ftemp],
((! SameQ[#1, w1]) && (! SameQ[#1, w2]))&];
If [ ! EmptyQ[varlist],
Message[ MagPhasePlot::unresolved, varlist];
fresp = fresp /. Map[(#1 -> 1)&, varlist] ];
(* Substitute dummy variables *)
fresp = fresp /. { w1 -> omega1, w2 -> omega2 };
(* Magnitude plotting *)
omitplot = SameQ[Replace[MagRangeScale, options], Null];
logrange = SameQ[Replace[MagRangeScale, options], Log];
plotlabel = If [ logrange,
"Magnitude Response (dB)",
"Magnitude Response" ];
plotrange = Replace[PlotRange, options];
If [ ! omitplot,
arglist = { fmag[w1, omega1, w2, omega2, fresp, logrange],
{ w1, wmin1, wmax1 },
{ w2, wmin2, wmax2 } } ~Join~
plotoptions ~Join~
{ AxesLabel -> { w1str, w2str, " " },
DisplayFunction -> disp,
PlotLabel -> plotlabel,
PlotRange -> plotrange };
magplot = Apply[ Plot3D, arglist ] ];
(* Phase plotting *)
omitplot = SameQ[Replace[PhaseRangeScale, options], Null];
degrees = SameQ[Replace[PhaseRangeScale, options], Degree];
plotlabel = If [ degrees,
"Phase Response (degrees)",
"Phase Response (radians)" ];
If [ ! omitplot,
arglist = { fphase[w1, omega1, w2, omega2, fresp, degrees],
{ w1, wmin1, wmax1 },
{ w2, wmin2, wmax2 } } ~Join~
plotoptions ~Join~
{ DisplayFunction -> disp,
PlotLabel -> plotlabel,
PlotRange -> All,
AxesLabel -> { w1str, w2str, " " } };
phaseplot = Apply[ Plot3D, arglist ] ];
(* Clean up and return the plots as graphics objects *)
On[Power::infy, Infinity::indt];
{ magplot, phaseplot } ]
(* Magnitude and phase plots for 1-D signals *)
MagPhasePlot[ freqresp_, {w_Symbol, wminimum_, wmaximum_}, op___ ] :=
Block [ {arglist, bogus, degrees, deltafuns, deltaplot, disp, fresp,
ftemp, i, list, logdomain, logrange, magfun,
magplot = NullPlot, max, maxdec, min, mindec, omega,
omitplot, options, period, phasefun, phaseplot = NullPlot,
plotlabel, plotoptions, plotrange, points, ticks,
varlist, varlist2, wmin, wmax, wstr, xlabel},
(* Check for errors in arguments *)
If [ N [ wminimum > wmaximum ],
Message[Plot::limits, {w, wminimum, wmaximum}];
Return[Null] ];
(* Set up for plotting *)
Off[Power::infy, Infinity::indt];
options = ToList[op] ~Join~ Options[MagPhasePlot];
wstr = StripPackage[w];
plotoptions = RemoveOptions[options, badoptions[MagPhasePlot]];
(* Adjust frequency response so that plot shows *)
(* periodicity (if any) in the freq. resp. *)
(* -- need to remove Aliasby operators -- *)
fresp = ( freqresp /. PeriodicRules );
(* separate delta funs from rest of freq. resp. *)
deltafuns = GetDeltaFunctions[fresp, w];
fresp = ( fresp /.
{ h_[ c_. Delta[a_. w + b_.] ] :> 0,
h_[ c_. Delta[a_. w + b_.] + x_ ] :> h[x],
Delta[a_. w + b_.] :> 0 } );
(* Convert expression to plottable/functional form *)
(* and replace constants like Pi with numerical value *)
fresp = N [ TheFunction[fresp] ];
(* If frequency response is zero, plot deltas and exit *)
If [ SameQ[fresp, 0],
arglist = Join[ plotoptions,
{ Axes -> axesdefaultvalue,
AxesLabel -> { wstr, " " },
PlotLabel -> "Magnitude Response",
PlotRange -> All,
Ticks -> Automatic } ];
PrependTo[ arglist,
DeltaPlot[ deltafuns, w, N[wminimum],
N[wmaximum], 0.0, 1.0 ] ];
magplot = Apply[ Show, arglist ];
Return[ { magplot, phaseplot } ] ];
(* Rename the frequency variable *)
fresp = fresp /. ( w -> omega );
(* Find valueless symbols other than frequency variable *)
(* The first parameter in Summation is a local symbol *)
ftemp = fresp /.
( Summation[n_, l_, u_, inc_][expr_] :>
bogus[l, u, inc, GetVariables[expr /. n -> l]] );
varlist = Select[GetVariables[ftemp], (! SameQ[#1, omega])&];
If [ ! EmptyQ[varlist],
Message[ MagPhasePlot::unresolved, Sort[ varlist ] ];
fresp = fresp /. Map[(#1 -> 1)&, varlist] ];
(* Set up for plotting *)
logrange = SameQ[Replace[MagRangeScale, options], Log];
plotlabel = If [ logrange,
"Magnitude Response (dB)",
"Magnitude Response" ];
logdomain = SameQ[Replace[DomainScale, options], Log];
If [ N[fresp == 0], logdomain = False ];
xlabel = If [ logdomain, loglabel[wstr], wstr ];
If [ logdomain,
Off[General::dby0];
mindec = Floor[ N[ Log[10, wminimum] ] ];
maxdec = Ceiling[ N[ Log[10, wmaximum] ] ];
wmin = 10^mindec;
wmax = 10^maxdec;
ticks = {};
For [ i = mindec, i < maxdec, i++,
ticks = ticks ~Join~
{i, i + 0.3, i + 0.7} ];
ticks = { ticks, Automatic };
On[General::dby0],
wmin = wminimum;
wmax = wmaximum;
ticks = Automatic ];
(* Magnitude plot with Delta functions (if any) *)
omitplot = SameQ[Replace[MagRangeScale, options], Null];
If [ ! omitplot,
arglist = Join[ { { waxis[w, logdomain],
fmag[w, omega, fresp, logrange] },
{ w, wmin, wmax },
DisplayFunction -> Identity },
plotoptions ];
magplot = Apply[ ParametricPlot, arglist ];
(* determine the range of the magnitude plot *)
plotrange = Replace[PlotRange, options];
list = Map[ Second, magplot[[1]] [[1]] [[1]] [[1]] ];
list = Select[list, NumberQ]; (* remove infinities *)
max = Max[list];
min = Min[list];
If [ logrange,
max = 20 Ceiling[max / 20];
min = 20 Floor[min / 20];
If [ ! SameQ[deltafuns, Null] && max < 20,
max = 20 ];
If [ ( max - min ) > 100,
min = max - 100 ];
If [ SameQ[plotrange, All],
plotrange = {min, max} ] ];
If [ ! logrange && max == 0, max = 1 ];
(* complete magnitude plot -- add in delta functions *)
magplot =
Apply[ Show,
Join[ { DeltaPlot[deltafuns, omega, wmin, wmax,
magplot, min, max] },
plotoptions,
{ PlotLabel -> plotlabel,
PlotRange -> plotrange,
Ticks -> ticks,
AxesLabel -> { xlabel, " " } } ] ] ];
(* Phase plot *)
omitplot = SameQ[Replace[PhaseRangeScale, options], Null];
If [ ! omitplot,
degrees = SameQ[Replace[PhaseRangeScale, options], Degree];
plotlabel = If [ degrees,
"Phase Response (degrees)",
"Phase Response (radians)" ];
arglist = Join[ { { waxis[w, logdomain],
fphase[w, omega, fresp, degrees] },
{ w, wmin, wmax } },
plotoptions,
{ PlotLabel -> plotlabel,
PlotRange -> All,
Ticks -> ticks,
AxesLabel -> { xlabel, " " } } ];
phaseplot = Apply[ ParametricPlot, arglist ] ];
(* Clean up and return the two plots as graphics objects *)
On[Power::infy, Infinity::indt];
{ magplot, phaseplot } ]
(* C. R O O T - L O C U S P L O T T I N G *)
polyorder[poly_, z_] := Exponent[poly, z] + Exponent[poly /. z -> z^-1, z]
Options[RootLocus] :=
{ Axes -> axesdefaultvalue, AxesLabel -> { "Re", "" },
DisplayFunction :> $DisplayFunction, PlotRange -> All }
getplotrange[ xvalues_ ] :=
Block [ {offset, xmin, xmax},
xmin = Min[xvalues];
xmax = Max[xvalues];
offset = ( xmax - xmin ) / 4; (* always >= 0 *)
xmin -= offset;
xmax += offset;
{ xmin, xmax } ]
localplot[ label_, oplist_, plots___ ] :=
Apply[ Show, Join[ { plots }, { PlotLabel -> label }, oplist ] ]
RootLocus[ poly_, z_Symbol, {k_Symbol, start_, end_, step_:1}, op___ ] :=
MyMessage[ RootLocus::constant,
{ NullPlot, NullPlot, NullPlot }, poly, z ] /;
FreeQ[poly, z]
RootLocus[ poly_, z_Symbol, {k_Symbol, start_, end_, step_:1}, op___ ] :=
Block [ {Epplot = NullPlot, Ezplot = NullPlot, Szplot = NullPlot,
Spplot = NullPlot, allcoords = {}, coords, denom, factor,
kmax, kmin, kstep, ktemp, nend, newpoly, nstart, numer,
options, order = 1, pointsize, pointthickness,
poles = {}, pplot, rootlist = {}, xmin, xmax,
ymin, ymax, zeroes = {}, zplot},
(* Determine the numeric range for the free param k *)
nstart = N[start];
nend = N[end];
options = Flatten[ToList[op]] ~Join~ Options[RootLocus];
kmin = Min[nstart, nend];
kmax = Max[nstart, nend];
kstep = Abs[N[step]];
(* Put the function over a common denominator *)
newpoly = Together[poly] /. k -> ktemp;
(* Get numerator and denominator-- expand them for speed *)
(* note: expand also hides errors in Solve under Mma 1.2 *)
numer = Expand[Numerator[newpoly]];
denom = Expand[Denominator[newpoly]];
(* Compute the thickness of the points in the locus *)
If [ MixedPolynomialQ[poly, z],
order = polyorder[Expand[numer], z] +
polyorder[Expand[denom], z] ];
factor = N[ Log[10, order (kmax - kmin) / kstep] ];
pointthickness = Max[0.024 / factor, 0.001];
pointsize = PointSize[pointthickness];
(* Compute the initial and final zero plot *)
If [ ! FreeQ[numer, z],
rootlist = GetRootList[ numer /. ktemp -> kmin, z ];
coords = ComplexTo2DCoordList[rootlist];
allcoords = Join[allcoords, coords];
Szplot = PointwisePlot[ coords, "O", "*O*" ];
rootlist = GetRootList[ numer /. ktemp -> kmax, z ];
coords = ComplexTo2DCoordList[rootlist];
allcoords = Join[allcoords, coords];
Ezplot = PointwisePlot[ coords, "O", "*O*" ] ];
(* Compute the initial and final pole plot *)
If [ ! FreeQ[denom, z],
rootlist = GetRootList[ denom /. ktemp -> kmin, z ];
coords = ComplexTo2DCoordList[rootlist];
allcoords = Join[allcoords, coords];
Spplot = PointwisePlot[ coords, "X", "*X*" ];
rootlist = GetRootList[ denom /. ktemp -> kmax, z ];
coords = ComplexTo2DCoordList[rootlist];
allcoords = Join[allcoords, coords];
Epplot = PointwisePlot[ coords, "X", "*X*" ] ];
(* Compute the zero locus *)
zplot = If [ FreeQ[numer, ktemp],
NullPlot,
zeroes =
Table[ ToCollection[GetRootList[numer, z]],
{ktemp, kmin, kmax, kstep} ];
coords = ComplexTo2DCoordList[zeroes];
Clear[zeroes];
allcoords = Join[allcoords, coords];
coords = Map[Point, coords];
PrependTo[coords, pointsize];
Graphics[coords] ];
(* Compute the pole locus *)
pplot = If [ FreeQ[denom, ktemp],
NullPlot,
poles =
Table[ ToCollection[GetRootList[denom, z]],
{ktemp, kmin, kmax, kstep} ];
coords = ComplexTo2DCoordList[poles];
Clear[poles];
allcoords = Join[allcoords, coords];
coords = Map[Point, coords];
PrependTo[coords, pointsize];
Graphics[coords] ];
(* Determine the true plot range *)
If [ SameQ[Replace[PlotRange, options], All],
{xmin, xmax} = getplotrange[ Map[First, allcoords] ];
{ymin, ymax} = getplotrange[ Map[Second, allcoords] ];
options = Prepend[ options,
PlotRange -> {{xmin, xmax},
{ymin, ymax}} ] ];
(* Free up memory before plots are combined *)
Clear[allcoords];
Clear[coords];
{ localplot[ "Pole-Zero Diagram at Start", options,
Szplot, Spplot],
localplot[ "Root Locus", options, Szplot, Spplot,
Ezplot, Epplot, zplot, pplot ],
localplot[ "Pole-Zero Diagram at End", options,
Ezplot, Epplot ] } ]
(* D. P O L E - Z E R O P L O T T I N G *)
(* 1. S u p p o r t i n g F u n c t i o n s *)
(* a. For graphing the region of convergence as an annulus *)
(* findline *)
findline[ycoord_, cos_, sin_] :=
Block [ {xcoord},
xcoord = N[ Sqrt[1 - ycoord^2] ];
{ rotate[xcoord, ycoord, cos, sin],
rotate[-xcoord, ycoord, cos, sin] } ]
makeannulus[{rm_, rp_}] :=
ToCollection[ CirclePS[rp, Dashing[{0.05,0.05}]],
ShadedAnnulus[rm, rp],
CirclePS[rm, Dashing[{0.05,0.05}]] ]
(* makefilllines *)
(* Fill lines for a unit circle: tilt in degrees, *)
(* and angular separation ~= 360 / numsamples *)
makefilllines[tilt_, numsamples_] :=
Block [ {costilt, maxi, sintilt},
costilt = N[ Cos[ tilt Pi / 180 ] ];
sintilt = N[ Sin[ tilt Pi / 180 ] ];
maxi = Floor[numsamples / 4] - 1;
Join[ Table[ findline[4 i / numsamples, costilt, sintilt],
{i, 0, maxi} ],
Table[ findline[-1 + 4 i / numsamples, costilt, sintilt],
{i, 1, maxi} ] ] ]
(* rotate *)
rotate[x_, y_, cos_, sin_] := { x cos - y sin, x sin + y cos }
(* splitline *)
splitline[ap_, rm_, slope_, {point1_, point2_}] :=
Block [ {b, sqrtterm, x1, x2, y1, y2},
b = Second[point1] - slope First[point1];
sqrtterm = Re[ N[ Sqrt[ rm^2 ap - b^2 ] ] ];
x1 = ( - slope b + sqrtterm ) / ap;
x2 = ( - slope b - sqrtterm ) / ap;
y1 = slope x1 + b;
y2 = slope x2 + b;
{{point1, {x1, y1}}, {{x2, y2}, point2}} ]
(* ShadedAnnulus *)
AngularTilt = 20 (* default values *)
NumSamples = 72
FillLines = makefilllines[AngularTilt, NumSamples]
FillLineSlope = N[ Tan[AngularTilt Pi / 180] ] (* slope of all fill lines *)
ShadedAnnulus[rm_, rp_] :=
ShadedAnnulus[ rm, rp, AngularTilt, FillLines, FillLineSlope ]
ShadedAnnulus[rm_, rp_, tilt_] :=
ShadedAnnulus[ rm, rp, tilt, makefilllines[tilt, NumSamples] ]
ShadedAnnulus[rm_, rp_, tilt_, filllines_] :=
ShadedAnnulus[ rm, rp, tilt, filllines, N[Tan[tilt Pi / 180]] ]
ShadedAnnulus[0, rp_, tilt_, filllines_, slope_] :=
Graphics[ Map[Line, rp filllines] ]
ShadedAnnulus[rm_, rp_, tilt_, filllines_, slope_] :=
Block [ {ap, fillpattern, i, intersections,
maxi, newlines, numlines, split, xnew},
fillpattern = (rp filllines);
numlines = Length[filllines];
(* Find intersection of lines with rminus circle *)
(* Equations: x^2 + y^2 = r^2 and y = m x + b *)
(* reduces to quadratic a' x^2 + b' x + c', with *)
(* a' = (m^2 + 1), b' = 2 m b, c' = b^2 - r^2 *)
ap = ( slope^2 + 1 );
(* Intersection of first line, zero intercept *)
xnew = N[ - rm / Sqrt[ap] ];
newlines = { { {xnew, slope xnew},
Second[ First[fillpattern] ] } };
fillpattern[[1]] = { First[ First[fillpattern] ],
{-xnew, - slope xnew} };
(* Find intersections of other lines *)
intersections = ( numlines + 2 ) rm / rp;
maxi = Floor[ intersections / 2];
For [ i = 1, i <= maxi, i++,
split = splitline[ ap, rm, slope, fillpattern[[i + 1]] ];
AppendTo[ newlines, First[split] ];
fillpattern[[i + 1]] = Second[split];
split = splitline[ ap, rm, slope,
fillpattern[[numlines - i + 1]] ];
AppendTo[ newlines, First[split] ];
fillpattern[[numlines - i + 1]] = Second[split] ];
Graphics[ Map[Line, Join[fillpattern, newlines]] ] ]
(* ROCPlot *)
Options[ROCPlot] :=
{ Axes -> axesdefaultvalue, AspectRatio -> 1,
DisplayFunction :> $DisplayFunction }
ROCPlot[{rm_, Infinity}, options___] :=
Block [ {oplist, plotrange, pmax, rmax},
oplist = ToList[options];
plotrange = Replace[PlotRange, oplist];
pmax = If [ SameQ[plotrange, PlotRange],
N[ Sqrt[2] rm ],
Max[ N[Sqrt[2] rm],
Abs[Select[N[Flatten[plotrange]], NumberQ]] ] ];
rmax = 1.75 pmax; (* pmax times some number > Sqrt[2] *)
Apply [ Show,
Join [ {makeannulus[{rm, rmax}]},
oplist,
Options[ROCPlot],
PlotRange -> {{-pmax, pmax}, {-pmax, pmax}} ] ] ]
ROCPlot[{rm_, rp_}, options___] :=
Apply[ Show,
Join[{makeannulus[{rm, rp}]},
ToList[options],
Options[ROCPlot]] ]
(* b. other supporting functions *)
(* InvalidZPoleQ *)
InvalidZPoleQ[rminus_, r_, rplus_] := TrueQ[rminus < Abs[N[r]] < rplus]
(* InvalidSPoleQ *)
InvalidSPoleQ[rminus_, r_, rplus_] := TrueQ[rminus < Re[N[r]] < rplus]
(* MyRationalPolyQ *)
MyRationalPolyQ[f_, z_] :=
MixedPolynomialQ[ Expand[Numerator[f]], z ] &&
MixedPolynomialQ[ Expand[Denominator[f]], z ]
(* pzOptions *)
pzOptions[ op___ ] :=
Join[ ToList[op],
ToList[ Options[PoleZeroPlot] ],
{ DisplayFunction :> $DisplayFunction } ]
(* SeparableQ *)
Separate[x_ y_, {var1_, var2_}] := { x, y } /; FreeQ[x, var2] && FreeQ[y, var1]
Separable[x_ y_, {var1_, var2_}] := True /; FreeQ[x, var2] && FreeQ[y, var1]
SeparableQ[p_, vars_] := TrueQ[ Separable[p, vars] ]
(* 2. D r i v e r *)
Options[ PoleZeroPlot ] :=
{ Dialogue -> True, DisplayFunction :> $DisplayFunction }
(* One-dimensional pole-zero plotting for z-transform objects *)
PoleZeroPlot[ ZTransData[f_, Rminus[rm_], Rplus[rp_], ZVariables[z_]], op___ ] :=
If [ Length[z] <= 2,
PoleZeroPlot[ f, z, rm, rp, True, op ],
Message[PoleZeroPlot::multidimensions] ]
(* One-dimensional pole-zero plotting for Laplace transform objects *)
PoleZeroPlot[ LTransData[f_, Rminus[rm_], Rplus[rp_], LVariables[s_]], op___ ] :=
If [ Length[s] <= 2,
PoleZeroPlot[ f, s, rm, rp, False, op ],
MyMessage[PoleZeroPlot::multidimensions] ]
(* Two-dimensional pole-zero plotting driver *)
(* -- separable transforms generate two separable pole-zero plots *)
(* -- symmetric transfer function f only requires one plot *)
(* -- otherwise, project z1 onto z2 and z2 onto z1 *)
PoleZeroPlot[f_, {z1_Symbol, z2_Symbol}, {rm1_, rm2_}, {rp1_, rp2_}, rest__ ] :=
Block [ {seplist, z1expr, z2expr},
If [ SeparableQ[f, {z1, z2}],
Block [ {poles, seplist, z1expr, z2expr},
seplist = Separate[f, {z1, z2}];
z1expr = seplist[[1]];
z2expr = seplist[[2]];
poles = PoleZeroPlot[z1expr, z1, rm1, rp1, rest];
poles = poles ~Join~
PoleZeroPlot[z2expr, z2, rm2, rp2, rest];
poles ],
If [ ! SameQ[f, f /. { z1 -> z2, z2 -> z1 }],
PoleZeroPlot[f, z2, z1, rm1, rp1, rest] ];
PoleZeroPlot[f, z1, z2, rm2, rp2, rest] ] ]
PoleZeroRootLocus[ ztrans_, z_Symbol, {w_Symbol, wmin_, wmax_, wstep_},
ucROCplot_, options_ ] :=
Block [ {denom, display, numer, polelist, rootmap, xlabel, ylabel,
zerolist, zstr},
numer = Numerator[ztrans];
denom = Denominator[ztrans];
display = Replace[DisplayFunction, options];
(* Get zeros and poles as expressions of zproj *)
zerolist = GetRootList[numer, z];
polelist = GetRootList[denom, z];
If [ InformUserQ[options],
Print[ " " ];
Print[ "The zeroes are: ", zerolist /. w -> Global`w ];
Print[ "The poles are: ", polelist /. w -> Global`w ] ];
(* Root map labels *)
zstr = StripPackage[z];
xlabel = "Re " ~StringJoin~ zstr;
ylabel = "Im " ~StringJoin~ zstr;
(* Zero root map *)
If [ ! FreeQ[numer, z],
rootmap = RootLocus[ numer, z,
{w, wmin, wmax, wstep},
AspectRatio -> Automatic,
DisplayFunction :> Identity,
PlotLabel -> "Zero Root Map",
AxesLabel -> { xlabel, ylabel } ];
Show[ { ucROCplot, rootmap[[2]] },
DisplayFunction -> display ] ];
(* Pole root map *)
If [ ! FreeQ[denom, z],
rootmap = RootLocus[ 1 / denom, z,
{w, wmin, wmax, wstep},
AspectRatio -> Automatic,
DisplayFunction :> Identity,
PlotLabel -> "Pole Root Map",
AxesLabel -> { xlabel, ylabel } ];
Show[ { ucROCplot, rootmap[[2]] },
DisplayFunction -> display ] ];
polelist ]
(* 3. P l o t t i n g R o u t i n e s f o r z - d o m a i n *)
PoleZeroPlot[f_, z_Symbol, rm_, rp_, True, op___] :=
Block [ {abspolelist, display, fillplot, invalidROClist, options,
polelist, poleplot, rmax, rminus, rplus, rminuscircle,
rpluscircle, unitcircle, xlabel, ylabel, zeroplot,
zstr, ztrans},
options = pzOptions[op];
display = Replace[DisplayFunction, options];
rmax = 1; (* default values *)
unitcircle = CirclePS[1];
rminus = N[rm]; (* numeric approximation *)
rplus = N[rp];
(* Analyze the z-transform function *)
ztrans = Together[f];
If [ ! MyRationalPolyQ[ztrans, z],
Message[ PoleZeroPlot::notrational ];
Message[ PoleZeroPlot::noplot ];
Return[{}] ];
(* Find the poles and zeroes *)
zerolist = GetRootList[Numerator[ztrans], z];
polelist = GetRootList[Denominator[ztrans], z];
If [ InformUserQ[options],
Print[ " " ];
Print[ "The zeroes are: ", zerolist ];
Print[ "The poles are: ", polelist ] ];
(* Check consistency of poles with region of convergence *)
(* Exit function if inconsistency encountered *)
invalidROClist = Select[ polelist,
InvalidZPoleQ[rminus, #1, rplus]& ];
If [ ! SameQ[invalidROClist, {}],
Message[PoleZeroPlot::invalidROC, invalidROClist, rm, rp];
Return[polelist] ];
(* Compute plot of poles and zeroes *)
(* Determine the plot's its circular extent *)
rmax = Max[ Join[ {rmax},
Map[Abs, ToList[N[zerolist]]],
Map[Abs, ToList[N[polelist]]] ] ];
zeroplot = If [ EmptyQ[zerolist],
NullPlot,
PointwisePlot[ComplexTo2DCoordList[zerolist],
"O", "*O*"] ];
poleplot = If [ EmptyQ[polelist],
NullPlot,
PointwisePlot[ComplexTo2DCoordList[polelist],
"X", "*X*"] ];
(* Build region of convergence plot *)
If [ NumberQ[rminus],
rminuscircle = CirclePS[ rminus, Dashing[{0.05,0.05}] ];
rmax = Max[rmax, Abs[rminus]],
rminuscircle = NullPlot ];
If [ NumberQ[rplus],
rpluscircle = CirclePS[ rplus, Dashing[{0.05,0.05}] ];
rmax = Max[rmax, Abs[rplus]],
rpluscircle = NullPlot ];
rplus = N[ If[ SameQ[rp, Infinity], 2 rmax, rp] ]; (* kludge *)
fillplot = If [ NumberQ[rminus] && NumberQ[rplus],
ShadedAnnulus[rminus, rplus],
NullPlot ];
(* Send the contour plot to the screen *)
rmax *= 1.25;
zstr = StripPackage[z];
xlabel = "Re " ~StringJoin~ zstr;
ylabel = "Im " ~StringJoin~ zstr;
Show [unitcircle, rminuscircle, rpluscircle,
fillplot, zeroplot, poleplot,
PlotRange -> {{-rmax, rmax}, {-rmax, rmax}},
AspectRatio -> Automatic,
AxesLabel -> { xlabel, ylabel },
Axes -> axesdefaultvalue,
DisplayFunction -> display ];
polelist ]
(* Root map two-dimensional pole-zero plotting *)
(* variable zproj is projected onto z-domain *)
PoleZeroPlot[f_, zproj_Symbol, z_Symbol, rm_, rp_, True, op___] :=
Block [ {graphfuns, options, plotstyle,
rminus, rplus, ucROCplot, w, ztrans},
options = pzOptions[op] ~Join~ { AspectRatio -> 1 };
(* Make sure that the transform is over common denominator *)
ztrans = Together[f];
If [ ! MyRationalPolyQ[ztrans, z],
Message[ PoleZeroPlot::notrational ];
Message[ PoleZeroPlot::noplot ];
Return[{}] ];
(* We will plot a unit circle and R- and possibly R+ *)
rminus = N[rm /. zproj -> w];
rplus = N[rp /. zproj -> w];
plotstyle = { Thickness[0.01] };
graphfuns = { { Cos[w], Sin[w] } };
If [ ! ZeroQ[rminus],
PrependTo[ graphfuns, { rminus Cos[w], rminus Sin[w] } ];
PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
If [ ! InfinityQ[rplus],
PrependTo[ graphfuns, { rplus Cos[w], rplus Sin[w] } ];
PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
ucROCplot = Apply[ ParametricPlot,
{ graphfuns,
{w, 0, 2 Pi},
AspectRatio -> 1,
DisplayFunction :> Identity,
PlotStyle -> plotstyle } ];
(* Project the transfer function and draw the root locus *)
ztrans = ztrans /. zproj -> Exp[I w];
PoleZeroRootLocus[ ztrans, z, {w, 0, 2 Pi, 2 Pi / 360},
ucROCplot, options ] ]
(* 4. P l o t t i n g R o u t i n e s f o r s - d o m a i n *)
(* One-dimensional pole-zero plotting *)
PoleZeroPlot[f_, s_Symbol, rm_, rp_, False, op___] :=
Block [ {display, invalidROClist, jj, ltrans, numericlist, options,
points, rminus, rminusline, rplus, rplusline, polelist,
poleplot, shading, skew, sstr, x1, x2, xlabel, xmin,
xmax, w, ylabel, ymin, ymax, ystep, zerolist, zeroplot},
options = pzOptions[op];
display = Replace[DisplayFunction, options];
ltrans = Together[f]; (* default values *)
rminusline = rplusline = shading = NullPlot;
rminus = N[rm]; (* numeric approximation *)
rplus = N[rp];
If [ ! RationalPolynomialQ[ltrans, s],
Message[ PoleZeroPlot::notrational ];
Message[ PoleZeroPlot::noplot ];
Return[{}] ];
zerolist = GetRootList[Numerator[ltrans], s];
polelist = GetRootList[Denominator[ltrans], s];
If [ InformUserQ[options],
Print[ " " ];
Print[ "The zeroes are: ", zerolist ];
Print[ "The poles are: ", polelist ] ];
(* Check consistency of poles with region of convergence *)
(* Exit function if inconsistency encountered *)
invalidROClist = Select[ polelist,
InvalidSPoleQ[rminus, #1, rplus]& ];
If [ ! SameQ[invalidROClist, {}],
Message[PoleZeroPlot::invalidROC, invalidROClist, rm, rp];
Return[polelist] ];
numericlist = ToList[N[zerolist]] ~Join~ ToList[N[polelist]];
(* Determine extent of pole-zero plot *)
(* Mathematica does not do this automatically *)
xmin = Min[Re[numericlist]];
xmax = Max[Re[numericlist]];
If [ NumberQ[rminus], xmin = Min[xmin, rminus] ];
If [ NumberQ[rplus], xmax = Max[xmax, rplus] ];
xmin = ( 1.0 - 0.25 Sign[xmin] ) xmin; (* move toward -Inf *)
xmax = ( 1.0 + 0.25 Sign[xmax] ) xmax; (* move toward +Inf *)
If [ ZeroQ[xmax], xmax = -.1 xmin ];
ymin = Min[Im[numericlist]];
ymax = Max[Im[numericlist]];
If [ SameQ[ymin, ymax],
ymin = ymin - 1;
ymax = ymax + 1,
ymin = ( 1.0 - 0.25 Sign[ymin] ) ymin;
ymax = ( 1.0 + 0.25 Sign[ymax] ) ymax ];
(* Build convergence boundaries (lines) *)
If [ NumberQ[rminus],
rminusline = Graphics[ { Dashing[{0.05, 0.05}],
Line[{{rminus, ymin},
{rminus, ymax}}] } ] ];
If [ NumberQ[rplus],
rplusline = Graphics[ { Dashing[{0.05, 0.05}],
Line[{{rplus, ymin},
{rplus, ymax}}] } ] ];
If [ (NumberQ[rminus] || SameQ[rm, -Infinity]) &&
(NumberQ[rplus] || SameQ[rp, Infinity]),
x1 = Max[rminus, xmin];
x2 = Min[rplus, xmax];
ystep = ( ymax - ymin ) / 15;
skew = 4 ystep;
points = Table[ Line[{{x1,jj}, {x2,jj+skew}}],
{jj, ymin - skew, ymax, ystep} ];
shading = Graphics[ points ] ];
(* Build plots of zeroes as O's and poles as X's *)
If [ EmptyQ[zerolist],
zeroplot = NullPlot,
zeroplot = PointwisePlot[ ComplexTo2DCoordList[zerolist],
"O", "*O*" ] ];
If [ EmptyQ[polelist],
poleplot = NullPlot,
poleplot = PointwisePlot[ ComplexTo2DCoordList[polelist],
"X", "*X*" ] ];
(* Plot pole-zero diagram, all together now *)
sstr = StripPackage[s];
xlabel = "Re " ~StringJoin~ sstr;
ylabel = "Im " ~StringJoin~ sstr;
Show [ rminusline, rplusline, shading, zeroplot, poleplot,
PlotRange -> {{xmin, xmax}, {ymin, ymax}},
AxesLabel -> { xlabel, ylabel },
Axes -> axesdefaultvalue,
DisplayFunction -> display ];
polelist
]
(* Root map two-dimensional pole-zero plotting *)
(* variable sproj is projected onto s-domain *)
PoleZeroPlot[f_, sproj_Symbol, s_Symbol, rm_, rp_, False, op___] :=
Block [ {graphfuns, ltrans, options, plotstyle, rminus,
rplus, ucROCplot, ucROCplot1, ucROCplot2, w},
options = pzOptions[op];
(* Make sure the transfer function is a ratio of two polys *)
ltrans = Together[f];
If [ ! MyRationalPolyQ[ltrans, s],
Message[ PoleZeroPlot::notrational ];
Message[ PoleZeroPlot::noplot ];
Return[{}] ];
(* We will possibly plot R- and R+ *)
rminusline = NullPlot;
rplusline = NullPlot;
rminus = N[rm /. sproj -> (I w)];
rplus = N[rp /. sproj -> (I w)];
plotstyle = { };
graphfuns = { };
If [ ! InfinityQ[rminus],
PrependTo[ graphfuns, {rminus, w} ];
PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
If [ ! InfinityQ[rplus],
PrependTo[ graphfuns, {rplus, w} ];
PrependTo[ plotstyle, Dashing[{0.05,0.05}] ] ];
ucROCplot1 = Apply[ ParametricPlot,
{ graphfuns,
{w, -10, 10},
DisplayFunction :> Identity,
PlotStyle -> plotstyle } ];
ucROCplot2 = Apply[ ParametricPlot,
{ graphfuns,
{w, -10000, 10000},
DisplayFunction :> Identity,
PlotStyle -> plotstyle } ];
ucROCplot = Show[ {ucROCplot1, ucROCplot2},
DisplayFunction :> Identity ];
(* Project the transfer function and draw the root locus *)
ltrans = ltrans /. sproj -> (I w);
PoleZeroRootLocus[ ltrans, s, {w, -10000, 10000, 50},
ucROCplot, options ] ]
(* E. S T A B I L I T Y A N A L Y S I S *)
(*
The trick here for m-D transforms is to replace any occurrence of
the transform variables in the ROC with the value of the boundary at
which the ROC becomes unstable. For the Laplace transform, this is
the line Re(s) = 0. For the z-transform, this is the m-D unit circle;
because the ROC in the z-transform are real numbers (absolute values
of complex-valued quantities), we let z1 = 1, z2 = 1, etc., instead
of z1 = Exp[I w1], z2 = Exp[I w2], etc.
*)
unknown[ condexpr_ ] :=
Block [ {simplified},
simplified = Simplify[condexpr] /. SPLessGreaterRules;
If [ N[simplified], True, False, simplified ] ]
LTransData/: Stable[ LTransData[t_, rm_, rp_, v_] ] :=
Block [ {inrange, vars, ltrans},
ltrans = LTransData[t, rm, rp, v];
inrange = InRange[GetRMinus[ltrans], 0, GetRPlus[ltrans],
-Infinity, Infinity, Less, Less];
vars = First[v];
If [ ListQ[vars], (* evaluate along imag. axes *)
inrange = inrange /. ( Map[Re, vars] ~ReplaceWith~ 0 ) ];
If [ inrange, True, False, unknown[inrange] ] ]
ZTransData/: Stable[ ZTransData[t_, rm_, rp_, v_] ] :=
Block [ {inrange, vars, ztrans},
ztrans = ZTransData[t, rm, rp, v];
inrange = InRange[GetRMinus[ztrans], 1, GetRPlus[ztrans],
0, Infinity, Less, Less];
vars = First[v];
If [ ListQ[vars], (* evaluate along unit circles *)
inrange = inrange /. ( vars ~ReplaceWith~ 1 ) ];
If [ inrange, True, False, unknown[inrange] ] ]
(* E N D P A C K A G E *)
End[]
EndPackage[]
(* A L I A S E S *)
Unprotect[ MagnitudePhasePlot ]
Clear[ MagnitudePhasePlot ]
MagnitudePhasePlot = MagPhasePlot
MagnitudePhasePlot::usage = MagPhasePlot::usage
Protect[ MagnitudePhasePlot ]
Unprotect[ RootLocusPlot ]
Clear[ RootLocusPlot ]
RootLocusPlot = RootLocus
RootLocusPlot::usage = RootLocus::usage
Protect[ RootLocusPlot ]
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 =
{ AddT, ConjT, ConvolveT,
DerivativeT, IntegrateT, LineImpulsemDT,
MagPhasePlot, MultiDInvTransform, MultiDTransform,
MultT, PoleZeroPlot, RootLocus,
ScaleT, Stable, SubstituteForT,
SubT, TransformFixUp };
Combine[ SPfunctions, newfuns ];
Apply[ Protect, newfuns ] ]
(* E N D I N G M E S S A G E *)
Print["Supporting routines for transform rule bases are loaded."]
Null
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.