Introducing fixed point math functions

libs/spandsp/src/Makefile.am

@@ -22,6 +22,7 @@ AM_LDFLAGS = $(COMP_VENDOR_LDFLAGS)
 MAINTAINERCLEANFILES = Makefile.in
 
 DISTCLEANFILES = $(srcdir)/at_interpreter_dictionary.h \
+                 $(srcdir)/math_fixed_tables.h \
                  $(srcdir)/v17_v32bis_rx_fixed_rrc.h \
                  $(srcdir)/v17_v32bis_rx_floating_rrc.h \
                  $(srcdir)/v17_v32bis_tx_fixed_rrc.h \
@@ -55,6 +56,7 @@ EXTRA_DIST = floating_fudge.h \
              libtiff.2008.vcproj \
              filter_tools.c \
              make_at_dictionary.c \
+             make_math_fixed_tables.c \
              make_modem_filter.c \
              msvc/config.h \
              msvc/Download_TIFF.2005.vcproj \
@@ -123,6 +125,7 @@ libspandsp_la_SOURCES = adsi.c \
                         lpc10_encode.c \
                         lpc10_placev.c \
                         lpc10_voicing.c \
+                        math_fixed.c \
                         modem_echo.c \
                         modem_connect_tones.c \
                         noise.c \
@@ -204,6 +207,7 @@ nobase_include_HEADERS = spandsp/adsi.h \
                          spandsp/image_translate.h \
                          spandsp/logging.h \
                          spandsp/lpc10.h \
+                         spandsp/math_fixed.h \
                          spandsp/modem_echo.h \
                          spandsp/modem_connect_tones.h \
                          spandsp/noise.h \
@@ -329,6 +333,9 @@ noinst_HEADERS =         faxfont.h \
 make_at_dictionary$(EXEEXT): $(top_srcdir)/src/make_at_dictionary.c
 	$(CC_FOR_BUILD) -o make_at_dictionary$(EXEEXT) $(top_srcdir)/src/make_at_dictionary.c  -DHAVE_CONFIG_H -I$(top_builddir)/src
 
+make_math_fixed_tables$(EXEEXT): $(top_srcdir)/src/make_math_fixed_tables.c
+	$(CC_FOR_BUILD) -o make_math_fixed_tables$(EXEEXT) $(top_srcdir)/src/make_math_fixed_tables.c  -DHAVE_CONFIG_H -I$(top_builddir)/src -lm
+
 make_modem_filter$(EXEEXT): $(top_srcdir)/src/make_modem_filter.c $(top_srcdir)/src/filter_tools.c
 	$(CC_FOR_BUILD) -o make_modem_filter$(EXEEXT) $(top_srcdir)/src/make_modem_filter.c $(top_srcdir)/src/filter_tools.c -DHAVE_CONFIG_H -I$(top_builddir)/src -lm
 
@@ -342,6 +349,13 @@ at_interpreter.lo: at_interpreter_dictionary.h
 at_interpreter_dictionary.h: make_at_dictionary$(EXEEXT)
 	./make_at_dictionary$(EXEEXT) >at_interpreter_dictionary.h
 
+math_fixed.$(OBJEXT): math_fixed_tables.h
+
+math_fixed.lo: math_fixed_tables.h
+
+math_fixed_tables.h: make_math_fixed_tables$(EXEEXT)
+	./make_math_fixed_tables$(EXEEXT) >math_fixed_tables.h
+
 t4_rx.$(OBJEXT): spandsp/version.h
 
 t4_rx.lo: spandsp/version.h

libs/spandsp/src/make_math_fixed_tables.c

+/*
+ * SpanDSP - a series of DSP components for telephony
+ *
+ * make_fixed_point_math_tables.c - Generate lookup tables for some of the
+ *                                  fixed point maths functions.
+ *
+ * Written by Steve Underwood <steveu@coppice.org>
+ *
+ * Copyright (C) 2010 Steve Underwood
+ *
+ * All rights reserved.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License version 2, as
+ * published by the Free Software Foundation.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+#if defined(HAVE_CONFIG_H)
+#include "config.h"
+#endif
+
+#include <inttypes.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <fcntl.h>
+#include <memory.h>
+#include <math.h>
+
+int main(int argc, char *argv[])
+{
+    int i;
+    double val;
+    int ival;
+
+    printf("static const uint16_t fixed_reciprocal_table[129] =\n");
+    printf("{\n");
+    for (i = 0;  i < 129;  i++)
+    {
+        val = 32768.0*128.0/(128 + i) + 0.5;
+        ival = val;
+        if (i < 128)
+            printf("    %6d,\n", ival);
+        else
+            printf("    %6d\n", ival);
+    }
+    printf("};\n\n");
+
+    printf("static const uint16_t fixed_sqrt_table[193] =\n");
+    printf("{\n");
+    for (i = 64;  i <= 256;  i++)
+    {
+        ival = sqrt(i/256.0)*65536.0 + 0.5;
+        if (ival > 65535)
+            ival = 65535;
+        if (i < 256)
+            printf("    %6d,\n", ival);
+        else
+            printf("    %6d\n", ival);
+    }
+    printf("};\n\n");
+
+    printf("static const int16_t fixed_log10_table[129] =\n");
+    printf("{\n");
+    for (i = 128;  i <= 256;  i++)
+    {
+        ival = log10(i/256.0)*32768.0 - 0.5;
+        if (i <= 255)
+            printf("    %6d,\n", ival);
+        else
+            printf("    %6d\n", ival);
+    }
+    printf("};\n\n");
+
+    printf("static const int16_t fixed_sine_table[257] =\n");
+    printf("{\n");
+    for (i = 0;  i <= 256;  i++)
+    {
+        val = sin(i*3.1415926535/512.0)*32768.0;
+        ival = val + 0.5;
+        if (ival > 32767)
+            ival = 32767;
+        if (i <= 255)
+            printf("    %6d,\n", ival);
+        else
+            printf("    %6d\n", ival);
+    }
+    printf("};\n\n");
+
+    printf("static const uint16_t fixed_arctan_table[257] =\n");
+    printf("{\n");
+    for (i = 0;  i <= 256;  i++)
+    {
+        val = atan(i/256.0)*65536.0/(2.0*3.1415926535);
+        ival = val + 0.5;
+        /* Nudge the result away from zero, so things sit consistently on
+           the correct side of the axes. */
+        if (ival == 0)
+            ival = 1;
+        if (i <= 255)
+            printf("    %6d,\n", ival);
+        else
+            printf("    %6d\n", ival);
+    }
+    printf("};\n\n");
+
+    return 0;
+}

libs/spandsp/src/math_fixed.c

+/*
+ * SpanDSP - a series of DSP components for telephony
+ *
+ * math_fixed.c
+ *
+ * Written by Steve Underwood <steveu@coppice.org>
+ *
+ * Copyright (C) 2010 Steve Underwood
+ *
+ * All rights reserved.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License version 2.1,
+ * as published by the Free Software Foundation.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this program; if not, write to the Free Software
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+/*! \file */
+
+#if defined(HAVE_CONFIG_H)
+#include "config.h"
+#endif
+
+#include <inttypes.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <fcntl.h>
+#include <string.h>
+#include <float.h>
+#if defined(HAVE_TGMATH_H)
+#include <tgmath.h>
+#endif
+#if defined(HAVE_MATH_H)
+#include <math.h>
+#endif
+#include "floating_fudge.h"
+#include <assert.h>
+
+#include "math_fixed_tables.h"
+
+#include "spandsp/telephony.h"
+#include "spandsp/bit_operations.h"
+#include "spandsp/math_fixed.h"
+
+#if defined(SPANDSP_USE_FIXED_POINT)
+SPAN_DECLARE(uint16_t) sqrtu32_u16(uint32_t x)
+{
+    uint16_t zz;
+    uint16_t z;
+    uint16_t i;
+
+    z = 0;
+    for (i = 0x8000;  i;  i >>= 1)
+    {
+        zz = z | i;
+        if (((int32_t) zz*zz) <= x)
+            z = zz;
+    }
+    return z;
+}
+/*- End of function --------------------------------------------------------*/
+#endif
+
+SPAN_DECLARE(uint16_t) fixed_reciprocal16(uint16_t x, int *shift)
+{
+    if (x == 0)
+    {
+        *shift = 0;
+        return 0xFFFF;
+    }
+    *shift = 15 - top_bit(x);
+    x <<= *shift;
+    return fixed_reciprocal_table[((x + 0x80) >> 8) - 128];
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(uint16_t) fixed_divide16(uint16_t y, uint16_t x)
+{
+    int shift;
+    uint32_t z;
+    uint16_t recip;
+
+    if (x == 0)
+        return 0xFFFF;
+    recip = fixed_reciprocal16(x, &shift);
+    z = (((uint32_t) y*recip) >> 15) << shift;
+    return z;
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(uint16_t) fixed_divide32(uint32_t y, uint16_t x)
+{
+    int shift;
+    uint32_t z;
+    uint16_t recip;
+
+    if (x == 0)
+        return 0xFFFF;
+    recip = fixed_reciprocal16(x, &shift);
+    z = (((uint32_t) y*recip) >> 15) << shift;
+    return z;
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(int16_t) fixed_log10_16(uint16_t x)
+{
+    int shift;
+
+    if (x == 0)
+        return 0;
+    shift = 14 - top_bit(x);
+    x <<= shift;
+    return (fixed_log10_table[((x + 0x40) >> 7) - 128] >> 3) - shift*1233;
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(int32_t) fixed_log10_32(uint32_t x)
+{
+    int shift;
+
+    if (x == 0)
+        return 0;
+    shift = 30 - top_bit(x);
+    x <<= shift;
+    return (fixed_log10_table[((x + 0x400000) >> 23) - 128] >> 3) - shift*1233;
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(uint16_t) fixed_sqrt16(uint16_t x)
+{
+    int shift;
+
+    if (x == 0)
+        return 0;
+    shift = 14 - (top_bit(x) & ~1);
+    x <<= shift;
+    //return fixed_sqrt_table[(((x + 0x80) >> 8) & 0xFF) - 64] >> (shift >> 1);
+    return fixed_sqrt_table[((x >> 8) & 0xFF) - 64] >> (shift >> 1);
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(uint16_t) fixed_sqrt32(uint32_t x)
+{
+    int shift;
+
+    if (x == 0)
+        return 0;
+    shift = 30 - (top_bit(x) & ~1);
+    x <<= shift;
+    //return fixed_sqrt_table[(((x + 0x800000) >> 24) & 0xFF) - 64] >> (shift >> 1);
+    return fixed_sqrt_table[((x >> 24) & 0xFF) - 64] >> (shift >> 1);
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(int16_t) fixed_sin(uint16_t x)
+{
+    int step;
+    int step_after;
+    int16_t frac;
+    int16_t z;
+
+    step = (x & 0x3FFF) >> 6;
+    frac = x & 0x3F;
+    if ((x & 0x4000))
+    {
+        step = 256 - step;
+        step_after = step - 1;
+    }
+    else
+    {
+        step_after = step + 1;
+    }
+    z = fixed_sine_table[step] + ((frac*(fixed_sine_table[step_after] - fixed_sine_table[step])) >> 6);
+    if ((x & 0x8000))
+    	z = -z;
+    return z;
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(int16_t) fixed_cos(uint16_t x)
+{
+    int step;
+    int step_after;
+    int16_t frac;
+    int16_t z;
+
+    x += 0x4000;
+    step = (x & 0x3FFF) >> 6;
+    frac = x & 0x3F;
+    if ((x & 0x4000))
+    {
+        step = 256 - step;
+        step_after = step - 1;
+    }
+    else
+    {
+        step_after = step + 1;
+    }
+    z = fixed_sine_table[step] + ((frac*(fixed_sine_table[step_after] - fixed_sine_table[step])) >> 6);
+    if ((x & 0x8000))
+    	z = -z;
+    return z;
+}
+/*- End of function --------------------------------------------------------*/
+
+SPAN_DECLARE(uint16_t) fixed_atan2(int16_t y, int16_t x)
+{
+    int16_t abs_x;
+    int16_t abs_y;
+    uint16_t angle;
+    uint16_t recip;
+    uint32_t z;
+    int step;
+    int shift;
+
+    if (y == 0)
+        return (x & 0x8000);
+    if (x == 0)
+        return ((y & 0x8000) | 0x4000);
+    abs_x = abs(x);
+    abs_y = abs(y);
+    
+    if (abs_y < abs_x)
+    {
+        recip = fixed_reciprocal16(abs_x, &shift);
+        z = (((uint32_t) recip*abs_y) >> 15) << shift;
+        step = z >> 7;
+        angle = fixed_arctan_table[step];
+    }
+    else
+    {
+        recip = fixed_reciprocal16(abs_y, &shift);
+        z = (((uint32_t) recip*abs_x) >> 15) << shift;
+        step = z >> 7;
+        angle = 0x4000 - fixed_arctan_table[step];
+    }
+    /* If we are in quadrant II or III, flip things around */
+    if (x < 0)
+        angle = 0x8000 - angle;
+    /* If we are in quadrant III or IV, negate to return an
+       answer in the full circle range. */
+    if (y < 0)
+        angle = -angle;
+    return (uint16_t) angle;
+}
+/*- End of function --------------------------------------------------------*/
+/*- End of file ------------------------------------------------------------*/

libs/spandsp/src/spandsp/math_fixed.h

+/*
+ * SpanDSP - a series of DSP components for telephony
+ *
+ * math_fixed.h
+ *
+ * Written by Steve Underwood <steveu@coppice.org>
+ *
+ * Copyright (C) 2010 Steve Underwood
+ *
+ * All rights reserved.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License version 2.1,
+ * as published by the Free Software Foundation.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this program; if not, write to the Free Software
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+#if !defined(_MATH_FIXED_H_)
+#define _MATH_FIXED_H_
+
+/*! \page math_fixed_page Fixed point math functions
+
+\section math_fixed_page_sec_1 What does it do?
+
+\section math_fixed_page_sec_2 How does it work?
+*/
+
+#if defined(__cplusplus)
+extern "C"
+{
+#endif
+
+#if defined(SPANDSP_USE_FIXED_POINT)
+SPAN_DECLARE(uint16_t) sqrtu32_u16(uint32_t x);
+#endif
+
+SPAN_DECLARE(uint16_t) fixed_reciprocal16(uint16_t x, int *shift);
+
+SPAN_DECLARE(uint16_t) fixed_divide16(uint16_t y, uint16_t x);
+
+SPAN_DECLARE(uint16_t) fixed_divide32(uint32_t y, uint16_t x);
+
+SPAN_DECLARE(int16_t) fixed_log10_16(uint16_t x);
+
+SPAN_DECLARE(int32_t) fixed_log10_32(uint32_t x);
+
+SPAN_DECLARE(uint16_t) fixed_sqrt16(uint16_t x);
+
+SPAN_DECLARE(uint16_t) fixed_sqrt32(uint32_t x);
+
+/*! Evaluate an approximate 16 bit fixed point sine.
+    \brief Evaluate an approximate 16 bit fixed point sine.
+    \param x A 16 bit unsigned angle, in 360/65536 degree steps.
+    \return sin(x)*32767. */
+SPAN_DECLARE(int16_t) fixed_sin(uint16_t x);
+
+/*! Evaluate an approximate 16 bit fixed point cosine.
+    \brief Evaluate an approximate 16 bit fixed point cosine.
+    \param x A 16 bit unsigned angle, in 360/65536 degree steps.
+    \return cos(x)*32767. */
+SPAN_DECLARE(int16_t) fixed_cos(uint16_t x);
+
+/*! Evaluate an approximate 16 bit fixed point sine.
+    \brief Evaluate an approximate 16 bit fixed point sine.
+    \param y .
+    \param x .
+    \return The 16 bit unsigned angle, in 360/65536 degree steps. */
+SPAN_DECLARE(uint16_t) fixed_atan2(int16_t y, int16_t x);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
+/*- End of file ------------------------------------------------------------*/

libs/spandsp/tests/Makefile.am

@@ -84,6 +84,7 @@ noinst_PROGRAMS =   adsi_tests \
                     line_model_tests \
                     logging_tests \
                     lpc10_tests \
+                    math_fixed_tests \
                     make_g168_css \
                     modem_connect_tones_tests \
                     modem_echo_tests \
@@ -236,6 +237,9 @@ logging_tests_LDADD = $(LIBDIR) -lspandsp
 lpc10_tests_SOURCES = lpc10_tests.c
 lpc10_tests_LDADD = -L$(top_builddir)/spandsp-sim -lspandsp-sim $(LIBDIR) -lspandsp
 
+math_fixed_tests_SOURCES = math_fixed_tests.c
+math_fixed_tests_LDADD = $(LIBDIR) -lspandsp
+
 make_g168_css_SOURCES = make_g168_css.c
 make_g168_css_LDADD = $(LIBDIR) -lspandsp
 
@@ -300,7 +304,7 @@ t38_core_tests_SOURCES = t38_core_tests.c
 t38_core_tests_LDADD = $(LIBDIR) -lspandsp
 
 t38_decode_SOURCES = t38_decode.c fax_utils.c pcap_parse.c udptl.c
-t38_decode_LDADD = $(LIBDIR) -lspandsp -lpcap
+t38_decode_LDADD = -L$(top_builddir)/spandsp-sim -lspandsp-sim $(LIBDIR) -lspandsp -lpcap
 
 t38_gateway_tests_SOURCES = t38_gateway_tests.c fax_utils.c media_monitor.cpp
 t38_gateway_tests_LDADD = -L$(top_builddir)/spandsp-sim -lspandsp-sim $(LIBDIR) -lspandsp