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

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