From bug-octave-request at che dot utexas dot edu  Mon Jun 20 17:44:54 1994
Subject: polynomial toolbox
From: amr at mpl dot UCSD dot EDU (Tony Richardson)
To: bug-octave at che dot utexas dot edu
Date: Mon, 20 Jun 94 15:44:49 pdt

John,
In response to the query on help-octave I cleaned up my polynomial
toolbox.  I've attached it to this message and you have my permission
to include it with octave if you'd like.

The included routines are: roots, poly, conv, deconv, filter, filter2,
compan, polyval, and polyvalm.

Comments:

 Conv, deconv, and filter are the basis of a signal processing toolbox
 also.  I'm continuing work in this area, primarily in the area of
 filter design.

 filter is the backbone routine of the polynomial routines (and also of
 any signal processing toolbox).  The speed of several of the routines
 would be greatly improved if filter were converted to a built-in routine.

 filter2 is an alternative implementation of filter.  It's slightly
 faster in the tests that I've run, but is not fully compatible with
 filter (it uses a different set of state variables).

 Matlab routines of the same name (all routines except filter2 have Matlab
 counterparts) are believed to be compatible.  The routines are in no way
 based (even in part) on Matlab routines however.  I do have access to
 Matlab, but I avoided looking at the Matlab code.  I only used Matlab to
 test for compatibility at the user level.

Regards and keep up the good work!

Tony Richardson
amr at mpl dot ucsd dot edu

#!/bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #!/bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	roots.m
#	poly.m
#	conv.m
#	deconv.m
#	compan.m
#	filter.m
#	filter2.m
#	polyval.m
#	polyvalm.m
# This archive created: Mon Jun 20 15:26:32 1994
export PATH; PATH=/bin:$PATH
if test -f 'roots.m'
then
	echo shar: over-writing existing file "'roots.m'"
fi
cat << \SHAR_EOF > 'roots.m'
function r = roots(c)
#Find the roots of a polynomial.
#
#In octave, a polynomial is represented by it's coefficients (arranged
#in descending order). For example, a vector c of length n+1 corresponds
#to the following nth order polynomial
#
#  p(x) = c(1) x^n + ... + c(n) x + c(n+1).
#
#roots(c) will return a vector r of length n such that
#
#  p(x) = c(1) [ (x - r(1)) * (x - r(2)) * ... * (x - r(n)) ]
#
#roots and poly are inverse functions to within a scaling factor.
#
#See also: poly, compan

  if(nargin != 1)
    error("usage: roots(vector)");
  endif

  if(!is_vector(c))
    error("roots: expecting a vector argument.");
  endif

  # Ensure that c is a row vector.
  if(rows(c) > 1)
    c = c.';
  endif

  # We could replace this with a call to compan, but it's faster to
  # just reproduce the code here.
  n = length(c);
  A = diag(ones(n-2,1),-1);
  A(1,:) = -c(2:n)/c(1);

  r = eig(A);

endfunction
SHAR_EOF
if test -f 'poly.m'
then
	echo shar: over-writing existing file "'poly.m'"
fi
cat << \SHAR_EOF > 'poly.m'
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, compan

  if(nargin != 1)
    error("usage: roots(argument)");
  endif

  if (is_scalar(r))
    c = [1 -r];
    return;
  elseif (is_square(r))
    r = eig(r);
  elseif (is_matrix(r))
    error("poly: matrix argument must be square.");
  endif

  l = length(r) + 1;
  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
SHAR_EOF
if test -f 'conv.m'
then
	echo shar: over-writing existing file "'conv.m'"
fi
cat << \SHAR_EOF > 'conv.m'
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: poly, roots, filter.

  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.
  a = reshape(a,1,la);
  b = reshape(b,1,lb);

  # Use the shortest vector as the coefficent vector to filter.
  if (la < lb)
    x = [b zeros(1,ly-lb)];
    y = filter(a,1,x);
  else
    x = [a zeros(1,ly-la)];
    y = filter(b,1,x);
  endif

endfunction
SHAR_EOF
if test -f 'deconv.m'
then
	echo shar: over-writing existing file "'deconv.m'"
fi
cat << \SHAR_EOF > 'deconv.m'
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, filter.

  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;

  b = filter(y,a,[1 zeros(1,ly-la)]);

  if (nargout == 2)
    r = y - conv(a,b);
  endif
endfunction
SHAR_EOF
if test -f 'compan.m'
then
	echo shar: over-writing existing file "'compan.m'"
fi
cat << \SHAR_EOF > 'compan.m'
function A = compan(c)
#compan (c)
#Compute the companion matrix corresponding to polynomial vector c.
#
#In octave a polynomial is represented by it's coefficients (arranged
#in descending order). For example a vector c of length n+1 corresponds
#to the following nth order polynomial
#
#  p(x) = c(1) x^n + ... + c(n) x + c(n+1).
#
#The corresponding companion matrix is
#        _                                                        _
#       |  -c(2)/c(1)   -c(3)/c(1)  ...  -c(n)/c(1)  -c(n+1)/c(1)  |
#       |       1            0      ...       0             0      |
#       |       0            1      ...       0             0      |
#   A = |       .            .   .            .             .      |
#       |       .            .       .        .             .      |
#       |       .            .           .    .             .      |
#       |_      0            0      ...       1             0     _|
#
#The eigenvalues of the companion matrix are equal to the roots of the
#polynomial.
#
#See also: roots, poly

  if(nargin != 1)
    error("usage: compan(vector)");
  endif

  if(!is_vector(c))
    error("compan: expecting a vector argument.");
  endif

  # Ensure that c is a row vector.
  if(rows(c) > 1)
    c = c.';
  endif

  n = length(c);
  A = diag(ones(n-2,1),-1);
  A(1,:) = -c(2:n)/c(1);

endfunction
SHAR_EOF
if test -f 'filter.m'
then
	echo shar: over-writing existing file "'filter.m'"
fi
cat << \SHAR_EOF > 'filter.m'
function [y, w] = filter(b,a,x,w)
#Filter a vector.
#y = filter(b,a,x) returns the solution to the following linear,
#time-invariant difference equation:
#
#   N                   M
#  sum a(k+1) y(n-k) + sum b(k+1) x(n-k) = 0  for 1<=n<=length(x)
#  k=0                 k=0
#
#where N=length(a)-1 and M=length(b)-1. An equivalent form of this
#equation is:
#
#          N                   M
#  y(n) = sum c(k+1) y(n-k) + sum d(k+1) x(n-k)  for 1<=n<=length(x)
#         k=1                 k=0
#
#where c = a/a(1) and d = b/a(1).
#
#In terms of the z-transform, y is the result of passing the discrete-
#time signal x through a system characterized by the following rational
#system function:
#
#             M
#            sum d(k+1) z^(-k)
#            k=0
#  H(z) = ----------------------
#               N
#          1 + sum c(k+1) z(-k)
#              k=1
#
#[y, sf] = filter(b,a,x,si) sets the initial state of the system, si,
#and returns the final state, sf.  The state vector is a column vector
#whose length is equal to the length of the longest coefficient vector.
#If si is not set, the initial state vector is set to all zeros.
#
#The particular algorithm employed is known as (in the signal
#processing literature) a transposed Direct Form II implementation.
#
#This routine is compatible with the Matlab routine of the same name.
#The filter2 algorithm is slightly faster, but uses a different set
#of state variables, hence the state vectors returned from filter2 are
#not the same as returned from this routine or from Matlab.

  if(nargin < 3 || nargin > 4)
    error("usage: [y, sf] = filter(b,a,x[,si])");
  endif

  if(is_matrix(a) || is_matrix(b) || is_matrix(x))
    error("Argument must be a vector.");
  endif

  N = length(a);
  M = length(b);
  L = length(x);

  MN = max([N M]);
  lw = MN - 1;

  # Ensure that all vectors have the assumed dimension.
  a = reshape(a,N,1);
  b = reshape(b,M,1);

  if(nargin == 3)
    # Set initial state to zero.
    w = zeros(lw,1);
  else
    if(is_matrix(w) || length(w) != lw)
      error("state vector has the wrong dimensions.");
    endif
    w = reshape(w,lw,1);
  endif

  # Allocate space for result.
  y = zeros(1,L);

  # It's convenient to pad both vectors to the same length.
  if(MN > M)
    b = [ b; zeros(MN-M,1) ];
  endif

  norm = a(1);
  if (norm != 1.)
    b = b/norm;
  endif

  # Distinguish between IIR and FIR cases.  The IIR code can easily be made
  # to  work for both cases, but the FIR code is slightly faster when it can
  # be used.

  if (N > 1)
    # IIR filter.
    # Reverse a coefficient vector and discard a(1) (it's always equal
    # to one after the normalization.
    if (MN > N)
      a = [ a; zeros(MN-N,1) ];
    endif
    if (norm != 1.)
      a = a/norm;
    endif
    for index = 1:L
      y(index) = w(1) + b(1)*x(index);
      # Update state vector
      w(1:(lw-1)) = w(2:lw) - a(2:lw)*y(index) + b(2:lw)*x(index);
      w(lw) = b(MN)*x(index) - a(MN) * y(index);
    endfor
  else
    for index = 1:L
      y(index) = w(1) + b(1)*x(index);
      # Update state vector
      w(1:lw-1) = w(2:lw) + b(2:lw)*x(index);
      w(lw) = b(MN)*x(index);
    endfor
  endif
endfunction
SHAR_EOF
if test -f 'filter2.m'
then
	echo shar: over-writing existing file "'filter2.m'"
fi
cat << \SHAR_EOF > 'filter2.m'
function [y, w] = filter2(b,a,x,w)
#Filter a vector.
#y = filter2(b,a,x) returns the solution to the following linear,
#time-invariant difference equation:
#
#   N                   M
#  sum a(k+1) y(n-k) + sum b(k+1) x(n-k) = 0  for 1<=n<=length(x)
#  k=0                 k=0
#
#where N=length(a)-1 and M=length(b)-1. An equivalent form of this
#equation is:
#
#          N                   M
#  y(n) = sum c(k+1) y(n-k) + sum d(k+1) x(n-k)  for 1<=n<=length(x)
#         k=1                 k=0
#
#where c = a/a(1) and d = b/a(1).
#
#In terms of the z-transform, y is the result of passing the discrete-
#time signal x through a system characterized by the following rational
#system function:
#
#             M
#            sum d(k+1) z^(-k)
#            k=0
#  H(z) = ----------------------
#               N
#          1 + sum c(k+1) z(-k)
#              k=1
#
#[y, sf] = filter2(b,a,x,si) sets the initial state of the system, si,
#and returns the final state, sf.  The state vector is a column vector
#whose length is equal to the length of the longest coefficient vector.
#If si is not set, the initial state vector is set to all zeros.
#
#The particular algorithm employed is known as (in the signal
#processing literature) a Direct Form II implementation.
#
#This routine is *not* fully compatible with the Matlab routine of the
#same name. It uses a different algorithm and a different set of state
#variables than the Matlab routine.  Refer to the filter routine for
#a fully compatible Matlab routine.
  if(nargin < 3 || nargin > 4)
    error("usage: [y, sf] = filter2(b,a,x[,si])");
  endif

  if(is_matrix(a) || is_matrix(b) || is_matrix(x))
    error("Argument must be a vector.");
  endif

  N = length(a);
  M = length(b);
  L = length(x);

  MN = max([N M]);

  # Ensure that all vectors have the assumed dimension.
  a = reshape(a,1,N);
  b = reshape(b,1,M);
  x = reshape(x,1,L);

  if(nargin == 3)
    # Set initial state to zero.
    w = zeros(MN,1);
  else
    if(is_matrix(w) || length(w) != MN)
      error("state vector has the wrong dimensions.");
    endif
    w = reshape(w,MN,1);
  endif

  # Allocate space for result.
  y = zeros(1,L);

  # It's convenient to work with the reversed coefficient vectors.
  b = b(M:-1:1);

  norm = a(1);

  if ( norm != 1. )
    b = b/norm;
  endif

  # Distinguish between IIR and FIR cases.  The IIR code can easily be made
  # to  work for both cases, but the FIR code is slightly faster when it can
  # be used.

  if (N > 1)
    # IIR filter.
    # Reverse a coefficient vector and discard a(1) (it's always equal
    # to one after the normalization.
    a = a(N:-1:2);
    if ( norm != 1. )
      a = a/norm;
    endif
    for index = 1:L
      # Shift state vector
      w(1:(MN-1)) = w(2:MN);
      # and update it.
      w(MN) = x(index) - a * w(1:(MN-1));
  
      # Calculate new output.
      y(index) = b * w((MN-M+1):MN);

    endfor
  else
    # FIR filter.
    xi = [ w; x.' ];
    for index = 1:L
      y(index) = b * xi((index+1):(index+M));
    endfor
    w = xi((index+1):(index+M));
  endif
endfunction
SHAR_EOF
if test -f 'polyval.m'
then
	echo shar: over-writing existing file "'polyval.m'"
fi
cat << \SHAR_EOF > 'polyval.m'
function y = polyval(c,x)
#Evaluate a polynomial.
#
#In octave, a polynomial is represented by it's coefficients (arranged
#in descending order). For example a vector c of length n+1 corresponds
#to the following nth order polynomial
#
#  p(x) = c(1) x^n + ... + c(n) x + c(n+1).
#
#polyval(c,x) will evaluate the polynomial at the specified value of x.
#
#If x is a vector or matrix, the polynomial is evaluated at each of the
#elements of x.
#
#See also: polyvalm, poly, roots, compan

  if(nargin != 2)
    error("usage: polyval(c,x)");
  endif

  if(is_matrix(c))
    error("poly: first argument must be a vector.");
  endif

  n = length(c);
  y = c(1)*ones(rows(x),columns(x));
  for index = 2:n
    y = c(index) + x .* y;
  endfor
endfunction
SHAR_EOF
if test -f 'polyvalm.m'
then
	echo shar: over-writing existing file "'polyvalm.m'"
fi
cat << \SHAR_EOF > 'polyvalm.m'
function y = polyvalm(c,x)
#Evaluate a polynomial in the matrix sense.
#
#In octave, a polynomial is represented by it's coefficients (arranged
#in descending order). For example a vector c of length n+1 corresponds
#to the following nth order polynomial
#
#  p(x) = c(1) x^n + ... + c(n) x + c(n+1).
#
#polyvalm(c,X) will evaluate the polynomial in the matrix sense, i.e. matrix
#multiplication is used instead of element by element multiplication as is
#used in polyval.
#
#X must be a square matrix.
#
#See also: polyval, poly, roots, compan

  if(nargin != 2)
    error("usage: polyvalm(c,x)");
  endif

  if(is_matrix(c))
    error("poly: first argument must be a vector.");
  endif

  if(!is_square(x))
    error("poly: second argument must be a square matrix.");
  endif

  [v, d] = eig(x);

  y = v * diag(polyval(c,diag(d))) * v';

endfunction
SHAR_EOF
#	End of shell archive
exit 0