diff options
| author | Marcus Müller <marcus.mueller@ettus.com> | 2017-06-27 21:19:04 +0200 |
|---|---|---|
| committer | Marcus Müller <marcus.mueller@ettus.com> | 2017-06-27 21:19:04 +0200 |
| commit | 071d27a8b9ce11594fc387be2f3fee58b6813d4d (patch) | |
| tree | 465d9c2ce2e66de8efbf4c436e5a486f92305fa6 /gr-fft | |
| parent | b57a37f7c676542f08a27d3f141f4d9ed2ab1132 (diff) | |
Fix: hazardous floating point sqrt(1.0 - (1.0+-epsilon)²)
https://github.com/gnuradio/gnuradio/issues/1348
sqrt(x<0) yields NaN; to avoid a situation where a variable would be
close to, but not necessarily exactly +-1, extracted the relevant
floating point corner cases from the loop.
Diffstat (limited to 'gr-fft')
| -rw-r--r-- | gr-fft/lib/window.cc | 11 |
1 files changed, 10 insertions, 1 deletions
diff --git a/gr-fft/lib/window.cc b/gr-fft/lib/window.cc index 126b28978d..cfbdb93dbd 100644 --- a/gr-fft/lib/window.cc +++ b/gr-fft/lib/window.cc @@ -261,10 +261,19 @@ namespace gr { double inm1 = 1.0/((double)(ntaps-1)); double temp; - for(int i = 0; i < ntaps; i++) { + /* extracting first and last element out of the loop, since + sqrt(1.0-temp*temp) might trigger unexpected floating point behaviour + if |temp| = 1.0+epsilon, which can happen for i==0 and + 1/i==1/(ntaps-1)==inm1 ; compare + https://github.com/gnuradio/gnuradio/issues/1348 . + In any case, the 0. Bessel function of first kind is 1 at point 0. + */ + taps[0] = IBeta; + for(int i = 1; i < ntaps-1; i++) { temp = 2*i*inm1 - 1; taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta; } + taps[ntaps-1] = IBeta; return taps; } |