##
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 |

`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 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

`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