Filter window calculation illegal calculation.

[cor-mt]

This bug get triggered on an "AMD FX(tm)-8150 Eight-Core Processor" when gnuradio is compiled with CFLAGS="-march=native -O2".

In window.cc line 265, with i == ntaps - 1 gives temp = 1.0000000000000000002 that cause sqrt (0) witch return "-nan" on the next line for the last value of the vector. This value get used for the gain calculation in the low_pass filter and thus screwing all the rest of the calculations.

[marcusmueller]

Thanks for documenting/reporting this here 👍

For future readers: http://lists.gnu.org/archive/html/discuss-gnuradio/2017-06/msg00277.html

[marcusmueller]

Just to be clear: this is a serious bug in the way we do floating point math, and I (hopefully) have a patch coming up. We need to be vigilant, though – this is most likely not the only mistake of this type.

[marcusmueller]

diff --git a/gr-fft/lib/window.cc b/gr-fft/lib/window.cc
index 126b28978..e80d469e9 100644
--- a/gr-fft/lib/window.cc
+++ b/gr-fft/lib/window.cc
@@ -30,9 +30,10 @@
 namespace gr {
   namespace fft {
 
-#define IzeroEPSILON 1E-21               /* Max error acceptable in Izero */
+#define BESSEL_EPSILON 1E-21               /* Max error acceptable in bessel_i_zero */
 
-    static double Izero(double x)
+    /*! \brief Iterative first kind Bessel function approximation */
+    static double bessel_i_zero(double x)
     {
       double sum, u, halfx, temp;
       int n;
@@ -45,7 +46,7 @@ namespace gr {
         temp *= temp;
         u *= temp;
         sum += u;
-      } while (u >= IzeroEPSILON*sum);
+      } while (u >= BESSEL_EPSILON*sum);
       return(sum);
     }
 
@@ -257,13 +258,17 @@ namespace gr {
 
       std::vector<float> taps(ntaps);
 
-      double IBeta = 1.0/Izero(beta);
+      double bessel_i_beta = 1.0/bessel_i_zero(beta);
       double inm1 = 1.0/((double)(ntaps-1));
       double temp;
 
-      for(int i = 0; i < ntaps; i++) {
+      /* dragging taps[0] out of the loop, because ((double)(-1))*((double)(-1))
+      might be 1.0 + epsilon and then we get sqrt() of something < 0, i.e. NaN.
+      https://github.com/gnuradio/gnuradio/issues/1348 */
+      taps[0] = bessel_i_zero(0) * bessel_i_beta;
+      for(int i = 1; i < ntaps; i++) {
         temp = 2*i*inm1 - 1;
-        taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta;
+        taps[i] = bessel_i_zero(beta*sqrt(1.0-temp*temp)) * bessel_i_beta;
       }
       return taps;
     }

[noc0lour]

Is this fixed now (: ?

[mbr0wn]

This was fixed in #1349. Closing.