*** ./liboctave/ChangeLog.orig 2004-01-30 15:24:31.000000000 +0100
--- ./liboctave/ChangeLog 2004-02-10 11:08:02.000000000 +0100
***************
*** 1,3 ****
--- 1,21 ----
+ 2004-02-10 David Bateman
+
+ * oct-fftw.cc: Add support for FFTW 3.x. Include the ability to
+ use the real to complex transform for fft's of real matrices
+ * oct-fftw.h: Decls updated for the above
+
+ * dMatrix.cc: fft's use real to complex transforms. 1D fft of a
+ matrix done as single call rather than loop. Update for FFTW 3.x
+ * CMatrix.cc: 1D fft of a matrix done as single call rather than
+ loop. Update for FFTW 3.x
+
+ * dNDArray.cc (fourier, ifourier, fourierNd, ifouriourNd: New
+ fourier transform functions for Nd arrays
+ * dNArray.h (fourier, ifourier, fourierNd, ifouriourNd: decls
+ * CNDArray.cc (fourier, ifourier, fourierNd, ifouriourNd: New
+ fourier transform functions for complex Nd arrays
+ * CNArray.h (fourier, ifourier, fourierNd, ifouriourNd: decls
+
2004-01-22 John W. Eaton
* Array.cc (Array::assign2, Array::assignN):
*** ./liboctave/dMatrix.cc.orig 2004-01-28 15:17:03.000000000 +0100
--- ./liboctave/dMatrix.cc 2004-02-10 11:54:20.000000000 +0100
***************
*** 52,58 ****
#include "oct-cmplx.h"
#include "quit.h"
! #ifdef HAVE_FFTW
#include "oct-fftw.h"
#endif
--- 52,58 ----
#include "oct-cmplx.h"
#include "quit.h"
! #if defined (HAVE_FFTW3)
#include "oct-fftw.h"
#endif
***************
*** 766,772 ****
}
}
! #ifdef HAVE_FFTW
ComplexMatrix
Matrix::fourier (void) const
--- 766,772 ----
}
}
! #if defined (HAVE_FFTW3)
ComplexMatrix
Matrix::fourier (void) const
***************
*** 789,804 ****
nsamples = nc;
}
! ComplexMatrix tmp (*this);
! Complex *in (tmp.fortran_vec ());
Complex *out (retval.fortran_vec ());
! for (size_t i = 0; i < nsamples; i++)
! {
! OCTAVE_QUIT;
!
! octave_fftw::fft (&in[npts * i], &out[npts * i], npts);
! }
return retval;
}
--- 789,798 ----
nsamples = nc;
}
! const double *in (fortran_vec ());
Complex *out (retval.fortran_vec ());
! octave_fftw::fft (in, out, npts, nsamples);
return retval;
}
***************
*** 828,839 ****
Complex *in (tmp.fortran_vec ());
Complex *out (retval.fortran_vec ());
! for (size_t i = 0; i < nsamples; i++)
! {
! OCTAVE_QUIT;
!
! octave_fftw::ifft (&in[npts * i], &out[npts * i], npts);
! }
return retval;
}
--- 822,828 ----
Complex *in (tmp.fortran_vec ());
Complex *out (retval.fortran_vec ());
! octave_fftw::ifft (in, out, npts, nsamples);
return retval;
}
***************
*** 841,853 ****
ComplexMatrix
Matrix::fourier2d (void) const
{
! int nr = rows ();
! int nc = cols ();
! ComplexMatrix retval (*this);
! // Note the order of passing the rows and columns to account for
! // column-major storage.
! octave_fftw::fft2d (retval.fortran_vec (), nc, nr);
return retval;
}
--- 830,840 ----
ComplexMatrix
Matrix::fourier2d (void) const
{
! dim_vector dv(rows (), cols ());
! const double *in = fortran_vec ();
! ComplexMatrix retval (rows (), cols ());
! octave_fftw::fftNd (in, retval.fortran_vec (), 2, dv);
return retval;
}
***************
*** 855,867 ****
ComplexMatrix
Matrix::ifourier2d (void) const
{
! int nr = rows ();
! int nc = cols ();
ComplexMatrix retval (*this);
! // Note the order of passing the rows and columns to account for
! // column-major storage.
! octave_fftw::ifft2d (retval.fortran_vec (), nc, nr);
return retval;
}
--- 842,853 ----
ComplexMatrix
Matrix::ifourier2d (void) const
{
! dim_vector dv(rows (), cols ());
ComplexMatrix retval (*this);
! Complex *out (retval.fortran_vec ());
!
! octave_fftw::ifftNd (out, out, 2, dv);
return retval;
}
*** ./liboctave/CMatrix.cc.orig 2004-01-28 15:17:09.000000000 +0100
--- ./liboctave/CMatrix.cc 2004-02-10 11:54:41.000000000 +0100
***************
*** 56,62 ****
#include "mx-inlines.cc"
#include "oct-cmplx.h"
! #ifdef HAVE_FFTW
#include "oct-fftw.h"
#endif
--- 56,62 ----
#include "mx-inlines.cc"
#include "oct-cmplx.h"
! #if defined (HAVE_FFTW3)
#include "oct-fftw.h"
#endif
***************
*** 1102,1108 ****
return retval;
}
! #ifdef HAVE_FFTW
ComplexMatrix
ComplexMatrix::fourier (void) const
--- 1102,1108 ----
return retval;
}
! #if defined (HAVE_FFTW3)
ComplexMatrix
ComplexMatrix::fourier (void) const
***************
*** 1128,1139 ****
const Complex *in (data ());
Complex *out (retval.fortran_vec ());
! for (size_t i = 0; i < nsamples; i++)
! {
! OCTAVE_QUIT;
!
! octave_fftw::fft (&in[npts * i], &out[npts * i], npts);
! }
return retval;
}
--- 1128,1134 ----
const Complex *in (data ());
Complex *out (retval.fortran_vec ());
! octave_fftw::fft (in, out, npts, nsamples);
return retval;
}
***************
*** 1162,1173 ****
const Complex *in (data ());
Complex *out (retval.fortran_vec ());
! for (size_t i = 0; i < nsamples; i++)
! {
! OCTAVE_QUIT;
!
! octave_fftw::ifft (&in[npts * i], &out[npts * i], npts);
! }
return retval;
}
--- 1157,1163 ----
const Complex *in (data ());
Complex *out (retval.fortran_vec ());
! octave_fftw::ifft (in, out, npts, nsamples);
return retval;
}
***************
*** 1175,1187 ****
ComplexMatrix
ComplexMatrix::fourier2d (void) const
{
! int nr = rows ();
! int nc = cols ();
! ComplexMatrix retval (*this);
! // Note the order of passing the rows and columns to account for
! // column-major storage.
! octave_fftw::fft2d (retval.fortran_vec (), nc, nr);
return retval;
}
--- 1165,1177 ----
ComplexMatrix
ComplexMatrix::fourier2d (void) const
{
! dim_vector dv(rows (), cols ());
! ComplexMatrix retval (rows (), cols ());
! const Complex *in (data ());
! Complex *out (retval.fortran_vec ());
!
! octave_fftw::fftNd (in, out, 2, dv);
return retval;
}
***************
*** 1189,1201 ****
ComplexMatrix
ComplexMatrix::ifourier2d (void) const
{
! int nr = rows ();
! int nc = cols ();
! ComplexMatrix retval (*this);
! // Note the order of passing the rows and columns to account for
! // column-major storage.
! octave_fftw::ifft2d (retval.fortran_vec (), nc, nr);
return retval;
}
--- 1179,1191 ----
ComplexMatrix
ComplexMatrix::ifourier2d (void) const
{
! dim_vector dv(rows (), cols ());
!
! ComplexMatrix retval (rows (), cols ());
! const Complex *in (data ());
! Complex *out (retval.fortran_vec ());
! octave_fftw::ifftNd (in, out, 2, dv);
return retval;
}
*** ./liboctave/oct-fftw.h.orig 2004-01-28 15:17:17.000000000 +0100
--- ./liboctave/oct-fftw.h 2004-02-10 11:39:44.000000000 +0100
***************
*** 22,45 ****
#define octave_oct_fftw_h 1
#include
!
! #if defined (HAVE_DFFTW_H)
! #include
! #else
! #include
! #endif
#include "oct-cmplx.h"
class
octave_fftw
{
public:
! static int fft (const Complex*, Complex *, size_t);
! static int ifft (const Complex*, Complex *, size_t);
!
! static int fft2d (Complex*, size_t, size_t);
! static int ifft2d (Complex*, size_t, size_t);
private:
octave_fftw ();
--- 22,48 ----
#define octave_oct_fftw_h 1
#include
! #include
#include "oct-cmplx.h"
+ #include "dim-vector.h"
class
octave_fftw
{
public:
! static int fft (const double *in, Complex *out, size_t npts,
! size_t nsamples = 1, int stride = 1, int dist = -1);
! static int fft (const Complex *in, Complex *out, size_t npts,
! size_t nsamples = 1, int stride = 1, int dist = -1);
! static int ifft (const Complex *in, Complex *out, size_t npts,
! size_t nsamples = 1, int stride = 1, int 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 &);
private:
octave_fftw ();
*** ./liboctave/oct-fftw.cc.orig 2004-01-28 15:17:22.000000000 +0100
--- ./liboctave/oct-fftw.cc 2004-02-12 15:23:35.000000000 +0100
***************
*** 22,39 ****
#include
#endif
! #ifdef HAVE_FFTW
#include "oct-fftw.h"
#include "lo-error.h"
!
// Helper class to create and cache fftw plans for both 1d and 2d. This
// implementation uses FFTW_ESTIMATE to create the plans, which in theory
! // is suboptimal, but provides quite reasonable performance. Future
! // enhancement will be to add a dynamically loadable interface ("fftw")
! // to manipulate fftw wisdom so that users may choose the appropriate
! // planner.
class
octave_fftw_planner
--- 22,47 ----
#include
#endif
! #if defined (HAVE_FFTW3)
#include "oct-fftw.h"
#include "lo-error.h"
! #include
// Helper class to create and cache fftw plans for both 1d and 2d. This
// implementation uses FFTW_ESTIMATE to create the plans, which in theory
! // is suboptimal, but provides quite reasonable performance.
!
! // Also note that if FFTW_ESTIMATE is not used the planner in FFTW3
! // destroys the input and output arrays. So with the form of the
! // current code we definitely want FFTW_ESTIMATE!! However, we use
! // any wsidom that is available, either in a FFTW3 system wide file
! // or as supplied by the user.
!
! // XXX FIXME XXX If we can ensure 16 byte alignment in Array ( *data)
! // the FFTW3 can use SIMD instructions for further acceleration.
!
! // Note that it is profitable to store the FFTW3 plans, for small ffts
class
octave_fftw_planner
***************
*** 41,57 ****
public:
octave_fftw_planner ();
! fftw_plan create_plan (fftw_direction, size_t);
! fftwnd_plan create_plan2d (fftw_direction, size_t, size_t);
private:
int plan_flags;
fftw_plan plan[2];
! fftwnd_plan plan2d[2];
!
! size_t n[2];
! size_t n2d[2][2];
};
octave_fftw_planner::octave_fftw_planner ()
--- 49,83 ----
public:
octave_fftw_planner ();
! fftw_plan create_plan (int dir, const int rank, const dim_vector dims,
! int howmany, int stride, int dist,
! const Complex *in, Complex *out);
! fftw_plan create_plan (const int rank, const dim_vector dims,
! int howmany, int stride, int dist,
! const double *in, Complex *out);
private:
int plan_flags;
+ // Plan for fft and ifft of complex values
fftw_plan plan[2];
! int d[2]; // dist
! int s[2]; // stride
! int r[2]; // rank
! int h[2]; // howmany
! dim_vector n[2]; // dims
! char ialign[2];
! char oalign[2];
!
! // Plan for fft of real values
! fftw_plan rplan;
! int rd; // dist
! int rs; // stride
! int rr; // rank
! int rh; // howmany
! dim_vector rn; // dims
! char rialign;
! char roalign;
};
octave_fftw_planner::octave_fftw_planner ()
***************
*** 59,89 ****
plan_flags = FFTW_ESTIMATE;
plan[0] = plan[1] = 0;
! plan2d[0] = plan2d[1] = 0;
!
! n[0] = n[1] = 0;
! n2d[0][0] = n2d[0][1] = n2d[1][0] = n2d[1][1] = 0;
}
fftw_plan
! octave_fftw_planner::create_plan (fftw_direction dir, size_t npts)
{
! size_t which = (dir == FFTW_FORWARD) ? 0 : 1;
fftw_plan *cur_plan_p = &plan[which];
bool create_new_plan = false;
! if (plan[which] == 0 || n[which] != npts)
! {
! create_new_plan = true;
! n[which] = npts;
! }
if (create_new_plan)
{
if (*cur_plan_p)
fftw_destroy_plan (*cur_plan_p);
! *cur_plan_p = fftw_create_plan (npts, dir, plan_flags);
if (*cur_plan_p == 0)
(*current_liboctave_error_handler) ("Error creating fftw plan");
--- 85,151 ----
plan_flags = FFTW_ESTIMATE;
plan[0] = plan[1] = 0;
! d[0] = d[1] = s[0] = s[1] = r[0] = r[1] = h[0] = h[1] = 0;
! ialign[0] = ialign[1] = oalign[0] = oalign[1] = 0;
! n[0] = n[1] = dim_vector();
!
! rplan = 0;
! rd = rs = rr = rh = 0;
! rialign = roalign = 0;
! rn = dim_vector ();
!
! // If we have a system wide wisdom file import it
! fftw_import_system_wisdom ( );
}
fftw_plan
! octave_fftw_planner::create_plan (int dir, const int rank,
! const dim_vector dims, int howmany,
! int stride, int dist,
! const Complex *in, Complex *out)
{
! int which = (dir == FFTW_FORWARD) ? 0 : 1;
fftw_plan *cur_plan_p = &plan[which];
bool create_new_plan = false;
+ char in_align = ((int) in) & 0xF;
+ char out_align = ((int) out) & 0xF;
! if (plan[which] == 0 || d[which] != dist || s[which] != stride ||
! r[which] != rank || h[which] != howmany ||
! ialign[which] != in_align || oalign[which] != out_align)
! 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) != n[which](i))
! {
! create_new_plan = true;
! break;
! }
if (create_new_plan)
{
+ d[which] = dist;
+ s[which] = stride;
+ r[which] = rank;
+ h[which] = howmany;
+ ialign[which] = in_align;
+ oalign[which] = out_align;
+ n[which] = dims;
+
if (*cur_plan_p)
fftw_destroy_plan (*cur_plan_p);
! // 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);
!
! *cur_plan_p =
! fftw_plan_many_dft (rank, tmp, howmany,
! reinterpret_cast (const_cast (in)),
! NULL, stride, dist, reinterpret_cast (out),
! NULL, stride, dist, dir, plan_flags);
if (*cur_plan_p == 0)
(*current_liboctave_error_handler) ("Error creating fftw plan");
***************
*** 92,123 ****
return *cur_plan_p;
}
! fftwnd_plan
! octave_fftw_planner::create_plan2d (fftw_direction dir,
! size_t nrows, size_t ncols)
{
! size_t which = (dir == FFTW_FORWARD) ? 0 : 1;
! fftwnd_plan *cur_plan_p = &plan2d[which];
bool create_new_plan = false;
! if (plan2d[which] == 0 || n2d[which][0] != nrows || n2d[which][1] != ncols)
! {
! create_new_plan = true;
!
! n2d[which][0] = nrows;
! n2d[which][1] = ncols;
! }
if (create_new_plan)
{
if (*cur_plan_p)
! fftwnd_destroy_plan (*cur_plan_p);
! *cur_plan_p = fftw2d_create_plan (nrows, ncols, dir,
! plan_flags | FFTW_IN_PLACE);
if (*cur_plan_p == 0)
! (*current_liboctave_error_handler) ("Error creating 2d fftw plan");
}
return *cur_plan_p;
--- 154,208 ----
return *cur_plan_p;
}
! fftw_plan
! octave_fftw_planner::create_plan (const int rank, const dim_vector dims,
! int howmany, int stride, int dist,
! const double *in, Complex *out)
{
! fftw_plan *cur_plan_p = &rplan;
bool create_new_plan = false;
+ char in_align = ((int) in) & 0xF;
+ char out_align = ((int) out) & 0xF;
! if (rplan == 0 || rd != dist || rs != stride ||
! rr != rank || rh != howmany ||
! rialign != in_align || roalign != out_align)
! 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) != rn(i))
! {
! create_new_plan = true;
! break;
! }
if (create_new_plan)
{
+ rd = dist;
+ rs = stride;
+ rr = rank;
+ rh = howmany;
+ rialign = in_align;
+ roalign = out_align;
+ rn = dims;
+
if (*cur_plan_p)
! fftw_destroy_plan (*cur_plan_p);
! // 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);
!
! *cur_plan_p =
! fftw_plan_many_dft_r2c (rank, tmp, howmany,
! (const_cast (in)),
! NULL, stride, dist, reinterpret_cast (out),
! NULL, stride, dist, plan_flags);
if (*cur_plan_p == 0)
! (*current_liboctave_error_handler) ("Error creating fftw plan");
}
return *cur_plan_p;
***************
*** 125,175 ****
static octave_fftw_planner fftw_planner;
int
! octave_fftw::fft (const Complex *in, Complex *out, size_t npts)
{
! fftw_one (fftw_planner.create_plan (FFTW_FORWARD, npts),
! reinterpret_cast (const_cast (in)),
! reinterpret_cast (out));
return 0;
}
int
! octave_fftw::ifft (const Complex *in, Complex *out, size_t npts)
{
! fftw_one (fftw_planner.create_plan (FFTW_BACKWARD, npts),
! reinterpret_cast (const_cast (in)),
! reinterpret_cast (out));
const Complex scale = npts;
! for (size_t i = 0; i < npts; i++)
! out[i] /= scale;
return 0;
}
int
! octave_fftw::fft2d (Complex *inout, size_t nr, size_t nc)
{
! fftwnd_one (fftw_planner.create_plan2d (FFTW_FORWARD, nr, nc),
! reinterpret_cast (inout),
! 0);
return 0;
}
int
! octave_fftw::ifft2d (Complex *inout, size_t nr, size_t nc)
{
! fftwnd_one (fftw_planner.create_plan2d (FFTW_BACKWARD, nr, nc),
! reinterpret_cast (inout),
! 0);
! const size_t npts = nr * nc;
const Complex scale = npts;
for (size_t i = 0; i < npts; i++)
! inout[i] /= scale;
return 0;
}
--- 210,370 ----
static octave_fftw_planner fftw_planner;
+ static inline void convert_packcomplex_1d (Complex *out, size_t nr,
+ size_t nc, int stride, int dist)
+ {
+ // Fill in the missing data
+ for (size_t i = 0; i < nr; i++)
+ for (size_t j = nc/2+1; j < nc; j++)
+ out[j*stride + i*dist] = conj(out[(nc - j)*stride + i*dist]);
+ }
+
+ static inline void convert_packcomplex_Nd (Complex *out,
+ const dim_vector &dv)
+ {
+ size_t nc = dv(0);
+ size_t nr = dv(1);
+ size_t np = (dv.length() > 2 ? dv.numel () / nc / nr : 1);
+ size_t nrp = nr * np;
+ Complex *ptr1, *ptr2;
+
+ // Create space for the missing elements
+ for (size_t i = 0; i < nrp; i++)
+ {
+ ptr1 = out + i * (nc/2 + 1) + nrp*((nc-1)/2);
+ ptr2 = out + i * nc;
+ for (size_t j = 0; j < nc/2+1; j++)
+ *ptr2++ = *ptr1++;
+ }
+
+ // Fill in the missing data
+ for (size_t i = 0; i < np; i++)
+ for (size_t j = 1; j < nr; j++)
+ for (size_t k = nc/2+1; k < nc; k++)
+ out[k + (j + i*nr)*nc] = conj(out[nc - k + ((i+1)*nr - j)*nc]);
+ for (size_t j = nc/2+1; j < nc; j++)
+ out[j] = conj(out[nc - j]);
+ }
+
int
! octave_fftw::fft (const double *in, Complex *out, size_t npts,
! size_t nsamples, int stride, int dist)
{
! dist = (dist < 0 ? npts : dist);
!
! dim_vector dv (npts);
! fftw_plan plan = fftw_planner.create_plan (1, dv, nsamples, stride, dist,
! in, out);
!
! fftw_execute_dft_r2c (plan, (const_cast(in)),
! reinterpret_cast (out));
!
! // Need to create other half of the transform
! convert_packcomplex_1d (out, nsamples, npts, stride, dist);
return 0;
}
int
! octave_fftw::fft (const Complex *in, Complex *out, size_t npts,
! size_t nsamples, int stride, int dist)
{
! dist = (dist < 0 ? npts : dist);
!
! dim_vector dv (npts);
! fftw_plan plan = fftw_planner.create_plan (FFTW_FORWARD, 1, dv, nsamples,
! stride, dist, in, out);
!
! fftw_execute_dft (plan,
! reinterpret_cast (const_cast(in)),
! reinterpret_cast (out));
!
! return 0;
! }
!
! int
! octave_fftw::ifft (const Complex *in, Complex *out, size_t npts,
! size_t nsamples, int stride, int dist)
! {
! dist = (dist < 0 ? npts : dist);
!
! dim_vector dv (npts);
! fftw_plan plan = fftw_planner.create_plan (FFTW_BACKWARD, 1, dv, nsamples,
! stride, dist, in, out);
!
! fftw_execute_dft (plan,
! reinterpret_cast (const_cast(in)),
! reinterpret_cast (out));
const Complex scale = npts;
! for (size_t j = 0; j < nsamples; j++)
! for (size_t i = 0; i < npts; i++)
! out[i*stride + j*dist] /= scale;
!
! return 0;
! }
!
! int
! octave_fftw::fftNd (const double *in, Complex *out, const int rank,
! const dim_vector &dv)
! {
! int dist = 1;
! for (int i = 0; i < rank; i++)
! dist *= dv(i);
!
! // Fool with the position of the start of the output matrix, so that
! // creating other half of the matrix won't cause cache problems
! int offset = (dv.numel () / dv(0)) * ((dv(0) - 1) / 2);
!
! fftw_plan plan = fftw_planner.create_plan (rank, dv, 1, 1, dist,
! in, out + offset);
!
! fftw_execute_dft_r2c (plan, (const_cast(in)),
! reinterpret_cast (out+ offset));
!
! // Need to create other half of the transform
! convert_packcomplex_Nd (out, dv);
return 0;
}
int
! octave_fftw::fftNd (const Complex *in, Complex *out, const int rank,
! const dim_vector &dv)
{
! int dist = 1;
! for (int i = 0; i < rank; i++)
! dist *= dv(i);
!
! fftw_plan plan = fftw_planner.create_plan (FFTW_FORWARD, rank, dv, 1, 1,
! dist, in, out);
!
! fftw_execute_dft (plan,
! reinterpret_cast (const_cast(in)),
! reinterpret_cast (out));
return 0;
}
int
! octave_fftw::ifftNd (const Complex *in, Complex *out, const int rank,
! const dim_vector &dv)
{
! int dist = 1;
! for (int i = 0; i < rank; i++)
! dist *= dv(i);
!
! fftw_plan plan = fftw_planner.create_plan (FFTW_BACKWARD, rank, dv, 1, 1,
! dist, in, out);
!
! fftw_execute_dft (plan,
! reinterpret_cast (const_cast(in)),
! reinterpret_cast (out));
! const size_t npts = dv.numel ();
const Complex scale = npts;
for (size_t i = 0; i < npts; i++)
! out[i] /= scale;
return 0;
}
*** ./liboctave/CNDArray.cc.orig 2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/CNDArray.cc 2004-02-10 18:00:10.000000000 +0100
***************
*** 37,42 ****
--- 37,63 ----
#include "lo-ieee.h"
#include "lo-mappers.h"
+ #if defined (HAVE_FFTW3)
+ # include "oct-fftw.h"
+ #else
+ extern "C"
+ {
+ // Note that the original complex fft routines were not written for
+ // double complex arguments. They have been modified by adding an
+ // implicit double precision (a-h,o-z) statement at the beginning of
+ // each subroutine.
+
+ F77_RET_T
+ F77_FUNC (cffti, CFFTI) (const int&, Complex*);
+
+ F77_RET_T
+ F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
+
+ F77_RET_T
+ F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
+ }
+ #endif
+
// XXX FIXME XXX -- could we use a templated mixed-type copy function
// here?
***************
*** 61,66 ****
--- 82,505 ----
elem (i) = a.elem (i);
}
+ #if defined (HAVE_FFTW3)
+ ComplexNDArray
+ ComplexNDArray::fourier (int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ int stride = 1;
+ int n = dv(dim);
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / dv (dim);
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
+ int dist = (stride == 1 ? n : 1);
+
+ const Complex *in (fortran_vec ());
+ ComplexNDArray retval (dv);
+ Complex *out (retval.fortran_vec ());
+
+ // Need to be careful here about the distance between fft's
+ for (int k = 0; k < nloop; k++)
+ octave_fftw::fft (in + k * stride * n, out + k * stride * n,
+ n, howmany, stride, dist);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::ifourier (const int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ int stride = 1;
+ int n = dv(dim);
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / dv (dim);
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
+ int dist = (stride == 1 ? n : 1);
+
+ const Complex *in (fortran_vec ());
+ ComplexNDArray retval (dv);
+ Complex *out (retval.fortran_vec ());
+
+ // Need to be careful here about the distance between fft's
+ for (int k = 0; k < nloop; k++)
+ octave_fftw::ifft (in + k * stride * n, out + k * stride * n,
+ n, howmany, stride, dist);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::fourier2d (void) const
+ {
+ dim_vector dv = dims();
+ if (dv.length () < 2)
+ return ComplexNDArray ();
+
+ dim_vector dv2(dv(0), dv(1));
+ const Complex *in = fortran_vec ();
+ ComplexNDArray retval (dv);
+ Complex *out = retval.fortran_vec ();
+ int howmany = numel() / dv(0) / dv(1);
+ int dist = dv(0) * dv(1);
+
+ for (int i=0; i < howmany; i++)
+ octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::ifourier2d (void) const
+ {
+ dim_vector dv = dims();
+ if (dv.length () < 2)
+ return ComplexNDArray ();
+
+ dim_vector dv2(dv(0), dv(1));
+ const Complex *in = fortran_vec ();
+ ComplexNDArray retval (dv);
+ Complex *out = retval.fortran_vec ();
+ int howmany = numel() / dv(0) / dv(1);
+ int dist = dv(0) * dv(1);
+
+ for (int i=0; i < howmany; i++)
+ octave_fftw::ifftNd (in + i*dist, out + i*dist, 2, dv2);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::fourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+
+ const Complex *in (fortran_vec ());
+ ComplexNDArray retval (dv);
+ Complex *out (retval.fortran_vec ());
+
+ octave_fftw::fftNd (in, out, rank, dv);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::ifourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+
+ const Complex *in (fortran_vec ());
+ ComplexNDArray retval (dv);
+ Complex *out (retval.fortran_vec ());
+
+ octave_fftw::ifftNd (in, out, rank, dv);
+
+ return retval;
+ }
+
+ #else
+ ComplexNDArray
+ ComplexNDArray::fourier (int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ ComplexNDArray retval (dv);
+ int npts = dv(dim);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
+
+ int stride = 1;
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int i = 0; i < npts; i++)
+ tmp[i] = elem((i + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftf, CFFTF) (npts, &tmp, pwsave);
+
+ for (int i = 0; i < npts; i++)
+ retval ((i + k*npts)*stride + j*dist) = tmp[i];
+ }
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::ifourier (int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ ComplexNDArray retval (dv);
+ int npts = dv(dim);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
+
+ int stride = 1;
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int i = 0; i < npts; i++)
+ tmp[i] = elem((i + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftb, CFFTB) (npts, &tmp, pwsave);
+
+ for (int i = 0; i < npts; i++)
+ retval ((i + k*npts)*stride + j*dist) = tmp[i];
+ }
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::fourier2d (void) const
+ {
+ dim_vector dv (dims(0), dims(1));
+ int rank = 2;
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::ifourier2d (void) const
+ {
+ dim_vector dv (dims(0), dims(1));
+ int rank = 2;
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::fourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ ComplexNDArray::ifourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ #endif
+
// unary operations
boolNDArray
*** ./liboctave/CNDArray.h.orig 2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/CNDArray.h 2004-02-10 17:24:25.000000000 +0100
***************
*** 88,93 ****
--- 88,102 ----
ComplexNDArray sum (int dim = -1) const;
ComplexNDArray sumsq (int dim = -1) const;
+ ComplexNDArray fourier (int dim = 1) const;
+ ComplexNDArray ifourier (int dim = 1) const;
+
+ ComplexNDArray fourier2d (void) const;
+ ComplexNDArray ifourier2d (void) const;
+
+ ComplexNDArray fourierNd (void) const;
+ ComplexNDArray ifourierNd (void) const;
+
NDArray abs (void) const;
ComplexMatrix matrix_value (void) const;
*** ./liboctave/dNDArray.cc.orig 2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/dNDArray.cc 2004-02-11 13:42:48.000000000 +0100
***************
*** 38,43 ****
--- 38,480 ----
#include "lo-ieee.h"
#include "lo-mappers.h"
+ #if defined (HAVE_FFTW3)
+ #include "oct-fftw.h"
+
+ ComplexNDArray
+ NDArray::fourier (int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ int stride = 1;
+ int n = dv(dim);
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / dv (dim);
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
+ int dist = (stride == 1 ? n : 1);
+
+ const double *in (fortran_vec ());
+ ComplexNDArray retval (dv);
+ Complex *out (retval.fortran_vec ());
+
+ // Need to be careful here about the distance between fft's
+ for (int k = 0; k < nloop; k++)
+ octave_fftw::fft (in + k * stride * n, out + k * stride * n,
+ n, howmany, stride, dist);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::ifourier (const int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ int stride = 1;
+ int n = dv(dim);
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / dv (dim);
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
+ int dist = (stride == 1 ? n : 1);
+
+ ComplexNDArray retval (*this);
+ Complex *out (retval.fortran_vec ());
+
+ // Need to be careful here about the distance between fft's
+ for (int k = 0; k < nloop; k++)
+ octave_fftw::ifft (out + k * stride * n, out + k * stride * n,
+ n, howmany, stride, dist);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::fourier2d (void) const
+ {
+ dim_vector dv = dims();
+ if (dv.length () < 2)
+ return ComplexNDArray ();
+
+ dim_vector dv2(dv(0), dv(1));
+ const double *in = fortran_vec ();
+ ComplexNDArray retval (dv);
+ Complex *out = retval.fortran_vec ();
+ int howmany = numel() / dv(0) / dv(1);
+ int dist = dv(0) * dv(1);
+
+ for (int i=0; i < howmany; i++)
+ octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::ifourier2d (void) const
+ {
+ dim_vector dv = dims();
+ if (dv.length () < 2)
+ return ComplexNDArray ();
+
+ dim_vector dv2(dv(0), dv(1));
+ ComplexNDArray retval (*this);
+ Complex *out = retval.fortran_vec ();
+ int howmany = numel() / dv(0) / dv(1);
+ int dist = dv(0) * dv(1);
+
+ for (int i=0; i < howmany; i++)
+ octave_fftw::ifftNd (out + i*dist, out + i*dist, 2, dv2);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::fourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+
+ const double *in (fortran_vec ());
+ ComplexNDArray retval (dv);
+ Complex *out (retval.fortran_vec ());
+
+ octave_fftw::fftNd (in, out, rank, dv);
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::ifourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+
+ ComplexNDArray tmp (*this);
+ Complex *in (tmp.fortran_vec ());
+ ComplexNDArray retval (dv);
+ Complex *out (retval.fortran_vec ());
+
+ octave_fftw::ifftNd (in, out, rank, dv);
+
+ return retval;
+ }
+
+ #else
+
+ extern "C"
+ {
+ // Note that the original complex fft routines were not written for
+ // double complex arguments. They have been modified by adding an
+ // implicit double precision (a-h,o-z) statement at the beginning of
+ // each subroutine.
+
+ F77_RET_T
+ F77_FUNC (cffti, CFFTI) (const int&, Complex*);
+
+ F77_RET_T
+ F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
+
+ F77_RET_T
+ F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
+ }
+
+ ComplexNDArray
+ NDArray::fourier (int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ ComplexNDArray retval (dv);
+ int npts = dv(dim);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
+
+ int stride = 1;
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int i = 0; i < npts; i++)
+ tmp[i] = elem((i + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftf, CFFTF) (npts, &tmp, pwsave);
+
+ for (int i = 0; i < npts; i++)
+ retval ((i + k*npts)*stride + j*dist) = tmp[i];
+ }
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::ifourier (int dim) const
+ {
+ dim_vector dv = dims ();
+
+ if (dim > dv.length () || dim < 0)
+ return ComplexNDArray ();
+
+ ComplexNDArray retval (dv);
+ int npts = dv(dim);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
+
+ int stride = 1;
+
+ for (int i = 0; i < dim; i++)
+ stride *= dv(i);
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int i = 0; i < npts; i++)
+ tmp[i] = elem((i + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftb, CFFTB) (npts, &tmp, pwsave);
+
+ for (int i = 0; i < npts; i++)
+ retval ((i + k*npts)*stride + j*dist) = tmp[i];
+ }
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::fourier2d (void) const
+ {
+ dim_vector dv (dims(0), dims(1));
+ int rank = 2;
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::ifourier2d (void) const
+ {
+ dim_vector dv (dims(0), dims(1));
+ int rank = 2;
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::fourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ ComplexNDArray
+ NDArray::ifourierNd (void) const
+ {
+ dim_vector dv = dims ();
+ int rank = dv.length ();
+ ComplexNDArray retval (*this);
+ int stride = 1;
+
+ for (int i = 0; i < rank; i++)
+ {
+ int npts = dv(i);
+ int nn = 4*npts+15;
+ Array wsave (nn);
+ Complex *pwsave = wsave.fortran_vec ();
+ Array row (npts);
+ Complex *prow = row.fortran_vec ();
+
+ int howmany = numel () / npts;
+ howmany = (stride == 1 ? howmany :
+ (howmany > stride ? stride : howmany));
+ int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+ int dist = (stride == 1 ? npts : 1);
+
+ F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+ for (int k = 0; k < nloop; k++)
+ {
+ for (int j = 0; j < howmany; j++)
+ {
+ OCTAVE_QUIT;
+
+ for (int l = 0; l < npts; l++)
+ prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+ F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
+
+ for (int l = 0; l < npts; l++)
+ retval ((l + k*npts)*stride + j*dist) = prow[l];
+ }
+ }
+
+ stride *= dv(i);
+ }
+
+ return retval;
+ }
+
+ #endif
+
NDArray::NDArray (const boolNDArray& a)
: MArrayN (a.dims ())
{
*** ./liboctave/dNDArray.h.orig 2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/dNDArray.h 2004-02-10 17:24:16.000000000 +0100
***************
*** 87,92 ****
--- 87,101 ----
NDArray sum (int dim = -1) const;
NDArray sumsq (int dim = -1) const;
+ ComplexNDArray fourier (int dim = 1) const;
+ ComplexNDArray ifourier (int dim = 1) const;
+
+ ComplexNDArray fourier2d (void) const;
+ ComplexNDArray ifourier2d (void) const;
+
+ ComplexNDArray fourierNd (void) const;
+ ComplexNDArray ifourierNd (void) const;
+
NDArray abs (void) const;
friend NDArray real (const ComplexNDArray& a);
*** ./src/DLD-FUNCTIONS/ifft.cc.orig 2004-01-30 22:46:35.000000000 +0100
--- ./src/DLD-FUNCTIONS/ifft.cc 2004-02-11 14:08:36.000000000 +0100
***************
*** 1,130 ****
- /*
-
- Copyright (C) 1996, 1997 John W. Eaton
-
- This file is part of Octave.
-
- Octave is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation; either version 2, or (at your option) any
- later version.
-
- Octave is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with Octave; see the file COPYING. If not, write to the Free
- Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-
- */
-
- #ifdef HAVE_CONFIG_H
- #include
- #endif
-
- #include "lo-mappers.h"
-
- #include "defun-dld.h"
- #include "error.h"
- #include "gripes.h"
- #include "oct-obj.h"
- #include "utils.h"
-
- // This function should be merged with Ffft.
-
- DEFUN_DLD (ifft, args, ,
- "-*- texinfo -*-\n\
- @deftypefn {Loadable Function} {} ifft (@var{a}, @var{n})\n\
- Compute the inverse FFT of @var{a} using subroutines from @sc{Fftpack}. If\n\
- @var{a} is a matrix, @code{fft} computes the inverse FFT for each column\n\
- of @var{a}.\n\
- \n\
- If called with two arguments, @var{n} is expected to be an integer\n\
- specifying the number of elements of @var{a} to use. If @var{a} is a\n\
- matrix, @var{n} specifies the number of rows of @var{a} to use. If\n\
- @var{n} is larger than the size of @var{a}, @var{a} is resized and\n\
- padded with zeros.\n\
- @end deftypefn")
- {
- octave_value retval;
-
- int nargin = args.length ();
-
- if (nargin < 1 || nargin > 2)
- {
- print_usage ("ifft");
- return retval;
- }
-
- octave_value arg = args(0);
-
- int n_points = arg.rows ();
- if (n_points == 1)
- n_points = arg.columns ();
-
- if (nargin == 2)
- {
- double dval = args(1).double_value ();
- if (xisnan (dval))
- error ("fft: NaN is invalid as the N_POINTS");
- else
- n_points = NINT (dval);
- }
-
- 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 octave_value (Matrix ());
-
- if (arg.is_real_type ())
- {
- Matrix m = arg.matrix_value ();
-
- if (! error_state)
- {
- if (m.rows () == 1)
- m.resize (1, n_points, 0.0);
- else
- m.resize (n_points, m.columns (), 0.0);
- retval = m.ifourier ();
- }
- }
- else if (arg.is_complex_type ())
- {
- ComplexMatrix m = arg.complex_matrix_value ();
-
- if (! error_state)
- {
- if (m.rows () == 1)
- m.resize (1, n_points, 0.0);
- else
- m.resize (n_points, m.columns (), 0.0);
- retval = m.ifourier ();
- }
- }
- else
- {
- gripe_wrong_type_arg ("ifft", arg);
- }
-
- return retval;
- }
-
- /*
- ;;; Local Variables: ***
- ;;; mode: C++ ***
- ;;; End: ***
- */
--- 0 ----
*** ./src/DLD-FUNCTIONS/fftn.cc.orig 2004-02-11 14:09:07.000000000 +0100
--- ./src/DLD-FUNCTIONS/fftn.cc 2004-02-11 14:12:44.000000000 +0100
***************
*** 0 ****
--- 1,165 ----
+ /*
+
+ Copyright (C) 2004 David Bateman
+
+ This file is part of Octave.
+
+ Octave is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2, or (at your option) any
+ later version.
+
+ Octave is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Octave; see the file COPYING. If not, write to the Free
+ Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+ */
+
+ #ifdef HAVE_CONFIG_H
+ #include
+ #endif
+
+ #include "lo-mappers.h"
+
+ #include "defun-dld.h"
+ #include "error.h"
+ #include "gripes.h"
+ #include "oct-obj.h"
+ #include "utils.h"
+
+ // This function should be merged with Fifft.
+
+ #if defined (HAVE_FFTW3)
+ #define FFTSRC "@sc{Fftw}"
+ #define WISDOM ", fft_wisdom"
+ #else
+ #define FFTSRC "@sc{Fftpack}"
+ #define WISDOM
+ #endif
+
+ static octave_value
+ do_fftn (const octave_value_list &args, const char *fcn, int type)
+ {
+ octave_value retval;
+
+ int nargin = args.length ();
+
+ if (nargin < 1 || nargin > 2)
+ {
+ print_usage (fcn);
+ return retval;
+ }
+
+ octave_value arg = args(0);
+ dim_vector dims = arg.dims ();
+
+ for (int i = 0; i < dims.length (); i++)
+ if (dims(i) < 0)
+ return retval;
+
+ if (nargin > 1)
+ {
+ Matrix val = args(1).matrix_value ();
+ if (val.rows () > val.columns ())
+ val.transpose ();
+
+ if (val.columns () != dims.length ())
+ error ("%s: invalid size of second argument", fcn);
+ else
+ {
+ for (int i = 0; i < dims.length (); i++)
+ {
+ if (xisnan (val(i,0)))
+ error ("%s: NaN is invalid as a dimension", fcn);
+ else if (NINT (val(i,0)) < 0)
+ error ("%s: all dimension must be greater than zero", fcn);
+ else
+ {
+ dims(i) = NINT(val(i,0));
+ }
+ }
+ }
+ }
+
+ if (error_state)
+ return retval;
+
+ if (dims.all_zero ())
+ return octave_value (Matrix ());
+
+ if (arg.is_real_type ())
+ {
+ NDArray nda = arg.array_value ();
+
+ if (! error_state)
+ {
+ nda.resize (dims, 0.0);
+ retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
+ }
+ }
+ else if (arg.is_complex_type ())
+ {
+ ComplexNDArray cnda = arg.complex_matrix_value ();
+
+ if (! error_state)
+ {
+ cnda.resize (dims, 0.0);
+ retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
+ }
+ }
+ else
+ {
+ gripe_wrong_type_arg (fcn, arg);
+ }
+
+ return retval;
+ }
+
+ DEFUN_DLD (fftn, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} fftn (@var{a}, @var{siz})\n\
+ Compute the N dimensional FFT of @var{a} using subroutines from\n"
+ FFTSRC
+ ". The optional vector argument @var{siz} may be used specify the\n\
+ dimensions of the array to be used. If an element of @var{siz} is\n\
+ smaller than the corresponding dimension, then the dimension is\n\
+ truncated prior to performing the FFT. Otherwise if an element\n\
+ of @var{siz} is larger than the corresponding dimension @var{a}\n\
+ is resized and padded with zeros.\n\
+ @end deftypefn\n\
+ @seealso {ifftn, fft, fft2"
+ WISDOM
+ "}")
+ {
+ return do_fftn(args, "fftn", 0);
+ }
+
+ DEFUN_DLD (ifftn, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} ifftn (@var{a}, @var{siz})\n\
+ Compute the invesre N dimensional FFT of @var{a} using subroutines from\n"
+ FFTSRC
+ ". The optional vector argument @var{siz} may be used specify the\n\
+ dimensions of the array to be used. If an element of @var{siz} is\n\
+ smaller than the corresponding dimension, then the dimension is\n\
+ truncated prior to performing the inverse FFT. Otherwise if an element\n\
+ of @var{siz} is larger than the corresponding dimension @var{a}\n\
+ is resized and padded with zeros.\n\
+ @end deftypefn\n\
+ @seealso {fftn, ifft, ifft2"
+ WISDOM
+ "}")
+ {
+ return do_fftn(args, "ifftn", 1);
+ }
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C++ ***
+ ;;; End: ***
+ */
*** ./src/DLD-FUNCTIONS/fft.cc.orig 2004-01-30 22:46:06.000000000 +0100
--- ./src/DLD-FUNCTIONS/fft.cc 2004-02-11 14:12:22.000000000 +0100
***************
*** 1,5 ****
--- 1,6 ----
/*
+ Copyright (C) 2004 David Bateman
Copyright (C) 1996, 1997 John W. Eaton
This file is part of Octave.
***************
*** 32,127 ****
#include "oct-obj.h"
#include "utils.h"
! // This function should be merged with Fifft.
! DEFUN_DLD (fft, args, ,
! "-*- texinfo -*-\n\
! @deftypefn {Loadable Function} {} fft (@var{a}, @var{n})\n\
! Compute the FFT of @var{a} using subroutines from @sc{Fftpack}. If @var{a}\n\
! is a matrix, @code{fft} computes the FFT for each column of @var{a}.\n\
! \n\
! If called with two arguments, @var{n} is expected to be an integer\n\
! specifying the number of elements of @var{a} to use. If @var{a} is a\n\
! matrix, @var{n} specifies the number of rows of @var{a} to use. If\n\
! @var{n} is larger than the size of @var{a}, @var{a} is resized and\n\
! padded with zeros.\n\
! @end deftypefn")
{
octave_value retval;
int nargin = args.length ();
! if (nargin < 1 || nargin > 2)
{
! print_usage ("fft");
return retval;
}
octave_value arg = args(0);
! int n_points = arg.rows ();
! if (n_points == 1)
! n_points = arg.columns ();
! if (nargin == 2)
{
! double dval = args(1).double_value ();
if (xisnan (dval))
! error ("fft: NaN is invalid as the N_POINTS");
else
! n_points = NINT (dval);
}
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 octave_value (Matrix ());
if (arg.is_real_type ())
{
! Matrix m = arg.matrix_value ();
if (! error_state)
{
! if (m.rows () == 1)
! m.resize (1, n_points, 0.0);
! else
! m.resize (n_points, m.columns (), 0.0);
! retval = m.fourier ();
}
}
else if (arg.is_complex_type ())
{
! ComplexMatrix m = arg.complex_matrix_value ();
if (! error_state)
{
! if (m.rows () == 1)
! m.resize (1, n_points, 0.0);
! else
! m.resize (n_points, m.columns (), 0.0);
! retval = m.fourier ();
}
}
else
{
! gripe_wrong_type_arg ("fft", arg);
}
return retval;
}
/*
;;; Local Variables: ***
;;; mode: C++ ***
--- 33,201 ----
#include "oct-obj.h"
#include "utils.h"
! #if defined (HAVE_FFTW3)
! #define FFTSRC "@sc{Fftw}"
! #define WISDOM ", fft_wisdom"
! #else
! #define FFTSRC "@sc{Fftpack}"
! #define WISDOM
! #endif
! static octave_value
! do_fft (const octave_value_list &args, const char *fcn, int type)
{
octave_value retval;
int nargin = args.length ();
! if (nargin < 1 || nargin > 3)
{
! print_usage (fcn);
return retval;
}
octave_value arg = args(0);
+ dim_vector dims = arg.dims ();
+ int n_points = -1;
+ int dim = -1;
+
+ if (nargin > 1)
+ {
+ if (! args(1).is_empty ())
+ {
+ double dval = args(1).double_value ();
+ if (xisnan (dval))
+ error ("%s: NaN is invalid as the N_POINTS", fcn);
+ else
+ {
+ n_points = NINT (dval);
+ if (n_points < 0)
+ error ("%s: number of points must be greater than zero", fcn);
+ }
+ }
+ }
! if (error_state)
! return retval;
! if (nargin > 2)
{
! double dval = args(2).double_value ();
if (xisnan (dval))
! error ("%s: NaN is invalid as the N_POINTS", fcn);
! else if (dval < 1 || dval > dims.length ())
! error ("%s: invalid dimension along which to perform fft", fcn);
else
! dim = NINT (dval) - 1;
}
if (error_state)
return retval;
! for (int i = 0; i < dims.length (); i++)
! if (dims(i) < 0)
return retval;
! if (dim < 0)
! for (int i = 0; i < dims.length (); i++)
! if ( dims(i) > 1)
! {
! dim = i;
! break;
! }
! if (n_points < 0)
! n_points = dims (dim);
! else
! dims (dim) = n_points;
!
! if (dims.all_zero () || n_points == 0)
return octave_value (Matrix ());
if (arg.is_real_type ())
{
! NDArray nda = arg.array_value ();
if (! error_state)
{
! nda.resize (dims, 0.0);
! retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
}
}
else if (arg.is_complex_type ())
{
! ComplexNDArray cnda = arg.complex_matrix_value ();
if (! error_state)
{
! cnda.resize (dims, 0.0);
! retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
}
}
else
{
! gripe_wrong_type_arg (fcn, arg);
}
return retval;
}
+
+ DEFUN_DLD (fft, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} fft (@var{a}, @var{n}, @var{dim})\n\
+ Compute the FFT of @var{a} using subroutines from\n"
+ FFTSRC
+ ". The FFT is calculated along the first non-singleton dimension of the\n\
+ array. Thus if @var{a} is a matrix, @code{fft (@var{a})} computes the\n\
+ FFT for each column of @var{a}.\n\
+ \n\
+ If called with two arguments, @var{n} is expected to be an integer\n\
+ specifying the number of elements of @var{a} to use, or an empty\n\
+ matrix to specify that its value should be ignored. If @var{n} is\n\
+ larger than the dimension along which the FFT is calculated, then\n\
+ @var{a} is resized and padded with zeros. Otherwise, address@hidden is\n\
+ smaller than the dimension along which the FFT is calculated, then\n\
+ @var{a} is truncated.\n\
+ \n\
+ If called with three agruments, @var{dim} is an integer specifying the\n\
+ dimension of the matrix along which the FFT is performed\n\
+ @end deftypefn\n\
+ @seealso {ifft, fft2, fftn"
+ WISDOM
+ "}")
+ {
+ return do_fft(args, "fft", 0);
+ }
+
+
+ DEFUN_DLD (ifft, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} ifft (@var{a}, @var{n}, @var{dim})\n\
+ Compute the inverse FFT of @var{a} using subroutines from\n"
+ FFTSRC
+ ". The inverse FFT is calculated along the first non-singleton dimension\n\
+ of the array. Thus if @var{a} is a matrix, @code{fft (@var{a})} computes\n\
+ the inverse FFT for each column of @var{a}.\n\
+ \n\
+ If called with two arguments, @var{n} is expected to be an integer\n\
+ specifying the number of elements of @var{a} to use, or an empty\n\
+ matrix to specify that its value should be ignored. If @var{n} is\n\
+ larger than the dimension along which the inverse FFT is calculated, then\n\
+ @var{a} is resized and padded with zeros. Otherwise, address@hidden is\n\
+ smaller than the dimension along which the inverse FFT is calculated,\n\
+ then @var{a} is truncated.\n\
+ \n\
+ If called with three agruments, @var{dim} is an integer specifying the\n\
+ dimension of the matrix along which the inverse FFT is performed\n\
+ @end deftypefn\n\
+ @seealso {fft, ifft2, ifftn"
+ WISDOM
+ "}")
+ {
+ return do_fft(args, "ifft", 1);
+ }
+
/*
;;; Local Variables: ***
;;; mode: C++ ***
*** ./src/DLD-FUNCTIONS/fft2.cc.orig 2004-01-30 22:49:01.000000000 +0100
--- ./src/DLD-FUNCTIONS/fft2.cc 2004-02-11 14:12:33.000000000 +0100
***************
*** 1,5 ****
--- 1,6 ----
/*
+ Copyright (C) 2004 David Bateman
Copyright (C) 1996, 1997 John W. Eaton
This file is part of Octave.
***************
*** 32,49 ****
#include "oct-obj.h"
#include "utils.h"
! // This function should be merged with Fifft2.
! DEFUN_DLD (fft2, args, ,
! "-*- texinfo -*-\n\
! @deftypefn {Loadable Function} {} fft2 (@var{a}, @var{n}, @var{m})\n\
! Compute the two dimensional FFT of @var{a}.\n\
! \n\
! The optional arguments @var{n} and @var{m} may be used specify the\n\
! number of rows and columns of @var{a} to use. If either of these is\n\
! larger than the size of @var{a}, @var{a} is resized and padded with\n\
! zeros.\n\
! @end deftypefn")
{
octave_value retval;
--- 33,50 ----
#include "oct-obj.h"
#include "utils.h"
! // This function should be merged with Fifft.
! #if defined (HAVE_FFTW3)
! #define FFTSRC "@sc{Fftw}"
! #define WISDOM ", fft_wisdom"
! #else
! #define FFTSRC "@sc{Fftpack}"
! #define WISDOM
! #endif
!
! static octave_value
! do_fft2 (const octave_value_list &args, const char *fcn, int type)
{
octave_value retval;
***************
*** 51,129 ****
if (nargin < 1 || nargin > 3)
{
! print_usage ("fft2");
return retval;
}
octave_value arg = args(0);
!
! int n_rows = arg.rows ();
if (nargin > 1)
{
double dval = args(1).double_value ();
if (xisnan (dval))
! error ("fft2: NaN is invalid as N_ROWS");
else
! n_rows = NINT (dval);
}
if (error_state)
return retval;
! int n_cols = arg.columns ();
if (nargin > 2)
{
double dval = args(2).double_value ();
if (xisnan (dval))
! error ("fft2: NaN is invalid as N_COLS");
else
! n_cols = NINT (dval);
}
if (error_state)
return retval;
! if (n_rows < 0 || n_cols < 0)
! {
! error ("fft2: number of points must be greater than zero");
return retval;
- }
! int arg_is_empty = empty_arg ("fft2", arg.rows (), arg.columns ());
! if (arg_is_empty < 0)
! return retval;
! else if (arg_is_empty || n_rows == 0 || n_cols == 0)
return octave_value (Matrix ());
if (arg.is_real_type ())
{
! Matrix m = arg.matrix_value ();
if (! error_state)
{
! m.resize (n_rows, n_cols, 0.0);
! retval = m.fourier2d ();
}
}
else if (arg.is_complex_type ())
{
! ComplexMatrix m = arg.complex_matrix_value ();
if (! error_state)
{
! m.resize (n_rows, n_cols, 0.0);
! retval = m.fourier2d ();
}
}
else
{
! gripe_wrong_type_arg ("fft2", arg);
}
return retval;
}
/*
;;; Local Variables: ***
;;; mode: C++ ***
--- 52,184 ----
if (nargin < 1 || nargin > 3)
{
! print_usage (fcn);
return retval;
}
octave_value arg = args(0);
! dim_vector dims = arg.dims ();
! int n_rows = -1;
!
if (nargin > 1)
{
double dval = args(1).double_value ();
if (xisnan (dval))
! error ("%s: NaN is invalid as the N_ROWS", fcn);
else
! {
! n_rows = NINT (dval);
! if (n_rows < 0)
! error ("%s: number of rows must be greater than zero", fcn);
! }
}
if (error_state)
return retval;
! int n_cols = -1;
if (nargin > 2)
{
double dval = args(2).double_value ();
if (xisnan (dval))
! error ("%s: NaN is invalid as the N_COLS", fcn);
else
! {
! n_cols = NINT (dval);
! if (n_cols < 0)
! error ("%s: number of columns must be greater than zero", fcn);
! }
}
if (error_state)
return retval;
! for (int i = 0; i < dims.length (); i++)
! if (dims(i) < 0)
return retval;
! if (n_rows < 0)
! n_rows = dims (0);
! else
! dims (0) = n_rows;
! if (n_cols < 0)
! n_cols = dims (1);
! else
! dims (1) = n_cols;
!
! if (dims.all_zero () || n_rows == 0 || n_cols == 0)
return octave_value (Matrix ());
if (arg.is_real_type ())
{
! NDArray nda = arg.array_value ();
if (! error_state)
{
! nda.resize (dims, 0.0);
! retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
}
}
else if (arg.is_complex_type ())
{
! ComplexNDArray cnda = arg.complex_matrix_value ();
if (! error_state)
{
! cnda.resize (dims, 0.0);
! retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
}
}
else
{
! gripe_wrong_type_arg (fcn, arg);
}
return retval;
}
+ DEFUN_DLD (fft2, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} fft2 (@var{a}, @var{n}, @var{m})\n\
+ Compute the two dimensional FFT of @var{a} using subroutines from\n"
+ FFTSRC
+ ". The optional arguments @var{n} and @var{m} may be used specify the\n\
+ number of rows and columns of @var{a} to use. If either of these is\n\
+ larger than the size of @var{a}, @var{a} is resized and padded with\n\
+ zeros.\n\
+ \n\
+ If @var{a} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
+ of @var{a} is treated seperately\n\
+ @end deftypefn\n\
+ @seealso {ifft2, fft, fftn"
+ WISDOM
+ "}")
+ {
+ return do_fft2(args, "fft2", 0);
+ }
+
+
+ DEFUN_DLD (ifft2, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} fft2 (@var{a}, @var{n}, @var{m})\n\
+ Compute the inverse two dimensional FFT of @var{a} using subroutines from\n"
+ FFTSRC
+ ". The optional arguments @var{n} and @var{m} may be used specify the\n\
+ number of rows and columns of @var{a} to use. If either of these is\n\
+ larger than the size of @var{a}, @var{a} is resized and padded with\n\
+ zeros.\n\
+ \n\
+ If @var{a} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
+ of @var{a} is treated seperately\n\
+ @end deftypefn\n\
+ @seealso {fft2, ifft, ifftn"
+ WISDOM
+ "}")
+ {
+ return do_fft2(args, "ifft2", 1);
+ }
+
/*
;;; Local Variables: ***
;;; mode: C++ ***
*** ./src/DLD-FUNCTIONS/ifft2.cc.orig 2004-01-30 22:50:46.000000000 +0100
--- ./src/DLD-FUNCTIONS/ifft2.cc 2004-02-11 14:08:32.000000000 +0100
***************
*** 1,131 ****
- /*
-
- Copyright (C) 1996, 1997 John W. Eaton
-
- This file is part of Octave.
-
- Octave is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation; either version 2, or (at your option) any
- later version.
-
- Octave is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with Octave; see the file COPYING. If not, write to the Free
- Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-
- */
-
- #ifdef HAVE_CONFIG_H
- #include
- #endif
-
- #include "lo-mappers.h"
-
- #include "defun-dld.h"
- #include "error.h"
- #include "gripes.h"
- #include "oct-obj.h"
- #include "utils.h"
-
- // This function should be merged with Ffft2.
-
- DEFUN_DLD (ifft2, args, ,
- "-*- texinfo -*-\n\
- @deftypefn {Loadable Function} {} ifft2 (@var{a}, @var{n}, @var{m})\n\
- Compute the two dimensional inverse FFT of @var{a}.\n\
- \n\
- The optional arguments @var{n} and @var{m} may be used specify the\n\
- number of rows and columns of @var{a} to use. If either of these is\n\
- larger than the size of @var{a}, @var{a} is resized and padded with\n\
- zeros.\n\
- @end deftypefn")
- {
- octave_value retval;
-
- int nargin = args.length ();
-
- if (nargin < 1 || nargin > 3)
- {
- print_usage ("ifft2");
- return retval;
- }
-
- octave_value arg = args(0);
-
- int n_rows = arg.rows ();
- if (nargin > 1)
- {
- double dval = args(1).double_value ();
- if (xisnan (dval))
- error ("fft2: NaN is invalid as N_ROWS");
- else
- n_rows = NINT (dval);
- }
-
- if (error_state)
- return retval;
-
- int n_cols = arg.columns ();
- if (nargin > 2)
- {
- double dval = args(2).double_value ();
- if (xisnan (dval))
- error ("fft2: NaN is invalid as N_COLS");
- else
- n_cols = NINT (dval);
- }
-
- if (error_state)
- return retval;
-
- if (n_rows < 0 || n_cols < 0)
- {
- error ("ifft2: number of points must be greater than zero");
- return retval;
- }
-
- int arg_is_empty = empty_arg ("ifft2", arg.rows (), arg.columns ());
-
- if (arg_is_empty < 0)
- return retval;
- else if (arg_is_empty || n_rows == 0 || n_cols == 0)
- return octave_value (Matrix ());
-
- if (arg.is_real_type ())
- {
- Matrix m = arg.matrix_value ();
-
- if (! error_state)
- {
- m.resize (n_rows, n_cols, 0.0);
- retval = m.ifourier2d ();
- }
- }
- else if (arg.is_complex_type ())
- {
- ComplexMatrix m = arg.complex_matrix_value ();
-
- if (! error_state)
- {
- m.resize (n_rows, n_cols, 0.0);
- retval = m.ifourier2d ();
- }
- }
- else
- {
- gripe_wrong_type_arg ("ifft2", arg);
- }
-
- return retval;
- }
-
- /*
- ;;; Local Variables: ***
- ;;; mode: C++ ***
- ;;; End: ***
- */
--- 0 ----
*** ./src/DLD-FUNCTIONS/fft_wisdom.cc.orig 2004-02-11 14:09:27.000000000 +0100
--- ./src/DLD-FUNCTIONS/fft_wisdom.cc 2004-02-12 14:17:18.000000000 +0100
***************
*** 0 ****
--- 1,155 ----
+ /*
+
+ Copyright (C) 2004 David Bateman
+
+ This file is part of Octave.
+
+ Octave is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2, or (at your option) any
+ later version.
+
+ Octave is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Octave; see the file COPYING. If not, write to the Free
+ Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+ */
+
+ #ifdef HAVE_CONFIG_H
+ #include
+ #endif
+
+ #include
+
+ #include "oct-env.h"
+ #include "defaults.h"
+ #include "lo-mappers.h"
+ #include "defun-dld.h"
+ #include "error.h"
+ #include "gripes.h"
+ #include "oct-obj.h"
+ #include "utils.h"
+
+ DEFUN_DLD (fft_wisdom, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} fft_wisdom (@var{file}, @var{ow})\n\
+ @deftypefnx {Loadable Function} {} fft_wisdom (@var{n})\n\
+ \n\
+ This function is used to manipulate the FFTW wisdom data. It can\n\
+ be used to load previously stored wisdom from a file, create wisdom\n\
+ needed by Octave and to save the current wisdom to a file. Wisdom\n\
+ data can be used to significantly accelerate the calculation of the FFTs,\n\
+ but is only profitable if the same FFT is called many times due to the\n\
+ overhead in calculating the wisdom data.\n\
+ \n\
+ Called with a single string argument, @code{fft_wisdom (@var{file})}\n\
+ will load the wisdom data found in @var{file}. If @var{file} does\n\
+ not exist, then the current wisdom used by FFTW is saved to this\n\
+ file. If the flag @var{ow} is non-zero, then even if @var{file}\n\
+ exists, it will be overwritten.\n\
+ \n\
+ It is assumed that each row of @var{n} represents the size of a FFT for\n\
+ which it is desired to pre-calculate the wisdom needed to accelerate it.\n\
+ That is\n\
+ \n\
+ @example\n\
+ fftw_wisdom ([102; 103]);\n\
+ a = fft (rand (1,102));\n\
+ b = fft (rand (1,103));\n\
+ fftw_wisdom ([102, 103, 105]);\n\
+ c = fftn (rand ([102, 103, 105]));\n\
+ @end example\n\
+ \n\
+ calculates the wisdom necessary to accelerate the 102- and 103-point\n\
+ FFT, and the 102x103x105 3d FFT. Note that calculated wisdom will be lost\n\
+ when restarting Octave. However, if it is saved as discussed above, it\n\
+ can be reloaded. Also, any system-wide wisdom file that has been found\n\
+ will also be used. Saved wisdom files should not be used on different\n\
+ platforms since they will not be efficient and the point of calculating\n\
+ the wisdom is lost.\n\
+ \n\
+ Note that the program @sc{fftw-wisdom} supplied with FFTW can equally be\n\
+ used to create a file containing wisdom that can be imported into Octave.\n\
+ @end deftypefn\n\
+ @seealso {fft, ifft, fft2, ifft2, fftn, ifftn}")
+ {
+ octave_value retval;
+ int nargin = args.length();
+
+ if (nargin < 1 || nargin > 2)
+ {
+ print_usage ("fft_wisdom");
+ return retval;
+ }
+
+ if (args(0).is_string ())
+ {
+ bool overwrite = false;
+
+ if (nargin != 1)
+ {
+ double dval = args (1).double_value ();
+ if (NINT (dval) != 0)
+ overwrite = true;
+ }
+
+ std::string wisdom = octave_env::make_absolute
+ (Vload_path_dir_path.find_first_of (args(0).string_value ().c_str ()),
+ octave_env::getcwd ());
+
+ if (wisdom.empty () || overwrite)
+ {
+ FILE *ofile = fopen (wisdom.c_str (), "wb");
+ fftw_export_wisdom_to_file (ofile);
+ fclose (ofile);
+ }
+ else
+ {
+ FILE *ifile = fopen (wisdom.c_str (), "r");
+ if (! fftw_import_wisdom_from_file (ifile))
+ error ("fft_wisdom: can not import wisdom from file");
+ fclose (ifile);
+ }
+
+ }
+ else
+ {
+ #if 1
+ error ("fft_wisdom: XXX FIXME XXX not yet implemented!!");
+ #else
+ Matrix m = args (0).matrix_value ();
+
+ if (error_state)
+ {
+ error ("fft_wisdom: argument must be a matrix or a string");
+ return retval;
+ }
+
+ rank = m.columns ();
+ fftw_plan plan;
+ OCTAVE_LOCAL_BUFFER (int, n, rank);
+
+ for (int k = 0; k < m.rows (); k++)
+ {
+ int fftsize = 1;
+
+ // Note reversal of dimensions for column major storage in FFTW
+ for (int i = 0, j = rank-1; i < rank; i++, j--)
+ {
+ n[i] = m(k,j);
+ fftsize *= n[i];
+ }
+
+ // XXX FIXME add the planning here to get the wisdom!!!
+
+ }
+ #endif
+ }
+
+ return retval;
+ }
*** ./src/Makefile.in.orig 2004-02-11 14:08:02.000000000 +0100
--- ./src/Makefile.in 2004-02-11 14:08:05.000000000 +0100
***************
*** 45,52 ****
DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc \
daspk.cc dasrt.cc dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc \
! filter.cc find.cc fsolve.cc gammainc.cc getgrent.cc \
! getpwent.cc getrusage.cc givens.cc hess.cc ifft.cc ifft2.cc \
inv.cc kron.cc lpsolve.cc lsode.cc lu.cc minmax.cc \
odessa.cc pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc \
sort.cc sqrtm.cc svd.cc syl.cc time.cc
--- 45,52 ----
DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc \
daspk.cc dasrt.cc dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc \
! fftn.cc fft_wisdom.cc filter.cc find.cc fsolve.cc gammainc.cc \
! getgrent.cc getpwent.cc getrusage.cc givens.cc hess.cc \
inv.cc kron.cc lpsolve.cc lsode.cc lu.cc minmax.cc \
odessa.cc pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc \
sort.cc sqrtm.cc svd.cc syl.cc time.cc
*** ./src/ChangeLog.orig 2004-01-30 23:14:31.000000000 +0100
--- ./src/ChangeLog 2004-02-11 14:02:10.000000000 +0100
***************
*** 1,3 ****
--- 1,15 ----
+ 2004-02-11 Dvaid Bateman
+
+ * DLD-FUNCTIONS/fft.cc: Adapt for Nd arrays, combine with ifft.cc
+ * DLD-FUNCTIONS/ifft.cc: remove and link to fft.cc
+ * DLD-FUNCTIONS/fft2.cc: Adapt for Nd arrays, combine with ifft.cc
+ * DLD-FUNCTIONS/ifft2.cc: remove and link to fft.cc
+ * DLD-FUNCTIONS/fftn.cc: new function for Nd FFT, combined with
+ ifftn.cc
+ * DLD-FUNCTIONS/ifftn.cc: link to fftn.cc
+ * DLD-FUNCTIONS/fft_wisdom.cc: New function to manipulate FFTW
+ wisdom
+
2004-01-22 John W. Eaton
* error.cc (pr_where): New arg, print_code with default value true.
*** ./configure.in.orig 2004-01-28 15:16:49.000000000 +0100
--- ./configure.in 2004-02-10 10:50:54.000000000 +0100
***************
*** 407,425 ****
with_fftw=$withval, with_fftw=yes)
if test "$with_fftw" = "yes"; then
! have_fftw_header=no
! AC_CHECK_HEADERS(dfftw.h fftw.h, [have_fftw_header=yes; break])
! if test "$have_fftw_header" = yes; then
! AC_CHECK_LIB(dfftw, fftw_create_plan, FFTW_LIBS="-ldfftw",
! [AC_CHECK_LIB(fftw, fftw_create_plan, FFTW_LIBS="-lfftw", with_fftw=no)])
! else
! with_fftw=no
fi
fi
! if test "$with_fftw" = yes; then
FFT_DIR=''
! AC_DEFINE(HAVE_FFTW, 1, [Define if the FFTW library is available.])
fi
WITH_MPI=true
--- 407,423 ----
with_fftw=$withval, with_fftw=yes)
if test "$with_fftw" = "yes"; then
! have_fftw3_header=no
! with_fftw3=no
! AC_CHECK_HEADER(fftw3.h, [have_fftw3_header=yes; break])
! if test "$have_fftw3_header" = yes; then
! AC_CHECK_LIB(fftw3, fftw_plan_dft_1d, [FFTW_LIBS="-lfftw3"; with_fftw3=yes])
fi
fi
! if test "$with_fftw" = yes && test "$with_fftw3" = yes; then
FFT_DIR=''
! AC_DEFINE(HAVE_FFTW3, 1, [Define if the FFTW3 library is used.])
fi
WITH_MPI=true
*** ./ChangeLog.orig 2004-01-30 15:25:15.000000000 +0100
--- ./ChangeLog 2004-01-30 15:26:32.000000000 +0100
***************
*** 1,3 ****
--- 1,8 ----
+ 2004-01-30 David Bateman
+
+ * configure.in: Test for the presence of FFTW 3.x and use it in
+ preference to FFTW 2.x. Define HAVE_FFTW3
+
2004-01-22 John W. Eaton
* octMakefile.in (maintainer-clean, distclean):