--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/conv.m	Tue Oct 11 00:34:39 1994 +0000
@@ -0,0 +1,55 @@
+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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/deconv.m	Tue Oct 11 00:34:39 1994 +0000
@@ -0,0 +1,52 @@
+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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/poly.m	Tue Oct 11 00:34:39 1994 +0000
@@ -0,0 +1,33 @@
+function y = poly (x)
+  #
+  # If A is a square n-by-n matrix, poly (A) is the row vector of 
+  # the coefficients of det (z * eye(n) - A), the characteristic
+  # polynomial of 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 (all (imag (x) == 0)))
+    y = real (y);
+  endif
+  
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/roots.m	Tue Oct 11 00:34:39 1994 +0000
@@ -0,0 +1,38 @@
+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
+  
+  [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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/signal/fftconv.m	Tue Oct 11 00:34:39 1994 +0000
@@ -0,0 +1,41 @@
+function c = fftconv (a, b, N)
+
+  # usage:  fftconv (a, b [, N])
+  #
+  # c = fftconv (a, b) returns the convolution of the vectors a and b,
+  # a vector with length equal to length (a) + length (b) - 1.  
+  # If a and b are the coefficient vectors of two polynomials, c is
+  # the coefficient vector of the product polynomial.
+  #
+  # The computation uses the FFT by calling fftfilt.  If the optional
+  # argument N is specified, an N-point FFT is used.
+
+  # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994
+  # Copyright Dept of Statistics and Probability Theory TU Wien
+  
+  if (nargin < 2 || nargin > 3)
+    error ("usage:  fftconv (b, x [, N])");
+  endif
+  
+  if (is_matrix (a) || is_matrix (b))
+    error ("fftconv:  both a and b should be vectors");
+  endif
+  la = length (a);
+  lb = length (b);
+  if ((la == 1) || (lb == 1))
+    c = a * b;
+  else
+    lc = la + lb - 1;
+    a(lc) = 0;
+    b(lc) = 0;
+    if (nargin == 2)
+      c = fftfilt (a, b);
+    else
+      if !(is_scalar (N))
+	error ("fftconv:  N has to be a scalar");
+      endif
+      c = fftfilt (a, b, N);
+    endif
+  endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/signal/fftfilt.m	Tue Oct 11 00:34:39 1994 +0000
@@ -0,0 +1,82 @@
+function y = fftfilt (b, x, N)
+
+  # usage:  fftfilt (b, x [, N])
+  #
+  # y = fftfilt (b, x) filters x with the FIR filter b using the FFT.
+  # y = fftfilt (b, x, N) uses the overlap-add method to filter x with
+  # b using an N-point FFT.
+
+  # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994
+  # Copyright Dept of Statistics and Probability Theory TU Wien
+
+  # Reference:  Oppenheim & Schafer (1989).  Discrete-time Signal
+  # Processing (Chapter 8).  Prentice-Hall.
+  
+  # If N is not specified explicitly, we do not use the overlap-add 
+  # method at all because loops are really slow.  Otherwise, we only
+  # ensure that the number of points in the FFT is the smallest power
+  # of two larger than N and length(b).  This could result in length
+  # one blocks, but if the user knows better ...
+  
+  if (nargin < 2 || nargin > 3)
+    error ("usage:  fftfilt (b, x [, N])");
+  endif
+  
+  [r_x, c_x] = size (x);
+  [r_b, c_b] = size (b);
+  if !( (min ([r_x c_x]) == 1) || (min ([r_b c_b]) == 1) )
+    error ("fftfilt:  both x and b should be vectors.");
+  endif
+  l_x  = r_x * c_x;
+  l_b  = r_b * c_b;
+
+  if ((l_x == 1) && (l_b == 1))
+    y = b * x;  
+    return;
+  endif
+  
+  x    = reshape (x, 1, l_x);
+  b    = reshape (b, 1, l_b);
+  
+  if (nargin == 2)
+    # use FFT with the smallest power of 2 which is >= length (x) +
+    # length (b) - 1 as number of points ...
+    N    = 2^(ceil (log (l_x + l_b - 1) / log(2)));
+    y    = ifft (fft (x, N) .* fft(b, N));
+  else
+    # use overlap-add method ...
+    if !(is_scalar (N))
+      error ("fftfilt:  N has to be a scalar");
+    endif
+    N    = 2^(ceil (log (max ([N l_b])) / log(2)));
+    L    = N - l_b + 1;
+    B    = fft (b, N);
+    R    = ceil (l_x / L);
+    y    = zeros (1, l_x);
+    for r=1:R;
+      lo  = (r - 1) * L + 1;
+      hi  = min (r * L, l_x);
+      tmp = ifft (fft (x(lo:hi), N) .* B);
+      hi  = min (lo+N-1, l_x);
+      y(lo:hi) = y(lo:hi) + tmp(1:(hi-lo+1));
+    endfor  
+  endif
+    
+  y    = reshape (y(1:l_x), r_x, c_x);
+
+  # final cleanups:  if both x and b are real respectively integer, y
+  # should also be so
+  if !(any (imag (x)) || any (imag (b)))
+    y = real (y);
+  endif
+  if !(any (x - round (x)) || any (b - round (b)))
+    y = round (y);
+  endif
+
+endfunction
+
+
+
+
+
+