*** ./DLD-FUNCTIONS/filter.cc.orig	2004-03-26 12:21:58.000000000 +0100
--- ./DLD-FUNCTIONS/filter.cc	2004-03-26 15:49:32.000000000 +0100
***************
*** 39,170 ****
  #include "oct-obj.h"
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&);
  
! extern MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&);
  #endif
  
  template <class T>
! MArray<T>
! filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si)
  {
!   MArray<T> y;
  
    int a_len  = a.length ();
    int b_len  = b.length ();
-   int x_len  = x.length ();
- 
-   int si_len = si.length ();
  
    int ab_len = a_len > b_len ? a_len : b_len;
  
    b.resize (ab_len, 0.0);
  
!   if (si.length () != ab_len - 1)
      {
!       error ("filter: si must be a vector of length max (length (a), length (b)) - 1");
        return y;
      }
  
!   T norm = a (0);
  
!   if (norm == 0.0)
      {
!       error ("filter: the first element of a must be non-zero");
        return y;
      }
  
!   y.resize (x_len, 0.0);
  
!   if (norm != 1.0)
!     b = b / norm;
  
!   if (a_len > 1)
      {
!       a.resize (ab_len, 0.0);
  
!       if (norm != 1.0)
! 	a = a / norm;
  
!       for (int i = 0; i < x_len; i++)
  	{
! 	  y (i) = si (0) + b (0) * x (i);
  
! 	  if (si_len > 1)
  	    {
! 	      for (int j = 0; j < si_len - 1; j++)
  		{
! 		  OCTAVE_QUIT;
  
! 		  si (j) = si (j+1) - a (j+1) * y (i)
! 		    + b (j+1) * x (i);
  		}
! 
! 	      si (si_len-1) = b (si_len) * x (i)
! 		- a (si_len) * y (i);
  	    }
- 	  else
- 	    si (0) = b (si_len) * x (i)
- 	      - a (si_len) * y (i);
  	}
!     }
!   else if (si_len > 0)
!     {
!       for (int i = 0; i < x_len; i++)
  	{
! 	  y (i) = si (0) + b (0) * x (i);
! 
! 	  if (si_len > 1)
  	    {
! 	      for (int j = 0; j < si_len - 1; j++)
  		{
! 		  OCTAVE_QUIT;
  
! 		  si (j) = si (j+1) + b (j+1) * x (i);
  		}
! 
! 	      si (si_len-1) = b (si_len) * x (i);
  	    }
- 	  else
- 	    si (0) = b (1) * x (i);
  	}
      }
!   else
!     y = b (0) * x;
! 
    return y;
  }
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&,
! 	MArray<double>&);
! 
! extern MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
! 	MArray<Complex>&);
  #endif
  
  template <class T>
! MArray<T>
! filter (MArray<T>& b, MArray<T>& a, MArray<T>& x)
  {
    int a_len = a.length ();
    int b_len = b.length ();
  
    int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
!   MArray<T> si (si_len, T (0.0));
! 
!   return filter (b, a, x, si);
  }
  
  DEFUN_DLD (filter, args, nargout,
    "-*- texinfo -*-\n\
  @deftypefn {Loadable Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\
  @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\
  Return the solution to the following linear, time-invariant difference\n\
  equation:\n\
  @iftex\n\
--- 39,252 ----
  #include "oct-obj.h"
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
  
! extern MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
  #endif
  
  template <class T>
! MArrayN<T>
! filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, MArrayN<T>& si, 
! 	int dim = 0)
  {
!   MArrayN<T> y;
  
    int a_len  = a.length ();
    int b_len  = b.length ();
  
    int ab_len = a_len > b_len ? a_len : b_len;
  
    b.resize (ab_len, 0.0);
+   if (a_len > 1)
+     a.resize (ab_len, 0.0);
  
!   T norm = a (0);
! 
!   if (norm == 0.0)
      {
!       error ("filter: the first element of a must be non-zero");
        return y;
      }
  
!   dim_vector x_dims = x.dims ();
!   if ((dim < 0) || (dim > x_dims.length ()))
!     {
!       error ("filter: filtering over invalid dimension");
!       return y;
!     }
  
!   int x_len = x_dims (dim);
! 
!   dim_vector si_dims = si.dims ();
!   int si_len = si_dims (0);
! 
!   if (si_len != ab_len - 1)
      {
!       error ("filter: first dimension of si must be of length max (length (a), length (b)) - 1");
        return y;
      }
  
!   if (si_dims.length () != x_dims.length ())
!     {
!       error ("filter: dimensionality of si and x must agree");
!       return y;
!     }
  
!   int si_dim = 0;
!   for (int i = 0; i < x_dims.length (); i++)
!     {
!       if (i == dim)
! 	continue;
!       
!       if (si_dims (++si_dim) != x_dims (i))
! 	{
! 	  error ("filter: dimensionality of si and x must agree");
! 	  return y;
! 	}
!     }
  
!   if (norm != 1.0)
      {
!       a = a / norm;
!       b = b / norm;
!     }
  
!   if ((a_len <= 1) && (si_len <= 1))
!     return b(0) * x;
  
!   y.resize (x_dims, 0.0);
! 
!   int x_stride = 1;
!   for (int i = 0; i < dim; i++)
!     x_stride *= x_dims(i);
! 
!   int x_num = x_dims.numel () / x_len;
!   for (int num = 0; num < x_num; num++)
!     {
!       int x_offset;
!       if (x_stride == 1)
! 	x_offset = num * x_len;
!       else
  	{
! 	  int x_offset2 = 0;
! 	  x_offset = num;
! 	  while (x_offset >= x_stride)
! 	    {
! 	      x_offset -= x_stride;
! 	      x_offset2++;
! 	    }
! 	  x_offset += x_offset2 * x_stride * x_len;
! 	}
!       int si_offset = num * si_len;
  
!       if (a_len > 1)
! 	{
! 	  for (int i = 0; i < x_len; i++)
  	    {
! 	      int idx = i * x_stride + x_offset; 
! 	      y (idx) = si (si_offset) + b (0) * x (idx);
! 
! 	      if (si_len > 1)
  		{
! 		  for (int j = 0; j < si_len - 1; j++)
! 		    {
! 		      OCTAVE_QUIT;
! 
! 		      si (j + si_offset) =  si (j + 1 + si_offset) - 
! 			a (j+1) * y (idx) + b (j+1) * x (idx);
! 		    }
  
! 		  si (si_len - 1 + si_offset) = b (si_len) * x (idx)
! 		    - a (si_len) * y (idx);
  		}
! 	      else
! 		si (si_offset) = b (si_len) * x (idx)
! 		  - a (si_len) * y (idx);
  	    }
  	}
!       else if (si_len > 0)
  	{
! 	  for (int i = 0; i < x_len; i++)
  	    {
! 	      int idx = i * x_stride + x_offset; 
! 	      y (idx) = si (si_offset) + b (0) * x (idx);
! 
! 	      if (si_len > 1)
  		{
! 		  for (int j = 0; j < si_len - 1; j++)
! 		    {
! 		      OCTAVE_QUIT;
! 
! 		      si (j + si_offset) = si (j + 1 + si_offset) + 
! 			b (j+1) * x (idx);
! 		    }
  
! 		  si (si_len - 1 + si_offset) = b (si_len) * x (idx);
  		}
! 	      else
! 		si (si_offset) = b (1) * x (idx);
  	    }
  	}
      }
!   
    return y;
  }
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
! 	MArrayN<double>&, int dim);
! 
! extern MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
! 	MArrayN<Complex>&, int dim);
  #endif
  
  template <class T>
! MArrayN<T>
! filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, int dim = -1)
  {
+   dim_vector x_dims = x.dims ();
+ 
+   if (dim < 0)
+     {
+       // Find first non-singleton dimension
+       while ((dim < x_dims.length()) && (x_dims (dim) <= 1))
+ 	dim++;
+   
+       // All dimensions singleton, pick first dimension
+       if (dim == x_dims.length ())
+ 	dim = 0;
+     }
+   else
+     if (dim < 0 || dim > x_dims.length ())
+       {
+ 	error ("filter: filtering over invalid dimension");
+ 	return MArrayN<T> ();
+       }
+ 
    int a_len = a.length ();
    int b_len = b.length ();
  
    int si_len = (a_len > b_len ? a_len : b_len) - 1;
+   dim_vector si_dims = x.dims ();
+   for (int i = dim; i > 0; i--)
+     si_dims (i) = si_dims (i-1);
+   si_dims (0) = si_len;
+   
+   MArrayN<T> si (si_dims, T (0.0));
  
!   return filter (b, a, x, si, dim);
  }
  
  DEFUN_DLD (filter, args, nargout,
    "-*- texinfo -*-\n\
  @deftypefn {Loadable Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\
  @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\
+ @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, [], @var{dim})\n\
+ @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si}, @var{dim})\n\
  Return the solution to the following linear, time-invariant difference\n\
  equation:\n\
  @iftex\n\
***************
*** 194,200 ****
   $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\
  @end tex\n\
  @end iftex\n\
! An equivalent form of this equation is:\n\
  @iftex\n\
  @tex\n\
  $$\n\
--- 276,283 ----
   $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\
  @end tex\n\
  @end iftex\n\
! over the first non-singleton dimension of @var{x} or over @var{dim} if\n\
! supplied. An equivalent form of this equation is:\n\
  @iftex\n\
  @tex\n\
  $$\n\
***************
*** 259,317 ****
  
    int nargin  = args.length ();
  
!   if (nargin < 3 || nargin > 4)
      {
        print_usage ("filter");
        return retval;
      }
  
!   const char *errmsg = "filter: arguments must be vectors";
  
!   bool x_is_row_vector = (args(2).rows () == 1);
  
!   bool si_is_row_vector = (nargin == 4 && args(3).rows () == 1);
  
    if (args(0).is_complex_type ()
        || args(1).is_complex_type ()
        || args(2).is_complex_type ()
!       || (nargin == 4 && args(3).is_complex_type ()))
      {
        ComplexColumnVector b (args(0).complex_vector_value ());
        ComplexColumnVector a (args(1).complex_vector_value ());
!       ComplexColumnVector x (args(2).complex_vector_value ());
  
        if (! error_state)
  	{
! 	  ComplexColumnVector si;
  
! 	  if (nargin == 3)
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      si.resize (si_len, 0.0);
  	    }
  	  else
! 	    si = ComplexColumnVector (args(3).complex_vector_value ());
  
  	  if (! error_state)
  	    {
! 	      ComplexColumnVector y (filter (b, a, x, si));
  
  	      if (nargout == 2)
! 		{
! 		  if (si_is_row_vector)
! 		    retval(1) = si.transpose ();
! 		  else
! 		    retval(1) = si;
! 		}
  
! 	      if (x_is_row_vector)
! 		retval(0) = y.transpose ();
! 	      else
! 		retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
--- 342,432 ----
  
    int nargin  = args.length ();
  
!   if (nargin < 3 || nargin > 5)
      {
        print_usage ("filter");
        return retval;
      }
  
!   const char *errmsg = "filter: arguments a and b must be vectors";
  
!   int dim;
!   dim_vector x_dims = args(2).dims ();
  
!   if (nargin == 5)
!     {
!       dim = args(4).nint_value() - 1;
!       if (dim < 0 || dim >= x_dims.length ())
! 	{
! 	  error ("filter: filtering over invalid dimension");
! 	  return retval;
! 	}
!     }
!   else
!     {
!       // Find first non-singleton dimension
!       dim = 0;
!       while ((dim < x_dims.length()) && (x_dims (dim) <= 1))
! 	dim++;
!   
!       // All dimensions singleton, pick first dimension
!       if (dim == x_dims.length ())
! 	dim = 0;
!     }
  
    if (args(0).is_complex_type ()
        || args(1).is_complex_type ()
        || args(2).is_complex_type ()
!       || (nargin >= 4 && args(3).is_complex_type ()))
      {
        ComplexColumnVector b (args(0).complex_vector_value ());
        ComplexColumnVector a (args(1).complex_vector_value ());
! 
!       ComplexNDArray x (args(2).complex_array_value ());
  
        if (! error_state)
  	{
! 	  ComplexNDArray si;
  
! 	  if (nargin == 3 || args(3).is_empty ())
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      dim_vector si_dims = x.dims ();
! 	      for (int i = dim; i > 0; i--)
! 		si_dims (i) = si_dims (i-1);
! 	      si_dims (0) = si_len;
! 
! 	      si.resize (si_dims, 0.0);
  	    }
  	  else
! 	    {
! 	      dim_vector si_dims = args (3).dims ();
! 	      bool si_is_vector = true;
! 	      for (int i=0; i < si_dims.length (); i++)
! 		if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ()))
! 		  {
! 		    si_is_vector = false;
! 		    break;
! 		  }
! 
! 	      if (si_is_vector)
! 		si = ComplexNDArray (args(3).complex_vector_value ());
! 	      else
! 		si = args(3).complex_array_value ();
! 	    }
  
  	  if (! error_state)
  	    {
! 	      ComplexNDArray y (filter (b, a, x, si, dim));
  
  	      if (nargout == 2)
! 		retval(1) = si;
  
! 	      retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
***************
*** 323,362 ****
      {
        ColumnVector b (args(0).vector_value ());
        ColumnVector a (args(1).vector_value ());
!       ColumnVector x (args(2).vector_value ());
  
        if (! error_state)
  	{
! 	  ColumnVector si;
  
! 	  if (nargin == 3)
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      si.resize (si_len, 0.0);
  	    }
  	  else
! 	    si = ColumnVector (args(3).vector_value ());
  
  	  if (! error_state)
  	    {
! 	      ColumnVector y (filter (b, a, x, si));
  
  	      if (nargout == 2)
! 		{
! 		  if (si_is_row_vector)
! 		    retval(1) = si.transpose ();
! 		  else
! 		    retval(1) = si;
! 		}
  
! 	      if (x_is_row_vector)
! 		retval(0) = y.transpose ();
! 	      else
! 		retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
--- 438,489 ----
      {
        ColumnVector b (args(0).vector_value ());
        ColumnVector a (args(1).vector_value ());
! 
!       NDArray x (args(2).array_value ());
  
        if (! error_state)
  	{
! 	  NDArray si;
  
! 	  if (nargin == 3 || args(3).is_empty ())
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      dim_vector si_dims = x.dims ();
! 	      for (int i = dim; i > 0; i--)
! 		si_dims (i) = si_dims (i-1);
! 	      si_dims (0) = si_len;
! 
! 	      si.resize (si_dims, 0.0);
  	    }
  	  else
! 	    {
! 	      dim_vector si_dims = args (3).dims ();
! 	      bool si_is_vector = true;
! 	      for (int i=0; i < si_dims.length (); i++)
! 		if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ()))
! 		  {
! 		    si_is_vector = false;
! 		    break;
! 		  }
! 
! 	      if (si_is_vector)
! 		si = NDArray (args(3).vector_value ());
! 	      else
! 		si = args(3).array_value ();
! 	    }
  
  	  if (! error_state)
  	    {
! 	      NDArray y (filter (b, a, x, si, dim));
  
  	      if (nargout == 2)
! 		retval(1) = si;
  
! 	      retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
***************
*** 368,386 ****
    return retval;
  }
  
! template MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&,
! 	MArray<double>&);
! 
! template MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&);
! 
! template MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
! 	MArray<Complex>&);
  
! template MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray <Complex>&);
  
  /*
  ;;; Local Variables: ***
--- 495,513 ----
    return retval;
  }
  
! template MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
! 	MArrayN<double>&, int dim);
! 
! template MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
! 
! template MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
! 	MArrayN<Complex>&, int dim);
  
! template MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
  
  /*
  ;;; Local Variables: ***
*** ./ChangeLog.orig	2004-03-26 15:55:00.000000000 +0100
--- ./ChangeLog	2004-03-26 15:54:52.000000000 +0100
***************
*** 1,3 ****
--- 1,8 ----
+ 2004-03-26 David Bateman <address@hidden>
+ 
+ 	* DLD-FUNCTIONS/filter.cc: Convert for NDArray, better Matlab
+ 	compatibility.
+ 	
  2004-03-12  John W. Eaton  <address@hidden>
  
  	* ov-cell.cc (octave_cell::save_hdf5): Handle empty cells.