# HG changeset patch
# User Chris Cannam
# Date 1382605296 -3600
# Node ID 3a7b70ecfb706ebadce4764c858a756394f29f66
# Parent  926f2986a5e00dae6cedbbc91991094956194cbc
Threshold for Kaiser terms

diff -r 926f2986a5e0 -r 3a7b70ecfb70 src/may/signal/window.yeti
--- a/src/may/signal/window.yeti	Thu Oct 24 09:18:15 2013 +0100
+++ b/src/may/signal/window.yeti	Thu Oct 24 10:01:36 2013 +0100
@@ -103,12 +103,16 @@
 calculateKaiser { length, beta } =
    (terms = 20;
     factorials = array (map maths.factorial [1..terms]);
-    bes0 x = 
-        1 + sum
-           (map do i:
-                f = factorials[i-1]; // this is i!, not (i-1)!
-                (maths.pow (x/2) (i*2)) / (f * f);
-                done [1..terms]);
+    threshold = 1e-12;
+    bes0term x i =
+       (f = factorials[i-1]; // this is i!, not (i-1)!
+        (maths.pow (x/2) (i*2)) / (f * f));
+    bes0sub x i acc =
+       (term = bes0term x i;
+        if i == terms or abs term < threshold then acc + term
+        else bes0sub x (i+1) (acc + term)
+        fi);
+    bes0 x = bes0sub x 1 1;
     denominator = bes0 beta;
     even = (length % 2 == 0);
     halflen = if even then length/2 else (length+1)/2 fi;