git.gag.com Git - debian/gnuradio/blob - gcell/lib/wrapper/gcp_fft_1d_r2.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 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 along
  18  * with this program; if not, write to the Free Software Foundation, Inc.,
  19  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  20  */
  21 #ifdef HAVE_CONFIG_H
  22 #include <config.h>
  23 #endif
  24 #include <gcell/gcp_fft_1d_r2.h>
  25 #include <stdint.h>
  26 #include <stdexcept>
  27 #include <math.h>
  28
  29 static void
  30 init_jd(gc_job_desc *jd,
  31         gc_proc_id_t proc_id,
  32         unsigned log2_fft_length,
  33         bool shift,
  34         std::complex<float> *out,
  35         const std::complex<float> *in,
  36         const std::complex<float> *twiddle,
  37         const float *window)
  38 {
  39   jd->proc_id = proc_id;
  40   jd->input.nargs = 2;
  41   jd->output.nargs = 0;
  42   jd->eaa.nargs = 4;
  43
  44   jd->input.arg[0].u32 = log2_fft_length;
  45   jd->input.arg[1].u32 = shift;
  46   unsigned int fft_length = 1 << log2_fft_length;
  47
  48   jd->eaa.arg[0].ea_addr = ptr_to_ea(out);
  49   jd->eaa.arg[0].direction = GCJD_DMA_PUT;
  50   jd->eaa.arg[0].put_size = sizeof(std::complex<float>) * fft_length;
  51
  52   jd->eaa.arg[1].ea_addr = ptr_to_ea(const_cast<std::complex<float>*>(in));
  53   jd->eaa.arg[1].direction = GCJD_DMA_GET;
  54   jd->eaa.arg[1].get_size = sizeof(std::complex<float>) * fft_length;
  55
  56   jd->eaa.arg[2].ea_addr = ptr_to_ea(const_cast<std::complex<float>*>(twiddle));
  57   jd->eaa.arg[2].direction = GCJD_DMA_GET;
  58   jd->eaa.arg[2].get_size = sizeof(std::complex<float>) * fft_length / 4;
  59
  60   jd->eaa.arg[3].ea_addr = ptr_to_ea(const_cast<float*>(window));
  61   jd->eaa.arg[3].direction = GCJD_DMA_GET;
  62   if (window == 0)
  63     jd->eaa.arg[3].get_size = 0;
  64   else
  65     jd->eaa.arg[3].get_size = sizeof(float) * fft_length;
  66 }
  67
  68
  69 gc_job_desc_sptr
  70 gcp_fft_1d_r2_submit(gc_job_manager_sptr mgr,
  71                      unsigned int log2_fft_length,
  72                      bool forward,
  73                      bool shift,
  74                      std::complex<float> *out,
  75                      const std::complex<float> *in,
  76                      const std::complex<float> *twiddle,
  77                      const float *window)
  78 {
  79   unsigned int fft_length = 1 << log2_fft_length;
  80   if (fft_length > 4096)
  81     throw std::invalid_argument("fft_length > 4096");
  82
  83   if ((intptr_t)out & 0xf)
  84     throw gc_bad_align("out");
  85   if ((intptr_t)in & 0xf)
  86     throw gc_bad_align("in");
  87   if ((intptr_t)twiddle & 0xf)
  88     throw gc_bad_align("twiddle");
  89   if ((intptr_t)window & 0xf)
  90     throw gc_bad_align("window");
  91
  92   std::string proc_name;
  93   if (forward)
  94     proc_name = "fwd_fft_1d_r2";
  95   else
  96     proc_name = "inv_fft_1d_r2";
  97
  98   gc_proc_id_t fft_id = mgr->lookup_proc(proc_name);
  99   gc_job_desc_sptr jd = gc_job_manager::alloc_job_desc(mgr);
 100   init_jd(jd.get(), fft_id, log2_fft_length, shift, out, in, twiddle, window);
 101   if (!mgr->submit_job(jd.get())){
 102     gc_job_status_t s = jd->status;
 103     throw gc_bad_submit(proc_name, s);
 104   }
 105   return jd;
 106 }
 107
 108 void
 109 gcp_fft_1d_r2_twiddle(unsigned int log2_fft_length, std::complex<float> *twiddle)
 110 {
 111   unsigned int n = 1 << log2_fft_length;
 112
 113   twiddle[0].real() = 1.0;
 114   twiddle[0].imag() = 0.0;
 115   for (unsigned i=1; i < n/4; i++){
 116     twiddle[i].real() =  cos(i * 2*M_PI/n);
 117     twiddle[n/4 - i].imag() = -twiddle[i].real();
 118   }
 119 }