View Raw SPL
/*****************************************************************************
* *
* ZP2CAS.SPL Copyright (C) 2015 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Zeros, poles and gain to cascade form *
* *
* Revisions: 15 Jun 2015 RRR Creation *
* *
*****************************************************************************/
#if @HELP_ZP2CAS
ZP2CAS
Purpose: Converts zeros, poles and gain to 2nd order cascade form.
Syntax: ZP2CAS(z, p, k, "flag")
z - A series, the zeros.
p - A series, the poles.
k - Optional. A scalar, the gain. Defaults to 1.0.
"flag" - Optional. A string, the stage ordering flag.
"up" : order stages such that the poles of each
stage are closer to the unit circle than
the previous stage (default)
"down" : order stages such that the poles of each
stage are closer to zero than the previous
stage
Returns: A series, the system coefficients in cascaded biquad form.
Example:
z = {0.0, 2.0};
p = {0.5, 0.2};
k = 1.0;
c = zp2cas(z, p, k);
c == {1, 1, -2, 0, -0.7, 0.1};
The cascade coefficients represent the system:
1 - 2z^(-1)
H(z) = (1) * --------------------------
1 - 0.7 z^(-1) + 0.1 z^(-2)
Remarks:
ZP2CAS converts the zeros, poles and gain of a discrete
system to cascade coefficient form where the coefficients
represent the following Z transform:
-1 -2 -1 -2
b10 + b11 z + b12 z b20 + b21 z + b22 z
H(z) = G * ________________________ * ________________________ ...
-1 -2 -1 -2
1 + a11 z + a12 z 1 + a21 z + a22 z
The cascade filter coefficients are represented by a single
column series in the following order:
G, b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...
To yield real coefficients, the zeros and poles must occur in
complex conjugate pairs or be real.
By default, the stages are ordered such that the poles of
each stage are closer to the unit circle than the previous
stage. The zeros of each stage are chosen to be closest to
the poles of the same stage.
If FLAG == "down", the stages are reversed.
See CAS2ZP to convert cascade coefficients to zeros, poles and
gain terms.
See Also:
Cascade
Cas2zp
Tf2zpk
Zp2tf
#endif
/* zeros, poles and gain to 2nd order cascade form */
zp2cas(z, p, k, flag)
{
local b, a, zp, pp, cas, nc, npr, nzp, ncz, ncp, dx = 1.0;
/* input args */
(z, p, k, flag) = zp2cas_parse_args(z, p, k, flag);
/* sample interval */
if (length(z) > 0)
{
dx = deltax(z);
}
else if (length(p) > 0)
{
dx = deltax(p);
}
/* zeros complex conjugate pairs */
zp = cplxpair(z);
/* poles complex conjugate pairs */
pp = cplxpair(p);
/* sort */
(z, p, nzr, npr) = zp2cas_sort(zp, pp);
ncz = numcols(z);
ncp = numcols(p);
if (ncz > ncp)
{
if (ncp == 0)
{
/* fir case */
cas = zp2cas_fir(zp, k);
}
else
{
error("zp2cas - number of zeros must be less than number of poles");
}
}
else
{
if (ncz < ncp)
{
/* more poles than zeros */
z = ravel(z, ones(1, ncp - ncz) * inf);
}
/* num columns for gain */
nc = max(ncz, ncp);
/* second order stage coefficients - unity gain */
(b, a) = zp2tf(z, p, ones(1, nc));
if (flag == "up")
{
/* order stages so succeeding stage poles closer to unit circle */
b = fliplr(b);
a = fliplr(a);
}
/* sizes */
b = extract(b, 1, 3);
a = extract(a, 2, 2);
/* combine stages */
cas = concat(b, a);
/* gain first term, unravel to single column series */
cas = {k, unravel(cas)};
}
/* comment and rate */
setcomment(cas, "Cascade");
setdeltax(cas, dx);
return(cas);
}
/* sort poles to unit circle, order zeros closest to poles */
zp2cas_sort(z, p)
{
local pc, pr, zc, zt, npr, npz, zzc, zzr, ppc, ppr, j, m, idx, n;
/* break into conjugate pairs and non-conjugate values */
(pc, pr) = zp2cas_conjreal(p, 1);
(zc, zr) = zp2cas_conjreal(z, 0);
npr = length(pr);
/* each column is a pole stage */
p = ravel(ravel((pc), 2), ravel((pr), 2));
zzc = {};
n = int(length(zc) / 2);
/* order conjugate pair zeros closest to poles */
loop (j = 1..n)
{
if (length(pc) > 0)
{
(m, idx) = min(mag(zc - pc[1]));
zzc = {zzc, zc[idx]};
zc[idx] = {};
(m, idx) = min(mag(zc - pc[2]));
zzc = {zzc, zc[idx]};
zc[idx] = {};
pc[{1, 2}] = {};
}
else if (length(pr) > 0)
{
(m, idx) = min(mag(zc - pr[1]));
zzc = {zzc, zc[idx], zc[idx + 1]};
zc[{idx, idx + 1}] = {};
pr[1] = {};
if (length(pr) > 0)
{
pr[1] = {};
}
}
else
{
zzc = {zzc, zc};
}
}
zzr = {};
/* order real zeros closest to poles */
loop (j = 1..length(zr))
{
if (length(pc) > 0)
{
(m, idx) = min(mag(zr - pc[1]));
zzr = {zzr, zr[idx]};
zr[idx] = {};
pc[1] = {};
}
else if (length(pr) > 0)
{
(m, idx) = min(mag(zr - pr[1]));
zzr = {zzr, zr[idx]};
zr[idx] = {};
pr[1] = {};
}
else
{
zzr = {zzr, zr};
}
}
npz = length(zzr);
/* each column is a zero stage */
z = ravel(ravel((zzc), 2), ravel((zzr), 2));
return(z, p, npz, npr);
}
/* break into complex conjugate pairs and non-conjugates */
zp2cas_conjreal(pz, sortf)
{
local idx, zc, zr, idx, ridx;
/* find real poles */
idx = find(mag(imag(pz)) > 0);
/* conjugate pairs */
zc = pz[idx];
if (sortf)
{
/* sort based on distance from unit circle */
ridx = grade(mag(zc - exp(i*phase(zc))), 1);
zc = zc[ridx];
}
/* non-conjugates */
ridx = 1..length(pz);
ridx[idx] = {};
zr = pz[ridx];
if (sortf)
{
/* sort based on magnitude, positive real part first */
ridx = grade(mag(zr - sign(zr)), 1);
zr = zr[ridx];
}
return(zc, zr);
}
/* handles all zeros case (fir) */
zp2cas_fir(broots, b0)
{
local n, k, m, i, brow, a1, b1, c;
/* sizes */
n = length(broots);
k = floor(n/2);
m = 2 * k;
b1 = zeros(k, 3);
a1 = zeros(k, 3);
/* find corresponding polynomials for poles and zeros */
a1[.., 1] = 1;
loop (i = 1..2..m)
{
brow = broots[i..1..i+1, ..];
brow = poly(brow);
b1[fix((i+1)/2), ..] = transpose(brow);
}
c = concat(transpose(b1), transpose(a1[.., 2..numcols(a1)]));
c = {b0, unravel(c)};
return(c);
}
zp2cas_parse_args(z, p, k, flag)
{
local ff, kk;
/* defaults */
ff = "up";
kk = 1.0;
if (not(isunspecified(flag)))
{
ff = tolower(caststring(flag));
}
if (not(isunspecified(k)))
{
if (isstring(k))
{
ff = k;
kk = 1.0;
}
else if (isarray(k))
{
kk = k[1];
}
else
{
kk = k;
}
}
if (isunspecified(p))
{
p = {};
}
if (isunspecified(z))
{
z = {};
}
/* checks */
if (not(ff == "up" || ff == "down"))
{
error('zp2cas - direction flag must be "up" or "down"');
}
if (not(isarray(z) && isarray(p)))
{
error("zp2cas - series required for poles and zeros");
}
return(z, p, kk, ff);
}