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 #include <boost/thread.hpp>
  32
  33 typedef boost::mutex::scoped_lock scoped_lock;
  34 static boost::mutex  s_planning_mutex;
  35
  36
  37 static char *
  38 wisdom_filename ()
  39 {
  40   static char *filename = ".gr_fftw_wisdom";
  41
  42   char  *home = getenv ("HOME");
  43   if (home){
  44     char *p = new char[strlen (home) + strlen (filename) + 2];
  45     strcpy (p, home);
  46     strcat (p, "/");
  47     strcat (p, filename);
  48     return p;
  49   }
  50   return 0;
  51 }
  52
  53 static void
  54 gri_fftw_import_wisdom ()
  55 {
  56   char *filename = wisdom_filename ();
  57   FILE *fp = fopen (filename, "r");
  58   if (fp != 0){
  59     int r = fftwf_import_wisdom_from_file (fp);
  60     fclose (fp);
  61     if (!r){
  62       fprintf (stderr, "gri_fftw: can't import wisdom from %s\n", filename);
  63     }
  64   }
  65   delete [] filename;
  66 }
  67
  68 static void
  69 gri_fftw_export_wisdom ()
  70 {
  71   char *filename = wisdom_filename ();
  72   FILE *fp = fopen (filename, "w");
  73   if (fp != 0){
  74     fftwf_export_wisdom_to_file (fp);
  75     fclose (fp);
  76   }
  77   else {
  78     fprintf (stderr, "gri_fftw: ");
  79     perror (filename);
  80   }
  81   delete [] filename;
  82 }
  83
  84 // ----------------------------------------------------------------
  85
  86 gri_fft_complex::gri_fft_complex (int fft_size, bool forward)
  87 {
  88   // Hold global mutex during plan construction and destruction.
  89   scoped_lock   lock(s_planning_mutex);
  90
  91   assert (sizeof (fftwf_complex) == sizeof (gr_complex));
  92
  93   if (fft_size <= 0)
  94     throw std::out_of_range ("gri_fftw: invalid fft_size");
  95
  96   d_fft_size = fft_size;
  97   d_inbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * inbuf_length ());
  98   if (d_inbuf == 0)
  99     throw std::runtime_error ("fftwf_malloc");
 100
 101   d_outbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * outbuf_length ());
 102   if (d_outbuf == 0){
 103     fftwf_free (d_inbuf);
 104     throw std::runtime_error ("fftwf_malloc");
 105   }
 106
 107   gri_fftw_import_wisdom ();    // load prior wisdom from disk
 108   d_plan = fftwf_plan_dft_1d (fft_size,
 109                               reinterpret_cast<fftwf_complex *>(d_inbuf),
 110                               reinterpret_cast<fftwf_complex *>(d_outbuf),
 111                               forward ? FFTW_FORWARD : FFTW_BACKWARD,
 112                               FFTW_MEASURE);
 113
 114   if (d_plan == NULL) {
 115     fprintf(stderr, "gri_fft_complex: error creating plan\n");
 116     throw std::runtime_error ("fftwf_plan_dft_1d failed");
 117   }
 118   gri_fftw_export_wisdom ();    // store new wisdom to disk
 119 }
 120
 121 gri_fft_complex::~gri_fft_complex ()
 122 {
 123   // Hold global mutex during plan construction and destruction.
 124   scoped_lock   lock(s_planning_mutex);
 125
 126   fftwf_destroy_plan ((fftwf_plan) d_plan);
 127   fftwf_free (d_inbuf);
 128   fftwf_free (d_outbuf);
 129 }
 130
 131 void
 132 gri_fft_complex::execute ()
 133 {
 134   fftwf_execute ((fftwf_plan) d_plan);
 135 }
 136
 137 // ----------------------------------------------------------------
 138
 139 gri_fft_real_fwd::gri_fft_real_fwd (int fft_size)
 140 {
 141   // Hold global mutex during plan construction and destruction.
 142   scoped_lock   lock(s_planning_mutex);
 143
 144   assert (sizeof (fftwf_complex) == sizeof (gr_complex));
 145
 146   if (fft_size <= 0)
 147     throw std::out_of_range ("gri_fftw: invalid fft_size");
 148
 149   d_fft_size = fft_size;
 150   d_inbuf = (float *) fftwf_malloc (sizeof (float) * inbuf_length ());
 151   if (d_inbuf == 0)
 152     throw std::runtime_error ("fftwf_malloc");
 153
 154   d_outbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * outbuf_length ());
 155   if (d_outbuf == 0){
 156     fftwf_free (d_inbuf);
 157     throw std::runtime_error ("fftwf_malloc");
 158   }
 159
 160   gri_fftw_import_wisdom ();    // load prior wisdom from disk
 161   d_plan = fftwf_plan_dft_r2c_1d (fft_size,
 162                                   d_inbuf,
 163                                   reinterpret_cast<fftwf_complex *>(d_outbuf),
 164                                   FFTW_MEASURE);
 165
 166   if (d_plan == NULL) {
 167     fprintf(stderr, "gri_fft_real_fwd: error creating plan\n");
 168     throw std::runtime_error ("fftwf_plan_dft_r2c_1d failed");
 169   }
 170   gri_fftw_export_wisdom ();    // store new wisdom to disk
 171 }
 172
 173 gri_fft_real_fwd::~gri_fft_real_fwd ()
 174 {
 175   // Hold global mutex during plan construction and destruction.
 176   scoped_lock   lock(s_planning_mutex);
 177
 178   fftwf_destroy_plan ((fftwf_plan) d_plan);
 179   fftwf_free (d_inbuf);
 180   fftwf_free (d_outbuf);
 181 }
 182
 183 void
 184 gri_fft_real_fwd::execute ()
 185 {
 186   fftwf_execute ((fftwf_plan) d_plan);
 187 }
 188
 189 // ----------------------------------------------------------------
 190
 191 gri_fft_real_rev::gri_fft_real_rev (int fft_size)
 192 {
 193   assert (sizeof (fftwf_complex) == sizeof (gr_complex));
 194
 195   if (fft_size <= 0)
 196     throw std::out_of_range ("gri_fftw: invalid fft_size");
 197
 198   d_fft_size = fft_size;
 199   d_inbuf = (gr_complex *) fftwf_malloc (sizeof (gr_complex) * inbuf_length ());
 200   if (d_inbuf == 0)
 201     throw std::runtime_error ("fftwf_malloc");
 202
 203   d_outbuf = (float *) fftwf_malloc (sizeof (float) * outbuf_length ());
 204   if (d_outbuf == 0){
 205     fftwf_free (d_inbuf);
 206     throw std::runtime_error ("fftwf_malloc");
 207   }
 208
 209   // FIXME If there's ever a chance that the planning functions
 210   // will be called in multiple threads, we've got to ensure single
 211   // threaded access.  They are not thread-safe.
 212
 213   gri_fftw_import_wisdom ();    // load prior wisdom from disk
 214   d_plan = fftwf_plan_dft_c2r_1d (fft_size,
 215                                   reinterpret_cast<fftwf_complex *>(d_inbuf),
 216                                   d_outbuf,
 217                                   FFTW_MEASURE);
 218
 219   if (d_plan == NULL) {
 220     fprintf(stderr, "gri_fft_real_rev: error creating plan\n");
 221     throw std::runtime_error ("fftwf_plan_dft_c2r_1d failed");
 222   }
 223   gri_fftw_export_wisdom ();    // store new wisdom to disk
 224 }
 225
 226 gri_fft_real_rev::~gri_fft_real_rev ()
 227 {
 228   fftwf_destroy_plan ((fftwf_plan) d_plan);
 229   fftwf_free (d_inbuf);
 230   fftwf_free (d_outbuf);
 231 }
 232
 233 void
 234 gri_fft_real_rev::execute ()
 235 {
 236   fftwf_execute ((fftwf_plan) d_plan);
 237 }
 238