# HG changeset patch
# User jwe
# Date 775168725 0
# Node ID faf108b99d21df4a19cbe55d1d80c9970cdf0667
# Parent  e6ae50d8d4eb6337bf0324201550914b478223b3
[project @ 1994-07-25 20:38:45 by jwe]
Initial revision

diff --git a/scripts/polynomial/conv-amr.m b/scripts/polynomial/conv-amr.m
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
diff --git a/scripts/polynomial/conv-tuwien.m b/scripts/polynomial/conv-tuwien.m
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
diff --git a/scripts/polynomial/deconv-amr.m b/scripts/polynomial/deconv-amr.m
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
diff --git a/scripts/polynomial/deconv-tuwien.m b/scripts/polynomial/deconv-tuwien.m
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
+    
+
+
+
+
+
+
+
+
diff --git a/scripts/polynomial/poly-amr.m b/scripts/polynomial/poly-amr.m
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
diff --git a/scripts/polynomial/poly-tuwien.m b/scripts/polynomial/poly-tuwien.m
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
diff --git a/scripts/polynomial/roots-tuwien.m b/scripts/polynomial/roots-tuwien.m
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