# HG changeset patch
# User jwe
# Date 775175528 0
# Node ID e79ff1f4df3cf863228f3cc815512f9e79aad96f
# Parent  34c713db72b684ebfdeea2a7e98f1731bd37424f
[project @ 1994-07-25 22:32:08 by jwe]
Initial revision

diff --git a/scripts/general/postpad.m b/scripts/general/postpad.m
new file mode 100644
--- /dev/null
+++ b/scripts/general/postpad.m
@@ -0,0 +1,42 @@
+function y = postpad(x,l,c)
+#postpad(x,l)
+#Appends zeros to the vector x until it is of length l.
+#postpad(x,l,c) appends the constant c instead of zero.
+#
+#If length(x) > l, elements from the end of x are removed
+#until a vector of length l is obtained.
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  if(nargin == 2)
+    c = 0;
+  elseif(nargin<2 || nargin>3)
+    error("usage: postpad(x,l) or postpad(x,l,c)");
+  endif
+
+  if(is_matrix(x))
+    error("first argument must be a vector");
+  elseif(!is_scalar(l))
+    error("second argument must be a scaler");
+  endif
+
+  if(l<0)
+    error("second argument must be non-negative");
+  endif
+
+  lx = length(x);
+
+  if(lx >= l)
+    y = x(1:l);
+  else
+    if(rows(x)>1)
+      y = [ x; c*ones(l-lx,1) ];
+    else
+      y = [ x c*ones(1,l-lx) ];
+    endif
+  endif
+
+endfunction
diff --git a/scripts/polynomial/compan.m b/scripts/polynomial/compan.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/compan.m
@@ -0,0 +1,48 @@
+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: poly, roots, residue, conv, deconv, polyval, polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  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
diff --git a/scripts/polynomial/polyderiv.m b/scripts/polynomial/polyderiv.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/polyderiv.m
@@ -0,0 +1,33 @@
+function p = polyderiv(p)
+#polyderiv(c)
+#Returns the coefficients of the derivative of the polynomial whose
+#coefficients are given by vector c.
+#
+#SEE ALSO: poly, polyinteg, polyreduce, roots, conv, deconv, residue,
+#          filter, polyval, polyvalm
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  if(nargin != 1)
+    error("usage: polyderiv(vector)");
+  endif
+
+  if(is_matrix(p))
+    error("argument must be a vector");
+  endif
+
+  lp = length(p);
+  if(lp == 1)
+    p = 0;
+    return;
+  elseif (lp == 0)
+    p = [];
+    return;
+  end
+
+  p = p(1:(lp-1)) .* [(lp-1):-1:1];
+
+endfunction
diff --git a/scripts/polynomial/polyreduce.m b/scripts/polynomial/polyreduce.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/polyreduce.m
@@ -0,0 +1,29 @@
+function p = polyreduce(p)
+#polyreduce(c)
+#Reduces a polynomial coefficient vector to a minimum number of terms,
+#i.e. it strips off any leading zeros.
+#
+#SEE ALSO: poly, roots, conv, deconv, residue, filter, polyval, polyvalm,
+#          polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  index = find(p==0);
+
+  index = find(index == 1:length(index));
+
+  if (length(index) == 0)
+    return;
+  endif
+
+  if(length(p)>1)
+    p = p(index(length(index))+1:length(p));
+  endif
+
+  if(length(p)==0)
+    p = 0;
+  endif
+endfunction
diff --git a/scripts/polynomial/polyval.m b/scripts/polynomial/polyval.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/polyval.m
@@ -0,0 +1,41 @@
+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, conv, deconv, residue, filter,
+#          polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  if(nargin != 2)
+    error("usage: polyval(c,x)");
+  endif
+
+  if(is_matrix(c))
+    error("poly: first argument must be a vector.");
+  endif
+
+  if(length(c) == 0)
+    y = c;
+    return;
+  endif
+
+  n = length(c);
+  y = c(1)*ones(rows(x),columns(x));
+  for index = 2:n
+    y = c(index) + x .* y;
+  endfor
+endfunction
diff --git a/scripts/polynomial/polyvalm.m b/scripts/polynomial/polyvalm.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/polyvalm.m
@@ -0,0 +1,40 @@
+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, conv, deconv, residue, filter,
+#          polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  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
diff --git a/scripts/signal/filter.m b/scripts/signal/filter.m
new file mode 100644
--- /dev/null
+++ b/scripts/signal/filter.m
@@ -0,0 +1,135 @@
+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
+#minus one.  If si is not set, the initial state vector is set to all
+#zeros.
+#
+#The particular algorithm employed is known as a transposed Direct
+#Form II implementation.
+#
+#SEE ALSO: poly, roots, conv, deconv, residue, polyval, polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  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.
+  if(columns(a) > 1)
+    a = reshape(a,N,1);
+  endif
+  if(columns(b) > 1)
+    b = reshape(b,M,1);
+  endif
+
+  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
+    if(columns(w) > 1)
+      w = reshape(w,lw,1);
+    endif
+  endif
+
+  # Allocate space for result.
+  y = zeros(1,L);
+
+  # It's convenient to pad the coefficient vectors to the same length.
+  b = postpad(b,MN);
+
+  norm = a(1);
+  if (norm == 0.)
+    error("First element in second argument must be non-zero.");
+  endif
+
+  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.
+    a = postpad(a,MN);
+    if (norm != 1.)
+      a = a/norm;
+    endif
+    for index = 1:L
+      y(index) = w(1) + b(1)*x(index);
+      # Update state vector
+      if(lw > 1)
+        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);
+      else
+        w(1) = b(MN)*x(index) - a(MN) * y(index);
+      endif
+    endfor
+  else
+    # FIR filter.
+    if(lw > 0)
+      for index = 1:L
+        y(index) = w(1) + b(1)*x(index);
+        if ( lw > 1)
+          # Update state vector
+          w(1:lw-1) = w(2:lw) + b(2:lw)*x(index);
+          w(lw) = b(MN)*x(index);
+        else
+          w(1) = b(2)*x(index);
+        endif
+      endfor
+    else
+      # Handle special case where there is no delay separately.
+      y = b(1)*x;
+    endif
+  endif
+endfunction