[Commit-gnuradio] r8200 - in gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper: . spu

[eb]

Author: eb
Date: 2008-04-15 12:50:47 -0600 (Tue, 15 Apr 2008)
New Revision: 8200

Modified:
   
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc
   
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.h
   
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc
   
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c
Log:
gcell work-in-progress

Modified: 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc
===================================================================
--- 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc
    2008-04-15 17:12:18 UTC (rev 8199)
+++ 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc
    2008-04-15 18:50:47 UTC (rev 8200)
@@ -30,16 +30,18 @@
 init_jd(gc_job_desc *jd,
        gc_proc_id_t proc_id,
        unsigned log2_fft_length,
+       bool forward,
        std::complex<float> *out,
        const std::complex<float> *in,
        const std::complex<float> *W)
 {
   jd->proc_id = proc_id;
-  jd->input.nargs = 1;
+  jd->input.nargs = 2;
   jd->output.nargs = 0;
   jd->eaa.nargs = 3;
 
   jd->input.arg[0].u32 = log2_fft_length;
+  jd->input.arg[1].u32 = forward;
   unsigned int fft_length = 1 << log2_fft_length;
 
   jd->eaa.arg[0].ea_addr = ptr_to_ea(out);
@@ -59,6 +61,7 @@
 gc_job_desc *
 gcp_fft_1d_r2_submit(gc_job_manager_sptr mgr,
                     unsigned int log2_fft_length,
+                    bool forward,
                     std::complex<float> *out,
                     const std::complex<float> *in,
                     const std::complex<float> *W)
@@ -76,7 +79,7 @@
 
   gc_proc_id_t fft_id = mgr->lookup_proc("fft_1d_r2");
   gc_job_desc *jd = mgr->alloc_job_desc();
-  init_jd(jd, fft_id, log2_fft_length, out, in, W);
+  init_jd(jd, fft_id, log2_fft_length, forward, out, in, W);
   if (!mgr->submit_job(jd)){
     gc_job_status_t s = jd->status;
     mgr->free_job_desc(jd);
@@ -102,6 +105,8 @@
 void
 gcp_fft_1d_r2_reverse_twiddle(unsigned int log2_fft_length, 
std::complex<float> *W)
 {
+  // FIXME this is wrong/insufficient.  inverse is still incorrect
+
   // reverse factors are the conjugate of the forward factors
   gcp_fft_1d_r2_forward_twiddle(log2_fft_length, W);
 

Modified: 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.h
===================================================================
--- 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.h 
    2008-04-15 17:12:18 UTC (rev 8199)
+++ 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/gcp_fft_1d_r2.h 
    2008-04-15 18:50:47 UTC (rev 8200)
@@ -29,6 +29,7 @@
  *
  * \param mgr is the job manager instance
  * \param log2_fft_length is the log2 of the fft_length (4 <= x <= 13).
+ * \param forward is true to compute the forward xform
  * \param out is the fft_length output from FFT (must be 16-byte aligned).
  * \param in is the fft_length input to FFT (must be 16-byte aligned).
  * \param W is fft_length/4 twiddle factor input to FFT (must be 16-byte 
aligned).
@@ -39,6 +40,7 @@
 gc_job_desc *
 gcp_fft_1d_r2_submit(gc_job_manager_sptr mgr,
                     unsigned int log2_fft_length,
+                    bool forward,
                     std::complex<float> *out,
                     const std::complex<float> *in,
                     const std::complex<float> *W);

Modified: 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc
===================================================================
--- 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc
 2008-04-15 17:12:18 UTC (rev 8199)
+++ 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc
 2008-04-15 18:50:47 UTC (rev 8200)
@@ -71,7 +71,7 @@
   opts.nspes = 1;
   gc_job_manager_sptr mgr = gc_make_job_manager(&opts);
 
-#if 0
+#if 1
   for (int log2_fft_size = 5; log2_fft_size <= 11; log2_fft_size++){
     test(mgr, log2_fft_size, true);
   }
@@ -89,7 +89,7 @@
   opts.nspes = 1;
   gc_job_manager_sptr mgr = gc_make_job_manager(&opts);
 
-#if 0
+#if 1
   for (int log2_fft_size = 5; log2_fft_size <= 11; log2_fft_size++){
     test(mgr, log2_fft_size, false);
   }
@@ -115,6 +115,22 @@
                  std::abs(x.imag()-y.imag()));
 }
 
+static float
+float_abs_rel_error(float ref, float actual)
+{
+  float delta = ref - actual;
+  if (std::abs(ref) < 1e-18)
+    ref = 1e-18;
+  return std::abs(delta/ref);
+}
+
+static float
+abs_rel_error(std::complex<float> ref, std::complex<float> actual)
+{
+  return std::max(float_abs_rel_error(ref.real(), actual.real()),
+                 float_abs_rel_error(ref.imag(), actual.imag()));
+}
+
 void 
 qa_gcp_fft_1d_r2::test(gc_job_manager_sptr mgr, int log2_fft_size, bool 
forward)
 {
@@ -165,7 +181,7 @@
 
   // ------------------------------------------------------------------------
   // compute the answer on the cell
-  gc_job_desc *jd = gcp_fft_1d_r2_submit(mgr, log2_fft_size,
+  gc_job_desc *jd = gcp_fft_1d_r2_submit(mgr, log2_fft_size, forward,
                                         cell_out, cell_in, cell_twiddle);
   if (!mgr->wait_job(jd)){
     fprintf(stderr, "wait_job failed: %s\n", 
gc_job_status_string(jd->status).c_str());
@@ -175,11 +191,11 @@
   mgr->free_job_desc(jd);
 
   // ------------------------------------------------------------------------
-  // compute the maximum of the abs diff of the real and imag components
-  float max_diff = 0.0;
+  // compute the maximum of the relative error
+  float max_rel = 0.0;
   for (int i = 0; i < fft_size; i++){
-    max_diff = std::max(max_diff, abs_diff(fftw_out[i], cell_out[i]));
-    if (1)
+    max_rel = std::max(max_rel, abs_rel_error(fftw_out[i], cell_out[i]));
+    if (0)
       printf("(%16.3f, %16.3fj)  (%16.3f, %16.3fj)  (%16.3f, %16.3fj)\n",
             fftw_out[i].real(), fftw_out[i].imag(),
             cell_out[i].real(), cell_out[i].imag(),
@@ -187,5 +203,9 @@
             fftw_out[i].imag() - cell_out[i].imag());
   }
 
-  fprintf(stdout, "%s fft_size = %4d  max_diff = %f\n", forward ? "fwd" : 
"rev", fft_size, max_diff);
+  fprintf(stdout, "%s fft_size = %4d  max_rel_error = %e\n",
+         forward ? "fwd" : "rev", fft_size, max_rel);
+
+  // CPPUNIT_ASSERT(max_rel <= 1e-4);
+
 }

Modified: 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c
===================================================================
--- 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c
 2008-04-15 17:12:18 UTC (rev 8199)
+++ 
gnuradio/branches/developers/eb/gcell-wip/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c
 2008-04-15 18:50:47 UTC (rev 8200)
@@ -31,6 +31,7 @@
   vector float *in = (vector float *) eaa->arg[1].ls_addr;
   vector float *W = (vector float *) eaa->arg[2].ls_addr;
   int log2_fft_length = input->arg[0].u32;
+  int forward = input->arg[1].u32;     // non-zero if forward xform (FIXME use)
 
   fft_1d_r2(out, in, W, log2_fft_length);
 }