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):^H.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:6
  % 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