Mercurial > octave-nkf

comparison scripts/polynomial/polyderiv.m @ 5217:e88886a6934d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-03-16 20:03:01 by jwe]
author jwe
date Wed, 16 Mar 2005 20:03:01 +0000
parents 5ed60b8b1ac4
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5216:5ed60b8b1ac4 5217:e88886a6934d
31 ## filter, polygcd, polyval, and polyvalm} 31 ## filter, polygcd, polyval, and polyvalm}
32 32
33 ## Author: Tony Richardson <arichard@stark.cc.oh.us> 33 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
34 ## Created: June 1994 34 ## Created: June 1994
35 ## Adapted-By: jwe 35 ## Adapted-By: jwe
36 ## Paul Kienzle <pkienzle@kienzle.powernet.co.uk>
37 ## handle b/a and b*a
38 36
39 function [q, r] = polyderiv (p, a) 37 function [q, r] = polyderiv (p, a)
40 38
41 if (nargin < 1 || nargin > 3) 39 if (nargin == 1 || nargin == 2)
42 usage ("q=polyderiv(p) or q=polyderiv(b,a) or [q, r]=polyderiv(b,a)"); 40 if (! isvector (p))
41 error ("polyderiv: argument must be a vector");
42 endif
43 if (nargin == 2)
44 if (! isvector (a))
45 error ("polyderiv: argument must be a vector");
46 endif
47 if (nargout == 1)
48 ## derivative of p*a returns a single polynomial
49 q = polyderiv (conv (p, a));
50 else
51 ## derivative of p/a returns numerator and denominator
52 r = conv (a, a);
53 if (numel (p) == 1)
54 q = -p * polyderiv (a);
55 elseif (numel (a) == 1)
56 q = a * polyderiv (p);
57 else
58 q = conv (polyderiv (p), a) - conv (p, polyderiv (a));
59 q = polyreduce (q);
60 endif
61
62 ## remove common factors from numerator and denominator
63 x = polygcd (q, r);
64 if (length(x) != 1)
65 q = deconv (q, x);
66 r = deconv (r, x);
67 endif
68
69 ## move all the gain into the numerator
70 q = q/r(1);
71 r = r/r(1);
72 endif
73 else
74 lp = numel (p);
75 if (lp == 1)
76 q = 0;
77 return;
78 elseif (lp == 0)
79 q = [];
80 return;
81 endif
82
83 ## Force P to be a row vector.
84 p = p(:).';
85
86 q = p(1:(lp-1)) .* [(lp-1):-1:1];
87 endif
88 else
89 usage ("q = polyderiv (p) or q = polyderiv (b, a) or [q, r] = polyderiv (b, a)");
43 endif 90 endif
44 91
45 if (! isvector (p))
46 error ("polyderiv: argument must be a vector");
47 endif
48
49 if (nargin == 2)
50 if (! isvector (a))
51 error ("polyderiv: argument must be a vector");
52 endif
53 if (nargout == 1)
54 ## derivative of p*a returns a single polynomial
55 q = polyderiv(conv(p,a));
56 else
57 ## derivative of p/a returns numerator and denominator
58 r = conv(a, a);
59 if numel(p) == 1
60 q = -p * polyderiv(a);
61 elseif numel(a) == 1
62 q = a * polyderiv(p);
63 else
64 q = conv(polyderiv(p),a) - conv(p,polyderiv(a));
65 q = polyreduce(q);
66 endif
67
68 ## remove common factors from numerator and denominator
69 x = polygcd(q,r);
70 if length(x)!=1
71 q=deconv(q,x);
72 r=deconv(r,x);
73 endif
74
75 ## move all the gain into the numerator
76 q=q/r(1);
77 r=r/r(1);
78 endif
79 else
80 lp = numel (p);
81 if (lp == 1)
82 q = 0;
83 return;
84 elseif (lp == 0)
85 q = [];
86 return;
87 end
88
89 ## Force P to be a row vector.
90 p = p(:).';
91
92 q = p (1:(lp-1)) .* [(lp-1):-1:1];
93 endif
94 92
95 endfunction 93 endfunction