Mercurial > hg > octave-terminal
changeset 558:faf108b99d21
[project @ 1994-07-25 20:38:45 by jwe]
Initial revision
author | jwe |
---|---|
date | Mon, 25 Jul 1994 20:38:45 +0000 |
parents | e6ae50d8d4eb |
children | 4e826edfbc56 |
files | scripts/polynomial/conv-amr.m scripts/polynomial/conv-tuwien.m scripts/polynomial/deconv-amr.m scripts/polynomial/deconv-tuwien.m scripts/polynomial/poly-amr.m scripts/polynomial/poly-tuwien.m scripts/polynomial/roots-tuwien.m |
diffstat | 7 files changed, 316 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/conv-amr.m @@ -0,0 +1,52 @@ +function y = conv(a,b) +#Convolve two vectors. +#y = conv(a,b) returns a vector of length equal to length(a)+length(b)-1. +#If a and b are polynomial coefficient vectors, conv returns the +#coefficients of the product polynomial. +# +#SEE ALSO: deconv, poly, roots, residue, polyval, polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 2) + error("usage: conv(a,b)"); + endif + + if(is_matrix(a) || is_matrix(b)) + error("conv: both arguments must be vectors"); + endif + + la = length(a); + lb = length(b); + + ly = la + lb - 1; + + # Ensure that both vectors are row vectors. + if(rows(a) > 1) + a = reshape(a,1,la); + endif + if(rows(b) > 1) + b = reshape(b,1,lb); + endif + + # Use the shortest vector as the coefficent vector to filter. + if (la < lb) + if(ly>lb) + x = [b zeros(1,ly-lb)]; + else + x = b; + endif + y = filter(a,1,x); + else + if(ly>la) + x = [a zeros(1,ly-la)]; + else + x = a; + endif + y = filter(b,1,x); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/conv-tuwien.m @@ -0,0 +1,43 @@ +function c=conv(a,b) +# +# usage: conv(a,b) +# +# Returns the convolution of vectors a and b. The resulting vector +# is of length(a)+length(b)-1. If a and b are polynomial coefficients +# conv(a,b) is equivalent to polynomial multiplication. + +# written by Gerhard Kircher on Aug 27, 1993 +# modified by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 23, 1993. + + l_a = length(a); + l_b = length(b); + return_row = 0; + if (l_a > l_b) + if (rows(a) == 1) + return_row = 1; + endif + else + if (rows(b) == 1) + return_row = 1; + endif + endif + a = reshape(a, l_a, 1); + b = reshape(b, l_b, 1); + if (l_a == 1 || l_b == 1) + c = a * b; + else + l_c = l_a + l_b - 1; + a(l_c) = 0; + b(l_c) = 0; + c = ifft(fft(a) .* fft(b)); + if !( any(imag(a)) || any(imag(b)) ) + c = real(c); + endif + if !( any(a-round(a)) || any(b-round(b)) ) + c = round(c); + endif + endif + if (return_row == 1) + c = c'; + end +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/deconv-amr.m @@ -0,0 +1,51 @@ +function [b, r] = deconv(y,a) +#Deconvolve two vectors. +# +#[b, r] = deconv(y,a) solves for b and r such that: +# +# y = conv(a,b) + r +# +#If y and a are polynomial coefficient vectors, b will contain the +#coefficients of the polynomial quotient and r will be a remander +#polynomial of lowest order. +# +#SEE ALSO: conv, poly, roots, residue, polyval, polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 2) + error("usage: deconv(y,a)"); + endif + + if(is_matrix(y) || is_matrix(a)) + error("conv: both arguments must be vectors"); + endif + + la = length(a); + ly = length(y); + + lb = ly - la + 1; + + if (ly>la) + b = filter(y,a,[1 zeros(1,ly-la)]); + elseif (ly == la) + b = filter(y,a,1); + else + b = 0; + endif + + b = polyreduce(b); + + lc = la + length(b) - 1; + if(ly == lc) + r = y - conv(a,b); + else + r = [ zeros(1,lc-ly) y] - conv(a,b); + endif + + r = polyreduce(r); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/deconv-tuwien.m @@ -0,0 +1,47 @@ +function [x, r] = deconv(a, b) +# +# Returns x and r such that a = conv(b, x) + r. +# If a and b are vectors of polynomial coefficients, x and r are the +# vectors of coefficients of quotient and remainder in the polynomial +# division of a by b. + +# written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 27, 1993 +# copyright Dept of Probability Theory and Statistics TU Wien + + if !(nargin == 2) + error("usage: deconv(a, b)"); + endif + f = find(b); + l_b = length(f); + if (l_b == 0) + error("deconv(a, b): b has to be nonzero"); + endif + l_a = length(a); + if (l_a < l_b) + x = 0; + r = a; + else + b = reshape(b(f(1):f(l_b)), 1, l_b); + a = reshape(a, 1, l_a); + # the stupid way: + if (l_b == 1) + x = a / b; + else + x(1) = a(1) ./ b(1); + for i = 2:l_a-l_b+1; + x(i) = (a(i) - b(i:-1:2) * x(1:i-1)) ./ b(1); + endfor + endif + r = a - conv(b, x); + endif + +endfunction + + + + + + + + +
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/poly-amr.m @@ -0,0 +1,54 @@ +function c = poly(r) +#Find the coefficients of a polynomial from its roots. +#Find the characteristic polynomial of a matrix. +# +#Given a vector r of length n containing the n roots of polynomial p(x) +# +# p(x) = (x - r(1)) * (x - r(2)) * ... * (x - r(n)) +# +#poly(r) will return a coefficient vector c of length n+1 such that +# +# p(x) = c(1) x^n + ... + c(n) x + c(n+1). +# +#and c(1) will always be equal to one. +# +#Given a matrix A, poly(A) will return a vector containing the coefficients +#of the characteristic polynomial of A, det(xI - A). +# +#poly and roots are inverse functions to within a scaling factor. +# +#SEE ALSO: roots, conv, deconv, residue, filter, polyval, polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 1) + error("usage: roots(argument)"); + endif + + l = length(r) + 1; + + if (is_scalar(r)) + c = [1 -r]; + return; + elseif(l == 1) + # r is an empty matrix. + # Matlab compatibility + c = 1; + return; + elseif (is_square(r)) + r = eig(r); + elseif (is_matrix(r)) + error("poly: matrix argument must be square."); + endif + + c = ones(1,l); + c(l) = -r(1); + for index = 2:(l-1) + m = l + 2 - index; + c((m-1):l) = [c(m:l) 0] - r(index)*[1 c(m:l)]; + endfor + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/poly-tuwien.m @@ -0,0 +1,32 @@ +function y = poly (x) +# +# If A is a square matrix, poly (A) is the row vector of coefficients of +# the characteristic polynomial det (z * eye(A) - A). +# If x is a vector, poly (x) is a vector of coefficients of the polynomial +# whose roots are the elements of x. + +# written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993 +# copyright Dept of Probability Theory and Statistics TU Wien + + m = min(size(x)); + n = max(size(x)); + if (m == 0) + y = 1; + elseif (m == 1) + v = x; + elseif (m == n) + v = eig(x); + else + error("usage: poly(x), where x is a vector or a square matrix"); + endif + + y = [ 1 zeros(1,n) ]; + for j = 1:n; + y(2:(j+1)) = y(2:(j+1)) - v(j) .* y(1:j); + endfor + + if all(imag(x) == 0) + y = real(y); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/roots-tuwien.m @@ -0,0 +1,37 @@ +function r = roots(v) +# +# For a vector v with n components, return the roots of the polynomial +# v(1)*z^(n-1) + ... + v(n-1) * z + v(n). + +# written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993 +# copyright Dept of Probability Theory and Statistics TU Wien +# modified by KH on Jan 10, 1994 + + [nr, nc] = size(v); + if !((nr == 1 && nc > 1) || (nc == 1 && nr > 1)) + error("usage: roots(v), where v is a nonzero vector"); + endif + n = nr + nc - 1; + v = reshape(v,1,n); + # If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the + # leading k zeros and n-k-l roots of the polynomial are zero. + f = find(v); + m = max(size(f)); + if (m > 0) + v = v(f(1):f(m)); + l = max(size(v)); + if (l > 1) + A = diag(ones(1, l-2), -1); + A(1,:) = -v(2:l) ./ v(1); + r = eig(A); + if (f(m) < n) + r = [r; zeros(n-f(m), 1)]; + endif + else + r = zeros(n-f(m), 1); + endif + else + error("usage: roots(v), where v is a nonzero vector"); + endif + +endfunction