From octave-sources-request at bevo dot che dot wisc dot edu Fri Jun 11 11:35:03 1999 Subject: Re: Use 'fftw-2.1' for fast fourier transformation if you want From: Thomas Walter To: walter at pctc dot chemie dot uni-erlangen dot de CC: octave-sources at bevo dot che dot wisc dot edu Date: Fri, 11 Jun 1999 18:35:30 +0200 Hello again, here is an update of my example for fast fourier transformation of vectors using FFTW. There are some bug fixes and improvements. The code in C++ is at the end. It also contains how to compile it. This description and the location where to find 'fftw-2.1.2.tar.gz' is part of the source below. As usual, hints are welcome. Nice weekend 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.18 1999/06/11 16:23: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; // Always remember: All function share the wisdom static int Oct_cmplx_to_fftw_cmplx (fftw_complex *cv_fftw, const ComplexColumnVector & cvo) { int counter; const 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; const 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; static int plan_has_changed = 0; // flag to avoid unnecessary saving of wisdom const int n_points = cvo_in.length (); static char wisdom_file_name[1024] = "\0"; // flag also to avoid reading of wisdom file more than once #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) if (wisdom_file_name[0] == '\0') { 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); } plan_has_changed = 1; // set flag } 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 backward plan\n"); return (1); } plan_has_changed = 1; // set flag } break; default: error ("do_fftw: Someone 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); // // There could be a major improvement if one can check if there is new wisdom // #if defined(GENERATE_WISDOM) if (plan_has_changed && (wisdom_flags & FFTW_MEASURE)) { FILE *wf = fopen (wisdom_file_name, "w"); if (wf == NULL) { warning ("do_fftw: error creating wisdom file. Wisdom not saved."); return (0); } fftw_export_wisdom_to_file (wf); fclose (wf); plan_has_changed = 0; // reset flag } #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); } // // Remarks: // Doing 'ifftw (fftw (z))' gives the original data. // The function 'ifftw()' does the scaling by '1/N' ! // DEFUN_DLD (fftw, args, /* nargout not used */, "fftw (z): fast fourier transformation (FORWARD) of a vector. Its opposite is 'ifftw(z)'.\n\tNOTE: Use 'fft (1)' to free the internal table used for calculation.") { const fftw_direction direction = FFTW_FORWARD; // FFTW_BACKWARD ComplexMatrix retval; int rv = 0; const int nargin = args.length (); const octave_value arg = args (0); // May be NULL, but the sanity check below prevents its use // Sanity checks if ((nargin != 1) || ((arg.rows () != 1) && (arg.columns () != 1))) { print_usage ("fftw"); return retval; } #if 1 // This gives a larger binary (a few bytes) const int n_points = (arg.rows () == 1) ? arg.columns () : arg.rows (); #else int n_points = arg.rows (); if (n_points == 1) { n_points = arg.columns (); } #endif const 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); } // Notes: // If you want // a = fftw ([1:5]) returning a row vector // and // a = fftw ([1:5]') returning a column vector // the result of this function must be a Matrix !!! if (arg.rows () > 1) { // return (octave_value (z_oct_out)); // This returns alway a column vector return (ComplexMatrix (z_oct_out)); // Return as Column Vector } else { // ComplexRowVector bad_choice = z_oct_out.transpose(); ComplexMatrix bad_choice = z_oct_out.transpose(); return ((bad_choice)); // Return as Row Vector // return (octave_value (z_oct_out.transpose() )); // return (octave_value ((double)arg.columns())); } } DEFUN_DLD (ifftw, args, /* nargout not used */, "ifftw (z): fast fourier transformation (INVERSE) of a vector. This is the opposite of 'fftw(z)'.\n\tNOTE: Use 'ifft (1)' to free the internal table used for calculation.") { const fftw_direction direction = FFTW_BACKWARD; // FFTW_FORWARD ComplexMatrix retval; int rv = 0; const int nargin = args.length (); const octave_value arg = args (0); // May be NULL, but the sanity check below prevents its use // Sanity checks if ((nargin != 1) || ((arg.rows () != 1) && (arg.columns () != 1))) { print_usage ("ifftw"); return retval; } #if 1 const int n_points = (arg.rows () == 1) ? arg.columns () : arg.rows (); #else int n_points = arg.rows (); if (n_points == 1) { n_points = arg.columns (); } #endif const 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. if (arg.rows () > 1) { return (ComplexMatrix(z_oct_out)); // Return as Column Vector } else { ComplexMatrix bad_choice = z_oct_out.transpose(); return (bad_choice); // Return as Row Vector } } ============= 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