# HG changeset patch # User jwe # Date 1076958153 0 # Node ID ccfbd6047a540b272b3ba1443d5bc76cc350f35e # Parent 9eed17b2c8d13d644dbe25891628b49f2a364ec6 [project @ 2004-02-16 19:02:32 by jwe] diff -r 9eed17b2c8d1 -r ccfbd6047a54 ChangeLog --- a/ChangeLog Mon Feb 16 17:57:34 2004 +0000 +++ b/ChangeLog Mon Feb 16 19:02:33 2004 +0000 @@ -1,3 +1,8 @@ +2004-02-16 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-02-16 John W. Eaton * mkoctfile.in (LINK_DEPS): Include $LIBS and $RLD_FLAG. diff -r 9eed17b2c8d1 -r ccfbd6047a54 configure.in --- a/configure.in Mon Feb 16 17:57:34 2004 +0000 +++ b/configure.in Mon Feb 16 19:02:33 2004 +0000 @@ -29,7 +29,7 @@ EXTERN_CXXFLAGS="$CXXFLAGS" AC_INIT -AC_REVISION($Revision: 1.445 $) +AC_REVISION($Revision: 1.446 $) AC_PREREQ(2.57) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -407,19 +407,17 @@ 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 + 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; then +if test "$with_fftw" = yes && test "$with_fftw3" = yes; then FFT_DIR='' - AC_DEFINE(HAVE_FFTW, 1, [Define if the FFTW library is available.]) + AC_DEFINE(HAVE_FFTW3, 1, [Define if the FFTW3 library is used.]) fi WITH_MPI=true diff -r 9eed17b2c8d1 -r ccfbd6047a54 doc/interpreter/signal.txi --- a/doc/interpreter/signal.txi Mon Feb 16 17:57:34 2004 +0000 +++ b/doc/interpreter/signal.txi Mon Feb 16 19:02:33 2004 +0000 @@ -19,6 +19,10 @@ @DOCSTRING(ifft2) +@DOCSTRING(fftn) + +@DOCSTRING(ifftn) + @DOCSTRING(fftconv) @DOCSTRING(fftfilt) diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/CMatrix.cc --- a/liboctave/CMatrix.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/CMatrix.cc Mon Feb 16 19:02:33 2004 +0000 @@ -56,7 +56,7 @@ #include "mx-inlines.cc" #include "oct-cmplx.h" -#ifdef HAVE_FFTW +#if defined (HAVE_FFTW3) #include "oct-fftw.h" #endif @@ -1102,7 +1102,7 @@ return retval; } -#ifdef HAVE_FFTW +#if defined (HAVE_FFTW3) ComplexMatrix ComplexMatrix::fourier (void) const @@ -1128,12 +1128,7 @@ 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); - } + octave_fftw::fft (in, out, npts, nsamples); return retval; } @@ -1162,12 +1157,7 @@ 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); - } + octave_fftw::ifft (in, out, npts, nsamples); return retval; } @@ -1175,13 +1165,13 @@ 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); + 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,13 +1179,13 @@ 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); + 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; } @@ -1332,8 +1322,8 @@ wsave.resize (nn); pwsave = wsave.fortran_vec (); - Array row (npts); - Complex *prow = row.fortran_vec (); + Array tmp (npts); + Complex *prow = tmp.fortran_vec (); F77_FUNC (cffti, CFFTI) (npts, pwsave); @@ -1401,8 +1391,8 @@ wsave.resize (nn); pwsave = wsave.fortran_vec (); - Array row (npts); - Complex *prow = row.fortran_vec (); + Array tmp (npts); + Complex *prow = tmp.fortran_vec (); F77_FUNC (cffti, CFFTI) (npts, pwsave); diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/CNDArray.cc --- a/liboctave/CNDArray.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/CNDArray.cc Mon Feb 16 19:02:33 2004 +0000 @@ -34,9 +34,31 @@ #include "Array-util.h" #include "CNDArray.h" #include "mx-base.h" +#include "f77-fcn.h" #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,6 +83,429 @@ 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] / + static_cast (npts); + } + } + + return retval; +} + +ComplexNDArray +ComplexNDArray::fourier2d (void) const +{ + dim_vector dv = dims (); + dim_vector dv2 (dv(0), dv(1)); + int rank = 2; + ComplexNDArray retval (*this); + int stride = 1; + + for (int i = 0; i < rank; i++) + { + int npts = dv2(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 *= dv2(i); + } + + return retval; +} + +ComplexNDArray +ComplexNDArray::ifourier2d (void) const +{ + dim_vector dv = dims(); + dim_vector dv2 (dv(0), dv(1)); + int rank = 2; + ComplexNDArray retval (*this); + int stride = 1; + + for (int i = 0; i < rank; i++) + { + int npts = dv2(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] / + static_cast (npts); + } + } + + stride *= dv2(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] / + static_cast (npts); + } + } + + stride *= dv(i); + } + + return retval; +} + +#endif + // unary operations boolNDArray diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/CNDArray.h --- a/liboctave/CNDArray.h Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/CNDArray.h Mon Feb 16 19:02:33 2004 +0000 @@ -94,6 +94,15 @@ NDArray abs (void) 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; + ComplexMatrix matrix_value (void) const; ComplexNDArray squeeze (void) const { return ArrayN::squeeze (); } diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/ChangeLog --- a/liboctave/ChangeLog Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/ChangeLog Mon Feb 16 19:02:33 2004 +0000 @@ -1,3 +1,31 @@ +2004-02-16 David Bateman + + * oct-fftw.cc (octave_fftw_planner::create_plan, octave_fftw::fftNd): + Add support for FFTW 3.x. Include the ability to + use the real to complex transform for fft's of real matrices + (octave_fftw_planner::create_plan2d): Delete. + (octave_fftw::fft2d): Delete. + (convert_packcomplex_1d, convert_packcomplex_Nd): + New static functions. + * oct-fftw.h: Update decls. + + * dMatrix.cc (Matrix::fourier, Matrix::ifourier, + Matrix::fourier2d, Matrix::ifourier2d): 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 (ComplexMatrix::fourier, ComplexMatrix::ifourier, + ComplexMatrix::fourier2d, ComplexMatrix::ifourier2d): 1D fft of a + matrix done as single call rather than loop. Update for FFTW 3.x. + + * dNDArray.cc (NDArray::fourier, NDArray::ifourier, + NDArray::fourierNd, NDArray::ifouriourNd): New fourier transform + functions for Nd arrays. + * dNArray.h Provide decls. + * CNDArray.cc (ComplexNDArray::fourier, ComplexNDArray::ifourier, + ComplexNDArray::fourierNd, ComplexNDArray::ifouriourNd): New + fourier transform functions for complex Nd arrays. + * CNArray.h: Provide decls. + 2004-02-15 Petter Risholm * Array.cc (Array::insert (const Array&, int, int)): diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/dMatrix.cc --- a/liboctave/dMatrix.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/dMatrix.cc Mon Feb 16 19:02:33 2004 +0000 @@ -52,7 +52,7 @@ #include "oct-cmplx.h" #include "quit.h" -#ifdef HAVE_FFTW +#if defined (HAVE_FFTW3) #include "oct-fftw.h" #endif @@ -766,7 +766,7 @@ } } -#ifdef HAVE_FFTW +#if defined (HAVE_FFTW3) ComplexMatrix Matrix::fourier (void) const @@ -789,16 +789,10 @@ nsamples = nc; } - ComplexMatrix tmp (*this); - Complex *in (tmp.fortran_vec ()); + const double *in (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); - } + octave_fftw::fft (in, out, npts, nsamples); return retval; } @@ -828,12 +822,7 @@ 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); - } + octave_fftw::ifft (in, out, npts, nsamples); return retval; } @@ -841,13 +830,11 @@ 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); + 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,13 +842,12 @@ ComplexMatrix Matrix::ifourier2d (void) const { - int nr = rows (); - int nc = cols (); + dim_vector dv(rows (), 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); + Complex *out (retval.fortran_vec ()); + + octave_fftw::ifftNd (out, out, 2, dv); return retval; } @@ -998,8 +984,8 @@ wsave.resize (nn); pwsave = wsave.fortran_vec (); - Array row (npts); - Complex *prow = row.fortran_vec (); + Array tmp (npts); + Complex *prow = tmp.fortran_vec (); F77_FUNC (cffti, CFFTI) (npts, pwsave); @@ -1067,8 +1053,8 @@ wsave.resize (nn); pwsave = wsave.fortran_vec (); - Array row (npts); - Complex *prow = row.fortran_vec (); + Array tmp (npts); + Complex *prow = tmp.fortran_vec (); F77_FUNC (cffti, CFFTI) (npts, pwsave); diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/dNDArray.cc --- a/liboctave/dNDArray.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/dNDArray.cc Mon Feb 16 19:02:33 2004 +0000 @@ -34,10 +34,453 @@ #include "Array-util.h" #include "dNDArray.h" #include "mx-base.h" +#include "f77-fcn.h" #include "lo-error.h" #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] / + static_cast (npts); + } + } + + return retval; +} + +ComplexNDArray +NDArray::fourier2d (void) const +{ + dim_vector dv = dims(); + dim_vector dv2 (dv(0), dv(1)); + int rank = 2; + ComplexNDArray retval (*this); + int stride = 1; + + for (int i = 0; i < rank; i++) + { + int npts = dv2(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 *= dv2(i); + } + + return retval; +} + +ComplexNDArray +NDArray::ifourier2d (void) const +{ + dim_vector dv = dims(); + dim_vector dv2 (dv(0), dv(1)); + int rank = 2; + ComplexNDArray retval (*this); + int stride = 1; + + for (int i = 0; i < rank; i++) + { + int npts = dv2(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] / + static_cast (npts); + } + } + + stride *= dv2(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] / + static_cast (npts); + } + } + + stride *= dv(i); + } + + return retval; +} + +#endif + NDArray::NDArray (const boolNDArray& a) : MArrayN (a.dims ()) { diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/dNDArray.h --- a/liboctave/dNDArray.h Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/dNDArray.h Mon Feb 16 19:02:33 2004 +0000 @@ -90,6 +90,15 @@ NDArray abs (void) 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; + friend NDArray real (const ComplexNDArray& a); friend NDArray imag (const ComplexNDArray& a); diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/oct-fftw.cc --- a/liboctave/oct-fftw.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/oct-fftw.cc Mon Feb 16 19:02:33 2004 +0000 @@ -22,18 +22,26 @@ #include #endif -#ifdef HAVE_FFTW +#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. Future -// enhancement will be to add a dynamically loadable interface ("fftw") -// to manipulate fftw wisdom so that users may choose the appropriate -// planner. +// 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,17 +49,35 @@ public: octave_fftw_planner (); - fftw_plan create_plan (fftw_direction, size_t); - fftwnd_plan create_plan2d (fftw_direction, size_t, size_t); + 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]; - fftwnd_plan plan2d[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]; - size_t n[2]; - size_t n2d[2][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,31 +85,67 @@ 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; + 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 (fftw_direction dir, size_t npts) +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) { - size_t which = (dir == FFTW_FORWARD) ? 0 : 1; + 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 || n[which] != npts) - { - create_new_plan = true; - n[which] = npts; - } + 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); - *cur_plan_p = fftw_create_plan (npts, dir, plan_flags); + // 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,32 +154,55 @@ return *cur_plan_p; } -fftwnd_plan -octave_fftw_planner::create_plan2d (fftw_direction dir, - size_t nrows, size_t ncols) +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) { - size_t which = (dir == FFTW_FORWARD) ? 0 : 1; - fftwnd_plan *cur_plan_p = &plan2d[which]; + fftw_plan *cur_plan_p = &rplan; bool create_new_plan = false; + char in_align = ((int) in) & 0xF; + char out_align = ((int) out) & 0xF; - 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 (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) - fftwnd_destroy_plan (*cur_plan_p); + fftw_destroy_plan (*cur_plan_p); - *cur_plan_p = fftw2d_create_plan (nrows, ncols, dir, - plan_flags | FFTW_IN_PLACE); + // 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 2d fftw plan"); + (*current_liboctave_error_handler) ("Error creating fftw plan"); } return *cur_plan_p; @@ -125,51 +210,184 @@ static octave_fftw_planner fftw_planner; -int -octave_fftw::fft (const Complex *in, Complex *out, size_t npts) +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) { - fftw_one (fftw_planner.create_plan (FFTW_FORWARD, npts), - reinterpret_cast (const_cast (in)), - reinterpret_cast (out)); + 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 the rank = 2 case directly for speed + 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 + i*nr*nc] = conj(out[(i*nr+1)*nc - j]); + } + + // Now do the permutations needed for rank > 2 cases + size_t jstart = dv(0) * dv(1); + size_t kstep = dv(0); + size_t nel = dv.numel (); + for (int inner = 2; inner < dv.length(); inner++) + { + size_t jmax = jstart * dv(inner); + for (size_t i = 0; i < nel; i+=jmax) + for (size_t j = jstart, jj = jmax-jstart; j < jj; + j+=jstart, jj-=jstart) + for (size_t k = 0; k < jstart; k+= kstep) + for (size_t l = nc/2+1; l < nc; l++) + { + Complex tmp = out[i+ j + k + l]; + out[i + j + k + l] = out[i + jj + k + l]; + out[i + jj + k + l] = tmp; + } + jstart = jmax; + } +} + +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::ifft (const Complex *in, Complex *out, size_t npts) +octave_fftw::fft (const Complex *in, Complex *out, size_t npts, + size_t nsamples, int stride, int dist) { - fftw_one (fftw_planner.create_plan (FFTW_BACKWARD, npts), - reinterpret_cast (const_cast (in)), - reinterpret_cast (out)); + 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 i = 0; i < npts; i++) - out[i] /= scale; + 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::fft2d (Complex *inout, size_t nr, size_t nc) +octave_fftw::fftNd (const double *in, Complex *out, const int rank, + const dim_vector &dv) { - fftwnd_one (fftw_planner.create_plan2d (FFTW_FORWARD, nr, nc), - reinterpret_cast (inout), - 0); + 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::ifft2d (Complex *inout, size_t nr, size_t nc) +octave_fftw::fftNd (const Complex *in, Complex *out, const int rank, + const dim_vector &dv) { - fftwnd_one (fftw_planner.create_plan2d (FFTW_BACKWARD, nr, nc), - reinterpret_cast (inout), - 0); + 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; +} - const size_t npts = nr * nc; +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++) - inout[i] /= scale; + out[i] /= scale; return 0; } diff -r 9eed17b2c8d1 -r ccfbd6047a54 liboctave/oct-fftw.h --- a/liboctave/oct-fftw.h Mon Feb 16 17:57:34 2004 +0000 +++ b/liboctave/oct-fftw.h Mon Feb 16 19:02:33 2004 +0000 @@ -22,24 +22,27 @@ #define octave_oct_fftw_h 1 #include - -#if defined (HAVE_DFFTW_H) -#include -#else -#include -#endif +#include #include "oct-cmplx.h" +#include "dim-vector.h" class octave_fftw { public: - static int fft (const Complex*, Complex *, size_t); - static int ifft (const Complex*, Complex *, size_t); + 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 fft2d (Complex*, size_t, size_t); - static int ifft2d (Complex*, size_t, size_t); + 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 (); diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/ChangeLog --- a/src/ChangeLog Mon Feb 16 17:57:34 2004 +0000 +++ b/src/ChangeLog Mon Feb 16 19:02:33 2004 +0000 @@ -1,3 +1,20 @@ +2004-02-16 David Bateman + + * DLD-FUNCTIONS/fft.cc: Adapt for Nd arrays, combine with ifft.cc. + * DLD-FUNCTIONS/ifft.cc: Delete. + * DLD-FUNCTIONS/fft2.cc: Adapt for Nd arrays, combine with ifft.cc. + * DLD-FUNCTIONS/ifft2.cc: Delete. + * DLD-FUNCTIONS/fftn.cc: New function for Nd FFT and inverse FFT. + * DLD-FUNCTIONS/fft_wisdom.cc: New function to manipulate FFTW + wisdom. + + * Makefile.in: Remove ifft.cc and ifft2.cc. Add fftn.cc and + ifftn.cc. Make building of fft-wisdom.cc conditional on the + value of FFTW_LIBS. + * defaults.cc (Vwisdom_prog): New variable + (set_default_wisdom_prog, wisdom_prog): New functions. + (symbols_of_defaults): Add DEFVAR for wisdom_prog. + 2004-02-16 John W. Eaton * ov-list.cc (octave_list::subsasgn): Call assign for Cell objects diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/DLD-FUNCTIONS/fft.cc --- a/src/DLD-FUNCTIONS/fft.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/src/DLD-FUNCTIONS/fft.cc Mon Feb 16 19:02:33 2004 +0000 @@ -1,5 +1,6 @@ /* +Copyright (C) 2004 David Bateman Copyright (C) 1996, 1997 John W. Eaton This file is part of Octave. @@ -32,96 +33,169 @@ #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 -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") +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 > 2) + if (nargin < 1 || nargin > 3) { - print_usage ("fft"); + print_usage (fcn); return retval; } octave_value arg = args(0); - - int n_points = arg.rows (); - if (n_points == 1) - n_points = arg.columns (); + 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 (nargin == 2) + if (error_state) + return retval; + + if (nargin > 2) { - double dval = args(1).double_value (); + double dval = args(2).double_value (); if (xisnan (dval)) - error ("fft: NaN is invalid as the N_POINTS"); + 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 - n_points = NINT (dval); + dim = NINT (dval) - 1; } if (error_state) return retval; - if (n_points < 0) - { - error ("fft: number of points must be greater than zero"); + for (int i = 0; i < dims.length (); i++) + if (dims(i) < 0) return retval; - } - int arg_is_empty = empty_arg ("fft", arg.rows (), arg.columns ()); + if (dim < 0) + for (int i = 0; i < dims.length (); i++) + if ( dims(i) > 1) + { + dim = i; + break; + } - if (arg_is_empty < 0) - return retval; - else if (arg_is_empty || n_points == 0) + 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 ()) { - Matrix m = arg.matrix_value (); + NDArray nda = arg.array_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 (); + nda.resize (dims, 0.0); + retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim)); } } else if (arg.is_complex_type ()) { - ComplexMatrix m = arg.complex_matrix_value (); + ComplexNDArray cnda = arg.complex_array_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 (); + cnda.resize (dims, 0.0); + retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim)); } } else { - gripe_wrong_type_arg ("fft", arg); + 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, if@var{n} 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, if@var{n} 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++ *** diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/DLD-FUNCTIONS/fft2.cc --- a/src/DLD-FUNCTIONS/fft2.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/src/DLD-FUNCTIONS/fft2.cc Mon Feb 16 19:02:33 2004 +0000 @@ -1,5 +1,6 @@ /* +Copyright (C) 2004 David Bateman Copyright (C) 1996, 1997 John W. Eaton This file is part of Octave. @@ -32,18 +33,18 @@ #include "oct-obj.h" #include "utils.h" -// This function should be merged with Fifft2. +// This function should be merged with Fifft. -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") +#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,79 +52,133 @@ if (nargin < 1 || nargin > 3) { - print_usage ("fft2"); + print_usage (fcn); return retval; } octave_value arg = args(0); - - int n_rows = arg.rows (); + dim_vector dims = arg.dims (); + int n_rows = -1; + if (nargin > 1) { double dval = args(1).double_value (); if (xisnan (dval)) - error ("fft2: NaN is invalid as N_ROWS"); + error ("%s: NaN is invalid as the N_ROWS", fcn); else - n_rows = NINT (dval); + { + 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 = arg.columns (); + int n_cols = -1; if (nargin > 2) { double dval = args(2).double_value (); if (xisnan (dval)) - error ("fft2: NaN is invalid as N_COLS"); + error ("%s: NaN is invalid as the N_COLS", fcn); else - n_cols = NINT (dval); + { + n_cols = NINT (dval); + if (n_cols < 0) + error ("%s: number of columns must be greater than zero", fcn); + } } if (error_state) return retval; - if (n_rows < 0 || n_cols < 0) - { - error ("fft2: number of points must be greater than zero"); + for (int i = 0; i < dims.length (); i++) + if (dims(i) < 0) return retval; - } - int arg_is_empty = empty_arg ("fft2", arg.rows (), arg.columns ()); + if (n_rows < 0) + n_rows = dims (0); + else + dims (0) = n_rows; - if (arg_is_empty < 0) - return retval; - else if (arg_is_empty || n_rows == 0 || n_cols == 0) + 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 ()) { - Matrix m = arg.matrix_value (); + NDArray nda = arg.array_value (); if (! error_state) { - m.resize (n_rows, n_cols, 0.0); - retval = m.fourier2d (); + nda.resize (dims, 0.0); + retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ()); } } else if (arg.is_complex_type ()) { - ComplexMatrix m = arg.complex_matrix_value (); + ComplexNDArray cnda = arg.complex_array_value (); if (! error_state) { - m.resize (n_rows, n_cols, 0.0); - retval = m.fourier2d (); + cnda.resize (dims, 0.0); + retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ()); } } else { - gripe_wrong_type_arg ("fft2", arg); + 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++ *** diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/DLD-FUNCTIONS/fft_wisdom.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/fft_wisdom.cc Mon Feb 16 19:02:33 2004 +0000 @@ -0,0 +1,190 @@ +/* + +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 "file-ops.h" +#include "lo-sstream.h" +#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" +#include "sighandlers.h" + +// Name of the FFTW wisdom program we'd like to use. +extern std::string Vwisdom_prog; + +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\ +Any value of the matrix that is less than 1, is assumed to indicate an\n\ +absent dimension. That is\n\ +\n\ +@example\n\ +fftw_wisdom ([102, 0, 0; 103, 103, 0; 102, 103, 105]);\n\ +a = fft (rand (1,102));\n\ +b = fft (rand (103,103));\n\ +c = fftn (rand ([102, 103, 105]));\n\ +@end example\n\ +\n\ +calculates the wisdom necessary to accelerate the 103, 102x102 and\n\ +the 102x103x105 FFTs. Note that calculated wisdom will be lost when\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 @code{fftw-wisdom} supplied with FFTW can equally\n\ +be used to create a file containing wisdom that can be imported into\n\ +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 + { + Matrix m = args (0).matrix_value (); + + if (error_state) + { + error ("fft_wisdom: argument must be a matrix or a string"); + return retval; + } + + std::string name = file_ops::tempnam ("", "oct-"); + + if (name.empty ()) + { + error ("fft_wisdom: can not open temporary file"); + return retval; + } + + OSSTREAM cmd_buf; + cmd_buf << Vwisdom_prog << " -n -o \"" << name << "\""; + + for (int k = 0; k < m.rows (); k++) + { + bool first = true; + + cmd_buf << " "; + + // Note reversal of dimensions for column major storage in FFTW + for (int j = m.columns()-1; j >= 0; j--) + if (NINT(m(k,j)) > 0) + { + if (first) + first = false; + else + cmd_buf << "x"; + cmd_buf << NINT(m(k,j)) ; + } + } + + cmd_buf << OSSTREAM_ENDS; + + volatile octave_interrupt_handler old_interrupt_handler + = octave_ignore_interrupts (); + + int status = system (OSSTREAM_C_STR (cmd_buf)); + + OSSTREAM_FREEZE (cmd_buf); + + octave_set_interrupt_handler (old_interrupt_handler); + + if (WIFEXITED (status)) + { + FILE *ifile = fopen (name.c_str (), "r"); + if (! fftw_import_wisdom_from_file (ifile)) + error ("fft_wisdom: can not import wisdom from temporary file"); + fclose (ifile); + } + else + error ("fft_wisdom: error running %s", Vwisdom_prog.c_str ()); + + } + + return retval; +} diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/DLD-FUNCTIONS/fftn.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/fftn.cc Mon Feb 16 19:02:33 2004 +0000 @@ -0,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 () || val.rows () != 1) + error ("%s: second argument must be a vector of length dim", 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_array_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: *** +*/ diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/DLD-FUNCTIONS/ifft.cc --- a/src/DLD-FUNCTIONS/ifft.cc Mon Feb 16 17:57:34 2004 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,130 +0,0 @@ -/* - -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: *** -*/ diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/DLD-FUNCTIONS/ifft2.cc --- a/src/DLD-FUNCTIONS/ifft2.cc Mon Feb 16 17:57:34 2004 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,131 +0,0 @@ -/* - -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: *** -*/ diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/Makefile.in --- a/src/Makefile.in Mon Feb 16 17:57:34 2004 +0000 +++ b/src/Makefile.in Mon Feb 16 19:02:33 2004 +0000 @@ -43,10 +43,14 @@ OPT_HANDLERS := DASPK-opts.cc DASRT-opts.cc DASSL-opts.cc \ LSODE-opts.cc NLEqn-opts.cc ODESSA-opts.cc Quad-opts.cc +ifneq ($(FFTW_LIBS), ) + DLD_WISDOM := fft_wisdom.cc +endif + 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 \ + fftn.cc $(DLD_WISDOM) 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 diff -r 9eed17b2c8d1 -r ccfbd6047a54 src/defaults.cc --- a/src/defaults.cc Mon Feb 16 17:57:34 2004 +0000 +++ b/src/defaults.cc Mon Feb 16 19:02:33 2004 +0000 @@ -92,6 +92,11 @@ std::string Vlocal_site_defaults_file; std::string Vsite_defaults_file; +#ifdef HAVE_FFTW3 +// Name of the FFTW wisdom program +std::string Vwisdom_prog; +#endif + // Each element of A and B should be directory names. For each // element of A not in the list B, execute SCRIPT_FILE in that // directory if it exists. @@ -310,6 +315,19 @@ Vinfo_prog = std::string (oct_info_prog); } +#ifdef HAVE_FFTW3 +static void +set_default_wisdom_prog (void) +{ + std::string oct_wisdom_prog = octave_env::getenv ("OCTAVE_WISDOM_PROGRAM"); + + if (oct_wisdom_prog.empty ()) + Vwisdom_prog = "fftw-wisdom"; + else + Vwisdom_prog = std::string (oct_wisdom_prog); +} +#endif + static void set_default_editor (void) { @@ -400,6 +418,10 @@ set_default_info_prog (); +#ifdef HAVE_FFTW3 + set_default_wisdom_prog (); +#endif + set_default_editor (); set_local_site_defaults_file (); @@ -484,6 +506,26 @@ return status; } +#ifdef HAVE_FFTW3 +static int +wisdom_prog (void) +{ + int status = 0; + + std::string s = builtin_string_variable ("WISDOM_PROGRAM"); + + if (s.empty ()) + { + gripe_invalid_value_specified ("WISDOM_PROGRAM"); + status = -1; + } + else + Vwisdom_prog = s; + + return status; +} +#endif + static int default_exec_path (void) { @@ -587,6 +629,18 @@ @code{\"emacs\"}.\n\ @end defvr"); +#ifdef HAVE_FFTW3 + DEFVAR (WISDOM_PROGRAM, Vwisdom_prog, wisdom_prog, + "-*- texinfo -*-\n\ +@defvr {Built-in Variable} WISDOM_PROGRAM\n\ +A string naming the FFTW wisdom program to use to create wisdom data\n\ +to accelerate Fourier transforms. If the environment variable\n\ +@code{OCTAVE_WISDOM_PROGRAM} is set when Octave starts, its value is used\n\ +as the default. Otherwise, @code{WISDOM_PROGRAM} is set to\n\ +@code{\"fftw-wisdom\"}.\n\ +@end defvr"); +#endif + DEFVAR (EXEC_PATH, Vexec_path, exec_path, "-*- texinfo -*-\n\ @defvr {Built-in Variable} EXEC_PATH\n\ @@ -621,7 +675,7 @@ substituted for leading, trailing, or doubled colons that appear in the\n\ built-in variable @code{EXEC_PATH}.\n\ @end defvr"); - + DEFVAR (LOADPATH, Vload_path, loadpath, "-*- texinfo -*-\n\ @defvr {Built-in Variable} LOADPATH\n\