Compare commits
75 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
b865a273d2 | ||
|
|
f906ae8d66 | ||
|
|
a26e271971 | ||
|
|
bd2314727e | ||
|
|
310b7759f1 | ||
|
|
de4b96e123 | ||
|
|
e35df26b18 | ||
|
|
59bd91d4fa | ||
|
|
1bde2f9824 | ||
|
|
1b34b65fb6 | ||
|
|
87c0e462dc | ||
|
|
22145f26fd | ||
|
|
4b187b1e02 | ||
|
|
661432388d | ||
|
|
46ef7cb1b9 | ||
|
|
72a3e5c580 | ||
|
|
e1baa7fc17 | ||
|
|
ee88fece29 | ||
|
|
fa68a811af | ||
|
|
cb01ae39ec | ||
|
|
ede7d68eb0 | ||
|
|
615ea94091 | ||
|
|
4b80222408 | ||
|
|
d29e08d06b | ||
|
|
6e647fcf6d | ||
|
|
9a0a271e4d | ||
|
|
0297e002d5 | ||
|
|
318e9b0f86 | ||
|
|
eeae494901 | ||
|
|
35c426bb12 | ||
|
|
5c9a4b5a88 | ||
|
|
b511fe864d | ||
|
|
3b92f27042 | ||
|
|
69bf4dc0b0 | ||
|
|
dc38e9f39b | ||
|
|
a2c709715b | ||
|
|
2b3b78acbd | ||
|
|
2e94a2f078 | ||
|
|
f0339fbc4b | ||
|
|
79d0802f39 | ||
|
|
0c7cd6018b | ||
|
|
47b28c620d | ||
|
|
e4a4d1a17d | ||
|
|
e2ff743cbe | ||
|
|
da84e602bc | ||
|
|
895e1a1f3b | ||
|
|
32738628ee | ||
|
|
ce2a416bda | ||
|
|
552034c74f | ||
|
|
cf869ef0c4 | ||
|
|
aad90aebfc | ||
|
|
f691347822 | ||
|
|
77982d9ee5 | ||
|
|
2cbc33003b | ||
|
|
d66c84f96c | ||
|
|
d146251571 | ||
|
|
e39fce2a16 | ||
|
|
effbfb463b | ||
|
|
04548c820a | ||
|
|
8d1d842c9b | ||
|
|
dd10c33a1a | ||
|
|
d68b19c1c8 | ||
|
|
0ed5135772 | ||
|
|
9ce0ad4830 | ||
|
|
4075ae5bc8 | ||
|
|
1c0644364a | ||
|
|
f7786d4dd5 | ||
|
|
5e5fe8ffa8 | ||
|
|
ebcc2f9860 | ||
|
|
4ce49a05e9 | ||
|
|
487de13131 | ||
|
|
66dc615ac8 | ||
|
|
2197daab4d | ||
|
|
74d5dd3b3f | ||
|
|
23d5075329 |
35
.gitignore
vendored
Normal file
35
.gitignore
vendored
Normal file
@@ -0,0 +1,35 @@
|
||||
#Folders Ignore
|
||||
Bin/
|
||||
Testing/
|
||||
|
||||
#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
|
||||
36
.travis.yml
Normal file
36
.travis.yml
Normal file
@@ -0,0 +1,36 @@
|
||||
sudo: false
|
||||
language: cpp
|
||||
os:
|
||||
- linux
|
||||
- osx
|
||||
|
||||
compiler:
|
||||
- gcc
|
||||
- clang
|
||||
|
||||
addons:
|
||||
apt:
|
||||
sources:
|
||||
- boost-latest
|
||||
- ubuntu-toolchain-r-test
|
||||
packages:
|
||||
- libboost1.55-all-dev
|
||||
|
||||
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
|
||||
51
CMakeLists.txt
Normal file
51
CMakeLists.txt
Normal file
@@ -0,0 +1,51 @@
|
||||
cmake_minimum_required(VERSION 2.8.0 FATAL_ERROR)
|
||||
|
||||
set(PROJECT_NAME wavelib)
|
||||
project(${PROJECT_NAME} CXX C)
|
||||
|
||||
|
||||
# src root path
|
||||
set(WAVELIB_SRC_ROOT ${PROJECT_SOURCE_DIR} CACHE PATH "Wavelib source root")
|
||||
# binary output by default
|
||||
set(COMMON_BIN_PATH ${CMAKE_BINARY_DIR}/Bin)
|
||||
set(LIBRARY_OUTPUT_PATH ${COMMON_BIN_PATH}/${CMAKE_BUILD_TYPE})
|
||||
set(EXECUTABLE_OUTPUT_PATH ${COMMON_BIN_PATH}/${CMAKE_BUILD_TYPE})
|
||||
# set where to find additional cmake modules if any
|
||||
set(CMAKE_MODULE_PATH ${WAVELIB_SRC_ROOT}/cmake ${CMAKE_MODULE_PATH})
|
||||
|
||||
|
||||
set(WAVELIB_VERSION "1.0.0" CACHE STRING "Wavelib version" FORCE)
|
||||
message(">>> Building Wavelib version: ${WAVELIB_VERSION}")
|
||||
message(">>> EXECUTABLE_OUTPUT_PATH = ${EXECUTABLE_OUTPUT_PATH}")
|
||||
|
||||
option(BUILD_UT "Enable Unit test" ON)
|
||||
|
||||
# cleanup prefix lib for Unix-like OSes
|
||||
set(CMAKE_SHARED_MODULE_PREFIX)
|
||||
|
||||
# install target to this folder by default
|
||||
set(WAVELIB_BINARY_DIR ${WAVELIB_SRC_ROOT}/bin)
|
||||
if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
|
||||
set(CMAKE_INSTALL_PREFIX "${WAVELIB_BINARY_DIR}" CACHE PATH "default install path" FORCE)
|
||||
endif()
|
||||
|
||||
# make include globaly visible
|
||||
set(PROJECT_WIDE_INCLUDE ${WAVELIB_SRC_ROOT}/include)
|
||||
include_directories(${PROJECT_WIDE_INCLUDE})
|
||||
|
||||
|
||||
|
||||
|
||||
include_directories(${COMMON_BIN_PATH})
|
||||
|
||||
if(BUILD_UT)
|
||||
include(CTest)
|
||||
enable_testing()
|
||||
add_subdirectory(unitTests)
|
||||
endif()
|
||||
|
||||
add_subdirectory(src)
|
||||
add_subdirectory(auxiliary)
|
||||
add_subdirectory(test)
|
||||
|
||||
install(DIRECTORY ${WAVELIB_SRC_ROOT}/include/ DESTINATION include FILES_MATCHING PATTERN "*.h")
|
||||
@@ -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:
|
||||
|
||||
12
README.md
12
README.md
@@ -1,9 +1,9 @@
|
||||
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.
|
||||
|
||||
@@ -11,6 +11,14 @@ SWT/ISWT Stationary Wavelet Transform. It works only for signal lengths that are
|
||||
|
||||
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.
|
||||
|
||||
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/
|
||||
|
||||
16
appveyor.yml
Normal file
16
appveyor.yml
Normal file
@@ -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
|
||||
18
auxiliary/CMakeLists.txt
Normal file
18
auxiliary/CMakeLists.txt
Normal file
@@ -0,0 +1,18 @@
|
||||
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
|
||||
|
||||
set(SOURCE_FILES denoise.c
|
||||
waux.c
|
||||
)
|
||||
|
||||
set(HEADER_FILES denoise.h
|
||||
waux.h
|
||||
)
|
||||
|
||||
add_library(wauxlib STATIC ${SOURCE_FILES} ${HEADER_FILES})
|
||||
|
||||
target_link_libraries(wauxlib wavelib)
|
||||
|
||||
set_property(TARGET wauxlib PROPERTY FOLDER "lib")
|
||||
|
||||
target_include_directories(wauxlib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
|
||||
|
||||
359
auxiliary/denoise.c
Normal file
359
auxiliary/denoise.c
Normal file
@@ -0,0 +1,359 @@
|
||||
|
||||
#include "denoise.h"
|
||||
|
||||
denoise_object denoise_init(int length, int J,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,"first");
|
||||
strcpy(obj->thresh,"soft");
|
||||
strcpy(obj->wmethod,"dwt");
|
||||
|
||||
return obj;
|
||||
}
|
||||
|
||||
void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,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,char *wname,char *method,char *ext,char *thresh,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;
|
||||
}
|
||||
printf(" \n");
|
||||
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 denoise(denoise_object obj, double *signal,double *denoised) {
|
||||
if(!strcmp(obj->dmethod,"sureshrink")) {
|
||||
sureshrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);
|
||||
} else if(!strcmp(obj->dmethod,"visushrink")) {
|
||||
visushrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);;
|
||||
} else {
|
||||
printf("Acceptable Denoising methods are - sureshrink and visushrink\n");
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
|
||||
void setDenoiseMethod(denoise_object obj, char *dmethod) {
|
||||
if (!strcmp(dmethod, "sureshrink")) {
|
||||
strcpy(obj->dmethod, "sureshrink");
|
||||
}
|
||||
else if (!strcmp(dmethod, "visushrink")) {
|
||||
strcpy(obj->dmethod, "visushrink");
|
||||
}
|
||||
else {
|
||||
printf("Acceptable Denoising methods are - sureshrink and visushrink\n");
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
|
||||
void setDenoiseWTMethod(denoise_object obj, char *wmethod) {
|
||||
if (!strcmp(wmethod, "dwt")) {
|
||||
strcpy(obj->wmethod, "dwt");
|
||||
}
|
||||
else if (!strcmp(wmethod, "swt")) {
|
||||
strcpy(obj->wmethod, "swt");
|
||||
}
|
||||
else {
|
||||
printf("Wavelet decomposition method can be either dwt or swt");
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
|
||||
void setDenoiseWTExtension(denoise_object obj, 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, char *thresh,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);
|
||||
}
|
||||
52
auxiliary/denoise.h
Normal file
52
auxiliary/denoise.h
Normal file
@@ -0,0 +1,52 @@
|
||||
/*
|
||||
Copyright (c) 2017, Rafat Hussain
|
||||
*/
|
||||
#ifndef DENOISE_H_
|
||||
#define DENOISE_H_
|
||||
|
||||
|
||||
#include "waux.h"
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
typedef struct denoise_set* denoise_object;
|
||||
|
||||
denoise_object denoise_init(int length, int J,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 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,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised);
|
||||
|
||||
void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised);
|
||||
|
||||
void denoise(denoise_object obj, double *signal,double *denoised);
|
||||
|
||||
void setDenoiseMethod(denoise_object obj, char *dmethod);
|
||||
|
||||
void setDenoiseWTMethod(denoise_object obj, char *wmethod);
|
||||
|
||||
void setDenoiseWTExtension(denoise_object obj, char *extension);
|
||||
|
||||
void setDenoiseParameters(denoise_object obj, char *thresh,char *level);
|
||||
|
||||
void denoise_free(denoise_object object);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#endif /* DENOISE_H_ */
|
||||
322
auxiliary/waux.c
Normal file
322
auxiliary/waux.c
Normal file
@@ -0,0 +1,322 @@
|
||||
#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(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(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(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);
|
||||
}
|
||||
|
||||
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,char *ctype,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(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(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;
|
||||
}
|
||||
|
||||
}
|
||||
52
auxiliary/waux.h
Normal file
52
auxiliary/waux.h
Normal file
@@ -0,0 +1,52 @@
|
||||
/*
|
||||
* 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(double* vec, int N);
|
||||
|
||||
double var(double* vec, int N);
|
||||
|
||||
double median(double *x, int N);
|
||||
|
||||
double mad(double *x, int N);
|
||||
|
||||
int minindex(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 getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, int J,double *lpr,
|
||||
double *hpr,int lf,int siglength,double *reccoeff);
|
||||
|
||||
void autocovar(double* vec,int N,double* acov, int M);
|
||||
|
||||
void autocorr(double* vec,int N,double* acorr, int M);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
#endif /* AUXILIARY_WAUX_H_ */
|
||||
56
header/wauxlib.h
Normal file
56
header/wauxlib.h
Normal file
@@ -0,0 +1,56 @@
|
||||
/*
|
||||
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,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 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,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised);
|
||||
|
||||
void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised);
|
||||
|
||||
void denoise(denoise_object obj, double *signal,double *denoised);
|
||||
|
||||
void setDenoiseMethod(denoise_object obj, char *dmethod);
|
||||
|
||||
void setDenoiseWTMethod(denoise_object obj, char *wmethod);
|
||||
|
||||
void setDenoiseWTExtension(denoise_object obj, char *extension);
|
||||
|
||||
void setDenoiseParameters(denoise_object obj, char *thresh,char *level);
|
||||
|
||||
void denoise_free(denoise_object object);
|
||||
|
||||
void getDWTRecCoeff(double *coeff,int *length,char *ctype,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_ */
|
||||
151
header/wavelib.h
151
header/wavelib.h
@@ -5,10 +5,25 @@
|
||||
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);
|
||||
@@ -83,21 +98,109 @@ struct wt_set{
|
||||
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 cfftset;
|
||||
int zpad;
|
||||
int length[102];
|
||||
double *output;
|
||||
double params[0];
|
||||
};
|
||||
|
||||
void dwt11(wt_object wt, double *Vin, int M, double *Wout,
|
||||
double *Vout);
|
||||
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];
|
||||
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"
|
||||
|
||||
int N; //
|
||||
int nodes;
|
||||
int cfftset;
|
||||
int zpad;
|
||||
int length[102];
|
||||
double *output;
|
||||
int *nodelength;
|
||||
int *coeflength;
|
||||
double params[0];
|
||||
};
|
||||
|
||||
typedef struct wpt_set* wpt_object;
|
||||
|
||||
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(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];
|
||||
};
|
||||
|
||||
|
||||
void dwt(wt_object wt, double *inp);
|
||||
|
||||
void idwt(wt_object wt, double *dwtop);
|
||||
|
||||
void wtree(wtree_object wt, double *inp);
|
||||
|
||||
void dwpt(wpt_object wt, double *inp);
|
||||
|
||||
void idwpt(wpt_object wt, double *dwtop);
|
||||
|
||||
void swt(wt_object wt, double *inp);
|
||||
|
||||
void iswt(wt_object wt, double *swtop);
|
||||
@@ -108,16 +211,54 @@ void imodwt(wt_object wt, double *dwtop);
|
||||
|
||||
void setDWTExtension(wt_object wt, char *extension);
|
||||
|
||||
void setWTREEExtension(wtree_object wt, char *extension);
|
||||
|
||||
void setDWPTExtension(wpt_object wt, char *extension);
|
||||
|
||||
void setDWPTEntropy(wpt_object wt, char *entropy, double eparam);
|
||||
|
||||
void setWTConv(wt_object wt, 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, char *type, int power);
|
||||
|
||||
void setCWTScaleVector(cwt_object wt, double *scale, int J, double s0, double dj);
|
||||
|
||||
void setCWTPadding(cwt_object wt, int pad);
|
||||
|
||||
void cwt(cwt_object wt, double *inp);
|
||||
|
||||
void icwt(cwt_object wt, double *cwtop);
|
||||
|
||||
int getCWTScaleLength(int N);
|
||||
|
||||
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 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);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
@@ -125,5 +266,3 @@ void wt_free(wt_object object);
|
||||
|
||||
|
||||
#endif /* WAVELIB_H_ */
|
||||
|
||||
|
||||
|
||||
32
src/CMakeLists.txt
Normal file
32
src/CMakeLists.txt
Normal file
@@ -0,0 +1,32 @@
|
||||
|
||||
|
||||
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
|
||||
|
||||
set(SOURCE_FILES conv.c
|
||||
cwt.c
|
||||
cwtmath.c
|
||||
hsfft.c
|
||||
real.c
|
||||
wavefilt.c
|
||||
wavefunc.c
|
||||
wavelib.c
|
||||
wtmath.c
|
||||
)
|
||||
|
||||
set(HEADER_FILES conv.h
|
||||
cwt.h
|
||||
cwtmath.h
|
||||
hsfft.h
|
||||
real.h
|
||||
wavefilt.h
|
||||
wavefunc.h
|
||||
wavelib.h
|
||||
wtmath.h
|
||||
)
|
||||
|
||||
add_library(wavelib STATIC ${SOURCE_FILES} ${HEADER_FILES})
|
||||
|
||||
set_property(TARGET wavelib PROPERTY FOLDER "lib")
|
||||
|
||||
target_include_directories(wavelib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
|
||||
|
||||
@@ -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);
|
||||
}
|
||||
|
||||
420
src/cwt.c
Executable file
420
src/cwt.c
Executable file
@@ -0,0 +1,420 @@
|
||||
/*
|
||||
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"
|
||||
|
||||
static int factorial2(int N) {
|
||||
int factorial,i;
|
||||
|
||||
factorial = 1;
|
||||
|
||||
for (i = 1; i <= N;++i) {
|
||||
factorial *= i;
|
||||
}
|
||||
|
||||
return factorial;
|
||||
}
|
||||
|
||||
static double factorial3(int N) {
|
||||
int i;
|
||||
double factorial;
|
||||
|
||||
factorial = 1;
|
||||
|
||||
for (i = 1; i <= N; ++i) {
|
||||
factorial *= i;
|
||||
}
|
||||
|
||||
return factorial;
|
||||
}
|
||||
|
||||
double factorial(int N) {
|
||||
if (N > 40) {
|
||||
printf("This program is only valid for N <= 40 \n");
|
||||
return -1.0;
|
||||
}
|
||||
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 };
|
||||
|
||||
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 / 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(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;
|
||||
|
||||
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) / gamma(0.5);
|
||||
*val = coeff * gamma(((double)m + 1.0) / 2.0) / sqrt(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;
|
||||
}
|
||||
}
|
||||
29
src/cwt.h
Executable file
29
src/cwt.h
Executable file
@@ -0,0 +1,29 @@
|
||||
#ifndef CWT_H_
|
||||
#define CWT_H_
|
||||
|
||||
#include "wavefunc.h"
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
void cwavelet(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_ */
|
||||
|
||||
|
||||
296
src/cwtmath.c
Executable file
296
src/cwtmath.c
Executable file
@@ -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 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;
|
||||
}
|
||||
22
src/cwtmath.h
Executable file
22
src/cwtmath.h
Executable file
@@ -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 gamma(double x);
|
||||
|
||||
int nint(double N);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#endif /* WAVELIB_H_ */
|
||||
@@ -106,6 +106,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);
|
||||
}
|
||||
|
||||
|
||||
7338
src/wavefilt.c
7338
src/wavefilt.c
File diff suppressed because it is too large
Load Diff
@@ -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
|
||||
|
||||
210
src/wavefunc.c
Executable file
210
src/wavefunc.c
Executable file
@@ -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(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]);
|
||||
}
|
||||
}
|
||||
24
src/wavefunc.h
Executable file
24
src/wavefunc.h
Executable file
@@ -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_ */
|
||||
4015
src/wavelib.c
4015
src/wavelib.c
File diff suppressed because it is too large
Load Diff
315
src/wavelib.h
315
src/wavelib.h
@@ -1,85 +1,230 @@
|
||||
#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_ */
|
||||
|
||||
|
||||
/*
|
||||
Copyright (c) 2014, Rafat Hussain
|
||||
*/
|
||||
#ifndef WAVELIB_H_
|
||||
#define WAVELIB_H_
|
||||
|
||||
#include "wtmath.h"
|
||||
#include "cwt.h"
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
#if defined(_MSC_VER)
|
||||
#pragma warning(disable : 4200)
|
||||
#pragma warning(disable : 4996)
|
||||
#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);
|
||||
|
||||
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];
|
||||
};
|
||||
|
||||
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];
|
||||
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"
|
||||
|
||||
int N; //
|
||||
int nodes;
|
||||
int cfftset;
|
||||
int zpad;
|
||||
int length[102];
|
||||
double *output;
|
||||
int *nodelength;
|
||||
int *coeflength;
|
||||
double params[0];
|
||||
};
|
||||
|
||||
typedef struct wpt_set* wpt_object;
|
||||
|
||||
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(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];
|
||||
};
|
||||
|
||||
|
||||
void dwt(wt_object wt, double *inp);
|
||||
|
||||
void idwt(wt_object wt, double *dwtop);
|
||||
|
||||
void wtree(wtree_object wt, double *inp);
|
||||
|
||||
void dwpt(wpt_object wt, double *inp);
|
||||
|
||||
void idwpt(wpt_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 setWTREEExtension(wtree_object wt, char *extension);
|
||||
|
||||
void setDWPTExtension(wpt_object wt, char *extension);
|
||||
|
||||
void setDWPTEntropy(wpt_object wt, char *entropy, double eparam);
|
||||
|
||||
void setWTConv(wt_object wt, 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, char *type, int power);
|
||||
|
||||
void setCWTScaleVector(cwt_object wt, double *scale, int J, double s0, double dj);
|
||||
|
||||
void setCWTPadding(cwt_object wt, int pad);
|
||||
|
||||
void cwt(cwt_object wt, double *inp);
|
||||
|
||||
void icwt(cwt_object wt, double *cwtop);
|
||||
|
||||
int getCWTScaleLength(int N);
|
||||
|
||||
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 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);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#endif /* WAVELIB_H_ */
|
||||
|
||||
96
src/wtmath.c
96
src/wtmath.c
@@ -1,3 +1,6 @@
|
||||
/*
|
||||
Copyright (c) 2014, Rafat Hussain
|
||||
*/
|
||||
#include "wtmath.h"
|
||||
|
||||
int upsamp(double *x, int lenx, int M, double *y) {
|
||||
@@ -236,4 +239,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;
|
||||
}
|
||||
|
||||
@@ -1,3 +1,6 @@
|
||||
/*
|
||||
Copyright (c) 2014, Rafat Hussain
|
||||
*/
|
||||
#ifndef WTMATH_H_
|
||||
#define WTMATH_H_
|
||||
|
||||
@@ -23,9 +26,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_ */
|
||||
|
||||
32
test/CMakeLists.txt
Normal file
32
test/CMakeLists.txt
Normal file
@@ -0,0 +1,32 @@
|
||||
add_executable(cwttest cwttest.c)
|
||||
|
||||
target_link_libraries(cwttest wavelib m)
|
||||
|
||||
add_executable(dwttest dwttest.c)
|
||||
|
||||
target_link_libraries(dwttest wavelib m)
|
||||
|
||||
add_executable(swttest swttest.c)
|
||||
|
||||
target_link_libraries(swttest wavelib m)
|
||||
|
||||
add_executable(modwttest modwttest.c)
|
||||
|
||||
target_link_libraries(modwttest wavelib m)
|
||||
|
||||
add_executable(dwpttest dwpttest.c)
|
||||
|
||||
target_link_libraries(dwpttest wavelib m)
|
||||
|
||||
add_executable(wtreetest wtreetest.c)
|
||||
|
||||
target_link_libraries(wtreetest wavelib m)
|
||||
|
||||
add_executable(denoisetest denoisetest.c)
|
||||
|
||||
target_link_libraries(denoisetest wauxlib wavelib m)
|
||||
|
||||
set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest
|
||||
PROPERTIES
|
||||
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/test"
|
||||
)
|
||||
1024
test/PieceRegular10.txt
Normal file
1024
test/PieceRegular10.txt
Normal file
File diff suppressed because it is too large
Load Diff
109
test/cwttest.c
Normal file
109
test/cwttest.c
Normal file
@@ -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;
|
||||
}
|
||||
135
test/denoisetest.c
Normal file
135
test/denoisetest.c
Normal file
@@ -0,0 +1,135 @@
|
||||
#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";
|
||||
char *ext = "sym";
|
||||
char *thresh = "soft";
|
||||
char *level = "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 is visushrink
|
||||
setDenoiseWTMethod(obj,method);// Default is dwt. the other option is swt
|
||||
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 first. The other option is all
|
||||
|
||||
denoise(obj,inp,oup);
|
||||
|
||||
// Alternative to denoise_object
|
||||
// Just use visushrink and sureshrink functions
|
||||
//visushrink(inp,N,J,wname,method,ext,thresh,level,oup);
|
||||
//sureshrink(inp,N,J,wname,method,ext,thresh,level,oup);
|
||||
|
||||
//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);
|
||||
return 0;
|
||||
}
|
||||
61
test/dwpttest.c
Normal file
61
test/dwpttest.c
Normal file
@@ -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;
|
||||
}
|
||||
2048
test/noisybumps.txt
Normal file
2048
test/noisybumps.txt
Normal file
File diff suppressed because it is too large
Load Diff
512
test/noisyheavisine.txt
Normal file
512
test/noisyheavisine.txt
Normal file
@@ -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
|
||||
1024
test/pieceregular1024.txt
Normal file
1024
test/pieceregular1024.txt
Normal file
File diff suppressed because it is too large
Load Diff
79926
test/s1.txt
Normal file
79926
test/s1.txt
Normal file
File diff suppressed because it is too large
Load Diff
504
test/sst_nino3.dat
Executable file
504
test/sst_nino3.dat
Executable file
@@ -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
|
||||
47
test/wtreetest.c
Normal file
47
test/wtreetest.c
Normal file
@@ -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;
|
||||
}
|
||||
3
unitTests/CMakeLists.txt
Normal file
3
unitTests/CMakeLists.txt
Normal file
@@ -0,0 +1,3 @@
|
||||
|
||||
add_subdirectory(wavelibBoostTests)
|
||||
|
||||
5
unitTests/wavelibBoostTests/BoostTest.cpp
Normal file
5
unitTests/wavelibBoostTests/BoostTest.cpp
Normal file
@@ -0,0 +1,5 @@
|
||||
|
||||
|
||||
#define BOOST_TEST_MODULE WaveLibTests
|
||||
|
||||
#include <boost/test/included/unit_test.hpp>
|
||||
10
unitTests/wavelibBoostTests/BoostTest.h
Normal file
10
unitTests/wavelibBoostTests/BoostTest.h
Normal file
@@ -0,0 +1,10 @@
|
||||
|
||||
|
||||
#ifndef SWALLOWING_BOOSTTEST_H_
|
||||
#define SWALLOWING_BOOSTTEST_H_
|
||||
|
||||
// we use the header only version of boost unit test
|
||||
#define BOOST_TEST_NO_LIB
|
||||
#include <boost/test/unit_test.hpp>
|
||||
|
||||
#endif // SWALLOWING_BOOSTTEST_H_
|
||||
27
unitTests/wavelibBoostTests/CMakeLists.txt
Normal file
27
unitTests/wavelibBoostTests/CMakeLists.txt
Normal file
@@ -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_SOURCE_DIR}/../../header
|
||||
)
|
||||
|
||||
|
||||
install(TARGETS wavelibLibTests
|
||||
RUNTIME DESTINATION bin
|
||||
LIBRARY DESTINATION tests
|
||||
ARCHIVE DESTINATION tests
|
||||
)
|
||||
661
unitTests/wavelibBoostTests/tst_dwt.cpp
Normal file
661
unitTests/wavelibBoostTests/tst_dwt.cpp
Normal file
@@ -0,0 +1,661 @@
|
||||
/*
|
||||
* Copyright (c) 2016 Holger Nahrstaedt (TU Berlin)
|
||||
*/
|
||||
|
||||
#include <sstream>
|
||||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
|
||||
#include <string>
|
||||
#include <cmath>
|
||||
#include "wavelib.h"
|
||||
|
||||
#include<vector>
|
||||
|
||||
namespace patch
|
||||
{
|
||||
template < typename T > std::string to_string( const T& n )
|
||||
{
|
||||
std::ostringstream stm ;
|
||||
stm << n ;
|
||||
return stm.str() ;
|
||||
}
|
||||
}
|
||||
|
||||
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 sum1(double *array, int N) {
|
||||
double sum;
|
||||
int i;
|
||||
|
||||
sum = 0.0;
|
||||
for (i = 0; i < N; ++i) {
|
||||
sum += array[i];
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
double sum2(double *array, int N) {
|
||||
double sum;
|
||||
int i;
|
||||
|
||||
sum = 0.0;
|
||||
for (i = 0; i < N; i+=2) {
|
||||
sum += array[i];
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
double sum3(double *array, int N) {
|
||||
double sum;
|
||||
int i;
|
||||
|
||||
sum = 0.0;
|
||||
for (i = 1; i < N; i += 2) {
|
||||
sum += array[i];
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
//np.sum(w[2*m:(2*N+2*m)]*w[0:2*N])
|
||||
double sum4(double *array, int N) {
|
||||
double sum;
|
||||
int i;
|
||||
|
||||
sum = 0.0;
|
||||
for (i = 0; i < N; i += 1) {
|
||||
sum += array[i] * array[i];
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
//np.sum(w[2 * m:(2 * N)] * w[0:2 * N - 2 * m])
|
||||
double sum5(double *array, int N,int m) {
|
||||
double sum;
|
||||
int i;
|
||||
|
||||
sum = 0.0;
|
||||
for (i = 2*m; i < N; i += 1) {
|
||||
sum += array[i] * array[i-2*m];
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
double RMS_Error(double *data, double *rec, int N) {
|
||||
int i;
|
||||
double sum = 0;
|
||||
for (i = 0; i < N; ++i) {
|
||||
sum += (data[i] - rec[i])*(data[i] - rec[i]);
|
||||
}
|
||||
return sqrt(sum/((double)N-1));
|
||||
}
|
||||
|
||||
double REL_Error(double *data, double *rec, int N) {
|
||||
int i;
|
||||
double sum1 = 0;
|
||||
double sum2 = 0;
|
||||
for (i = 0; i < N; ++i) {
|
||||
sum1 += (data[i] - rec[i])*(data[i] - rec[i]);
|
||||
sum2 += data[i] * data[i];
|
||||
}
|
||||
return sqrt(sum1)/sqrt(sum2);
|
||||
}
|
||||
|
||||
void ReconstructionTest()
|
||||
{
|
||||
|
||||
wave_object obj;
|
||||
wt_object wt;
|
||||
double *inp,*out;
|
||||
int N, i,J;
|
||||
double epsilon = 1e-15;
|
||||
char *type = (char*) "dwt";
|
||||
|
||||
N = 79926;
|
||||
|
||||
//N = 256;
|
||||
|
||||
inp = (double*)malloc(sizeof(double)* N);
|
||||
out = (double*)malloc(sizeof(double)* N);
|
||||
//wmean = mean(temp, N);
|
||||
|
||||
for (i = 0; i < N; ++i) {
|
||||
inp[i] = (rand() / (double)(RAND_MAX));
|
||||
}
|
||||
std::vector<std::string > waveletNames;
|
||||
|
||||
for (unsigned int j = 0; j < 36; j++)
|
||||
{
|
||||
waveletNames.push_back(std::string("db") + patch::to_string(j + 1));
|
||||
}
|
||||
for (unsigned int j = 0; j < 17; j++)
|
||||
{
|
||||
waveletNames.push_back(std::string("coif") + patch::to_string(j + 1));
|
||||
}
|
||||
for (unsigned int j = 1; j < 20; j++)
|
||||
{
|
||||
waveletNames.push_back(std::string("sym") + patch::to_string(j + 1));
|
||||
}
|
||||
|
||||
waveletNames.push_back("bior1.1");
|
||||
waveletNames.push_back("bior1.3");
|
||||
waveletNames.push_back("bior1.5");
|
||||
waveletNames.push_back("bior2.2");
|
||||
waveletNames.push_back("bior2.4");
|
||||
waveletNames.push_back("bior2.6");
|
||||
waveletNames.push_back("bior2.8");
|
||||
waveletNames.push_back("bior3.1");
|
||||
waveletNames.push_back("bior3.3");
|
||||
waveletNames.push_back("bior3.5");
|
||||
waveletNames.push_back("bior3.7");
|
||||
waveletNames.push_back("bior3.9");
|
||||
waveletNames.push_back("bior4.4");
|
||||
waveletNames.push_back("bior5.5");
|
||||
waveletNames.push_back("bior6.8");
|
||||
|
||||
waveletNames.push_back("rbior1.1");
|
||||
waveletNames.push_back("rbior1.3");
|
||||
waveletNames.push_back("rbior1.5");
|
||||
waveletNames.push_back("rbior2.2");
|
||||
waveletNames.push_back("rbior2.4");
|
||||
waveletNames.push_back("rbior2.6");
|
||||
waveletNames.push_back("rbior2.8");
|
||||
waveletNames.push_back("rbior3.1");
|
||||
waveletNames.push_back("rbior3.3");
|
||||
waveletNames.push_back("rbior3.5");
|
||||
waveletNames.push_back("rbior3.7");
|
||||
waveletNames.push_back("rbior3.9");
|
||||
waveletNames.push_back("rbior4.4");
|
||||
waveletNames.push_back("rbior5.5");
|
||||
waveletNames.push_back("rbior6.8");
|
||||
|
||||
for (unsigned int direct_fft = 0; direct_fft < 2; direct_fft++)
|
||||
{
|
||||
for (unsigned int sym_per = 0; sym_per < 2; sym_per++)
|
||||
{
|
||||
for (unsigned int j = 0; j < waveletNames.size(); j++)
|
||||
{
|
||||
char * name = new char[waveletNames[j].size() + 1];
|
||||
memcpy(name, waveletNames[j].c_str(), waveletNames[j].size() + 1);
|
||||
obj = wave_init(name);// Initialize the wavelet
|
||||
for (J = 1; J < 3; J++)
|
||||
{
|
||||
//J = 3;
|
||||
|
||||
wt = wt_init(obj,(char*) "dwt", N, J);// Initialize the wavelet transform object
|
||||
if (sym_per == 0)
|
||||
setDWTExtension(wt, (char*) "sym");// Options are "per" and "sym". Symmetric is the default option
|
||||
else
|
||||
setDWTExtension(wt, (char*) "per");
|
||||
if (direct_fft == 0)
|
||||
setWTConv(wt, (char*) "direct");
|
||||
else
|
||||
setWTConv(wt, (char*) "fft");
|
||||
|
||||
dwt(wt, inp);// Perform DWT
|
||||
|
||||
idwt(wt, out);// Perform IDWT (if needed)
|
||||
// Test Reconstruction
|
||||
|
||||
if (direct_fft == 0)
|
||||
epsilon = 1e-8;
|
||||
else
|
||||
epsilon = 1e-10;
|
||||
//BOOST_CHECK_SMALL(RMS_Error(out, inp, wt->siglength), epsilon); // If Reconstruction succeeded then the output should be a small value.
|
||||
|
||||
//printf("%g ",RMS_Error(out, inp, wt->siglength));
|
||||
if (RMS_Error(out, inp, wt->siglength) > epsilon) {
|
||||
printf("\n ERROR : DWT Reconstruction Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
wt_free(wt);
|
||||
}
|
||||
wave_free(obj);
|
||||
delete[] name;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
free(out);
|
||||
free(inp);
|
||||
}
|
||||
|
||||
void DWPTReconstructionTest()
|
||||
{
|
||||
|
||||
wave_object obj;
|
||||
wpt_object wt;
|
||||
double *inp,*out;
|
||||
int N, i,J;
|
||||
double epsilon = 1e-8;
|
||||
|
||||
N = 79926;
|
||||
|
||||
//N = 256;
|
||||
|
||||
inp = (double*)malloc(sizeof(double)* N);
|
||||
out = (double*)malloc(sizeof(double)* N);
|
||||
//wmean = mean(temp, N);
|
||||
|
||||
for (i = 0; i < N; ++i) {
|
||||
inp[i] = (rand() / (double)(RAND_MAX));
|
||||
}
|
||||
std::vector<std::string > waveletNames;
|
||||
|
||||
for (unsigned int j = 0; j < 36; j++)
|
||||
{
|
||||
waveletNames.push_back(std::string("db") + patch::to_string(j + 1));
|
||||
}
|
||||
for (unsigned int j = 0; j < 17; j++)
|
||||
{
|
||||
waveletNames.push_back(std::string("coif") + patch::to_string(j + 1));
|
||||
}
|
||||
for (unsigned int j = 1; j < 20; j++)
|
||||
{
|
||||
waveletNames.push_back(std::string("sym") + patch::to_string(j + 1));
|
||||
}
|
||||
|
||||
waveletNames.push_back("bior1.1");
|
||||
waveletNames.push_back("bior1.3");
|
||||
waveletNames.push_back("bior1.5");
|
||||
waveletNames.push_back("bior2.2");
|
||||
waveletNames.push_back("bior2.4");
|
||||
waveletNames.push_back("bior2.6");
|
||||
waveletNames.push_back("bior2.8");
|
||||
waveletNames.push_back("bior3.1");
|
||||
waveletNames.push_back("bior3.3");
|
||||
waveletNames.push_back("bior3.5");
|
||||
waveletNames.push_back("bior3.7");
|
||||
waveletNames.push_back("bior3.9");
|
||||
waveletNames.push_back("bior4.4");
|
||||
waveletNames.push_back("bior5.5");
|
||||
waveletNames.push_back("bior6.8");
|
||||
|
||||
waveletNames.push_back("rbior1.1");
|
||||
waveletNames.push_back("rbior1.3");
|
||||
waveletNames.push_back("rbior1.5");
|
||||
waveletNames.push_back("rbior2.2");
|
||||
waveletNames.push_back("rbior2.4");
|
||||
waveletNames.push_back("rbior2.6");
|
||||
waveletNames.push_back("rbior2.8");
|
||||
waveletNames.push_back("rbior3.1");
|
||||
waveletNames.push_back("rbior3.3");
|
||||
waveletNames.push_back("rbior3.5");
|
||||
waveletNames.push_back("rbior3.7");
|
||||
waveletNames.push_back("rbior3.9");
|
||||
waveletNames.push_back("rbior4.4");
|
||||
waveletNames.push_back("rbior5.5");
|
||||
waveletNames.push_back("rbior6.8");
|
||||
|
||||
for (unsigned int ent = 0; ent < 2; ent++)
|
||||
{
|
||||
for (unsigned int sym_per = 0; sym_per < 2; sym_per++)
|
||||
{
|
||||
for (unsigned int j = 0; j < waveletNames.size(); j++)
|
||||
{
|
||||
char * name = new char[waveletNames[j].size() + 1];
|
||||
memcpy(name, waveletNames[j].c_str(), waveletNames[j].size() + 1);
|
||||
obj = wave_init(name);// Initialize the wavelet
|
||||
for (J = 1; J < 3; J++)
|
||||
{
|
||||
//J = 3;
|
||||
|
||||
wt = wpt_init(obj, N, J);// Initialize the wavelet transform object
|
||||
if (sym_per == 0)
|
||||
setDWPTExtension(wt, (char*) "sym");// Options are "per" and "sym". Symmetric is the default option
|
||||
else
|
||||
setDWPTExtension(wt, (char*) "per");
|
||||
|
||||
if (ent == 0)
|
||||
setDWPTEntropy(wt, (char*) "shannon",0);
|
||||
else
|
||||
setDWPTEntropy(wt, (char*) "logenergy",0);
|
||||
|
||||
dwpt(wt, inp);// Perform DWT
|
||||
|
||||
idwpt(wt, out);// Perform IDWT (if needed)
|
||||
// Test Reconstruction
|
||||
|
||||
|
||||
//BOOST_CHECK_SMALL(RMS_Error(out, inp, wt->siglength), epsilon); // If Reconstruction succeeded then the output should be a small value.
|
||||
|
||||
//printf("%s %g \n",name,RMS_Error(out, inp, wt->siglength));
|
||||
|
||||
if (RMS_Error(out, inp, wt->siglength) > epsilon) {
|
||||
printf("\n ERROR : DWPT Reconstruction Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
wpt_free(wt);
|
||||
}
|
||||
wave_free(obj);
|
||||
delete[] name;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
free(out);
|
||||
free(inp);
|
||||
}
|
||||
|
||||
void CWTReconstructionTest() {
|
||||
int i, j,N, J,subscale,a0,iter;
|
||||
double *inp,*oup;
|
||||
double dt, dj,s0, pi,t;
|
||||
double val, epsilon;
|
||||
int it1,it2;
|
||||
cwt_object wt;
|
||||
|
||||
|
||||
char *wave[3];
|
||||
wave[0] = (char*) "morl";
|
||||
wave[1] =(char*) "paul";
|
||||
wave[2] = (char*) "dog";
|
||||
double param[30] = {4.5,5,5.5,6,6.5,8,10,13,17,20,
|
||||
4,5,7,8,10,12,13,14,17,20,2,4,6,8,10,12,14,16,18,20};
|
||||
char *type = (char*) "pow";
|
||||
|
||||
epsilon = 0.01;
|
||||
N = 2048;
|
||||
dt = 0.000125;
|
||||
subscale = 20;
|
||||
dj = 1.0 / (double) subscale;
|
||||
s0 = dt/32;
|
||||
J = 32 * subscale;
|
||||
a0 = 2;//power
|
||||
|
||||
|
||||
inp = (double*)malloc(sizeof(double)* N);
|
||||
oup = (double*)malloc(sizeof(double)* N);
|
||||
|
||||
pi = 4.0 * atan(1.0);
|
||||
|
||||
|
||||
for (i = 0; i < N; ++i) {
|
||||
t = dt * i;
|
||||
inp[i] = sin(2 * pi * 500 * t) + sin(2 * pi * 1000 * t) + 0.1 * sin(2 * pi * 8 * t);
|
||||
if (i == 1200 || i ==1232) {
|
||||
inp[i] += 5.0;
|
||||
}
|
||||
}
|
||||
|
||||
for(it1 = 0; it1 < 3;++it1) {
|
||||
for(it2 = 0; it2 < 10;++it2) {
|
||||
|
||||
wt = cwt_init(wave[it1], param[it1*10+it2], N,dt, J);
|
||||
|
||||
setCWTScales(wt, s0, dj, type, a0);
|
||||
|
||||
cwt(wt, inp);
|
||||
|
||||
|
||||
icwt(wt, oup);
|
||||
|
||||
|
||||
//printf("\nWavelet : %s Parameter %g Error %g \n", wave[it1],param[it1*10+it2],REL_Error(inp,oup, wt->siglength));
|
||||
if (REL_Error(inp,oup, wt->siglength) > epsilon) {
|
||||
printf("\n ERROR : DWPT Reconstruction Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
cwt_free(wt);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
free(inp);
|
||||
free(oup);
|
||||
|
||||
}
|
||||
|
||||
|
||||
void DBCoefTests()
|
||||
{
|
||||
wave_object obj;
|
||||
double epsilon = 1e-15;
|
||||
double t1,t2,t3,t4,t5;
|
||||
std::vector<std::string > waveletNames;
|
||||
waveletNames.resize(38);
|
||||
for (unsigned int i = 0; i < waveletNames.size();i++)
|
||||
{
|
||||
waveletNames[i] = std::string("db") + patch::to_string(i+1);
|
||||
}
|
||||
|
||||
for (unsigned int j = 0; j < waveletNames.size(); j++)
|
||||
{
|
||||
char * name = new char[waveletNames[j].size() + 1];
|
||||
memcpy(name, waveletNames[j].c_str(), waveletNames[j].size() + 1);
|
||||
obj = wave_init(name);// Initialize the wavelet
|
||||
t1 = sum1(obj->lpr, obj->lpr_len) - sqrt(2.0);
|
||||
t2 = sum2(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t3 = sum3(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t4 = sum4(obj->lpr, obj->lpr_len) - 1.;
|
||||
|
||||
if (fabs(t1) > epsilon || fabs(t2) > epsilon || fabs(t3) > epsilon || fabs(t4) > epsilon) {
|
||||
printf("\n ERROR : DB Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
for (int m = 1; m < (obj->lpr_len / 2) - 1;m++) {
|
||||
t5 = sum5(obj->lpr, obj->lpr_len, m);
|
||||
if (fabs(t5) > epsilon) {
|
||||
printf("\n ERROR : DB Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
wave_free(obj);
|
||||
delete[] name;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void CoifCoefTests()
|
||||
{
|
||||
wave_object obj;
|
||||
double epsilon = 1e-15;
|
||||
double t1,t2,t3,t4,t5;
|
||||
std::vector<std::string > waveletNames;
|
||||
waveletNames.resize(17);
|
||||
for (unsigned int i = 0; i < waveletNames.size(); i++)
|
||||
{
|
||||
waveletNames[i] = std::string("coif") + patch::to_string(i + 1);
|
||||
}
|
||||
|
||||
for (unsigned int j = 0; j < waveletNames.size(); j++)
|
||||
{
|
||||
char * name = new char[waveletNames[j].size() + 1];
|
||||
memcpy(name, waveletNames[j].c_str(), waveletNames[j].size() + 1);
|
||||
obj = wave_init(name);// Initialize the wavelet
|
||||
t1 = sum1(obj->lpr, obj->lpr_len) - sqrt(2.0);
|
||||
t2 = sum2(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t3 = sum3(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t4 = sum4(obj->lpr, obj->lpr_len) - 1.;
|
||||
|
||||
if (fabs(t1) > epsilon || fabs(t2) > epsilon || fabs(t3) > epsilon || fabs(t4) > epsilon) {
|
||||
printf("\n ERROR : Coif Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
for (int m = 1; m < (obj->lpr_len / 2) - 1;m++) {
|
||||
t5 = sum5(obj->lpr, obj->lpr_len, m);
|
||||
if (fabs(t5) > epsilon) {
|
||||
printf("\n ERROR : Coif Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
wave_free(obj);
|
||||
delete[] name;
|
||||
}
|
||||
}
|
||||
|
||||
void SymCoefTests()
|
||||
{
|
||||
wave_object obj;
|
||||
double epsilon = 1e-10;
|
||||
double t1,t2,t3,t4,t5;
|
||||
std::vector<std::string > waveletNames;
|
||||
for (unsigned int i = 1; i < 20; i++)
|
||||
{
|
||||
waveletNames.push_back(std::string("sym") + patch::to_string(i + 1));
|
||||
}
|
||||
|
||||
for (unsigned int j = 0; j < waveletNames.size(); j++)
|
||||
{
|
||||
char * name = new char[waveletNames[j].size() + 1];
|
||||
memcpy(name, waveletNames[j].c_str(), waveletNames[j].size() + 1);
|
||||
obj = wave_init(name);// Initialize the wavelet
|
||||
t1 = sum1(obj->lpr, obj->lpr_len) - sqrt(2.0);
|
||||
t2 = sum2(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t3 = sum3(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t4 = sum4(obj->lpr, obj->lpr_len) - 1.;
|
||||
|
||||
if (fabs(t1) > epsilon || fabs(t2) > epsilon || fabs(t3) > epsilon || fabs(t4) > epsilon) {
|
||||
printf("\n ERROR : Sym Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
for (int m = 1; m < (obj->lpr_len / 2) - 1;m++) {
|
||||
t5 = sum5(obj->lpr, obj->lpr_len, m);
|
||||
if (fabs(t5) > epsilon) {
|
||||
printf("\n ERROR : Sym Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
wave_free(obj);
|
||||
delete[] name;
|
||||
}
|
||||
}
|
||||
|
||||
void BiorCoefTests()
|
||||
{
|
||||
wave_object obj;
|
||||
double epsilon = 1e-10;
|
||||
double t1,t2,t3,t4,t5,t6;
|
||||
std::vector<std::string > waveletNames;
|
||||
waveletNames.push_back("bior1.1");
|
||||
waveletNames.push_back("bior1.3");
|
||||
waveletNames.push_back("bior1.5");
|
||||
waveletNames.push_back("bior2.2");
|
||||
waveletNames.push_back("bior2.4");
|
||||
waveletNames.push_back("bior2.6");
|
||||
waveletNames.push_back("bior2.8");
|
||||
waveletNames.push_back("bior3.1");
|
||||
waveletNames.push_back("bior3.3");
|
||||
waveletNames.push_back("bior3.5");
|
||||
waveletNames.push_back("bior3.7");
|
||||
waveletNames.push_back("bior3.9");
|
||||
waveletNames.push_back("bior4.4");
|
||||
waveletNames.push_back("bior5.5");
|
||||
waveletNames.push_back("bior6.8");
|
||||
|
||||
for (unsigned int j = 0; j < waveletNames.size(); j++)
|
||||
{
|
||||
char * name = new char[waveletNames[j].size() + 1];
|
||||
memcpy(name, waveletNames[j].c_str(), waveletNames[j].size() + 1);
|
||||
obj = wave_init(name);// Initialize the wavelet
|
||||
|
||||
t1 = sum1(obj->lpr, obj->lpr_len) - sqrt(2.0);
|
||||
t2 = sum1(obj->lpd, obj->lpd_len) - sqrt(2.0);
|
||||
|
||||
t3 = sum2(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t4 = sum2(obj->lpd, obj->lpd_len) - 1. / sqrt(2.0);
|
||||
|
||||
t5 = sum3(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t6 = sum3(obj->lpd, obj->lpd_len) - 1. / sqrt(2.0);
|
||||
|
||||
if (fabs(t1) > epsilon || fabs(t2) > epsilon || fabs(t3) > epsilon || fabs(t4) > epsilon || fabs(t5) > epsilon || fabs(t6) > epsilon ) {
|
||||
printf("\n ERROR : Bior Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
|
||||
wave_free(obj);
|
||||
delete[] name;
|
||||
}
|
||||
}
|
||||
|
||||
void RBiorCoefTests()
|
||||
{
|
||||
wave_object obj;
|
||||
double epsilon = 1e-10;
|
||||
double t1,t2,t3,t4,t5,t6;
|
||||
std::vector<std::string > waveletNames;
|
||||
waveletNames.push_back("rbior1.1");
|
||||
waveletNames.push_back("rbior1.3");
|
||||
waveletNames.push_back("rbior1.5");
|
||||
waveletNames.push_back("rbior2.2");
|
||||
waveletNames.push_back("rbior2.4");
|
||||
waveletNames.push_back("rbior2.6");
|
||||
waveletNames.push_back("rbior2.8");
|
||||
waveletNames.push_back("rbior3.1");
|
||||
waveletNames.push_back("rbior3.3");
|
||||
waveletNames.push_back("rbior3.5");
|
||||
waveletNames.push_back("rbior3.7");
|
||||
waveletNames.push_back("rbior3.9");
|
||||
waveletNames.push_back("rbior4.4");
|
||||
waveletNames.push_back("rbior5.5");
|
||||
waveletNames.push_back("rbior6.8");
|
||||
|
||||
for (unsigned int j = 0; j < waveletNames.size(); j++)
|
||||
{
|
||||
char * name = new char[waveletNames[j].size() + 1];
|
||||
memcpy(name, waveletNames[j].c_str(), waveletNames[j].size() + 1);
|
||||
obj = wave_init(name);// Initialize the wavelet
|
||||
|
||||
t1 = sum1(obj->lpr, obj->lpr_len) - sqrt(2.0);
|
||||
t2 = sum1(obj->lpd, obj->lpd_len) - sqrt(2.0);
|
||||
|
||||
t3 = sum2(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t4 = sum2(obj->lpd, obj->lpd_len) - 1. / sqrt(2.0);
|
||||
|
||||
t5 = sum3(obj->lpr, obj->lpr_len) - 1. / sqrt(2.0);
|
||||
t6 = sum3(obj->lpd, obj->lpd_len) - 1. / sqrt(2.0);
|
||||
|
||||
if (fabs(t1) > epsilon || fabs(t2) > epsilon || fabs(t3) > epsilon || fabs(t4) > epsilon || fabs(t5) > epsilon || fabs(t6) > epsilon ) {
|
||||
printf("\n ERROR : RBior Coefficients Unit Test Failed. Exiting. \n");
|
||||
exit(-1);
|
||||
}
|
||||
wave_free(obj);
|
||||
delete[] name;
|
||||
}
|
||||
}
|
||||
|
||||
int main() {
|
||||
printf("Running Unit Tests : \n \n");
|
||||
printf("Running DBCoefTests ... ");
|
||||
DBCoefTests();
|
||||
printf("DONE \n");
|
||||
printf("Running CoifCoefTests ... ");
|
||||
CoifCoefTests();
|
||||
printf("DONE \n");
|
||||
printf("Running SymCoefTests ... ");
|
||||
SymCoefTests();
|
||||
printf("DONE \n");
|
||||
printf("Running BiorCoefTests ... ");
|
||||
BiorCoefTests();
|
||||
printf("DONE \n");
|
||||
printf("Running RBiorCoefTests ... ");
|
||||
RBiorCoefTests();
|
||||
printf("DONE \n");
|
||||
printf("Running DWT ReconstructionTests ... ");
|
||||
ReconstructionTest();
|
||||
printf("DONE \n");
|
||||
printf("Running DWPT ReconstructionTests ... ");
|
||||
DWPTReconstructionTest();
|
||||
printf("DONE \n");
|
||||
printf("Running CWT ReconstructionTests ... ");
|
||||
CWTReconstructionTest();
|
||||
printf("DONE \n");
|
||||
printf("\n\nUnit Tests Successful\n\n");
|
||||
return 0;
|
||||
|
||||
}
|
||||
|
||||
BIN
wavelib-doc.pdf
Normal file
BIN
wavelib-doc.pdf
Normal file
Binary file not shown.
Reference in New Issue
Block a user