add installer

commit : UWP build error fixed

.gitignore

@@ -0,0 +1,42 @@
+#Folders Ignore
+Bin/
+Testing/
+.vscode/
+
+#cmake ignore
+CMakeCache.txt
+CMakeFiles/
+CMakeScripts/
+Makefile/
+cmake_install.cmake
+install_manifest.txt
+CTestTestfile.cmake
+Makefile
+
+#yml
+*.yml
+
+#Compiled
+*.a
+*.o
+*.lib
+*.so
+*.exe
+*.dll
+
+#Eclipse
+.project
+.cproject
+language.settings.xml
+
+#tcl
+*.tcl
+
+#wavelib-specific
+denoised.txt
+
+test/
+build/
+.vscode/
+.DS_Store
+*.sh

.travis.yml

@@ -0,0 +1,33 @@
+sudo: false
+language: c
+os:
+    - linux
+    - osx
+
+compiler:
+  - gcc
+  - clang
+ 
+addons:
+  apt:
+    sources:
+        - ubuntu-toolchain-r-test
+
+matrix:
+  allow_failures:
+    - compiler: clang
+
+before_install:
+    # linux prereqisite packages
+    - if [ "$TRAVIS_OS_NAME" == "linux" ]; then wget --no-check-certificate https://www.cmake.org/files/v3.2/cmake-3.2.3-Linux-x86_64.tar.gz; fi
+    - if [ "$TRAVIS_OS_NAME" == "linux" ]; then tar -xzvf cmake-3.2.3-Linux-x86_64.tar.gz; fi
+    - if [ "$TRAVIS_OS_NAME" == "linux" ]; then export PATH=$PWD/cmake-3.2.3-Linux-x86_64/bin:$PATH; fi
+
+before_script:
+
+
+script:
+    - mkdir build.ci && cd build.ci
+    - cmake .. -DBUILD_UT=ON -DCMAKE_BUILD_TYPE=$BUILD_CONFIG -DUSE_STATIC_BOOST=YES 
+    - cmake --build .
+    - ctest -VV

CMakeLists.txt

@@ -0,0 +1,15 @@
+cmake_minimum_required(VERSION 3.15.2)
+# 设置工程名称和语言
+project(WaveLib VERSION 1.0)
+# 添加配置配件编写的函数
+include(CMakePackageConfigHelpers)
+
+message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME})
+# CMake默认的安装路径 Windows下为C:/Program\ Files/${Project_Name} Linux/Unix下为/usr/local
+message(STATUS "Install prefix: " ${CMAKE_INSTALL_PREFIX})
+# CMake默认的变异类型为空
+message(STATUS "Build type: " ${CMAKE_BUILD_TYPE})
+
+# 添加源文件地址
+add_subdirectory(src)
+add_subdirectory(auxiliary)

COPYRIGHT

@@ -1,4 +1,5 @@
 Copyright (c) 2014, Rafat Hussain
+Copyright (c) 2016, Holger Nahrstaedt
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

README.md

@@ -1,19 +1,29 @@
+[![Build Status](https://travis-ci.org/rafat/wavelib.svg?branch=master)](https://travis-ci.org/rafat/wavelib)
+
 wavelib
 =======
 
-C Implementation of Wavelet Transform (DWT,SWT and MODWT)
+C Implementation of Discrete Wavelet Transform (DWT,SWT and MODWT), Continuous Wavelet transform (CWT) and Discrete Packet Transform ( Full Tree Decomposition and Best Basis DWPT).
 
-Methods Implemented
+Discrete Wavelet Transform Methods Implemented
 
-DWT/IDWT A decimated Discrete Wavelet Transform implementation using implicit signal extension and up/downsampling so it is a fast implementation. A FFT based implementation is optional but will not be usually needed. Both periodic and symmetric options are available.
+DWT/IDWT and DWT2/IDWT2 A decimated Discrete Wavelet Transform implementation using implicit signal extension and up/downsampling so it is a fast implementation. A FFT based implementation is optional but will not be usually needed. Both periodic and symmetric options are available.
 
-SWT/ISWT Stationary Wavelet Transform. It works only for signal lengths that are multiples of 2^J where J is the number of decomposition levels. For signals of other lengths see MODWT implementation.
+SWT/ISWT and SWT2/ISWT2 Stationary Wavelet Transform. It works only for signal lengths that are multiples of 2^J where J is the number of decomposition levels. For signals of other lengths see MODWT implementation.
 
-MODWT/IMODWT Maximal Overlap Discrete Wavelet Transform is another undecimated transform. It is implemented for signals of any length but only orthogonal wavelets (Daubechies, Symlets and Coiflets) can be deployed. This implementation is based on the method laid out in "Wavelet Methods For Wavelet Analysis" by Donald Percival and Andrew Walden. 
+MODWT/IMODWT and MODWT2/IMODWT2 Maximal Overlap Discrete Wavelet Transform is another undecimated transform. It is implemented for signals of any length but only orthogonal wavelets (Daubechies, Symlets and Coiflets) can be deployed. This implementation is based on the method laid out in "Wavelet Methods For Wavelet Analysis" by Donald Percival and Andrew Walden. 
+
+Discrete Wavelet Packet Transform Methods Implemented
+
+WTREE A Fully Decimated Wavelet Tree Decomposition. This is a highly redundant transform and retains all coefficients at each node. This is not recommended for compression and denoising applications.
+
+DWPT/IDWPT Is a derivative of WTREE method which retains coefficients based on entropy methods. This is a non-redundant transform and output length is of the same order as the input.
+
+CWT/ICWT C translation ( with some modifications) of Continuous Wavelet  Transform Software provided by C. Torrence and G. Compo, and is available at URL: http://atoc.colorado.edu/research/wavelets/'. A generalized Inverse Transform with approximate reconstruction is also added.
 
 Documentation Available at - https://github.com/rafat/wavelib/wiki
 
-Live Demo (Emscripten) - http://rafat.github.io/wavelib/
+Live Demo of 1D DWT and 1D CWT (Emscripten) - http://rafat.github.io/wavelib/
 
 License - BSD 3-Clause
 

WaveLibConfig.cmake.in

@@ -0,0 +1,19 @@
+@PACKAGE_INIT@
+
+set(@PROJECT_NAME@_VERSION "@PROJECT_VERSION@")
+set_and_check(@PROJECT_NAME@_INSTALL_PREFIX "${PACKAGE_PREFIX_DIR}")
+set_and_check(@PROJECT_NAME@_INC_DIR "${PACKAGE_PREFIX_DIR}/include")
+set_and_check(@PROJECT_NAME@_INCULDE_DIR "${PACKAGE_PREFIX_DIR}/include")
+set_and_check(@PROJECT_NAME@_LIB_DIR "${PACKAGE_PREFIX_DIR}/lib")
+set_and_check(@PROJECT_NAME@_LIBRARY_DIR "${PACKAGE_PREFIX_DIR}/lib")
+
+set(@PROJECT_NAME@_LIB wavelib)
+set(@PROJECT_NAME@_LIBRARY wavelib)
+
+set(WauxLib_LIB wauxlib)
+set(WauxLib_LIBRARY wauxlib)
+
+# include target information
+include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake")
+
+check_required_components(@PROJECT_NAME@)

appveyor.yml

@@ -0,0 +1,16 @@
+os: Visual Studio 2015
+
+platform: x64
+
+environment: 
+  BOOST_ROOT: C:\Libraries\boost_1_59_0 
+  BOOST_LIBRARYDIR: C:\Libraries\boost_1_59_0\lib64-msvc-14.0
+
+build_script: 
+  - mkdir build 
+  - cd build 
+  - cmake -G "Visual Studio 14 2015 Win64" -DUSE_STATIC_BOOST=NO  ..
+  - cmake --build . --config Debug
+
+test_script: 
+  - ctest -VV -C Debug 

auxiliary/CMakeLists.txt

@@ -0,0 +1,55 @@
+# 设定源文件文件夹
+aux_source_directory(. WAUX_SRC)
+
+# 以下部分为库的编译
+# 注意目标名必须唯一 所以不能直接生成相同名称的动态库与静态库
+# 注意此处不必为目标名称添加lib前缀和相应后缀，cmake会自行添加
+add_library(wauxlib SHARED ${WAUX_SRC})
+# 首先添加静态库的生成命令
+add_library(wauxlib_static STATIC ${WAUX_SRC})
+# 设置静态库的输出名称从而获得与动态库名称相同的静态库
+set_target_properties(wauxlib_static PROPERTIES OUTPUT_NAME "wauxlib")
+# 设置输出目标属性以同时输出动态库与静态库
+set_target_properties(wauxlib PROPERTIES CLEAN_DIRECT_OUTPUT 1)
+set_target_properties(wauxlib_static PROPERTIES CLEAN_DIRECT_OUTPUT 1)
+# 设置动态库的版本号
+set_target_properties(wauxlib PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR})
+# 设置动态库的运行搜索地址
+set_target_properties(wauxlib PROPERTIES INSTALL_RPATH /usr/local/lib)
+set_target_properties(wauxlib_static PROPERTIES INSTALL_RPATH /usr/local/lib)
+
+target_link_libraries(wauxlib PUBLIC wavelib)
+target_link_libraries(wauxlib_static PUBLIC wavelib_static)
+
+# 设置库文件的输出地址
+set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
+
+# 设置编译选项
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3")
+
+set(CONFIG_FILE_PATH lib/cmake/${PROJECT_NAME})
+
+configure_package_config_file(${PROJECT_SOURCE_DIR}/${PROJECT_NAME}Config.cmake.in 
+    ${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake
+    INSTALL_DESTINATION ${CONFIG_FILE_PATH})
+
+write_basic_package_version_file(${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake
+        VERSION ${PROJECT_VERSION}
+        COMPATIBILITY SameMajorVersion)
+
+# 库的安装命令
+if(WIN32)
+    install(TARGETS wauxlib DESTINATION lib)
+    install(TARGETS wauxlib_static DESTINATION lib)
+else()
+    install(TARGETS wauxlib wauxlib_static
+        EXPORT ${PROJECT_NAME}Targets
+        LIBRARY DESTINATION lib
+        ARCHIVE DESTINATION lib)
+    install(EXPORT ${PROJECT_NAME}Targets
+            DESTINATION ${CONFIG_FILE_PATH})
+    install(FILES
+        ${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake
+        ${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake
+        DESTINATION ${CONFIG_FILE_PATH})
+endif()

auxiliary/denoise.c

@@ -0,0 +1,489 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "waux.h"
+#include "../header/wauxlib.h"
+
+denoise_object denoise_init(int length, int J,const char* wname) {
+	denoise_object obj = NULL;
+
+	obj = (denoise_object)malloc(sizeof(struct denoise_set) +sizeof(double));
+
+	obj->N = length;
+	obj->J = J;
+
+	strcpy(obj->wname,wname);
+
+	//Set Default Values
+	strcpy(obj->dmethod,"sureshrink");
+	strcpy(obj->ext,"sym");
+	strcpy(obj->level,"all");
+	strcpy(obj->thresh,"soft");
+	strcpy(obj->wmethod,"dwt");
+	strcpy(obj->cmethod,"direct");
+
+	return obj;
+}
+
+void visushrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised) {
+	int filt_len,iter,i,dlen,dwt_len,sgn, MaxIter,it;
+	double sigma,td,tmp;
+	wave_object wave;
+	wt_object wt;
+	double *dout,*lnoise;
+
+	wave = wave_init(wname);
+	
+	filt_len = wave->filtlength;
+	
+	MaxIter = (int) (log((double)N / ((double)filt_len - 1.0)) / log(2.0));
+
+	if (J > MaxIter) {
+		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter);
+		exit(-1);
+	}
+
+	wt = wt_init(wave,method,N,J);
+	if(!strcmp(method,"dwt")) {
+		setDWTExtension(wt,ext);
+		dwt(wt,signal);
+	} else if(!strcmp(method,"swt")) {
+		swt(wt,signal);
+	} else {
+		printf("Acceptable WT methods are - dwt,swt and modwt\n");
+		exit(-1);
+	}
+
+	lnoise = (double*)malloc(sizeof(double) * J);
+
+	//Set sigma
+
+	iter = wt->length[0];
+	dlen = wt->length[J];
+
+	dout = (double*)malloc(sizeof(double) * dlen);
+
+	if(!strcmp(level,"first")) {
+		for (i = 1; i < J; ++i) {
+			iter += wt->length[i];
+		}
+
+		for(i = 0; i < dlen;++i) {
+			dout[i] = fabs(wt->output[iter+i]);
+		}
+
+		sigma = median(dout,dlen) / 0.6745;
+		for(it = 0; it < J;++it) {
+			lnoise[it] = sigma;
+		}
+	} else if(!strcmp(level,"all")){
+		for(it = 0; it < J;++it) {
+			dlen = wt->length[it+1];
+			for(i = 0; i < dlen;++i) {
+				dout[i] = fabs(wt->output[iter+i]);
+			}
+			sigma = median(dout,dlen) / 0.6745;
+			lnoise[it] = sigma;
+			iter += dlen;
+		}
+
+	} else {
+		printf("Acceptable Noise estimation level values are - first and all \n");
+		exit(-1);
+	}
+
+	dwt_len = wt->outlength;
+	iter = wt->length[0];
+	for(it = 0; it < J;++it) {
+		sigma = lnoise[it];
+		dlen = wt->length[it+1];
+		td = sqrt(2.0 * log(dwt_len)) * sigma;
+
+		if(!strcmp(thresh,"hard")) {
+			for(i = 0; i < dlen;++i) {
+				if (fabs(wt->output[iter+i]) < td) {
+					wt->output[iter+i] = 0;
+				}
+			}
+		} else if(!strcmp(thresh,"soft")) {
+			for(i = 0; i < dlen;++i) {
+					if (fabs(wt->output[iter + i]) < td) {
+						wt->output[iter+i] = 0;
+					} else {
+						sgn = wt->output[iter+i] >= 0 ? 1 : -1;
+						tmp = sgn * (fabs(wt->output[iter+i]) - td);
+						wt->output[iter+i] = tmp;
+					}
+				}
+		}
+
+		iter += wt->length[it+1];
+	}
+
+	if(!strcmp(method,"dwt")) {
+		idwt(wt,denoised);
+	} else if(!strcmp(method,"swt")) {
+		iswt(wt,denoised);
+	}
+
+	free(dout);
+	free(lnoise);
+	wave_free(wave);
+	wt_free(wt);
+}
+
+void sureshrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised) {
+	int filt_len,i,it,len,dlen,dwt_len,min_index,sgn, MaxIter,iter;
+	double sigma,norm,td,tv,te,ct,thr,temp,x_sum;
+	wave_object wave;
+	wt_object wt;
+	double *dout,*risk,*dsum,*lnoise;
+
+	wave = wave_init(wname);
+
+	filt_len = wave->filtlength;
+
+	MaxIter = (int) (log((double)N / ((double)filt_len - 1.0)) / log(2.0));
+	// Depends on J
+	if (J > MaxIter) {
+		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter);
+		exit(-1);
+	}
+
+	wt = wt_init(wave,method,N,J);
+
+	if(!strcmp(method,"dwt")) {
+		setDWTExtension(wt,ext);
+		dwt(wt,signal);
+	} else if(!strcmp(method,"swt")) {
+		swt(wt,signal);
+	} else {
+		printf("Acceptable WT methods are - dwt and swt\n");
+		exit(-1);
+	}
+
+	len = wt->length[0];
+	dlen = wt->length[J];
+
+	dout = (double*)malloc(sizeof(double) * dlen);
+	risk = (double*)malloc(sizeof(double) * dlen);
+	dsum = (double*)malloc(sizeof(double) * dlen);
+	lnoise = (double*)malloc(sizeof(double) * J);
+
+	iter = wt->length[0];
+
+	if(!strcmp(level,"first")) {
+		for (i = 1; i < J; ++i) {
+			iter += wt->length[i];
+		}
+
+		for(i = 0; i < dlen;++i) {
+			dout[i] = fabs(wt->output[iter+i]);
+		}
+
+		sigma = median(dout,dlen) / 0.6745;
+		for(it = 0; it < J;++it) {
+			lnoise[it] = sigma;
+		}
+	} else if(!strcmp(level,"all")){
+		for(it = 0; it < J;++it) {
+			dlen = wt->length[it+1];
+			for(i = 0; i < dlen;++i) {
+				dout[i] = fabs(wt->output[iter+i]);
+			}
+			sigma = median(dout,dlen) / 0.6745;
+			lnoise[it] = sigma;
+			iter += dlen;
+		}
+
+	} else {
+		printf("Acceptable Noise estimation level values are - first and all \n");
+		exit(-1);
+	}
+
+
+	for(it = 0; it < J;++it) {
+		dwt_len = wt->length[it+1];
+		sigma = lnoise[it];
+
+		if ( sigma < 0.00000001) {
+			td = 0;
+		} else {
+			tv = sqrt(2.0 * log(dwt_len));
+			norm = 0.0;
+			for(i = 0; i < dwt_len;++i) {
+				norm += (wt->output[len+i] *wt->output[len+i] /(sigma*sigma));
+			}
+			te =(norm - (double) dwt_len)/(double) dwt_len;
+			ct = pow(log((double) dwt_len)/log(2.0),1.5)/sqrt((double) dwt_len);
+
+			if (te < ct) {
+				td = tv;
+			} else {
+				x_sum = 0.0;
+
+				for(i = 0; i < dwt_len;++i) {
+					dout[i] = fabs(wt->output[len+i]/sigma);
+				}
+
+				qsort(dout, dwt_len, sizeof(double), compare_double);
+				for(i = 0; i < dwt_len;++i) {
+					dout[i] = (dout[i]*dout[i]);
+					x_sum += dout[i];
+					dsum[i] = x_sum;
+				}
+
+				for(i = 0;i < dwt_len;++i) {
+					risk[i] = ((double)dwt_len - 2 * ((double)i + 1) +dsum[i] +
+							dout[i]*((double)dwt_len - 1 -(double) i))/(double)dwt_len;
+				}
+				min_index = minindex(risk,dwt_len);
+				thr = sqrt(dout[min_index]);
+				td = thr < tv ? thr : tv;
+			}
+		}
+
+		td = td * sigma;
+
+		if(!strcmp(thresh,"hard")) {
+			for(i = 0; i < dwt_len;++i) {
+				if (fabs(wt->output[len+i]) < td) {
+					wt->output[len+i] = 0;
+				}
+			}
+		} else if(!strcmp(thresh,"soft")) {
+			for(i = 0; i < dwt_len;++i) {
+					if (fabs(wt->output[len + i]) < td) {
+						wt->output[len+i] = 0;
+					} else {
+						sgn = wt->output[len+i] >= 0 ? 1 : -1;
+						temp = sgn * (fabs(wt->output[len+i]) - td);
+						wt->output[len+i] = temp;
+					}
+				}
+		}
+
+		len += wt->length[it+1];
+	}
+
+	if(!strcmp(method,"dwt")) {
+		idwt(wt,denoised);
+	} else if(!strcmp(method,"swt")) {
+		iswt(wt,denoised);
+	}
+	free(dout);
+	free(dsum);
+	free(risk);
+	free(lnoise);
+	wave_free(wave);
+	wt_free(wt);
+}
+
+void modwtshrink(double *signal, int N, int J, const char *wname, const char *cmethod, const char *ext, const char *thresh, double *denoised) {
+	int filt_len, iter, i, dlen, sgn, MaxIter, it;
+	double sigma, td, tmp, M, llen;
+	wave_object wave;
+	wt_object wt;
+	double *dout, *lnoise;
+
+	wave = wave_init(wname);
+
+	filt_len = wave->filtlength;
+
+	MaxIter = (int)(log((double)N / ((double)filt_len - 1.0)) / log(2.0));
+
+	if (J > MaxIter) {
+		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
+		exit(-1);
+	}
+
+	wt = wt_init(wave, "modwt", N, J);
+
+	if (!strcmp(ext, "sym") && !strcmp(cmethod,"fft")) {
+		setWTConv(wt, "fft");
+		setDWTExtension(wt, "sym");
+	}
+	else if (!strcmp(ext, "sym") && !strcmp(cmethod, "direct")) {
+		printf("Symmetric Extension is not available for direct method");
+		exit(-1);
+	}
+	else if (!strcmp(ext, "per") && !strcmp(cmethod, "direct"))  {
+		setWTConv(wt, "direct");
+		setDWTExtension(wt, "per");
+	}
+	else if (!strcmp(ext, "per") && !strcmp(cmethod, "fft"))  {
+		setWTConv(wt, "fft");
+		setDWTExtension(wt, "per");
+	}
+	else {
+		printf("Signal extension can be either per or sym");
+		exit(-1);
+	}
+
+	modwt(wt, signal);
+
+	lnoise = (double*)malloc(sizeof(double)* J);
+
+	//Set sigma
+
+	iter = wt->length[0];
+	dlen = wt->length[J];
+	dout = (double*)malloc(sizeof(double)* dlen);
+
+	for (it = 0; it < J; ++it) {
+		dlen = wt->length[it + 1];
+		for (i = 0; i < dlen; ++i) {
+			dout[i] = fabs(wt->output[iter + i]);
+		}
+
+		sigma = sqrt(2.0) * median(dout, dlen) / 0.6745;
+		lnoise[it] = sigma;
+		iter += dlen;
+	}
+
+	M = pow(2.0,J);
+	llen = log((double)wt->modwtsiglength);
+	// Thresholding
+
+	iter = wt->length[0];
+	for (it = 0; it < J; ++it) {
+		sigma = lnoise[it];
+		dlen = wt->length[it + 1];
+		td = sqrt(2.0 * llen / M) * sigma;
+
+		if (!strcmp(thresh, "hard")) {
+			for (i = 0; i < dlen; ++i) {
+				if (fabs(wt->output[iter + i]) < td) {
+					wt->output[iter + i] = 0;
+				}
+			}
+		}
+		else if (!strcmp(thresh, "soft")) {
+			for (i = 0; i < dlen; ++i) {
+				if (fabs(wt->output[iter + i]) < td) {
+					wt->output[iter + i] = 0;
+				}
+				else {
+					sgn = wt->output[iter + i] >= 0 ? 1 : -1;
+					tmp = sgn * (fabs(wt->output[iter + i]) - td);
+					wt->output[iter + i] = tmp;
+				}
+			}
+		}
+
+		iter += wt->length[it + 1];
+		M /= 2.0;
+	}
+
+	imodwt(wt, denoised);
+
+	free(dout);
+	free(lnoise);
+	wave_free(wave);
+	wt_free(wt);
+}
+
+
+void denoise(denoise_object obj, double *signal,double *denoised) {
+	if(!strcmp(obj->dmethod,"sureshrink")) {
+		if (!strcmp(obj->wmethod, "modwt")) {
+			printf("sureshrink method only works with swt and dwt. Please use setDenoiseWTMethod to set the correct method\n");
+			exit(-1);
+		}
+		sureshrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);
+	} else if(!strcmp(obj->dmethod,"visushrink")) {
+		if (!strcmp(obj->wmethod, "modwt")) {
+			printf("visushrink method only works with swt and dwt. Please use setDenoiseWTMethod to set the correct method\n");
+			exit(-1);
+		}
+		visushrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);;
+	} else if(!strcmp(obj->dmethod,"modwtshrink")) {
+		if (strcmp(obj->wmethod, "modwt")) {
+			printf("modwtshrink method only works with modwt. Please use setDenoiseWTMethod to set the correct method\n");
+			exit(-1);
+		}
+		modwtshrink(signal,obj->N,obj->J,obj->wname,obj->cmethod,obj->ext,obj->thresh,denoised);;
+	} else {
+		printf("Acceptable Denoising methods are - sureshrink and visushrink\n");
+		exit(-1);
+	}
+}
+
+void setDenoiseMethod(denoise_object obj, const char *dmethod) {
+	if (!strcmp(dmethod, "sureshrink")) {
+		strcpy(obj->dmethod, "sureshrink");
+	}
+	else if (!strcmp(dmethod, "visushrink")) {
+		strcpy(obj->dmethod, "visushrink");
+	}
+	 else if (!strcmp(dmethod, "modwtshrink")) {
+		strcpy(obj->dmethod, "modwtshrink");
+	}
+	else {
+		printf("Acceptable Denoising methods are - sureshrink, visushrink and modwtshrink\n");
+		exit(-1);
+	}
+}
+
+void setDenoiseWTMethod(denoise_object obj, const char *wmethod) {
+	if (!strcmp(wmethod, "dwt")) {
+		strcpy(obj->wmethod, "dwt");
+	}
+	else if (!strcmp(wmethod, "swt")) {
+		strcpy(obj->wmethod, "swt");
+	}
+	 else if (!strcmp(wmethod, "modwt")) {
+		strcpy(obj->wmethod, "modwt");
+	}
+	else {
+		printf("Wavelet decomposition method can be one of dwt, modwt or swt.\n");
+		exit(-1);
+	}
+}
+
+void setDenoiseWTExtension(denoise_object obj, const char *extension) {
+	if (!strcmp(extension, "sym")) {
+		strcpy(obj->ext, "sym");
+	}
+	else if (!strcmp(extension, "per")) {
+		strcpy(obj->ext, "per");
+	}
+	else {
+		printf("Signal extension can be either per or sym");
+		exit(-1);
+	}
+}
+
+void setDenoiseParameters(denoise_object obj, const char *thresh,const char *level) {
+
+	//Set thresholding
+	if (!strcmp(thresh, "soft")) {
+		strcpy(obj->thresh, "soft");
+	}
+	else if (!strcmp(thresh, "hard")) {
+		strcpy(obj->thresh, "hard");
+	}
+	else {
+		printf("Thresholding Method - soft or hard");
+		exit(-1);
+	}
+
+	// Set Noise estimation at the first level or at all levels
+
+	if (!strcmp(level, "first")) {
+		strcpy(obj->level, "first");
+	}
+	else if (!strcmp(level, "all")) {
+		strcpy(obj->level, "all");
+	}
+	else {
+		printf("Noise Estimation at level - first or all");
+		exit(-1);
+	}
+}
+
+void denoise_free(denoise_object object) {
+	free(object);
+}

auxiliary/waux.c

@@ -0,0 +1,324 @@
+#include "../header/wauxlib.h"
+#include "waux.h"
+
+int compare_double(const void* a, const void* b)
+{
+    double arg1 = *(const double*)a;
+    double arg2 = *(const double*)b;
+
+    if (arg1 < arg2) return -1;
+    if (arg1 > arg2) return 1;
+    return 0;
+
+}
+
+double mean(const double* vec, int N) {
+	int i;
+	double m;
+	m = 0.0;
+
+	for (i = 0; i < N; ++i) {
+		m+= vec[i];
+	}
+	m = m / N;
+	return m;
+}
+
+double var(const double* vec, int N) {
+	double v,temp,m;
+	int i;
+	v = 0.0;
+	m = mean(vec,N);
+
+	for (i = 0; i < N; ++i) {
+		temp = vec[i] - m;
+		v+= temp*temp;
+	}
+
+	v = v / N;
+
+	return v;
+
+}
+
+double median(double *x, int N) {
+	double sigma;
+
+    qsort(x, N, sizeof(double), compare_double);
+
+    if ((N % 2) == 0) {
+		sigma = (x[N/2 - 1] + x[N/2] ) / 2.0;
+	} else {
+		sigma = x[N/2];
+	}
+
+	return sigma;
+}
+
+double mad(double *x, int N) {
+	double sigma;
+	int i;
+
+	sigma = median(x,N);
+
+	for(i = 0; i < N;++i) {
+		x[i] = (x[i] - sigma) > 0 ? (x[i] - sigma) : -(x[i] - sigma);
+	}
+
+	sigma = median(x,N);
+
+	return sigma;
+}
+
+int minindex(const double *arr, int N) {
+	double min;
+	int index,i;
+
+	min = DBL_MAX;
+	index = 0;
+	for(i = 0; i < N;++i) {
+		if (arr[i] < min) {
+			min = arr[i];
+			index = i;
+		}
+	}
+
+	return index;
+
+}
+
+void getDWTAppx(wt_object wt, double *appx,int N) {
+	/*
+	Wavelet decomposition is stored as
+	[A(J) D(J) D(J-1) ..... D(1)] in wt->output vector
+
+	Length of A(J) , N = wt->length[0]
+	*/
+	int i;
+
+	for (i = 0; i < N; ++i) {
+		appx[i] = wt->output[i];
+	}
+}
+
+void getDWTDetail(wt_object wt, double *detail, int N, int level) {
+	/*
+	returns Detail coefficents at the jth level where j = J,J-1,...,1
+	and Wavelet decomposition is stored as
+	[A(J) D(J) D(J-1) ..... D(1)] in wt->output vector
+	Use getDWTAppx() to get A(J)
+	Level 1 : Length of D(J), ie N, is stored in wt->length[1]
+	Level 2 :Length of D(J-1), ie N, is stored in wt->length[2]
+	....
+	Level J : Length of D(1), ie N, is stored in wt->length[J]
+	*/
+	int i, iter, J;
+	J = wt->J;
+
+	if (level > J || level < 1) {
+		printf("The decomposition only has 1,..,%d levels", J);
+		exit(-1);
+	}
+
+	iter = wt->length[0];
+
+	for (i = 1; i < J-level; ++i) {
+		iter += wt->length[i];
+	}
+
+	for (i = 0; i < N; ++i) {
+		detail[i] = wt->output[i + iter];
+	}
+}
+
+void getDWTRecCoeff(double *coeff,int *length,const char *ctype,const char *ext, int level, int J,double *lpr,
+		double *hpr,int lf,int siglength,double *reccoeff) {
+
+	int i,j,k,det_len,N,l,m,n,v,t,l2;
+	double *out,*X_lp,*filt;
+	out = (double*)malloc(sizeof(double)* (siglength + 1));
+	l2 = lf / 2;
+	m = -2;
+	n = -1;
+	if (!strcmp(ext, "per")) {
+		if (!strcmp((ctype), "appx")) {
+			det_len = length[0];
+		} else {
+			det_len = length[J - level + 1];
+		}
+
+		N = 2 * length[J];
+
+		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+
+		for (i = 0; i < det_len; ++i) {
+			out[i] = coeff[i];
+		}
+
+		for (j = level; j > 0; --j) {
+
+			//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
+
+			if (!strcmp((ctype), "det") && j == level) {
+				filt = hpr;
+			} else {
+				filt = lpr;
+			}
+
+			//idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp);
+			m = -2;
+			n = -1;
+
+			for (i = 0; i < det_len + l2 - 1; ++i) {
+				m += 2;
+				n += 2;
+				X_lp[m] = 0.0;
+				X_lp[n] = 0.0;
+				for (l = 0; l < l2; ++l) {
+					t = 2 * l;
+					if ((i - l) >= 0 && (i - l) < det_len) {
+						X_lp[m] += filt[t] * out[i - l];
+						X_lp[n] += filt[t + 1] * out[i - l];
+					}
+					else if ((i - l) >= det_len && (i-l) < det_len + lf - 1) {
+						X_lp[m] += filt[t] * out[i - l - det_len];
+						X_lp[n] += filt[t + 1] * out[i - l - det_len];
+					}
+					else if ((i - l) < 0 && (i-l) > -l2) {
+						X_lp[m] += filt[t] * out[det_len + i - l] ;
+						X_lp[n] += filt[t + 1] * out[det_len + i - l];
+					}
+				}
+			}
+
+			for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
+				out[k - lf/2 + 1] = X_lp[k];
+			}
+
+			if (j != 1) {
+				det_len = length[J - j + 2];
+			}
+		}
+
+		free(X_lp);
+
+	}
+	else if (!strcmp(ext, "sym")) {
+		if (!strcmp((ctype), "appx")) {
+			det_len = length[0];
+		} else {
+			det_len = length[J - level + 1];
+		}
+
+		N = 2 * length[J] - 1;
+
+		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+
+		for (i = 0; i < det_len; ++i) {
+			out[i] = coeff[i];
+		}
+
+		for (j = level; j > 0; --j) {
+
+			//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
+
+			if (!strcmp((ctype), "det") && j == level) {
+				filt = hpr;
+			} else {
+				filt = lpr;
+			}
+
+			//idwt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp);
+
+			m = -2;
+			n = -1;
+
+			for (v = 0; v < det_len; ++v) {
+				i = v;
+				m += 2;
+				n += 2;
+				X_lp[m] = 0.0;
+				X_lp[n] = 0.0;
+				for (l = 0; l < lf / 2; ++l) {
+					t = 2 * l;
+					if ((i - l) >= 0 && (i - l) < det_len) {
+						X_lp[m] += filt[t] * out[i - l];
+						X_lp[n] += filt[t + 1] * out[i - l];
+					}
+				}
+			}
+
+			for (k = lf-2; k < 2 * det_len; ++k) {
+				out[k - lf + 2] = X_lp[k];
+			}
+
+
+			if (j != 1) {
+				det_len = length[J - j + 2];
+			}
+		}
+
+		free(X_lp);
+
+	}
+	else {
+			printf("Signal extension can be either per or sym");
+			exit(-1);
+	}
+
+	for (i = 0; i < siglength; ++i) {
+		reccoeff[i] = out[i];
+	}
+
+	free(out);
+
+}
+
+
+void autocovar(const double* vec,int N, double* acov,int M) {
+	double m,temp1,temp2;
+	int i,t;
+	m = mean(vec,N);
+
+	if ( M > N) {
+		M = N-1;
+		printf("\n Lag is greater than the length N of the input vector. It is automatically set to length N - 1.\n");
+		printf("\n The Output Vector only contains N calculated values.");
+	} else if ( M < 0) {
+		M = 0;
+	}
+
+	for(i = 0; i < M;i++) {
+		acov[i] = 0.0;
+		for (t = 0; t < N-i;t++) {
+			temp1 = vec[t] - m;
+			temp2 = vec[t+i] - m;
+			acov[i]+= temp1*temp2;
+		}
+		acov[i] = acov[i] / N;
+
+	}
+
+
+}
+
+void autocorr(const double* vec,int N,double* acorr, int M) {
+	double var;
+	int i;
+	if (M > N) {
+		M = N - 1;
+		printf("\n Lag is greater than the length N of the input vector. It is automatically set to length N - 1.\n");
+		printf("\n The Output Vector only contains N calculated values.");
+	}
+	else if (M < 0) {
+		M = 0;
+	}
+	autocovar(vec,N,acorr,M);
+	var = acorr[0];
+	acorr[0] = 1.0;
+
+	for(i = 1; i < M; i++) {
+		acorr[i] = acorr[i]/var;
+	}
+
+}

auxiliary/waux.h

@@ -0,0 +1,47 @@
+/*
+ * waux.h
+ *
+ *  Created on: Aug 22, 2017
+ *      Author: Rafat Hussain
+ */
+
+#ifndef AUXILIARY_WAUX_H_
+#define AUXILIARY_WAUX_H_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int compare_double(const void* a, const void* b);
+
+double mean(const double* vec, int N);
+
+double var(const double* vec, int N);
+
+double median(double *x, int N);
+
+int minindex(const double *arr, int N);
+
+void getDWTAppx(wt_object wt, double *appx,int N);
+
+void getDWTDetail(wt_object wt, double *detail, int N, int level);
+
+void autocovar(const double* vec,int N,double* acov, int M);
+
+void autocorr(const double* vec,int N,double* acorr, int M);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+
+#endif /* AUXILIARY_WAUX_H_ */

header/wauxlib.h

@@ -0,0 +1,60 @@
+/*
+Copyright (c) 2017, Rafat Hussain
+*/
+#ifndef WAUXLIB_H_
+#define WAUXLIB_H_
+
+#include "wavelib.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct denoise_set* denoise_object;
+
+denoise_object denoise_init(int length, int J,const char* wname);
+
+struct denoise_set{
+	int N; //signal length
+	int J; // Levels of Wavelet decomposition
+	char wname[10]; //Wavelet name
+	char wmethod[10]; //Wavelet decomposition method - dwt or swt
+	char cmethod[10]; //Cnvolution Method - direct or fft . Available only for modwt.
+	// SWT and DWT only use direct method. 
+	char ext[10]; // Signal Extension - sym or per
+	char thresh[10]; // thresholding - soft or hard
+	char level[10]; // Noise Estimation level - first or all
+	char dmethod[20]; //Denoising Method -sureshrink or visushrink
+	//double params[0];
+};
+
+void visushrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised);
+
+void sureshrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised);
+
+void modwtshrink(double *signal, int N, int J, const char *wname, const char *cmethod, const char *ext, const char *thresh, double *denoised);
+
+void denoise(denoise_object obj, double *signal,double *denoised);
+
+void setDenoiseMethod(denoise_object obj, const char *dmethod);
+
+void setDenoiseWTMethod(denoise_object obj, const char *wmethod);
+
+void setDenoiseWTExtension(denoise_object obj, const char *extension);
+
+void setDenoiseParameters(denoise_object obj, const char *thresh,const char *level);
+
+void denoise_free(denoise_object object);
+
+void getDWTRecCoeff(double *coeff,int *length,const char *ctype,const char *ext, int level, int J,double *lpr,
+		double *hpr,int lf,int siglength,double *reccoeff);
+
+double mad(double *x, int N);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAUXLIB_H_ */

header/wavelib.h

@@ -5,13 +5,28 @@
 extern "C" {
 #endif
 
+#if defined(_MSC_VER)
+#pragma warning(disable : 4200)
+#pragma warning(disable : 4996)
+#endif
+
 #ifndef fft_type
 #define fft_type double
 #endif
 
+#ifndef cplx_type
+#define cplx_type double
+#endif
+
+
+typedef struct cplx_t {
+	cplx_type re;
+	cplx_type im;
+} cplx_data;
+
 typedef struct wave_set* wave_object;
 
-wave_object wave_init(char* wname);
+wave_object wave_init(const char* wname);
 
 struct wave_set{
 	char wname[50];
@@ -68,9 +83,35 @@ struct conv_set{
 
 typedef struct wt_set* wt_object;
 
-wt_object wt_init(wave_object wave,char* method, int siglength, int J);
+wt_object wt_init(wave_object wave,const char* method, int siglength, int J);
 
 struct wt_set{
+	wave_object wave;
+	conv_object cobj;
+	char method[10];
+	int siglength;// Length of the original signal.
+	int modwtsiglength; // Modified signal length for MODWT
+	int outlength;// Length of the output DWT vector
+	int lenlength;// Length of the Output Dimension Vector "length"
+	int J; // Number of decomposition Levels
+	int MaxIter;// Maximum Iterations J <= MaxIter
+	int even;// even = 1 if signal is of even length. even = 0 otherwise
+	char ext[10];// Type of Extension used - "per" or "sym"
+	char cmethod[10]; // Convolution Method - "direct" or "FFT"
+
+	int N; //
+	int cfftset;
+	int zpad;
+	int length[102];
+	double *output;
+	double params[0];
+};
+
+typedef struct wtree_set* wtree_object;
+
+wtree_object wtree_init(wave_object wave, int siglength, int J);
+
+struct wtree_set{
 	wave_object wave;
 	conv_object cobj;
 	char method[10];
@@ -81,43 +122,192 @@ struct wt_set{
 	int MaxIter;// Maximum Iterations J <= MaxIter
 	int even;// even = 1 if signal is of even length. even = 0 otherwise
 	char ext[10];// Type of Extension used - "per" or "sym"
-	char cmethod[10]; // Convolution Method - "direct" or "FFT"
 
-	int N; // 
-	int cfftset; 
+	int N; //
+	int nodes;
+	int cfftset;
 	int zpad;
 	int length[102];
 	double *output;
+	int *nodelength;
+	int *coeflength;
 	double params[0];
 };
 
-void dwt11(wt_object wt, double *Vin, int M, double *Wout,
-	double *Vout);
+typedef struct wpt_set* wpt_object;
 
-void dwt(wt_object wt, double *inp);
+wpt_object wpt_init(wave_object wave, int siglength, int J);
+
+struct wpt_set{
+	wave_object wave;
+	conv_object cobj;
+	int siglength;// Length of the original signal.
+	int outlength;// Length of the output DWT vector
+	int lenlength;// Length of the Output Dimension Vector "length"
+	int J; // Number of decomposition Levels
+	int MaxIter;// Maximum Iterations J <= MaxIter
+	int even;// even = 1 if signal is of even length. even = 0 otherwise
+	char ext[10];// Type of Extension used - "per" or "sym"
+	char entropy[20];
+	double eparam;
+
+	int N; //
+	int nodes;
+	int length[102];
+	double *output;
+	double *costvalues;
+	double *basisvector;
+	int *nodeindex;
+	int *numnodeslevel;
+	int *coeflength;
+	double params[0];
+};
+
+
+typedef struct cwt_set* cwt_object;
+
+cwt_object cwt_init(const char* wave, double param, int siglength,double dt, int J);
+
+struct cwt_set{
+	char wave[10];// Wavelet - morl/morlet,paul,dog/dgauss
+	int siglength;// Length of Input Data
+	int J;// Total Number of Scales
+	double s0;// Smallest scale. It depends on the sampling rate. s0 <= 2 * dt for most wavelets
+	double dt;// Sampling Rate
+	double dj;// Separation between scales. eg., scale = s0 * 2 ^ ( [0:N-1] *dj ) or scale = s0 *[0:N-1] * dj
+	char type[10];// Scale Type - Power or Linear
+	int pow;// Base of Power in case type = pow. Typical value is pow = 2
+	int sflag;
+	int pflag;
+	int npad;
+	int mother;
+	double m;// Wavelet parameter param
+	double smean;// Input Signal mean
+
+	cplx_data *output;
+	double *scale;
+	double *period;
+	double *coi;
+	double params[0];
+};
+
+typedef struct wt2_set* wt2_object;
+
+wt2_object wt2_init(wave_object wave, const char* method, int rows, int cols, int J);
+
+struct wt2_set{
+	wave_object wave;
+	char method[10];
+	int rows;// Matrix Number of rows
+	int cols; // Matrix Number of columns
+	int outlength;// Length of the output DWT vector
+	int J; // Number of decomposition Levels
+	int MaxIter;// Maximum Iterations J <= MaxIter
+	char ext[10];// Type of Extension used - "per" or "sym"
+	int coeffaccesslength;
+
+	int N; //
+	int *dimensions;
+	int *coeffaccess;
+	int params[0];
+};
+
+void dwt(wt_object wt, const double *inp);
 
 void idwt(wt_object wt, double *dwtop);
 
-void swt(wt_object wt, double *inp);
+double *getDWTmra(wt_object wt, double *wavecoeffs);
+
+void wtree(wtree_object wt, const double *inp);
+
+void dwpt(wpt_object wt, const double *inp);
+
+void idwpt(wpt_object wt, double *dwtop);
+
+void swt(wt_object wt, const double *inp);
 
 void iswt(wt_object wt, double *swtop);
 
-void modwt(wt_object wt, double *inp);
+double *getSWTmra(wt_object wt, double *wavecoeffs);
+
+void modwt(wt_object wt, const double *inp);
 
 void imodwt(wt_object wt, double *dwtop);
 
-void setDWTExtension(wt_object wt, char *extension);
+double* getMODWTmra(wt_object wt, double *wavecoeffs);
 
-void setWTConv(wt_object wt, char *cmethod);
+void setDWTExtension(wt_object wt, const char *extension);
+
+void setWTREEExtension(wtree_object wt, const char *extension);
+
+void setDWPTExtension(wpt_object wt, const char *extension);
+
+void setDWT2Extension(wt2_object wt, const char *extension);
+
+void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam);
+
+void setWTConv(wt_object wt, const char *cmethod);
+
+int getWTREENodelength(wtree_object wt, int X);
+
+void getWTREECoeffs(wtree_object wt, int X, int Y, double *coeffs, int N);
+
+int getDWPTNodelength(wpt_object wt, int X);
+
+void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N);
+
+void setCWTScales(cwt_object wt, double s0, double dj, const char *type, int power);
+
+void setCWTScaleVector(cwt_object wt, const double *scale, int J, double s0, double dj);
+
+void setCWTPadding(cwt_object wt, int pad);
+
+void cwt(cwt_object wt, const double *inp);
+
+void icwt(cwt_object wt, double *cwtop);
+
+int getCWTScaleLength(int N);
+
+double* dwt2(wt2_object wt, double *inp);
+
+void idwt2(wt2_object wt,double *wavecoeff, double *oup);
+
+double* swt2(wt2_object wt, double *inp);
+
+void iswt2(wt2_object wt, double *wavecoeffs, double *oup);
+
+double* modwt2(wt2_object wt, double *inp);
+
+void imodwt2(wt2_object wt, double *wavecoeff, double *oup);
+
+double* getWT2Coeffs(wt2_object wt,double* wcoeffs, int level,char *type, int *rows, int *cols);
+
+void dispWT2Coeffs(double *A, int row, int col);
 
 void wave_summary(wave_object obj);
 
 void wt_summary(wt_object wt);
 
+void wtree_summary(wtree_object wt);
+
+void wpt_summary(wpt_object wt);
+
+void cwt_summary(cwt_object wt);
+
+void wt2_summary(wt2_object wt);
+
 void wave_free(wave_object object);
 
 void wt_free(wt_object object);
 
+void wtree_free(wtree_object object);
+
+void wpt_free(wpt_object object);
+
+void cwt_free(cwt_object object);
+
+void wt2_free(wt2_object wt);
+
 
 #ifdef __cplusplus
 }
@@ -125,5 +315,3 @@ void wt_free(wt_object object);
 
 
 #endif /* WAVELIB_H_ */
-
-

installer

@@ -0,0 +1,54 @@
+#!/bin/bash
+
+if [[ $# == 0 || ${1} == "help" ]]; then
+	echo "Compiles executables/libraries and maintains installed files. Two tools 'Cmake' and 'stow' are empolyed here. For more information, see https://cmake.org and https://www.gnu.org/software/stow/."
+	echo ""
+	echo "School of Earth Sciences, Zhejiang University"
+	echo "Yi Zhang (yizhang-geo@zju.edu.cn)"
+	echo ""
+	echo "Usage: ./config.sh [option] [Cmake options]"
+	echo ""
+	echo "Options:"
+	echo "(1) configure: Configure Cmake project(s). This option could take extra Cmake options as in <option>=<value>."
+	echo "(2)     build: Build executables/libraries."
+	echo "(3)   install: Install executables/libraries to the directory of CMAKE_INSTALL_PREFIX and sym-links them to the target address. This offers a quick and clean remove of the installed files."
+	echo "(4)     clean: Clean build/ folder(s)."
+	echo "(5) uninstall: Delete the installed files and sym-links."
+	echo "(6)      info: Print out current setups."
+	echo "(7)      help: Show help information."
+	exit 0
+fi
+
+package=wavelib
+address=/opt/stow
+taress=/usr/local
+option="-DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=${address}/${package}"
+
+if [[ $# -gt 1 ]]; then
+	for opt in "$@"; do
+		if [[ ${opt} != "configure" ]]; then
+			option="${option} -D${opt}"
+		fi
+	done
+fi
+
+if [[ ${1} == "configure" && ! -d "build/" ]]; then
+	mkdir build && cd build && cmake .. ${option}
+elif [[ ${1} == "configure" ]]; then
+	cd build && rm -rf * && cmake .. ${option}
+elif [[ ${1} == "build" ]]; then
+	cd build && make
+elif [[ ${1} == "install" ]]; then
+	cd build && sudo make install
+	sudo stow --dir=${address} --target=${taress} -S ${package}
+elif [[ ${1} == "clean" ]]; then
+	rm -rf build/
+elif [[ ${1} == "uninstall" ]]; then
+	sudo stow --dir=${address} --target=${taress} -D ${package}
+	sudo rm -rf ${address}/${package}
+elif [[ ${1} == "info" ]]; then
+	echo "package name:" ${package}
+	echo "stow address:" ${address}
+	echo "target address:" ${taress}
+	echo "Cmake options:" ${option}
+fi

src/CMakeLists.txt

@@ -0,0 +1,55 @@
+# 设定源文件文件夹
+aux_source_directory(. WAVE_SRC)
+
+# 以下部分为库的编译
+# 注意目标名必须唯一 所以不能直接生成相同名称的动态库与静态库
+# 注意此处不必为目标名称添加lib前缀和相应后缀，cmake会自行添加
+add_library(wavelib SHARED ${WAVE_SRC})
+# 首先添加静态库的生成命令
+add_library(wavelib_static STATIC ${WAVE_SRC})
+# 设置静态库的输出名称从而获得与动态库名称相同的静态库
+set_target_properties(wavelib_static PROPERTIES OUTPUT_NAME "wavelib")
+# 设置输出目标属性以同时输出动态库与静态库
+set_target_properties(wavelib PROPERTIES CLEAN_DIRECT_OUTPUT 1)
+set_target_properties(wavelib_static PROPERTIES CLEAN_DIRECT_OUTPUT 1)
+# 设置动态库的版本号
+set_target_properties(wavelib PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR})
+# 设置动态库的运行搜索地址
+set_target_properties(wavelib PROPERTIES INSTALL_RPATH /usr/local/lib)
+set_target_properties(wavelib_static PROPERTIES INSTALL_RPATH /usr/local/lib)
+
+# 设置库文件的输出地址
+set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
+
+# 设置编译选项
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3")
+
+set(CONFIG_FILE_PATH lib/cmake/${PROJECT_NAME})
+
+configure_package_config_file(${PROJECT_SOURCE_DIR}/${PROJECT_NAME}Config.cmake.in 
+    ${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake
+    INSTALL_DESTINATION ${CONFIG_FILE_PATH})
+
+write_basic_package_version_file(${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake
+        VERSION ${PROJECT_VERSION}
+        COMPATIBILITY SameMajorVersion)
+
+# 库的安装命令
+if(WIN32)
+    install(TARGETS wavelib DESTINATION lib)
+    install(TARGETS wavelib_static DESTINATION lib)
+else()
+    install(TARGETS wavelib wavelib_static
+        EXPORT ${PROJECT_NAME}Targets
+        LIBRARY DESTINATION lib
+        ARCHIVE DESTINATION lib)
+    install(EXPORT ${PROJECT_NAME}Targets
+            DESTINATION ${CONFIG_FILE_PATH})
+    install(FILES
+        ${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake
+        ${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake
+        DESTINATION ${CONFIG_FILE_PATH})
+endif()
+# 头文件安装命令
+install(FILES ../header/wavelib.h DESTINATION include)
+install(FILES ../header/wauxlib.h DESTINATION include)

src/conv.c

@@ -204,5 +204,7 @@ void conv_fft(const conv_object obj,fft_type *inp1,fft_type *inp2,fft_type *oup)
 
 
 void free_conv(conv_object object) {
+	free_real_fft(object->fobj);
+	free_real_fft(object->iobj);
 	free(object);
 }

src/conv.h

@@ -17,18 +17,8 @@ extern "C" {
 #define MIN(a,b) (((a)<(b))?(a):(b))
 #define MAX(a,b) (((a)>(b))?(a):(b))
 
-typedef struct conv_set* conv_object;
-
 conv_object conv_init(int N, int L);
 
-struct conv_set{
-	fft_real_object fobj;
-	fft_real_object iobj;
-	int ilen1;
-	int ilen2;
-	int clen;
-};
-
 int factorf(int M);
 
 int findnext(int M);

src/cwt.c

@@ -0,0 +1,398 @@
+/*
+  Copyright (c) 2015, Rafat Hussain
+*/
+/*
+This code is a C translation ( with some modifications) of Wavelet Software provided by
+C. Torrence and G. Compo, and is available at URL: http://atoc.colorado.edu/research/wavelets/''.
+*/
+
+#include "cwt.h"
+
+double factorial(int N) {
+	static const double fact[41] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000,
+		20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000, 51090942171709440000.0, 1124000727777607680000.0,
+		25852016738884976640000.0, 620448401733239439360000.0, 15511210043330985984000000.0, 403291461126605635584000000.0, 10888869450418352160768000000.0,
+		304888344611713860501504000000.0, 8841761993739701954543616000000.0, 265252859812191058636308480000000.0, 8222838654177922817725562880000000.0,
+		263130836933693530167218012160000000.0, 8683317618811886495518194401280000000.0, 295232799039604140847618609643520000000.0, 10333147966386144929666651337523200000000.0,
+		371993326789901217467999448150835200000000.0, 13763753091226345046315979581580902400000000.0, 523022617466601111760007224100074291200000000.0,
+		20397882081197443358640281739902897356800000000.0, 815915283247897734345611269596115894272000000000.0 };
+
+	if (N > 40 || N < 0) {
+		printf("This program is only valid for 0 <= N <= 40 \n");
+		return -1.0;
+	}
+
+	return fact[N];
+
+}
+static void wave_function(int nk, double dt,int mother, double param,double scale1, double *kwave, double pi,double *period1,
+	double *coi1, fft_data *daughter) {
+
+	double norm, expnt, fourier_factor;
+	int k, m;
+	double temp;
+	int sign,re;
+
+
+	if (mother == 0) {
+		//MORLET
+		if (param < 0.0) {
+			param = 6.0;
+		}
+		norm = sqrt(2.0*pi*scale1 / dt)*pow(pi,-0.25);
+
+		for (k = 1; k <= nk / 2 + 1; ++k) {
+			temp = (scale1*kwave[k-1] - param);
+			expnt = -0.5 * temp * temp;
+			daughter[k - 1].re = norm * exp(expnt);
+			daughter[k - 1].im = 0.0;
+		}
+		for (k = nk / 2 + 2; k <= nk; ++k) {
+			daughter[k - 1].re = daughter[k - 1].im = 0.0;
+		}
+		fourier_factor = (4.0*pi) / (param + sqrt(2.0 + param*param));
+		*period1 = scale1*fourier_factor;
+		*coi1 = fourier_factor / sqrt(2.0);
+	}
+	else if (mother == 1) {
+		// PAUL
+		if (param < 0.0) {
+			param = 4.0;
+		}
+		m = (int)param;
+		norm = sqrt(2.0*pi*scale1 / dt)*(pow(2.0,(double)m) / sqrt((double)(m*factorial(2 * m - 1))));
+		for (k = 1; k <= nk / 2 + 1; ++k) {
+			temp = scale1 * kwave[k - 1];
+			expnt = - temp;
+			daughter[k - 1].re = norm * pow(temp,(double)m) * exp(expnt);
+			daughter[k - 1].im = 0.0;
+		}
+		for (k = nk / 2 + 2; k <= nk; ++k) {
+			daughter[k - 1].re = daughter[k - 1].im = 0.0;
+		}
+		fourier_factor = (4.0*pi) / (2.0 * m + 1.0);
+		*period1 = scale1*fourier_factor;
+		*coi1 = fourier_factor * sqrt(2.0);
+	}
+	else if (mother == 2) {
+		if (param < 0.0) {
+			param = 2.0;
+		}
+		m = (int)param;
+
+		if (m % 2 == 0) {
+			re = 1;
+		}
+		else {
+			re = 0;
+		}
+
+		if (m % 4 == 0 || m % 4 == 1) {
+			sign = -1;
+		}
+		else {
+			sign = 1;
+		}
+
+
+		norm = sqrt(2.0*pi*scale1 / dt)*sqrt(1.0 / cwt_gamma(m + 0.50));
+		norm *= sign;
+
+		if (re == 1) {
+			for (k = 1; k <= nk; ++k) {
+				temp = scale1 * kwave[k - 1];
+				daughter[k - 1].re = norm*pow(temp,(double)m)*exp(-0.50*pow(temp,2.0));
+				daughter[k - 1].im = 0.0;
+			}
+		}
+		else if (re == 0) {
+			for (k = 1; k <= nk; ++k) {
+				temp = scale1 * kwave[k - 1];
+				daughter[k - 1].re = 0.0;
+				daughter[k - 1].im = norm*pow(temp, (double)m)*exp(-0.50*pow(temp, 2.0));
+			}
+		}
+		fourier_factor = (2.0*pi) * sqrt(2.0 / (2.0 * m + 1.0));
+		*period1 = scale1*fourier_factor;
+		*coi1 = fourier_factor / sqrt(2.0);
+	}
+}
+
+void cwavelet(const double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad,
+	double *wave, double *scale, double *period, double *coi) {
+
+	int i, j, k, iter;
+	double ymean, freq1, pi, period1, coi1;
+	double tmp1, tmp2;
+	double scale1;
+	double *kwave;
+	fft_object obj, iobj;
+	fft_data *ypad, *yfft,*daughter;
+
+	(void)s0; (void)dj; /* yes, we need these parameters unused */
+
+	pi = 4.0 * atan(1.0);
+
+	if (npad < N) {
+		printf("npad must be >= N \n");
+		exit(-1);
+	}
+
+	obj = fft_init(npad, 1);
+	iobj = fft_init(npad, -1);
+
+	ypad = (fft_data*)malloc(sizeof(fft_data)* npad);
+	yfft = (fft_data*)malloc(sizeof(fft_data)* npad);
+	daughter = (fft_data*)malloc(sizeof(fft_data)* npad);
+	kwave = (double*)malloc(sizeof(double)* npad);
+
+	ymean = 0.0;
+
+	for (i = 0; i < N; ++i) {
+		ymean += y[i];
+	}
+
+	ymean /= N;
+
+	for (i = 0; i < N; ++i) {
+		ypad[i].re = y[i] - ymean;
+		ypad[i].im = 0.0;
+	}
+
+	for (i = N; i < npad; ++i) {
+		ypad[i].re = ypad[i].im = 0.0;
+	}
+
+
+	// Find FFT of the input y (ypad)
+
+	fft_exec(obj, ypad, yfft);
+
+	for (i = 0; i < npad; ++i) {
+		yfft[i].re /= (double) npad;
+		yfft[i].im /= (double) npad;
+	}
+
+
+	//Construct the wavenumber array
+
+	freq1 = 2.0*pi / ((double)npad*dt);
+	kwave[0] = 0.0;
+
+	for (i = 1; i < npad / 2 + 1; ++i) {
+		kwave[i] = i * freq1;
+	}
+
+	for (i = npad / 2 + 1; i < npad; ++i) {
+		kwave[i] = -kwave[npad - i ];
+	}
+
+	
+	// Main loop
+
+	for (j = 1; j <= jtot; ++j) {
+		scale1 = scale[j - 1];// = s0*pow(2.0, (double)(j - 1)*dj);
+		wave_function(npad, dt, mother, param, scale1, kwave, pi,&period1,&coi1, daughter);
+		period[j - 1] = period1;
+		for (k = 0; k < npad; ++k) {
+			tmp1 = daughter[k].re * yfft[k].re - daughter[k].im * yfft[k].im;
+			tmp2 = daughter[k].re * yfft[k].im + daughter[k].im * yfft[k].re;
+			daughter[k].re = tmp1;
+			daughter[k].im = tmp2;
+		}
+		fft_exec(iobj, daughter, ypad);
+		iter = 2 * (j - 1) * N;
+		for (i = 0; i < N; ++i) {
+			wave[iter + 2 * i] = ypad[i].re;
+			wave[iter + 2 * i + 1] = ypad[i].im;
+		}
+	}
+	
+
+	for (i = 1; i <= (N + 1) / 2; ++i) {
+		coi[i - 1] = coi1 * dt * ((double)i - 1.0);
+		coi[N - i] = coi[i - 1];
+	}
+	
+
+	
+	free(kwave);
+	free(ypad);
+	free(yfft);
+	free(daughter);
+
+	free_fft(obj);
+	free_fft(iobj);
+
+}
+
+void psi0(int mother, double param,double *val,int *real) {
+	double pi,coeff;
+	int m,sign;
+
+	m = (int)param;
+	pi = 4.0 * atan(1.0);
+
+	if (mother == 0) {
+		// Morlet
+		*val = 1.0 / sqrt(sqrt(pi));
+		*real = 1;
+	}
+	else if (mother == 1) {
+		//Paul
+		if (m % 2 == 0) {
+			*real = 1;
+		}
+		else {
+			*real = 0;
+		}
+
+		if (m % 4 == 0 || m % 4 == 1) {
+			sign = 1;
+		}
+		else {
+			sign = -1;
+		}
+		*val = sign * pow(2.0, (double)m) * factorial(m) / (sqrt(pi * factorial(2 * m)));
+
+	}
+	else if (mother == 2) {
+		// D.O.G
+		*real = 1;
+
+		if (m % 2 == 0) {
+			if (m % 4 == 0) {
+				sign = -1;
+			}
+			else {
+				sign = 1;
+			}
+			coeff = sign * pow(2.0, (double)m / 2) / cwt_gamma(0.5);
+			*val = coeff * cwt_gamma(((double)m + 1.0) / 2.0) / sqrt(cwt_gamma(m + 0.50));
+		}
+		else {
+			*val = 0;
+		}
+	}
+}
+
+static int maxabs(double *array,int N) {
+	double maxval,temp;
+	int i,index;
+	maxval = 0.0;
+	index = -1;
+
+	for (i = 0; i < N; ++i) {
+		temp = fabs(array[i]);
+		if (temp >= maxval) {
+			maxval = temp;
+			index = i;
+		}
+	}
+
+	return index;
+}
+
+
+double cdelta(int mother, double param, double psi0 ) {
+	int N,i,j,iter;
+	double *delta, *scale,*period,*wave,*coi,*mval;
+	double den,cdel;
+	double subscale,dt,dj,s0;
+	int jtot;
+	int maxarr;
+
+	subscale = 8.0;
+	dt = 0.25;
+	if (mother == 0) {
+		N = 16;
+		s0 = dt/4;
+	}
+	else if (mother == 1) {
+		N = 16;
+		s0 = dt / 4.0;
+	}
+	else if (mother == 2)
+	{
+		s0 = dt/8.0;
+		N = 256;
+		if (param == 2.0) {
+			subscale = 16.0;
+			s0 = dt / 16.0;
+			N = 2048;
+		}
+	}
+
+	dj = 1.0 / subscale;
+	jtot = 16 * (int) subscale;
+
+	delta = (double*)malloc(sizeof(double)* N);
+	wave = (double*)malloc(sizeof(double)* 2 * N * jtot);
+	coi = (double*)malloc(sizeof(double)* N);
+	scale = (double*)malloc(sizeof(double)* jtot);
+	period = (double*)malloc(sizeof(double)* jtot);
+	mval = (double*)malloc(sizeof(double)* N);
+
+
+	delta[0] = 1;
+
+	for (i = 1; i < N; ++i) {
+		delta[i] = 0;
+	}
+
+	for (i = 0; i < jtot; ++i) {
+		scale[i] = s0*pow(2.0, (double)(i)*dj);
+	}
+
+	cwavelet(delta, N, dt, mother, param, s0, dj, jtot, N, wave, scale, period, coi);
+
+	for (i = 0; i < N; ++i) {
+		mval[i] = 0;
+	}
+
+	for (j = 0; j < jtot; ++j) {
+		iter = 2 * j * N;
+		den = sqrt(scale[j]);
+		for (i = 0; i < N; ++i) {
+			mval[i] += wave[iter + 2 * i]/den;
+		}
+	}
+
+
+	maxarr = maxabs(mval, N);
+
+	cdel = sqrt(dt) * dj * mval[maxarr] / psi0;
+
+	free(delta);
+	free(wave);
+
+	free(scale);
+	free(period);
+	free(coi);
+	free(mval);
+
+	return cdel;
+}
+
+void icwavelet(double *wave, int N, double *scale,int jtot,double dt,double dj,double cdelta,double psi0,double *oup) {
+	int i, j,iter;
+	double den, coeff;
+
+	coeff = sqrt(dt) * dj / (cdelta *psi0);
+
+	for (i = 0; i < N; ++i) {
+		oup[i] = 0.0;
+	}
+
+	for (j = 0; j < jtot; ++j) {
+		iter = 2 * j * N;
+		den = sqrt(scale[j]);
+		for (i = 0; i < N; ++i) {
+			oup[i] += wave[iter + 2 * i] / den;
+		}
+	}
+
+	for (i = 0; i < N; ++i) {
+		oup[i] *= coeff;
+	}
+}

src/cwt.h

@@ -0,0 +1,29 @@
+#ifndef CWT_H_
+#define CWT_H_
+
+#include "wavefunc.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void cwavelet(const double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad,
+		double *wave, double *scale, double *period, double *coi);
+
+void psi0(int mother, double param, double *val, int *real);
+
+double factorial(int N);
+
+double cdelta(int mother, double param, double psi0);
+
+void icwavelet(double *wave, int N, double *scale, int jtot, double dt, double dj, double cdelta, double psi0, double *oup);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */
+
+

src/cwtmath.c

@@ -0,0 +1,296 @@
+#include "cwtmath.h"
+
+static void nsfft_fd(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *w) {
+	int M,N,i,j,L;
+	double delta,den,theta,tempr,tempi,plb;
+	double *temp1,*temp2;
+
+	N = obj->N;
+	L = N/2;
+	//w = (double*)malloc(sizeof(double)*N);
+	
+	M = divideby(N, 2);
+	
+	if (M == 0) {
+		printf("The Non-Standard FFT Length must be a power of 2");
+		exit(1);
+	}
+	
+	temp1 = (double*)malloc(sizeof(double)*L);
+	temp2 = (double*)malloc(sizeof(double)*L);
+	
+	delta = (ub - lb)/ N;
+	j = -N;
+	den = 2 * (ub-lb);
+	
+	for(i = 0; i < N;++i) {
+		w[i] = (double)j/den;
+		j += 2;
+	}
+	
+	fft_exec(obj,inp,oup);
+	
+
+	for (i = 0; i < L; ++i) {
+		temp1[i] = oup[i].re;
+		temp2[i] = oup[i].im;
+	}
+
+	for (i = 0; i < N - L; ++i) {
+		oup[i].re = oup[i + L].re;
+		oup[i].im = oup[i + L].im;
+	}
+
+	for (i = 0; i < L; ++i) {
+		oup[N - L + i].re = temp1[i];
+		oup[N - L + i].im = temp2[i];
+	}
+	
+	plb = PI2 * lb;
+	
+	for(i = 0; i < N;++i) {
+		tempr = oup[i].re;
+		tempi = oup[i].im;
+		theta = w[i] * plb;
+		
+		oup[i].re = delta * (tempr*cos(theta) + tempi*sin(theta));
+		oup[i].im = delta * (tempi*cos(theta) - tempr*sin(theta));
+	}
+
+	
+	//free(w);
+	free(temp1);
+	free(temp2);
+}
+
+static void nsfft_bk(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *t) {
+	int M,N,i,j,L;
+	double *w;
+	double delta,den,plb,theta;
+	double *temp1,*temp2;
+	fft_data *inpt;
+
+	N = obj->N;
+	L = N/2;
+	
+	M = divideby(N, 2);
+	
+	if (M == 0) {
+		printf("The Non-Standard FFT Length must be a power of 2");
+		exit(1);
+	}
+	
+	temp1 = (double*)malloc(sizeof(double)*L);
+	temp2 = (double*)malloc(sizeof(double)*L);
+	w = (double*)malloc(sizeof(double)*N);
+	inpt = (fft_data*) malloc (sizeof(fft_data) * N);
+	
+	delta = (ub - lb)/ N;
+	j = -N;
+	den = 2 * (ub-lb);
+	
+	for(i = 0; i < N;++i) {
+		w[i] = (double)j/den;
+		j += 2;
+	}
+	
+	plb = PI2 * lb;
+	
+	for(i = 0; i < N;++i) {
+		theta = w[i] * plb;
+		
+		inpt[i].re = (inp[i].re*cos(theta) - inp[i].im*sin(theta))/delta;
+		inpt[i].im = (inp[i].im*cos(theta) + inp[i].re*sin(theta))/delta;
+	}
+	
+	for (i = 0; i < L; ++i) {
+		temp1[i] = inpt[i].re;
+		temp2[i] = inpt[i].im;
+	}
+
+	for (i = 0; i < N - L; ++i) {
+		inpt[i].re = inpt[i + L].re;
+		inpt[i].im = inpt[i + L].im;
+	}
+
+	for (i = 0; i < L; ++i) {
+		inpt[N - L + i].re = temp1[i];
+		inpt[N - L + i].im = temp2[i];
+	}
+	
+	fft_exec(obj,inpt,oup);
+	
+	for(i = 0; i < N;++i) {
+		t[i] = lb + i*delta;
+	}
+	
+	free(w);
+	free(temp1);
+	free(temp2);
+	free(inpt);
+}
+
+void nsfft_exec(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *w) {
+	if (obj->sgn == 1) {
+		nsfft_fd(obj,inp,oup,lb,ub,w);
+	} else if (obj->sgn == -1) {
+		nsfft_bk(obj,inp,oup,lb,ub,w);
+	}
+}
+
+static double fix(double x) {
+	// Rounds to the integer nearest to zero 
+	if (x >= 0.) {
+		return floor(x);
+	} else {
+		return ceil(x);
+	}
+}
+
+int nint(double N) {
+	int i;
+
+	i = (int)(N + 0.49999);
+
+	return i;
+}
+
+double cwt_gamma(double x) {
+	/*
+	 * This C program code is based on  W J Cody's fortran code.
+	 * http://www.netlib.org/specfun/gamma
+	 * 
+	 * References:
+   "An Overview of Software Development for Special Functions",
+	W. J. Cody, Lecture Notes in Mathematics, 506,
+	Numerical Analysis Dundee, 1975, G. A. Watson (ed.),
+	Springer Verlag, Berlin, 1976.
+
+   Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
+   */ 
+   
+	// numerator and denominator coefficients for 1 <= x <= 2
+	
+	double y,oup,fact,sum,y2,yi,z,nsum,dsum;
+	int swi,n,i;
+	
+	double spi = 0.9189385332046727417803297;
+    double pi  = 3.1415926535897932384626434;
+	double xmax = 171.624e+0;
+	double xinf = 1.79e308;
+	double eps = 2.22e-16;
+	double xninf = 1.79e-308;
+	
+	double num[8] = { -1.71618513886549492533811e+0,
+                        2.47656508055759199108314e+1,
+                       -3.79804256470945635097577e+2,
+                        6.29331155312818442661052e+2,
+                        8.66966202790413211295064e+2,
+                       -3.14512729688483675254357e+4,
+                       -3.61444134186911729807069e+4,
+                        6.64561438202405440627855e+4 };
+                       
+    double den[8] = { -3.08402300119738975254353e+1,
+                        3.15350626979604161529144e+2,
+                       -1.01515636749021914166146e+3,
+                       -3.10777167157231109440444e+3,
+                        2.25381184209801510330112e+4,
+                        4.75584627752788110767815e+3,
+                       -1.34659959864969306392456e+5,
+                       -1.15132259675553483497211e+5 }; 
+                       
+    // Coefficients for Hart's Minimax approximation x >= 12       
+    
+    
+	double c[7] = { -1.910444077728e-03,
+                        8.4171387781295e-04,
+                       -5.952379913043012e-04,
+                        7.93650793500350248e-04,
+                       -2.777777777777681622553e-03,
+                        8.333333333333333331554247e-02,
+                        5.7083835261e-03 };            
+    
+    y = x;             
+    swi = 0; 
+    fact = 1.0;
+    n = 0;      
+    
+    
+    if ( y < 0.) {
+		// Negative x
+		y = -x;
+		yi = fix(y);
+		oup = y - yi;
+		
+		if (oup != 0.0) {
+			if (yi != fix(yi * .5) * 2.) {
+				swi = 1;
+			}
+			fact = -pi / sin(pi * oup);
+			y += 1.;
+		} else {
+			return xinf;
+		}
+	}      
+	
+	if (y < eps) {
+		if (y >= xninf) {
+			oup = 1.0/y;
+		} else {
+			return xinf;
+		}
+		
+	} else if (y < 12.) {
+		yi = y;
+		if ( y < 1.) {
+			z = y;
+			y += 1.;
+		} else {
+			n = ( int ) y - 1;
+            y -= ( double ) n;
+            z = y - 1.0;
+		}
+		nsum = 0.;
+		dsum = 1.;
+		for (i = 0; i < 8; ++i) {
+			nsum = (nsum + num[i]) * z;
+			dsum = dsum * z + den[i];
+		}
+		oup = nsum / dsum + 1.;
+		
+		if (yi < y) {
+	   
+			oup /= yi;
+		} else if (yi > y) {
+	    
+			for (i = 0; i < n; ++i) {
+				oup *= y;
+				y += 1.;
+			}
+		
+		}      
+	
+	} else {
+		if (y <= xmax) {
+			y2 = y * y;
+			sum = c[6];
+			for (i = 0; i < 6; ++i) {
+				sum = sum / y2 + c[i];
+			}
+			sum = sum / y - y + spi;
+			sum += (y - .5) * log(y);
+			oup = exp(sum);
+		} else {
+			return(xinf);
+		}
+	}
+	
+	if (swi) {
+		oup = -oup;
+	}
+    if (fact != 1.) {
+		oup = fact / oup; 
+	}       
+	
+	return oup;
+}

src/cwtmath.h

@@ -0,0 +1,22 @@
+#ifndef CWTMATH_H_
+#define CWTMATH_H_
+
+#include "wtmath.h"
+#include "hsfft.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void nsfft_exec(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *w);// lb -lower bound, ub - upper bound, w - time or frequency grid (Size N)
+
+double cwt_gamma(double x);
+
+int nint(double N);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */

src/hsfft.c

@@ -18,7 +18,7 @@ fft_object fft_init(int N, int sgn) {
 	if (out == 1) {
 		obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (N-1));
 		obj->lf = factors(N,obj->factors);
-		longvectorN(obj->twiddle,N,obj->factors,obj->lf);
+		longvectorN(obj->twiddle,obj->factors,obj->lf);
 		twi_len = N;
 		obj->lt = 0;
 	} else {
@@ -32,7 +32,7 @@ fft_object fft_init(int N, int sgn) {
 		}
 		obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (M-1));
 		obj->lf = factors(M,obj->factors);
-		longvectorN(obj->twiddle,M,obj->factors,obj->lf);
+		longvectorN(obj->twiddle,obj->factors,obj->lf);
 		obj->lt = 1;
 		twi_len = M;
 	}
@@ -1831,7 +1831,7 @@ void twiddle(fft_data *vec,int N, int radix) {
 
 }
 
-void longvectorN(fft_data *sig,int N, int *array, int tx) {
+void longvectorN(fft_data *sig, int *array, int tx) {
 	int L,i,Ls,ct,j,k;
 	fft_type theta;
 	L = 1;

src/hsfft.h

@@ -12,6 +12,8 @@
 #include <math.h>
 #include <string.h>
 
+#include "../header/wavelib.h"
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -22,11 +24,6 @@ extern "C" {
 #define fft_type double
 #endif
 
-
-typedef struct fft_t {
-  fft_type re;
-  fft_type im;
-} fft_data;
 /*
 #define SADD(a,b) ((a)+(b))
 
@@ -35,19 +32,8 @@ typedef struct fft_t {
 #define SMUL(a,b) ((a)*(b))
 */
 
-typedef struct fft_set* fft_object;
-
 fft_object fft_init(int N, int sgn);
 
-struct fft_set{
-	int N;
-	int sgn;
-	int factors[64];
-	int lf;
-	int lt;
-	fft_data twiddle[1];
-};
-
 void fft_exec(fft_object obj,fft_data *inp,fft_data *oup);
 
 int divideby(int M,int d);
@@ -60,7 +46,7 @@ int factors(int M, int* arr);
 
 void twiddle(fft_data *sig,int N, int radix);
 
-void longvectorN(fft_data *sig,int N, int *array, int M);
+void longvectorN(fft_data *sig, int *array, int M);
 
 void free_fft(fft_object object);
 

src/real.c

@@ -9,11 +9,9 @@
 
 fft_real_object fft_real_init(int N, int sgn) {
 	fft_real_object obj = NULL;
-	fft_type PI, theta;
+	fft_type theta;
 	int k;
 	
-	PI = 3.1415926535897932384626433832795;
-	
 	obj = (fft_real_object) malloc (sizeof(struct fft_real_set) + sizeof(fft_data)* (N/2));
 	
 	obj->cobj = fft_init(N/2,sgn);
@@ -79,10 +77,9 @@ void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup) {
 	
 	fft_data* cinp;
 	fft_data* coup;
-	int i,N2,N;
+	int i,N2;
 	fft_type temp1,temp2;
 	N2 = obj->cobj->N;
-	N = N2*2;
 	
 	cinp = (fft_data*) malloc (sizeof(fft_data) * N2);
 	coup = (fft_data*) malloc (sizeof(fft_data) * N2);
@@ -106,6 +103,7 @@ void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup) {
 }
 
 void free_real_fft(fft_real_object object) {
+	free_fft(object->cobj);
 	free(object);
 }
 

src/real.h

@@ -14,15 +14,8 @@
 extern "C" {
 #endif
 
-typedef struct fft_real_set* fft_real_object;
-
 fft_real_object fft_real_init(int N, int sgn);
 
-struct fft_real_set{
-	fft_object cobj;
-	fft_data twiddle2[1];
-};
-
 void fft_r2c_exec(fft_real_object obj,fft_type *inp,fft_data *oup);
 
 void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup);

src/wavefilt.h

@@ -1,19 +1,28 @@
+/*
+Copyright (c) 2014, Rafat Hussain
+Copyright (c) 2016, Holger Nahrstaedt
+*/
 #ifndef WAVEFILT_H_
 #define WAVEFILT_H_
 
 #include <stdio.h>
 #include "conv.h"
+#define _USE_MATH_DEFINES
+#include "math.h"
 
 #ifdef __cplusplus
 extern "C" {
 #endif
 
 
-int filtlength(char* name);
-
-int filtcoef(char* name, double *lp1, double *hp1, double *lp2, double *hp2);
+int filtlength(const char* name);
 
+int filtcoef(const char* name, double *lp1, double *hp1, double *lp2, double *hp2);
 
+void copy_reverse(const double *in, int N, double *out);
+void qmf_even(const double *in, int N, double *out);
+void qmf_wrev(const double *in, int N, double *out);
+void copy(const double *in, int N, double *out);
 #ifdef __cplusplus
 }
 #endif

src/wavefunc.c

@@ -0,0 +1,210 @@
+#include "wavefunc.h"
+
+void meyer(int N,double lb,double ub,double *phi,double *psi,double *tgrid) {
+	int M,i;
+	double *w;
+	double delta,j;
+	double theta,x,x2,x3,x4,v,cs,sn;
+	double wf;
+	fft_data *phiw,*psiw,*oup;
+	fft_object obj;
+	
+	M = divideby(N, 2);
+	
+	if (M == 0) {
+		printf("Size of Wavelet must be a power of 2");
+		exit(1);
+	}
+	if (lb >= ub) {
+		printf("upper bound must be greater than lower bound");
+		exit(1);
+	}
+	
+	obj = fft_init(N,-1);
+	w = (double*)malloc(sizeof(double)*N);
+	phiw = (fft_data*) malloc (sizeof(fft_data) * N);
+	psiw = (fft_data*) malloc (sizeof(fft_data) * N);
+	oup = (fft_data*) malloc (sizeof(fft_data) * N);
+	
+	delta = 2 * (ub-lb) / PI2;
+	
+	j = (double) N;
+	j *= -1.0;
+	
+	for(i = 0; i < N;++i) {
+		w[i] = j / delta;
+		j += 2.0; 
+		psiw[i].re = psiw[i].im = 0.0;
+		phiw[i].re = phiw[i].im = 0.0;
+	}
+	
+	
+	for(i = 0; i < N;++i) {
+		wf = fabs(w[i]);
+		if (wf <= PI2/3.0) {
+			phiw[i].re = 1.0;
+		}
+		if (wf > PI2/3.0 && wf <= 2 * PI2 / 3.0) {
+			x = (3 * wf / PI2) - 1.0;
+			x2 = x*x;
+			x3 = x2 * x;
+			x4 = x3 *x;
+			v = x4 *(35 - 84*x + 70*x2 - 20*x3);
+			theta = v * PI2 / 4.0;
+			cs = cos(theta);
+			sn = sin(theta);
+			
+			phiw[i].re = cs;
+			psiw[i].re = cos(w[i]/2.0) * sn;
+			psiw[i].im = sin(w[i]/2.0) * sn;
+		}
+		if (wf > 2.0 * PI2/3.0 && wf <= 4 * PI2 / 3.0) {
+			x = (1.5 * wf / PI2) - 1.0;
+			x2 = x*x;
+			x3 = x2 * x;
+			x4 = x3 *x;
+			v = x4 *(35 - 84*x + 70*x2 - 20*x3);
+			theta = v * PI2 / 4.0;
+			cs = cos(theta);
+			
+			psiw[i].re = cos(w[i]/2.0) * cs;
+			psiw[i].im = sin(w[i]/2.0) * cs;
+		}
+	}
+	
+	
+	nsfft_exec(obj,phiw,oup,lb,ub,tgrid);
+	
+	
+	for(i = 0; i < N;++i) {
+		phi[i] = oup[i].re/N;
+	}
+	
+	nsfft_exec(obj,psiw,oup,lb,ub,tgrid);
+	
+	
+	for(i = 0; i < N;++i) {
+		psi[i] = oup[i].re/N;
+	}
+	
+	
+	
+	free(oup);
+	free(phiw);
+	free(psiw);
+	free(w);
+}
+
+void gauss(int N,int p,double lb,double ub,double *psi,double *t) {
+  double delta,num,den,t2,t4;
+	int i;
+	
+	if (lb >= ub) {
+		printf("upper bound must be greater than lower bound");
+		exit(1);
+	}
+	
+	t[0] = lb;
+	t[N-1] = ub;
+	delta = (ub - lb) / (N-1);
+	for(i = 1; i < N-1;++i) {
+		t[i] = lb + delta * i;
+	}
+	
+	den = sqrt(cwt_gamma(p+0.5));
+	
+	if ((p+1)%2 == 0) {
+		num = 1.0;
+	} else {
+		num = -1.0;
+	}
+	
+	num /= den;
+	
+	//printf("\n%g\n",num);
+	
+	if (p == 1) {
+		for(i = 0; i < N;++i) {
+			psi[i] = -t[i] * exp(- t[i] * t[i]/2.0) * num; 
+		}
+	} else if (p == 2) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = (-1.0 + t2) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 3) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = t[i] * (3.0 - t2) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 4) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = (t2 * t2 - 6.0 * t2 + 3.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 5) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = t[i] * (-t2 * t2 + 10.0 * t2 - 15.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 6) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = (t2 * t2 * t2 - 15.0 * t2 * t2 + 45.0 * t2 - 15.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 7) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = t[i] * (-t2 * t2 * t2 + 21.0 * t2 * t2 - 105.0 * t2 + 105.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 8) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			t4 = t2 * t2;
+			psi[i] = (t4 * t4 - 28.0 * t4 * t2 + 210.0 * t4 - 420.0 * t2 + 105.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 9) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			t4 = t2 * t2;
+			psi[i] = t[i] * (- t4 * t4 + 36.0 * t4 * t2 - 378.0 * t4 + 1260.0 * t2 - 945.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 10) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			t4 = t2 * t2;
+			psi[i] = (t4 * t4 * t2 - 45.0 * t4 * t4 + 630.0 * t4 * t2 - 3150.0 * t4 + 4725.0 * t2 - 945.0) * exp(- t2/2.0) * num; 
+		}
+	} else {
+	  printf("\n The Gaussian Derivative Wavelet is only available for Derivatives 1 to 10");
+	  exit(1);
+	}
+
+	
+	
+}
+
+void mexhat(int N,double lb,double ub,double *psi,double *t) {
+  gauss(N,2,lb,ub,psi,t);
+}
+
+void morlet(int N,double lb,double ub,double *psi,double *t) {
+       int i;
+       double delta;
+
+  	if (lb >= ub) {
+		printf("upper bound must be greater than lower bound");
+		exit(1);
+	}
+	
+	t[0] = lb;
+	t[N-1] = ub;
+	delta = (ub - lb) / (N-1);
+	for(i = 1; i < N-1;++i) {
+		t[i] = lb + delta * i;
+	}
+
+	for(i = 0; i < N;++i) {
+	  psi[i] = exp(- t[i] * t[i] / 2.0) * cos(5 * t[i]);
+	}
+}

src/wavefunc.h

@@ -0,0 +1,24 @@
+#ifndef WAVEFUNC_H_
+#define WAVEFUNC_H_
+
+#include "cwtmath.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void meyer(int N,double lb,double ub,double *phi,double *psi,double *tgrid);
+
+void gauss(int N,int p,double lb,double ub,double *psi,double *t);
+
+void mexhat(int N,double lb,double ub,double *psi,double *t);
+
+void morlet(int N,double lb,double ub,double *psi,double *t);  
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVEFUNC_H_ */

src/wavelib.h

@@ -1,85 +0,0 @@
-#ifndef WAVELIB_H_
-#define WAVELIB_H_
-
-#include "wtmath.h"
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-typedef struct wave_set* wave_object;
-
-wave_object wave_init(char* wname);
-
-struct wave_set{
-	char wname[50];
-	int filtlength;// When all filters are of the same length. [Matlab uses zero-padding to make all filters of the same length]
-	int lpd_len;// Default filtlength = lpd_len = lpr_len = hpd_len = hpr_len
-	int hpd_len;
-	int lpr_len;
-	int hpr_len;
-	double *lpd;
-	double *hpd;
-	double *lpr;
-	double *hpr;
-	double params[0];
-};
-
-typedef struct wt_set* wt_object;
-
-wt_object wt_init(wave_object wave,char* method, int siglength, int J);
-
-struct wt_set{
-	wave_object wave;
-	conv_object cobj;
-	char method[10];
-	int siglength;// Length of the original signal.
-	int outlength;// Length of the output DWT vector
-	int lenlength;// Length of the Output Dimension Vector "length"
-	int J; // Number of decomposition Levels
-	int MaxIter;// Maximum Iterations J <= MaxIter
-	int even;// even = 1 if signal is of even length. even = 0 otherwise
-	char ext[10];// Type of Extension used - "per" or "sym"
-	char cmethod[10]; // Convolution Method - "direct" or "FFT"
-
-	int N; // 
-	int cfftset; 
-	int zpad;
-	int length[102];
-	double *output;
-	double params[0];
-};
-
-void dwt(wt_object wt, double *inp);
-
-void idwt(wt_object wt, double *dwtop);
-
-void swt(wt_object wt, double *inp);
-
-void iswt(wt_object wt, double *swtop);
-
-void modwt(wt_object wt, double *inp);
-
-void imodwt(wt_object wt, double *dwtop);
-
-void setDWTExtension(wt_object wt, char *extension);
-
-void setWTConv(wt_object wt, char *cmethod);
-
-void wave_summary(wave_object obj);
-
-void wt_summary(wt_object wt);
-
-void wave_free(wave_object object);
-
-void wt_free(wt_object object);
-
-
-#ifdef __cplusplus
-}
-#endif
-
-
-#endif /* WAVELIB_H_ */
-
-

src/wtmath.c

@@ -1,5 +1,355 @@
+/*
+Copyright (c) 2018, Rafat Hussain
+*/
 #include "wtmath.h"
 
+void dwt_per_stride(double *inp, int N, double *lpd,double*hpd,int lpd_len,double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int l, l2, isodd, i, t, len_avg,is,os;
+
+	len_avg = lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + l2;
+		os = i *ostride;
+		cA[os] = 0.0;
+		cD[os] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= l2 && (t - l) < N) {
+				is = (t - l) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < l2 && (t - l) >= 0) {
+				is = (t - l) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < 0 && isodd == 0) {
+				is = (t - l + N) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < 0 && isodd == 1) {
+				if ((t - l) != -1) {
+					is = (t - l + N + 1) * istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+				else {
+					is = (N - 1) * istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+			}
+			else if ((t - l) >= N && isodd == 0) {
+				is = (t - l - N) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) >= N && isodd == 1) {
+				is = (t - l - (N + 1)) * istride;
+				if (t - l != N) {
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+				else {
+					is = (N - 1) * istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+			}
+
+		}
+		
+	}
+
+}
+
+void dwt_sym_stride(double *inp, int N, double *lpd, double*hpd, int lpd_len, double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int i, l, t, len_avg;
+	int is, os;
+	len_avg = lpd_len;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + 1;
+		os = i *ostride;
+		cA[os] = 0.0;
+		cD[os] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= 0 && (t - l) < N) {
+				is = (t - l) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < 0) {
+				is = (-t + l - 1) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) >= N) {
+				is = (2 * N - t + l - 1) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+
+		}
+	}
+
+
+}
+
+void modwt_per_stride(int M, double *inp, int N, double *filt, int lpd_len, double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int l, i, t, len_avg;
+	int is, os;
+	len_avg = lpd_len;
+
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i;
+		os = i *ostride;
+		is = t *istride;
+		cA[os] = filt[0] * inp[is];
+		cD[os] = filt[len_avg] * inp[is];
+		for (l = 1; l < len_avg; l++) {
+			t -= M;
+			while (t >= len_cA) {
+				t -= len_cA;
+			}
+			while (t < 0) {
+				t += len_cA;
+			}
+			os = i * ostride;
+			is = t * istride;
+			cA[os] += filt[l] * inp[is];
+			cD[os] += filt[len_avg + l] * inp[is];
+
+		}
+	}
+	
+}
+
+void swt_per_stride(int M, double *inp, int N, double *lpd, double*hpd, int lpd_len, double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int l, l2, isodd, i, t, len_avg, j;
+	int is, os;
+	len_avg = M * lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i + l2;
+		os = i *ostride;
+		cA[os] = 0.0;
+		cD[os] = 0.0;
+		l = -1;
+		for (j = 0; j < len_avg; j += M) {
+			l++;
+			while (j >= len_cA) {
+				j -= len_cA;
+			}
+			if ((t - j) >= l2 && (t - j) < N) {
+				is = (t - j)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) < l2 && (t - j) >= 0) {
+				is = (t - j)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) < 0) {
+				is = (t - j + N)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) >= N && isodd == 0) {
+				is = (t - j - N)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) >= N && isodd == 1) {
+				if (t - l != N) {
+					is = (t - j - (N + 1))*istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+				else {
+					is = (N - 1)*istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[N - 1];
+				}
+			}
+
+		}
+	}
+
+}
+
+void idwt_per_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, int lpr_len, double *X, int istride, int ostride) {
+	int len_avg, i, l, m, n, t, l2;
+	int is, ms, ns;
+
+	len_avg = lpr_len;
+	l2 = len_avg / 2;
+	m = -2;
+	n = -1;
+
+	for (i = 0; i < len_cA + l2 - 1; ++i) {
+		m += 2;
+		n += 2;
+		ms = m * ostride;
+		ns = n * ostride;
+		X[ms] = 0.0;
+		X[ns] = 0.0;
+		for (l = 0; l < l2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				is = (i - l) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+			else if ((i - l) >= len_cA && (i - l) < len_cA + len_avg - 1) {
+				is = (i - l - len_cA) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+			else if ((i - l) < 0 && (i - l) > -l2) {
+				is = (len_cA + i - l) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+		}
+	}
+}
+
+void idwt_sym_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, int lpr_len, double *X, int istride, int ostride) {
+	int len_avg, i, l, m, n, t, v;
+	int ms, ns, is;
+	len_avg = lpr_len;
+	m = -2;
+	n = -1;
+
+	for (v = 0; v < len_cA; ++v) {
+		i = v;
+		m += 2;
+		n += 2;
+		ms = m * ostride;
+		ns = n * ostride;
+		X[ms] = 0.0;
+		X[ns] = 0.0;
+		for (l = 0; l < len_avg / 2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				is = (i - l) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+		}
+	}
+}
+
+void imodwt_per_stride(int M, double *cA, int len_cA, double *cD, double *filt,int lf,double *X,int istride, int ostride) {
+	int len_avg, i, l, t;
+	int is, os;
+	
+	len_avg = lf;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i;
+		os = i * ostride;
+		is = t *istride;
+		X[os] = (filt[0] * cA[is]) + (filt[len_avg] * cD[is]);
+		for (l = 1; l < len_avg; l++) {
+			t += M;
+			while (t >= len_cA) {
+				t -= len_cA;
+			}
+			while (t < 0) {
+				t += len_cA;
+			}
+			is = t *istride;
+			X[os] += (filt[l] * cA[is]) + (filt[len_avg + l] * cD[is]);
+
+		}
+	}
+
+}
+
+void idwt2_shift(int shift, int rows, int cols, double *lpr, double *hpr, int lf, double *A,double *H, double *V,double *D, double *oup) {
+	int i, k, N, ir, ic, J, dim1, dim2;
+	int istride, ostride;
+	double *cL, *cH, *X_lp;
+
+
+	N = rows > cols ? 2 * rows : 2 * cols;
+
+	J = 1;
+	dim1 = 2 * rows;
+	dim2 = 2 * cols;
+
+	X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+	cL = (double*)calloc(dim1*dim2, sizeof(double));
+	cH = (double*)calloc(dim1*dim2, sizeof(double));
+
+	ir = rows;
+	ic = cols;
+	istride = ic;
+	ostride = 1;
+	for (i = 0; i < ic; ++i) {
+		idwt_per_stride(A+i, ir, H+i, lpr, hpr, lf, X_lp, istride, ostride);
+
+		for (k = lf / 2 - 1; k < 2 * ir + lf / 2 - 1; ++k) {
+			cL[(k - lf / 2 + 1)*ic + i] = X_lp[k];
+		}
+
+		idwt_per_stride(V+i, ir, D+i, lpr, hpr, lf, X_lp, istride, ostride);
+
+		for (k = lf / 2 - 1; k < 2 * ir + lf / 2 - 1; ++k) {
+			cH[(k - lf / 2 + 1)*ic + i] = X_lp[k];
+		}
+	}
+	ir *= 2;
+	istride = 1;
+	ostride = 1;
+
+	for (i = 0; i < ir; ++i) {
+		idwt_per_stride(cL + i*ic, ic, cH + i*ic, lpr, hpr, lf, X_lp, istride, ostride);
+
+		for (k = lf / 2 - 1; k < 2 * ic + lf / 2 - 1; ++k) {
+			oup[(k - lf / 2 + 1) + i*ic * 2] = X_lp[k];
+		}
+	}
+	ic *= 2;
+
+
+	if (shift == -1) {
+		//Save the last column
+		for (i = 0; i < ir; ++i) {
+			cL[i] = oup[(i + 1)*ic - 1];
+		}
+		// Save the last row
+		memcpy(cH, oup + (ir - 1)*ic, sizeof(double)*ic);
+		for (i = ir - 1; i > 0; --i) {
+			memcpy(oup + i*ic + 1, oup + (i - 1)*ic, sizeof(double)*(ic - 1));
+		}
+		oup[0] = cL[ir - 1];
+		for (i = 1; i < ir; ++i) {
+			oup[i*ic] = cL[i - 1];
+		}
+
+		for (i = 1; i < ic; ++i) {
+			oup[i] = cH[i - 1];
+		}
+	}
+
+
+	free(X_lp);
+	free(cL);
+	free(cH);
+
+}
+
 int upsamp(double *x, int lenx, int M, double *y) {
 	int N, i, j, k;
 
@@ -236,4 +586,95 @@ int wmaxiter(int sig_len, int filt_len) {
 	lev = (int)temp;
 
 	return lev;
-}
+}
+
+static double entropy_s(double *x,int N) {
+  int i;
+  double val,x2;
+
+  val = 0.0;
+
+  for(i = 0; i < N; ++i) {
+    if (x[i] != 0) {
+      x2 = x[i] * x[i];
+      val -= x2 * log(x2);
+    }
+  }
+  return val;
+}
+
+static double entropy_t(double *x,int N, double t) {
+  int i;
+  double val,x2;
+  if (t < 0) {
+    printf("Threshold value must be >= 0");
+    exit(1);
+  }
+  val = 0.0;
+
+  for(i = 0; i < N; ++i) {
+    x2 = fabs(x[i]);
+    if (x2 > t) {
+      val += 1;
+    }
+
+  }
+
+  return val;
+
+}
+
+static double entropy_n(double *x,int N,double p) {
+  int i;
+  double val,x2;
+  if (p < 1) {
+    printf("Norm power value must be >= 1");
+    exit(1);
+  }
+  val = 0.0;
+  for(i = 0; i < N; ++i) {
+    x2 = fabs(x[i]);
+    val += pow(x2,(double)p);
+
+  }
+
+  return val;
+}
+
+static double entropy_l(double *x,int N) {
+  int i;
+  double val,x2;
+
+  val = 0.0;
+
+  for(i = 0; i < N; ++i) {
+    if (x[i] != 0) {
+      x2 = x[i] * x[i];
+      val += log(x2);
+    }
+  }
+  return val;
+}
+
+double costfunc(double *x, int N ,char *entropy,double p) {
+	double val;
+
+	if (!strcmp(entropy, "shannon")) {
+		val = entropy_s(x, N);
+	}
+	else if (!strcmp(entropy, "threshold")) {
+		val = entropy_t(x, N,p);
+	}
+	else if (!strcmp(entropy, "norm")) {
+		val = entropy_n(x, N,p);
+	}
+	else if (!strcmp(entropy, "logenergy") || !strcmp(entropy, "log energy") || !strcmp(entropy, "energy")) {
+		val = entropy_l(x, N);
+	}
+	else {
+		printf("Entropy must be one of shannon, threshold, norm or energy");
+		exit(-1);
+	}
+
+	return val;
+}

src/wtmath.h

@@ -1,3 +1,6 @@
+/*
+Copyright (c) 2014, Rafat Hussain
+*/
 #ifndef WTMATH_H_
 #define WTMATH_H_
 
@@ -7,6 +10,30 @@
 extern "C" {
 #endif
 
+void dwt_per_stride(double *inp, int N, double *lpd,double*hpd,int lpd_len,
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void dwt_sym_stride(double *inp, int N, double *lpd, double*hpd, int lpd_len,
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void modwt_per_stride(int M, double *inp, int N, double *filt, int lpd_len, 
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void swt_per_stride(int M, double *inp, int N, double *lpd, double*hpd, int lpd_len, 
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void idwt_per_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, 
+int lpr_len, double *X, int istride, int ostride);
+
+void idwt_sym_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, 
+int lpr_len, double *X, int istride, int ostride);
+
+void imodwt_per_stride(int M, double *cA, int len_cA, double *cD, double *filt,
+int lf,double *X,int istride, int ostride);
+
+void idwt2_shift(int shift, int rows, int cols, double *lpr, double *hpr, int lf, 
+double *A,double *H, double *V,double *D, double *oup);
+
 int upsamp(double *x, int lenx, int M, double *y);
 
 int upsamp2(double *x, int lenx, int M, double *y);
@@ -23,9 +50,11 @@ int testSWTlength(int N, int J);
 
 int wmaxiter(int sig_len, int filt_len);
 
+double costfunc(double *x, int N, char *entropy, double p);
+
 #ifdef __cplusplus
 }
 #endif
 
 
-#endif /* WAVELIB_H_ */
+#endif /* WAVELIB_H_ */

test/CMakeLists.txt

@@ -0,0 +1,48 @@
+add_executable(cwttest cwttest.c)
+
+target_link_libraries(cwttest wavelib)
+
+add_executable(dwttest dwttest.c)
+
+target_link_libraries(dwttest wavelib)
+
+add_executable(swttest swttest.c)
+
+target_link_libraries(swttest wavelib)
+
+add_executable(modwttest modwttest.c)
+
+target_link_libraries(modwttest wavelib)
+
+add_executable(dwpttest dwpttest.c)
+
+target_link_libraries(dwpttest wavelib)
+
+add_executable(wtreetest wtreetest.c)
+
+target_link_libraries(wtreetest wavelib)
+
+add_executable(denoisetest denoisetest.c)
+
+target_link_libraries(denoisetest wauxlib wavelib)
+
+add_executable(modwtdenoisetest modwtdenoisetest.c)
+
+target_link_libraries(modwtdenoisetest wauxlib wavelib)
+
+add_executable(dwt2test dwt2test.c)
+
+target_link_libraries(dwt2test wavelib)
+
+add_executable(swt2test swt2test.c)
+
+target_link_libraries(swt2test wavelib)
+
+add_executable(modwt2test modwt2test.c)
+
+target_link_libraries(modwt2test wavelib)
+
+set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest modwtdenoisetest dwt2test swt2test modwt2test
+        PROPERTIES
+        RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/test"
+        )

test/cwttest.c

@@ -0,0 +1,109 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+int main() {
+	int i, N, J,subscale,a0,iter,nd,k;
+	double *inp,*oup;
+	double dt, dj,s0, param,mn;
+	double td,tn,den, num, recon_mean, recon_var;
+	cwt_object wt;
+
+	FILE *ifp;
+	double temp[1200];
+
+	char *wave = "morlet";// Set Morlet wavelet. Other options "paul" and "dog"
+	char *type = "pow";
+
+	N = 504;
+	param = 6.0;
+	subscale = 4;
+	dt = 0.25;
+	s0 = dt;
+	dj = 1.0 / (double)subscale;
+	J = 11 * subscale; // Total Number of scales
+	a0 = 2;//power
+
+	ifp = fopen("sst_nino3.dat", "r");
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+
+	fclose(ifp);
+
+	wt = cwt_init(wave, param, N,dt, J);
+
+	inp = (double*)malloc(sizeof(double)* N);
+	oup = (double*)malloc(sizeof(double)* N);
+
+	for (i = 0; i < N; ++i) {
+		inp[i] = temp[i] ;
+	}
+
+	setCWTScales(wt, s0, dj, type, a0);
+
+	cwt(wt, inp);
+
+	printf("\n MEAN %g \n", wt->smean);
+	
+	mn = 0.0;
+
+	for (i = 0; i < N; ++i) {
+		mn += sqrt(wt->output[i].re * wt->output[i].re + wt->output[i].im * wt->output[i].im);
+	}
+	
+	cwt_summary(wt);
+
+	printf("\n abs mean %g \n", mn / N);
+	
+	printf("\n\n");
+	printf("Let CWT w = w(j, n/2 - 1) where n = %d\n\n", N);
+	nd = N/2 - 1;
+	
+	printf("%-15s%-15s%-15s%-15s \n","j","Scale","Period","ABS(w)^2");
+	for(k = 0; k < wt->J;++k) {
+		iter = nd + k * N;
+		printf("%-15d%-15lf%-15lf%-15lf \n",k,wt->scale[k],wt->period[k],
+		wt->output[iter].re * wt->output[iter].re + wt->output[iter].im * wt->output[iter].im);
+	}
+
+	icwt(wt, oup);
+
+	num = den = recon_var = recon_mean = 0.0;
+	printf("\n\n");
+	printf("Signal Reconstruction\n");
+	printf("%-15s%-15s%-15s \n","i","Input(i)","Output(i)");
+
+	for (i = N - 10; i < N; ++i) {
+		printf("%-15d%-15lf%-15lf \n", i,inp[i] , oup[i]);
+	}
+
+	for (i = 0; i < N; ++i) {
+		//printf("%g %g \n", oup[i] ,inp[i] - wt->smean);
+		td = inp[i] ;
+		tn = oup[i] - td;
+		num += (tn * tn);
+		den += (td * td);
+		recon_mean += oup[i];
+	}
+
+	recon_var = sqrt(num / N);
+	recon_mean /= N;
+
+	printf("\nRMS Error %g \n", sqrt(num) / sqrt(den));
+	printf("\nVariance %g \n", recon_var);
+	printf("\nMean %g \n", recon_mean);
+
+	free(inp);
+	free(oup);
+	cwt_free(wt);
+	return 0;
+}

test/denoisetest.c

@@ -0,0 +1,139 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "../header/wauxlib.h"
+
+static double rmse(int N,double *x,double *y) {
+	double rms;
+	int i;
+
+	rms = 0.0;
+
+	for(i = 0; i < N;++i) {
+		rms += (x[i] - y[i]) * (x[i] - y[i]);
+	}
+
+	rms = sqrt(rms/(double)N);
+
+	return rms;
+}
+
+static double corrcoef(int N,double *x,double *y) {
+	double cc,xm,ym,tx,ty,num,den1,den2;
+	int i;
+	xm = ym = 0.0;
+	for(i = 0; i < N;++i) {
+		xm += x[i];
+		ym += y[i];
+	}
+
+	xm = xm/N;
+	ym = ym / N;
+	num = den1 = den2 = 0.0;
+
+	for(i = 0; i < N;++i) {
+		tx = x[i] - xm;
+		ty = y[i] - ym;
+		num += (tx*ty);
+		den1 += (tx*tx);
+		den2 += (ty*ty);
+	}
+
+	cc = num / sqrt(den1*den2);
+
+	return cc;
+}
+
+int main() {
+	// gcc -Wall -I../header -L../Bin denoisetest.c -o denoise -lwauxlib -lwavelib -lm
+	double *sig,*inp,*oup;
+	int i,N,J;
+	FILE *ifp;
+
+	denoise_object obj;
+	double temp[2400];
+
+	char *wname = "db5";
+	char *method = "dwt";// Available - dwt, swt and modwt. modwt works only with modwtshrink. The other two methods work with
+	// visushrink and sureshrink
+	char *ext = "sym";// sym and per work with dwt. swt and modwt only use per extension when called through denoise.
+	// You can use sy extension if you directly call modwtshrink with cmethod set to fft. See modwtdenoisetest.c file 
+	char *thresh = "soft";// soft or hard
+	char *level = "all"; // noise estimation at "first" or "all" levels. modwt only has the option of "all" 
+
+	ifp = fopen("pieceregular1024.txt", "r");
+	i = 0;
+	if (!ifp) {
+	    printf("Cannot Open File");
+	    exit(100);
+	}
+
+	while (!feof(ifp)) {
+	    fscanf(ifp, "%lf \n", &temp[i]);
+	    i++;
+	}
+
+	fclose(ifp);
+
+	N = i;
+	J = 4;
+
+	sig = (double*)malloc(sizeof(double)* N);
+	inp = (double*)malloc(sizeof(double)* N);
+	oup = (double*)malloc(sizeof(double)* N);
+
+	for(i = 0; i < N;++i) {
+		sig[i] = temp[i];
+	}
+
+	ifp = fopen("PieceRegular10.txt", "r");
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
+
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+
+	fclose(ifp);
+
+	for(i = 0; i < N;++i) {
+		inp[i] = temp[i];
+	}
+	obj = denoise_init(N,J,wname);
+	setDenoiseMethod(obj,"visushrink");// sureshrink is also the default. The other option with dwt and swt is visushrink.
+	// modwt works only with modwtshrink method
+	setDenoiseWTMethod(obj,method);// Default is dwt. the other options are swt and modwt
+	setDenoiseWTExtension(obj,ext);// Default is sym. the other option is per
+	setDenoiseParameters(obj,thresh,level);// Default for thresh is soft. Other option is hard
+	// Default for level is all. The other option is first
+
+	denoise(obj,inp,oup);
+
+	// Alternative to denoise_object
+	// Just use visushrink, modwtshrink and sureshrink functions
+	//visushrink(inp,N,J,wname,method,ext,thresh,level,oup);
+	//sureshrink(inp,N,J,wname,method,ext,thresh,level,oup);
+	// modwtshrink(sig,N,J,wname,cmethod,ext,thresh,oup); See modwtdenoisetest.c
+	//ofp = fopen("denoiseds.txt", "w");
+
+
+	printf("Signal - Noisy Signal Stats \n");
+	printf("RMSE %g\n",rmse(N,sig,inp));
+	printf("Corr Coeff %g\n",corrcoef(N,sig,inp));
+
+	printf("Signal - DeNoised Signal Stats \n");
+	printf("RMSE %g\n",rmse(N,sig,oup));
+		printf("Corr Coeff %g\n",corrcoef(N,sig,oup));
+
+	free(sig);
+	free(inp);
+	denoise_free(obj);
+	free(oup);
+	return 0;
+}

test/dwpttest.c

@@ -0,0 +1,61 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+int main() {
+	int i, J, N;
+	wave_object obj;
+	wpt_object wt;
+	double *inp, *oup, *diff;
+
+	char *name = "db4";
+	obj = wave_init(name);// Initialize the wavelet
+	N = 788 + 23;
+	inp = (double*)malloc(sizeof(double)* N);
+	oup = (double*)malloc(sizeof(double)* N);
+	diff = (double*)malloc(sizeof(double)* N);
+	for (i = 1; i < N + 1; ++i) {
+		//inp[i - 1] = -0.25*i*i*i + 25 * i *i + 10 * i;
+		inp[i - 1] = i;
+	}
+	J = 4;
+
+	wt = wpt_init(obj, N, J);// Initialize the wavelet transform Tree object
+	setDWPTExtension(wt, "per");// Options are "per" and "sym". Symmetric is the default option
+	setDWPTEntropy(wt, "logenergy", 0);
+
+	dwpt(wt, inp); // Discrete Wavelet Packet Transform
+
+	idwpt(wt, oup); // Inverse Discrete Wavelet Packet Transform
+
+	for (i = 0; i < N; ++i) {
+		diff[i] = (inp[i] - oup[i])/inp[i];
+	}
+
+	wpt_summary(wt); // Tree Summary
+
+	printf("\n MAX %g \n", absmax(diff, wt->siglength)); // If Reconstruction succeeded then the output should be a small value.
+
+	free(inp);
+	free(oup);
+	free(diff);
+	wave_free(obj);
+	wpt_free(wt);
+	return 0;
+}

test/dwt2test.c

@@ -0,0 +1,87 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+double generate_rnd() {
+	double rnd;
+
+	rnd = (double) (rand() % 100 + 1);
+
+	return rnd;
+}
+
+int main() {
+	wave_object obj;
+	wt2_object wt;
+	int i, k, J, rows, cols,N;
+	double *inp, *wavecoeffs,*oup,*diff;
+	double *cLL;
+	int ir, ic;
+	double amax;
+
+	rows = 32;
+	cols = 30;
+    N = rows*cols;
+
+	char *name = "db2";
+	obj = wave_init(name);// Initialize the wavelet
+	srand(time(0));
+	inp = (double*)calloc(N, sizeof(double));
+	oup = (double*)calloc(N, sizeof(double));
+	diff = (double*)calloc(N, sizeof(double));
+    
+	J = 3;
+
+	wt = wt2_init(obj, "dwt", rows,cols, J);
+
+	for (i = 0; i < rows; ++i) {
+		for (k = 0; k < cols; ++k) {
+			//inp[i*cols + k] = i*cols + k;
+			inp[i*cols + k] = generate_rnd();
+			oup[i*cols + k] = 0.0;
+		}
+	}
+
+	wavecoeffs = dwt2(wt, inp);
+
+	cLL = getWT2Coeffs(wt, wavecoeffs, 1, "D", &ir, &ic);
+
+	dispWT2Coeffs(cLL, ir, ic);
+	
+	idwt2(wt, wavecoeffs, oup);
+
+	for (i = 0; i < rows*cols; ++i) {
+		diff[i] = oup[i] - inp[i];
+	}
+
+	amax = absmax(diff, rows*cols);
+
+	wt2_summary(wt);
+
+	printf("Abs Max %g \n", amax);
+
+	wave_free(obj);
+	wt2_free(wt);
+	free(inp);
+	free(wavecoeffs);
+	free(oup);
+	free(diff);
+    return 0;
+}

test/dwttest.c

@@ -41,6 +41,8 @@ int main() {
 		i++;
 	}
 	N = 256;
+	
+	fclose(ifp);
 
 	inp = (double*)malloc(sizeof(double)* N);
 	out = (double*)malloc(sizeof(double)* N);

test/modwt2test.c

@@ -0,0 +1,82 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+double generate_rnd() {
+	double rnd;
+
+	rnd = (double) (rand() % 100 + 1);
+
+	return rnd;
+}
+
+int main() {
+	wave_object obj;
+	wt2_object wt;
+	int i, k, J, rows, cols,N,ir,ic;
+	double *inp, *wavecoeffs, *oup,*cLL,*diff;
+	double amax;
+	rows = 51;
+	cols = 40;
+	N = rows*cols;
+	char *name = "db2";
+	obj = wave_init(name);// Initialize the wavelet
+
+	inp = (double*)calloc(N, sizeof(double));
+	oup = (double*)calloc(N, sizeof(double));
+	diff = (double*)calloc(N, sizeof(double));
+	J = 2;
+
+	wt = wt2_init(obj, "modwt", rows, cols, J);
+
+	for (i = 0; i < rows; ++i) {
+		for (k = 0; k < cols; ++k) {
+			//inp[i*cols + k] = i*cols + k;
+            inp[i*cols + k] = generate_rnd();
+			oup[i*cols + k] = 0.0;
+		}
+	}
+
+	wavecoeffs = modwt2(wt, inp);
+
+	cLL = getWT2Coeffs(wt, wavecoeffs, J, "A", &ir, &ic);
+
+	//dispWT2Coeffs(cLL, ir, ic);
+	
+	imodwt2(wt, wavecoeffs, oup);
+
+	for (i = 0; i < N; ++i) {
+		diff[i] = oup[i] - inp[i];
+	}
+
+	amax = absmax(diff, N);
+
+	wt2_summary(wt);
+
+	printf("Abs Max %g \n", amax);
+	
+	wave_free(obj);
+	wt2_free(wt);
+	free(inp);
+	free(wavecoeffs);
+	free(oup);
+	free(diff);
+    return 0;
+}

test/modwtdenoisetest.c

@@ -0,0 +1,123 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "../header/wauxlib.h"
+
+static double rmse(int N,double *x,double *y) {
+	double rms;
+	int i;
+
+	rms = 0.0;
+
+	for(i = 0; i < N;++i) {
+		rms += (x[i] - y[i]) * (x[i] - y[i]);
+	}
+
+	rms = sqrt(rms/(double)N);
+
+	return rms;
+}
+
+static double corrcoef(int N,double *x,double *y) {
+	double cc,xm,ym,tx,ty,num,den1,den2;
+	int i;
+	xm = ym = 0.0;
+	for(i = 0; i < N;++i) {
+		xm += x[i];
+		ym += y[i];
+	}
+
+	xm = xm/N;
+	ym = ym / N;
+	num = den1 = den2 = 0.0;
+
+	for(i = 0; i < N;++i) {
+		tx = x[i] - xm;
+		ty = y[i] - ym;
+		num += (tx*ty);
+		den1 += (tx*tx);
+		den2 += (ty*ty);
+	}
+
+	cc = num / sqrt(den1*den2);
+
+	return cc;
+}
+
+int main() {
+	// gcc -Wall -I../header -L../Bin modwtdenoisetest.c -o modwtdenoise -lwauxlib -lwavelib -lm
+    /*
+    modwtshrink can also be called from the denoise object. See denoisetest.c for more information
+    */
+	double *sig,*inp,*oup;
+	int i,N,J;
+	FILE *ifp;
+
+	double temp[2400];
+
+	char *wname = "db5";
+	char *ext = "per";// The other option sym is only available with "fft" cmethod
+	char *thresh = "soft";
+	char *cmethod = "direct";// The other option is "fft"
+
+	ifp = fopen("pieceregular1024.txt", "r");
+	i = 0;
+	if (!ifp) {
+	    printf("Cannot Open File");
+	    exit(100);
+	}
+
+	while (!feof(ifp)) {
+	    fscanf(ifp, "%lf \n", &temp[i]);
+	    i++;
+	}
+
+	fclose(ifp);
+
+	N = i;
+	J = 4;
+
+	sig = (double*)malloc(sizeof(double)* N);
+	inp = (double*)malloc(sizeof(double)* N);
+	oup = (double*)malloc(sizeof(double)* N);
+
+	for(i = 0; i < N;++i) {
+		sig[i] = temp[i];
+	}
+
+	ifp = fopen("PieceRegular10.txt", "r");
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
+
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+
+	fclose(ifp);
+
+	for(i = 0; i < N;++i) {
+		inp[i] = temp[i];
+	}
+	
+    modwtshrink(sig,N,J,wname,cmethod,ext,thresh,oup);
+
+
+	printf("Signal - Noisy Signal Stats \n");
+	printf("RMSE %g\n",rmse(N,sig,inp));
+	printf("Corr Coeff %g\n",corrcoef(N,sig,inp));
+
+	printf("Signal - DeNoised Signal Stats \n");
+	printf("RMSE %g\n",rmse(N,sig,oup));
+		printf("Corr Coeff %g\n",corrcoef(N,sig,oup));
+
+	free(sig);
+	free(inp);
+	free(oup);
+	return 0;
+}

test/modwttest.c

@@ -42,6 +42,8 @@ int main() {
 		i++;
 	}
 	N = 177;
+	
+	fclose(ifp);
 
 	inp = (double*)malloc(sizeof(double)* N);
 	out = (double*)malloc(sizeof(double)* N);

test/noisyheavisine.txt

@@ -0,0 +1,512 @@
+1.0657
+    1.0006
+    0.6117
+    0.4044
+    0.7049
+    0.8592
+    0.4489
+    0.6764
+    1.4804
+    1.3729
+    1.1876
+    0.8359
+    1.0592
+    1.6929
+    1.5832
+    1.4987
+    1.9193
+    0.9968
+    2.2814
+    1.8615
+    1.8153
+    2.2991
+    2.3962
+    2.6282
+    2.6752
+    2.8719
+    2.3616
+    2.2559
+    2.3763
+    2.5818
+    3.2529
+    2.2118
+    2.7430
+    3.0733
+    3.3577
+    2.7022
+    3.1119
+    3.2637
+    3.2952
+    3.4738
+    3.3380
+    3.1535
+    3.3264
+    3.5341
+    4.1180
+    3.7819
+    3.6842
+    3.8126
+    3.9201
+    4.4759
+    3.9019
+    4.4292
+    3.9813
+    4.2000
+    4.1341
+    4.6312
+    3.5677
+    3.5918
+    4.1196
+    4.3429
+    4.0506
+    3.9147
+    4.5915
+    3.7525
+    4.3211
+    4.7630
+    4.1468
+    3.6378
+    4.9381
+    3.6221
+    4.0528
+    4.2498
+    4.3916
+    3.8548
+    4.1120
+    3.1226
+    3.8218
+    3.4635
+    3.0588
+    4.1349
+    4.2391
+    3.7100
+    3.1359
+    4.4892
+    3.3203
+    3.7948
+    3.3253
+    3.8491
+    3.3987
+    3.0231
+    3.1837
+    3.0502
+    3.5930
+    3.2456
+    3.0435
+    2.5884
+    2.4189
+    2.9470
+    3.1980
+    2.2149
+    1.9529
+    3.0684
+    2.3438
+    1.9030
+    1.6420
+    2.5492
+    2.5800
+    2.5930
+    2.0497
+    1.7434
+    2.4729
+    2.0337
+    0.7426
+    1.3418
+    2.2373
+    0.8947
+    1.1703
+    0.8237
+    1.4013
+    0.5204
+    0.7086
+    0.2940
+    0.4183
+    0.8239
+    0.3488
+    0.1901
+   -0.2755
+   -0.5884
+    0.0459
+   -0.0061
+   -0.1662
+   -1.1395
+   -0.2187
+   -0.4940
+   -0.1904
+   -0.7303
+   -0.8049
+   -1.2128
+   -1.4796
+   -1.2787
+   -0.6750
+   -0.6300
+   -1.9705
+   -1.3766
+   -1.6300
+   -1.6662
+   -1.4749
+   -2.2873
+   -2.3715
+   -1.9438
+   -2.0246
+   -2.3207
+   -2.3904
+   -4.0250
+   -4.7166
+   -4.6094
+   -4.3269
+   -4.5663
+   -4.1415
+   -4.8460
+   -4.9088
+   -5.1159
+   -5.2494
+   -4.7209
+   -5.6772
+   -5.3329
+   -6.3053
+   -5.0153
+   -5.1394
+   -5.0556
+   -5.8880
+   -5.5547
+   -5.5403
+   -6.3194
+   -6.3660
+   -5.9584
+   -5.1940
+   -4.9157
+   -5.7317
+   -6.5066
+   -5.7450
+   -5.7231
+   -5.9420
+   -5.8529
+   -6.6728
+   -6.5548
+   -5.6438
+   -6.0741
+   -6.6387
+   -6.1218
+   -6.3158
+   -5.7250
+   -6.0155
+   -5.8662
+   -5.7875
+   -6.3902
+   -5.9303
+   -6.0030
+   -5.6667
+   -5.1734
+   -5.7733
+   -5.9180
+   -5.8427
+   -6.0721
+   -6.4873
+   -5.5756
+   -5.9103
+   -5.5415
+   -5.6358
+   -5.8095
+   -5.4756
+   -5.2417
+   -5.4192
+   -5.3777
+   -5.7800
+   -4.8058
+   -4.7930
+   -6.2389
+   -5.9839
+   -4.9383
+   -5.3716
+   -5.4538
+   -3.8454
+   -5.1885
+   -5.2452
+   -4.5655
+   -4.9033
+   -4.9928
+   -5.0235
+   -4.6184
+   -4.0967
+   -4.8166
+   -4.1745
+   -4.0614
+   -4.1093
+   -3.4929
+   -3.5424
+   -2.5478
+   -4.1180
+   -3.4682
+   -3.1256
+   -3.5773
+   -3.0447
+   -2.4956
+   -2.7483
+   -2.6201
+   -2.9657
+   -2.6621
+   -2.8913
+   -2.6488
+   -2.5289
+   -1.9951
+   -2.1213
+   -2.2065
+   -2.2494
+   -2.0965
+   -2.3657
+   -1.5025
+   -1.2423
+   -2.0154
+   -0.8329
+   -1.6098
+   -1.2474
+   -1.0787
+   -1.2216
+   -1.0861
+   -1.3985
+   -0.8476
+   -0.4991
+    0.0904
+   -0.5278
+    0.1709
+   -0.5306
+   -0.8072
+   -0.4898
+   -0.3393
+    0.2191
+   -0.4752
+    0.0910
+   -0.2168
+   -0.7928
+    0.4831
+    0.1193
+    0.9896
+    0.4941
+    1.1458
+    1.1746
+    1.6752
+    0.6359
+    0.5090
+    1.4067
+    0.9310
+    1.0004
+    1.4047
+    1.4470
+    1.4776
+    1.8183
+    1.7719
+    1.0112
+    1.6877
+    1.3403
+    1.2260
+    1.7027
+    1.7228
+    1.5210
+    1.9816
+    2.0695
+    2.0422
+    1.6521
+    1.3538
+    1.6597
+    1.6981
+    1.9754
+    2.2320
+    2.8194
+    1.9796
+    1.9535
+    1.8937
+    1.6508
+    2.1684
+    1.9457
+    2.2100
+    2.3376
+    1.4828
+    2.3156
+    1.6363
+    1.6415
+    1.6262
+    1.7795
+    1.2742
+    2.1842
+    1.5837
+    2.1802
+    2.5516
+    1.8494
+    1.5392
+    1.8861
+    1.1615
+    1.5972
+    1.5326
+    1.4134
+    1.1573
+    0.9850
+    1.3061
+    1.5567
+    1.1001
+    0.5861
+    1.2758
+    1.4634
+    0.5481
+   -0.2347
+    1.2253
+    0.7886
+   -0.0569
+    0.3684
+    1.0031
+    0.2320
+    0.2774
+    0.3051
+    0.2066
+   -0.0612
+   -0.4045
+    0.2544
+    0.1755
+    0.1436
+    0.6782
+   -0.3352
+   -0.4587
+    1.8259
+    1.3455
+    1.8159
+    1.8610
+    1.4192
+    1.4261
+    1.0369
+    0.8564
+    0.4077
+    0.5913
+    0.0495
+    1.1516
+    0.2284
+   -0.0953
+   -0.2963
+    0.3560
+    0.0803
+    0.1577
+    0.1330
+   -0.4338
+    0.1264
+   -0.5193
+   -0.3638
+   -1.4667
+   -0.8071
+   -1.1646
+   -1.3581
+   -2.0099
+   -1.9754
+   -1.3684
+   -1.4739
+   -2.0044
+   -1.9212
+   -1.3331
+   -1.8712
+   -1.9120
+   -1.6113
+   -1.4759
+   -2.5851
+   -1.5004
+   -2.2432
+   -2.4955
+   -1.8040
+   -2.2723
+   -2.7506
+   -2.7914
+   -3.0147
+   -3.1889
+   -2.6117
+   -2.9667
+   -4.1494
+   -3.1516
+   -2.9235
+   -2.9130
+   -3.3368
+   -3.5575
+   -3.2338
+   -3.6494
+   -3.2499
+   -4.3063
+   -3.3651
+   -2.9785
+   -3.3652
+   -3.4743
+   -4.0558
+   -3.9807
+   -3.2774
+   -4.0199
+   -4.5528
+   -4.2490
+   -3.5356
+   -3.9068
+   -3.7764
+   -3.7189
+   -3.2039
+   -3.6964
+   -4.9097
+   -3.4455
+   -3.4451
+   -4.1807
+   -4.2489
+   -3.5878
+   -4.1839
+   -4.1409
+   -3.4127
+   -3.8450
+   -3.1923
+   -4.2415
+   -4.1260
+   -3.1998
+   -4.1118
+   -4.3941
+   -4.0991
+   -3.7035
+   -3.5813
+   -3.6244
+   -3.8968
+   -4.0114
+   -3.0996
+   -2.5770
+   -2.5784
+   -2.5148
+   -3.4869
+   -3.1257
+   -3.3178
+   -3.2136
+   -3.3256
+   -3.1696
+   -2.6366
+   -2.7773
+   -3.4404
+   -2.7195
+   -1.7045
+   -2.7076
+   -2.4246
+   -3.3656
+   -2.7804
+   -2.5646
+   -2.2261
+   -1.8682
+   -2.7736
+   -2.1846
+   -2.2518
+   -2.1819
+   -1.6506
+   -1.1380
+   -1.4379
+   -1.2677
+   -0.9920
+   -0.9576
+   -1.7788
+   -1.1704
+   -1.0133
+   -0.0132
+   -0.5174
+   -0.7500
+   -0.5398
+   -1.4065
+   -1.4180
+   -0.5397
+    0.2176
+    0.0255
+   -0.1699
+   -0.0142

test/sst_nino3.dat

@@ -0,0 +1,504 @@
+-0.15
+-0.30
+-0.14
+-0.41
+-0.46
+-0.66
+-0.50
+-0.80
+-0.95
+-0.72
+-0.31
+-0.71
+-1.04
+-0.77
+-0.86
+-0.84
+-0.41
+-0.49
+-0.48
+-0.72
+-1.21
+-0.80
+ 0.16
+ 0.46
+ 0.40
+ 1.00
+ 2.17
+ 2.50
+ 2.34
+ 0.80
+ 0.14
+-0.06
+-0.34
+-0.71
+-0.34
+-0.73
+-0.48
+-0.11
+ 0.22
+ 0.51
+ 0.51
+ 0.25
+-0.10
+-0.33
+-0.42
+-0.23
+-0.53
+-0.44
+-0.30
+ 0.15
+ 0.09
+ 0.19
+-0.06
+ 0.25
+ 0.30
+ 0.81
+ 0.26
+ 0.10
+ 0.34
+ 1.01
+-0.31
+-0.90
+-0.73
+-0.92
+-0.73
+-0.31
+-0.03
+ 0.12
+ 0.37
+ 0.82
+ 1.22
+ 1.83
+ 1.60
+ 0.34
+-0.72
+-0.87
+-0.85
+-0.40
+-0.39
+-0.65
+ 0.07
+ 0.67
+ 0.39
+ 0.03
+-0.17
+-0.76
+-0.87
+-1.36
+-1.10
+-0.99
+-0.78
+-0.93
+-0.87
+-0.44
+-0.34
+-0.50
+-0.39
+-0.04
+ 0.42
+ 0.62
+ 0.17
+ 0.23
+ 1.03
+ 1.54
+ 1.09
+ 0.01
+ 0.12
+-0.27
+-0.47
+-0.41
+-0.37
+-0.36
+-0.39
+ 0.43
+ 1.05
+ 1.58
+ 1.25
+ 0.86
+ 0.60
+ 0.21
+ 0.19
+-0.23
+-0.29
+ 0.18
+ 0.12
+ 0.71
+ 1.42
+ 1.59
+ 0.93
+-0.25
+-0.66
+-0.95
+-0.47
+ 0.06
+ 0.70
+ 0.81
+ 0.78
+ 1.43
+ 1.22
+ 1.05
+ 0.44
+-0.35
+-0.67
+-0.84
+-0.66
+-0.45
+-0.12
+-0.20
+-0.16
+-0.47
+-0.52
+-0.79
+-0.80
+-0.62
+-0.86
+-1.29
+-1.04
+-1.05
+-0.75
+-0.81
+-0.90
+-0.25
+ 0.62
+ 1.22
+ 0.96
+ 0.21
+-0.11
+-0.25
+-0.24
+-0.43
+ 0.23
+ 0.67
+ 0.78
+ 0.41
+ 0.98
+ 1.28
+ 1.45
+ 1.02
+ 0.03
+-0.59
+-1.34
+-0.99
+-1.49
+-1.74
+-1.33
+-0.55
+-0.51
+-0.36
+-0.99
+ 0.32
+ 1.04
+ 1.41
+ 0.99
+ 0.66
+ 0.50
+ 0.22
+ 0.71
+-0.16
+ 0.38
+ 0.00
+-1.11
+-1.04
+ 0.05
+-0.64
+-0.34
+-0.50
+-1.85
+-0.94
+-0.78
+ 0.29
+ 0.27
+ 0.69
+-0.06
+-0.83
+-0.80
+-1.02
+-0.96
+-0.09
+ 0.62
+ 0.87
+ 1.03
+ 0.70
+-0.10
+-0.31
+ 0.04
+-0.46
+ 0.04
+ 0.24
+-0.08
+-0.28
+ 0.06
+ 0.05
+-0.31
+ 0.11
+ 0.27
+ 0.26
+ 0.04
+ 0.12
+ 1.11
+ 1.53
+ 1.23
+ 0.17
+-0.18
+-0.56
+ 0.05
+ 0.41
+ 0.22
+ 0.04
+-0.19
+-0.46
+-0.65
+-1.06
+-0.54
+ 0.14
+ 0.25
+-0.21
+-0.73
+-0.43
+ 0.48
+ 0.26
+ 0.05
+ 0.11
+-0.27
+-0.08
+-0.10
+ 0.29
+-0.15
+-0.28
+-0.55
+-0.44
+-1.40
+-0.55
+-0.69
+ 0.58
+ 0.37
+ 0.42
+ 1.83
+ 1.23
+ 0.65
+ 0.41
+ 1.03
+ 0.64
+-0.07
+ 0.98
+ 0.36
+-0.30
+-1.33
+-1.39
+-0.94
+ 0.34
+-0.00
+-0.15
+ 0.06
+ 0.39
+ 0.36
+-0.49
+-0.53
+ 0.35
+ 0.07
+-0.24
+ 0.20
+-0.22
+-0.68
+-0.44
+ 0.02
+-0.22
+-0.30
+-0.59
+ 0.10
+-0.02
+-0.27
+-0.60
+-0.48
+-0.37
+-0.53
+-1.35
+-1.22
+-0.99
+-0.34
+-0.79
+-0.24
+ 0.02
+ 0.69
+ 0.78
+ 0.17
+-0.17
+-0.29
+-0.27
+ 0.31
+ 0.44
+ 0.38
+ 0.24
+-0.13
+-0.89
+-0.76
+-0.71
+-0.37
+-0.59
+-0.63
+-1.47
+-0.40
+-0.18
+-0.37
+-0.43
+-0.06
+ 0.61
+ 1.33
+ 1.19
+ 1.13
+ 0.31
+ 0.14
+ 0.03
+ 0.21
+ 0.15
+-0.22
+-0.02
+ 0.03
+-0.17
+ 0.12
+-0.35
+-0.06
+ 0.38
+-0.45
+-0.32
+-0.33
+-0.49
+-0.14
+-0.56
+-0.18
+ 0.46
+ 1.09
+ 1.04
+ 0.23
+-0.99
+-0.59
+-0.92
+-0.28
+ 0.52
+ 1.31
+ 1.45
+ 0.61
+-0.11
+-0.18
+-0.39
+-0.39
+-0.36
+-0.50
+-0.81
+-1.10
+-0.29
+ 0.57
+ 0.68
+ 0.78
+ 0.78
+ 0.63
+ 0.98
+ 0.49
+-0.42
+-1.34
+-1.20
+-1.18
+-0.65
+-0.42
+-0.97
+-0.28
+ 0.77
+ 1.77
+ 2.22
+ 1.05
+-0.67
+-0.99
+-1.52
+-1.17
+-0.22
+-0.04
+-0.45
+-0.46
+-0.75
+-0.70
+-1.38
+-1.15
+-0.01
+ 0.97
+ 1.10
+ 0.68
+-0.02
+-0.04
+ 0.47
+ 0.30
+-0.55
+-0.51
+-0.09
+-0.01
+ 0.34
+ 0.61
+ 0.58
+ 0.33
+ 0.38
+ 0.10
+ 0.18
+-0.30
+-0.06
+-0.28
+ 0.12
+ 0.58
+ 0.89
+ 0.93
+ 2.39
+ 2.44
+ 1.92
+ 0.64
+-0.24
+ 0.27
+-0.13
+-0.16
+-0.54
+-0.13
+-0.37
+-0.78
+-0.22
+ 0.03
+ 0.25
+ 0.31
+ 1.03
+ 1.10
+ 1.05
+ 1.11
+ 1.28
+ 0.57
+-0.55
+-1.16
+-0.99
+-0.38
+ 0.01
+-0.29
+ 0.09
+ 0.46
+ 0.57
+ 0.24
+ 0.39
+ 0.49
+ 0.86
+ 0.51
+ 0.95
+ 1.25
+ 1.33
+-0.00
+ 0.34
+ 0.66
+ 1.11
+ 0.34
+ 0.48
+ 0.56
+ 0.39
+-0.17
+ 1.04
+ 0.77
+ 0.12
+-0.35
+-0.22
+ 0.08
+-0.08
+-0.18
+-0.06

test/swt2test.c

@@ -0,0 +1,83 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+double generate_rnd() {
+	double rnd;
+
+	rnd = (double) (rand() % 100 + 1);
+
+	return rnd;
+}
+
+int main() {
+	wave_object obj;
+	wt2_object wt;
+	int i, k, J, rows, cols;
+	double *inp, *wavecoeffs, *oup, *cLL,*diff;
+	double amax;
+	int ir, ic, N;
+	rows = 64;
+	cols = 48;
+
+	char *name = "bior3.1";
+	obj = wave_init(name);// Initialize the wavelet
+	N = rows*cols;
+	inp = (double*)calloc(N, sizeof(double));
+	oup = (double*)calloc(N, sizeof(double));
+	diff = (double*)calloc(N, sizeof(double));
+	J = 2;
+
+	wt = wt2_init(obj, "swt", rows, cols, J);
+
+	for (i = 0; i < rows; ++i) {
+		for (k = 0; k < cols; ++k) {
+			//inp[i*cols + k] = i*cols + k;
+            inp[i*cols + k] = generate_rnd();
+			oup[i*cols + k] = 0.0;
+		}
+	}
+	
+	wavecoeffs = swt2(wt, inp);
+
+	cLL = getWT2Coeffs(wt, wavecoeffs, J, "A", &ir, &ic);
+
+	dispWT2Coeffs(cLL, ir, ic);
+	
+	iswt2(wt, wavecoeffs, oup);
+
+	for (i = 0; i < N; ++i) {
+		diff[i] = oup[i] - inp[i];
+	}
+
+	amax = absmax(diff, N);
+
+	wt2_summary(wt);
+
+	printf("Abs Max %g \n", amax);
+	
+	wave_free(obj);
+	wt2_free(wt);
+	free(inp);
+	free(wavecoeffs);
+	free(oup);
+	free(diff);
+    return 0;
+}

test/swttest.c

@@ -41,7 +41,7 @@ int main() {
 		i++;
 	}
 	N = 256;
-
+	fclose(ifp);
 	inp = (double*)malloc(sizeof(double)* N);
 	out = (double*)malloc(sizeof(double)* N);
 	diff = (double*)malloc(sizeof(double)* N);

test/wtreetest.c

@@ -0,0 +1,47 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+int main() {
+	int i, J, N, len;
+	int X, Y;
+	wave_object obj;
+	wtree_object wt;
+	double *inp, *oup;
+
+	char *name = "db3";
+	obj = wave_init(name);// Initialize the wavelet
+	N = 147;
+	inp = (double*)malloc(sizeof(double)* N);
+	for (i = 1; i < N + 1; ++i) {
+		inp[i - 1] = -0.25*i*i*i + 25 * i *i + 10 * i;
+	}
+	J = 3;
+
+	wt = wtree_init(obj, N, J);// Initialize the wavelet transform object
+	setWTREEExtension(wt, "sym");// Options are "per" and "sym". Symmetric is the default option
+
+	wtree(wt, inp);
+	wtree_summary(wt);
+	X = 3;
+	Y = 5;
+	len = getWTREENodelength(wt, X);
+	printf("\n %d", len);
+	printf("\n");
+	oup = (double*)malloc(sizeof(double)* len);
+
+	printf("Node [%d %d] Coefficients : \n",X,Y);
+	getWTREECoeffs(wt, X, Y, oup, len);
+	for (i = 0; i < len; ++i) {
+		printf("%g ", oup[i]);
+	}
+	printf("\n");
+
+	free(inp);
+	free(oup);
+	wave_free(obj);
+	wtree_free(wt);
+	return 0;
+}

unitTests/CMakeLists.txt

@@ -0,0 +1,3 @@
+
+add_subdirectory(wavelibBoostTests) 
+

unitTests/wavelibBoostTests/CMakeLists.txt

@@ -0,0 +1,27 @@
+                                
+set(SOURCE_FILES
+		tst_dwt.cpp
+
+	)
+
+
+
+
+
+add_executable(wavelibLibTests ${SOURCE_FILES} )
+
+add_test(NAME wavelibLibTests WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/test  COMMAND wavelibLibTests)
+
+add_dependencies(wavelibLibTests wavelib)
+target_link_libraries(wavelibLibTests wavelib)
+
+target_include_directories(wavelibLibTests PUBLIC
+							${CMAKE_CURRENT_SOURCE_DIR}/../../header
+                            )
+                            
+
+install(TARGETS wavelibLibTests
+    RUNTIME DESTINATION bin
+    LIBRARY DESTINATION tests
+    ARCHIVE DESTINATION tests
+    )