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
Index: gr_nco.h
===================================================================
RCS file: /cvsroot/gnuradio/gnuradio-core/src/lib/general/gr_nco.h,v
retrieving revision 1.3
diff -u -r1.3 gr_nco.h
--- src/lib/general/gr_nco.h 31 Jul 2004 22:23:25 -0000 1.3
+++ src/lib/general/gr_nco.h 15 Jun 2005 22:07:48 -0000
@@ -42,22 +42,22 @@
virtual ~gr_nco () {}
// radians
- void set_phase (float angle) {
+ void set_phase (double angle) {
phase = angle;
}
- void adjust_phase (float delta_phase) {
+ void adjust_phase (double delta_phase) {
phase += delta_phase;
}
// angle_rate is in radians / step
- void set_freq (float angle_rate){
+ void set_freq (double angle_rate){
phase_inc = angle_rate;
}
// angle_rate is a delta in radians / step
- void adjust_freq (float delta_angle_rate)
+ void adjust_freq (double delta_angle_rate)
{
phase_inc += delta_angle_rate;
}
@@ -77,11 +77,22 @@
}
}
- void step (int n) { phase += phase_inc * n; }
+ void step (int n)
+ {
+ phase += phase_inc * n;
+ if (fabs (phase) > M_PI){
+
+ while (phase > M_PI)
+ phase -= 2*M_PI;
+
+ while (phase < -M_PI)
+ phase += 2*M_PI;
+ }
+ }
// units are radians / step
- float get_phase () const { return phase; }
- float get_freq () const { return phase_inc; }
+ double get_phase () const { return phase; }
+ double get_freq () const { return phase_inc; }
// compute sin and cos for current phase angle
void sincos (float *sinx, float *cosx) const;
@@ -95,9 +104,9 @@
float sin () const { return std::sin (phase); }
protected:
- float phase;
- float phase_inc;
+ double phase;
+ double phase_inc;
};
template<class o_type, class i_type>
Index: gr_fxpt_nco.h
===================================================================
RCS file: /cvsroot/gnuradio/gnuradio-core/src/lib/general/gr_fxpt_nco.h,v
retrieving revision 1.1
diff -u -r1.1 gr_fxpt_nco.h
--- gr_fxpt_nco.h 19 Dec 2004 05:48:39 -0000 1.1
+++ gr_fxpt_nco.h 15 Jun 2005 22:57:18 -0000
@@ -28,32 +28,64 @@
* \brief Numerically Controlled Oscillator (NCO)
*/
class gr_fxpt_nco {
+ gr_int32 d_phase_origin;
gr_int32 d_phase;
gr_int32 d_phase_inc;
+ gr_int32 d_step_error;
+ gr_uint32 d_period;
+ gr_uint32 d_step;
+ gr_uint32 d_period_frac;
public:
- gr_fxpt_nco () : d_phase (0), d_phase_inc (0) {}
+ gr_fxpt_nco () : d_phase_origin (0), d_phase (0), d_phase_inc (0),
d_step_error (0), d_period (0), d_step (0), d_period_frac (0) {}
~gr_fxpt_nco () {}
// radians
void set_phase (float angle) {
- d_phase = gr_fxpt::float_to_fixed (angle);
+ d_phase_origin = gr_fxpt::float_to_fixed (angle);
+ d_phase = d_phase_origin;
}
void adjust_phase (float delta_phase) {
d_phase += gr_fxpt::float_to_fixed (delta_phase);
}
+#define STEP_CHECK (1U<<8)
+
+ // step_rate is in steps / period
+ void set_period (float step_rate){
+ float phase_inc_error;
+ float step_rate_ceil;
+ float phase_inc;
+
+ step_rate_ceil = ceilf(step_rate);
+ d_period = (gr_int32)step_rate_ceil;
+
+ phase_inc = 2*TWO_TO_THE_31 / step_rate;
+ d_phase_inc = (gr_int32)phase_inc;
+
+ d_period_frac = (gr_int32)(phase_inc*(step_rate_ceil - step_rate));
+
+ phase_inc_error = (2*TWO_TO_THE_31 / step_rate) - (float)d_phase_inc;
+ d_step_error = (gr_int32)(phase_inc_error * (float)STEP_CHECK);
+
+ d_phase_origin += d_step_error/2;
+ d_phase = d_phase_origin;
+
+ }
+
+
// angle_rate is in radians / step
void set_freq (float angle_rate){
- d_phase_inc = gr_fxpt::float_to_fixed (angle_rate);
+ set_period (1/(angle_rate/(2 * M_PI)));
}
// angle_rate is a delta in radians / step
void adjust_freq (float delta_angle_rate)
{
d_phase_inc += gr_fxpt::float_to_fixed (delta_angle_rate);
+ // FIXME: calc phase_inc error
}
// increment current phase angle
@@ -61,11 +93,33 @@
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;
+ }
+ }
}
void step (int n)
{
d_phase += d_phase_inc * n;
+
+ d_step += n;
+ if (d_step >= d_period) {
+ d_phase_origin -= d_period_frac;
+ d_phase = d_phase_origin + d_step_error * ((d_step %
d_period)/STEP_CHECK);
+ d_step %= d_period;
+ } else {
+ if (n >= (int)STEP_CHECK || (d_step % STEP_CHECK) == 0) {
+ d_phase += d_step_error;
+ }
+ }
}
// units are radians / step
_______________________________________________
Discuss-gnuradio mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/discuss-gnuradio
