r372 - forth/lib

Author: wmb
Date: 2007-05-12 03:15:35 +0200 (Sat, 12 May 2007)
New Revision: 372

Added:
   forth/lib/isin.fth
Log:
New file - Integer sin() calculation.

Added: forth/lib/isin.fth
===================================================================
--- forth/lib/isin.fth	                        (rev 0)
+++ forth/lib/isin.fth	2007-05-12 01:15:35 UTC (rev 372)
@@ -0,0 +1,92 @@
+purpose: Compute a fixed-point sin table
+\ See license at end of file
+
+\ Calculate the sin function using fractional integer arithmetic
+\ where the scale factor is 2^15.  1.0 is represented by 32767,
+\ pi is 102943.  The argument to isin ( index -- frac ) is a
+\ sample index from 0 to fs/freq/4 .
+
+d# 16000 value fs
+0 value x
+0 value xsq
+
+d#  32767 constant one
+d# 102943 constant pi
+
+0 value freq
+0 value fstep
+0 value #cycle
+0 value #half-cycle
+0 value #quarter-cycle
+
+: set-freq  ( freq -- )
+   dup to freq
+   pi * to fstep
+   fs freq /  dup  to #cycle
+   2/         dup  to #half-cycle
+   2/              to #quarter-cycle
+;
+: set-period  ( quarter-cycle -- )
+   dup to #quarter-cycle       ( quarter-cycle )
+   2* dup to #half-cycle       ( half-cycle )
+   2* dup to #cycle            ( cycle )
+   fs over / to freq           ( period )
+   pi swap /  fs *  to fstep   ( )
+;
+
+\ Multiply two fractional numbers where the scale factor is 2^15
+\ : times  ( n1 n2 -- n3 )  d# 32767 */  ;   \ Insignificantly more accurate, but slower
+: times  ( n1 n2 -- n3 )  *  d# 15 rshift   ;
+
+\ Computes  (1 - (x^2 / divisor) * last)
+: sin-step  ( last divisor -- next )  xsq  swap /  times  one min  one swap -  ;
+
+\ Taylor series expansion of sin, calculated as
+\ x * (1 - (x^2/(2*3)) * (1 - (x^2/(4*5)) * (1 - (x^2/(6*7)) * (1 - (x^2/(8*9))))))
+\ This is good for the first quadrant only, i.e. 0 <= index <= fs / freq / 4
+: isin  ( index -- frac )
+   fstep *  fs 2/  /  to x
+   x dup times to xsq
+   one  d# 72 sin-step d# 42 sin-step  d# 20 sin-step  6 sin-step  x times  one max
+;
+
+: one-cycle  ( adr -- )
+   #quarter-cycle 1+  0  do   ( adr )
+      i isin
+      2dup  swap  i wa+ w!                ( adr isin )
+      2dup  swap  #half-cycle i - wa+ w!  ( adr isin )
+      negate                              ( adr -isin )
+      2dup  swap  #half-cycle i + wa+ w!  ( adr -isin )
+      over  #cycle i - wa+ w!             ( adr )
+   loop                                   ( adr )
+   drop
+;
+
+\ d# 16000 to fs
+\ d# 150 set-freq
+\ here one-cycle
+\ here #cycle 2*  wdump
+
+\ LICENSE_BEGIN
+\ Copyright (c) 2006 FirmWorks
+\ 
+\ Permission is hereby granted, free of charge, to any person obtaining
+\ a copy of this software and associated documentation files (the
+\ "Software"), to deal in the Software without restriction, including
+\ without limitation the rights to use, copy, modify, merge, publish,
+\ distribute, sublicense, and/or sell copies of the Software, and to
+\ permit persons to whom the Software is furnished to do so, subject to
+\ the following conditions:
+\ 
+\ The above copyright notice and this permission notice shall be
+\ included in all copies or substantial portions of the Software.
+\ 
+\ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+\ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+\ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+\ NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+\ LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+\ OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+\ WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+\
+\ LICENSE_END