##
Partial Fraction Expansion: `residuez.m`

Figure J.8 gives a listing of a matlab function for computing
a ``left-justified'' partial fraction expansion (PFE) of an IIR
digital filter
as described in §6.8 (and
below). This function, along with its ``right justified''
counterpart, `residued`, are included in the
`octave-forge` matlab library for Octave.^{J.1}

function [r, p, f, m] = residuez(B, A, tol) if nargin<3, tol=0.001; end NUM = B(:)'; DEN = A(:)'; % Matlab's residue does not return m (implied by p): [r,p,f,m]=residue(conj(fliplr(NUM)),conj(fliplr(DEN)),tol); p = 1 ./ p; r = r .* ((-p) .^m); if f, f = conj(fliplr(f)); end |

This code was written for Octave, but it also runs in Matlab if the
'`m`' outputs (pole multiplicity counts) are omitted (two places).
The input arguments are compatible with the existing
`residuez` function in the Matlab Signal Processing Toolbox.

### Method

As can be seen from the code listing, this implementation of
`residuez` simply calls `residue`, which was written to
carry out the partial fraction expansions of -plane
(continuous-time) transfer functions :

where is the ``quotient'' and is the ``remainder'' in the PFE:

where is the order of the quotient polynomial in , and is the

*multiplicity*of the th pole. (When all poles are distinct, we have for all .) For , we define .

In the discrete-time case, we have the -plane transfer function

For compatibility with Matlab's

`residuez`, we need a PFE of the form such that

where .

We see that the -plane case formally does what we desire if we
treat -plane polynomials as polynomials in instead of
. From Eq.(J.2), we see that this requires reversing the
coefficient-order of `B` and `A` in the call to
`residue`. In the returned result, we obtain terms such as

### Example with Repeated Poles

The following Matlab code performs a partial fraction expansion of a
filter having three pairs of repeated roots (one real and two
complex):^{J.2}

N = 1000; % number of time samples to compute A = [ 1 0 0 1 0 0 0.25]; B = [ 1 0 0 0 0 0 0]; % Compute "trusted" impulse response: h_tdl = filter(B, A, [1 zeros(1, N-1)]); % Compute partial fraction expansion (PFE): [R,P,K] = residuez(B, A); % PFE impulse response: n = [0:N-1]; h_pfe = zeros(1,N); for i = 1:size(R) % repeated roots are not detected exactly: if i>1 && abs(P(i)-P(i-1))<1E-7 h_pfe = h_pfe + (R(i) * (n+1) .* P(i).^n); disp(sprintf('Pole %d is a repeat of pole %d',i,i-1)); % if i>2 && abs(P(i)-P(i-2))<1E-7 ... else h_pfe = h_pfe + (R(i) * P(i).^n); end end err = abs(max(h_pfe-h_tdl)) % should be about 5E-8

**Next Section:**

Partial Fraction Expansion: residued.m

**Previous Section:**

Frequency Response Plots: plotfr.m