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