View Raw SPL
/*****************************************************************************
* *
* SPLANE.SPL Copyright (C) 2015 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Pole-Zero plot of an S transform *
* *
* Revisions: 18 Jun 2015 RRR Creation - from zplane.spl *
* *
*****************************************************************************/
#if @HELP_SPLANE
SPLANE
Purpose: Displays a Pole-Zero plot of an S transform.
Syntax: SPLANE(b, a)
b - A series, the numerator coefficients in descending
powers of s.
a - A series, the denominator coefficients in descending
powers of s.
Alternate Syntax: SPLANE(z, p, gain)
z - A series, the zeros of the S transform.
p - A series, the poles of the S transform.
gain - A scalar, the gain of the system.
Alternate Syntax: SPLANE(c)
c - A series, the system coefficients in cascaded
biquad form. If c contains 2 columns, c is
assumed to be in direct form, where the first
column is b and the second column is a.
Returns: An XY series, the poles are plotted as light red "x"s and the
zeros are plotted as "o"s. The unit circle is displayed as
a dashed circle.
Example:
s
Given H(s) = -------
s + 0.5
W1: splane({1}, {s, 0.5})
Displays the s plane with a single pole at s = -0.5.
The input is given as a system function in descending terms
of s.
Example:
W1: splane({0}, {-0.5}, 1)
Displays the same plot as above. The input is given in terms
of zeros and poles.
Example:
s^2 + 2s
Given H(s) = ----------------
s^2 + 0.7s + 0.1
z = roots({1, 2, 0})
p = roots({1, 0.7, 0.1})
splane(z, p, 1)
Displays two real poles and two zeros in the current Window.
Remarks:
For SPLANE(b, a), the input series represent the
terms of the rational polynomial H(s) = b(s)/a(s) where:
N = length(b) and M = length(a):
b(s) b[1] s^(N-1) + b[1] s^(N-1} + ... + b[N]
H(s) = ------ = ----------------------------------------
a(s) a[1] s^(M-1) + a[1] s^(M-1} + ... + a[M]
For SPLANE(z, p, gain), the gain term must be present, but it
does not effect the resulting plot.
The aspect ratio of the window is not altered.
The poles are plotted as light red x's and the zeros as black
circles. Solid lines are drawn through the origin.
Multiple zeros and poles are labeled with a multiplicity
number to the upper right of the symbol.
For SPLANE(c), if the input c is a single column, the
coefficients are assumed to be in cascaded bi-quad form.
This is the output format of IIR filters designed by
DADiSP/Filters. If c contains 2 columns, the coefficients
are assumed to be in direct form, where the first column
is b and the second column is a.
See Also:
Residue
Sfreq
Zplane
#endif
static splane_pmode = -1;
/* S domain pole-zero plot */
splane(b, a, gain)
{
local z, p, t, x, y, tol, ymin, ymax, xmin, xmax, h, ext;
local pxy, zxy, i, v1, v2, fac, len, m, item, lstyle;
if (argc < 1)
{
error("splane - input coefficient series requied");
}
/* we do not change the aspect ratio */
if (argc == 3)
{
/* zeros, poles, gain */
z = {b};
p = {a};
}
else
{
z = p = {};
if (argc < 2)
{
if (numcols(b) == 2)
{
/* raveled B(z), A(z) */
a = col(b, 2);
b = col(b, 1);
}
else
{
/* assume cascaded bi-quad coefficients ala Filters Module */
(z, p, k) = cas2zp(b);
a = {};
b = {};
}
}
else
{
/* cast to series so we can enter a scalar */
a = {a};
b = {b};
}
if (length(a) > 0 || length(b) > 0)
{
(z, p, k) = tf2zp(b, a);
}
}
/* plotting off */
splane_pmode = plotmode();
plotmode(0);
/* convert zeros and poles to XY series */
if (length(z) > 0)
{
zxy = xy(real(z), imag(z));
setcomment(zxy, "Zeros");
}
if (length(p) > 0)
{
pxy = xy(real(p), imag(p));
setcomment(pxy, "Poles");
}
/* expansion factor for axes */
tol = 0.15;
ymax = 0.0;
ymin = 0.0;
xmax = 0.0;
xmin = 0.0;
/* range */
if (length(z) > 0)
{
/* get Y coordinates */
(ymin, ymax) = splane_coords(yvals(zxy), ymin, ymax, tol);
/* get X coordinates */
(xmin, xmax) = splane_coords(xvals(zxy), xmin, xmax, tol);
}
/* handle poles */
if (length(p) > 0)
{
/* get Y coordinates */
(ymin, ymax) = splane_coords(yvals(pxy), ymin, ymax, tol);
/* get X coordinates */
(xmin, xmax) = splane_coords(xvals(pxy), xmin, xmax, tol);
}
xfac = (xmax - xmin) / 10;
yfac = (ymax - ymin) / 10;
/* expand */
xmin -= xfac;
xmax += xfac;
ymin -= yfac;
ymax += yfac;
/* draw circle as primary series */
t = 0..(2 * pi / 100)..(2 * pi);
x = cos(t) * (xmax - xmin);
y = sin(t) * (ymax - ymin);
xy(x, y);
setcomment("S Plane");
/*
* transparent color for unit circle data, draw ellipse for the
* actual displayed circle
*/
setcolor(-1);
/*
* solid line style for axes
*/
lstyle = 1;
setline(lstyle);
/* delete old lines */
h = findshape(w0, "tag", "splane_shape");
if (length(h) > 0)
{
deletehandle(h);
}
/* labels */
setxlabel("Real Part");
setylabel("Imaginary Part");
label("Pole-Zero Plot");
/* zeros and poles will be overplots - start at item 2 */
item = 2;
/* handle zeros */
if (length(z) > 0)
{
/* overplot zeros in black */
overp(zxy, sys_text_color);
/* point plot for zeros and poles */
setplotstyle(1, item);
/* set symbols - circle for zeros */
setsymbol(15, item);
item++;
}
/* handle poles */
if (length(p) > 0)
{
/* overplot poles in lred */
overp(pxy, lred);
/* point plot for poles */
setplotstyle(1, item);
/* set symbols - cross for poles */
setsymbol(5, item);
}
/* set axes */
setx(xmin, xmax);
sety(ymin, ymax);
/* set autoscale coordinates */
setxauto(getxl(), getxr());
setyauto(getyb(), getyt());
/* remove old text */
h = findtext(w0, "tag", "splane_text");
if (length(h) > 0)
{
deletehandle(h);
}
/* anotate multiple poles */
if (length(p) > 0)
{
(p, m) = splane_mpoles(p);
if ((len = length(p)) > 0)
{
loop (i = 1..len)
{
x = real(p[i]);
y = imag(p[i]);
splane_text(x, y, lred, sprintf("%d", (m[i])));
}
}
}
/* anotate multiple zeros */
if (length(z) > 0)
{
(z, m) = splane_mpoles(z);
if ((len = length(z)) > 0)
{
loop (i = 1..len)
{
x = real(z[i]);
y = imag(z[i]);
splane_text(x, y, black, sprintf("%d", (m[i])));
}
}
}
/*
* draw crosshair through the origin
*/
/* x line */
splane_line(lstyle, -inf, 0.0, inf, 0.0);
/* y line */
splane_line(lstyle, 0.0, -inf, 0.0, inf);
/* add axes labels */
scales(2);
/* now plot it */
plotmode(splane_pmode);
splane_pmode = -1;
}
/* determine axis coords */
splane_coords(vals, valmin, valmax, tol)
{
local vmax, vmin, v1, v2;
/* see if max|min * (1 +- tol) > valmax|min */
vmax = max(vals);
vmin = min(vals);
v1 = vmax + abs(vmax * tol);
v2 = vmin - abs(vmin * tol);
if (v1 > valmax) valmax = v1;
if (v2 < valmin) valmin = v2;
return(valmin, valmax);
}
/* returns poles and multiplicity */
splane_mpoles(p, tol)
{
local i, j, pp, pm, mcnt, idx;
if (argc < 2) tol = 1e-3;
tol2 = 1e-3;
/* sort in descending order */
p = sort(p, 0);
/* search for repeated poles and record multiplicity */
if (length(p) > 1)
{
idx = find(mag(real(p)) < tol2 && mag(imag(p)) < tol2);
p[idx] = real(p[idx]);
/* initial poles and multiplicity */
pm = ones(length(p), 1);
pp = ones(length(p), 1);
j = 1;
mcnt = 0;
loop (i = 2..length(p))
{
if (mag(p[i] - p[i-1]) <= (tol * mag(p[i-1])))
{
/* repeated root - accumulate multiplicity */
if (pm[j] == 1)
{
pp[j] = p[i];
mcnt++;
}
pm[j]++;
}
else if (pm[j] > 1)
{
j++;
}
}
if (mcnt > 0)
{
pm = extract(pm, 1, mcnt);
pp = extract(pp, 1, mcnt);
}
else
{
pm = {};
pp = {};
}
}
else
{
pm = {};
pp = {};
}
return(pp, pm);
}
/* text */
splane_text(x, y, color, str)
{
local h;
h = text(castreal(x), castreal(y), PAPER, color, -1, 1, str);
/* mark text */
h.tag = "splane_text";
/* don't show tag */
h.showtag = 0;
/* no border */
h.box = 0;
/* lock it down */
h.lock = 1;
}
/* line */
splane_line(lstyle, x1, y1, x2, y2)
{
local h;
/* draw line */
h = linedraw(getwcolor(-1), lstyle, PAPER, 1, 1, 1, castreal(x1), castreal(y1), castreal(x2), castreal(y2));
/* tag shape */
h.tag = "splane_shape";
/* no tag */
h.showtag = 0;
/* lock it down */
h.lock = 1;
}
/* error handler */
splane_error(errno, errmes)
{
if (splane_pmode >= 0)
{
/* make sure plotting is restored */
plotmode(splane_pmode);
splane_pmode = -1;
}
error(errmes);
}