git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/general/gri_fft.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2003,2008 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 #include <gri_fft.h>
  24 #include <fftw3.h>
  25 #include <gr_complex.h>
  26 #include <stdlib.h>
  27 #include <string.h>
  28 #include <stdio.h>
  29 #include <cassert>
  30 #include <stdexcept>
  31
  32
  33 boost::mutex &
  34 gri_fft_planner::mutex()
  35 {
  36   static boost::mutex  s_planning_mutex;
  37
  38   return s_planning_mutex;
  39 }
  40
  41 static char *
  42 wisdom_filename ()
  43 {
  44   static const char *filename = ".gr_fftw_wisdom";
  45
  46   char  *home = getenv ("HOME");
  47   if (home){
  48     char *p = new char[strlen (home) + strlen (filename) + 2];
  49     strcpy (p, home);
  50     strcat (p, "/");
  51     strcat (p, filename);
  52     return p;
  53   }
  54   return 0;
  55 }
  56
  57 static void
  58 gri_fftw_import_wisdom ()
  59 {
  60   char *filename = wisdom_filename ();
  61   FILE *fp = fopen (filename, "r");
  62   if (fp != 0){
  63     int r = fftwf_import_wisdom_from_file (fp);
  64     fclose (fp);
  65     if (!r){
  66       fprintf (stderr, "gri_fftw: can't import wisdom from %s\n", filename);
  67     }
  68   }
  69   delete [] filename;
  70 }
  71
  72 static void
  73 gri_fftw_export_wisdom ()
  74 {
  75   char *filename = wisdom_filename ();
  76   FILE *fp = fopen (filename, "w");
  77   if (fp != 0){
  78     fftwf_export_wisdom_to_file (fp);
  79     fclose (fp);
  80   }
  81   else {
  82     fprintf (stderr, "gri_fftw: ");
  83     perror (filename);
  84   }
  85   delete [] filename;
  86 }
  87
  88 // ----------------------------------------------------------------
  89
  90 gri_fft_complex::gri_fft_complex (int fft_size, bool forward)
  91 {
  92   // Hold global mutex during plan construction and destruction.
  93   gri_fft_planner::scoped_lock  lock(gri_fft_planner::mutex());
  94
  95   assert (sizeof (fftwf_complex) == sizeof (gr_complex));
  96
  97   if (fft_size <= 0)
  98     throw std::out_of_range ("gri_fftw: invalid fft_size");
  99
 100   d_fft_size = fft_size;
 101   d_inbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * inbuf_length ());
 102   if (d_inbuf == 0)
 103     throw std::runtime_error ("fftwf_malloc");
 104
 105   d_outbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * outbuf_length ());
 106   if (d_outbuf == 0){
 107     fftwf_free (d_inbuf);
 108     throw std::runtime_error ("fftwf_malloc");
 109   }
 110
 111   gri_fftw_import_wisdom ();    // load prior wisdom from disk
 112   d_plan = fftwf_plan_dft_1d (fft_size,
 113                               reinterpret_cast<fftwf_complex *>(d_inbuf),
 114                               reinterpret_cast<fftwf_complex *>(d_outbuf),
 115                               forward ? FFTW_FORWARD : FFTW_BACKWARD,
 116                               FFTW_MEASURE);
 117
 118   if (d_plan == NULL) {
 119     fprintf(stderr, "gri_fft_complex: error creating plan\n");
 120     throw std::runtime_error ("fftwf_plan_dft_1d failed");
 121   }
 122   gri_fftw_export_wisdom ();    // store new wisdom to disk
 123 }
 124
 125 gri_fft_complex::~gri_fft_complex ()
 126 {
 127   // Hold global mutex during plan construction and destruction.
 128   gri_fft_planner::scoped_lock  lock(gri_fft_planner::mutex());
 129
 130   fftwf_destroy_plan ((fftwf_plan) d_plan);
 131   fftwf_free (d_inbuf);
 132   fftwf_free (d_outbuf);
 133 }
 134
 135 void
 136 gri_fft_complex::execute ()
 137 {
 138   fftwf_execute ((fftwf_plan) d_plan);
 139 }
 140
 141 // ----------------------------------------------------------------
 142
 143 gri_fft_real_fwd::gri_fft_real_fwd (int fft_size)
 144 {
 145   // Hold global mutex during plan construction and destruction.
 146   gri_fft_planner::scoped_lock  lock(gri_fft_planner::mutex());
 147
 148   assert (sizeof (fftwf_complex) == sizeof (gr_complex));
 149
 150   if (fft_size <= 0)
 151     throw std::out_of_range ("gri_fftw: invalid fft_size");
 152
 153   d_fft_size = fft_size;
 154   d_inbuf = (float *) fftwf_malloc (sizeof (float) * inbuf_length ());
 155   if (d_inbuf == 0)
 156     throw std::runtime_error ("fftwf_malloc");
 157
 158   d_outbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * outbuf_length ());
 159   if (d_outbuf == 0){
 160     fftwf_free (d_inbuf);
 161     throw std::runtime_error ("fftwf_malloc");
 162   }
 163
 164   gri_fftw_import_wisdom ();    // load prior wisdom from disk
 165   d_plan = fftwf_plan_dft_r2c_1d (fft_size,
 166                                   d_inbuf,
 167                                   reinterpret_cast<fftwf_complex *>(d_outbuf),
 168                                   FFTW_MEASURE);
 169
 170   if (d_plan == NULL) {
 171     fprintf(stderr, "gri_fft_real_fwd: error creating plan\n");
 172     throw std::runtime_error ("fftwf_plan_dft_r2c_1d failed");
 173   }
 174   gri_fftw_export_wisdom ();    // store new wisdom to disk
 175 }
 176
 177 gri_fft_real_fwd::~gri_fft_real_fwd ()
 178 {
 179   // Hold global mutex during plan construction and destruction.
 180   gri_fft_planner::scoped_lock  lock(gri_fft_planner::mutex());
 181
 182   fftwf_destroy_plan ((fftwf_plan) d_plan);
 183   fftwf_free (d_inbuf);
 184   fftwf_free (d_outbuf);
 185 }
 186
 187 void
 188 gri_fft_real_fwd::execute ()
 189 {
 190   fftwf_execute ((fftwf_plan) d_plan);
 191 }
 192
 193 // ----------------------------------------------------------------
 194
 195 gri_fft_real_rev::gri_fft_real_rev (int fft_size)
 196 {
 197   // Hold global mutex during plan construction and destruction.
 198   gri_fft_planner::scoped_lock  lock(gri_fft_planner::mutex());
 199
 200   assert (sizeof (fftwf_complex) == sizeof (gr_complex));
 201
 202   if (fft_size <= 0)
 203     throw std::out_of_range ("gri_fftw: invalid fft_size");
 204
 205   d_fft_size = fft_size;
 206   d_inbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * inbuf_length ());
 207   if (d_inbuf == 0)
 208     throw std::runtime_error ("fftwf_malloc");
 209
 210   d_outbuf = (float *) fftwf_malloc (sizeof (float) * outbuf_length ());
 211   if (d_outbuf == 0){
 212     fftwf_free (d_inbuf);
 213     throw std::runtime_error ("fftwf_malloc");
 214   }
 215
 216   // FIXME If there's ever a chance that the planning functions
 217   // will be called in multiple threads, we've got to ensure single
 218   // threaded access.  They are not thread-safe.
 219
 220   gri_fftw_import_wisdom ();    // load prior wisdom from disk
 221   d_plan = fftwf_plan_dft_c2r_1d (fft_size,
 222                                   reinterpret_cast<fftwf_complex *>(d_inbuf),
 223                                   d_outbuf,
 224                                   FFTW_MEASURE);
 225
 226   if (d_plan == NULL) {
 227     fprintf(stderr, "gri_fft_real_rev: error creating plan\n");
 228     throw std::runtime_error ("fftwf_plan_dft_c2r_1d failed");
 229   }
 230   gri_fftw_export_wisdom ();    // store new wisdom to disk
 231 }
 232
 233 gri_fft_real_rev::~gri_fft_real_rev ()
 234 {
 235   fftwf_destroy_plan ((fftwf_plan) d_plan);
 236   fftwf_free (d_inbuf);
 237   fftwf_free (d_outbuf);
 238 }
 239
 240 void
 241 gri_fft_real_rev::execute ()
 242 {
 243   fftwf_execute ((fftwf_plan) d_plan);
 244 }
 245