# HG changeset patch
# User Markus Meisinger <chloros2@gmx.de>
# Date 1711305268 -3600
#      Sun Mar 24 19:34:28 2024 +0100
# Node ID 342e2e9bff3d07d76f20b9cba25475e894e2b320
# Parent  399be7cc310f80f94558fadaa3eceda7bd9d57b6
ifft, ifft2, ifftn: Implement "symmetric" flag for conjugate-symmetric input (bug #43742).

* liboctave/numeric/oct-fftw.cc:  Macros for SIMD alignment calculation
(SIMD_PAD, SIMD_ALIGN, SIMD_OFFSET).  Added double and float FFTW planner for
complex to real "C2R" transformations (fftw_planner::do_create_plan).  Added
interface functions for FFTW execution (ifftsym, ifftNdsym).  Added functions
to enforce conjugate-symmetry of input data (makesymmetric_Nd,
make_subspace_symmetric).  Added "C2R" input data conversion functions
(c2r_convert_1d, c2r_convert_Nd).

* liboctave/numeric/oct-fftw.h:  Plannner attributes and prototypes.

* libinterp/corefcn/fft.cc, libinterp/corefcn/fft2.cc,
libinterp/corefcn/fftn.cc:  Detect symflag in import arguments.  Added new
tests.  Update of the documentation.

* liboctave/array/CNDArray.cc, liboctave/array/CNDArray.h,
liboctave/array/dNDArray.cc, liboctave/array/dNDArray.h,
liboctave/array/fCNDArray.cc, liboctave/array/fCNDArray.h,
liboctave/array/fNDArray.cc, liboctave/array/fNDArray.h:  Ifft interface for
returning a real value (ifouriersym, ifourier2dsym, ifourierNdsym).

diff -r 399be7cc310f -r 342e2e9bff3d libinterp/corefcn/fft.cc
--- a/libinterp/corefcn/fft.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/libinterp/corefcn/fft.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -42,7 +42,7 @@
 {
   int nargin = args.length ();
 
-  if (nargin < 1 || nargin > 3)
+  if (nargin < 1 || nargin > 4)
     print_usage ();
 
   octave_value retval;
@@ -51,7 +51,23 @@
   dim_vector dims = arg.dims ();
   int ndims = dims.ndims ();
   int dim = -1;
+  bool ifft_symmetric = false;
 
+  // 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);
+      nargin--;
+    }
+  else if (nargin == 4)
+    print_usage ();
+
+  // Get numeric arguments.
   if (nargin > 1)
     {
       if (! args(1).isempty ())
@@ -115,24 +131,34 @@
         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 fnda = arg.float_array_value ();
+          fnda.resize (dims, 0.0);
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
+          if (ifft_symmetric)
+            retval = fnda.ifouriersym (dim);
+          else
+            retval = (type != 0 ? fnda.ifourier (dim) : fnda.fourier (dim));
         }
       else
         {
-          FloatComplexNDArray cnda = arg.float_complex_array_value ();
+          FloatComplexNDArray fcnda = arg.float_complex_array_value ();
+          fcnda.resize (dims, 0.0);
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
+          if (ifft_symmetric)
+            retval = fcnda.ifouriersym (dim);
+          else
+            retval = (type != 0 ? fcnda.ifourier (dim) : fcnda.fourier (dim));
         }
     }
   else
@@ -140,16 +166,22 @@
       if (arg.isreal ())
         {
           NDArray nda = arg.array_value ();
+          nda.resize (dims, 0.0);
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
+          if (ifft_symmetric)
+            retval = nda.ifouriersym (dim);
+          else
+            retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
         }
       else if (arg.iscomplex ())
         {
           ComplexNDArray cnda = arg.complex_array_value ();
+          cnda.resize (dims, 0.0);
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
+          if (ifft_symmetric)
+            retval = cnda.ifouriersym (dim);
+          else
+            retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
         }
       else
         err_wrong_type_arg (fcn, arg);
@@ -227,6 +259,7 @@
 @deftypefn  {} {@var{x} =} ifft (@var{y})
 @deftypefnx {} {@var{x} =} ifft (@var{y}, @var{n})
 @deftypefnx {} {@var{x} =} ifft (@var{y}, @var{n}, @var{dim})
+@deftypefnx {} {@var{x} =} ifft (@dots{}, @var{symflag})
 Compute the inverse discrete Fourier transform of @var{y}
 using a Fast Fourier Transform (FFT) algorithm.
 
@@ -244,6 +277,12 @@
 
 If called with three arguments, @var{dim} is an integer specifying the
 dimension of the matrix along which the inverse FFT is performed.
+
+If @var{y} is conjugate-symmetric, the inverse Fourier transform @var{x} is
+real-valued.  In this case, the keyword @code{@var{symmetric}} can be specified
+to ensure a real-valued output.  Asymmetries in @var{y} are eliminated by
+using only the numbers associated with the positive frequencies for the
+transformation.
 @seealso{fft, ifft2, ifftn, fftw}
 @end deftypefn */)
 {
@@ -251,6 +290,129 @@
 }
 
 /*
+## Input validation.
+
+%!testif HAVE_FFTW
+%!error <unknown symflag> ifft (1, "unknown")
+
+%!testif HAVE_FFTW
+%! assert (ifft (1i), 1i, 0);
+%! assert (ifft (1i, "nonsymmetric"), 1i, 0);
+%! assert (ifft (1i, "symmetric"), 0, 0);
+
+## Padding with zero.
+%!testif HAVE_FFTW
+%! assert (ifft ([2, 0, 2], 4), [1, 0, 1, 0], 0);
+%! assert (ifft ([2, 0, 2], 4, "symmetric"), [1, 0, 1, 0], 0);
+
+## Truncate. Cases N=1 and N=2 are handled by different code.
+%!testif HAVE_FFTW
+%! y = [10+6i, 12+2i; 12-2i, -10+14i; -11+4i, 4-4i];
+%! assert (ifft (y, 1), [10+6i, 12+2i], 0);
+%! assert (ifft (y, 1, "symmetric"), [10, 12], 0);
+%! assert (ifft (y, 2), [11+2i, 1+8i; -1+4i, 11-6i], 0);
+%! assert (ifft (y, 2, "symmetric"), [11, 1; -1, 11], 0);
+
+## Double precision data types.
+
+## dNDArray.
+%!testif HAVE_FFTW
+%! y = [20, -8, -40, -8];
+%! x1 = ifft (y);
+%! x2 = ifft (y, "symmetric");
+%! assert (x1, [-9, 15, -1, 15], 0);
+%! assert (x2, [-9, 15, -1, 15], 0);
+
+## CNDArray with conjugate-symmetric input.
+%!testif HAVE_FFTW
+%! y = [4, 2+2i, 12, 2-2i];
+%! x1 = ifft (y);
+%! x2 = ifft (y, "symmetric");
+%! assert (x1, [5, -3, 3, -1], 0);
+%! assert (x2, [5, -3, 3, -1], 0);
+
+## CNDArray non-symmetric input.
+%!testif HAVE_FFTW
+%! y = [20+92i, -18+12i, 28+100i, 10+60i];
+%! x1 = ifft (y);
+%! x2 = ifft (y, "symmetric");
+%! assert (x1, [10+66i, 10-9i, 14+30i, -14+5i], 0);
+%! assert (x2, [3, -8, 21, 4], 0);
+
+## dim=0; stride=1; howmany=8; dist=3; n=3; nloop=1.
+%!testif HAVE_FFTW
+%! x = reshape (mod (1:24, 5), 3, 4, 2);
+%! y = fft (x, [], 1);
+%! x1 = ifft (y, [], 1);
+%! x2 = ifft (y, [], 1, "symmetric");
+%! assert (x1, x, 2*eps);
+%! assert (x2, x, 2*eps);
+
+## dim=1; stride=3; howmany=3; dist=1; n=4; nloop=2.
+%!testif HAVE_FFTW
+%! x = reshape (mod (1:24, 5), 3, 4, 2);
+%! y = fft (x, [], 2);
+%! x1 = ifft (y, [], 2);
+%! x2 = ifft (y, [], 2, "symmetric");
+%! assert (x1, x, 0);
+%! assert (x2, x, 0);
+
+## dim=2; stride=12; howmany=12; dist=1; n=2; nloop=1.
+%!testif HAVE_FFTW
+%! x = reshape (mod (1:24, 5), 3, 4, 2);
+%! y = fft (x, [], 3);
+%! x1 = ifft (y, [], 3);
+%! x2 = ifft (y, [], 3, "symmetric");
+%! assert (x1, x, 0);
+%! assert (x2, x, 0);
+
+## Single precision data types.
+
+%!testif HAVE_FFTW
+%! y = single ([20, -8, -40, -8]);
+%! x1 = ifft (y);
+%! x2 = ifft (y, "symmetric");
+%! assert (x1, single ([-9, 15, -1, 15]), 0);
+%! assert (x2, single ([-9, 15, -1, 15]), 0);
+
+%!testif HAVE_FFTW
+%! y = single ([4, 2+2i, 12, 2-2i]);
+%! x1 = ifft (y);
+%! x2 = ifft (y, "symmetric");
+%! assert (x1, single ([5, -3, 3, -1]), 0);
+%! assert (x2, single ([5, -3, 3, -1]), 0);
+
+%!testif HAVE_FFTW
+%! y = single ([20+92i, -18+12i, 28+100i, 10+60i]);
+%! x1 = ifft (y);
+%! x2 = ifft (y, "symmetric");
+%! assert (x1, single ([10+66i, 10-9i, 14+30i, -14+5i]), 0);
+%! assert (x2, single ([3, -8, 21, 4]), 0);
+
+%!testif HAVE_FFTW
+%! x = single (reshape (mod (1:24, 5), 3, 4, 2));
+%! y = fft (x, [], 1);
+%! x1 = ifft (y, [], 1);
+%! x2 = ifft (y, [], 1, "symmetric");
+%! assert (x1, x, 2*eps ("single"));
+%! assert (x2, x, 2*eps ("single"));
+
+%!testif HAVE_FFTW
+%! x = single (reshape (mod (1:24, 5), 3, 4, 2));
+%! y = fft (x, [], 2);
+%! x1 = ifft (y, [], 2);
+%! x2 = ifft (y, [], 2, "symmetric");
+%! assert (x1, x, 0);
+%! assert (x2, x, 0);
+
+%!testif HAVE_FFTW
+%! x = single (reshape (mod (1:24, 5), 3, 4, 2));
+%! y = fft (x, [], 3);
+%! x1 = ifft (y, [], 3);
+%! x2 = ifft (y, [], 3, "symmetric");
+%! assert (x1, x, 0);
+%! assert (x2, x, 0);
+
 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
 ##         Comalco Research and Technology
 ##         02 May 2000
diff -r 399be7cc310f -r 342e2e9bff3d libinterp/corefcn/fft2.cc
--- a/libinterp/corefcn/fft2.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/libinterp/corefcn/fft2.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -44,14 +44,30 @@
 {
   int nargin = args.length ();
 
-  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;
 
+  // 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);
+      nargin--;
+    }
+  else if (nargin == 4)
+    print_usage ();
+
+  // Get numeric arguments.
   if (nargin > 1)
     {
       double dval = args(1).double_value ();
@@ -101,17 +117,23 @@
     {
       if (arg.isreal ())
         {
-          FloatNDArray nda = arg.float_array_value ();
+          FloatNDArray fnda = arg.float_array_value ();
+          fnda.resize (dims, 0.0);
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
+          if (ifft2_symmetric)
+            retval = fnda.ifourier2dsym ();
+          else
+            retval = (type != 0 ? fnda.ifourier2d () : fnda.fourier2d ());
         }
       else
         {
-          FloatComplexNDArray cnda = arg.float_complex_array_value ();
+          FloatComplexNDArray fcnda = arg.float_complex_array_value ();
+          fcnda.resize (dims, 0.0);
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
+          if (ifft2_symmetric)
+            retval = fcnda.ifourier2dsym ();
+          else
+            retval = (type != 0 ? fcnda.ifourier2d () : fcnda.fourier2d ());
         }
     }
   else
@@ -119,16 +141,22 @@
       if (arg.isreal ())
         {
           NDArray nda = arg.array_value ();
+          nda.resize (dims, 0.0);
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
+          if (ifft2_symmetric)
+            retval = nda.ifourier2dsym ();
+          else
+            retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
         }
       else if (arg.iscomplex ())
         {
           ComplexNDArray cnda = arg.complex_array_value ();
+          cnda.resize (dims, 0.0);
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
+          if (ifft2_symmetric)
+            retval = cnda.ifourier2dsym ();
+          else
+            retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
         }
       else
         err_wrong_type_arg (fcn, arg);
@@ -175,6 +203,7 @@
        doc: /* -*- texinfo -*-
 @deftypefn  {} {@var{A} =} ifft2 (@var{B})
 @deftypefnx {} {@var{A} =} ifft2 (@var{B}, @var{m}, @var{n})
+@deftypefnx {} {@var{A} =} ifft2 (@dots{}, @var{symflag})
 Compute the inverse two-dimensional discrete Fourier transform of @var{B}
 using a Fast Fourier Transform (FFT) algorithm.
 
@@ -184,6 +213,11 @@
 
 If @var{B} is a multi-dimensional matrix, each two-dimensional sub-matrix
 of @var{B} is treated separately.
+
+If @var{B} is conjugate-symmetric, the inverse Fourier transform @var{A} is
+real-valued.  In this case, the keyword @code{@var{symmetric}} can be specified
+to ensure a real-valued output.  Asymmetries in @var{B} are eliminated by using
+a non-redundant part of the spectrum including the positive frequency locations.
 @seealso{fft2, ifft, ifftn, fftw}
 @end deftypefn */)
 {
@@ -191,6 +225,16 @@
 }
 
 /*
+## Input validation.
+
+%!testif HAVE_FFTW
+%!error <unknown symflag> ifft2 (1, "unknown")
+
+%!testif HAVE_FFTW
+%! assert (ifft2 (1i), 1i, 0);
+%! assert (ifft2 (1i, "nonsymmetric"), 1i, 0);
+%! assert (ifft2 (1i, "symmetric"), 0, 0);
+
 %!testif HAVE_FFTW <*65414>
 %! sz = size (ifft2 (ones (2, 0, 3)));
 %! assert (sz, [2, 0, 3]);
@@ -199,6 +243,92 @@
 %! sz = size (ifft2 (ones (5, 4, 3), 2, 0));
 %! assert (sz, [2, 0, 3]);
 
+## Truncate rows and columns.
+%!testif HAVE_FFTW
+%! B = [-1,-7-4i,35,-7+4i; -3-10i,23-24i,-13,-11-14i; -33,7-6i,39,7+6i; -3+10i,-11+14i,-13,23+24i];
+%! assert (ifft2 (B, 3, 2), ifft2 (B(1:3, 1:2)), 0);
+%! assert (ifft2 (B, 3, 2, "symmetric"), ifft2 (B(1:3, 1:2), "symmetric"), 0);
+
+## Double precision data types.
+
+## dNDArray.
+%!testif HAVE_FFTW
+%! A = reshape (mod (1:24, 5), 3, 2, 4);
+%! B = fft2 (A);
+%! A1 = ifft2 (B);
+%! A2 = ifft2 (B, "symmetric");
+%! assert (A1, A, 10*eps);
+%! assert (A2, A, 10*eps);
+
+## CNDArray with conjugate-symmetric input.
+%!testif HAVE_FFTW
+%! A = [1,1,1,1,1; 16,8,4,2,1; 81,27,9,3,1; 256,64,16,4,1];
+%! B = fft2 (A);
+%! A1 = ifft2 (B);
+%! A2 = ifft2 (B, "symmetric");
+%! assert (A1, A, 60*eps);
+%! assert (A2, A, 60*eps);
+
+## Check "makesymmetric" for non-symmetric input.
+## howmany=1; dist=16.
+%!testif HAVE_FFTW
+%! A1 = [34+10i,-16i,-8-8i,-16; -4i,0,0,0; -2-2i,0,0,0; -4,0,0,0];
+%! A2 = [25,-9,-1,-7; -1,6,-1,2; -1,1,1,-1; -3,-4,1,-4];
+%! B = 4 * (reshape (1:16, 4, 4) + 1i * reshape ((1:16) - 6, 4, 4));
+%! A3 = ifft2 (B);
+%! A4 = ifft2 (B, "symmetric");
+%! assert (A3, A1, 0);
+%! assert (A4, A2, 0);
+
+## Check "makesymmetric" for multivariate, non-symmetric input: The ifft2 of
+## individual slices in a multivariate data set must be equal to the ifft2 over
+## the entire data set.
+## howmany=2; dist=16.
+%!testif HAVE_FFTW
+%! B(:, :, 1) = reshape (1:16, 4, 4) + 1i * reshape ((1:16) - 5, 4, 4);
+%! B(:, :, 2) = reshape (17:32, 4, 4) + 1i * reshape ((17:32) - 5, 4, 4);
+%! A1 = ifft2 (B, "symmetric");
+%! A2 = ifft2 (B(:, :, 1), "symmetric");
+%! A3 = ifft2 (B(:, :, 2), "symmetric");
+%! assert (A2, A1(:, :, 1));
+%! assert (A3, A1(:, :, 2));
+
+## Single precision data types.
+
+%!testif HAVE_FFTW
+%! A = single (reshape (mod (1:24, 5), 3, 2, 4));
+%! B = fft2 (A);
+%! A1 = ifft2 (B);
+%! A2 = ifft2 (B, "symmetric");
+%! assert (A1, A, 10*eps ("single"));
+%! assert (A2, A, 10*eps ("single"));
+
+%!testif HAVE_FFTW
+%! A = single ([1,1,1,1,1; 16,8,4,2,1; 81,27,9,3,1; 256,64,16,4,1]);
+%! B = fft2 (A);
+%! A1 = ifft2 (B);
+%! A2 = ifft2 (B, "symmetric");
+%! assert (A1, A, 130*eps ("single"));
+%! assert (A2, A, 130*eps ("single"));
+
+%!testif HAVE_FFTW
+%! A1 = single ([34+10i,-16i,-8-8i,-16; -4i,0,0,0; -2-2i,0,0,0; -4,0,0,0]);
+%! A2 = single ([25,-9,-1,-7; -1,6,-1,2; -1,1,1,-1; -3,-4,1,-4]);
+%! B = single (4 * (reshape (1:16, 4, 4) + 1i * reshape ((1:16) - 6, 4, 4)));
+%! A3 = ifft2 (B);
+%! A4 = ifft2 (B, "symmetric");
+%! assert (A3, A1, 0);
+%! assert (A4, A2, 0);
+
+%!testif HAVE_FFTW
+%! B(:, :, 1) = single (reshape (1:16, 4, 4) + 1i * reshape ((1:16) - 5, 4, 4));
+%! B(:, :, 2) = single (reshape (17:32, 4, 4) + 1i * reshape ((17:32) - 5, 4, 4));
+%! A1 = ifft2 (B, "symmetric");
+%! A2 = ifft2 (B(:, :, 1), "symmetric");
+%! A3 = ifft2 (B(:, :, 2), "symmetric");
+%! assert (A2, A1(:, :, 1));
+%! assert (A3, A1(:, :, 2));
+
 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
 ##         Comalco Research and Technology
 ##         02 May 2000
diff -r 399be7cc310f -r 342e2e9bff3d libinterp/corefcn/fftn.cc
--- a/libinterp/corefcn/fftn.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/libinterp/corefcn/fftn.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -44,17 +44,33 @@
 {
   int nargin = args.length ();
 
-  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;
 
+  // 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);
+      nargin--;
+    }
+  else if (nargin == 3)
+    print_usage ();
+
+  // Get numeric arguments.
   if (nargin > 1)
     {
       Matrix val = args(1).xmatrix_value ("%s: SIZE must be a vector of length dim", fcn);
@@ -88,17 +104,23 @@
     {
       if (arg.isreal ())
         {
-          FloatNDArray nda = arg.float_array_value ();
+          FloatNDArray fnda = arg.float_array_value ();
+          fnda.resize (dims, 0.0);
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
+          if (ifftn_symmetric)
+            retval = fnda.ifourierNdsym ();
+          else
+            retval = (type != 0 ? fnda.ifourierNd () : fnda.fourierNd ());
         }
       else
         {
-          FloatComplexNDArray cnda = arg.float_complex_array_value ();
+          FloatComplexNDArray fcnda = arg.float_complex_array_value ();
+          fcnda.resize (dims, 0.0);
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
+          if (ifftn_symmetric)
+            retval = fcnda.ifourierNdsym ();
+          else
+            retval = (type != 0 ? fcnda.ifourierNd () : fcnda.fourierNd ());
         }
     }
   else
@@ -106,16 +128,22 @@
       if (arg.isreal ())
         {
           NDArray nda = arg.array_value ();
+          nda.resize (dims, 0.0);
 
-          nda.resize (dims, 0.0);
-          retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
+          if (ifftn_symmetric)
+            retval = nda.ifourierNdsym ();
+          else
+            retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
         }
       else if (arg.iscomplex ())
         {
           ComplexNDArray cnda = arg.complex_array_value ();
+          cnda.resize (dims, 0.0);
 
-          cnda.resize (dims, 0.0);
-          retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
+          if (ifftn_symmetric)
+            retval = cnda.ifourierNdsym ();
+          else
+            retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
         }
       else
         err_wrong_type_arg (fcn, arg);
@@ -157,6 +185,7 @@
        doc: /* -*- texinfo -*-
 @deftypefn  {} {@var{A} =} ifftn (@var{B})
 @deftypefnx {} {@var{A} =} ifftn (@var{B}, @var{size})
+@deftypefnx {} {@var{A} =} ifftn (@dots{}, @var{symflag})
 Compute the inverse N-dimensional discrete Fourier transform of @var{B}
 using a Fast Fourier Transform (FFT) algorithm.
 
@@ -166,6 +195,11 @@
 truncated prior to performing the inverse FFT@.  Otherwise, if an element of
 @var{size} is larger than the corresponding dimension then @var{B} is
 resized and padded with zeros.
+
+If @var{B} is conjugate-symmetric, the inverse Fourier transform @var{A} is
+real-valued.  In this case, the keyword @code{@var{symmetric}} can be specified
+to ensure a real-valued output.  Asymmetries in @var{B} are eliminated by using
+a non-redundant part of the spectrum including the positive frequency locations.
 @seealso{fftn, ifft, ifft2, fftw}
 @end deftypefn */)
 {
@@ -173,6 +207,16 @@
 }
 
 /*
+## Input validation.
+
+%!testif HAVE_FFTW
+%!error <unknown symflag> ifftn (1, "unknown")
+
+%!testif HAVE_FFTW
+%! assert (ifftn (1i), 1i, 0);
+%! assert (ifftn (1i, "nonsymmetric"), 1i, 0);
+%! assert (ifftn (1i, "symmetric"), 0, 0);
+
 %!testif HAVE_FFTW <*65414>
 %! sz = size (ifftn (ones (2, 0, 3)));
 %! assert (sz, [2, 0, 3]);
@@ -180,6 +224,62 @@
 %!testif HAVE_FFTW <*65414>
 %! sz = size (ifftn (ones (5, 4, 3), [2, 0, 3]));
 %! assert (sz, [2, 0, 3]);
+
+## Double precision data types.
+
+## dNDArray.
+%!testif HAVE_FFTW
+%! A(:, :, 1) = [5, 8; 6, 9];
+%! A(:, :, 2) = [5, 7; 1, 5];
+%! B(:, :, 1) = [46, -12; 4, 2];
+%! B(:, :, 2) = [10, 0; -8, -2];
+%! A1 = ifftn (B);
+%! A2 = ifftn (B, "symmetric");
+%! assert (A1, A, 0);
+%! assert (A2, A, 0);
+
+## Check "makesymmetric" for non-symmetric input.
+%!testif HAVE_FFTW
+%! A1(:, :, 1) = [65+25i,-8i,-4-4i,-8; -2i,0,0,0; -1-1i,0,0,0; -2,0,0,0];
+%! A1(:, :, 2) = [-32i,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0];
+%! A1(:, :, 3) = [-16-16i,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0];
+%! A1(:, :, 4) = [-32,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0];
+%! A2(:, :, 1) = [187,-13,-9,-19; -3,6,-3,2; -3,1,1,-1; -5,-4,3,-4];
+%! A2(:, :, 2) = [-82,48,8,-32; 20,0,16,0; 2,0,0,0; -16,0,-16,0];
+%! A2(:, :, 3) = [-9,-13,7,13; -3,6,-3,2; 1,1,1,-1; 3,-4,3,-4];
+%! A2(:, :, 4) = [-46,16,-8,-32; 12,0,16,0; -2,0,0,0; -16,0,-16,0];
+%! B = reshape (1:64, 4, 4, 4) + 1i * reshape ((1:64) - 20, 4, 4, 4);
+%! A3 = ifftn (B);
+%! A4 = ifftn (B, "symmetric");
+%! assert (2 * A3, A1, 0);
+%! assert (8 * A4, A2, 0);
+
+## Single precision data types.
+
+%!testif HAVE_FFTW
+%! A(:, :, 1) = single ([5, 8; 6, 9]);
+%! A(:, :, 2) = single ([5, 7; 1, 5]);
+%! B(:, :, 1) = single ([46, -12; 4, 2]);
+%! B(:, :, 2) = single ([10, 0; -8, -2]);
+%! A1 = ifftn (B);
+%! A2 = ifftn (B, "symmetric");
+%! assert (A1, A, 0);
+%! assert (A2, A, 0);
+
+%!testif HAVE_FFTW
+%! A1(:, :, 1) = single ([65+25i,-8i,-4-4i,-8; -2i,0,0,0; -1-1i,0,0,0; -2,0,0,0]);
+%! A1(:, :, 2) = single ([-32i,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0]);
+%! A1(:, :, 3) = single ([-16-16i,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0]);
+%! A1(:, :, 4) = single ([-32,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0]);
+%! A2(:, :, 1) = single ([187,-13,-9,-19; -3,6,-3,2; -3,1,1,-1; -5,-4,3,-4]);
+%! A2(:, :, 2) = single ([-82,48,8,-32; 20,0,16,0; 2,0,0,0; -16,0,-16,0]);
+%! A2(:, :, 3) = single ([-9,-13,7,13; -3,6,-3,2; 1,1,1,-1; 3,-4,3,-4]);
+%! A2(:, :, 4) = single ([-46,16,-8,-32; 12,0,16,0; -2,0,0,0; -16,0,-16,0]);
+%! B = single (reshape (1:64, 4, 4, 4) + 1i * reshape ((1:64) - 20, 4, 4, 4));
+%! A3 = ifftn (B);
+%! A4 = ifftn (B, "symmetric");
+%! assert (2 * A3, A1, 0);
+%! assert (8 * A4, A2, 0);
 */
 
 OCTAVE_END_NAMESPACE(octave)
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/CNDArray.cc
--- a/liboctave/array/CNDArray.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/CNDArray.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -116,6 +116,37 @@
   return retval;
 }
 
+NDArray
+ComplexNDArray::ifouriersym (int dim) const
+{
+  const dim_vector& dv = dims ();
+
+  if (dim > dv.ndims () || dim < 0)
+    return NDArray ();
+
+  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 () / n;
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
+  octave_idx_type dist = (stride == 1 ? n : 1);
+
+  const Complex *in (data ());
+  NDArray retval (dv);
+  double *out (retval.rwdata ());
+
+  // Need to be careful here about the distance between fft's
+  for (octave_idx_type k = 0; k < nloop; k++)
+    octave::fftw::ifftsym (in + k * stride * n, out + k * stride * n,
+                           n, howmany, stride, dist);
+
+  return retval;
+}
+
 ComplexNDArray
 ComplexNDArray::fourier2d () const
 {
@@ -156,6 +187,26 @@
   return retval;
 }
 
+NDArray
+ComplexNDArray::ifourier2dsym () const
+{
+  const dim_vector& dv = dims ();
+  if (dv.ndims () < 2)
+    return NDArray ();
+
+  dim_vector dv2 (dv(0), dv(1));
+  const Complex *in = data ();
+  NDArray retval (dv);
+  double *out = retval.rwdata ();
+  octave_idx_type howmany = numel () / dv(0) / dv(1);
+  octave_idx_type dist = dv(0) * dv(1);
+
+  for (octave_idx_type i=0; i < howmany; i++)
+    octave::fftw::ifftNdsym (in + i*dist, out + i*dist, 2, dv2);
+
+  return retval;
+}
+
 ComplexNDArray
 ComplexNDArray::fourierNd () const
 {
@@ -186,6 +237,21 @@
   return retval;
 }
 
+NDArray
+ComplexNDArray::ifourierNdsym () const
+{
+  const dim_vector& dv = dims ();
+  int rank = dv.ndims ();
+
+  const Complex *in (data ());
+  NDArray retval (dv);
+  double *out (retval.rwdata ());
+
+  octave::fftw::ifftNdsym (in, out, rank, dv);
+
+  return retval;
+}
+
 #else
 
 ComplexNDArray
@@ -210,6 +276,17 @@
   return ComplexNDArray ();
 }
 
+NDArray
+ComplexNDArray::ifouriersym (int dim) const
+{
+  octave_unused_parameter (dim);
+
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return NDArray ();
+}
+
 ComplexNDArray
 ComplexNDArray::fourier2d () const
 {
@@ -228,6 +305,15 @@
   return ComplexNDArray ();
 }
 
+NDArray
+ComplexNDArray::ifourier2dsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return NDArray ();
+}
+
 ComplexNDArray
 ComplexNDArray::fourierNd () const
 {
@@ -246,6 +332,15 @@
   return ComplexNDArray ();
 }
 
+NDArray
+ComplexNDArray::ifourierNdsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return NDArray ();
+}
+
 #endif
 
 // unary operations
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/CNDArray.h
--- a/liboctave/array/CNDArray.h	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/CNDArray.h	Sun Mar 24 19:34:28 2024 +0100
@@ -121,12 +121,15 @@
 
   OCTAVE_API ComplexNDArray fourier (int dim = 1) const;
   OCTAVE_API ComplexNDArray ifourier (int dim = 1) const;
+  OCTAVE_API NDArray ifouriersym (int dim = 1) const;
 
   OCTAVE_API ComplexNDArray fourier2d () const;
   OCTAVE_API ComplexNDArray ifourier2d () const;
+  OCTAVE_API NDArray ifourier2dsym () const;
 
   OCTAVE_API ComplexNDArray fourierNd () const;
   OCTAVE_API ComplexNDArray ifourierNd () const;
+  OCTAVE_API NDArray ifourierNdsym () const;
 
   ComplexNDArray squeeze () const { return MArray<Complex>::squeeze (); }
 
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/dNDArray.cc
--- a/liboctave/array/dNDArray.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/dNDArray.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -157,6 +157,38 @@
   return retval;
 }
 
+NDArray
+NDArray::ifouriersym (int dim) const
+{
+  const dim_vector& dv = dims ();
+
+  if (dim > dv.ndims () || dim < 0)
+    return NDArray ();
+
+  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 () / n;
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
+  octave_idx_type dist = (stride == 1 ? n : 1);
+
+  ComplexNDArray tmp (*this);
+  Complex *in (tmp.rwdata ());
+  NDArray retval (dv);
+  double *out (retval.rwdata ());
+
+  // Need to be careful here about the distance between fft's
+  for (octave_idx_type k = 0; k < nloop; k++)
+    octave::fftw::ifftsym (in + k * stride * n, out + k * stride * n,
+                           n, howmany, stride, dist);
+
+  return retval;
+}
+
 ComplexNDArray
 NDArray::fourier2d () const
 {
@@ -196,6 +228,27 @@
   return retval;
 }
 
+NDArray
+NDArray::ifourier2dsym () const
+{
+  const dim_vector& dv = dims ();
+  if (dv.ndims () < 2)
+    return NDArray ();
+
+  dim_vector dv2 (dv(0), dv(1));
+  ComplexNDArray tmp (*this);
+  Complex *in (tmp.rwdata ());
+  NDArray retval (dv);
+  double *out (retval.rwdata ());
+  octave_idx_type howmany = numel () / dv(0) / dv(1);
+  octave_idx_type dist = dv(0) * dv(1);
+
+  for (octave_idx_type i=0; i < howmany; i++)
+    octave::fftw::ifftNdsym (in + i*dist, out + i*dist, 2, dv2);
+
+  return retval;
+}
+
 ComplexNDArray
 NDArray::fourierNd () const
 {
@@ -227,6 +280,22 @@
   return retval;
 }
 
+NDArray
+NDArray::ifourierNdsym () const
+{
+  const dim_vector& dv = dims ();
+  int rank = dv.ndims ();
+
+  ComplexNDArray tmp (*this);
+  Complex *in (tmp.rwdata ());
+  NDArray retval (dv);
+  double *out (retval.rwdata ());
+
+  octave::fftw::ifftNdsym (in, out, rank, dv);
+
+  return retval;
+}
+
 #else
 
 ComplexNDArray
@@ -251,6 +320,17 @@
   return ComplexNDArray ();
 }
 
+NDArray
+NDArray::ifouriersym (int dim) const
+{
+  octave_unused_parameter (dim);
+
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return NDArray ();
+}
+
 ComplexNDArray
 NDArray::fourier2d () const
 {
@@ -269,6 +349,15 @@
   return ComplexNDArray ();
 }
 
+NDArray
+NDArray::ifourier2dsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return NDArray ();
+}
+
 ComplexNDArray
 NDArray::fourierNd () const
 {
@@ -287,6 +376,15 @@
   return ComplexNDArray ();
 }
 
+NDArray
+NDArray::ifourierNdsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return NDArray ();
+}
+
 #endif
 
 // unary operations
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/dNDArray.h
--- a/liboctave/array/dNDArray.h	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/dNDArray.h	Sun Mar 24 19:34:28 2024 +0100
@@ -128,12 +128,15 @@
 
   OCTAVE_API ComplexNDArray fourier (int dim = 1) const;
   OCTAVE_API ComplexNDArray ifourier (int dim = 1) const;
+  OCTAVE_API NDArray ifouriersym (int dim = 1) const;
 
   OCTAVE_API ComplexNDArray fourier2d () const;
   OCTAVE_API ComplexNDArray ifourier2d () const;
+  OCTAVE_API NDArray ifourier2dsym () const;
 
   OCTAVE_API ComplexNDArray fourierNd () const;
   OCTAVE_API ComplexNDArray ifourierNd () const;
+  OCTAVE_API NDArray ifourierNdsym () const;
 
   friend OCTAVE_API NDArray real (const ComplexNDArray& a);
   friend OCTAVE_API NDArray imag (const ComplexNDArray& a);
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/fCNDArray.cc
--- a/liboctave/array/fCNDArray.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/fCNDArray.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -116,6 +116,37 @@
   return retval;
 }
 
+FloatNDArray
+FloatComplexNDArray::ifouriersym (int dim) const
+{
+  const dim_vector& dv = dims ();
+
+  if (dim > dv.ndims () || dim < 0)
+    return FloatNDArray ();
+
+  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 () / n;
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
+  octave_idx_type dist = (stride == 1 ? n : 1);
+
+  const FloatComplex *in (data ());
+  FloatNDArray retval (dv);
+  float *out (retval.rwdata ());
+
+  // Need to be careful here about the distance between fft's
+  for (octave_idx_type k = 0; k < nloop; k++)
+    octave::fftw::ifftsym (in + k * stride * n, out + k * stride * n,
+                           n, howmany, stride, dist);
+
+  return retval;
+}
+
 FloatComplexNDArray
 FloatComplexNDArray::fourier2d () const
 {
@@ -156,6 +187,26 @@
   return retval;
 }
 
+FloatNDArray
+FloatComplexNDArray::ifourier2dsym () const
+{
+  const dim_vector& dv = dims ();
+  if (dv.ndims () < 2)
+    return FloatNDArray ();
+
+  dim_vector dv2 (dv(0), dv(1));
+  const FloatComplex *in = data ();
+  FloatNDArray retval (dv);
+  float *out = retval.rwdata ();
+  octave_idx_type howmany = numel () / dv(0) / dv(1);
+  octave_idx_type dist = dv(0) * dv(1);
+
+  for (octave_idx_type i=0; i < howmany; i++)
+    octave::fftw::ifftNdsym (in + i*dist, out + i*dist, 2, dv2);
+
+  return retval;
+}
+
 FloatComplexNDArray
 FloatComplexNDArray::fourierNd () const
 {
@@ -186,6 +237,21 @@
   return retval;
 }
 
+FloatNDArray
+FloatComplexNDArray::ifourierNdsym () const
+{
+  const dim_vector& dv = dims ();
+  int rank = dv.ndims ();
+
+  const FloatComplex *in (data ());
+  FloatNDArray retval (dv);
+  float *out (retval.rwdata ());
+
+  octave::fftw::ifftNdsym (in, out, rank, dv);
+
+  return retval;
+}
+
 #else
 
 FloatComplexNDArray
@@ -210,6 +276,17 @@
   return FloatComplexNDArray ();
 }
 
+FloatNDArray
+FloatComplexNDArray::ifouriersym (int dim) const
+{
+  octave_unused_parameter (dim);
+
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return FloatNDArray ();
+}
+
 FloatComplexNDArray
 FloatComplexNDArray::fourier2d () const
 {
@@ -228,6 +305,15 @@
   return FloatComplexNDArray ();
 }
 
+FloatNDArray
+FloatComplexNDArray::ifourier2dsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return FloatNDArray ();
+}
+
 FloatComplexNDArray
 FloatComplexNDArray::fourierNd () const
 {
@@ -246,6 +332,15 @@
   return FloatComplexNDArray ();
 }
 
+FloatNDArray
+FloatComplexNDArray::ifourierNdsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return FloatNDArray ();
+}
+
 #endif
 
 // unary operations
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/fCNDArray.h
--- a/liboctave/array/fCNDArray.h	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/fCNDArray.h	Sun Mar 24 19:34:28 2024 +0100
@@ -123,12 +123,15 @@
 
   OCTAVE_API FloatComplexNDArray fourier (int dim = 1) const;
   OCTAVE_API FloatComplexNDArray ifourier (int dim = 1) const;
+  OCTAVE_API FloatNDArray ifouriersym (int dim = 1) const;
 
   OCTAVE_API FloatComplexNDArray fourier2d () const;
   OCTAVE_API FloatComplexNDArray ifourier2d () const;
+  OCTAVE_API FloatNDArray ifourier2dsym () const;
 
   OCTAVE_API FloatComplexNDArray fourierNd () const;
   OCTAVE_API FloatComplexNDArray ifourierNd () const;
+  OCTAVE_API FloatNDArray ifourierNdsym () const;
 
   FloatComplexNDArray squeeze () const
   { return MArray<FloatComplex>::squeeze (); }
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/fNDArray.cc
--- a/liboctave/array/fNDArray.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/fNDArray.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -115,6 +115,38 @@
   return retval;
 }
 
+FloatNDArray
+FloatNDArray::ifouriersym (int dim) const
+{
+  const dim_vector& dv = dims ();
+
+  if (dim > dv.ndims () || dim < 0)
+    return FloatNDArray ();
+
+  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 () / n;
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  octave_idx_type nloop = (stride == 1 ? 1 : numel () / n / stride);
+  octave_idx_type dist = (stride == 1 ? n : 1);
+
+  FloatComplexNDArray tmp (*this);
+  FloatComplex *in (tmp.rwdata ());
+  FloatNDArray retval (dv);
+  float *out (retval.rwdata ());
+
+  // Need to be careful here about the distance between fft's
+  for (octave_idx_type k = 0; k < nloop; k++)
+    octave::fftw::ifftsym (in + k * stride * n, out + k * stride * n,
+                           n, howmany, stride, dist);
+
+  return retval;
+}
+
 FloatComplexNDArray
 FloatNDArray::fourier2d () const
 {
@@ -154,6 +186,27 @@
   return retval;
 }
 
+FloatNDArray
+FloatNDArray::ifourier2dsym () const
+{
+  const dim_vector& dv = dims ();
+  if (dv.ndims () < 2)
+    return FloatNDArray ();
+
+  dim_vector dv2 (dv(0), dv(1));
+  FloatComplexNDArray tmp (*this);
+  FloatComplex *in (tmp.rwdata ());
+  FloatNDArray retval (dv);
+  float *out (retval.rwdata ());
+  octave_idx_type howmany = numel () / dv(0) / dv(1);
+  octave_idx_type dist = dv(0) * dv(1);
+
+  for (octave_idx_type i=0; i < howmany; i++)
+    octave::fftw::ifftNdsym (in + i*dist, out + i*dist, 2, dv2);
+
+  return retval;
+}
+
 FloatComplexNDArray
 FloatNDArray::fourierNd () const
 {
@@ -185,6 +238,22 @@
   return retval;
 }
 
+FloatNDArray
+FloatNDArray::ifourierNdsym () const
+{
+  const dim_vector& dv = dims ();
+  int rank = dv.ndims ();
+
+  FloatComplexNDArray tmp (*this);
+  FloatComplex *in (tmp.rwdata ());
+  FloatNDArray retval (dv);
+  float *out (retval.rwdata ());
+
+  octave::fftw::ifftNdsym (in, out, rank, dv);
+
+  return retval;
+}
+
 #else
 
 FloatComplexNDArray
@@ -208,6 +277,16 @@
   return FloatComplexNDArray ();
 }
 
+FloatNDArray
+FloatNDArray::ifouriersym (int dim) const
+{
+  octave_unused_parameter (dim);
+
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+  return FloatNDArray ();
+}
+
 FloatComplexNDArray
 FloatNDArray::fourier2d () const
 {
@@ -226,6 +305,15 @@
   return FloatComplexNDArray ();
 }
 
+FloatNDArray
+FloatNDArray::ifourier2dsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return FloatNDArray ();
+}
+
 FloatComplexNDArray
 FloatNDArray::fourierNd () const
 {
@@ -244,6 +332,15 @@
   return FloatComplexNDArray ();
 }
 
+FloatNDArray
+FloatNDArray::ifourierNdsym () const
+{
+  (*current_liboctave_error_handler)
+    ("support for FFTW was unavailable or disabled when liboctave was built");
+
+  return FloatNDArray ();
+}
+
 #endif
 
 // unary operations
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/array/fNDArray.h
--- a/liboctave/array/fNDArray.h	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/array/fNDArray.h	Sun Mar 24 19:34:28 2024 +0100
@@ -128,12 +128,15 @@
 
   OCTAVE_API FloatComplexNDArray fourier (int dim = 1) const;
   OCTAVE_API FloatComplexNDArray ifourier (int dim = 1) const;
+  OCTAVE_API FloatNDArray ifouriersym (int dim = 1) const;
 
   OCTAVE_API FloatComplexNDArray fourier2d () const;
   OCTAVE_API FloatComplexNDArray ifourier2d () const;
+  OCTAVE_API FloatNDArray ifourier2dsym () const;
 
   OCTAVE_API FloatComplexNDArray fourierNd () const;
   OCTAVE_API FloatComplexNDArray ifourierNd () const;
+  OCTAVE_API FloatNDArray ifourierNdsym () const;
 
   friend OCTAVE_API FloatNDArray real (const FloatComplexNDArray& a);
   friend OCTAVE_API FloatNDArray imag (const FloatComplexNDArray& a);
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/numeric/oct-fftw.cc
--- a/liboctave/numeric/oct-fftw.cc	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/numeric/oct-fftw.cc	Sun Mar 24 19:34:28 2024 +0100
@@ -75,6 +75,11 @@
   m_inplace[0] = m_inplace[1] = false;
   m_n[0] = m_n[1] = dim_vector ();
 
+  m_cplan = nullptr;
+  m_cd = m_cs = m_cr = m_ch = 0;
+  m_csimd_align = false;
+  m_cn = dim_vector ();
+
 #if defined (HAVE_FFTW3_THREADS)
   int init_ret = fftw_init_threads ();
   if (! init_ret)
@@ -105,6 +110,10 @@
   if (*plan_p)
     fftw_destroy_plan (*plan_p);
 
+  plan_p = reinterpret_cast<fftw_plan *> (&m_cplan);
+  if (*plan_p)
+    fftw_destroy_plan (*plan_p);
+
   plan_p = reinterpret_cast<fftw_plan *> (&m_plan[0]);
   if (*plan_p)
     fftw_destroy_plan (*plan_p);
@@ -138,6 +147,7 @@
       fftw_plan_with_nthreads (nt);
       // Clear the current plans.
       s_instance->m_rplan = nullptr;
+      s_instance->m_cplan = nullptr;
       s_instance->m_plan[0] = s_instance->m_plan[1] = nullptr;
     }
 #else
@@ -148,6 +158,14 @@
 #endif
 }
 
+//FIXME: Macros used for C2R but can be also used for C2C and R2C.
+#define SIMD_ALIGN_BYTES (16)
+#define SIMD_PAD(type, size) (((size) + (sizeof(type) - 1)) / (sizeof(type)))
+#define SIMD_ALIGN(value, size) ((reinterpret_cast<std::ptrdiff_t>   \
+                                 (value) + ((size) - 1)) & ~((size) - 1))
+#define SIMD_OFFSET(value, size) ((reinterpret_cast<std::ptrdiff_t>  \
+                                 (value)) & ((size) - 1))
+
 #define CHECK_SIMD_ALIGNMENT(x)                         \
   (((reinterpret_cast<std::ptrdiff_t> (x)) & 0xF) == 0)
 
@@ -392,6 +410,122 @@
   return *cur_plan_p;
 }
 
+void *
+fftw_planner::do_create_plan (const int rank, const dim_vector& dims,
+                              octave_idx_type howmany,
+                              octave_idx_type stride,
+                              octave_idx_type dist,
+                              const Complex *in, double *out)
+{
+  fftw_plan *cur_plan_p = reinterpret_cast<fftw_plan *> (&m_cplan);
+  bool create_new_plan = false;
+  bool ioalign = CHECK_SIMD_ALIGNMENT (in) && CHECK_SIMD_ALIGNMENT (out);
+
+  // Don't create a new plan if we have a non SIMD plan already but
+  // can do SIMD.  This prevents endlessly recreating plans if we
+  // change the alignment.
+
+  if (m_cplan == nullptr || m_cd != dist || m_cs != stride || m_cr != rank
+      || m_ch != howmany || ((ioalign != m_csimd_align) ? ! ioalign : false))
+    create_new_plan = true;
+  else
+    {
+      // We still might not have the same shape of array.
+
+      for (int i = 0; i < rank; i++)
+        if (dims(i) != m_cn(i))
+          {
+            create_new_plan = true;
+            break;
+          }
+    }
+
+  if (create_new_plan)
+    {
+      m_cd = dist;
+      m_cs = stride;
+      m_cr = rank;
+      m_ch = howmany;
+      m_csimd_align = ioalign;
+      m_cn = dims;
+
+      // Note reversal of dimensions for column major storage in FFTW.
+      OCTAVE_LOCAL_BUFFER (int, tmp, rank);
+      for (int i = 0, j = rank-1; i < rank; i++, j--)
+        tmp[i] = dims(j);
+
+      octave_idx_type nn = dims(0) / 2 + 1;
+      for (int i = 1; i < rank; i++)
+        nn *= dims(i);
+
+      int plan_flags = 0;
+      bool plan_destroys_in = true;
+
+      switch (m_meth)
+        {
+        case UNKNOWN:
+        case ESTIMATE:
+          plan_flags |= FFTW_ESTIMATE;
+          plan_destroys_in = false;
+          break;
+        case MEASURE:
+          plan_flags |= FFTW_MEASURE;
+          break;
+        case PATIENT:
+          plan_flags |= FFTW_PATIENT;
+          break;
+        case EXHAUSTIVE:
+          plan_flags |= FFTW_EXHAUSTIVE;
+          break;
+        case HYBRID:
+          if (nn < 8193)
+            plan_flags |= FFTW_MEASURE;
+          else
+            {
+              plan_flags |= FFTW_ESTIMATE;
+              plan_destroys_in = false;
+            }
+          break;
+        }
+
+      if (ioalign)
+        plan_flags &= ~FFTW_UNALIGNED;
+      else
+        plan_flags |= FFTW_UNALIGNED;
+
+      if (*cur_plan_p)
+        fftw_destroy_plan (*cur_plan_p);
+
+      OCTAVE_SCOPED_BUFFER_ANCHOR (Complex, itmp);
+      itmp = const_cast<Complex *> (in);
+      double *otmp = out;
+
+      if (plan_destroys_in)
+        {
+          // Create matrix with the same size and alignment as input
+          octave_idx_type in_place = reinterpret_cast<Complex *> (out) == in;
+          OCTAVE_SCOPED_BUFFER (Complex, itmp, nn * howmany +
+                                2 * SIMD_PAD(Complex, SIMD_ALIGN_BYTES));
+          itmp = reinterpret_cast<Complex *> (SIMD_ALIGN(itmp, SIMD_ALIGN_BYTES) +
+                                              SIMD_OFFSET(in, SIMD_ALIGN_BYTES));
+
+          if (in_place)
+            otmp = reinterpret_cast<double *> (itmp);
+        }
+
+      *cur_plan_p
+        = fftw_plan_many_dft_c2r (rank, tmp, howmany,
+                                  reinterpret_cast<fftw_complex *> (itmp),
+                                  nullptr, 1, nn, otmp,
+                                  nullptr, stride, dist, plan_flags);
+
+      if (*cur_plan_p == nullptr)
+        (*current_liboctave_error_handler) ("Error creating FFTW plan");
+    }
+
+  return *cur_plan_p;
+}
+
 fftw_planner::FftwMethod
 fftw_planner::do_method ()
 {
@@ -411,11 +545,13 @@
           m_meth = _meth;
           if (m_rplan)
             fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_rplan));
+          if (m_cplan)
+            fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_cplan));
           if (m_plan[0])
             fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_plan[0]));
           if (m_plan[1])
             fftw_destroy_plan (reinterpret_cast<fftw_plan> (m_plan[1]));
-          m_rplan = m_plan[0] = m_plan[1] = nullptr;
+          m_rplan = m_cplan = m_plan[0] = m_plan[1] = nullptr;
         }
     }
   else
@@ -435,6 +571,11 @@
   m_inplace[0] = m_inplace[1] = false;
   m_n[0] = m_n[1] = dim_vector ();
 
+  m_cplan = nullptr;
+  m_cd = m_cs = m_cr = m_ch = 0;
+  m_csimd_align = false;
+  m_cn = dim_vector ();
+
 #if defined (HAVE_FFTW3F_THREADS)
   int init_ret = fftwf_init_threads ();
   if (! init_ret)
@@ -460,6 +601,10 @@
   if (*plan_p)
     fftwf_destroy_plan (*plan_p);
 
+  plan_p = reinterpret_cast<fftwf_plan *> (&m_cplan);
+  if (*plan_p)
+    fftwf_destroy_plan (*plan_p);
+
   plan_p = reinterpret_cast<fftwf_plan *> (&m_plan[0]);
   if (*plan_p)
     fftwf_destroy_plan (*plan_p);
@@ -493,6 +638,7 @@
       fftwf_plan_with_nthreads (nt);
       // Clear the current plans.
       s_instance->m_rplan = nullptr;
+      s_instance->m_cplan = nullptr;
       s_instance->m_plan[0] = s_instance->m_plan[1] = nullptr;
     }
 #else
@@ -745,6 +891,122 @@
   return *cur_plan_p;
 }
 
+void *
+float_fftw_planner::do_create_plan (const int rank, const dim_vector& dims,
+                                    octave_idx_type howmany,
+                                    octave_idx_type stride,
+                                    octave_idx_type dist,
+                                    const FloatComplex *in, float *out)
+{
+  fftwf_plan *cur_plan_p = reinterpret_cast<fftwf_plan *> (&m_cplan);
+  bool create_new_plan = false;
+  bool ioalign = CHECK_SIMD_ALIGNMENT (in) && CHECK_SIMD_ALIGNMENT (out);
+
+  // Don't create a new plan if we have a non SIMD plan already but
+  // can do SIMD.  This prevents endlessly recreating plans if we
+  // change the alignment.
+
+  if (m_cplan == nullptr || m_cd != dist || m_cs != stride || m_cr != rank
+      || m_ch != howmany || ((ioalign != m_csimd_align) ? ! ioalign : false))
+    create_new_plan = true;
+  else
+    {
+      // We still might not have the same shape of array.
+
+      for (int i = 0; i < rank; i++)
+        if (dims(i) != m_cn(i))
+          {
+            create_new_plan = true;
+            break;
+          }
+    }
+
+  if (create_new_plan)
+    {
+      m_cd = dist;
+      m_cs = stride;
+      m_cr = rank;
+      m_ch = howmany;
+      m_csimd_align = ioalign;
+      m_cn = dims;
+
+      // Note reversal of dimensions for column major storage in FFTW.
+      OCTAVE_LOCAL_BUFFER (int, tmp, rank);
+      for (int i = 0, j = rank-1; i < rank; i++, j--)
+        tmp[i] = dims(j);
+
+      octave_idx_type nn = dims(0) / 2 + 1;
+      for (int i = 1; i < rank; i++)
+        nn *= dims(i);
+
+      int plan_flags = 0;
+      bool plan_destroys_in = true;
+
+      switch (m_meth)
+        {
+        case UNKNOWN:
+        case ESTIMATE:
+          plan_flags |= FFTW_ESTIMATE;
+          plan_destroys_in = false;
+          break;
+        case MEASURE:
+          plan_flags |= FFTW_MEASURE;
+          break;
+        case PATIENT:
+          plan_flags |= FFTW_PATIENT;
+          break;
+        case EXHAUSTIVE:
+          plan_flags |= FFTW_EXHAUSTIVE;
+          break;
+        case HYBRID:
+          if (nn < 8193)
+            plan_flags |= FFTW_MEASURE;
+          else
+            {
+              plan_flags |= FFTW_ESTIMATE;
+              plan_destroys_in = false;
+            }
+          break;
+        }
+
+      if (ioalign)
+        plan_flags &= ~FFTW_UNALIGNED;
+      else
+        plan_flags |= FFTW_UNALIGNED;
+
+      if (*cur_plan_p)
+        fftwf_destroy_plan (*cur_plan_p);
+
+      OCTAVE_SCOPED_BUFFER_ANCHOR (FloatComplex, itmp);
+      itmp = const_cast<FloatComplex *> (in);
+      float *otmp = out;
+
+      if (plan_destroys_in)
+        {
+          // Create matrix with the same size and alignment as input
+          octave_idx_type in_place = reinterpret_cast<FloatComplex *> (out) == in;
+          OCTAVE_SCOPED_BUFFER (FloatComplex, itmp, nn * howmany +
+                                2 * SIMD_PAD(FloatComplex, SIMD_ALIGN_BYTES));
+          itmp = reinterpret_cast<FloatComplex *> (SIMD_ALIGN(itmp, SIMD_ALIGN_BYTES) +
+                                                   SIMD_OFFSET(in, SIMD_ALIGN_BYTES));
+
+          if (in_place)
+            otmp = reinterpret_cast<float *> (itmp);
+        }
+
+      *cur_plan_p
+        = fftwf_plan_many_dft_c2r (rank, tmp, howmany,
+                                   reinterpret_cast<fftwf_complex *> (itmp),
+                                   nullptr, 1, nn, otmp,
+                                   nullptr, stride, dist, plan_flags);
+
+      if (*cur_plan_p == nullptr)
+        (*current_liboctave_error_handler) ("Error creating FFTW plan");
+    }
+
+  return *cur_plan_p;
+}
+
 float_fftw_planner::FftwMethod
 float_fftw_planner::do_method ()
 {
@@ -764,11 +1026,13 @@
           m_meth = _meth;
           if (m_rplan)
             fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (m_rplan));
+          if (m_cplan)
+            fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (m_cplan));
           if (m_plan[0])
             fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (m_plan[0]));
           if (m_plan[1])
             fftwf_destroy_plan (reinterpret_cast<fftwf_plan> (m_plan[1]));
-          m_rplan = m_plan[0] = m_plan[1] = nullptr;
+          m_rplan = m_cplan = m_plan[0] = m_plan[1] = nullptr;
         }
     }
   else
@@ -778,6 +1042,178 @@
 
 template <typename T>
 static inline void
+c2r_convert_1d (const T *in, T *out, const std::size_t nsamples,
+                const std::size_t npts, const octave_idx_type stride,
+                const octave_idx_type dist)
+{
+  T *ptr = out;
+
+  octave_quit ();
+
+  for (std::size_t i = 0; i < nsamples; i++)
+    for (std::size_t j = 0; j < npts/2+1; j++)
+      *ptr++ = in[j*stride + i*dist];
+
+  octave_quit ();
+}
+
+template <typename T>
+static inline void
+c2r_convert_Nd (T *Y, const dim_vector& dv)
+{
+  const std::size_t nc = dv(0);
+  const std::size_t nr = dv(1);
+  const std::size_t np = (dv.ndims () > 2 ? dv.numel () / nc / nr : 1);
+  T *ptr = Y;
+
+  octave_quit ();
+
+  for (std::size_t i = 0; i < np; i++)
+    for (std::size_t j = 0; j < nr; j++)
+      for (std::size_t k = 0; k < nc/2+1; k++)
+        *ptr++ = Y[k + (j + i*nr)*nc];
+
+  octave_quit ();
+}
+
+// Note this is a recursive function.
+template <typename T>
+static void
+make_subspace_symmetric (const T *in, T *out, const dim_vector& dv,
+                         std::size_t *sub_dims, std::size_t *sub_n,
+                         std::size_t nsubdims, std::size_t ofs,
+                         bool lastdim_odd)
+{
+  const std::size_t ndim = dv.ndims ();
+
+  // The conjugate values from the lower part of the subspace are
+  // cross mirrored along the highest dimension into the upper part.
+
+  OCTAVE_LOCAL_BUFFER (std::size_t, tmp, ndim);
+  tmp[0] = 1;
+  std::size_t nn = 1;
+  for (std::size_t i = 0; i < ndim-1; i++)
+    {
+      nn *= dv(i);
+      tmp[i+1] = nn;
+    }
+
+  OCTAVE_LOCAL_BUFFER (std::size_t, mul, nsubdims);
+  OCTAVE_LOCAL_BUFFER (std::size_t, coord, nsubdims);
+  std::size_t nel = 1;
+  for (std::size_t k = 0; k < nsubdims; k++)
+    {
+      mul[k] = tmp[sub_dims[k]];
+      coord[k] = 1;
+      nel *= sub_n[k];
+    }
+
+  // Main loop with copy action.
+  for (std::size_t i = 0; i < nel; i++)
+    {
+      std::size_t id = 0;
+      std::size_t idm = 0;
+      for (std::size_t k = 0; k < nsubdims; k++)
+        {
+          id += mul[k] * coord[k];
+          idm += mul[k] * (dv(sub_dims[k]) - coord[k]);
+        }
+      // Make the mirror point conjugate complex.
+      out[idm + ofs] = conj (in[id + ofs]);
+      out[id + ofs] = in[id + ofs];
+
+      // Update coordinates associated with current point.
+      std::size_t dim = 0;
+      while (dim < nsubdims)
+        {
+          if (coord[dim] < sub_n[dim])
+            {
+              coord[dim]++;
+              break;
+            }
+          else
+            {
+              // At end of the current dimension, overflow into the next dimension.
+              coord[dim] = 1;
+              dim++;
+            }
+        }
+    }
+
+  // If the last subspace size is odd, there is a remainder and a
+  // sub-subspace is created.
+
+  if ((nsubdims == 1) && lastdim_odd)
+    {
+      // Terminate when a one-dimensional subspace is reached.
+      const std::size_t id = mul[0] * (1 + sub_n[0]);
+      out[id + ofs] = real (in[id + ofs]);
+    }
+  else if (lastdim_odd)
+    {
+      // Move the sub-subspace origin to the correct position.
+      ofs += mul[nsubdims-1] * (1 + sub_n[nsubdims-1]);
+
+      nsubdims--;
+
+      const std::size_t n = dv(sub_dims[nsubdims-1]) - 1;
+      const std::size_t n2 = n / 2;
+      lastdim_odd = static_cast<bool> (n - 2 * n2);
+
+      for (std::size_t k = 0; k < nsubdims-1; k++)
+        sub_n[k] = dv(sub_dims[k]) - 1;
+      sub_n[nsubdims-1] = n2;
+
+      make_subspace_symmetric (in, out, dv, sub_dims, sub_n, nsubdims,
+                               ofs, lastdim_odd);
+    }
+}
+
+template <typename T>
+static inline void
+makesymmetric_Nd (const T *in, T *out, const dim_vector& dv)
+{
+  const std::size_t ndim = dv.ndims ();
+
+  OCTAVE_LOCAL_BUFFER (std::size_t, sub_dims, ndim);
+  OCTAVE_LOCAL_BUFFER (std::size_t, sub_n, ndim);
+
+  octave_quit ();
+
+  out[0] = real (in[0]);
+
+  // For all subspace permutations.
+  const std::size_t bit = 1;
+  for (std::size_t i = 1; i < (bit << ndim); i++)
+    {
+      std::size_t nsubdims = 0;
+      for (std::size_t j = 0; j < ndim; j++)
+        if (i & (bit << j))
+          {
+            sub_dims[nsubdims] = j;
+            nsubdims++;
+          }
+
+      const std::size_t n = dv(sub_dims[nsubdims-1]) - 1;
+      const std::size_t n2 = n / 2;
+      const bool lastdim_odd = static_cast<bool> (n - 2 * n2);
+
+      // The last subspace size is truncated.
+      for (std::size_t k = 0; k < nsubdims-1; k++)
+        sub_n[k] = dv(sub_dims[k]) - 1;
+      sub_n[nsubdims-1] = n2;
+
+      octave_quit ();
+
+      make_subspace_symmetric (in, out, dv, sub_dims, sub_n, nsubdims,
+                               0, lastdim_odd);
+    }
+
+  octave_quit ();
+}
+
+template <typename T>
+static inline void
 convert_packcomplex_1d (T *out, std::size_t nr, std::size_t nc,
                         octave_idx_type stride, octave_idx_type dist)
 {
@@ -922,6 +1358,35 @@
 }
 
 int
+fftw::ifftsym (const Complex *in, double *out, std::size_t npts,
+               std::size_t nsamples, octave_idx_type stride,
+               octave_idx_type dist)
+{
+  dist = (dist < 0 ? npts : dist);
+
+  dim_vector dv (npts, 1);
+
+  OCTAVE_LOCAL_BUFFER (Complex, in_c2r, nsamples * (npts / 2 + 1) +
+                       SIMD_PAD(Complex, SIMD_ALIGN_BYTES));
+  in_c2r = reinterpret_cast<Complex *> (SIMD_ALIGN(in_c2r, SIMD_ALIGN_BYTES));
+  c2r_convert_1d (in, in_c2r, nsamples, npts, stride, dist);
+
+  void *vplan = fftw_planner::create_plan (1, dv, nsamples, stride, dist, in_c2r, out);
+  fftw_plan m_plan = reinterpret_cast<fftw_plan> (vplan);
+
+  fftw_execute_dft_c2r (m_plan,
+                        reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in_c2r)),
+                        reinterpret_cast<double *> (out));
+
+  const double scale = npts;
+  for (std::size_t j = 0; j < nsamples; j++)
+    for (std::size_t i = 0; i < npts; i++)
+      out[i*stride + j*dist] /= scale;
+
+  return 0;
+}
+
+int
 fftw::fftNd (const double *in, Complex *out, const int rank,
              const dim_vector& dv)
 {
@@ -992,6 +1457,35 @@
 }
 
 int
+fftw::ifftNdsym (const Complex *in, double *out, const int rank,
+                 const dim_vector& dv)
+{
+  octave_idx_type dist = 1;
+  for (int i = 0; i < rank; i++)
+    dist *= dv(i);
+
+  const std::size_t npts = dv.numel ();
+
+  OCTAVE_LOCAL_BUFFER (Complex, in_c2r, npts + SIMD_PAD(Complex, SIMD_ALIGN_BYTES));
+  in_c2r = reinterpret_cast<Complex *> (SIMD_ALIGN(in_c2r, SIMD_ALIGN_BYTES));
+  makesymmetric_Nd (in, in_c2r, dv);
+  c2r_convert_Nd (in_c2r, dv);
+
+  void *vplan = fftw_planner::create_plan (rank, dv, 1, 1, dist, in_c2r, out);
+  fftw_plan m_plan = reinterpret_cast<fftw_plan> (vplan);
+
+  fftw_execute_dft_c2r (m_plan,
+                        reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in_c2r)),
+                        reinterpret_cast<double *> (out));
+
+  const double scale = npts;
+  for (std::size_t i = 0; i < npts; i++)
+    out[i] /= scale;
+
+  return 0;
+}
+
+int
 fftw::fft (const float *in, FloatComplex *out, std::size_t npts,
            std::size_t nsamples, octave_idx_type stride,
            octave_idx_type dist)
@@ -1059,6 +1553,35 @@
 }
 
 int
+fftw::ifftsym (const FloatComplex *in, float *out, std::size_t npts,
+               std::size_t nsamples, octave_idx_type stride,
+               octave_idx_type dist)
+{
+  dist = (dist < 0 ? npts : dist);
+
+  dim_vector dv (npts, 1);
+
+  OCTAVE_LOCAL_BUFFER (FloatComplex, in_c2r, nsamples * (npts / 2 + 1) +
+                       SIMD_PAD(FloatComplex, SIMD_ALIGN_BYTES));
+  in_c2r = reinterpret_cast<FloatComplex *> (SIMD_ALIGN(in_c2r, SIMD_ALIGN_BYTES));
+  c2r_convert_1d (in, in_c2r, nsamples, npts, stride, dist);
+
+  void *vplan = float_fftw_planner::create_plan (1, dv, nsamples, stride, dist, in_c2r, out);
+  fftwf_plan m_plan = reinterpret_cast<fftwf_plan> (vplan);
+
+  fftwf_execute_dft_c2r (m_plan,
+                         reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *> (in_c2r)),
+                         reinterpret_cast<float *> (out));
+
+  const float scale = npts;
+  for (std::size_t j = 0; j < nsamples; j++)
+    for (std::size_t i = 0; i < npts; i++)
+      out[i*stride + j*dist] /= scale;
+
+  return 0;
+}
+
+int
 fftw::fftNd (const float *in, FloatComplex *out, const int rank,
              const dim_vector& dv)
 {
@@ -1128,6 +1651,35 @@
   return 0;
 }
 
+int
+fftw::ifftNdsym (const FloatComplex *in, float *out, const int rank,
+                 const dim_vector& dv)
+{
+  octave_idx_type dist = 1;
+  for (int i = 0; i < rank; i++)
+    dist *= dv(i);
+
+  const std::size_t npts = dv.numel ();
+
+  OCTAVE_LOCAL_BUFFER (FloatComplex, in_c2r, npts + SIMD_PAD(FloatComplex, SIMD_ALIGN_BYTES));
+  in_c2r = reinterpret_cast<FloatComplex *> (SIMD_ALIGN(in_c2r, SIMD_ALIGN_BYTES));
+  makesymmetric_Nd (in, in_c2r, dv);
+  c2r_convert_Nd (in_c2r, dv);
+
+  void *vplan = float_fftw_planner::create_plan (rank, dv, 1, 1, dist, in_c2r, out);
+  fftwf_plan m_plan = reinterpret_cast<fftwf_plan> (vplan);
+
+  fftwf_execute_dft_c2r (m_plan,
+                         reinterpret_cast<fftwf_complex *> (const_cast<FloatComplex *> (in_c2r)),
+                         reinterpret_cast<float *> (out));
+
+  const float scale = npts;
+  for (std::size_t i = 0; i < npts; i++)
+    out[i] /= scale;
+
+  return 0;
+}
+
 #endif
 
 std::string
diff -r 399be7cc310f -r 342e2e9bff3d liboctave/numeric/oct-fftw.h
--- a/liboctave/numeric/oct-fftw.h	Sat Mar 23 13:54:13 2024 -0400
+++ b/liboctave/numeric/oct-fftw.h	Sun Mar 24 19:34:28 2024 +0100
@@ -84,6 +84,17 @@
            : nullptr;
   }
 
+  static void *
+  create_plan (const int rank, const dim_vector& dims,
+               octave_idx_type howmany, octave_idx_type stride,
+               octave_idx_type dist, const Complex *in, double *out)
+  {
+    return instance_ok ()
+           ? s_instance->do_create_plan (rank, dims, howmany, stride, dist,
+                                         in, out)
+           : nullptr;
+  }
+
   static FftwMethod method ()
   {
     static FftwMethod dummy;
@@ -123,6 +134,11 @@
                   octave_idx_type howmany, octave_idx_type stride,
                   octave_idx_type dist, const double *in, Complex *out);
 
+  void *
+  do_create_plan (const int rank, const dim_vector& dims,
+                  octave_idx_type howmany, octave_idx_type stride,
+                  octave_idx_type dist, const Complex *in, double *out);
+
   FftwMethod do_method ();
 
   FftwMethod do_method (FftwMethod meth);
@@ -152,6 +168,26 @@
   bool m_simd_align[2];
   bool m_inplace[2];
 
+  // Plan for ifft symmetric
+  void *m_cplan;
+
+  // dist
+  octave_idx_type m_cd;
+
+  // stride
+  octave_idx_type m_cs;
+
+  // rank
+  int m_cr;
+
+  // howmany
+  octave_idx_type m_ch;
+
+  // dims
+  dim_vector m_cn;
+
+  bool m_csimd_align;
+
   // Plan for fft of real values
   void *m_rplan;
 
@@ -224,6 +260,17 @@
            : nullptr;
   }
 
+  static void *
+  create_plan (const int rank, const dim_vector& dims,
+               octave_idx_type howmany, octave_idx_type stride,
+               octave_idx_type dist, const FloatComplex *in, float *out)
+  {
+    return instance_ok ()
+           ? s_instance->do_create_plan (rank, dims, howmany, stride, dist,
+                                         in, out)
+           : nullptr;
+  }
+
   static FftwMethod method ()
   {
     static FftwMethod dummy;
@@ -263,6 +310,11 @@
                   octave_idx_type howmany, octave_idx_type stride,
                   octave_idx_type dist, const float *in, FloatComplex *out);
 
+  void *
+  do_create_plan (const int rank, const dim_vector& dims,
+                  octave_idx_type howmany, octave_idx_type stride,
+                  octave_idx_type dist, const FloatComplex *in, float *out);
+
   FftwMethod do_method ();
 
   FftwMethod do_method (FftwMethod meth);
@@ -292,6 +344,26 @@
   bool m_simd_align[2];
   bool m_inplace[2];
 
+  // Plan for ifft symmetric
+  void *m_cplan;
+
+  // dist
+  octave_idx_type m_cd;
+
+  // stride
+  octave_idx_type m_cs;
+
+  // rank
+  int m_cr;
+
+  // howmany
+  octave_idx_type m_ch;
+
+  // dims
+  dim_vector m_cn;
+
+  bool m_csimd_align;
+
   // Plan for fft of real values
   void *m_rplan;
 
@@ -332,12 +404,17 @@
   static int ifft (const Complex *in, Complex *out, std::size_t npts,
                    std::size_t nsamples = 1, octave_idx_type stride = 1,
                    octave_idx_type dist = -1);
+  static int ifftsym (const Complex *in, double *out, std::size_t npts,
+                      std::size_t nsamples = 1, octave_idx_type stride = 1,
+                      octave_idx_type dist = -1);
 
   static int fftNd (const double *, Complex *, const int, const dim_vector&);
   static int fftNd (const Complex *, Complex *, const int,
                     const dim_vector&);
   static int ifftNd (const Complex *, Complex *, const int,
                      const dim_vector&);
+  static int ifftNdsym (const Complex *, double *, const int,
+                        const dim_vector&);
 
   static int fft (const float *in, FloatComplex *out, std::size_t npts,
                   std::size_t nsamples = 1, octave_idx_type stride = 1,
@@ -348,6 +425,9 @@
   static int ifft (const FloatComplex *in, FloatComplex *out, std::size_t npts,
                    std::size_t nsamples = 1, octave_idx_type stride = 1,
                    octave_idx_type dist = -1);
+  static int ifftsym (const FloatComplex *in, float *out, std::size_t npts,
+                      std::size_t nsamples = 1, octave_idx_type stride = 1,
+                      octave_idx_type dist = -1);
 
   static int fftNd (const float *, FloatComplex *, const int,
                     const dim_vector&);
@@ -355,6 +435,9 @@
                     const dim_vector&);
   static int ifftNd (const FloatComplex *, FloatComplex *, const int,
                      const dim_vector&);
+  static int ifftNdsym (const FloatComplex *, float *, const int,
+                        const dim_vector&);
+
 };
 
 extern OCTAVE_API std::string fftw_version ();
