From octave-sources-request at bevo dot che dot wisc dot edu Thu May 20 08:30:25 1999 Subject: Use 'fftw-2.1' for fast fourier transformation if you want From: Thomas Walter To: octave-sources at bevo dot che dot wisc dot edu Date: Thu, 20 May 1999 15:30:30 +0200 Hello, I implemented a test interface to use FFTW for fast fourier transformations from octave. The source is C++ in 'fftw.cc' and implements currently 1D transformations only. The source it at the end and contains the call of 'mkoctfile' too. You need fftw compiled as shared library. For more about FFTW see http://theory.lcs.mit.edu/~fftw/homepage.html Hints are welcome Bye Thomas =========== fftw.cc start ==================== // Hi emacs! My mode is -*- C++ -*- // // Use 'fftw-2.1' for fast fourier transformation of 1 dimensional data. // This library can do: // fft any size without zero-padding // fft of higher dimensions (multi-dimensions) // can use multi processor/threaded architectures. // // See home site of FFTW at // // http://theory.lcs.mit.edu/~fftw/homepage.html // // To get a '.oct' file use: // env SH_LDFLAGS='-shared -Wl,-rpath,/usr/local/lib' mkoctfile [-v] [-s] -lfftw fftw.cc // followes by: // ln [-s] fftw.oct ifftw.oct // // Note: '-s' strips symbols // // Copyright (c) 1999 Thomas Walter // $Id: fftw.cc,v 1.8 1999/05/20 13:21:47 twalter Exp $ /******************************************************************************** ## Copyright (C) 1999 Thomas Walter ## ## This program is free software ## ## 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. ## ## Author: Thomas Walter ## Keywords: math fourier fft ifft FFTW ## Created: May 1999 ********************************************************************************/ #include #include /* Global variables */ static int wisdom_flags; static int Oct_cmplx_to_fftw_cmplx (fftw_complex *cv_fftw, const ComplexColumnVector & cvo) { int counter; int n_points = cvo.length (); for (counter = 0; counter < n_points; counter++) { c_re(cv_fftw[counter]) = real (cvo(counter)); c_im(cv_fftw[counter]) = imag (cvo(counter)); } return (0); } static int fftw_cmplx_to_Oct_cmplx (ComplexColumnVector & cvo, const fftw_complex *cv_fftw) { int counter; int n_points = cvo.length (); for (counter = 0; counter < n_points; counter++) { cvo(counter) = Complex (c_re(cv_fftw[counter]), c_im(cv_fftw[counter])); } return (0); } // Always generate wisdom. Uncomment this line if you don't want it. #define GENERATE_WISDOM 1 static int do_fftw (ComplexColumnVector & cvo_out, const ComplexColumnVector & cvo_in, const fftw_direction direction) { int rv; fftw_complex *z_fftw_in, *z_fftw_out; static fftw_plan plan_forward, plan_backward; static int plan_forward_size = -1; static int plan_backward_size = -1; int n_points = cvo_in.length (); char wisdom_file_name[1024]; #if !defined(GENERATE_WISDOM) wisdom_flags = FFTW_ESTIMATE | FFTW_USE_WISDOM; /* Always use wisdom */ // To generate wisdom replace FFTW_ESTIMATE with FFTW_MEASURE #else // defined(GENERATE_WISDOM) wisdom_flags = FFTW_MEASURE | FFTW_USE_WISDOM; /* Always use wisdom */ #endif // defined(GENERATE_WISDOM) { int be_quiet = 1; char *users_home = getenv("HOME"); sprintf (wisdom_file_name, "%s/%s", users_home, "wisdom"); FILE *wf = fopen (wisdom_file_name, "r"); if (wf == NULL) { if (!be_quiet) { printf ("Couldn't open wisdom file \"%s\".\n", wisdom_file_name); printf ("This file will be created upon completion.\n"); } } else { if (!(FFTW_SUCCESS == fftw_import_wisdom_from_file (wf))) { warning ("do_fftw: illegal wisdom file format. Use of wisdom switched off."); wisdom_flags = FFTW_ESTIMATE; fclose (wf); } } } // // Create plan specific to the direction // switch (direction) { case FFTW_FORWARD: if (plan_forward_size != n_points) { if (plan_forward_size > 0) { fftw_destroy_plan (plan_forward); plan_forward_size = -1; } plan_forward = fftw_create_plan (n_points, direction, wisdom_flags); plan_forward_size = n_points; if (plan_forward == NULL) { error ("do_fftw: Error creating forward plan\n"); return (1); } } break; case FFTW_BACKWARD: if (plan_backward_size != n_points) { if (plan_backward_size > 0) { fftw_destroy_plan (plan_backward); plan_backward_size = -1; } plan_backward = fftw_create_plan (n_points, direction, wisdom_flags); plan_backward_size = n_points; if (plan_backward == NULL) { error ("do_fftw: Error creating forward plan\n"); return (1); } } break; default: error ("do_fftw: Some requested an unknown direction\n"); return (1); break; } z_fftw_in = (fftw_complex *) malloc (n_points * sizeof (fftw_complex)); z_fftw_out = (fftw_complex *) malloc (n_points * sizeof (fftw_complex)); if (z_fftw_in == NULL) { return (1); } if (z_fftw_out == NULL) { return (1); } rv = Oct_cmplx_to_fftw_cmplx (z_fftw_in, cvo_in); /* * do the transformation */ switch (direction) { case FFTW_FORWARD: fftw_one (plan_forward, z_fftw_in, z_fftw_out); /* fftw (plan, 1, in, 1, 1, out, 1, 1); */ break; case FFTW_BACKWARD: fftw_one (plan_backward, z_fftw_in, z_fftw_out); /* fftw (plan, 1, in, 1, 1, out, 1, 1); */ } /* * destroy the plan, I do no further transformations */ #if defined (FFTW_ALWAYS_DESTROY_PLAN) switch (direction) { case FFTW_FORWARD: fftw_destroy_plan (plan_forward); plan_forward_size = -1; break; case FFTW_BACKWARD: fftw_destroy_plan (plan_backward); plan_backward_size = -1; break; } #endif // defined (FFTW_ALWAYS_DESTROY_PLAN) rv = fftw_cmplx_to_Oct_cmplx (cvo_out, z_fftw_out); free (z_fftw_out); free (z_fftw_in); #if defined(GENERATE_WISDOM) if (wisdom_flags & FFTW_MEASURE) { FILE *wf = fopen (wisdom_file_name, "w"); CHECK (wf != 0,"error creating wisdom file"); fftw_export_wisdom_to_file (wf); fclose (wf); } #endif // defined(GENERATE_WISDOM) return (0); } static int do_scale (const ComplexColumnVector & cvo_in, ComplexColumnVector & cvo_out, const double scale) { int counter; for (counter = 0; counter < cvo_in.length (); counter++) { cvo_out(counter) = cvo_in(counter) / scale; } return (0); } DEFUN_DLD (fftw, args, /* nargout not used */, "fftw (z): fast fourier transformation of a column vector") { fftw_direction direction = FFTW_FORWARD; // FFTW_BACKWARD octave_value_list retval; int rv = 0; int nargin = args.length (); if (nargin != 1) { print_usage ("fftw (z)"); return retval; } octave_value arg = args (0); int n_points = arg.rows (); if (n_points == 1) { n_points = arg.columns (); } ComplexColumnVector z_oct_in = arg.complex_vector_value (); ComplexColumnVector z_oct_out; z_oct_out.resize (n_points, 1); rv = do_fftw (z_oct_out, z_oct_in, direction); if (rv) { error ("fftw: transformation failed"); return (retval); } return (octave_value (z_oct_out)); } DEFUN_DLD (ifftw, args, /* nargout not used */, "ifftw (z): fast fourier transformation of a column vector") { fftw_direction direction = FFTW_BACKWARD; // FFTW_FORWARD octave_value_list retval; int rv = 0; int nargin = args.length (); if (nargin != 1) { print_usage ("ifftw (z)"); return retval; } octave_value arg = args (0); int n_points = arg.rows (); if (n_points == 1) { n_points = arg.columns (); } ComplexColumnVector z_oct_in = arg.complex_vector_value (); ComplexColumnVector z_oct_out; z_oct_out.resize (n_points, 1); rv = do_fftw (z_oct_out, z_oct_in, direction); if (rv) { error ("ifftw: transformation failed"); return (retval); } // Inverse direction needs scaling. Devide with 'n_points'. do_scale (z_oct_out, z_oct_out, n_points); // it is save to use the same variable. return (ComplexMatrix(z_oct_out)); } =========== fftw.cc end ==================== -- Platzangst: Der dauerhafte Zustand eines Luftballons. ---------------------------------------------- Dipl. Phys. Thomas Walter Inst. f. Physiklische Chemie II Egerlandstr. 3 Tel.: ++9131-85 27326 / 27330 91058 Erlangen, Germany email: walter at pctc dot chemie dot uni-erlangen dot de