# HG changeset patch
# User jwe
# Date 778999439 0
# Node ID 01da6806197be9ff073071a580320ca7b74c413e
# Parent  d31fde825c536708a276892b4f436b2e4464c25a
[project @ 1994-09-08 04:43:13 by jwe]

diff -r d31fde825c53 -r 01da6806197b liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc	Thu Sep 08 04:42:58 1994 +0000
+++ b/liboctave/CMatrix.cc	Thu Sep 08 04:43:59 1994 +0000
@@ -1021,6 +1021,115 @@
   return ComplexMatrix (tmp_data, nr, nc);
 }
 
+ComplexMatrix
+ComplexMatrix::fourier2d (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  int npts, nsamples;
+  if (nr == 1 || nc == 1)
+    {
+      npts = nr > nc ? nr : nc;
+      nsamples = 1;
+    }
+  else
+    {
+      npts = nr;
+      nsamples = nc;
+    }
+
+  int nn = 4*npts+15;
+  Complex *wsave = new Complex [nn];
+  Complex *tmp_data = dup (data (), length ());
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (int j = 0; j < nsamples; j++)
+    F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave);
+
+  delete [] wsave;
+
+  npts = nc;
+  nsamples = nr;
+  nn = 4*npts+15;
+  wsave = new Complex [nn];
+  Complex *row = new Complex[npts];
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (j = 0; j < nsamples; j++)
+    {
+      for (int i = 0; i < npts; i++)
+	row[i] = tmp_data[i*nr + j];
+
+      F77_FCN (cfftf) (&npts, row, wsave);
+
+      for (i = 0; i < npts; i++)
+	tmp_data[i*nr + j] = row[i];
+    }
+
+  delete [] wsave;
+  delete [] row;
+
+  return ComplexMatrix (tmp_data, nr, nc);
+}
+
+ComplexMatrix
+ComplexMatrix::ifourier2d (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  int npts, nsamples;
+  if (nr == 1 || nc == 1)
+    {
+      npts = nr > nc ? nr : nc;
+      nsamples = 1;
+    }
+  else
+    {
+      npts = nr;
+      nsamples = nc;
+    }
+
+  int nn = 4*npts+15;
+  Complex *wsave = new Complex [nn];
+  Complex *tmp_data = dup (data (), length ());
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (int j = 0; j < nsamples; j++)
+    F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave);
+
+  delete [] wsave;
+
+  for (j = 0; j < npts*nsamples; j++)
+    tmp_data[j] = tmp_data[j] / (double) npts;
+
+  npts = nc;
+  nsamples = nr;
+  nn = 4*npts+15;
+  wsave = new Complex [nn];
+  Complex *row = new Complex[npts];
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (j = 0; j < nsamples; j++)
+    {
+      for (int i = 0; i < npts; i++)
+	row[i] = tmp_data[i*nr + j];
+
+      F77_FCN (cfftb) (&npts, row, wsave);
+
+      for (i = 0; i < npts; i++)
+	tmp_data[i*nr + j] = row[i] / (double) npts;
+    }
+
+  delete [] wsave;
+  delete [] row;
+
+  return ComplexMatrix (tmp_data, nr, nc);
+}
+
 ComplexDET
 ComplexMatrix::determinant (void) const
 {
diff -r d31fde825c53 -r 01da6806197b liboctave/CMatrix.h
--- a/liboctave/CMatrix.h	Thu Sep 08 04:42:58 1994 +0000
+++ b/liboctave/CMatrix.h	Thu Sep 08 04:43:59 1994 +0000
@@ -137,6 +137,9 @@
   ComplexMatrix fourier (void) const;
   ComplexMatrix ifourier (void) const;
 
+  ComplexMatrix fourier2d (void) const;
+  ComplexMatrix ifourier2d (void) const;
+
   ComplexDET determinant (void) const;
   ComplexDET determinant (int& info) const;
   ComplexDET determinant (int& info, double& rcond) const;
diff -r d31fde825c53 -r 01da6806197b liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc	Thu Sep 08 04:42:58 1994 +0000
+++ b/liboctave/dMatrix.cc	Thu Sep 08 04:43:59 1994 +0000
@@ -670,6 +670,115 @@
   return ComplexMatrix (tmp_data, nr, nc);
 }
 
+ComplexMatrix
+Matrix::fourier2d (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  int npts, nsamples;
+  if (nr == 1 || nc == 1)
+    {
+      npts = nr > nc ? nr : nc;
+      nsamples = 1;
+    }
+  else
+    {
+      npts = nr;
+      nsamples = nc;
+    }
+
+  int nn = 4*npts+15;
+  Complex *wsave = new Complex [nn];
+  Complex *tmp_data = make_complex (data (), length ());
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (int j = 0; j < nsamples; j++)
+    F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave);
+
+  delete [] wsave;
+
+  npts = nc;
+  nsamples = nr;
+  nn = 4*npts+15;
+  wsave = new Complex [nn];
+  Complex *row = new Complex[npts];
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (j = 0; j < nsamples; j++)
+    {
+      for (int i = 0; i < npts; i++)
+	row[i] = tmp_data[i*nr + j];
+
+      F77_FCN (cfftf) (&npts, row, wsave);
+
+      for (i = 0; i < npts; i++)
+	tmp_data[i*nr + j] = row[i];
+    }
+
+  delete [] wsave;
+  delete [] row;
+
+  return ComplexMatrix (tmp_data, nr, nc);
+}
+
+ComplexMatrix
+Matrix::ifourier2d (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  int npts, nsamples;
+  if (nr == 1 || nc == 1)
+    {
+      npts = nr > nc ? nr : nc;
+      nsamples = 1;
+    }
+  else
+    {
+      npts = nr;
+      nsamples = nc;
+    }
+
+  int nn = 4*npts+15;
+  Complex *wsave = new Complex [nn];
+  Complex *tmp_data = make_complex (data (), length ());
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (int j = 0; j < nsamples; j++)
+    F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave);
+
+  delete [] wsave;
+
+  for (j = 0; j < npts*nsamples; j++)
+    tmp_data[j] = tmp_data[j] / (double) npts;
+
+  npts = nc;
+  nsamples = nr;
+  nn = 4*npts+15;
+  wsave = new Complex [nn];
+  Complex *row = new Complex[npts];
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (j = 0; j < nsamples; j++)
+    {
+      for (int i = 0; i < npts; i++)
+	row[i] = tmp_data[i*nr + j];
+
+      F77_FCN (cfftb) (&npts, row, wsave);
+
+      for (i = 0; i < npts; i++)
+	tmp_data[i*nr + j] = row[i] / (double) npts;
+    }
+
+  delete [] wsave;
+  delete [] row;
+
+  return ComplexMatrix (tmp_data, nr, nc);
+}
+
 DET
 Matrix::determinant (void) const
 {
diff -r d31fde825c53 -r 01da6806197b liboctave/dMatrix.h
--- a/liboctave/dMatrix.h	Thu Sep 08 04:42:58 1994 +0000
+++ b/liboctave/dMatrix.h	Thu Sep 08 04:43:59 1994 +0000
@@ -112,6 +112,9 @@
   ComplexMatrix fourier (void) const;
   ComplexMatrix ifourier (void) const;
 
+  ComplexMatrix fourier2d (void) const;
+  ComplexMatrix ifourier2d (void) const;
+
   DET determinant (void) const;
   DET determinant (int& info) const;
   DET determinant (int& info, double& rcond) const;
diff -r d31fde825c53 -r 01da6806197b src/fft.cc
--- a/src/fft.cc	Thu Sep 08 04:42:58 1994 +0000
+++ b/src/fft.cc	Thu Sep 08 04:43:59 1994 +0000
@@ -36,12 +36,16 @@
 #include "help.h"
 #include "defun-dld.h"
 
-DEFUN_DLD ("fft", Ffft, Sfft, 2, 1,
-  "fft (X): fast fourier transform of a vector")
+// This function should be merged with Fifft.
+
+DEFUN_DLD ("fft", Ffft, Sfft, 3, 1,
+  "fft (X [, N]): fast fourier transform of a vector")
 {
   Octave_object retval;
 
-  if (args.length () != 2)
+  int nargin = args.length ();
+
+  if (nargin < 2 || nargin > 3)
     {
       print_usage ("fft");
       return retval;
@@ -49,17 +53,34 @@
 
   tree_constant arg = args(1);
 
-  if (empty_arg ("fft", arg.rows (), arg.columns ()) < 0)
+  int n_points = arg.rows ();
+  if (nargin == 3)
+    n_points = NINT (args(2).double_value ());
+
+  if (error_state)
     return retval;
 
+  if (n_points < 0)
+    {
+      error ("fft: number of points must be greater than zero");
+      return retval;
+    }
+
+  int arg_is_empty = empty_arg ("fft", arg.rows (), arg.columns ());
+
+  if (arg_is_empty < 0)
+    return retval;
+  else if (arg_is_empty || n_points == 0)
+    return Matrix ();
+
   if (arg.is_real_type ())
     {
       Matrix m = arg.matrix_value ();
 
       if (! error_state)
 	{
-	  ComplexMatrix mfft = m.fourier ();
-	  retval = mfft;
+	  m.resize (n_points, m.columns (), 0.0);
+	  retval = m.fourier ();
 	}
     }
   else if (arg.is_complex_type ())
@@ -68,8 +89,8 @@
 
       if (! error_state)
 	{
-	  ComplexMatrix mfft = m.fourier ();
-	  retval = mfft;
+	  m.resize (n_points, m.columns (), 0.0);
+	  retval = m.fourier ();
 	}
     }
   else
diff -r d31fde825c53 -r 01da6806197b src/ifft.cc
--- a/src/ifft.cc	Thu Sep 08 04:42:58 1994 +0000
+++ b/src/ifft.cc	Thu Sep 08 04:43:59 1994 +0000
@@ -36,12 +36,16 @@
 #include "help.h"
 #include "defun-dld.h"
 
-DEFUN_DLD ("ifft", Fifft, Sifft,2, 1,
-  "ifft (X): inverse fast fourier transform of a vector")
+// This function should be merged with Ffft.
+
+DEFUN_DLD ("ifft", Fifft, Sifft, 3, 1,
+  "ifft (X [, N]): inverse fast fourier transform of a vector")
 {
   Octave_object retval;
 
-  if (args.length () != 2)
+  int nargin = args.length ();
+
+  if (nargin < 2 || nargin > 3)
     {
       print_usage ("ifft");
       return retval;
@@ -49,17 +53,34 @@
 
   tree_constant arg = args(1);
     
-  if (empty_arg ("ifft", arg.rows (), arg.columns ()) < 0)
+  int n_points = arg.rows ();
+  if (nargin == 3)
+    n_points = NINT (args(2).double_value ());
+
+  if (error_state)
     return retval;
 
+  if (n_points < 0)
+    {
+      error ("ifft: number of points must be greater than zero");
+      return retval;
+    }
+
+  int arg_is_empty = empty_arg ("ifft", arg.rows (), arg.columns ());
+
+  if (arg_is_empty < 0)
+    return retval;
+  else if (arg_is_empty || n_points == 0)
+    return Matrix ();
+
   if (arg.is_real_type ())
     {
       Matrix m = arg.matrix_value ();
 
       if (! error_state)
 	{
-	  ComplexMatrix mifft = m.ifourier ();
-	  retval = mifft;
+	  m.resize (n_points, m.columns (), 0.0);
+	  retval = m.ifourier ();
 	}
     }
   else if (arg.is_complex_type ())
@@ -68,8 +89,8 @@
 
       if (! error_state)
 	{
-	  ComplexMatrix mifft = m.ifourier ();
-	  retval = mifft;
+	  m.resize (n_points, m.columns (), 0.0);
+	  retval = m.ifourier ();
 	}
     }
   else