diff -r 933b074f600e libinterp/corefcn/fft.cc
--- a/libinterp/corefcn/fft.cc	Sun Nov 20 21:42:27 2016 -0800
+++ b/libinterp/corefcn/fft.cc	Wed Nov 30 04:01:02 2016 +0800
@@ -59,26 +59,33 @@ do_fft (const octave_value_list &args, c
       if (! args(1).is_empty ())
         {
           double dval = args(1).double_value ();
-          if (octave::math::isnan (dval))
-            error ("%s: number of points (N) cannot be NaN", fcn);
+          if (octave::math::isnan (dval) || !octave::math::finite (dval))
+            error ("%s: number of points (N) cannot be NaN or Inf", fcn);
 
           n_points = octave::math::nint_big (dval);
           if (n_points < 0)
-            error ("%s: number of points (N) must be greater than zero", fcn);
+            error ("%s: number of points (N) must be non-negative", fcn);
         }
     }
 
   if (nargin > 2)
     {
       double dval = args(2).double_value ();
-      if (octave::math::isnan (dval))
-        error ("%s: DIM cannot be NaN", fcn);
-      else if (dval < 1 || dval > ndims)
+      if (octave::math::isnan (dval) || !octave::math::finite (dval))
+        error ("%s: DIM cannot be NaN or Inf", fcn);
+      else if (dval < 1)
         error ("%s: DIM must be a valid dimension along which to perform FFT",
                fcn);
       else
         // to be safe, cast it back to int since dim is an int
         dim = octave::math::nint (dval) - 1;
+
+      // Allow dim >= ndims, for degenerated cases, e.g. fft(a(:,:,1:1), n, 3)
+      if (dim >= ndims)
+        {
+          dims = dims.redim (dim + 1);
+          ndims = dims.ndims ();
+        }
     }
 
   // FIXME: This seems strange and unnecessary (10/21/16).
@@ -110,16 +117,6 @@ do_fft (const octave_value_list &args, c
         return octave_value (NDArray (dims));
     }
 
-  if (n_points == 1)
-    {
-      octave_value_list idx (ndims);
-      for (octave_idx_type i = 0; i < ndims; i++)
-        idx(i) = idx_vector::colon;
-      idx(dim) = idx_vector (static_cast<octave_idx_type> (0));
-
-      return arg.do_index_op (idx);
-    }
-
   if (arg.is_single_type ())
     {
       if (arg.is_real_type ())
@@ -169,6 +166,11 @@ do_fft (const octave_value_list &args, c
 %!assert (fft (ones (2,2)), [2,2; 0,0])
 %!assert (fft (eye (2,2)), [1,1; 1,-1])
 
+%!assert (fft ([1 2;3 4], 0, 3), [1 2;3 4](:,:,1:0))
+%!assert (fft ([1 2;3 4], 1, 3), [1 2;3 4])
+%!assert (fft ([1 2;3 4], 2, 3), repmat ([1 2;3 4], [1 1 2]))
+%!assert (fft (1, 3, 3), ones (1,1,3))
+
 %!assert (fft (single ([])), single ([]))
 %!assert (fft (zeros (10,0,"single")), zeros (10,0,"single"))
 %!assert (fft (zeros (0,10,"single")), zeros (0,10,"single"))
@@ -177,7 +179,15 @@ do_fft (const octave_value_list &args, c
 %!assert (fft (ones (2,2,"single")), single ([2,2; 0,0]))
 %!assert (fft (eye (2,2,"single")), single ([1,1; 1,-1]))
 
+%!assert (fft (single ([1 2;3 4]), 0, 3), single ([1 2;3 4](:,:,1:0)))
+%!assert (fft (single ([1 2;3 4]), 1, 3), single ([1 2;3 4]))
+%!assert (fft (single ([1 2;3 4]), 2, 3), single (repmat ([1 2;3 4], [1 1 2])))
+%!assert (fft (single (1), 3, 3), ones (1,1,3,"single"))
+
+
 %!error (fft ())
+%!error (fft (1,2,Inf))
+%!error (fft (1,Inf,1))
 */
 
 
diff -r 933b074f600e liboctave/array/CNDArray.cc
--- a/liboctave/array/CNDArray.cc	Sun Nov 20 21:42:27 2016 -0800
+++ b/liboctave/array/CNDArray.cc	Wed Nov 30 04:01:02 2016 +0800
@@ -59,18 +59,26 @@ ComplexNDArray::fourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   const Complex *in (fortran_vec ());
@@ -90,18 +98,26 @@ ComplexNDArray::ifourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   const Complex *in (fortran_vec ());
@@ -195,9 +211,17 @@ ComplexNDArray::fourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   ComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
@@ -243,9 +267,17 @@ ComplexNDArray::ifourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   ComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
diff -r 933b074f600e liboctave/array/dNDArray.cc
--- a/liboctave/array/dNDArray.cc	Sun Nov 20 21:42:27 2016 -0800
+++ b/liboctave/array/dNDArray.cc	Wed Nov 30 04:01:02 2016 +0800
@@ -101,18 +101,26 @@ NDArray::fourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   const double *in (fortran_vec ());
@@ -132,18 +140,26 @@ NDArray::ifourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   ComplexNDArray retval (*this);
@@ -236,9 +252,17 @@ NDArray::fourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   ComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
@@ -284,9 +308,17 @@ NDArray::ifourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return ComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      ComplexNDArray retval (*this);
+      return retval;
+    }
+
   ComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
diff -r 933b074f600e liboctave/array/fCNDArray.cc
--- a/liboctave/array/fCNDArray.cc	Sun Nov 20 21:42:27 2016 -0800
+++ b/liboctave/array/fCNDArray.cc	Wed Nov 30 04:01:02 2016 +0800
@@ -59,18 +59,26 @@ FloatComplexNDArray::fourier (int dim) c
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   const FloatComplex *in (fortran_vec ());
@@ -90,18 +98,26 @@ FloatComplexNDArray::ifourier (int dim) 
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   const FloatComplex *in (fortran_vec ());
@@ -195,9 +211,17 @@ FloatComplexNDArray::fourier (int dim) c
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   FloatComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
@@ -242,9 +266,17 @@ FloatComplexNDArray::ifourier (int dim) 
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   FloatComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
diff -r 933b074f600e liboctave/array/fNDArray.cc
--- a/liboctave/array/fNDArray.cc	Sun Nov 20 21:42:27 2016 -0800
+++ b/liboctave/array/fNDArray.cc	Wed Nov 30 04:01:02 2016 +0800
@@ -59,18 +59,26 @@ FloatNDArray::fourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   const float *in (fortran_vec ());
@@ -90,18 +98,26 @@ FloatNDArray::ifourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
-  octave_idx_type howmany = numel () / dv(dim);
+  octave_idx_type howmany = numel () / n;
   howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
-  octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
   octave_idx_type dist = (stride == 1 ? n : 1);
 
   FloatComplexNDArray retval (*this);
@@ -194,9 +210,17 @@ FloatNDArray::fourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the Fourier transform result as if dims () has been extended
+      // to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   FloatComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
@@ -241,9 +265,17 @@ FloatNDArray::ifourier (int dim) const
 {
   dim_vector dv = dims ();
 
-  if (dim > dv.ndims () || dim < 0)
+  if (dim < 0)
     return FloatComplexNDArray ();
 
+  if (dim >= dv.ndims ())
+    {
+      // Return the inverse Fourier transform result as if dims () has been
+      // extended to length dim.
+      FloatComplexNDArray retval (*this);
+      return retval;
+    }
+
   FloatComplexNDArray retval (dv);
   octave_idx_type npts = dv(dim);
   octave_idx_type nn = 4*npts+15;
