From 3cee0589e5d4738c59ee5f86befc711c7e72807b Mon Sep 17 00:00:00 2001 From: Steve Underwood Date: Sat, 2 Jul 2011 22:04:29 +0800 Subject: [PATCH] Introducing fixed point math functions --- libs/spandsp/src/Makefile.am | 14 + libs/spandsp/src/make_math_fixed_tables.c | 116 ++++++ libs/spandsp/src/math_fixed.c | 255 +++++++++++++ libs/spandsp/src/spandsp/math_fixed.h | 83 +++++ libs/spandsp/tests/Makefile.am | 6 +- libs/spandsp/tests/math_fixed_tests.c | 431 ++++++++++++++++++++++ 6 files changed, 904 insertions(+), 1 deletion(-) create mode 100644 libs/spandsp/src/make_math_fixed_tables.c create mode 100644 libs/spandsp/src/math_fixed.c create mode 100644 libs/spandsp/src/spandsp/math_fixed.h create mode 100644 libs/spandsp/tests/math_fixed_tests.c diff --git a/libs/spandsp/src/Makefile.am b/libs/spandsp/src/Makefile.am index 7e8af50f5f..5c1f216586 100644 --- a/libs/spandsp/src/Makefile.am +++ b/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 diff --git a/libs/spandsp/src/make_math_fixed_tables.c b/libs/spandsp/src/make_math_fixed_tables.c new file mode 100644 index 0000000000..2bddd12ced --- /dev/null +++ b/libs/spandsp/src/make_math_fixed_tables.c @@ -0,0 +1,116 @@ +/* + * 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 + * + * 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 +#include +#include +#include +#include +#include + +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; +} diff --git a/libs/spandsp/src/math_fixed.c b/libs/spandsp/src/math_fixed.c new file mode 100644 index 0000000000..b627fe1719 --- /dev/null +++ b/libs/spandsp/src/math_fixed.c @@ -0,0 +1,255 @@ +/* + * SpanDSP - a series of DSP components for telephony + * + * math_fixed.c + * + * Written by Steve Underwood + * + * 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 +#include +#include +#include +#include +#include +#if defined(HAVE_TGMATH_H) +#include +#endif +#if defined(HAVE_MATH_H) +#include +#endif +#include "floating_fudge.h" +#include + +#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 ------------------------------------------------------------*/ diff --git a/libs/spandsp/src/spandsp/math_fixed.h b/libs/spandsp/src/spandsp/math_fixed.h new file mode 100644 index 0000000000..5bc94382ed --- /dev/null +++ b/libs/spandsp/src/spandsp/math_fixed.h @@ -0,0 +1,83 @@ +/* + * SpanDSP - a series of DSP components for telephony + * + * math_fixed.h + * + * Written by Steve Underwood + * + * 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 ------------------------------------------------------------*/ diff --git a/libs/spandsp/tests/Makefile.am b/libs/spandsp/tests/Makefile.am index 150977f39c..429ed34e98 100644 --- a/libs/spandsp/tests/Makefile.am +++ b/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 diff --git a/libs/spandsp/tests/math_fixed_tests.c b/libs/spandsp/tests/math_fixed_tests.c new file mode 100644 index 0000000000..5111fb4527 --- /dev/null +++ b/libs/spandsp/tests/math_fixed_tests.c @@ -0,0 +1,431 @@ +/* + * SpanDSP - a series of DSP components for telephony + * + * math_fixed_tests.c - Test the fixed point math functions. + * + * Written by Steve Underwood + * + * 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. + */ + +/*! \file */ + +/*! \page math_fixed_tests_page Fixed point math function tests +\section math_fixed_tests_page_sec_1 What does it do? +*/ + +#if defined(HAVE_CONFIG_H) +#include "config.h" +#endif + +#include +#include +#include +#include +#include +#include + +//#if defined(WITH_SPANDSP_INTERNALS) +#define SPANDSP_EXPOSE_INTERNAL_STRUCTURES +//#endif + +#include "spandsp.h" + +static void fixed_reciprocal16_tests(void) +{ + int x; + uint16_t yu16; + uint32_t yu32; + double z; + double ratio; + double max; + double min; + int shift; + + /* The reciprocal should be with about 0.4%, except for 0, which is obviously a + special case. */ + printf("fixed_reciprocal16() function tests\n"); + if (fixed_reciprocal16(0, &shift) != 0xFFFF || shift != 0) + { + printf("Test failed\n"); + exit(2); + } + min = 999999.0; + max = -999999.0; + for (x = 1; x < 65536; x++) + { + yu16 = fixed_reciprocal16(x, &shift); + yu32 = ((uint32_t) yu16) << shift; + z = 32768.0*32768.0/x; + ratio = z/yu32; + //printf("%6x %15d %f %f\n", x, yu32, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < 0.996 || ratio > 1.004) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_divide16_tests(void) +{ + int x; + int y; + uint16_t yu16; + double z; + double ratio; + double max; + double min; + + printf("fixed_divide16() function tests\n"); + if (fixed_divide16(12345, 0) != 0xFFFF) + { + printf("Test failed\n"); + exit(2); + } + min = 999999.0; + max = -999999.0; + for (y = 32; y < 65536; y += 16) + { + for (x = y*16; x < 65536; x += 16) + { + yu16 = fixed_divide16(y, x); + z = 32768.0*y/x; + ratio = z/yu16; + //printf("%6d %6d %6d %f %f\n", x, y, yu16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < 0.996 || ratio > 1.07) + { + printf("Test failed\n"); + exit(2); + } + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_divide32_tests(void) +{ + uint32_t xu32; + uint16_t yu16; + uint32_t yu32; + double z; + double ratio; + double max; + double min; + + printf("fixed_divide32() function tests\n"); + if (fixed_divide32(12345, 0) != 0xFFFF) + { + printf("Test failed\n"); + exit(2); + } + min = 999999.0; + max = -999999.0; + for (yu32 = 32; yu32 < 65536; yu32 += 16) + { + for (xu32 = yu32*16; xu32 < 65535; xu32 += 16) + { + yu16 = fixed_divide32(yu32, xu32); + z = 32768.0*yu32/xu32; + ratio = z/yu16; + //printf("%6u %6u %6u %f %f\n", xu32, yu32, yu16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < 0.996 || ratio > 1.07) + { + printf("Test failed\n"); + exit(2); + } + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_log10_16_tests(void) +{ + int x; + int16_t yi16; + double z; + double ratio; + double max; + double min; + + printf("Log10 16 bit function tests\n"); + min = 999999.0; + max = -999999.0; + for (x = 1; x < 32500; x++) + { + yi16 = fixed_log10_16(x); + z = 4096.0*log10(x/32768.0); + ratio = z - yi16; + //printf("%6d %15d %f %f\n", x, yi16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < -8.0 || ratio > 8.0) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_log10_32_tests(void) +{ + int x; + int32_t yi32; + double z; + double ratio; + double max; + double min; + + printf("fixed_log10_32() function tests\n"); + min = 999999.0; + max = -999999.0; + for (x = 1; x < 32767*65536; x += 0x4000) + { + yi32 = fixed_log10_32(x); + z = 4096.0*log10(x/(32768.0*65536.0)); + ratio = z - yi32; + //printf("%6d %15d %f %f\n", x, yi32, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < -8.0 || ratio > 8.0) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_sqrt16_tests(void) +{ + int x; + uint16_t yu16; + double z; + double ratio; + double max; + double min; + + printf("fixed_sqrt16() function tests\n"); + min = 999999.0; + max = -999999.0; + for (x = 1; x < 65536; x++) + { + yu16 = fixed_sqrt16(x); + z = sqrt(x)*256.0; + ratio = z/yu16; + //printf("%6d %6d %f %f\n", x, yu16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < 0.999 || ratio > 1.008) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_sqrt32_tests(void) +{ + uint32_t xu32; + uint16_t yu16; + double z; + double ratio; + double max; + double min; + + printf("fixed_sqrt32() function tests\n"); + min = 999999.0; + max = -999999.0; + for (xu32 = 20000; xu32 < 0xFFFF0000; xu32 += 10000) + { + yu16 = fixed_sqrt32(xu32); + z = sqrt(xu32); + ratio = z/yu16; + //printf("%10u %6d %f %f\n", xu32, yu16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < 0.999 || ratio > 1.009) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_sin_tests(void) +{ + int x; + int16_t yi16; + double z; + double ratio; + double max; + double min; + + printf("fixed_sin() function tests\n"); + min = 999999.0; + max = -999999.0; + for (x = 0; x < 65536; x++) + { + yi16 = fixed_sin(x); + z = sin(2.0*3.1415926535*x/65536.0)*32768.0; + ratio = z - yi16; + //printf("%6d %6d %f %f\n", x, yi16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < -2.0 || ratio > 2.0) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_cos_tests(void) +{ + int x; + int16_t yi16; + double z; + double ratio; + double max; + double min; + + printf("fixed_cos() function tests\n"); + min = 999999.0; + max = -999999.0; + for (x = 0; x < 65536; x++) + { + yi16 = fixed_cos(x); + z = cos(2.0*3.1415926535*x/65536.0)*32768.0; + ratio = z - yi16; + //printf("%6d %6d %f %f\n", x, yi16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < -2.0 || ratio > 2.0) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +static void fixed_atan2_tests(void) +{ + int i; + int x; + int y; + uint16_t yu16; + double z; + double ratio; + double max; + double min; + + printf("fixed_atan2() function tests\n"); + min = 999999.0; + max = -999999.0; + for (i = 0; i < 65536; i++) + { + x = 16384.0*cos(i*2.0*3.1415926535/65536.0); + y = 16384.0*sin(i*2.0*3.1415926535/65536.0); + yu16 = fixed_atan2(y, x); + z = atan2(y/32768.0, x/32768.0)*65536.0/(2.0*3.1415926535); + if (z < 0.0) + z += 65536.0; + ratio = z - yu16; + //printf("%6d %6d %6d %6d %f %f\n", i, x, y, yu16, z, ratio); + if (ratio < min) + min = ratio; + if (ratio > max) + max = ratio; + if (ratio < -43.0 || ratio > 43.0) + { + printf("Test failed\n"); + exit(2); + } + } + printf("min %f, max %f\n", min, max); + printf("Test passed\n"); +} +/*- End of function --------------------------------------------------------*/ + +int main(int argc, char *argv[]) +{ + fixed_reciprocal16_tests(); + fixed_divide16_tests(); + fixed_divide32_tests(); + fixed_log10_16_tests(); + fixed_log10_32_tests(); + fixed_sqrt16_tests(); + fixed_sqrt32_tests(); + fixed_sin_tests(); + fixed_cos_tests(); + fixed_atan2_tests(); + + printf("Tests passed\n"); + + return 0; +} +/*- End of function --------------------------------------------------------*/ +/*- End of file ------------------------------------------------------------*/