FFT inverse tranform is now correct.

git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@8211 221aa14e-8319-0410-a670-987f0aec2ac5

4 files changed, 27 insertions, 35 deletions

diff --git a/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc b/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc
index d639dad45..f92ee42c8 100644
--- a/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc
+++ b/gcell/src/lib/wrapper/gcp_fft_1d_r2.cc
@@ -89,7 +89,7 @@ gcp_fft_1d_r2_submit(gc_job_manager_sptr mgr,
 }
 
 void
-gcp_fft_1d_r2_forward_twiddle(unsigned int log2_fft_length, std::complex<float> *W)
+gcp_fft_1d_r2_twiddle(unsigned int log2_fft_length, std::complex<float> *W)
 {
   unsigned int n = 1 << log2_fft_length;
 
@@ -100,17 +100,3 @@ gcp_fft_1d_r2_forward_twiddle(unsigned int log2_fft_length, std::complex<float>
     W[n/4 - i].imag() = -W[i].real();
   }
 }
-
-
-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);
-
-  unsigned int n = 1 << log2_fft_length;
-  for (unsigned i=0; i < n/4; i++)
-    W[i] = conj(W[i]);
-}
diff --git a/gcell/src/lib/wrapper/gcp_fft_1d_r2.h b/gcell/src/lib/wrapper/gcp_fft_1d_r2.h
index be1440fd4..ed1d9e783 100644
--- a/gcell/src/lib/wrapper/gcp_fft_1d_r2.h
+++ b/gcell/src/lib/wrapper/gcp_fft_1d_r2.h
@@ -28,7 +28,7 @@
  * \brief Submit a job that computes the forward or reverse FFT.
  *
  * \param mgr is the job manager instance
- * \param log2_fft_length is the log2 of the fft_length (4 <= x <= 13).
+ * \param log2_fft_length is the log2 of the fft_length (4 <= x <= 12).
  * \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).
@@ -46,21 +46,12 @@ gcp_fft_1d_r2_submit(gc_job_manager_sptr mgr,
 		     const std::complex<float> *W);
 
 /*!
- * \brief Compute twiddle factors for forward transform.
+ * \brief Compute twiddle factors 
  *
  * \param log2_fft_length is the log2 of the fft_length.
  * \param W is fft_length/4 twiddle factor output (must be 16-byte aligned).
  */
 void
-gcp_fft_1d_r2_forward_twiddle(unsigned int log2_fft_length, std::complex<float> *W);
-
-/*!
- * \brief Compute twiddle factors for reverse transform.
- *
- * \param log2_fft_length is the log2 of the fft_length.
- * \param W is fft_length/4 twiddle factor output (must be 16-byte aligned).
- */
-void
-gcp_fft_1d_r2_reverse_twiddle(unsigned int log2_fft_length, std::complex<float> *W);
+gcp_fft_1d_r2_twiddle(unsigned int log2_fft_length, std::complex<float> *W);
 
 #endif /* INCLUDED_GCP_FFT_1D_R2_H */
diff --git a/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc b/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc
index 1bb676ac2..b177edede 100644
--- a/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc
+++ b/gcell/src/lib/wrapper/qa_gcp_fft_1d_r2.cc
@@ -150,10 +150,7 @@ qa_gcp_fft_1d_r2::test(gc_job_manager_sptr mgr, int log2_fft_size, bool forward)
   std::complex<float> *cell_out = (std::complex<float> *) cell_out_void.get();
   std::complex<float> *cell_twiddle = (std::complex<float> *) cell_twiddle_void.get();
 
-  if (forward)
-    gcp_fft_1d_r2_forward_twiddle(log2_fft_size, cell_twiddle);
-  else
-    gcp_fft_1d_r2_reverse_twiddle(log2_fft_size, cell_twiddle);
+  gcp_fft_1d_r2_twiddle(log2_fft_size, cell_twiddle);
 
   srandom(1);		// we want reproducibility
 
@@ -206,6 +203,5 @@ qa_gcp_fft_1d_r2::test(gc_job_manager_sptr mgr, int log2_fft_size, bool forward)
   fprintf(stdout, "%s fft_size = %4d  max_rel_error = %e\n",
 	  forward ? "fwd" : "rev", fft_size, max_rel);
 
-  // CPPUNIT_ASSERT(max_rel <= 1e-4);
-
+  CPPUNIT_ASSERT(max_rel <= 5e-3);
 }
diff --git a/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c b/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c
index 36bd878ed..bf4bdfd20 100644
--- a/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c
+++ b/gcell/src/lib/wrapper/spu/gcs_fft_1d_r2.c
@@ -22,6 +22,18 @@
 #include <gc_declare_proc.h>
 #include <libfft.h>
 
+/*
+ * v is really vector complex<float>
+ */
+static void
+conjugate_vector(vector float *v, int nelements)
+{
+  vector float k = {1, -1, 1, -1};
+  int i;
+  for (i = 0; i < nelements; i++)
+    v[i] *= k;
+}
+
 static void
 gcs_fft_1d_r2(const gc_job_direct_args_t *input,
 	      gc_job_direct_args_t *output __attribute__((unused)),
@@ -31,9 +43,16 @@ gcs_fft_1d_r2(const gc_job_direct_args_t *input,
   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)
+  int forward = input->arg[1].u32;	// non-zero if forward xform
 
-  fft_1d_r2(out, in, W, log2_fft_length);
+  if (forward){
+    fft_1d_r2(out, in, W, log2_fft_length);
+  }
+  else {
+    conjugate_vector(in, 1 << (log2_fft_length - 1));
+    fft_1d_r2(out, in, W, log2_fft_length);
+    conjugate_vector(out, 1 << (log2_fft_length - 1));
+  }
 }
 
 GC_DECLARE_PROC(gcs_fft_1d_r2, "fft_1d_r2");