1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
|
// Copyright (C) 2018 - IIT Bombay - FOSSEE
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// Author:[insert name]
// Organization: FOSSEE, IIT Bombay
// Email: toolbox@scilab.in
function [r, p, f, m] = residuez(B, A, tol)
// RESIDUEZ - return residues, poles, and FIR part of B(z)/A(z)
//
// Let nb = length(b), na = length(a), and N=na-1 = no. of poles.
// If nb<na, then f will be empty, and the returned filter is
//
// r(1) r(N)
// H(z) = ---------------- + ... + ----------------- = R(z)
// [ 1-p(1)/z ]^m(1) [ 1-p(N)/z ]^m(N)
//
// If, on the other hand, nb >= na, the FIR part f will not be empty.
// Let M = nb-na+1 = order of f = length(f)-1). Then the returned filter is
//
// H(z) = f(1) + f(2)/z + f(3)/z^2 + ... + f(M+1)/z^M + R(z)
//
// where R(z) is the parallel one-pole filter bank defined above.
// Note, in particular, that the impulse-response of the one-pole
// filter bank is in parallel with that of the the FIR part. This can
// be wasteful when matching the initial impulse response is important,
// since F(z) can already match the first N terms of the impulse
// response. To obtain a decomposition in which the impulse response of
// the IIR part R(z) starts after that of the FIR part F(z), use RESIDUED.
//
//NOTE that the polynomials 'b' and 'a' should have real coefficients(because of the function 'filter' used in polyval)
//Testcase
//B=[1 1 1]; A=[1 -2 1];
//[r,p,f,m] = residuez(B,A)
//OUTPUT:
//r=[0;3]
//p=[1;1]
//f=1
//e=[1;2]
[nargout,nargin]=argn();
if nargin==3
warning("tolerance ignored");
end
NUM = B(:)';
DEN = A(:)';
// Matlab's residue does not return m (since it is implied by p):
[r,p,f,m]=residue(conj(mtlb_fliplr(NUM)),conj(mtlb_fliplr(DEN)));
p = 1 ./ p;
r = r .* ((-p) .^m);
if f
f = conj(mtlb_fliplr(f));
end
endfunction
|
