git.gag.com Git - debian/gnuradio/blob - gcell/lib/wrapper/qa_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
  22 #include "qa_gcp_fft_1d_r2.h"
  23 #include <cppunit/TestAssert.h>
  24 #include <gcell/gcp_fft_1d_r2.h>
  25 #include <fftw3.h>
  26 #include <stdio.h>
  27 #include <stdlib.h>     // random, posix_memalign
  28 #include <algorithm>
  29 #include <string.h>
  30
  31 typedef boost::shared_ptr<void> void_sptr;
  32
  33 // handle to embedded SPU executable
  34 extern spe_program_handle_t gcell_all_spx;
  35
  36 /*
  37  * Return pointer to cache-aligned chunk of storage of size size bytes.
  38  * Throw if can't allocate memory.  The storage should be freed
  39  * with "free" when done.  The memory is initialized to zero.
  40  */
  41 static void *
  42 aligned_alloc(size_t size, size_t alignment = 128)
  43 {
  44   void *p = 0;
  45   if (posix_memalign(&p, alignment, size) != 0){
  46     perror("posix_memalign");
  47     throw std::runtime_error("memory");
  48   }
  49   memset(p, 0, size);           // zero the memory
  50   return p;
  51 }
  52
  53 class free_deleter {
  54 public:
  55   void operator()(void *p) {
  56     free(p);
  57   }
  58 };
  59
  60 static boost::shared_ptr<void>
  61 aligned_alloc_sptr(size_t size, size_t alignment = 128)
  62 {
  63   return boost::shared_ptr<void>(aligned_alloc(size, alignment), free_deleter());
  64 }
  65
  66 // test forward FFT
  67 void
  68 qa_gcp_fft_1d_r2::t1()
  69 {
  70   gc_jm_options opts;
  71   opts.program_handle = gc_program_handle_from_address(&gcell_all_spx);
  72   opts.nspes = 1;
  73   gc_job_manager_sptr mgr = gc_make_job_manager(&opts);
  74
  75 #if 1
  76   for (int log2_fft_size = 5; log2_fft_size <= 12; log2_fft_size++){
  77     test(mgr, log2_fft_size, true);
  78   }
  79 #else
  80   test(mgr, 5, true);
  81 #endif
  82 }
  83
  84 // test inverse FFT
  85 void
  86 qa_gcp_fft_1d_r2::t2()
  87 {
  88   gc_jm_options opts;
  89   opts.program_handle = gc_program_handle_from_address(&gcell_all_spx);
  90   opts.nspes = 1;
  91   gc_job_manager_sptr mgr = gc_make_job_manager(&opts);
  92
  93 #if 1
  94   for (int log2_fft_size = 5; log2_fft_size <= 12; log2_fft_size++){
  95     test(mgr, log2_fft_size, false);
  96   }
  97 #else
  98   test(mgr, 5, false);
  99 #endif
 100 }
 101
 102 void
 103 qa_gcp_fft_1d_r2::t3()
 104 {
 105   // FIXME Test fwd and inv with windowing option
 106 }
 107
 108 void
 109 qa_gcp_fft_1d_r2::t4()
 110 {
 111   // FIXME Test fwd and inv with shift option
 112 }
 113
 114 static inline float
 115 abs_diff(std::complex<float> x, std::complex<float> y)
 116 {
 117   return std::max(std::abs(x.real()-y.real()),
 118                   std::abs(x.imag()-y.imag()));
 119 }
 120
 121 static float
 122 float_abs_rel_error(float ref, float actual)
 123 {
 124   float delta = ref - actual;
 125   if (std::abs(ref) < 1e-18)
 126     ref = 1e-18;
 127   return std::abs(delta/ref);
 128 }
 129
 130 static float
 131 abs_rel_error(std::complex<float> ref, std::complex<float> actual)
 132 {
 133   return std::max(float_abs_rel_error(ref.real(), actual.real()),
 134                   float_abs_rel_error(ref.imag(), actual.imag()));
 135 }
 136
 137 void
 138 qa_gcp_fft_1d_r2::test(gc_job_manager_sptr mgr, int log2_fft_size, bool forward)
 139 {
 140   int fft_size = 1 << log2_fft_size;
 141
 142   // allocate aligned buffers with boost shared_ptr's
 143   void_sptr fftw_in_void = aligned_alloc_sptr(fft_size * sizeof(std::complex<float>), 128);
 144   void_sptr fftw_out_void = aligned_alloc_sptr(fft_size * sizeof(std::complex<float>), 128);
 145   void_sptr cell_in_void = aligned_alloc_sptr(fft_size * sizeof(std::complex<float>), 128);
 146   void_sptr cell_out_void = aligned_alloc_sptr(fft_size * sizeof(std::complex<float>), 128);
 147   void_sptr cell_twiddle_void = aligned_alloc_sptr(fft_size/4 * sizeof(std::complex<float>), 128);
 148
 149   // cast them to the type we really want
 150   std::complex<float> *fftw_in = (std::complex<float> *) fftw_in_void.get();
 151   std::complex<float> *fftw_out = (std::complex<float> *) fftw_out_void.get();
 152   std::complex<float> *cell_in = (std::complex<float> *) cell_in_void.get();
 153   std::complex<float> *cell_out = (std::complex<float> *) cell_out_void.get();
 154   std::complex<float> *cell_twiddle = (std::complex<float> *) cell_twiddle_void.get();
 155
 156   gcp_fft_1d_r2_twiddle(log2_fft_size, cell_twiddle);
 157
 158   srandom(1);           // we want reproducibility
 159
 160   // initialize the input buffers
 161   for (int i = 0; i < fft_size; i++){
 162     std::complex<float> t((float) (random() & 0xfffff), (float) (random() & 0xfffff));
 163     fftw_in[i] = t;
 164     cell_in[i] = t;
 165   }
 166
 167   // ------------------------------------------------------------------------
 168   // compute the reference answer
 169   fftwf_plan plan = fftwf_plan_dft_1d (fft_size,
 170                                        reinterpret_cast<fftwf_complex *>(fftw_in),
 171                                        reinterpret_cast<fftwf_complex *>(fftw_out),
 172                                        forward ? FFTW_FORWARD : FFTW_BACKWARD,
 173                                        FFTW_ESTIMATE);
 174   if (plan == 0){
 175     fprintf(stderr, "qa_gcp_fft_1d_r2: error creating FFTW plan\n");
 176     throw std::runtime_error ("fftwf_plan_dft_r2c_1d failed");
 177   }
 178
 179   fftwf_execute(plan);
 180   fftwf_destroy_plan(plan);
 181
 182   // ------------------------------------------------------------------------
 183   // compute the answer on the cell
 184   gc_job_desc_sptr jd = gcp_fft_1d_r2_submit(mgr, log2_fft_size, forward, false,
 185                                              cell_out, cell_in, cell_twiddle, 0);
 186   if (!mgr->wait_job(jd.get())){
 187     fprintf(stderr, "wait_job failed: %s\n", gc_job_status_string(jd->status).c_str());
 188     CPPUNIT_ASSERT(0);
 189   }
 190
 191   // ------------------------------------------------------------------------
 192   // compute the maximum of the relative error
 193   float max_rel = 0.0;
 194   for (int i = 0; i < fft_size; i++){
 195     max_rel = std::max(max_rel, abs_rel_error(fftw_out[i], cell_out[i]));
 196     if (0)
 197       printf("(%16.3f, %16.3fj)  (%16.3f, %16.3fj)  (%16.3f, %16.3fj)\n",
 198              fftw_out[i].real(), fftw_out[i].imag(),
 199              cell_out[i].real(), cell_out[i].imag(),
 200              fftw_out[i].real() - cell_out[i].real(),
 201              fftw_out[i].imag() - cell_out[i].imag());
 202   }
 203
 204   fprintf(stdout, "%s fft_size = %4d  max_rel_error = %e\n",
 205           forward ? "fwd" : "rev", fft_size, max_rel);
 206
 207   CPPUNIT_ASSERT(max_rel <= 5e-3);
 208 }