diff options
Diffstat (limited to 'gsl-1.9/wavelet')
| -rw-r--r-- | gsl-1.9/wavelet/ChangeLog | 14 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/Makefile.am | 17 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/Makefile.in | 546 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/TODO | 11 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/bspline.c | 591 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/daubechies.c | 458 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/dwt.c | 390 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/gsl_wavelet.h | 108 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/gsl_wavelet2d.h | 107 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/haar.c | 81 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/test.c | 270 | ||||
| -rw-r--r-- | gsl-1.9/wavelet/wavelet.c | 132 |
12 files changed, 2725 insertions, 0 deletions
diff --git a/gsl-1.9/wavelet/ChangeLog b/gsl-1.9/wavelet/ChangeLog new file mode 100644 index 0000000..9d5f240 --- /dev/null +++ b/gsl-1.9/wavelet/ChangeLog @@ -0,0 +1,14 @@ +2006-03-16 Brian Gough <bjg@network-theory.co.uk> + + * changed to gsl_wavelet_forward and gsl_wavelet_backward enums + throughout internally instead of forward and backward. + +2004-12-29 Brian Gough <bjg@network-theory.co.uk> + + * gsl_wavelet.h: added missing includes, use GSL_VAR instead of + extern + +2004-07-23 Brian Gough <bjg@network-theory.co.uk> + + * added wavelet directory from Ivo Alxneit. + diff --git a/gsl-1.9/wavelet/Makefile.am b/gsl-1.9/wavelet/Makefile.am new file mode 100644 index 0000000..fba86b4 --- /dev/null +++ b/gsl-1.9/wavelet/Makefile.am @@ -0,0 +1,17 @@ +noinst_LTLIBRARIES = libgslwavelet.la + +pkginclude_HEADERS = gsl_wavelet.h gsl_wavelet2d.h + +INCLUDES= -I$(top_builddir) + +libgslwavelet_la_SOURCES = dwt.c wavelet.c bspline.c daubechies.c haar.c + +check_PROGRAMS = test + +TESTS = $(check_PROGRAMS) + +test_LDADD = libgslwavelet.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la + +test_SOURCES = test.c + + diff --git a/gsl-1.9/wavelet/Makefile.in b/gsl-1.9/wavelet/Makefile.in new file mode 100644 index 0000000..ea12ea8 --- /dev/null +++ b/gsl-1.9/wavelet/Makefile.in @@ -0,0 +1,546 @@ +# Makefile.in generated by automake 1.9.6 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, +# 2003, 2004, 2005 Free Software Foundation, Inc. +# This Makefile.in is free software; the Free Software Foundation +# gives unlimited permission to copy and/or distribute it, +# with or without modifications, as long as this notice is preserved. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, to the extent permitted by law; without +# even the implied warranty of MERCHANTABILITY or FITNESS FOR A +# PARTICULAR PURPOSE. + +@SET_MAKE@ + + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ +pkgdatadir = $(datadir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +top_builddir = .. +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +INSTALL = @INSTALL@ +install_sh_DATA = $(install_sh) -c -m 644 +install_sh_PROGRAM = $(install_sh) -c +install_sh_SCRIPT = $(install_sh) -c +INSTALL_HEADER = $(INSTALL_DATA) +transform = $(program_transform_name) +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +build_triplet = @build@ +host_triplet = @host@ +check_PROGRAMS = test$(EXEEXT) +subdir = wavelet +DIST_COMMON = $(pkginclude_HEADERS) $(srcdir)/Makefile.am \ + $(srcdir)/Makefile.in ChangeLog TODO +ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 +am__aclocal_m4_deps = $(top_srcdir)/configure.ac +am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \ + $(ACLOCAL_M4) +mkinstalldirs = $(SHELL) $(top_srcdir)/mkinstalldirs +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +LTLIBRARIES = $(noinst_LTLIBRARIES) +libgslwavelet_la_LIBADD = +am_libgslwavelet_la_OBJECTS = dwt.lo wavelet.lo bspline.lo \ + daubechies.lo haar.lo +libgslwavelet_la_OBJECTS = $(am_libgslwavelet_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslwavelet.la ../blas/libgslblas.la \ + ../cblas/libgslcblas.la ../matrix/libgslmatrix.la \ + ../vector/libgslvector.la ../block/libgslblock.la \ + ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la \ + ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir) +depcomp = +am__depfiles_maybe = +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) \ + $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \ + $(AM_CFLAGS) $(CFLAGS) +CCLD = $(CC) +LINK = $(LIBTOOL) --tag=CC --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \ + $(AM_LDFLAGS) $(LDFLAGS) -o $@ +SOURCES = $(libgslwavelet_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslwavelet_la_SOURCES) $(test_SOURCES) +am__vpath_adj_setup = srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; +am__vpath_adj = case $$p in \ + $(srcdir)/*) f=`echo "$$p" | sed "s|^$$srcdirstrip/||"`;; \ + *) f=$$p;; \ + esac; +am__strip_dir = `echo $$p | sed -e 's|^.*/||'`; +am__installdirs = "$(DESTDIR)$(pkgincludedir)" +pkgincludeHEADERS_INSTALL = $(INSTALL_HEADER) +HEADERS = $(pkginclude_HEADERS) +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AR = @AR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFS = @DEFS@ +ECHO = @ECHO@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +GSL_CFLAGS = @GSL_CFLAGS@ +GSL_LIBS = @GSL_LIBS@ +GSL_LT_CBLAS_VERSION = @GSL_LT_CBLAS_VERSION@ +GSL_LT_VERSION = @GSL_LT_VERSION@ +HAVE_AIX_IEEE_INTERFACE = @HAVE_AIX_IEEE_INTERFACE@ +HAVE_DARWIN86_IEEE_INTERFACE = @HAVE_DARWIN86_IEEE_INTERFACE@ +HAVE_DARWIN_IEEE_INTERFACE = @HAVE_DARWIN_IEEE_INTERFACE@ +HAVE_EXTENDED_PRECISION_REGISTERS = @HAVE_EXTENDED_PRECISION_REGISTERS@ +HAVE_FREEBSD_IEEE_INTERFACE = @HAVE_FREEBSD_IEEE_INTERFACE@ +HAVE_GNUM68K_IEEE_INTERFACE = @HAVE_GNUM68K_IEEE_INTERFACE@ +HAVE_GNUPPC_IEEE_INTERFACE = @HAVE_GNUPPC_IEEE_INTERFACE@ +HAVE_GNUSPARC_IEEE_INTERFACE = @HAVE_GNUSPARC_IEEE_INTERFACE@ +HAVE_GNUX86_IEEE_INTERFACE = @HAVE_GNUX86_IEEE_INTERFACE@ +HAVE_HPUX11_IEEE_INTERFACE = @HAVE_HPUX11_IEEE_INTERFACE@ +HAVE_HPUX_IEEE_INTERFACE = @HAVE_HPUX_IEEE_INTERFACE@ +HAVE_IEEE_COMPARISONS = @HAVE_IEEE_COMPARISONS@ +HAVE_IEEE_DENORMALS = @HAVE_IEEE_DENORMALS@ +HAVE_INLINE = @HAVE_INLINE@ +HAVE_IRIX_IEEE_INTERFACE = @HAVE_IRIX_IEEE_INTERFACE@ +HAVE_NETBSD_IEEE_INTERFACE = @HAVE_NETBSD_IEEE_INTERFACE@ +HAVE_OPENBSD_IEEE_INTERFACE = @HAVE_OPENBSD_IEEE_INTERFACE@ +HAVE_OS2EMX_IEEE_INTERFACE = @HAVE_OS2EMX_IEEE_INTERFACE@ +HAVE_PRINTF_LONGDOUBLE = @HAVE_PRINTF_LONGDOUBLE@ +HAVE_SOLARIS_IEEE_INTERFACE = @HAVE_SOLARIS_IEEE_INTERFACE@ +HAVE_SUNOS4_IEEE_INTERFACE = @HAVE_SUNOS4_IEEE_INTERFACE@ +HAVE_TRU64_IEEE_INTERFACE = @HAVE_TRU64_IEEE_INTERFACE@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LDFLAGS = @LDFLAGS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +MAINT = @MAINT@ +MAINTAINER_MODE_FALSE = @MAINTAINER_MODE_FALSE@ +MAINTAINER_MODE_TRUE = @MAINTAINER_MODE_TRUE@ +MAKEINFO = @MAKEINFO@ +OBJEXT = @OBJEXT@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +RANLIB = @RANLIB@ +RELEASED = @RELEASED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STRIP = @STRIP@ +VERSION = @VERSION@ +ac_ct_AR = @ac_ct_AR@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_RANLIB = @ac_ct_RANLIB@ +ac_ct_STRIP = @ac_ct_STRIP@ +am__leading_dot = @am__leading_dot@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build = @build@ +build_alias = @build_alias@ +build_cpu = @build_cpu@ +build_os = @build_os@ +build_vendor = @build_vendor@ +datadir = @datadir@ +exec_prefix = @exec_prefix@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +noinst_LTLIBRARIES = libgslwavelet.la +pkginclude_HEADERS = gsl_wavelet.h gsl_wavelet2d.h +INCLUDES = -I$(top_builddir) +libgslwavelet_la_SOURCES = dwt.c wavelet.c bspline.c daubechies.c haar.c +TESTS = $(check_PROGRAMS) +test_LDADD = libgslwavelet.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +test_SOURCES = test.c +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .lo .o .obj +$(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \ + && exit 0; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu --ignore-deps wavelet/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps wavelet/Makefile +.PRECIOUS: Makefile +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + @case '$?' in \ + *config.status*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \ + *) \ + echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: @MAINTAINER_MODE_TRUE@ $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): @MAINTAINER_MODE_TRUE@ $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +clean-noinstLTLIBRARIES: + -test -z "$(noinst_LTLIBRARIES)" || rm -f $(noinst_LTLIBRARIES) + @list='$(noinst_LTLIBRARIES)'; for p in $$list; do \ + dir="`echo $$p | sed -e 's|/[^/]*$$||'`"; \ + test "$$dir" != "$$p" || dir=.; \ + echo "rm -f \"$${dir}/so_locations\""; \ + rm -f "$${dir}/so_locations"; \ + done +libgslwavelet.la: $(libgslwavelet_la_OBJECTS) $(libgslwavelet_la_DEPENDENCIES) + $(LINK) $(libgslwavelet_la_LDFLAGS) $(libgslwavelet_la_OBJECTS) $(libgslwavelet_la_LIBADD) $(LIBS) + +clean-checkPROGRAMS: + @list='$(check_PROGRAMS)'; for p in $$list; do \ + f=`echo $$p|sed 's/$(EXEEXT)$$//'`; \ + echo " rm -f $$p $$f"; \ + rm -f $$p $$f ; \ + done +test$(EXEEXT): $(test_OBJECTS) $(test_DEPENDENCIES) + @rm -f test$(EXEEXT) + $(LINK) $(test_LDFLAGS) $(test_OBJECTS) $(test_LDADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +.c.o: + $(COMPILE) -c $< + +.c.obj: + $(COMPILE) -c `$(CYGPATH_W) '$<'` + +.c.lo: + $(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + +distclean-libtool: + -rm -f libtool +uninstall-info-am: +install-pkgincludeHEADERS: $(pkginclude_HEADERS) + @$(NORMAL_INSTALL) + test -z "$(pkgincludedir)" || $(mkdir_p) "$(DESTDIR)$(pkgincludedir)" + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + if test -f "$$p"; then d=; else d="$(srcdir)/"; fi; \ + f=$(am__strip_dir) \ + echo " $(pkgincludeHEADERS_INSTALL) '$$d$$p' '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + $(pkgincludeHEADERS_INSTALL) "$$d$$p" "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +uninstall-pkgincludeHEADERS: + @$(NORMAL_UNINSTALL) + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + f=$(am__strip_dir) \ + echo " rm -f '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + rm -f "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$tags $$unique; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + test -z "$(CTAGS_ARGS)$$tags$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$tags $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && cd $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) $$here + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +check-TESTS: $(TESTS) + @failed=0; all=0; xfail=0; xpass=0; skip=0; \ + srcdir=$(srcdir); export srcdir; \ + list='$(TESTS)'; \ + if test -n "$$list"; then \ + for tst in $$list; do \ + if test -f ./$$tst; then dir=./; \ + elif test -f $$tst; then dir=; \ + else dir="$(srcdir)/"; fi; \ + if $(TESTS_ENVIRONMENT) $${dir}$$tst; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xpass=`expr $$xpass + 1`; \ + failed=`expr $$failed + 1`; \ + echo "XPASS: $$tst"; \ + ;; \ + *) \ + echo "PASS: $$tst"; \ + ;; \ + esac; \ + elif test $$? -ne 77; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xfail=`expr $$xfail + 1`; \ + echo "XFAIL: $$tst"; \ + ;; \ + *) \ + failed=`expr $$failed + 1`; \ + echo "FAIL: $$tst"; \ + ;; \ + esac; \ + else \ + skip=`expr $$skip + 1`; \ + echo "SKIP: $$tst"; \ + fi; \ + done; \ + if test "$$failed" -eq 0; then \ + if test "$$xfail" -eq 0; then \ + banner="All $$all tests passed"; \ + else \ + banner="All $$all tests behaved as expected ($$xfail expected failures)"; \ + fi; \ + else \ + if test "$$xpass" -eq 0; then \ + banner="$$failed of $$all tests failed"; \ + else \ + banner="$$failed of $$all tests did not behave as expected ($$xpass unexpected passes)"; \ + fi; \ + fi; \ + dashes="$$banner"; \ + skipped=""; \ + if test "$$skip" -ne 0; then \ + skipped="($$skip tests were not run)"; \ + test `echo "$$skipped" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$skipped"; \ + fi; \ + report=""; \ + if test "$$failed" -ne 0 && test -n "$(PACKAGE_BUGREPORT)"; then \ + report="Please report to $(PACKAGE_BUGREPORT)"; \ + test `echo "$$report" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$report"; \ + fi; \ + dashes=`echo "$$dashes" | sed s/./=/g`; \ + echo "$$dashes"; \ + echo "$$banner"; \ + test -z "$$skipped" || echo "$$skipped"; \ + test -z "$$report" || echo "$$report"; \ + echo "$$dashes"; \ + test "$$failed" -eq 0; \ + else :; fi + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \ + list='$(DISTFILES)'; for file in $$list; do \ + case $$file in \ + $(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \ + $(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \ + esac; \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test "$$dir" != "$$file" && test "$$dir" != "."; then \ + dir="/$$dir"; \ + $(mkdir_p) "$(distdir)$$dir"; \ + else \ + dir=''; \ + fi; \ + if test -d $$d/$$file; then \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \ + fi; \ + cp -pR $$d/$$file $(distdir)$$dir || exit 1; \ + else \ + test -f $(distdir)/$$file \ + || cp -p $$d/$$file $(distdir)/$$file \ + || exit 1; \ + fi; \ + done +check-am: all-am + $(MAKE) $(AM_MAKEFLAGS) $(check_PROGRAMS) + $(MAKE) $(AM_MAKEFLAGS) check-TESTS +check: check-am +all-am: Makefile $(LTLIBRARIES) $(HEADERS) +installdirs: + for dir in "$(DESTDIR)$(pkgincludedir)"; do \ + test -z "$$dir" || $(mkdir_p) "$$dir"; \ + done +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + `test -z '$(STRIP)' || \ + echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install +mostlyclean-generic: + +clean-generic: + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES mostlyclean-am + +distclean: distclean-am + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-libtool distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +info: info-am + +info-am: + +install-data-am: install-pkgincludeHEADERS + +install-exec-am: + +install-info: install-info-am + +install-man: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic \ + mostlyclean-libtool + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: uninstall-info-am uninstall-pkgincludeHEADERS + +.PHONY: CTAGS GTAGS all all-am check check-TESTS check-am clean \ + clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES ctags distclean distclean-compile \ + distclean-generic distclean-libtool distclean-tags distdir dvi \ + dvi-am html html-am info info-am install install-am \ + install-data install-data-am install-exec install-exec-am \ + install-info install-info-am install-man \ + install-pkgincludeHEADERS install-strip installcheck \ + installcheck-am installdirs maintainer-clean \ + maintainer-clean-generic mostlyclean mostlyclean-compile \ + mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \ + tags uninstall uninstall-am uninstall-info-am \ + uninstall-pkgincludeHEADERS + +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/gsl-1.9/wavelet/TODO b/gsl-1.9/wavelet/TODO new file mode 100644 index 0000000..ea5347b --- /dev/null +++ b/gsl-1.9/wavelet/TODO @@ -0,0 +1,11 @@ +* Implement more wavelet types: + Check the literature to find out what are commonly used types. Candidates + could be coiflets and symmlets. +** Coefficients for coiflets and symmlets found so far are only with a + precision of about eight digts. This is probaly insufficient. + +* Wavelet packet transform: + Should include utility functions for selecting the coefficients according + to "dwt-type", "best basis" or "best level". + +* Continuous wavelet transform. diff --git a/gsl-1.9/wavelet/bspline.c b/gsl-1.9/wavelet/bspline.c new file mode 100644 index 0000000..7ac1411 --- /dev/null +++ b/gsl-1.9/wavelet/bspline.c @@ -0,0 +1,591 @@ +/* wavelet/bspline.c + * + * Copyright (C) 2004 Ivo Alxneit + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Coefficients are from A. Cohen, I. Daubechies, and J.-C. Feauveau; + * "Biorthogonal Bases of Compactly Supported Wavelets", Communications + * on Pure and Applied Mathematics, 45 (1992) 485--560 (table 6.1). + * + * Note the following errors in table 1: + * + * N = 2, N~ = 4, m0~ + * the second term in z^-1 (45/64 z^-1) should be left out. + * + * N = 3, N~ = 7, m0~ + * the term 336z^-3 should read 363z^-3. + */ + +#include <config.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_wavelet.h> + +static const double h1_103[6] = { -0.0883883476483184405501055452631, + 0.0883883476483184405501055452631, + M_SQRT1_2, + M_SQRT1_2, + 0.0883883476483184405501055452631, + -0.0883883476483184405501055452631 +}; + +static const double g2_103[6] = { -0.0883883476483184405501055452631, + -0.0883883476483184405501055452631, + M_SQRT1_2, + -(M_SQRT1_2), + 0.0883883476483184405501055452631, + 0.0883883476483184405501055452631 +}; + +static const double h1_105[10] = { 0.0165728151840597076031447897368, + -0.0165728151840597076031447897368, + -0.1215339780164378557563951247368, + 0.1215339780164378557563951247368, + M_SQRT1_2, + M_SQRT1_2, + 0.1215339780164378557563951247368, + -0.1215339780164378557563951247368, + -0.0165728151840597076031447897368, + 0.0165728151840597076031447897368 +}; + +static const double g2_105[10] = { 0.0165728151840597076031447897368, + 0.0165728151840597076031447897368, + -0.1215339780164378557563951247368, + -0.1215339780164378557563951247368, + M_SQRT1_2, + -(M_SQRT1_2), + 0.1215339780164378557563951247368, + 0.1215339780164378557563951247368, + -0.0165728151840597076031447897368, + -0.0165728151840597076031447897368 +}; + +static const double g1_1[10] = { 0.0, 0.0, 0.0, 0.0, + M_SQRT1_2, + -(M_SQRT1_2), + 0.0, 0.0, 0.0, 0.0 +}; + +static const double h2_1[10] = { 0.0, 0.0, 0.0, 0.0, + M_SQRT1_2, + M_SQRT1_2, + 0.0, 0.0, 0.0, 0.0 +}; + +static const double h1_202[6] = { -0.1767766952966368811002110905262, + 0.3535533905932737622004221810524, + 1.0606601717798212866012665431573, + 0.3535533905932737622004221810524, + -0.1767766952966368811002110905262, + 0.0 +}; + +static const double g2_202[6] = { 0.0, + -0.1767766952966368811002110905262, + -0.3535533905932737622004221810524, + 1.0606601717798212866012665431573, + -0.3535533905932737622004221810524, + -0.1767766952966368811002110905262 +}; + +static const double h1_204[10] = { 0.0331456303681194152062895794737, + -0.0662912607362388304125791589473, + -0.1767766952966368811002110905262, + 0.4198446513295125926130013399998, + 0.9943689110435824561886873842099, + 0.4198446513295125926130013399998, + -0.1767766952966368811002110905262, + -0.0662912607362388304125791589473, + 0.0331456303681194152062895794737, + 0.0 +}; + +static const double g2_204[10] = { 0.0, + 0.0331456303681194152062895794737, + 0.0662912607362388304125791589473, + -0.1767766952966368811002110905262, + -0.4198446513295125926130013399998, + 0.9943689110435824561886873842099, + -0.4198446513295125926130013399998, + -0.1767766952966368811002110905262, + 0.0662912607362388304125791589473, + 0.0331456303681194152062895794737 +}; + +static const double h1_206[14] = { -0.0069053396600248781679769957237, + 0.0138106793200497563359539914474, + 0.0469563096881691715422435709210, + -0.1077232986963880994204411332894, + -0.1698713556366120029322340948025, + 0.4474660099696121052849093228945, + 0.9667475524034829435167794013152, + 0.4474660099696121052849093228945, + -0.1698713556366120029322340948025, + -0.1077232986963880994204411332894, + 0.0469563096881691715422435709210, + 0.0138106793200497563359539914474, + -0.0069053396600248781679769957237, + 0.0 +}; + +static const double g2_206[14] = { 0.0, + -0.0069053396600248781679769957237, + -0.0138106793200497563359539914474, + 0.0469563096881691715422435709210, + 0.1077232986963880994204411332894, + -0.1698713556366120029322340948025, + -0.4474660099696121052849093228945, + 0.9667475524034829435167794013152, + -0.4474660099696121052849093228945, + -0.1698713556366120029322340948025, + 0.1077232986963880994204411332894, + 0.0469563096881691715422435709210, + -0.0138106793200497563359539914474, + -0.0069053396600248781679769957237, +}; + +static const double h1_208[18] = { 0.0015105430506304420992449678146, + -0.0030210861012608841984899356291, + -0.0129475118625466465649568669819, + 0.0289161098263541773284036695929, + 0.0529984818906909399392234421792, + -0.1349130736077360572068505539514, + -0.1638291834340902345352542235443, + 0.4625714404759165262773590010400, + 0.9516421218971785225243297231697, + 0.4625714404759165262773590010400, + -0.1638291834340902345352542235443, + -0.1349130736077360572068505539514, + 0.0529984818906909399392234421792, + 0.0289161098263541773284036695929, + -0.0129475118625466465649568669819, + -0.0030210861012608841984899356291, + 0.0015105430506304420992449678146, + 0.0 +}; + +static const double g2_208[18] = { 0.0, + 0.0015105430506304420992449678146, + 0.0030210861012608841984899356291, + -0.0129475118625466465649568669819, + -0.0289161098263541773284036695929, + 0.0529984818906909399392234421792, + 0.1349130736077360572068505539514, + -0.1638291834340902345352542235443, + -0.4625714404759165262773590010400, + 0.9516421218971785225243297231697, + -0.4625714404759165262773590010400, + -0.1638291834340902345352542235443, + 0.1349130736077360572068505539514, + 0.0529984818906909399392234421792, + -0.0289161098263541773284036695929, + -0.0129475118625466465649568669819, + 0.0030210861012608841984899356291, + 0.0015105430506304420992449678146, +}; + +static const double h2_2[18] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.3535533905932737622004221810524, + 0.7071067811865475244008443621048, + 0.3535533905932737622004221810524, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +}; + +static const double g1_2[18] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + -0.3535533905932737622004221810524, + 0.7071067811865475244008443621048, + -0.3535533905932737622004221810524, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +}; + +static const double h1_301[4] = { -0.3535533905932737622004221810524, + 1.0606601717798212866012665431573, + 1.0606601717798212866012665431573, + -0.3535533905932737622004221810524 +}; + +static const double g2_301[4] = { 0.3535533905932737622004221810524, + 1.0606601717798212866012665431573, + -1.0606601717798212866012665431573, + -0.3535533905932737622004221810524 +}; + +static const double h1_303[8] = { 0.0662912607362388304125791589473, + -0.1988737822087164912377374768420, + -0.1546796083845572709626847042104, + 0.9943689110435824561886873842099, + 0.9943689110435824561886873842099, + -0.1546796083845572709626847042104, + -0.1988737822087164912377374768420, + 0.0662912607362388304125791589473 +}; + +static const double g2_303[8] = { -0.0662912607362388304125791589473, + -0.1988737822087164912377374768420, + 0.1546796083845572709626847042104, + 0.9943689110435824561886873842099, + -0.9943689110435824561886873842099, + -0.1546796083845572709626847042104, + 0.1988737822087164912377374768420, + 0.0662912607362388304125791589473 +}; + +static const double h1_305[12] = { -0.0138106793200497563359539914474, + 0.0414320379601492690078619743421, + 0.0524805814161890740766251675000, + -0.2679271788089652729175074340788, + -0.0718155324642587329469607555263, + 0.9667475524034829435167794013152, + 0.9667475524034829435167794013152, + -0.0718155324642587329469607555263, + -0.2679271788089652729175074340788, + 0.0524805814161890740766251675000, + 0.0414320379601492690078619743421, + -0.0138106793200497563359539914474 +}; + +static const double g2_305[12] = { 0.0138106793200497563359539914474, + 0.0414320379601492690078619743421, + -0.0524805814161890740766251675000, + -0.2679271788089652729175074340788, + 0.0718155324642587329469607555263, + 0.9667475524034829435167794013152, + -0.9667475524034829435167794013152, + -0.0718155324642587329469607555263, + 0.2679271788089652729175074340788, + 0.0524805814161890740766251675000, + -0.0414320379601492690078619743421, + -0.0138106793200497563359539914474 +}; + +static const double h1_307[16] = { 0.0030210861012608841984899356291, + -0.0090632583037826525954698068873, + -0.0168317654213106405344439270765, + 0.0746639850740189951912512662623, + 0.0313329787073628846871956180962, + -0.3011591259228349991008967259990, + -0.0264992409453454699696117210896, + 0.9516421218971785225243297231697, + 0.9516421218971785225243297231697, + -0.0264992409453454699696117210896, + -0.3011591259228349991008967259990, + 0.0313329787073628846871956180962, + 0.0746639850740189951912512662623, + -0.0168317654213106405344439270765, + -0.0090632583037826525954698068873, + 0.0030210861012608841984899356291 +}; + +static const double g2_307[16] = { -0.0030210861012608841984899356291, + -0.0090632583037826525954698068873, + 0.0168317654213106405344439270765, + 0.0746639850740189951912512662623, + -0.0313329787073628846871956180962, + -0.3011591259228349991008967259990, + 0.0264992409453454699696117210896, + 0.9516421218971785225243297231697, + -0.9516421218971785225243297231697, + -0.0264992409453454699696117210896, + 0.3011591259228349991008967259990, + 0.0313329787073628846871956180962, + -0.0746639850740189951912512662623, + -0.0168317654213106405344439270765, + 0.0090632583037826525954698068873, + 0.0030210861012608841984899356291 +}; + +static const double h1_309[20] = { -0.0006797443727836989446602355165, + 0.0020392331183510968339807065496, + 0.0050603192196119810324706421788, + -0.0206189126411055346546938106687, + -0.0141127879301758447558029850103, + 0.0991347824942321571990197448581, + 0.0123001362694193142367090236328, + -0.3201919683607785695513833204624, + 0.0020500227115698857061181706055, + 0.9421257006782067372990864259380, + 0.9421257006782067372990864259380, + 0.0020500227115698857061181706055, + -0.3201919683607785695513833204624, + 0.0123001362694193142367090236328, + 0.0991347824942321571990197448581, + -0.0141127879301758447558029850103, + -0.0206189126411055346546938106687, + 0.0050603192196119810324706421788, + 0.0020392331183510968339807065496, + -0.0006797443727836989446602355165 +}; + +static const double g2_309[20] = { 0.0006797443727836989446602355165, + 0.0020392331183510968339807065496, + -0.0050603192196119810324706421788, + -0.0206189126411055346546938106687, + 0.0141127879301758447558029850103, + 0.0991347824942321571990197448581, + -0.0123001362694193142367090236328, + -0.3201919683607785695513833204624, + -0.0020500227115698857061181706055, + 0.9421257006782067372990864259380, + -0.9421257006782067372990864259380, + 0.0020500227115698857061181706055, + 0.3201919683607785695513833204624, + 0.0123001362694193142367090236328, + -0.0991347824942321571990197448581, + -0.0141127879301758447558029850103, + 0.0206189126411055346546938106687, + 0.0050603192196119810324706421788, + -0.0020392331183510968339807065496, + -0.0006797443727836989446602355165 +}; + +static const double h2_3[20] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.1767766952966368811002110905262, + 0.5303300858899106433006332715786, + 0.5303300858899106433006332715786, + 0.1767766952966368811002110905262, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +}; + +static const double g1_3[20] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + -0.1767766952966368811002110905262, + 0.5303300858899106433006332715786, + -0.5303300858899106433006332715786, + 0.1767766952966368811002110905262, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +}; + +static int +bspline_init (const double **h1, const double **g1, + const double **h2, const double **g2, size_t * nc, + size_t * offset, size_t member) +{ + switch (member) + { + case 103: + *nc = 6; + *h1 = h1_103; + *g1 = &g1_1[2]; + *h2 = &h2_1[2]; + *g2 = g2_103; + break; + + case 105: + *nc = 10; + *h1 = h1_105; + *g1 = g1_1; + *h2 = h2_1; + *g2 = g2_105; + break; + + case 202: + *nc = 6; + *h1 = h1_202; + *g1 = &g1_2[6]; + *h2 = &h2_2[6]; + *g2 = g2_202; + break; + + case 204: + *nc = 10; + *h1 = h1_204; + *g1 = &g1_2[4]; + *h2 = &h2_2[4]; + *g2 = g2_204; + break; + + case 206: + *nc = 14; + *h1 = h1_206; + *g1 = &g1_2[2]; + *h2 = &h2_2[2]; + *g2 = g2_206; + break; + + case 208: + *nc = 18; + *h1 = h1_208; + *g1 = g1_2; + *h2 = h2_2; + *g2 = g2_208; + break; + + case 301: + *nc = 4; + *h1 = h1_301; + *g1 = &g1_3[8]; + *h2 = &h2_3[8]; + *g2 = g2_301; + break; + + case 303: + *nc = 8; + *h1 = h1_303; + *g1 = &g1_3[6]; + *h2 = &h2_3[6]; + *g2 = g2_303; + break; + + case 305: + *nc = 12; + *h1 = h1_305; + *g1 = &g1_3[4]; + *h2 = &h2_3[4]; + *g2 = g2_305; + break; + + case 307: + *nc = 16; + *h1 = h1_307; + *g1 = &g1_3[2]; + *h2 = &h2_3[2]; + *g2 = g2_307; + break; + + case 309: + *nc = 20; + *h1 = h1_309; + *g1 = g1_3; + *h2 = h2_3; + *g2 = g2_309; + break; + + default: + return GSL_FAILURE; + } + + *offset = 0; + + return GSL_SUCCESS; +} + +static int +bspline_centered_init (const double **h1, const double **g1, + const double **h2, const double **g2, size_t * nc, + size_t * offset, size_t member) +{ + switch (member) + { + case 103: + *nc = 6; + *h1 = h1_103; + *g1 = &g1_1[2]; + *h2 = &h2_1[2]; + *g2 = g2_103; + break; + + case 105: + *nc = 10; + *h1 = h1_105; + *g1 = g1_1; + *h2 = h2_1; + *g2 = g2_105; + break; + + case 202: + *nc = 6; + *h1 = h1_202; + *g1 = &g1_2[6]; + *h2 = &h2_2[6]; + *g2 = g2_202; + break; + + case 204: + *nc = 10; + *h1 = h1_204; + *g1 = &g1_2[4]; + *h2 = &h2_2[4]; + *g2 = g2_204; + break; + + case 206: + *nc = 14; + *h1 = h1_206; + *g1 = &g1_2[2]; + *h2 = &h2_2[2]; + *g2 = g2_206; + break; + + case 208: + *nc = 18; + *h1 = h1_208; + *g1 = g1_2; + *h2 = h2_2; + *g2 = g2_208; + break; + + case 301: + *nc = 4; + *h1 = h1_301; + *g1 = &g1_3[8]; + *h2 = &h2_3[8]; + *g2 = g2_301; + break; + + case 303: + *nc = 8; + *h1 = h1_303; + *g1 = &g1_3[6]; + *h2 = &h2_3[6]; + *g2 = g2_303; + break; + + case 305: + *nc = 12; + *h1 = h1_305; + *g1 = &g1_3[4]; + *h2 = &h2_3[4]; + *g2 = g2_305; + break; + + case 307: + *nc = 16; + *h1 = h1_307; + *g1 = &g1_3[2]; + *h2 = &h2_3[2]; + *g2 = g2_307; + break; + + case 309: + *nc = 20; + *h1 = h1_309; + *g1 = g1_3; + *h2 = h2_3; + *g2 = g2_309; + break; + + default: + return GSL_FAILURE; + } + + *offset = ((*nc) >> 1); + + return GSL_SUCCESS; +} + +static const gsl_wavelet_type bspline_type = { + "bspline", + &bspline_init +}; + +static const gsl_wavelet_type bspline_centered_type = { + "bspline-centered", + &bspline_centered_init +}; + +const gsl_wavelet_type *gsl_wavelet_bspline = &bspline_type; +const gsl_wavelet_type *gsl_wavelet_bspline_centered = &bspline_centered_type; diff --git a/gsl-1.9/wavelet/daubechies.c b/gsl-1.9/wavelet/daubechies.c new file mode 100644 index 0000000..67b7daa --- /dev/null +++ b/gsl-1.9/wavelet/daubechies.c @@ -0,0 +1,458 @@ +/* wavelet/daubechies.c + * + * Copyright (C) 2004 Ivo Alxneit + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* + * Coefficients for Daubechies wavelets of extremal phase are from + * I. Daubechies, "Orthonormal Bases of Compactly Supported Wavelets", + * Communications on Pure and Applied Mathematics, 41 (1988) 909--996 + * (table 1). + * Additional digits have been obtained using the Mathematica package + * Daubechies.m by Tong Chen & Meng Xu available at + * http://www.cwp.mines.edu/wavelets/. + */ + +#include <config.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_wavelet.h> + +static const double h_4[4] = { 0.48296291314453414337487159986, + 0.83651630373780790557529378092, + 0.22414386804201338102597276224, + -0.12940952255126038117444941881 +}; + +static const double g_4[4] = { -0.12940952255126038117444941881, + -0.22414386804201338102597276224, + 0.83651630373780790557529378092, + -0.48296291314453414337487159986 +}; + +static const double h_6[6] = { 0.33267055295008261599851158914, + 0.80689150931109257649449360409, + 0.45987750211849157009515194215, + -0.13501102001025458869638990670, + -0.08544127388202666169281916918, + 0.03522629188570953660274066472 +}; + +static const double g_6[6] = { 0.03522629188570953660274066472, + 0.08544127388202666169281916918, + -0.13501102001025458869638990670, + -0.45987750211849157009515194215, + 0.80689150931109257649449360409, + -0.33267055295008261599851158914 +}; + +static const double h_8[8] = { 0.23037781330889650086329118304, + 0.71484657055291564708992195527, + 0.63088076792985890788171633830, + -0.02798376941685985421141374718, + -0.18703481171909308407957067279, + 0.03084138183556076362721936253, + 0.03288301166688519973540751355, + -0.01059740178506903210488320852 +}; + +static const double g_8[8] = { -0.01059740178506903210488320852, + -0.03288301166688519973540751355, + 0.03084138183556076362721936253, + 0.18703481171909308407957067279, + -0.02798376941685985421141374718, + -0.63088076792985890788171633830, + 0.71484657055291564708992195527, + -0.23037781330889650086329118304 +}; + +static const double h_10[10] = { 0.16010239797419291448072374802, + 0.60382926979718967054011930653, + 0.72430852843777292772807124410, + 0.13842814590132073150539714634, + -0.24229488706638203186257137947, + -0.03224486958463837464847975506, + 0.07757149384004571352313048939, + -0.00624149021279827427419051911, + -0.01258075199908199946850973993, + 0.00333572528547377127799818342 +}; + +static const double g_10[10] = { 0.00333572528547377127799818342, + 0.01258075199908199946850973993, + -0.00624149021279827427419051911, + -0.07757149384004571352313048939, + -0.03224486958463837464847975506, + 0.24229488706638203186257137947, + 0.13842814590132073150539714634, + -0.72430852843777292772807124410, + 0.60382926979718967054011930653, + -0.16010239797419291448072374802 +}; + +static const double h_12[12] = { 0.11154074335010946362132391724, + 0.49462389039845308567720417688, + 0.75113390802109535067893449844, + 0.31525035170919762908598965481, + -0.22626469396543982007631450066, + -0.12976686756726193556228960588, + 0.09750160558732304910234355254, + 0.02752286553030572862554083950, + -0.03158203931748602956507908070, + 0.00055384220116149613925191840, + 0.00477725751094551063963597525, + -0.00107730108530847956485262161 +}; + +static const double g_12[12] = { -0.00107730108530847956485262161, + -0.00477725751094551063963597525, + 0.00055384220116149613925191840, + 0.03158203931748602956507908070, + 0.02752286553030572862554083950, + -0.09750160558732304910234355254, + -0.12976686756726193556228960588, + 0.22626469396543982007631450066, + 0.31525035170919762908598965481, + -0.75113390802109535067893449844, + 0.49462389039845308567720417688, + -0.11154074335010946362132391724 +}; + +static const double h_14[14] = { 0.07785205408500917901996352196, + 0.39653931948191730653900039094, + 0.72913209084623511991694307034, + 0.46978228740519312247159116097, + -0.14390600392856497540506836221, + -0.22403618499387498263814042023, + 0.07130921926683026475087657050, + 0.08061260915108307191292248036, + -0.03802993693501441357959206160, + -0.01657454163066688065410767489, + 0.01255099855609984061298988603, + 0.00042957797292136652113212912, + -0.00180164070404749091526826291, + 0.00035371379997452024844629584 +}; + +static const double g_14[14] = { 0.00035371379997452024844629584, + 0.00180164070404749091526826291, + 0.00042957797292136652113212912, + -0.01255099855609984061298988603, + -0.01657454163066688065410767489, + 0.03802993693501441357959206160, + 0.08061260915108307191292248036, + -0.07130921926683026475087657050, + -0.22403618499387498263814042023, + 0.14390600392856497540506836221, + 0.46978228740519312247159116097, + -0.72913209084623511991694307034, + 0.39653931948191730653900039094, + -0.07785205408500917901996352196 +}; + +static const double h_16[16] = { 0.05441584224310400995500940520, + 0.31287159091429997065916237551, + 0.67563073629728980680780076705, + 0.58535468365420671277126552005, + -0.01582910525634930566738054788, + -0.28401554296154692651620313237, + 0.00047248457391328277036059001, + 0.12874742662047845885702928751, + -0.01736930100180754616961614887, + -0.04408825393079475150676372324, + 0.01398102791739828164872293057, + 0.00874609404740577671638274325, + -0.00487035299345157431042218156, + -0.00039174037337694704629808036, + 0.00067544940645056936636954757, + -0.00011747678412476953373062823 +}; + +static const double g_16[16] = { -0.00011747678412476953373062823, + -0.00067544940645056936636954757, + -0.00039174037337694704629808036, + 0.00487035299345157431042218156, + 0.00874609404740577671638274325, + -0.01398102791739828164872293057, + -0.04408825393079475150676372324, + 0.01736930100180754616961614887, + 0.12874742662047845885702928751, + -0.00047248457391328277036059001, + -0.28401554296154692651620313237, + 0.01582910525634930566738054788, + 0.58535468365420671277126552005, + -0.67563073629728980680780076705, + 0.31287159091429997065916237551, + -0.05441584224310400995500940520 +}; + +static const double h_18[18] = { 0.03807794736387834658869765888, + 0.24383467461259035373204158165, + 0.60482312369011111190307686743, + 0.65728807805130053807821263905, + 0.13319738582500757619095494590, + -0.29327378327917490880640319524, + -0.09684078322297646051350813354, + 0.14854074933810638013507271751, + 0.03072568147933337921231740072, + -0.06763282906132997367564227483, + 0.00025094711483145195758718975, + 0.02236166212367909720537378270, + -0.00472320475775139727792570785, + -0.00428150368246342983449679500, + 0.00184764688305622647661912949, + 0.00023038576352319596720521639, + -0.00025196318894271013697498868, + 0.00003934732031627159948068988 +}; + +static const double g_18[18] = { 0.00003934732031627159948068988, + 0.00025196318894271013697498868, + 0.00023038576352319596720521639, + -0.00184764688305622647661912949, + -0.00428150368246342983449679500, + 0.00472320475775139727792570785, + 0.02236166212367909720537378270, + -0.00025094711483145195758718975, + -0.06763282906132997367564227483, + -0.03072568147933337921231740072, + 0.14854074933810638013507271751, + 0.09684078322297646051350813354, + -0.29327378327917490880640319524, + -0.13319738582500757619095494590, + 0.65728807805130053807821263905, + -0.60482312369011111190307686743, + 0.24383467461259035373204158165, + -0.03807794736387834658869765888 +}; + +static const double h_20[20] = { 0.02667005790055555358661744877, + 0.18817680007769148902089297368, + 0.52720118893172558648174482796, + 0.68845903945360356574187178255, + 0.28117234366057746074872699845, + -0.24984642432731537941610189792, + -0.19594627437737704350429925432, + 0.12736934033579326008267723320, + 0.09305736460357235116035228984, + -0.07139414716639708714533609308, + -0.02945753682187581285828323760, + 0.03321267405934100173976365318, + 0.00360655356695616965542329142, + -0.01073317548333057504431811411, + 0.00139535174705290116578931845, + 0.00199240529518505611715874224, + -0.00068585669495971162656137098, + -0.00011646685512928545095148097, + 0.00009358867032006959133405013, + -0.00001326420289452124481243668 +}; + +static const double g_20[20] = { -0.00001326420289452124481243668, + -0.00009358867032006959133405013, + -0.00011646685512928545095148097, + 0.00068585669495971162656137098, + 0.00199240529518505611715874224, + -0.00139535174705290116578931845, + -0.01073317548333057504431811411, + -0.00360655356695616965542329142, + 0.03321267405934100173976365318, + 0.02945753682187581285828323760, + -0.07139414716639708714533609308, + -0.09305736460357235116035228984, + 0.12736934033579326008267723320, + 0.19594627437737704350429925432, + -0.24984642432731537941610189792, + -0.28117234366057746074872699845, + 0.68845903945360356574187178255, + -0.52720118893172558648174482796, + 0.18817680007769148902089297368, + -0.02667005790055555358661744877 +}; + +static int +daubechies_init (const double **h1, const double **g1, const double **h2, + const double **g2, size_t * nc, size_t * offset, + size_t member) +{ + switch (member) + { + case 4: + *h1 = h_4; + *g1 = g_4; + *h2 = h_4; + *g2 = g_4; + break; + + case 6: + *h1 = h_6; + *g1 = g_6; + *h2 = h_6; + *g2 = g_6; + break; + + case 8: + *h1 = h_8; + *g1 = g_8; + *h2 = h_8; + *g2 = g_8; + break; + + case 10: + *h1 = h_10; + *g1 = g_10; + *h2 = h_10; + *g2 = g_10; + break; + + case 12: + *h1 = h_12; + *g1 = g_12; + *h2 = h_12; + *g2 = g_12; + break; + + case 14: + *h1 = h_14; + *g1 = g_14; + *h2 = h_14; + *g2 = g_14; + break; + + case 16: + *h1 = h_16; + *g1 = g_16; + *h2 = h_16; + *g2 = g_16; + break; + + case 18: + *h1 = h_18; + *g1 = g_18; + *h2 = h_18; + *g2 = g_18; + break; + + case 20: + *h1 = h_20; + *g1 = g_20; + *h2 = h_20; + *g2 = g_20; + break; + + default: + return GSL_FAILURE; + } + + *nc = member; + *offset = 0; + + return GSL_SUCCESS; +} + +static int +daubechies_centered_init (const double **h1, const double **g1, + const double **h2, const double **g2, size_t * nc, + size_t * offset, size_t member) +{ + switch (member) + { + case 4: + *h1 = h_4; + *g1 = g_4; + *h2 = h_4; + *g2 = g_4; + break; + + case 6: + *h1 = h_6; + *g1 = g_6; + *h2 = h_6; + *g2 = g_6; + break; + + case 8: + *h1 = h_8; + *g1 = g_8; + *h2 = h_8; + *g2 = g_8; + break; + + case 10: + *h1 = h_10; + *g1 = g_10; + *h2 = h_10; + *g2 = g_10; + break; + + case 12: + *h1 = h_12; + *g1 = g_12; + *h2 = h_12; + *g2 = g_12; + break; + + case 14: + *h1 = h_14; + *g1 = g_14; + *h2 = h_14; + *g2 = g_14; + break; + + case 16: + *h1 = h_16; + *g1 = g_16; + *h2 = h_16; + *g2 = g_16; + break; + + case 18: + *h1 = h_18; + *g1 = g_18; + *h2 = h_18; + *g2 = g_18; + break; + + case 20: + *h1 = h_20; + *g1 = g_20; + *h2 = h_20; + *g2 = g_20; + break; + + default: + return GSL_FAILURE; + } + + *nc = member; + *offset = (member >> 1); + + return GSL_SUCCESS; +} + +static const gsl_wavelet_type daubechies_type = { + "daubechies", + &daubechies_init +}; + +static const gsl_wavelet_type daubechies_centered_type = { + "daubechies-centered", + &daubechies_centered_init +}; + +const gsl_wavelet_type *gsl_wavelet_daubechies = &daubechies_type; +const gsl_wavelet_type *gsl_wavelet_daubechies_centered = + &daubechies_centered_type; diff --git a/gsl-1.9/wavelet/dwt.c b/gsl-1.9/wavelet/dwt.c new file mode 100644 index 0000000..7d621eb --- /dev/null +++ b/gsl-1.9/wavelet/dwt.c @@ -0,0 +1,390 @@ +/* wavelet/dwt.c + * + * Copyright (C) 2004 Ivo Alxneit + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* function dwt_step is based on the public domain function pwt.c + * available from http://www.numerical-recipes.com + */ + +#include <config.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_wavelet.h> +#include <gsl/gsl_wavelet2d.h> + +#define ELEMENT(a,stride,i) ((a)[(stride)*(i)]) + +static int binary_logn (const size_t n); +static void dwt_step (const gsl_wavelet * w, double *a, size_t stride, size_t n, gsl_wavelet_direction dir, gsl_wavelet_workspace * work); + +static int +binary_logn (const size_t n) +{ + size_t ntest; + size_t logn = 0; + size_t k = 1; + + while (k < n) + { + k *= 2; + logn++; + } + + ntest = (1 << logn); + + if (n != ntest) + { + return -1; /* n is not a power of 2 */ + } + + return logn; +} + +static void +dwt_step (const gsl_wavelet * w, double *a, size_t stride, size_t n, + gsl_wavelet_direction dir, gsl_wavelet_workspace * work) +{ + double ai, ai1; + size_t i, ii; + size_t jf; + size_t k; + size_t n1, ni, nh, nmod; + + for (i = 0; i < work->n; i++) + { + work->scratch[i] = 0.0; + } + + nmod = w->nc * n; + nmod -= w->offset; /* center support */ + + n1 = n - 1; + nh = n >> 1; + + if (dir == gsl_wavelet_forward) + { + for (ii = 0, i = 0; i < n; i += 2, ii++) + { + ni = i + nmod; + for (k = 0; k < w->nc; k++) + { + jf = n1 & (ni + k); + work->scratch[ii] += w->h1[k] * ELEMENT (a, stride, jf); + work->scratch[ii + nh] += w->g1[k] * ELEMENT (a, stride, jf); + } + } + } + else + { + for (ii = 0, i = 0; i < n; i += 2, ii++) + { + ai = ELEMENT (a, stride, ii); + ai1 = ELEMENT (a, stride, ii + nh); + ni = i + nmod; + for (k = 0; k < w->nc; k++) + { + jf = (n1 & (ni + k)); + work->scratch[jf] += (w->h2[k] * ai + w->g2[k] * ai1); + } + } + } + + for (i = 0; i < n; i++) + { + ELEMENT (a, stride, i) = work->scratch[i]; + } +} + +int +gsl_wavelet_transform (const gsl_wavelet * w, + double *data, size_t stride, size_t n, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work) +{ + size_t i; + + if (work->n < n) + { + GSL_ERROR ("not enough workspace provided", GSL_EINVAL); + } + + if (binary_logn (n) == -1) + { + GSL_ERROR ("n is not a power of 2", GSL_EINVAL); + } + + if (n < 2) + { + return GSL_SUCCESS; + } + + if (dir == gsl_wavelet_forward) + { + for (i = n; i >= 2; i >>= 1) + { + dwt_step (w, data, stride, i, dir, work); + } + } + else + { + for (i = 2; i <= n; i <<= 1) + { + dwt_step (w, data, stride, i, dir, work); + } + } + + return GSL_SUCCESS; +} + +int +gsl_wavelet_transform_forward (const gsl_wavelet * w, + double *data, size_t stride, size_t n, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet_transform (w, data, stride, n, gsl_wavelet_forward, work); +} + +int +gsl_wavelet_transform_inverse (const gsl_wavelet * w, + double *data, size_t stride, size_t n, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet_transform (w, data, stride, n, gsl_wavelet_backward, work); +} + + +/* Leaving this out for now BJG */ +#if 0 +int +gsl_dwt_vector (const gsl_wavelet * w, gsl_vector *v, gsl_wavelet_direction + dir, gsl_wavelet_workspace * work) +{ + return gsl_dwt (w, v->data, v->stride, v->size, dir, work); +} +#endif + +int +gsl_wavelet2d_transform (const gsl_wavelet * w, + double *data, size_t tda, size_t size1, + size_t size2, gsl_wavelet_direction dir, + gsl_wavelet_workspace * work) +{ + size_t i; + + if (size1 != size2) + { + GSL_ERROR ("2d dwt works only with square matrix", GSL_EINVAL); + } + + if (work->n < size1) + { + GSL_ERROR ("not enough workspace provided", GSL_EINVAL); + } + + if (binary_logn (size1) == -1) + { + GSL_ERROR ("n is not a power of 2", GSL_EINVAL); + } + + if (size1 < 2) + { + return GSL_SUCCESS; + } + + if (dir == gsl_wavelet_forward) + { + for (i = 0; i < size1; i++) /* for every row j */ + { + gsl_wavelet_transform (w, &ELEMENT(data, tda, i), 1, size1, dir, work); + } + for (i = 0; i < size2; i++) /* for every column j */ + { + gsl_wavelet_transform (w, &ELEMENT(data, 1, i), tda, size2, dir, work); + } + } + else + { + for (i = 0; i < size2; i++) /* for every column j */ + { + gsl_wavelet_transform (w, &ELEMENT(data, 1, i), tda, size2, dir, work); + } + for (i = 0; i < size1; i++) /* for every row j */ + { + gsl_wavelet_transform (w, &ELEMENT(data, tda, i), 1, size1, dir, work); + } + } + + return GSL_SUCCESS; +} + +int +gsl_wavelet2d_nstransform (const gsl_wavelet * w, + double *data, size_t tda, size_t size1, + size_t size2, gsl_wavelet_direction dir, + gsl_wavelet_workspace * work) +{ + size_t i, j; + + if (size1 != size2) + { + GSL_ERROR ("2d dwt works only with square matrix", GSL_EINVAL); + } + + if (work->n < size1) + { + GSL_ERROR ("not enough workspace provided", GSL_EINVAL); + } + + if (binary_logn (size1) == -1) + { + GSL_ERROR ("n is not a power of 2", GSL_EINVAL); + } + + if (size1 < 2) + { + return GSL_SUCCESS; + } + + if (dir == gsl_wavelet_forward) + { + for (i = size1; i >= 2; i >>= 1) + { + for (j = 0; j < i; j++) /* for every row j */ + { + dwt_step (w, &ELEMENT(data, tda, j), 1, i, dir, work); + } + for (j = 0; j < i; j++) /* for every column j */ + { + dwt_step (w, &ELEMENT(data, 1, j), tda, i, dir, work); + } + } + } + else + { + for (i = 2; i <= size1; i <<= 1) + { + for (j = 0; j < i; j++) /* for every column j */ + { + dwt_step (w, &ELEMENT(data, 1, j), tda, i, dir, work); + } + for (j = 0; j < i; j++) /* for every row j */ + { + dwt_step (w, &ELEMENT(data, tda, j), 1, i, dir, work); + } + } + } + + return GSL_SUCCESS; +} + + +int +gsl_wavelet2d_transform_forward (const gsl_wavelet * w, + double *data, size_t tda, size_t size1, + size_t size2, gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_transform (w, data, tda, size1, size2, gsl_wavelet_forward, work); +} + +int +gsl_wavelet2d_transform_inverse (const gsl_wavelet * w, + double *data, size_t tda, size_t size1, + size_t size2, gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_transform (w, data, tda, size1, size2, gsl_wavelet_backward, work); +} + +int +gsl_wavelet2d_nstransform_forward (const gsl_wavelet * w, + double *data, size_t tda, size_t size1, + size_t size2, gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_nstransform (w, data, tda, size1, size2, gsl_wavelet_forward, work); +} + +int +gsl_wavelet2d_nstransform_inverse (const gsl_wavelet * w, + double *data, size_t tda, size_t size1, + size_t size2, gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_nstransform (w, data, tda, size1, size2, gsl_wavelet_backward, work); +} + + +int +gsl_wavelet2d_transform_matrix (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_transform (w, a->data, + a->tda, a->size1, a->size2, + dir, work); +} + +int +gsl_wavelet2d_transform_matrix_forward (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_transform (w, a->data, + a->tda, a->size1, a->size2, + gsl_wavelet_forward, work); +} + +int +gsl_wavelet2d_transform_matrix_inverse (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_transform (w, a->data, + a->tda, a->size1, a->size2, + gsl_wavelet_backward, work); +} + +int +gsl_wavelet2d_nstransform_matrix (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_nstransform (w, a->data, + a->tda, a->size1, a->size2, + dir, work); +} + +int +gsl_wavelet2d_nstransform_matrix_forward (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_nstransform (w, a->data, + a->tda, a->size1, a->size2, + gsl_wavelet_forward, work); +} + +int +gsl_wavelet2d_nstransform_matrix_inverse (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work) +{ + return gsl_wavelet2d_nstransform (w, a->data, + a->tda, a->size1, a->size2, + gsl_wavelet_backward, work); +} + + diff --git a/gsl-1.9/wavelet/gsl_wavelet.h b/gsl-1.9/wavelet/gsl_wavelet.h new file mode 100644 index 0000000..fde012c --- /dev/null +++ b/gsl-1.9/wavelet/gsl_wavelet.h @@ -0,0 +1,108 @@ +/* wavelet/gsl_wavelet.h + * + * Copyright (C) 2004 Ivo Alxneit + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_WAVELET_H__ +#define __GSL_WAVELET_H__ +#include <stdlib.h> +#include <gsl/gsl_types.h> +#include <gsl/gsl_errno.h> + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +#ifndef GSL_DISABLE_DEPRECATED +typedef enum { + forward = 1, backward = -1, + gsl_wavelet_forward = 1, gsl_wavelet_backward = -1 +} +gsl_wavelet_direction; +#else +typedef enum { + gsl_wavelet_forward = 1, gsl_wavelet_backward = -1 +} +gsl_wavelet_direction; +#endif + +typedef struct +{ + const char *name; + int (*init) (const double **h1, const double **g1, + const double **h2, const double **g2, size_t * nc, + size_t * offset, size_t member); +} +gsl_wavelet_type; + +typedef struct +{ + const gsl_wavelet_type *type; + const double *h1; + const double *g1; + const double *h2; + const double *g2; + size_t nc; + size_t offset; +} +gsl_wavelet; + +typedef struct +{ + double *scratch; + size_t n; +} +gsl_wavelet_workspace; + +GSL_VAR const gsl_wavelet_type *gsl_wavelet_daubechies; +GSL_VAR const gsl_wavelet_type *gsl_wavelet_daubechies_centered; +GSL_VAR const gsl_wavelet_type *gsl_wavelet_haar; +GSL_VAR const gsl_wavelet_type *gsl_wavelet_haar_centered; +GSL_VAR const gsl_wavelet_type *gsl_wavelet_bspline; +GSL_VAR const gsl_wavelet_type *gsl_wavelet_bspline_centered; + +gsl_wavelet *gsl_wavelet_alloc (const gsl_wavelet_type * T, size_t k); +void gsl_wavelet_free (gsl_wavelet * w); +const char *gsl_wavelet_name (const gsl_wavelet * w); + +gsl_wavelet_workspace *gsl_wavelet_workspace_alloc (size_t n); +void gsl_wavelet_workspace_free (gsl_wavelet_workspace * work); + +int gsl_wavelet_transform (const gsl_wavelet * w, + double *data, size_t stride, size_t n, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work); + +int gsl_wavelet_transform_forward (const gsl_wavelet * w, + double *data, size_t stride, size_t n, + gsl_wavelet_workspace * work); + +int gsl_wavelet_transform_inverse (const gsl_wavelet * w, + double *data, size_t stride, size_t n, + gsl_wavelet_workspace * work); + +__END_DECLS + +#endif /* __GSL_WAVELET_H__ */ diff --git a/gsl-1.9/wavelet/gsl_wavelet2d.h b/gsl-1.9/wavelet/gsl_wavelet2d.h new file mode 100644 index 0000000..1365755 --- /dev/null +++ b/gsl-1.9/wavelet/gsl_wavelet2d.h @@ -0,0 +1,107 @@ +/* wavelet/gsl_wavelet.h + * + * Copyright (C) 2004 Ivo Alxneit + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_WAVELET2D_H__ +#define __GSL_WAVELET2D_H__ +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_vector_double.h> +#include <gsl/gsl_matrix_double.h> +#include <gsl/gsl_wavelet.h> + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +int gsl_wavelet2d_transform (const gsl_wavelet * w, + double *data, + size_t tda, size_t size1, size_t size2, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work); + +int gsl_wavelet2d_transform_forward (const gsl_wavelet * w, + double *data, + size_t tda, size_t size1, size_t size2, + gsl_wavelet_workspace * work); + +int gsl_wavelet2d_transform_inverse (const gsl_wavelet * w, + double *data, + size_t tda, size_t size1, size_t size2, + gsl_wavelet_workspace * work); + +int gsl_wavelet2d_nstransform (const gsl_wavelet * w, + double *data, + size_t tda, size_t size1, size_t size2, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work); + +int gsl_wavelet2d_nstransform_forward (const gsl_wavelet * w, + double *data, + size_t tda, size_t size1, size_t size2, + gsl_wavelet_workspace * work); + +int gsl_wavelet2d_nstransform_inverse (const gsl_wavelet * w, + double *data, + size_t tda, size_t size1, size_t size2, + gsl_wavelet_workspace * work); + +int +gsl_wavelet2d_transform_matrix (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work); + +int +gsl_wavelet2d_transform_matrix_forward (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work); + +int +gsl_wavelet2d_transform_matrix_inverse (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work); + + +int +gsl_wavelet2d_nstransform_matrix (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_direction dir, + gsl_wavelet_workspace * work); + +int +gsl_wavelet2d_nstransform_matrix_forward (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work); + +int +gsl_wavelet2d_nstransform_matrix_inverse (const gsl_wavelet * w, + gsl_matrix * a, + gsl_wavelet_workspace * work); + +__END_DECLS + +#endif /* __GSL_WAVELET2D_H__ */ diff --git a/gsl-1.9/wavelet/haar.c b/gsl-1.9/wavelet/haar.c new file mode 100644 index 0000000..8591cfd --- /dev/null +++ b/gsl-1.9/wavelet/haar.c @@ -0,0 +1,81 @@ +/* wavelet/haar.c + * + * Copyright (C) 2004 Ivo Alxneit + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_wavelet.h> + +static const double ch_2[2] = { M_SQRT1_2, M_SQRT1_2 }; +static const double cg_2[2] = { M_SQRT1_2, -(M_SQRT1_2) }; + +static int +haar_init (const double **h1, const double **g1, const double **h2, + const double **g2, size_t * nc, size_t * offset, + const size_t member) +{ + if (member != 2) + { + return GSL_FAILURE; + } + + *h1 = ch_2; + *g1 = cg_2; + *h2 = ch_2; + *g2 = cg_2; + + *nc = 2; + *offset = 0; + + return GSL_SUCCESS; +} + +static int +haar_centered_init (const double **h1, const double **g1, const double **h2, + const double **g2, size_t * nc, size_t * offset, + const size_t member) +{ + if (member != 2) + { + return GSL_FAILURE; + } + + *h1 = ch_2; + *g1 = cg_2; + *h2 = ch_2; + *g2 = cg_2; + + *nc = 2; + *offset = 1; + + return GSL_SUCCESS; +} + +static const gsl_wavelet_type haar_type = { + "haar", + &haar_init +}; + +static const gsl_wavelet_type haar_centered_type = { + "haar-centered", + &haar_centered_init +}; + +const gsl_wavelet_type *gsl_wavelet_haar = &haar_type; +const gsl_wavelet_type *gsl_wavelet_haar_centered = &haar_centered_type; diff --git a/gsl-1.9/wavelet/test.c b/gsl-1.9/wavelet/test.c new file mode 100644 index 0000000..997e4aa --- /dev/null +++ b/gsl-1.9/wavelet/test.c @@ -0,0 +1,270 @@ +#include <config.h> +#include <stdio.h> +#include <math.h> + +#include <gsl/gsl_errno.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_test.h> + +#include <gsl/gsl_wavelet.h> +#include <gsl/gsl_wavelet2d.h> + +#define N_BS 11 + +double urand (void); + +double +urand (void) +{ + static unsigned long int x = 1; + x = (1103515245 * x + 12345) & 0x7fffffffUL; + return x / 2147483648.0; +} + +const size_t member[N_BS] = + { 309, 307, 305, 303, 301, 208, 206, 204, 202, 105, 103 }; + +void +test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member); + +void +test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type); + +int +main (int argc, char **argv) +{ + size_t i, N, stride, tda; + const int S = 1, NS = 2; /* Standard & Non-standard transforms */ + + /* One-dimensional tests */ + + for (N = 1; N <= 16384; N *= 2) + { + for (stride = 1; stride <= 5; stride++) + { + for (i = 0; i < N_BS; i++) + { + test_1d (N, stride, gsl_wavelet_bspline, member[i]); + test_1d (N, stride, gsl_wavelet_bspline_centered, member[i]); + } + + for (i = 4; i <= 20; i += 2) + { + test_1d (N, stride, gsl_wavelet_daubechies, i); + test_1d (N, stride, gsl_wavelet_daubechies_centered, i); + } + + test_1d (N, stride, gsl_wavelet_haar, 2); + test_1d (N, stride, gsl_wavelet_haar_centered, 2); + } + } + + /* Two-dimensional tests */ + + for (N = 1; N <= 64; N *= 2) + { + for (tda = N; tda <= N + 5; tda++) + { + for (i = 0; i < N_BS; i++) + { + test_2d (N, tda, gsl_wavelet_bspline, member[i], S); + test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], S); + + test_2d (N, tda, gsl_wavelet_bspline, member[i], NS); + test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], NS); + } + + for (i = 4; i <= 20; i += 2) + { + test_2d (N, tda, gsl_wavelet_daubechies, i, S); + test_2d (N, tda, gsl_wavelet_daubechies_centered, i, S); + + test_2d (N, tda, gsl_wavelet_daubechies, i, NS); + test_2d (N, tda, gsl_wavelet_daubechies_centered, i, NS); + } + + test_2d (N, tda, gsl_wavelet_haar, 2, S); + test_2d (N, tda, gsl_wavelet_haar_centered, 2, S); + + test_2d (N, tda, gsl_wavelet_haar, 2, NS); + test_2d (N, tda, gsl_wavelet_haar_centered, 2, NS); + } + } + + exit (gsl_test_summary ()); +} + + +void +test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member) +{ + gsl_wavelet_workspace *work; + gsl_vector *v1, *v2, *vdelta; + gsl_vector_view v; + gsl_wavelet *w; + + size_t i; + double *data = malloc (N * stride * sizeof (double)); + + for (i = 0; i < N * stride; i++) + data[i] = 12345.0 + i; + + v = gsl_vector_view_array_with_stride (data, stride, N); + v1 = &(v.vector); + + for (i = 0; i < N; i++) + { + gsl_vector_set (v1, i, urand ()); + } + + v2 = gsl_vector_alloc (N); + gsl_vector_memcpy (v2, v1); + + vdelta = gsl_vector_alloc (N); + + work = gsl_wavelet_workspace_alloc (N); + + w = gsl_wavelet_alloc (T, member); + + gsl_wavelet_transform_forward (w, v2->data, v2->stride, v2->size, work); + gsl_wavelet_transform_inverse (w, v2->data, v2->stride, v2->size, work); + + for (i = 0; i < N; i++) + { + double x1 = gsl_vector_get (v1, i); + double x2 = gsl_vector_get (v2, i); + gsl_vector_set (vdelta, i, fabs (x1 - x2)); + } + + { + double x1, x2; + i = gsl_vector_max_index (vdelta); + x1 = gsl_vector_get (v1, i); + x2 = gsl_vector_get (v2, i); + + gsl_test (fabs (x2 - x1) > N * 1e-15, + "%s(%d), n = %d, stride = %d, maxerr = %g", + gsl_wavelet_name (w), member, N, stride, fabs (x2 - x1)); + } + + if (stride > 1) + { + int status = 0; + + for (i = 0; i < N * stride; i++) + { + if (i % stride == 0) + continue; + + status |= (data[i] != (12345.0 + i)); + } + + gsl_test (status, "%s(%d) other data untouched, n = %d, stride = %d", + gsl_wavelet_name (w), member, N, stride); + } + + gsl_wavelet_workspace_free (work); + gsl_wavelet_free (w); + gsl_vector_free (vdelta); + gsl_vector_free (v2); + free (data); +} + + +void +test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type) +{ + gsl_wavelet_workspace *work; + gsl_matrix *m2; + gsl_wavelet *w; + gsl_matrix *m1; + gsl_matrix *mdelta; + gsl_matrix_view m; + size_t i; + size_t j; + + double *data = malloc (N * tda * sizeof (double)); + + const char * typename; + + typename = (type == 1) ? "standard" : "nonstd" ; + + for (i = 0; i < N * tda; i++) + data[i] = 12345.0 + i; + + m = gsl_matrix_view_array_with_tda (data, N, N, tda); + m1 = &(m.matrix); + + for (i = 0; i < N; i++) + { + for (j = 0; j < N; j++) + { + gsl_matrix_set (m1, i, j, urand()); + } + } + + m2 = gsl_matrix_alloc (N, N); + gsl_matrix_memcpy (m2, m1); + + mdelta = gsl_matrix_alloc (N, N); + + work = gsl_wavelet_workspace_alloc (N); + + w = gsl_wavelet_alloc (T, member); + + switch (type) + { + case 1: + gsl_wavelet2d_transform_matrix_forward (w, m2, work); + gsl_wavelet2d_transform_matrix_inverse (w, m2, work); + break; + case 2: + gsl_wavelet2d_nstransform_matrix_forward (w, m2, work); + gsl_wavelet2d_nstransform_matrix_inverse (w, m2, work); + break; + } + + for (i = 0; i < N; i++) + { + for (j = 0; j < N; j++) + { + double x1 = gsl_matrix_get (m1, i, j); + double x2 = gsl_matrix_get (m2, i, j ); + gsl_matrix_set (mdelta, i, j, fabs (x1 - x2)); + } + } + + { + double x1, x2; + gsl_matrix_max_index (mdelta, &i, &j); + x1 = gsl_matrix_get (m1, i, j); + x2 = gsl_matrix_get (m2, i, j); + + gsl_test (fabs (x2 - x1) > N * 1e-15, + "%s(%d)-2d %s, n = %d, tda = %d, maxerr = %g", + gsl_wavelet_name (w), member, typename, N, tda, fabs (x2 - x1)); + } + + if (tda > N) + { + int status = 0; + + for (i = 0; i < N ; i++) + { + for (j = N; j < tda; j++) + { + status |= (data[i*tda+j] != (12345.0 + (i*tda+j))); + } + } + + gsl_test (status, "%s(%d)-2d %s other data untouched, n = %d, tda = %d", + gsl_wavelet_name (w), member, typename, N, tda); + } + + free (data); + gsl_wavelet_workspace_free (work); + gsl_wavelet_free (w); + gsl_matrix_free (m2); + gsl_matrix_free (mdelta); +} diff --git a/gsl-1.9/wavelet/wavelet.c b/gsl-1.9/wavelet/wavelet.c new file mode 100644 index 0000000..fb51874 --- /dev/null +++ b/gsl-1.9/wavelet/wavelet.c @@ -0,0 +1,132 @@ +/* wavelet/wavelet.c + * + * Copyright (C) 2004 Ivo Alxneit + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_wavelet.h> + +gsl_wavelet * +gsl_wavelet_alloc (const gsl_wavelet_type * T, size_t k) +{ + int status; + + gsl_wavelet *w = (gsl_wavelet *) malloc (sizeof (gsl_wavelet)); + + if (w == NULL) + { + GSL_ERROR_VAL ("failed to allocate space for wavelet struct", + GSL_ENOMEM, 0); + }; + + w->type = T; + + status = (T->init) (&(w->h1), &(w->g1), &(w->h2), &(w->g2), + &(w->nc), &(w->offset), k); + + if (status) + { + free (w); + GSL_ERROR_VAL ("invalid wavelet member", GSL_EINVAL, 0); + } + + return w; +} + +void +gsl_wavelet_free (gsl_wavelet * w) +{ + free (w); +} + +const char * +gsl_wavelet_name (const gsl_wavelet * w) +{ + return w->type->name; +} + + +/* Let's not export this for now (BJG) */ +#if 0 +void +gsl_wavelet_print (const gsl_wavelet * w) +{ + size_t n = w->nc; + size_t i; + + printf ("Wavelet type: %s\n", w->type->name); + + printf + (" h1(%d):%12.8f g1(%d):%12.8f h2(%d):%12.8f g2(%d):%12.8f\n", + 0, w->h1[0], 0, w->g1[0], 0, w->h2[0], 0, w->g2[0]); + + for (i = 1; i < (n < 10 ? n : 10); i++) + { + printf + (" h1(%d):%12.8f g1(%d):%12.8f h2(%d):%12.8f g2(%d):%12.8f\n", + i, w->h1[i], i, w->g1[i], i, w->h2[i], i, w->g2[i]); + } + + for (; i < n; i++) + { + printf + ("h1(%d):%12.8f g1(%d):%12.8f h2(%d):%12.8f g2(%d):%12.8f\n", + i, w->h1[i], i, w->g1[i], i, w->h2[i], i, w->g2[i]); + } +} +#endif + +gsl_wavelet_workspace * +gsl_wavelet_workspace_alloc (size_t n) +{ + gsl_wavelet_workspace *work; + + if (n == 0) + { + GSL_ERROR_VAL ("length n must be positive integer", GSL_EDOM, 0); + } + + work = (gsl_wavelet_workspace *) malloc (sizeof (gsl_wavelet_workspace)); + + if (work == NULL) + { + GSL_ERROR_VAL ("failed to allocate struct", GSL_ENOMEM, 0); + } + + work->n = n; + work->scratch = (double *) malloc (n * sizeof (double)); + + if (work->scratch == NULL) + { + /* error in constructor, prevent memory leak */ + free (work); + GSL_ERROR_VAL ("failed to allocate scratch space", GSL_ENOMEM, 0); + } + + return work; +} + +void +gsl_wavelet_workspace_free (gsl_wavelet_workspace * work) +{ + /* release scratch space */ + free (work->scratch); + work->scratch = NULL; + free (work); +} |