Passar função de transferência para forma Paralela

Eis aqui um script de matlab para passar uma função de transferência H(z) = b(z)/a(z) para a forma paralela, utilizando blocos de segunda ordem . A função residue retorna simplesmente valores de primeira ordem, com coeficientes complexos. Utilizando esta função, o valor retornado já é da forma (bk1z+bk2)/(z^2+ak1z+ak2).

1. function [num_paralelo,den_paralelo]=formaParalela(b,a)
3. % forma paralela dado um polinomio H=b/a
4. % input: formaParalela(b,a)
5. % resposta do tipo
6. %
7. % [0 bk1 bk2 1 ak1 ak2]
8. %
9. % onde a funcao correspondente e':
10. %
11. % bk1z+bk2
12. %-----------
13. %z^2+ak1z+ak2
14. %
15. % fonte: Diniz, Digital Signal Processing 2nd Ed (p. 265)
17. [Residuos,Polos,Constante]=residue(b,a);
18. Nparalelo=length(Polos);
19. polos_reais=cplxpair(Polos);
20. residuos_reais=zeros(Nparalelo,1);
21. for i=1:Nparalelo,
22.    residuos_reais(find(polos_reais==Polos(i)),1)=Residuos(i);
23. end;
24. num_blocks=ceil(Nparalelo/2);
25. Hparalelo=zeros(num_blocks,6);
26. for count_p =1:num_blocks,
27.   if length(polos_reais)~=1,
28.     Hparalelo(count_p,2)=residuos_reais(1)+residuos_reais(2);
29.     Hparalelo(count_p,3)=-residuos_reais(1)*polos_reais(2)-polos_reais(1)*residuos_reais(2);
30.     Hparalelo(count_p,5)=-polos_reais(1)-polos_reais(2);
31.     Hparalelo(count_p,6)=polos_reais(1)*polos_reais(2);
32.     residuos_reais(1:2) = []; polos_reais(1:2) = [];
33.   else,
34.     Hparalelo(count_p,2) = residuos_reais(1);
35.     Hparalelo(count_p,5) = -polos_reais(1);
36.   end;
37. Hparalelo(count_p,4) = 1;
38. end;
39. Hparalelo= real(Hparalelo)

O arquivo deve ser salvo com o nome formaParalela.m e, caso precise de ajuda, o cabeçalho comentado é acessado através de help formaParalela.