[Discuss-gnuradio] NCO accuracy

[Stephane Fillod]

Fellow hackers,

Everybody can make copy/paste bugs, like the one fixed by the piece of 
patch (against current CVS) hereafter:


RCS file: /cvsroot/gnuradio/gnuradio-core/src/lib/general/qa_gr_fxpt_nco.cc,v
retrieving revision 1.1
diff -u -r1.1 qa_gr_fxpt_nco.cc
--- src/lib/general/qa_gr_fxpt_nco.cc   19 Dec 2004 05:48:39 -0000      1.1
+++ src/lib/general/qa_gr_fxpt_nco.cc   15 Jun 2005 21:45:03 -0000
@@ -42,11 +42,11 @@
 
   for (int i = 0; i < 100000; i++){
     float ref_sin = ref_nco.sin ();
-    float new_sin = ref_nco.sin ();
+    float new_sin = new_nco.sin ();
     CPPUNIT_ASSERT_DOUBLES_EQUAL (ref_sin, new_sin, SIN_COS_TOLERANCE);
 
     float ref_cos = ref_nco.cos ();
-    float new_cos = ref_nco.cos ();
+    float new_cos = new_nco.cos ();
     CPPUNIT_ASSERT_DOUBLES_EQUAL (ref_cos, new_cos, SIN_COS_TOLERANCE);
 
     ref_nco.step ();


Brown-bag thing, isn't it?
In fact, this is an interesting bug, worth it because it allows some 
hacking fun and get to understand new stuff in the end of the day (well,
several days in my case :-).

After fixing it, I reran the test_general, and this time it failed 
in qa_gr_fxpt_nco::t0. So what, the fxpt code is not accurate?
Wrong! the fxpt code is accurate (to some degree). This is what
I found after some investigation. The equality assertion actually 
failed because of a float type. Indeed, the phase type in gr_nco.h
is a float, and after several additions, rounding errors drove
the supposed linear phase out of its way (as seen with Gnuplot).

So it looks like gr_nco.h does need to have its phase and phase_inc 
variables of double type. The output types of sin/cos stays float.
Of course, the move from float to double does not slow down benchmark_nco
on systems with single&double FPU (x87 and alike).
CVS commiters will find a patch attached which addresses this (gr_nco.patch).

Pulling the magnifying glasses, the following patch is to be prefered
over the previous on qa_gr_fxpt_nco.cc.

Index: qa_gr_fxpt_nco.cc
===================================================================
RCS file: /cvsroot/gnuradio/gnuradio-core/src/lib/general/qa_gr_fxpt_nco.cc,v
retrieving revision 1.1
diff -u -r1.1 qa_gr_fxpt_nco.cc
--- src/lib/qa_gr_fxpt_nco.cc   19 Dec 2004 05:48:39 -0000      1.1
+++ src/lib/qa_gr_fxpt_nco.cc   15 Jun 2005 22:28:02 -0000
@@ -31,27 +31,45 @@
 
 static const float SIN_COS_TOLERANCE = 1e-5;
 
+static double max_d(double a, double b)
+{
+  return fabs(a) > fabs(b) ? a : b;
+}
+
 void
 qa_gr_fxpt_nco::t0 ()
 {
   gr_nco<float,float>  ref_nco;
   gr_fxpt_nco          new_nco;
+  double max_error = 0, max_phase_error = 0;
 
   ref_nco.set_freq (2 * M_PI / 5003);
   new_nco.set_freq (2 * M_PI / 5003);
 
+  CPPUNIT_ASSERT_DOUBLES_EQUAL (ref_nco.get_freq(), new_nco.get_freq(), 
SIN_COS_TOLERANCE);
+
   for (int i = 0; i < 100000; i++){
     float ref_sin = ref_nco.sin ();
-    float new_sin = ref_nco.sin ();
+    float new_sin = new_nco.sin ();
     CPPUNIT_ASSERT_DOUBLES_EQUAL (ref_sin, new_sin, SIN_COS_TOLERANCE);
 
+    max_error = max_d (max_error, ref_sin-new_sin);
+
     float ref_cos = ref_nco.cos ();
-    float new_cos = ref_nco.cos ();
+    float new_cos = new_nco.cos ();
     CPPUNIT_ASSERT_DOUBLES_EQUAL (ref_cos, new_cos, SIN_COS_TOLERANCE);
 
+    max_error = max_d (max_error, ref_cos-new_cos);
+
     ref_nco.step ();
     new_nco.step ();
+
+    CPPUNIT_ASSERT_DOUBLES_EQUAL (ref_nco.get_phase(), new_nco.get_phase(), 
SIN_COS_TOLERANCE);
+
+    max_phase_error = max_d (max_phase_error, 
ref_nco.get_phase()-new_nco.get_phase());
   }
+
+  printf ("Fxpt  max error %.9f, max phase error %.9f\n", max_error, 
max_phase_error);
 }
 
 void


With the new patch on qa_gr_fxpt_nco.cc and the a good ref (gr_nco with
phase on double), it appears fxpt_nco now has room for improvement.
What if fxpt suffers the same problems gr_nco did?
Indeed, not only fxpt has a bit of jitter due to linear approx which we
could live with it, fxpt_nco can seriously drift from the expected
frequency due to the fixed point phase. Aha! again the phase! I told you
it was an interesting bug "thread" :-)
Again with fxpt_nco, the error add up over step(), and little by 
little, it's not so little. I haven't computed the ppm instability,
but IMHO, the float_nco and fxpt_nco without a fix should not be used
as stable Numerically Controlled Oscillators in GNU Radio.

The float nco seems to have its fix (phase of double type).
I've tried to find something for fxpt_nco along a step counter,
with cheap fixup and a return back to the phase origin every period.
This is in the attached file gr_fxpt_nco.patch. The angle_rate
and friends should be double instead of float to stand the comparison
with float_nco. Rem: angle_rate's not rounding well give interesting 
cases. The qa_gr_fxpt_nco.cc should probably also test over
longer times (more periods).

Please comment on it, I am no expert. New ideas are also welcome (I did 
not bother googling or read specific books). BTW, is this nco quest worth
it?  Am I chasing precision which is useless unless for people measuring 
uHz signals?

Anyway, it was fun and refreshing to remind the IEEE-754 format and 
its caveats (stuff like http://www.cs.unc.edu/~ibr/projects/paranoia/ )
and try to devise some workaround. The NCO fxpt might eventually
be accelerated one day (expect at least 4x speedup over current fxpt).

Cheers,
-- 
Stephane

gr_nco.patch
Description: Text document

gr_fxpt_nco.patch
Description: Text document