# HG changeset patch
# User Markus Meisinger <chloros2@gmx.de>
# Date 1641552146 -3600
#      Fri Jan 07 11:42:26 2022 +0100
# Node ID 0bc62bc9cce608144ab66f2f102fa208d870a8d1
# Parent  83f9f8bda88362cf47a4cc10f126543d362e2e32
ifft, ifft2, ifftn: Return the real part if "symmetric" is specified (bug #43742).

* fft.cc, fft2.cc, fftn.cc: Separate the symflag and return the real part.

diff -r 83f9f8bda883 -r 0bc62bc9cce6 libinterp/corefcn/fft.cc
--- a/libinterp/corefcn/fft.cc	Tue Dec 28 18:59:33 2021 -0500
+++ b/libinterp/corefcn/fft.cc	Fri Jan 07 11:42:26 2022 +0100
@@ -41,8 +41,9 @@
 do_fft (const octave_value_list& args, const char *fcn, int type)
 {
   int nargin = args.length ();
+  int numeric_argin = nargin;
 
-  if (nargin < 1 || nargin > 3)
+  if (nargin < 1 || nargin > 4)
     print_usage ();
 
   octave_value retval;
@@ -51,8 +52,28 @@
   dim_vector dims = arg.dims ();
   int ndims = dims.ndims ();
   int dim = -1;
+  bool ifft_symmetric = false;
 
-  if (nargin > 1)
+  // Separate the symflag
+  if (type && (nargin > 1) && args(nargin - 1).is_string ())
+    {
+      std::string strval = args(nargin - 1).
+                           xstring_value ("%s: string expected", fcn);
+      if (strval == "symmetric")
+        {
+          ifft_symmetric = true;
+        }
+      else if (strval != "nonsymmetric")
+        {
+          error ("%s: unknown symflag", fcn);
+        }
+      numeric_argin = nargin - 1;
+    }
+  else if (nargin == 4)
+    print_usage ();
+
+  // Get numeric arguments
+  if (numeric_argin > 1)
     {
       if (! args(1).isempty ())
         {
@@ -66,7 +87,7 @@
         }
     }
 
-  if (nargin > 2)
+  if (numeric_argin > 2)
     {
       double dval = args(2).double_value ();
       if (math::isnan (dval))
@@ -115,41 +136,69 @@
         idx(i) = idx_vector::colon;
       idx(dim) = idx_vector (0);
 
-      return arg.index_op (idx);
+      retval = arg.index_op (idx);
+      if (ifft_symmetric)
+        return retval.real();
+      else
+        return retval;
     }
 
   if (arg.is_single_type ())
     {
       if (arg.isreal ())
         {
-          FloatNDArray nda = arg.float_array_value ();
+          FloatNDArray in = arg.float_array_value ();
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
+          in.resize (dims, 0.0);
+          FloatComplexNDArray out =
+                           (type != 0 ? in.ifourier (dim) : in.fourier (dim));
+
+          if (ifft_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else
         {
-          FloatComplexNDArray cnda = arg.float_complex_array_value ();
+          FloatComplexNDArray in = arg.float_complex_array_value ();
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
+          in.resize (dims, 0.0);
+          FloatComplexNDArray out =
+                           (type != 0 ? in.ifourier (dim) : in.fourier (dim));
+
+          if (ifft_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
     }
   else
     {
       if (arg.isreal ())
         {
-          NDArray nda = arg.array_value ();
+          NDArray in = arg.array_value ();
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
+          in.resize (dims, 0.0);
+          ComplexNDArray out =
+                           (type != 0 ? in.ifourier (dim) : in.fourier (dim));
+
+          if (ifft_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else if (arg.iscomplex ())
         {
-          ComplexNDArray cnda = arg.complex_array_value ();
+          ComplexNDArray in = arg.complex_array_value ();
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
+          in.resize (dims, 0.0);
+          ComplexNDArray out =
+                           (type != 0 ? in.ifourier (dim) : in.fourier (dim));
+
+          if (ifft_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else
         err_wrong_type_arg (fcn, arg);
@@ -227,6 +276,7 @@
 @deftypefn  {} {} ifft (@var{x})
 @deftypefnx {} {} ifft (@var{x}, @var{n})
 @deftypefnx {} {} ifft (@var{x}, @var{n}, @var{dim})
+@deftypefnx {} {} ifft (@dots{}, @var{symflag})
 Compute the inverse discrete Fourier transform of @var{x}
 using a Fast Fourier Transform (FFT) algorithm.
 
@@ -244,6 +294,10 @@
 
 If called with three arguments, @var{dim} is an integer specifying the
 dimension of the matrix along which the inverse FFT is performed.
+
+If the Fourier transform satisfies Hermitian redundancy, it's inverse is
+real-valued. To force a real-valued output, the keyword
+@code{@var{symmetric}} can be passed as an argument.
 @seealso{fft, ifft2, ifftn, fftw}
 @end deftypefn */)
 {
@@ -251,6 +305,41 @@
 }
 
 /*
+%!testif HAVE_FFTW
+%! assert (isreal (ifft (i, 'symmetric')), true);
+
+%!testif HAVE_FFTW
+%! assert (iscomplex (ifft (i, 'nonsymmetric')), true);
+
+%!testif HAVE_FFTW
+%! assert (iscomplex (ifft ([4, 2+2i, 4, 2-2i] + 1e-3i)), true);
+
+%!testif HAVE_FFTW
+%! s = ifft ([4, 2+2i, 4, 2-2i] + 1e-3i, 'symmetric');
+%! assert (s, [3 -1 1 1], 1e-3);
+%! assert (isreal (s), true);
+
+%!testif HAVE_FFTW
+%! s = ifft ([4, 2+2i, 4, 2-2i] + 1e-3i, 'nonsymmetric');
+%! assert (s, [3 -1 1 1], 1e-3);
+%! assert (iscomplex (s), true);
+
+%!testif HAVE_FFTW
+%! s = ifft (single ([4, 2+2i, 4, 2-2i] + 1e-3i), 'symmetric');
+%! assert (s, [3 -1 1 1], 1e-3);
+%! assert (isreal (s), true);
+
+%!testif HAVE_FFTW
+%! s = ifft (single ([4, 2+2i, 4, 2-2i] + 1e-3i), 'nonsymmetric');
+%! assert (s, [3 -1 1 1], 1e-3);
+%! assert (iscomplex (s), true);
+
+%!testif HAVE_FFTW 
+%!error <unknown symflag> ifft (1, 'string');
+
+%!testif HAVE_FFTW
+%! assert (numel (ifft ([4, 2+2i, 4, 2-2i], [])), 4);
+
 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
 ##         Comalco Research and Technology
 ##         02 May 2000
diff -r 83f9f8bda883 -r 0bc62bc9cce6 libinterp/corefcn/fft2.cc
--- a/libinterp/corefcn/fft2.cc	Tue Dec 28 18:59:33 2021 -0500
+++ b/libinterp/corefcn/fft2.cc	Fri Jan 07 11:42:26 2022 +0100
@@ -43,16 +43,37 @@
 do_fft2 (const octave_value_list& args, const char *fcn, int type)
 {
   int nargin = args.length ();
+  int numeric_argin = nargin;
 
-  if (nargin < 1 || nargin > 3)
+  if (nargin < 1 || nargin > 4)
     print_usage ();
 
   octave_value retval;
   octave_value arg = args(0);
   dim_vector dims = arg.dims ();
   octave_idx_type n_rows = -1;
+  bool ifft2_symmetric = false;
 
-  if (nargin > 1)
+  // Separate the symflag
+  if (type && (nargin > 1) && args(nargin - 1).is_string ())
+    {
+      std::string strval = args(nargin - 1).
+                                   xstring_value ("%s: string expected", fcn);
+      if (strval == "symmetric")
+        {
+          ifft2_symmetric = true;
+        }
+      else if (strval != "nonsymmetric")
+        {
+          error ("%s: unknown symflag", fcn);
+        }
+      numeric_argin = nargin - 1;
+    }
+  else if (nargin == 4)
+    print_usage ();
+
+  // Get numeric arguments
+  if (numeric_argin > 1)
     {
       double dval = args(1).double_value ();
       if (math::isnan (dval))
@@ -64,7 +85,7 @@
     }
 
   octave_idx_type n_cols = -1;
-  if (nargin > 2)
+  if (numeric_argin > 2)
     {
       double dval = args(2).double_value ();
       if (math::isnan (dval))
@@ -101,34 +122,58 @@
     {
       if (arg.isreal ())
         {
-          FloatNDArray nda = arg.float_array_value ();
+          FloatNDArray in = arg.float_array_value ();
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
+          in.resize (dims, 0.0);
+          FloatComplexNDArray out =
+                             (type != 0 ? in.ifourier2d () : in.fourier2d ());
+
+          if (ifft2_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else
         {
-          FloatComplexNDArray cnda = arg.float_complex_array_value ();
+          FloatComplexNDArray in = arg.float_complex_array_value ();
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
+          in.resize (dims, 0.0);
+          FloatComplexNDArray out =
+                             (type != 0 ? in.ifourier2d () : in.fourier2d ());
+
+          if (ifft2_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
     }
   else
     {
       if (arg.isreal ())
         {
-          NDArray nda = arg.array_value ();
+          NDArray in = arg.array_value ();
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
+          in.resize (dims, 0.0);
+          ComplexNDArray out =
+                             (type != 0 ? in.ifourier2d () : in.fourier2d ());
+
+          if (ifft2_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else if (arg.iscomplex ())
         {
-          ComplexNDArray cnda = arg.complex_array_value ();
+          ComplexNDArray in = arg.complex_array_value ();
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
+          in.resize (dims, 0.0);
+          ComplexNDArray out =
+                             (type != 0 ? in.ifourier2d () : in.fourier2d ());
+
+          if (ifft2_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else
         err_wrong_type_arg (fcn, arg);
@@ -161,6 +206,7 @@
        doc: /* -*- texinfo -*-
 @deftypefn  {} {} ifft2 (@var{A})
 @deftypefnx {} {} ifft2 (@var{A}, @var{m}, @var{n})
+@deftypefnx {} {} ifft2 (@dots{}, @var{symflag})
 Compute the inverse two-dimensional discrete Fourier transform of @var{A}
 using a Fast Fourier Transform (FFT) algorithm.
 
@@ -170,6 +216,10 @@
 
 If @var{A} is a multi-dimensional matrix, each two-dimensional sub-matrix
 of @var{A} is treated separately.
+
+If the Fourier transform satisfies Hermitian redundancy, it's inverse is
+real-valued. To force a real-valued output, the keyword
+@code{@var{symmetric}} can be passed as an argument.
 @seealso{fft2, ifft, ifftn, fftw}
 @end deftypefn */)
 {
@@ -177,6 +227,37 @@
 }
 
 /*
+%!testif HAVE_FFTW
+%!error <unknown symflag> ifft2 ([1, 3; 3, 1], 'test')
+
+%!testif HAVE_FFTW
+%! s1 = vander (1:4) + 1e-3i;
+%! S = fft2 (s1);
+%!
+%! s2 = ifft2 (S, 'symmetric');
+%! assert (isreal (s2), true);
+%! assert (s1, s2, 1e-3);
+%!
+%! s3 = ifft2 (S, 'nonsymmetric');
+%! assert (s1, s3, 0);
+%!
+%! s4 = ifft2 (S);
+%! assert (s1, s4, 0);
+
+%!testif HAVE_FFTW
+%! s1 = single (vander (1:4) + 1e-3i);
+%! S = fft2 (s1);
+%!
+%! s2 = ifft2 (S, 'symmetric');
+%! assert (isreal (s2), true);
+%! assert (s1, s2, 1e-3);
+%!
+%! s3 = ifft2 (S, 'nonsymmetric');
+%! assert (s1, s3, 0);
+%!
+%! s4 = ifft2 (S);
+%! assert (s1, s4, 0);
+
 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
 ##         Comalco Research and Technology
 ##         02 May 2000
diff -r 83f9f8bda883 -r 0bc62bc9cce6 libinterp/corefcn/fftn.cc
--- a/libinterp/corefcn/fftn.cc	Tue Dec 28 18:59:33 2021 -0500
+++ b/libinterp/corefcn/fftn.cc	Fri Jan 07 11:42:26 2022 +0100
@@ -43,19 +43,40 @@
 do_fftn (const octave_value_list& args, const char *fcn, int type)
 {
   int nargin = args.length ();
+  int numeric_argin = nargin;
 
-  if (nargin < 1 || nargin > 2)
+  if (nargin < 1 || nargin > 3)
     print_usage ();
 
   octave_value retval;
   octave_value arg = args(0);
   dim_vector dims = arg.dims ();
+  bool ifftn_symmetric = false;
 
   for (int i = 0; i < dims.ndims (); i++)
     if (dims(i) < 0)
       return retval;
 
-  if (nargin > 1)
+  // Separate the symflag
+  if (type && (nargin > 1) && args(nargin - 1).is_string ())
+    {
+      std::string strval = args(nargin - 1).
+                                   xstring_value ("%s: string expected", fcn);
+      if (strval == "symmetric")
+        {
+          ifftn_symmetric = true;
+        }
+      else if (strval != "nonsymmetric")
+        {
+          error ("%s: unknown symflag", fcn);
+        }
+      numeric_argin = nargin - 1;
+    }
+  else if (nargin == 3)
+    print_usage ();
+
+  // Get numeric arguments
+  if (numeric_argin > 1)
     {
       Matrix val = args(1).xmatrix_value ("%s: SIZE must be a vector of length dim", fcn);
 
@@ -88,34 +109,58 @@
     {
       if (arg.isreal ())
         {
-          FloatNDArray nda = arg.float_array_value ();
+          FloatNDArray in = arg.float_array_value ();
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
+          in.resize (dims, 0.0);
+          FloatComplexNDArray out =
+                             (type != 0 ? in.ifourierNd () : in.fourierNd ());
+
+          if (ifftn_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else
         {
-          FloatComplexNDArray cnda = arg.float_complex_array_value ();
+          FloatComplexNDArray in = arg.float_complex_array_value ();
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
+          in.resize (dims, 0.0);
+          FloatComplexNDArray out =
+                             (type != 0 ? in.ifourierNd () : in.fourierNd ());
+
+          if (ifftn_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
     }
   else
     {
       if (arg.isreal ())
         {
-          NDArray nda = arg.array_value ();
+          NDArray in = arg.array_value ();
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
+          in.resize (dims, 0.0);
+          ComplexNDArray out =
+                             (type != 0 ? in.ifourierNd () : in.fourierNd ());
+
+          if (ifftn_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else if (arg.iscomplex ())
         {
-          ComplexNDArray cnda = arg.complex_array_value ();
+          ComplexNDArray in = arg.complex_array_value ();
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
+          in.resize (dims, 0.0);
+          ComplexNDArray out =
+                             (type != 0 ? in.ifourierNd () : in.fourierNd ());
+
+          if (ifftn_symmetric)
+            retval = real (out);
+          else
+            retval = out;
         }
       else
         err_wrong_type_arg (fcn, arg);
@@ -147,6 +192,7 @@
        doc: /* -*- texinfo -*-
 @deftypefn  {} {} ifftn (@var{A})
 @deftypefnx {} {} ifftn (@var{A}, @var{size})
+@deftypefnx {} {} ifftn (@dots{}, @var{symflag})
 Compute the inverse N-dimensional discrete Fourier transform of @var{A}
 using a Fast Fourier Transform (FFT) algorithm.
 
@@ -156,10 +202,47 @@
 truncated prior to performing the inverse FFT@.  Otherwise, if an element of
 @var{size} is larger than the corresponding dimension then @var{A} is
 resized and padded with zeros.
+
+If the Fourier transform satisfies Hermitian redundancy, it's inverse is
+real-valued. To force a real-valued output, the keyword
+@code{@var{symmetric}} can be passed as an argument.
 @seealso{fftn, ifft, ifft2, fftw}
 @end deftypefn */)
 {
   return do_fftn (args, "ifftn", 1);
 }
 
+/*
+%!testif HAVE_FFTW
+%!error <unknown symflag> ifftn (1, 'test');
+
+%!testif HAVE_FFTW
+%! s1 = reshape (1:8,2,2,2) + 1e-3i;
+%! S = fftn (s1);
+%!
+%! s2 = ifftn (S, 'symmetric');
+%! assert (isreal (s2), true);
+%! assert (s1, s2, 1e-3);
+%!
+%! s3 = ifftn (S, 'nonsymmetric');
+%! assert (s1, s3, 0);
+%!
+%! s4 = ifftn (S);
+%! assert (s1, s4, 0);
+
+%!testif HAVE_FFTW
+%! s1 = single (reshape (1:8,2,2,2) + 1e-3i);
+%! S = fftn (s1);
+%!
+%! s2 = ifftn (S, 'symmetric');
+%! assert (isreal (s2), true);
+%! assert (s1, s2, 1e-3);
+%!
+%! s3 = ifftn (S, 'nonsymmetric');
+%! assert (s1, s3, 0);
+%!
+%! s4 = ifftn (S);
+%! assert (s1, s4, 0);
+*/
+
 OCTAVE_NAMESPACE_END
