View Raw SPL
/*****************************************************************************
* *
* RESIDUEZ.SPL Copyright (C) 2003-2018 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Partial fraction expansion of a Z transform *
* *
* Revisions: 1 Oct 2003 RRR Creation *
* 1 Nov 2018 RRR complex fix *
* *
*****************************************************************************/
#include
#if @HELP_RESIDUEZ
RESIDUEZ
Purpose: Finds the partial fraction expansion of a Z transform.
Syntax: RESIDUEZ(b, a)
(r, p, k) = RESIDUEZ(b, a)
b - A series, the numerator coefficients in ascending
powers of z^(-1).
a - A series, the denominator coefficients in ascending
powers of z^(-1).
Alternative Syntax:
RESIDUEZ(r, p, k)
(b, a, k) = RESIDUEZ(r, p, k)
r - A series, the residues representing the numerator terms
of the partial fraction expansion.
p - A series, the poles of the partial fraction expansion.
k - A series, the numerator coefficients for the direct
terms of the partial fraction expansion.
Returns: (r, p, k) = RESIDUEZ(b, a) returns the partial fraction expansion
of the polynomial.
R, p and k are series where r represents the residues of
the partial fraction expansion, p are the pole locations and k
represents the direct terms (if any).
RESIDUEZ(b, a) with one or zero output arguments returns r, p
and k in one series of three columns where r == COL(1),
p == COL(2) and k == COL(3).
(b, a, c) = RESIDUEZ(r, p, k) returns inverse partial fraction
expansion, converting the partial fraction expansion back into
b(z)/a(z) form. Series b and a are the numerator and denominator
terms of the rational polynomial with the partial fraction
expansion represented by r, p, and k. Series c represents the
remainder terms (if any).
RESIDUEZ(r, p, k) with one or zero output arguments returns b
and a in one series of two columns where b == COL(1) and
a == COL(2).
RESIDUEZ(f) or (b, a) = RESIDUEZ(f) assumes f is a three
column series with r, p and k as each of the columns.
Thus RESIDUEZ(RESIDUEZ(b, a)) == {{b/a[1], a/a[1]}}.
Example:
1 - z^-1 1 - z^-1
Given H(z) = ------------------- = ----------------------
1 - 5z^-1 + 6z^- 2 (1 - 3z^-1)(1 - 2z^-1)
(r, p, k) = residuez({1, -1}, {1, -5, 6})
r == {2, -1}
p == {3, 2}
k == {}
Representing the partial fraction expansion:
2 -1
H(z) = ----------- + -----------
(1 - 3z^-1) (1 - 2z^-1)
From this partial fraction expansion, the impulse response
can be derived by inspection: h[n] = 2 * 3^n - 2^n. Because the
poles are outside the unit circle, the system is unstable.
Now, performing the inverse transform:
(b, a, c) = residuez(r, p, k)
b == {1, -1}
a == {1, -5, 6}
c == {}
The series b and a represent the numerator and denominator
terms of the original rational polynomial.
Example:
2 + 3z^-1 + 4z^-2
Given H(z) = ------------------------
1 + 3z^-1 + 3z^-2 + z^-3
(r, p, k) = residuez({2, 3, 4}, {1, 3, 3, 1})
r == {4, -5, 3}
p == {-1, -1, -1}
k == {}
Since H(z) contains 3 repeated poles, the resulting partial
fraction expansion becomes:
4 -5 3
H(z) = -------- + ------------ + ------------
1 + z^-1 (1 + z^-1)^2 (1 + z^-1)^3
Now, performing the inverse transform:
(b, a, c) = residuez(r, p, k)
b == {2, 3, 4}
a == {1, 3, 3, 1}
c == {}
The series b and a represent the numerator and denominator
terms of the original rational polynomial.
Example:
z^3 - 10z^2 - 4z + 4
H(z) = --------------------
2z^2 - 2z - 4
1 - 10z^-1 - 4z^-2 + 4z^-3
or H(z) = --------------------------
0 + 2z^-1 - 2z^-2 - 4z^-3
Since residuez expects H(z) to be in terms of z^-1 and
the first denominator term cannot be zero, we form:
G(z) = H(z)/z
1 - 10z^-1 - 4z^-2 + 4z^-3
G(z) = H(z)/z = --------------------------
2 - 2z^-1 - 4z^-2
Now we find the partial fraction expansion for G(z):
(r, p, k) = residuez({1, -10, -4, 4}, {2, -2, -4})
r == {-1.5, 0.5}
p == {2, -1}
k == {1.5, -1}
Representing the partial fraction expansion:
-1.5 0.5
G(z) = --------- + -------- + 1.5 - z^-1
1 - 2z^-1 1 + z^-1
Since H(z) = zG(z), we have:
-1.5z^2 0.5z^2
H(z) = ------- + -------- + 1.5z - 1
z - 2 z + 1
Remarks:
Given the rational polynomial H(z) = b(z)/a(z) where:
M = length(b) and N = length(a):
b(z) b[1] + b[2] z^(-1) + ... + b[M] z^(-M+1)
H(z) = ------ = ----------------------------------------
a(z) a[1] + a[2] z^(-1) + ... + a[N] z^(-N+1)
If a[1] does not equal 1, the numerator and denominator terms
are normalized by dividing each coefficient by a[1].
If there are no repeated roots, the partial fraction expansion
of the rational polynomial is of the form:
r[1] r[N-1]
H(z) = ------------ + ... + -------------- + k[1] + k[2]z^(-1) ...
1-p[1]z^(-1) 1-p[N-1]z^(-1)
If there are K repeated roots (closer than 1.0e-3), then
the partial fraction expansion includes terms such as:
r[i] r[i+1] r[i+K-1]
-------------- + ---------------- + ... + ----------------
(1-p[i]z^(-1)) (1-p[i]z^(-1))^2 (1-p[i]z^(-1))^K
See Also:
Impz
Poly
Residue
Roots
Zfreq
Zplane
#endif
/* partial fraction expansion of H(z) */
residuez(b, a, k)
{
local lb, la, lm, acnt, r, p, c, rpole, pm, pp, idx;
/* parse args */
if (argc == 1)
{
if (numcols(b) > 1)
{
if (numcols(b) > 2)
{
k = col(b, 3);
}
else
{
k = {};
}
a = col(b, 2);
}
else
{
b = {b};
a = {1};
k = {};
}
b = col(b, 1);
acnt = 3;
}
else
{
b = {b};
a = {a};
acnt = argc;
}
lb = length(b);
la = length(a);
lm = max(lb, la);
if (acnt == 3)
{
/*
* forming original expansion from partial fraction
*/
/* check repeated poles or leading zero coefs in denominator */
if (reproot(a, 0) || (lb >= la) || iscomplex(a))
{
/*
* convert H(s) residues & poles to H(z) residues & poles
* Rz = Rs / (-Ps)^m where m is multiplicity
* Pz = 1 / Ps
*/
rpole = 1;
pm = residuez_mpole(a);
/* residues */
b = b / ((-pm[.., 1]) ^ pm[.., 2]);
/* poles */
a = 1 / a;
k = rev(k);
}
else
{
rpole = 0;
la = length(a);
if (length(k) > 0)
{
b = {b, k};
k = {};
}
a = extract(a, 1, length(b));
}
/* use normal residue to find polynomial coefficients */
(r, p, c) = residue(b, a, k);
if (rpole)
{
r = rev(r);
p = rev(p);
c = {};
r /= p[1];
p /= p[1];
}
else
{
p = extract(p, 1, la + 1);
}
/* dress it up */
setcomment(r, "Num");
setcomment(p, "Den");
setcomment(c, "Rem");
}
else
{
/* find partial fraction expansion */
if (reproot(a, 1) || (lb > la) || iscomplex(a))
{
/* have repeated poles or leading zero coeffs in denominator */
rpole = 1;
/* find H(s) for s = z^-1 */
b = reverse({b});
a = reverse({a});
}
else
{
rpole = 0;
/* form H(z)/z */
b = extract({b}, 1, lm);
a = extract({a}, 1, lm + 1);
}
/* now use standard residue function */
(r, p, c) = residue(b, a);
/* check last pole */
if (la <= length(p) && not(rpole))
{
idx = grade(p, "descend");
pp = p[idx];
if (pp[la] == conj(pp[la-1]))
{
p = pp;
r = r[idx];
p[la-1] *= 2;
r[la-1] *= 2;
}
}
if (rpole)
{
/*
* convert H(s) residues & poles to H(z) residues & poles
* Rz = Rs / (-Ps)^m where m is multiplicity
* Pz = 1 / Ps
*/
/* to sort the resulting pole order */
idx = grade(p, 1);
/* get poles and multiplicity */
pm = residuez_mpole(p);
/* residues */
r = r / ((-pm[.., 1]) ^ pm[.., 2]);
/* poles */
p = 1 / p;
/* sort */
r = r[idx];
p = p[idx];
c = rev(c);
}
else
{
/* reform to proper H(z) */
if (la > lb)
{
c = {};
}
else
{
c = extract(r, la, -1);
}
}
r = extract(r, 1, la - 1);
p = extract(p, 1, la - 1);
if (length(c) > 0)
{
/* see residue.spl */
c = residue_check_imag_tol(c);
}
/* dress it up */
setcomment(r, "Residues");
setcomment(p, "Poles");
setcomment(c, "Direct Terms");
}
/* check for imaginary parts near zero (residue.spl) */
r = residue_check_imag_tol(r);
p = residue_check_imag_tol(p);
if (outargc <= 1)
{
/* place result in current window */
r = ravel(r, p, c);
/* set to tableview */
setplotstyle(r, 4);
return(r);
}
else
{
return(r, p, c);
}
}
/* returns poles and multiplicity */
residuez_mpole(p, tol)
{
local i, pm;
if (argc < 2) tol = 1e-3;
/* initial poles and multiplicity */
pm = ravel(p, ones(length(p), 1));
/* search for repeated poles and accumulate multiplicity */
if (length(p) > 1)
{
loop (i = 2..length(p))
{
if (mag(pm[i, 1] - pm[i-1, 1]) <= tol)
{
/* repeated root - accumulate multiplicity */
pm[i, 2] += pm[i-1, 2];
}
}
}
return(pm);
}