Re: [Discuss-gnuradio] NCO accuracy

[Eric Blossom]

On Thu, Jun 16, 2005 at 01:28:31AM +0200, Stephane Fillod wrote:
> 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 ();

Thanks!

> 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).

Committed to CVS.  

Note that on AMD64, default floating point code generation uses scalar
SSE operations which are 32-bit for 32-bit ops.  Also, note that
single precision is twice as fast as double precision on AMD64 using
scalar SSE.  AMD64 can perform 3.8 FLOPS per cycle!  (~ 2 FADDs + ~ 2
FMULTs).  Using SSE gets you twice as many registers as x87 too.

> 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" :-)

This is all good.  What we're running into is an issue of how much
precision to carry in the phase accumulator and what is the "right
answer".  In an ideal world, we wouldn't be using floats, we'd 
use exact rationals.  Note that gr_fxpt_nco is using 32-bits of
precision vs 24-bits for single precision floating point (IEEE
mantissa size).

Now, after changing gr_nco to double, it's carrying more precision
than gr_fxpt_nco.  Of course we could change gr_fxpt_nco to use 64-bit
ints and then it would have more precision...

The code below is a too hairy for this "fast path" code.  It turns
what was a single cycle operation into something on the order of 100
cycles and also includes conditional branching.  Also note that the
real problem is that the frequency arg to set_freq is float not
double.  Our existing 32-bit fixed point accumlator can represent the
32-bit float freq value with no loss of precision.

>    void step () 
>    { 
>      d_phase += d_phase_inc;
> +
> +    ++d_step;
> +    if (d_step >= d_period) {
> +         d_phase_origin -= d_period_frac;
> +         d_phase = d_phase_origin;
> +         d_step = 0;
> +    } else {
> +         if ((d_step % STEP_CHECK) == 0) {
> +             d_phase += d_step_error;
> +         }
> +    }
>    }

Good work shaking all this down.

If you're so inclined, I would entertain a 64-bit fixed point version.
The 64-bit phase accumulator update will be on the order of 1 to 3
cycles depending on architecture.

Eric