|
2847
|
1 ## Copyright (C) 1996, 1997 John W. Eaton |
|
6964
|
2 ## Copyright (C) 2007 Ben Abbott |
|
2313
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
9 ## any later version. |
|
|
10 ## |
|
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
|
17 ## along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
18 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
19 ## 02110-1301, USA. |
|
904
|
20 |
|
3368
|
21 ## -*- texinfo -*- |
|
6964
|
22 ## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}] =} residue (@var{b}, @var{a}) |
|
|
23 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) |
|
|
24 ## In the instance of two inputs, they are assumed to correspond to vectors |
|
|
25 ## of polynomial coefficients, @var{b} and @var{a}. From these polynomial |
|
|
26 ## coefficients, the function residue calculates the partial fraction |
|
|
27 ## expansion corresponding to the ratio of the two polynomials. |
|
3368
|
28 ## @cindex partial fraction expansion |
|
3426
|
29 ## |
|
6964
|
30 ## In this instance, the function @code{residue} returns @var{r}, @var{p}, |
|
|
31 ## and @var{k}, where the vector @var{r} contains the residue terms, |
|
|
32 ## @var{p} contains the pole values, @var{k} contains the coefficients of a |
|
|
33 ## direct polynomial term (if it exists). |
|
|
34 ## |
|
|
35 ## In the instance of three inputs, the function 'residue' performs the |
|
|
36 ## reciprocal task. Meaning the partial fraction expansion is reconstituted |
|
|
37 ## into the corresponding pair of polynomial coefficients. |
|
3426
|
38 ## |
|
3368
|
39 ## Assuming @var{b} and @var{a} represent polynomials |
|
|
40 ## @iftex |
|
|
41 ## @tex |
|
|
42 ## $P(s)$ and $Q(s)$ |
|
|
43 ## @end tex |
|
|
44 ## @end iftex |
|
|
45 ## @ifinfo |
|
|
46 ## P (s) and Q(s) |
|
|
47 ## @end ifinfo |
|
|
48 ## we have: |
|
|
49 ## @iftex |
|
|
50 ## @tex |
|
|
51 ## $$ |
|
|
52 ## {P(s)\over Q(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
|
|
53 ## + \sum_{i=1}^N k_i s^{N-i}. |
|
|
54 ## $$ |
|
|
55 ## @end tex |
|
|
56 ## @end iftex |
|
|
57 ## @ifinfo |
|
3426
|
58 ## |
|
3368
|
59 ## @example |
|
2311
|
60 ## P(s) M r(m) N |
|
3368
|
61 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
|
|
62 ## Q(s) m=1 (s-p(m))^e(m) i=1 |
|
|
63 ## @end example |
|
|
64 ## @end ifinfo |
|
3426
|
65 ## |
|
3368
|
66 ## @noindent |
|
3499
|
67 ## where @math{M} is the number of poles (the length of the @var{r}, |
|
|
68 ## @var{p}, and @var{e} vectors) and @math{N} is the length of the |
|
6964
|
69 ## @var{k} vector. The @var{e} vector specifies the multiplicity of the |
|
|
70 ## mth residue's pole. |
|
3426
|
71 ## |
|
6964
|
72 ## For example, |
|
3426
|
73 ## |
|
3368
|
74 ## @example |
|
|
75 ## @group |
|
6964
|
76 ## b = [1, 1, 1]; |
|
|
77 ## a = [1, -5, 8, -4]; |
|
|
78 ## [r, p, k] = residue (b, a); |
|
3368
|
79 ## @result{} r = [-2, 7, 3] |
|
|
80 ## @result{} p = [2, 2, 1] |
|
|
81 ## @result{} k = [](0x0) |
|
|
82 ## @end group |
|
|
83 ## @end example |
|
3426
|
84 ## |
|
3368
|
85 ## @noindent |
|
|
86 ## which implies the following partial fraction expansion |
|
|
87 ## @iftex |
|
|
88 ## @tex |
|
|
89 ## $$ |
|
|
90 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} |
|
|
91 ## $$ |
|
|
92 ## @end tex |
|
|
93 ## @end iftex |
|
|
94 ## @ifinfo |
|
3426
|
95 ## |
|
3368
|
96 ## @example |
|
|
97 ## s^2 + s + 1 -2 7 3 |
|
2311
|
98 ## ------------------- = ----- + ------- + ----- |
|
|
99 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) |
|
3368
|
100 ## @end example |
|
6964
|
101 ## |
|
3368
|
102 ## @end ifinfo |
|
6964
|
103 ## A similar, but reciprocal example, where the fraction's polynomials are |
|
|
104 ## reconstituted from the residues, poles, and direct term is |
|
|
105 ## |
|
|
106 ## @example |
|
|
107 ## @group |
|
|
108 ## r = [-2, 7, 3]; |
|
|
109 ## p = [2, 2, 1]; |
|
|
110 ## k = [1 0]; |
|
|
111 ## [b, a] = residue (r, p, k); |
|
|
112 ## @result{} b = [1, -5, 9, -3, 1] |
|
|
113 ## @result{} a = [1, -5, 8, 4] |
|
|
114 ## @end group |
|
|
115 ## @end example |
|
|
116 ## |
|
|
117 ## @noindent |
|
|
118 ## which implies the following partial fraction expansion |
|
|
119 ## @iftex |
|
|
120 ## @tex |
|
|
121 ## $$ |
|
|
122 ## {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s |
|
|
123 ## $$ |
|
|
124 ## @end tex |
|
|
125 ## @end iftex |
|
|
126 ## @ifinfo |
|
|
127 ## |
|
|
128 ## @example |
|
|
129 ## s^4 - 5s^3 + 9s^2 - 3s + 1 -2 7 3 |
|
|
130 ## -------------------------- = ----- + ------- + ----- + s |
|
|
131 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) |
|
|
132 ## @end example |
|
|
133 ## @end ifinfo |
|
|
134 ## @seealso{poly, roots, conv, deconv, mpoles, polyval, polyderiv, polyinteg} |
|
3368
|
135 ## @end deftypefn |
|
1025
|
136 |
|
3202
|
137 ## Author: Tony Richardson <arichard@stark.cc.oh.us> |
|
6964
|
138 ## Author: Ben Abbott <bpabbott@mac.com> |
|
2312
|
139 ## Created: June 1994 |
|
|
140 ## Adapted-By: jwe |
|
559
|
141 |
|
6964
|
142 function [r, p, k] = residue (b, a, varargin) |
|
559
|
143 |
|
1025
|
144 if (nargin < 2 || nargin > 3) |
|
6046
|
145 print_usage (); |
|
559
|
146 endif |
|
|
147 |
|
6964
|
148 toler = .001; |
|
|
149 |
|
|
150 if (nargin == 3) |
|
|
151 ## The inputs are the residue, pole, and direct part. Solve for the |
|
|
152 ## corresponding numerator and denominator polynomials |
|
|
153 [r, p] = rresidue (b, a, varargin{1}, toler); |
|
|
154 return |
|
|
155 end |
|
559
|
156 |
|
2303
|
157 ## Make sure both polynomials are in reduced form. |
|
1025
|
158 |
|
|
159 a = polyreduce (a); |
|
|
160 b = polyreduce (b); |
|
559
|
161 |
|
1025
|
162 b = b / a(1); |
|
|
163 a = a / a(1); |
|
559
|
164 |
|
1025
|
165 la = length (a); |
|
|
166 lb = length (b); |
|
559
|
167 |
|
2303
|
168 ## Handle special cases here. |
|
1025
|
169 |
|
|
170 if (la == 0 || lb == 0) |
|
559
|
171 k = r = p = e = []; |
|
|
172 return; |
|
|
173 elseif (la == 1) |
|
1025
|
174 k = b / a; |
|
559
|
175 r = p = e = []; |
|
|
176 return; |
|
|
177 endif |
|
|
178 |
|
2303
|
179 ## Find the poles. |
|
1025
|
180 |
|
|
181 p = roots (a); |
|
|
182 lp = length (p); |
|
559
|
183 |
|
5464
|
184 ## Determine if the poles are (effectively) zero. |
|
5463
|
185 |
|
6964
|
186 small = max (abs (p)); |
|
|
187 small = max ([small, 1] ) * 1e-8 * (1 + numel (p))^2; |
|
|
188 p(abs (p) < small) = 0; |
|
|
189 |
|
|
190 ## Determine if the poles are (effectively) real, or imaginary. |
|
5463
|
191 |
|
6964
|
192 index = (abs (imag (p)) < small); |
|
|
193 p(index) = real (p(index)); |
|
|
194 index = (abs (real (p)) < small); |
|
|
195 p(index) = 1i * imag (p(index)); |
|
1025
|
196 |
|
6964
|
197 ## Sort poles so that multiplicity loop will work. |
|
|
198 |
|
|
199 [e, indx] = mpoles (p, toler, 1); |
|
|
200 p = p (indx); |
|
559
|
201 |
|
2303
|
202 ## Find the direct term if there is one. |
|
1025
|
203 |
|
|
204 if (lb >= la) |
|
6964
|
205 ## Also return the reduced numerator. |
|
1025
|
206 [k, b] = deconv (b, a); |
|
|
207 lb = length (b); |
|
559
|
208 else |
|
|
209 k = []; |
|
|
210 endif |
|
|
211 |
|
1025
|
212 if (lp == 1) |
|
|
213 r = polyval (b, p); |
|
559
|
214 return; |
|
|
215 endif |
|
|
216 |
|
6964
|
217 ## Determine the order of the denominator and remaining numerator. |
|
|
218 ## With the direct term removed the potential order of the numerator |
|
|
219 ## is one less than the order of the denominator. |
|
1025
|
220 |
|
6964
|
221 aorder = numel (a) - 1; |
|
|
222 border = aorder - 1; |
|
1025
|
223 |
|
6964
|
224 ## Construct a system of equations relating the individual |
|
|
225 ## contributions from each residue to the complete numerator. |
|
559
|
226 |
|
6964
|
227 A = zeros (border+1, border+1); |
|
|
228 B = prepad (reshape (b, [numel(b), 1]), border+1, 0); |
|
|
229 for ip = 1:numel(p) |
|
|
230 ri = zeros (size (p)); |
|
|
231 ri(ip) = 1; |
|
|
232 A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).'; |
|
|
233 endfor |
|
559
|
234 |
|
2303
|
235 ## Solve for the residues. |
|
1025
|
236 |
|
|
237 r = A \ B; |
|
559
|
238 |
|
|
239 endfunction |
|
6964
|
240 |
|
|
241 function [pnum, pden] = rresidue (r, p, k, toler) |
|
|
242 |
|
|
243 ## Reconstitute the numerator and denominator polynomials from the |
|
|
244 ## residues, poles, and direct term. |
|
|
245 |
|
|
246 if (nargin < 2 || nargin > 4) |
|
|
247 print_usage (); |
|
|
248 endif |
|
|
249 |
|
|
250 if (nargin < 4) |
|
|
251 toler = []; |
|
|
252 endif |
|
|
253 |
|
|
254 if (nargin < 3) |
|
|
255 k = []; |
|
|
256 endif |
|
|
257 |
|
|
258 [multp, indx] = mpoles (p, toler, 0); |
|
|
259 |
|
|
260 p = p (indx); |
|
|
261 r = r (indx); |
|
|
262 |
|
|
263 indx = 1:numel(p); |
|
|
264 |
|
|
265 for n = indx |
|
|
266 pn = [1, -p(n)]; |
|
|
267 if n == 1 |
|
|
268 pden = pn; |
|
|
269 else |
|
|
270 pden = conv (pden, pn); |
|
|
271 endif |
|
|
272 endfor |
|
|
273 |
|
|
274 ## D is the order of the denominator |
|
|
275 ## K is the order of the direct polynomial |
|
|
276 ## N is the order of the resulting numerator |
|
|
277 ## pnum(1:(N+1)) is the numerator's polynomial |
|
|
278 ## pden(1:(D+1)) is the denominator's polynomial |
|
|
279 ## pm is the multible pole for the nth residue |
|
|
280 ## pn is the numerator contribution for the nth residue |
|
|
281 |
|
|
282 D = numel (pden) - 1; |
|
|
283 K = numel (k) - 1; |
|
|
284 N = K + D; |
|
|
285 pnum = zeros (1, N+1); |
|
|
286 for n = indx(abs(r)>0) |
|
|
287 p1 = [1, -p(n)]; |
|
|
288 for m = 1:multp(n) |
|
|
289 if m == 1 |
|
|
290 pm = p1; |
|
|
291 else |
|
|
292 pm = conv (pm, p1); |
|
|
293 endif |
|
|
294 endfor |
|
|
295 pn = deconv (pden, pm); |
|
|
296 pn = r(n) * pn; |
|
|
297 pnum = pnum + prepad ( pn, N+1, 0); |
|
|
298 endfor |
|
|
299 |
|
|
300 ## Add the direct term. |
|
|
301 |
|
|
302 if (numel (k)) |
|
|
303 pnum = pnum + conv (pden, k); |
|
|
304 endif |
|
|
305 |
|
|
306 ## Check for leading zeros and trim the polynomial coefficients. |
|
|
307 |
|
|
308 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; |
|
|
309 |
|
|
310 pnum (abs (pnum) < small) = 0; |
|
|
311 pden (abs (pden) < small) = 0; |
|
|
312 |
|
|
313 pnum = polyreduce (pnum); |
|
|
314 pden = polyreduce (pden); |
|
|
315 |
|
|
316 endfunction |
|
6968
|
317 |
|
|
318 %% test/octave.test/poly/residue-1.m |
|
|
319 %!test |
|
|
320 %! b = [1, 1, 1]; |
|
|
321 %! a = [1, -5, 8, -4]; |
|
|
322 %! [r, p, k, e] = residue (b, a); |
|
|
323 %! assert((abs (r - [-2; 7; 3]) < 1e-6 |
|
|
324 %! && abs (p - [2; 2; 1]) < 1e-7 |
|
|
325 %! && isempty (k) |
|
|
326 %! && e == [1; 2; 1])); |
