|
7017
|
1 ## Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2004, 2005 |
|
|
2 ## 2006, 2007 John W. Eaton |
|
6964
|
3 ## Copyright (C) 2007 Ben Abbott |
|
2313
|
4 ## |
|
|
5 ## This file is part of Octave. |
|
|
6 ## |
|
|
7 ## Octave is free software; you can redistribute it and/or modify it |
|
|
8 ## under the terms of the GNU General Public License as published by |
|
7016
|
9 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
10 ## your option) any later version. |
|
2313
|
11 ## |
|
|
12 ## Octave is distributed in the hope that it will be useful, but |
|
|
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
15 ## General Public License for more details. |
|
|
16 ## |
|
|
17 ## You should have received a copy of the GNU General Public License |
|
7016
|
18 ## along with Octave; see the file COPYING. If not, see |
|
|
19 ## <http://www.gnu.org/licenses/>. |
|
904
|
20 |
|
3368
|
21 ## -*- texinfo -*- |
|
6978
|
22 ## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}, @var{e}] =} residue (@var{b}, @var{a}) |
|
|
23 ## Compute the partial fraction expansion for the quotient of the |
|
|
24 ## polynomials, @var{b} and @var{a}. |
|
3426
|
25 ## |
|
3368
|
26 ## @iftex |
|
|
27 ## @tex |
|
|
28 ## $$ |
|
6978
|
29 ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
|
3368
|
30 ## + \sum_{i=1}^N k_i s^{N-i}. |
|
|
31 ## $$ |
|
|
32 ## @end tex |
|
|
33 ## @end iftex |
|
|
34 ## @ifinfo |
|
3426
|
35 ## |
|
3368
|
36 ## @example |
|
6978
|
37 ## B(s) M r(m) N |
|
3368
|
38 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
|
6978
|
39 ## A(s) m=1 (s-p(m))^e(m) i=1 |
|
3368
|
40 ## @end example |
|
|
41 ## @end ifinfo |
|
3426
|
42 ## |
|
3368
|
43 ## @noindent |
|
3499
|
44 ## where @math{M} is the number of poles (the length of the @var{r}, |
|
6978
|
45 ## @var{p}, and @var{e}), the @var{k} vector is a polynomial of order @math{N-1} |
|
|
46 ## representing the direct contribution, and the @var{e} vector specifies |
|
|
47 ## the multiplicity of the mth residue's pole. |
|
3426
|
48 ## |
|
6964
|
49 ## For example, |
|
3426
|
50 ## |
|
3368
|
51 ## @example |
|
|
52 ## @group |
|
6964
|
53 ## b = [1, 1, 1]; |
|
|
54 ## a = [1, -5, 8, -4]; |
|
7011
|
55 ## [r, p, k, e] = residue (b, a); |
|
|
56 ## @result{} r = [-2; 7; 3] |
|
|
57 ## @result{} p = [2; 2; 1] |
|
3368
|
58 ## @result{} k = [](0x0) |
|
7011
|
59 ## @result{} e = [1; 2; 1] |
|
3368
|
60 ## @end group |
|
|
61 ## @end example |
|
3426
|
62 ## |
|
3368
|
63 ## @noindent |
|
6978
|
64 ## which represents the following partial fraction expansion |
|
3368
|
65 ## @iftex |
|
|
66 ## @tex |
|
|
67 ## $$ |
|
|
68 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} |
|
|
69 ## $$ |
|
|
70 ## @end tex |
|
|
71 ## @end iftex |
|
|
72 ## @ifinfo |
|
3426
|
73 ## |
|
3368
|
74 ## @example |
|
|
75 ## s^2 + s + 1 -2 7 3 |
|
2311
|
76 ## ------------------- = ----- + ------- + ----- |
|
|
77 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) |
|
3368
|
78 ## @end example |
|
6964
|
79 ## |
|
3368
|
80 ## @end ifinfo |
|
6978
|
81 ## |
|
|
82 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) |
|
7011
|
83 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e}) |
|
6978
|
84 ## Compute the reconstituted quotient of polynomials, |
|
7398
|
85 ## @var{b}(s)/@var{a}(s), from the partial fraction expansion; |
|
6978
|
86 ## represented by the residues, poles, and a direct polynomial specified |
|
7011
|
87 ## by @var{r}, @var{p} and @var{k}, and the pole multiplicity @var{e}. |
|
|
88 ## |
|
|
89 ## If the multiplicity, @var{e}, is not explicitly specified the multiplicity is |
|
|
90 ## determined by the script mpoles.m. |
|
6978
|
91 ## |
|
|
92 ## For example, |
|
6964
|
93 ## |
|
|
94 ## @example |
|
|
95 ## @group |
|
7011
|
96 ## r = [-2; 7; 3]; |
|
|
97 ## p = [2; 2; 1]; |
|
|
98 ## k = [1, 0]; |
|
6964
|
99 ## [b, a] = residue (r, p, k); |
|
|
100 ## @result{} b = [1, -5, 9, -3, 1] |
|
6978
|
101 ## @result{} a = [1, -5, 8, -4] |
|
7011
|
102 ## |
|
|
103 ## where mpoles.m is used to determine e = [1; 2; 1] |
|
|
104 ## |
|
|
105 ## @end group |
|
|
106 ## @end example |
|
|
107 ## |
|
|
108 ## Alternatively the multiplicity may be defined explicitly, for example, |
|
|
109 ## |
|
|
110 ## @example |
|
|
111 ## @group |
|
|
112 ## r = [7; 3; -2]; |
|
|
113 ## p = [2; 1; 2]; |
|
|
114 ## k = [1, 0]; |
|
|
115 ## e = [2; 1; 1]; |
|
|
116 ## [b, a] = residue (r, p, k, e); |
|
|
117 ## @result{} b = [1, -5, 9, -3, 1] |
|
|
118 ## @result{} a = [1, -5, 8, -4] |
|
6964
|
119 ## @end group |
|
|
120 ## @end example |
|
|
121 ## |
|
|
122 ## @noindent |
|
6978
|
123 ## which represents the following partial fraction expansion |
|
6964
|
124 ## @iftex |
|
|
125 ## @tex |
|
|
126 ## $$ |
|
6978
|
127 ## {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s = {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} |
|
6964
|
128 ## $$ |
|
|
129 ## @end tex |
|
|
130 ## @end iftex |
|
|
131 ## @ifinfo |
|
|
132 ## |
|
|
133 ## @example |
|
6978
|
134 ## -2 7 3 s^4 - 5s^3 + 9s^2 - 3s + 1 |
|
|
135 ## ----- + ------- + ----- + s = -------------------------- |
|
|
136 ## (s-2) (s-2)^2 (s-1) s^3 - 5s^2 + 8s - 4 |
|
6964
|
137 ## @end example |
|
|
138 ## @end ifinfo |
|
|
139 ## @seealso{poly, roots, conv, deconv, mpoles, polyval, polyderiv, polyinteg} |
|
3368
|
140 ## @end deftypefn |
|
1025
|
141 |
|
3202
|
142 ## Author: Tony Richardson <arichard@stark.cc.oh.us> |
|
6964
|
143 ## Author: Ben Abbott <bpabbott@mac.com> |
|
2312
|
144 ## Created: June 1994 |
|
|
145 ## Adapted-By: jwe |
|
559
|
146 |
|
6978
|
147 function [r, p, k, e] = residue (b, a, varargin) |
|
559
|
148 |
|
7011
|
149 if (nargin < 2 || nargin > 4) |
|
6046
|
150 print_usage (); |
|
559
|
151 endif |
|
|
152 |
|
6964
|
153 toler = .001; |
|
|
154 |
|
7011
|
155 if (nargin >= 3) |
|
|
156 if (nargin >= 4) |
|
|
157 e = varargin{2}; |
|
|
158 else |
|
|
159 e = []; |
|
|
160 endif |
|
6964
|
161 ## The inputs are the residue, pole, and direct part. Solve for the |
|
|
162 ## corresponding numerator and denominator polynomials |
|
7011
|
163 [r, p] = rresidue (b, a, varargin{1}, toler, e); |
|
6964
|
164 return |
|
7011
|
165 endif |
|
559
|
166 |
|
2303
|
167 ## Make sure both polynomials are in reduced form. |
|
1025
|
168 |
|
|
169 a = polyreduce (a); |
|
|
170 b = polyreduce (b); |
|
559
|
171 |
|
1025
|
172 b = b / a(1); |
|
|
173 a = a / a(1); |
|
559
|
174 |
|
1025
|
175 la = length (a); |
|
|
176 lb = length (b); |
|
559
|
177 |
|
2303
|
178 ## Handle special cases here. |
|
1025
|
179 |
|
|
180 if (la == 0 || lb == 0) |
|
559
|
181 k = r = p = e = []; |
|
|
182 return; |
|
|
183 elseif (la == 1) |
|
1025
|
184 k = b / a; |
|
559
|
185 r = p = e = []; |
|
|
186 return; |
|
|
187 endif |
|
|
188 |
|
2303
|
189 ## Find the poles. |
|
1025
|
190 |
|
|
191 p = roots (a); |
|
|
192 lp = length (p); |
|
559
|
193 |
|
6964
|
194 ## Sort poles so that multiplicity loop will work. |
|
|
195 |
|
|
196 [e, indx] = mpoles (p, toler, 1); |
|
|
197 p = p (indx); |
|
559
|
198 |
|
7398
|
199 ## For each group of pole multiplicity, set the value of each |
|
|
200 ## pole to the average of the group. This reduces the error in |
|
|
201 ## the resulting poles. |
|
|
202 |
|
|
203 p_group = cumsum (e == 1); |
|
|
204 for ng = 1:p_group(end) |
|
|
205 m = find (p_group == ng); |
|
|
206 p(m) = mean (p(m)); |
|
|
207 endfor |
|
|
208 |
|
2303
|
209 ## Find the direct term if there is one. |
|
1025
|
210 |
|
|
211 if (lb >= la) |
|
6964
|
212 ## Also return the reduced numerator. |
|
1025
|
213 [k, b] = deconv (b, a); |
|
|
214 lb = length (b); |
|
559
|
215 else |
|
|
216 k = []; |
|
|
217 endif |
|
|
218 |
|
7398
|
219 ## Determine if the poles are (effectively) zero. |
|
|
220 |
|
|
221 small = max (abs (p)); |
|
|
222 small = max ([small, 1]) * eps*1e4 * (1 + numel (p))^2; |
|
|
223 p(abs (p) < small) = 0; |
|
|
224 |
|
|
225 ## Determine if the poles are (effectively) real, or imaginary. |
|
|
226 |
|
|
227 index = (abs (imag (p)) < small); |
|
|
228 p(index) = real (p(index)); |
|
|
229 index = (abs (real (p)) < small); |
|
|
230 p(index) = 1i * imag (p(index)); |
|
|
231 |
|
|
232 ## The remainder determines the residues. The case of one pole |
|
|
233 ## is trivial. |
|
|
234 |
|
1025
|
235 if (lp == 1) |
|
|
236 r = polyval (b, p); |
|
559
|
237 return; |
|
|
238 endif |
|
|
239 |
|
6964
|
240 ## Determine the order of the denominator and remaining numerator. |
|
|
241 ## With the direct term removed the potential order of the numerator |
|
|
242 ## is one less than the order of the denominator. |
|
1025
|
243 |
|
6964
|
244 aorder = numel (a) - 1; |
|
|
245 border = aorder - 1; |
|
1025
|
246 |
|
6964
|
247 ## Construct a system of equations relating the individual |
|
|
248 ## contributions from each residue to the complete numerator. |
|
559
|
249 |
|
6964
|
250 A = zeros (border+1, border+1); |
|
|
251 B = prepad (reshape (b, [numel(b), 1]), border+1, 0); |
|
|
252 for ip = 1:numel(p) |
|
|
253 ri = zeros (size (p)); |
|
|
254 ri(ip) = 1; |
|
|
255 A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).'; |
|
|
256 endfor |
|
559
|
257 |
|
2303
|
258 ## Solve for the residues. |
|
1025
|
259 |
|
|
260 r = A \ B; |
|
559
|
261 |
|
|
262 endfunction |
|
6964
|
263 |
|
7011
|
264 function [pnum, pden, e] = rresidue (r, p, k, toler, e) |
|
6964
|
265 |
|
|
266 ## Reconstitute the numerator and denominator polynomials from the |
|
|
267 ## residues, poles, and direct term. |
|
|
268 |
|
7011
|
269 if (nargin < 2 || nargin > 5) |
|
6964
|
270 print_usage (); |
|
|
271 endif |
|
|
272 |
|
7011
|
273 if (nargin < 5) |
|
|
274 e = []; |
|
|
275 endif |
|
|
276 |
|
6964
|
277 if (nargin < 4) |
|
|
278 toler = []; |
|
|
279 endif |
|
|
280 |
|
|
281 if (nargin < 3) |
|
|
282 k = []; |
|
|
283 endif |
|
7011
|
284 |
|
|
285 if numel (e) |
|
|
286 indx = 1:numel(p); |
|
|
287 else |
|
|
288 [e, indx] = mpoles (p, toler, 0); |
|
|
289 p = p (indx); |
|
|
290 r = r (indx); |
|
|
291 endif |
|
6964
|
292 |
|
|
293 indx = 1:numel(p); |
|
|
294 |
|
|
295 for n = indx |
|
|
296 pn = [1, -p(n)]; |
|
|
297 if n == 1 |
|
|
298 pden = pn; |
|
|
299 else |
|
|
300 pden = conv (pden, pn); |
|
|
301 endif |
|
|
302 endfor |
|
|
303 |
|
|
304 ## D is the order of the denominator |
|
|
305 ## K is the order of the direct polynomial |
|
|
306 ## N is the order of the resulting numerator |
|
|
307 ## pnum(1:(N+1)) is the numerator's polynomial |
|
|
308 ## pden(1:(D+1)) is the denominator's polynomial |
|
|
309 ## pm is the multible pole for the nth residue |
|
|
310 ## pn is the numerator contribution for the nth residue |
|
|
311 |
|
|
312 D = numel (pden) - 1; |
|
|
313 K = numel (k) - 1; |
|
|
314 N = K + D; |
|
|
315 pnum = zeros (1, N+1); |
|
7011
|
316 for n = indx(abs (r) > 0) |
|
6964
|
317 p1 = [1, -p(n)]; |
|
7011
|
318 for m = 1:e(n) |
|
|
319 if (m == 1) |
|
6964
|
320 pm = p1; |
|
|
321 else |
|
|
322 pm = conv (pm, p1); |
|
|
323 endif |
|
|
324 endfor |
|
|
325 pn = deconv (pden, pm); |
|
|
326 pn = r(n) * pn; |
|
7183
|
327 pnum = pnum + prepad (pn, N+1, 0, 2); |
|
6964
|
328 endfor |
|
|
329 |
|
|
330 ## Add the direct term. |
|
|
331 |
|
|
332 if (numel (k)) |
|
|
333 pnum = pnum + conv (pden, k); |
|
|
334 endif |
|
|
335 |
|
|
336 ## Check for leading zeros and trim the polynomial coefficients. |
|
|
337 |
|
|
338 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; |
|
|
339 |
|
7011
|
340 pnum(abs (pnum) < small) = 0; |
|
|
341 pden(abs (pden) < small) = 0; |
|
6964
|
342 |
|
|
343 pnum = polyreduce (pnum); |
|
|
344 pden = polyreduce (pden); |
|
|
345 |
|
|
346 endfunction |
|
6968
|
347 |
|
|
348 %!test |
|
|
349 %! b = [1, 1, 1]; |
|
|
350 %! a = [1, -5, 8, -4]; |
|
|
351 %! [r, p, k, e] = residue (b, a); |
|
7398
|
352 %! assert (abs (r - [-2; 7; 3]) < 1e-12 |
|
|
353 %! && abs (p - [2; 2; 1]) < 1e-12 |
|
6998
|
354 %! && isempty (k) |
|
|
355 %! && e == [1; 2; 1]); |
|
6994
|
356 %! k = [1 0]; |
|
7011
|
357 %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); |
|
|
358 %! a = a; |
|
|
359 %! [br, ar] = residue (r, p, k); |
|
|
360 %! assert ((abs (br - b) < 1e-12 |
|
|
361 %! && abs (ar - a) < 1e-12)); |
|
|
362 %! [br, ar] = residue (r, p, k, e); |
|
|
363 %! assert ((abs (br - b) < 1e-12 |
|
|
364 %! && abs (ar - a) < 1e-12)); |
|
6994
|
365 |
|
|
366 %!test |
|
|
367 %! b = [1, 0, 1]; |
|
|
368 %! a = [1, 0, 18, 0, 81]; |
|
7398
|
369 %! [r, p, k, e] = residue (b, a); |
|
6994
|
370 %! r1 = [-5i; 12; +5i; 12]/54; |
|
|
371 %! p1 = [+3i; +3i; -3i; -3i]; |
|
7398
|
372 %! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12 |
|
6998
|
373 %! && isempty (k) |
|
|
374 %! && e == [1; 2; 1; 2]); |
|
6994
|
375 %! [br, ar] = residue (r, p, k); |
|
|
376 %! assert ((abs (br - b) < 1e-12 |
|
6998
|
377 %! && abs (ar - a) < 1e-12)); |
|
7011
|
378 |
|
|
379 %!test |
|
|
380 %! r = [7; 3; -2]; |
|
|
381 %! p = [2; 1; 2]; |
|
|
382 %! k = [1 0]; |
|
|
383 %! e = [2; 1; 1]; |
|
|
384 %! [b, a] = residue (r, p, k, e); |
|
|
385 %! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12 |
|
|
386 %! && abs (a - [1, -5, 8, -4]) < 1e-12)); |
|
|
387 %! [rr, pr, kr, er] = residue (b, a); |
|
|
388 %! [jnk, n] = mpoles(p); |
|
7398
|
389 %! assert ((abs (rr - r(n)) < 1e-12 |
|
|
390 %! && abs (pr - p(n)) < 1e-12 |
|
7011
|
391 %! && abs (kr - k) < 1e-12 |
|
|
392 %! && abs (er - e(n)) < 1e-12)); |
|
|
393 |
|
7188
|
394 %!test |
|
|
395 %! b = [1]; |
|
|
396 %! a = [1, 10, 25]; |
|
7398
|
397 %! [r, p, k, e] = residue (b, a); |
|
7188
|
398 %! r1 = [0; 1]; |
|
|
399 %! p1 = [-5; -5]; |
|
7398
|
400 %! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12 |
|
7188
|
401 %! && isempty (k) |
|
|
402 %! && e == [1; 2]); |
|
|
403 %! [br, ar] = residue (r, p, k); |
|
|
404 %! assert ((abs (br - b) < 1e-12 |
|
|
405 %! && abs (ar - a) < 1e-12)); |
