Mercurial > hg > may

changeset 448:3a7b70ecfb70

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Threshold for Kaiser terms
author Chris Cannam
date Thu, 24 Oct 2013 10:01:36 +0100
parents 926f2986a5e0
children ea95ee1cb3a6
files src/may/signal/window.yeti
diffstat 1 files changed, 10 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- 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;