diff --git a/auxiliary/CMakeLists.txt b/auxiliary/CMakeLists.txt index ba87b7b..51c836d 100644 --- a/auxiliary/CMakeLists.txt +++ b/auxiliary/CMakeLists.txt @@ -4,9 +4,7 @@ set(SOURCE_FILES denoise.c waux.c ) -set(HEADER_FILES denoise.h - waux.h -) +set(HEADER_FILES waux.h) add_library(wauxlib STATIC ${SOURCE_FILES} ${HEADER_FILES}) diff --git a/auxiliary/denoise.c b/auxiliary/denoise.c index 1cd10a7..5354252 100644 --- a/auxiliary/denoise.c +++ b/auxiliary/denoise.c @@ -1,7 +1,12 @@ +#include +#include +#include +#include -#include "denoise.h" +#include "waux.h" +#include "wauxlib.h" -denoise_object denoise_init(int length, int J,char* wname) { +denoise_object denoise_init(int length, int J,const char* wname) { denoise_object obj = NULL; obj = (denoise_object)malloc(sizeof(struct denoise_set) +sizeof(double)); @@ -21,7 +26,7 @@ denoise_object denoise_init(int length, int J,char* wname) { return obj; } -void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised) { +void visushrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised) { int filt_len,iter,i,dlen,dwt_len,sgn, MaxIter,it; double sigma,td,tmp; wave_object wave; @@ -128,7 +133,7 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch wt_free(wt); } -void sureshrink(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,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised) { int filt_len,i,it,len,dlen,dwt_len,min_index,sgn, MaxIter,iter; double sigma,norm,td,tv,te,ct,thr,temp,x_sum; wave_object wave; @@ -286,7 +291,7 @@ void denoise(denoise_object obj, double *signal,double *denoised) { } } -void setDenoiseMethod(denoise_object obj, char *dmethod) { +void setDenoiseMethod(denoise_object obj, const char *dmethod) { if (!strcmp(dmethod, "sureshrink")) { strcpy(obj->dmethod, "sureshrink"); } @@ -299,7 +304,7 @@ void setDenoiseMethod(denoise_object obj, char *dmethod) { } } -void setDenoiseWTMethod(denoise_object obj, char *wmethod) { +void setDenoiseWTMethod(denoise_object obj, const char *wmethod) { if (!strcmp(wmethod, "dwt")) { strcpy(obj->wmethod, "dwt"); } @@ -312,7 +317,7 @@ void setDenoiseWTMethod(denoise_object obj, char *wmethod) { } } -void setDenoiseWTExtension(denoise_object obj, char *extension) { +void setDenoiseWTExtension(denoise_object obj, const char *extension) { if (!strcmp(extension, "sym")) { strcpy(obj->ext, "sym"); } @@ -325,7 +330,7 @@ void setDenoiseWTExtension(denoise_object obj, char *extension) { } } -void setDenoiseParameters(denoise_object obj, char *thresh,char *level) { +void setDenoiseParameters(denoise_object obj, const char *thresh,const char *level) { //Set thresholding if (!strcmp(thresh, "soft")) { diff --git a/auxiliary/denoise.h b/auxiliary/denoise.h deleted file mode 100644 index 11e0e99..0000000 --- a/auxiliary/denoise.h +++ /dev/null @@ -1,52 +0,0 @@ -/* -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_ */ diff --git a/auxiliary/waux.c b/auxiliary/waux.c index 10ed389..1b22143 100644 --- a/auxiliary/waux.c +++ b/auxiliary/waux.c @@ -1,3 +1,4 @@ +#include "wauxlib.h" #include "waux.h" int compare_double(const void* a, const void* b) @@ -11,7 +12,7 @@ int compare_double(const void* a, const void* b) } -double mean(double* vec, int N) { +double mean(const double* vec, int N) { int i; double m; m = 0.0; @@ -23,7 +24,7 @@ double mean(double* vec, int N) { return m; } -double var(double* vec, int N) { +double var(const double* vec, int N) { double v,temp,m; int i; v = 0.0; @@ -69,7 +70,7 @@ double mad(double *x, int N) { return sigma; } -int minindex(double *arr, int N) { +int minindex(const double *arr, int N) { double min; int index,i; @@ -130,7 +131,7 @@ 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, +void getDWTRecCoeff(double *coeff,int *length,const char *ctype,const char *ext, int level, int J,double *lpr, double *hpr,int lf,int siglength,double *reccoeff) { int i,j,k,det_len,N,l,m,n,v,t,l2; @@ -274,7 +275,7 @@ void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, } -void autocovar(double* vec,int N, double* acov,int M) { +void autocovar(const double* vec,int N, double* acov,int M) { double m,temp1,temp2; int i,t; m = mean(vec,N); @@ -301,7 +302,7 @@ void autocovar(double* vec,int N, double* acov,int M) { } -void autocorr(double* vec,int N,double* acorr, int M) { +void autocorr(const double* vec,int N,double* acorr, int M) { double var; int i; if (M > N) { diff --git a/auxiliary/waux.h b/auxiliary/waux.h index d0a3dbb..f4d99d7 100644 --- a/auxiliary/waux.h +++ b/auxiliary/waux.h @@ -21,26 +21,21 @@ extern "C" { int compare_double(const void* a, const void* b); -double mean(double* vec, int N); +double mean(const double* vec, int N); -double var(double* vec, int N); +double var(const double* vec, int N); double median(double *x, int N); -double mad(double *x, int N); - -int minindex(double *arr, int N); +int minindex(const double *arr, int N); void getDWTAppx(wt_object wt, double *appx,int N); void getDWTDetail(wt_object wt, double *detail, int N, int level); -void 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(const double* vec,int N,double* acov, int M); -void autocovar(double* vec,int N,double* acov, int M); - -void autocorr(double* vec,int N,double* acorr, int M); +void autocorr(const double* vec,int N,double* acorr, int M); #ifdef __cplusplus diff --git a/header/wauxlib.h b/header/wauxlib.h index 85ca443..53e4647 100644 --- a/header/wauxlib.h +++ b/header/wauxlib.h @@ -12,7 +12,7 @@ extern "C" { typedef struct denoise_set* denoise_object; -denoise_object denoise_init(int length, int J,char* wname); +denoise_object denoise_init(int length, int J,const char* wname); struct denoise_set{ int N; //signal length @@ -26,23 +26,23 @@ struct denoise_set{ //double params[0]; }; -void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); +void visushrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised); -void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); +void sureshrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised); void denoise(denoise_object obj, double *signal,double *denoised); -void setDenoiseMethod(denoise_object obj, char *dmethod); +void setDenoiseMethod(denoise_object obj, const char *dmethod); -void setDenoiseWTMethod(denoise_object obj, char *wmethod); +void setDenoiseWTMethod(denoise_object obj, const char *wmethod); -void setDenoiseWTExtension(denoise_object obj, char *extension); +void setDenoiseWTExtension(denoise_object obj, const char *extension); -void setDenoiseParameters(denoise_object obj, char *thresh,char *level); +void setDenoiseParameters(denoise_object obj, const char *thresh,const char *level); void denoise_free(denoise_object object); -void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, int J,double *lpr, +void getDWTRecCoeff(double *coeff,int *length,const char *ctype,const char *ext, int level, int J,double *lpr, double *hpr,int lf,int siglength,double *reccoeff); double mad(double *x, int N); diff --git a/header/wavelib.h b/header/wavelib.h index 2399d8c..e36096b 100644 --- a/header/wavelib.h +++ b/header/wavelib.h @@ -26,7 +26,7 @@ typedef struct cplx_t { typedef struct wave_set* wave_object; -wave_object wave_init(char* wname); +wave_object wave_init(const char* wname); struct wave_set{ char wname[50]; @@ -83,7 +83,7 @@ struct conv_set{ typedef struct wt_set* wt_object; -wt_object wt_init(wave_object wave,char* method, int siglength, int J); +wt_object wt_init(wave_object wave,const char* method, int siglength, int J); struct wt_set{ wave_object wave; @@ -165,7 +165,7 @@ struct wpt_set{ typedef struct cwt_set* cwt_object; -cwt_object cwt_init(char* wave, double param, int siglength,double dt, int J); +cwt_object cwt_init(const char* wave, double param, int siglength,double dt, int J); struct cwt_set{ char wave[10];// Wavelet - morl/morlet,paul,dog/dgauss @@ -191,33 +191,33 @@ struct cwt_set{ }; -void dwt(wt_object wt, double *inp); +void dwt(wt_object wt, const double *inp); void idwt(wt_object wt, double *dwtop); -void wtree(wtree_object wt, double *inp); +void wtree(wtree_object wt, const double *inp); -void dwpt(wpt_object wt, double *inp); +void dwpt(wpt_object wt, const double *inp); void idwpt(wpt_object wt, double *dwtop); -void swt(wt_object wt, double *inp); +void swt(wt_object wt, const double *inp); void iswt(wt_object wt, double *swtop); -void modwt(wt_object wt, double *inp); +void modwt(wt_object wt, const double *inp); void imodwt(wt_object wt, double *dwtop); -void setDWTExtension(wt_object wt, char *extension); +void setDWTExtension(wt_object wt, const char *extension); -void setWTREEExtension(wtree_object wt, char *extension); +void setWTREEExtension(wtree_object wt, const char *extension); -void setDWPTExtension(wpt_object wt, char *extension); +void setDWPTExtension(wpt_object wt, const char *extension); -void setDWPTEntropy(wpt_object wt, char *entropy, double eparam); +void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam); -void setWTConv(wt_object wt, char *cmethod); +void setWTConv(wt_object wt, const char *cmethod); int getWTREENodelength(wtree_object wt, int X); @@ -227,13 +227,13 @@ 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 setCWTScales(cwt_object wt, double s0, double dj, const char *type, int power); -void setCWTScaleVector(cwt_object wt, double *scale, int J, double s0, double dj); +void setCWTScaleVector(cwt_object wt, const double *scale, int J, double s0, double dj); void setCWTPadding(cwt_object wt, int pad); -void cwt(cwt_object wt, double *inp); +void cwt(cwt_object wt, const double *inp); void icwt(cwt_object wt, double *cwtop); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 758f5c3..2781764 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -20,7 +20,6 @@ set(HEADER_FILES conv.h real.h wavefilt.h wavefunc.h - wavelib.h wtmath.h ) @@ -28,5 +27,5 @@ add_library(wavelib STATIC ${SOURCE_FILES} ${HEADER_FILES}) set_property(TARGET wavelib PROPERTY FOLDER "lib") -target_include_directories(wavelib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +target_include_directories(wavelib PUBLIC ${CMAKE_SOURCE_DIR}/header) diff --git a/src/conv.h b/src/conv.h index 40b2548..ffeba29 100644 --- a/src/conv.h +++ b/src/conv.h @@ -17,18 +17,8 @@ extern "C" { #define MIN(a,b) (((a)<(b))?(a):(b)) #define MAX(a,b) (((a)>(b))?(a):(b)) -typedef struct conv_set* conv_object; - conv_object conv_init(int N, int L); -struct conv_set{ - fft_real_object fobj; - fft_real_object iobj; - int ilen1; - int ilen2; - int clen; -}; - int factorf(int M); int findnext(int M); diff --git a/src/cwt.c b/src/cwt.c index 6e9efb1..1695e5a 100755 --- a/src/cwt.c +++ b/src/cwt.c @@ -8,37 +8,8 @@ C. Torrence and G. Compo, and is available at URL: http://atoc.colorado.edu/rese #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, + static const double fact[41] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000, 20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000, 51090942171709440000.0, 1124000727777607680000.0, 25852016738884976640000.0, 620448401733239439360000.0, 15511210043330985984000000.0, 403291461126605635584000000.0, 10888869450418352160768000000.0, 304888344611713860501504000000.0, 8841761993739701954543616000000.0, 265252859812191058636308480000000.0, 8222838654177922817725562880000000.0, @@ -46,6 +17,11 @@ double factorial(int N) { 371993326789901217467999448150835200000000.0, 13763753091226345046315979581580902400000000.0, 523022617466601111760007224100074291200000000.0, 20397882081197443358640281739902897356800000000.0, 815915283247897734345611269596115894272000000000.0 }; + if (N > 40 || N < 0) { + printf("This program is only valid for 0 <= N <= 40 \n"); + return -1.0; + } + return fact[N]; } @@ -142,7 +118,7 @@ static void wave_function(int nk, double dt,int mother, double param,double scal } } -void cwavelet(double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad, +void cwavelet(const double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad, double *wave, double *scale, double *period, double *coi) { int i, j, k, iter; @@ -153,6 +129,8 @@ void cwavelet(double *y, int N, double dt, int mother, double param, double s0, fft_object obj, iobj; fft_data *ypad, *yfft,*daughter; + (void)s0; (void)dj; /* yes, we need these parameters unused */ + pi = 4.0 * atan(1.0); if (npad < N) { diff --git a/src/cwt.h b/src/cwt.h index ed9ea9a..d1e4a03 100755 --- a/src/cwt.h +++ b/src/cwt.h @@ -7,7 +7,7 @@ extern "C" { #endif -void cwavelet(double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad, +void cwavelet(const double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad, double *wave, double *scale, double *period, double *coi); void psi0(int mother, double param, double *val, int *real); diff --git a/src/hsfft.c b/src/hsfft.c index bf4c650..3919c67 100644 --- a/src/hsfft.c +++ b/src/hsfft.c @@ -18,7 +18,7 @@ fft_object fft_init(int N, int sgn) { if (out == 1) { obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (N-1)); obj->lf = factors(N,obj->factors); - longvectorN(obj->twiddle,N,obj->factors,obj->lf); + longvectorN(obj->twiddle,obj->factors,obj->lf); twi_len = N; obj->lt = 0; } else { @@ -32,7 +32,7 @@ fft_object fft_init(int N, int sgn) { } obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (M-1)); obj->lf = factors(M,obj->factors); - longvectorN(obj->twiddle,M,obj->factors,obj->lf); + longvectorN(obj->twiddle,obj->factors,obj->lf); obj->lt = 1; twi_len = M; } @@ -1831,7 +1831,7 @@ void twiddle(fft_data *vec,int N, int radix) { } -void longvectorN(fft_data *sig,int N, int *array, int tx) { +void longvectorN(fft_data *sig, int *array, int tx) { int L,i,Ls,ct,j,k; fft_type theta; L = 1; diff --git a/src/hsfft.h b/src/hsfft.h index 95d6b71..52297dd 100644 --- a/src/hsfft.h +++ b/src/hsfft.h @@ -12,6 +12,8 @@ #include #include +#include "wavelib.h" + #ifdef __cplusplus extern "C" { #endif @@ -22,11 +24,6 @@ extern "C" { #define fft_type double #endif - -typedef struct fft_t { - fft_type re; - fft_type im; -} fft_data; /* #define SADD(a,b) ((a)+(b)) @@ -35,19 +32,8 @@ typedef struct fft_t { #define SMUL(a,b) ((a)*(b)) */ -typedef struct fft_set* fft_object; - fft_object fft_init(int N, int sgn); -struct fft_set{ - int N; - int sgn; - int factors[64]; - int lf; - int lt; - fft_data twiddle[1]; -}; - void fft_exec(fft_object obj,fft_data *inp,fft_data *oup); int divideby(int M,int d); @@ -60,7 +46,7 @@ int factors(int M, int* arr); void twiddle(fft_data *sig,int N, int radix); -void longvectorN(fft_data *sig,int N, int *array, int M); +void longvectorN(fft_data *sig, int *array, int M); void free_fft(fft_object object); diff --git a/src/real.c b/src/real.c index c990f19..37f69d6 100644 --- a/src/real.c +++ b/src/real.c @@ -9,11 +9,9 @@ fft_real_object fft_real_init(int N, int sgn) { fft_real_object obj = NULL; - fft_type PI, theta; + fft_type theta; int k; - PI = 3.1415926535897932384626433832795; - obj = (fft_real_object) malloc (sizeof(struct fft_real_set) + sizeof(fft_data)* (N/2)); obj->cobj = fft_init(N/2,sgn); @@ -79,10 +77,9 @@ void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup) { fft_data* cinp; fft_data* coup; - int i,N2,N; + int i,N2; fft_type temp1,temp2; N2 = obj->cobj->N; - N = N2*2; cinp = (fft_data*) malloc (sizeof(fft_data) * N2); coup = (fft_data*) malloc (sizeof(fft_data) * N2); diff --git a/src/real.h b/src/real.h index 42d30eb..8632305 100644 --- a/src/real.h +++ b/src/real.h @@ -14,15 +14,8 @@ extern "C" { #endif -typedef struct fft_real_set* fft_real_object; - fft_real_object fft_real_init(int N, int sgn); -struct fft_real_set{ - fft_object cobj; - fft_data twiddle2[1]; -}; - void fft_r2c_exec(fft_real_object obj,fft_type *inp,fft_data *oup); void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup); diff --git a/src/wavefilt.c b/src/wavefilt.c index 53fb9c4..cf3e041 100644 --- a/src/wavefilt.c +++ b/src/wavefilt.c @@ -3311,7 +3311,6 @@ void copy_reverse(const double *in, int N,double *out) void qmf_wrev(const double *in, int N, double *out) { - int count = 0; double *sigOutTemp; sigOutTemp = (double*)malloc(N*sizeof(double)); diff --git a/src/wavelib.c b/src/wavelib.c index 2936572..aaf028c 100644 --- a/src/wavelib.c +++ b/src/wavelib.c @@ -2,9 +2,16 @@ Copyright (c) 2014, Rafat Hussain */ -#include "wavelib.h" +#include +#include +#include +#include -wave_object wave_init(char* wname) { +#include "cwt.h" +#include "wavelib.h" +#include "wtmath.h" + +wave_object wave_init(const char* wname) { wave_object obj = NULL; int retval; retval = 0; @@ -30,7 +37,7 @@ wave_object wave_init(char* wname) { return obj; } -wt_object wt_init(wave_object wave,char* method, int siglength,int J) { +wt_object wt_init(wave_object wave,const char* method, int siglength,int J) { int size,i,MaxIter; wt_object obj = NULL; @@ -258,13 +265,13 @@ wpt_object wpt_init(wave_object wave, int siglength, int J) { return obj; } -cwt_object cwt_init(char* wave, double param,int siglength, double dt, int J) { +cwt_object cwt_init(const char* wave, double param,int siglength, double dt, int J) { cwt_object obj = NULL; int N, i,nj2,ibase2,mother; double s0, dj; double t1; int m, odd; - char *pdefault = "pow"; + const char *pdefault = "pow"; m = (int)param; odd = 1; @@ -370,7 +377,7 @@ static void wconv(wt_object wt, double *sig, int N, double *filt, int L, double } -static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) { int l,l2,isodd,i,t,len_avg; len_avg = wt->wave->lpd_len; @@ -426,7 +433,7 @@ static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, do } -static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD) { int l, l2, isodd, i, t, len_avg; len_avg = wt->wave->lpd_len; @@ -482,7 +489,7 @@ static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_c } -static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) { int l, l2, isodd, i, t, len_avg; len_avg = wt->wave->lpd_len; @@ -538,7 +545,7 @@ static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, } -static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) { int i,l, t, len_avg; len_avg = wt->wave->lpd_len; @@ -567,7 +574,7 @@ static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, do } -static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD) { int i, l, t, len_avg; len_avg = wt->wave->lpd_len; @@ -596,7 +603,7 @@ static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_c } -static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) { int i, l, t, len_avg; len_avg = wt->wave->lpd_len; @@ -625,7 +632,7 @@ static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, } -static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, int len_cA, double *cD, int len_cD) { +static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, double *cD) { int len_avg,D,lf; double *signal,*cA_undec; len_avg = (wt->wave->lpd_len + wt->wave->hpd_len) / 2; @@ -700,7 +707,7 @@ static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, int len_cA, d free(cA_undec); } -void dwt(wt_object wt,double *inp) { +void dwt(wt_object wt,const double *inp) { int i,J,temp_len,iter,N,lp; int len_cA; double *orig,*orig2; @@ -753,10 +760,10 @@ void dwt(wt_object wt,double *inp) { len_cA = wt->length[J - iter]; N -= len_cA; if ( !strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT") ) { - dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + dwt1(wt, orig, temp_len, orig2, wt->params + N); } else { - dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N); } temp_len = wt->length[J - iter]; if (iter == J - 1) { @@ -789,10 +796,10 @@ void dwt(wt_object wt,double *inp) { len_cA = wt->length[J - iter]; N -= len_cA; if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) { - dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + dwt1(wt, orig, temp_len, orig2, wt->params + N); } else { - dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N); } temp_len = wt->length[J - iter]; @@ -817,7 +824,7 @@ void dwt(wt_object wt,double *inp) { free(orig2); } -void wtree(wtree_object wt,double *inp) { +void wtree(wtree_object wt,const double *inp) { int i,J,temp_len,iter,N,lp,p2,k,N2,Np; int len_cA,t,t2,it1; double *orig; @@ -871,9 +878,9 @@ void wtree(wtree_object wt,double *inp) { N = N2; for(k = 0; k < p2;++k) { if (iter == 0) { - wtree_per(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); + wtree_per(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA); } else { - wtree_per(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); + wtree_per(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA); } N += 2 * len_cA; } @@ -906,9 +913,9 @@ void wtree(wtree_object wt,double *inp) { N = N2; for(k = 0; k < p2;++k) { if (iter == 0) { - wtree_sym(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); + wtree_sym(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA); } else { - wtree_sym(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); + wtree_sym(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA); } N += 2 * len_cA; } @@ -959,7 +966,7 @@ static int ipow2(int n) { return p; } -void dwpt(wpt_object wt, double *inp) { +void dwpt(wpt_object wt, const double *inp) { int i, J, temp_len, iter, N, lp, p2, k, N2, Np; int temp, elength, temp2,size,nodes,llb,n1,j; double eparam,v1,v2; @@ -1026,10 +1033,10 @@ void dwpt(wpt_object wt, double *inp) { N = N2; for (k = 0; k < p2; ++k) { if (iter == 0) { - dwpt_per(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); + dwpt_per(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA); } else { - dwpt_per(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); + dwpt_per(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA); } wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam); it2++; @@ -1066,10 +1073,10 @@ void dwpt(wpt_object wt, double *inp) { N = N2; for (k = 0; k < p2; ++k) { if (iter == 0) { - dwpt_sym(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); + dwpt_sym(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA); } else { - dwpt_sym(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); + dwpt_sym(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA); } wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam); it2++; @@ -1327,7 +1334,7 @@ int getCWTScaleLength(int N) { return J; } -void setCWTScales(cwt_object wt, double s0, double dj,char *type,int power) { +void setCWTScales(cwt_object wt, double s0, double dj,const char *type,int power) { int i; strcpy(wt->type,type); //s0*pow(2.0, (double)(j - 1)*dj); @@ -1353,7 +1360,7 @@ void setCWTScales(cwt_object wt, double s0, double dj,char *type,int power) { wt->dj = dj; } -void setCWTScaleVector(cwt_object wt, double *scale, int J,double s0,double dj) { +void setCWTScaleVector(cwt_object wt, const double *scale, int J,double s0,double dj) { int i; if (J != wt->J) { @@ -1378,7 +1385,7 @@ void setCWTPadding(cwt_object wt, int pad) { } } -void cwt(cwt_object wt, double *inp) { +void cwt(cwt_object wt, const double *inp) { int i, N, npad,nj2,j,j2; N = wt->siglength; if (wt->sflag == 0) { @@ -1406,7 +1413,7 @@ void cwt(cwt_object wt, double *inp) { } wt->smean /= N; - cwavelet(inp, N, wt->dt, wt->mother, wt->m, wt->s0, wt->dj, wt->J,npad,wt->params, wt->params+nj2, wt->params+nj2+j, wt->params+nj2+j2); + cwavelet(inp, N, wt->dt, wt->mother, wt->m, wt->s0,wt->dj,wt->J,npad,wt->params, wt->params+nj2, wt->params+nj2+j, wt->params+nj2+j2); } @@ -1477,7 +1484,7 @@ static void idwt1(wt_object wt,double *temp, double *cA_up,double *cA, int len_c } -static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { +static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, double *X) { int len_avg,i,l,m,n,t,l2; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; @@ -1508,7 +1515,7 @@ static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_c } } -static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { +static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, double *X) { int len_avg, i, l, m, n, t, v; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; @@ -1592,7 +1599,7 @@ void idwt(wt_object wt, double *dwtop) { //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); - idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp); + idwt_per(wt,out, det_len, wt->output + iter, X_lp); for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) { out[k - lf/2 + 1] = X_lp[k]; } @@ -1621,7 +1628,7 @@ void idwt(wt_object wt, double *dwtop) { //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); - idwt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp); + idwt_sym(wt, out, det_len, wt->output + iter, X_lp); for (k = lf-2; k < 2 * det_len; ++k) { out[k - lf + 2] = X_lp[k]; } @@ -1695,7 +1702,7 @@ void idwt(wt_object wt, double *dwtop) { } -static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { +static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, double *X) { int len_avg, i, l, m, n, t, l2; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; @@ -1726,7 +1733,7 @@ static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, int len } } -static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { +static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, double *X) { int len_avg, i, l, m, n, t, v; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; @@ -1750,13 +1757,12 @@ static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, int len } void idwpt(wpt_object wt, double *dwtop) { - int J, U, i, lf, N, k,p,l; + int J, i, lf, k,p,l; int app_len, det_len, index, n1, llb, index2, index3, index4,indexp,xlen; double *X_lp, *X, *out, *out2; int *prep,*ptemp; J = wt->J; - U = 2; app_len = wt->length[0]; p = ipow2(J); lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; @@ -1772,251 +1778,250 @@ void idwpt(wpt_object wt, double *dwtop) { llb = 1; index2 = xlen / p; indexp = 0; - if (wt->basisvector[0] == 1) { - for (i = 0; i < wt->siglength; ++i) { - dwtop[i] = wt->output[i]; - } - - } - else { - for (i = 0; i < J; ++i) { - llb *= 2; - n1 += llb; - } - - for (i = 0; i < xlen; ++i) { - X[i] = 0.0; - } - - for (i = 0; i < llb; ++i) { - prep[i] = (int)wt->basisvector[n1 - llb + i]; - ptemp[i] = 0; - } - - if (!strcmp(wt->ext, "per")) { - app_len = wt->length[0]; - det_len = wt->length[1]; - index = 0; - - - for (i = 0; i < J; ++i) { - p = ipow2(J - i - 1); - det_len = wt->length[i + 1]; - index2 *= 2; - index3 = 0; - index4 = 0; - //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); - n1 -= llb; - for (l = 0; l < llb; ++l) { - if (ptemp[l] != 2) { - prep[l] = (int)wt->basisvector[n1 + l]; - } - else { - prep[l] = ptemp[l]; - } - ptemp[l] = 0; - } - - - for (l = 0; l < p; ++l) { - if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) { - for (k = 0; k < det_len; ++k) { - out[k] = wt->output[index + k]; - out2[k] = wt->output[index + det_len + k]; - } - idwpt_per(wt, out, det_len, out2, det_len, X_lp); - for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { - X[index3 + k - lf / 2 + 1] = X_lp[k]; - } - index += 2 * det_len; - index3 += index2; - index4 += 2 * indexp; - ptemp[l] = 2; - } - else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) { - index4 += indexp; - for (k = 0; k < det_len; ++k) { - out[k] = wt->output[index + k]; - out2[k] = X[index4 + k]; - } - idwpt_per(wt, out, det_len, out2, det_len, X_lp); - for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { - X[index3 + k - lf / 2 + 1] = X_lp[k]; - } - index += det_len; - index3 += index2; - index4 += indexp; - ptemp[l] = 2; - } - else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) { - for (k = 0; k < det_len; ++k) { - out[k] = X[index4 + k]; - out2[k] = wt->output[index + k]; - } - idwpt_per(wt, out, det_len, out2, det_len, X_lp); - for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { - X[index3 + k - lf / 2 + 1] = X_lp[k]; - } - index += det_len; - index3 += index2; - index4 += 2 * indexp; - ptemp[l] = 2; - } - else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) { - for (k = 0; k < det_len; ++k) { - out[k] = X[index4 + k]; - out2[k] = X[index4 + indexp + k]; - } - idwpt_per(wt, out, det_len, out2, det_len, X_lp); - for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { - X[index3 + k - lf / 2 + 1] = X_lp[k]; - } - index4 += 2 * indexp; - index3 += index2; - ptemp[l] = 2; - } - else { - index3 += index2; - index4 += 2 * indexp; - } - - } - - - /* - idwt_per(wt, out, det_len, wt->output + iter, det_len, X_lp); - for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { - out[k - lf / 2 + 1] = X_lp[k]; - } - - iter += det_len; - det_len = wt->length[i + 2]; - */ - llb /= 2; - indexp = index2; - } - - //free(X_lp); - - } - else if (!strcmp(wt->ext, "sym")) { - app_len = wt->length[0]; - det_len = wt->length[1]; - N = 2 * wt->length[J] - 1; - - //X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); - index = 0; - - for (i = 0; i < J; ++i) { - p = ipow2(J - i - 1); - det_len = wt->length[i + 1]; - index2 *= 2; - index3 = 0; - index4 = 0; - //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); - n1 -= llb; - for (l = 0; l < llb; ++l) { - if (ptemp[l] != 2) { - prep[l] = (int)wt->basisvector[n1 + l]; - } - else { - prep[l] = ptemp[l]; - } - ptemp[l] = 0; - } - - - for (l = 0; l < p; ++l) { - if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) { - for (k = 0; k < det_len; ++k) { - out[k] = wt->output[index + k]; - out2[k] = wt->output[index + det_len + k]; - } - idwpt_sym(wt, out, det_len, out2, det_len, X_lp); - for (k = lf - 2; k < 2 * det_len; ++k) { - X[index3 + k - lf + 2] = X_lp[k]; - } - index += 2 * det_len; - index3 += index2; - index4 += 2 * indexp; - ptemp[l] = 2; - } - else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) { - index4 += indexp; - for (k = 0; k < det_len; ++k) { - out[k] = wt->output[index + k]; - out2[k] = X[index4 + k]; - } - idwpt_sym(wt, out, det_len, out2, det_len, X_lp); - for (k = lf - 2; k < 2 * det_len; ++k) { - X[index3 + k - lf + 2] = X_lp[k]; - } - index += det_len; - index3 += index2; - index4 += indexp; - ptemp[l] = 2; - } - else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) { - for (k = 0; k < det_len; ++k) { - out[k] = X[index4 + k]; - out2[k] = wt->output[index + k]; - } - idwpt_sym(wt, out, det_len, out2, det_len, X_lp); - for (k = lf - 2; k < 2 * det_len; ++k) { - X[index3 + k - lf + 2] = X_lp[k]; - } - index += det_len; - index3 += index2; - index4 += 2 * indexp; - ptemp[l] = 2; - } - else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) { - for (k = 0; k < det_len; ++k) { - out[k] = X[index4 + k]; - out2[k] = X[index4 + indexp + k]; - } - idwpt_sym(wt, out, det_len, out2, det_len, X_lp); - for (k = lf - 2; k < 2 * det_len; ++k) { - X[index3 + k - lf + 2] = X_lp[k]; - } - index4 += 2 * indexp; - index3 += index2; - ptemp[l] = 2; - } - else { - index3 += index2; - index4 += 2 * indexp; - } - - } - - //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); - /* - idwpt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp); - for (k = lf - 2; k < 2 * det_len; ++k) { - out[k - lf + 2] = X_lp[k]; - } - - iter += det_len; - det_len = wt->length[i + 2]; - */ - llb /= 2; - indexp = index2; - } - - //free(X_lp); - - } - else { - printf("Signal extension can be either per or sym"); - exit(-1); - } - - for (i = 0; i < wt->siglength; ++i) { - //printf("%g ", X[i]); - dwtop[i] = X[i]; - } - + if (wt->basisvector[0] == 1) { + for (i = 0; i < wt->siglength; ++i) { + dwtop[i] = wt->output[i]; + } + + } + else { + for (i = 0; i < J; ++i) { + llb *= 2; + n1 += llb; + } + + for (i = 0; i < xlen; ++i) { + X[i] = 0.0; + } + + for (i = 0; i < llb; ++i) { + prep[i] = (int)wt->basisvector[n1 - llb + i]; + ptemp[i] = 0; + } + + if (!strcmp(wt->ext, "per")) { + app_len = wt->length[0]; + det_len = wt->length[1]; + index = 0; + + + for (i = 0; i < J; ++i) { + p = ipow2(J - i - 1); + det_len = wt->length[i + 1]; + index2 *= 2; + index3 = 0; + index4 = 0; + //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); + n1 -= llb; + for (l = 0; l < llb; ++l) { + if (ptemp[l] != 2) { + prep[l] = (int)wt->basisvector[n1 + l]; + } + else { + prep[l] = ptemp[l]; + } + ptemp[l] = 0; + } + + + for (l = 0; l < p; ++l) { + if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) { + for (k = 0; k < det_len; ++k) { + out[k] = wt->output[index + k]; + out2[k] = wt->output[index + det_len + k]; + } + idwpt_per(wt, out, det_len, out2, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index += 2 * det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) { + index4 += indexp; + for (k = 0; k < det_len; ++k) { + out[k] = wt->output[index + k]; + out2[k] = X[index4 + k]; + } + idwpt_per(wt, out, det_len, out2, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = wt->output[index + k]; + } + idwpt_per(wt, out, det_len, out2, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = X[index4 + indexp + k]; + } + idwpt_per(wt, out, det_len, out2, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index4 += 2 * indexp; + index3 += index2; + ptemp[l] = 2; + } + else { + index3 += index2; + index4 += 2 * indexp; + } + + } + + + /* + idwt_per(wt, out, det_len, wt->output + iter, det_len, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + out[k - lf / 2 + 1] = X_lp[k]; + } + + iter += det_len; + det_len = wt->length[i + 2]; + */ + llb /= 2; + indexp = index2; + } + + //free(X_lp); + + } + else if (!strcmp(wt->ext, "sym")) { + app_len = wt->length[0]; + det_len = wt->length[1]; + + //X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); + index = 0; + + for (i = 0; i < J; ++i) { + p = ipow2(J - i - 1); + det_len = wt->length[i + 1]; + index2 *= 2; + index3 = 0; + index4 = 0; + //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); + n1 -= llb; + for (l = 0; l < llb; ++l) { + if (ptemp[l] != 2) { + prep[l] = (int)wt->basisvector[n1 + l]; + } + else { + prep[l] = ptemp[l]; + } + ptemp[l] = 0; + } + + + for (l = 0; l < p; ++l) { + if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) { + for (k = 0; k < det_len; ++k) { + out[k] = wt->output[index + k]; + out2[k] = wt->output[index + det_len + k]; + } + idwpt_sym(wt, out, det_len, out2, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index += 2 * det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) { + index4 += indexp; + for (k = 0; k < det_len; ++k) { + out[k] = wt->output[index + k]; + out2[k] = X[index4 + k]; + } + idwpt_sym(wt, out, det_len, out2, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = wt->output[index + k]; + } + idwpt_sym(wt, out, det_len, out2, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = X[index4 + indexp + k]; + } + idwpt_sym(wt, out, det_len, out2, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index4 += 2 * indexp; + index3 += index2; + ptemp[l] = 2; + } + else { + index3 += index2; + index4 += 2 * indexp; + } + + } + + //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); + /* + idwpt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + out[k - lf + 2] = X_lp[k]; + } + + iter += det_len; + det_len = wt->length[i + 2]; + */ + llb /= 2; + indexp = index2; + } + + //free(X_lp); + + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } + + for (i = 0; i < wt->siglength; ++i) { + //printf("%g ", X[i]); + dwtop[i] = X[i]; + } + } @@ -2029,7 +2034,7 @@ void idwpt(wpt_object wt, double *dwtop) { } -static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD) { int l, l2, isodd, i, t, len_avg,j; len_avg = M * wt->wave->lpd_len; @@ -2079,7 +2084,7 @@ static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_ } -static void swt_fft(wt_object wt, double *inp) { +static void swt_fft(wt_object wt, const double *inp) { int i, J, temp_len, iter, M, N, len_filt; int lenacc; double *low_pass, *high_pass,*sig,*cA,*cD; @@ -2164,16 +2169,15 @@ static void swt_fft(wt_object wt, double *inp) { free(cD); } -static void swt_direct(wt_object wt, double *inp) { - int i, J, temp_len, iter, M, N; - int lenacc, len_filt; +static void swt_direct(wt_object wt, const double *inp) { + int i, J, temp_len, iter, M; + int lenacc; double *cA, *cD; temp_len = wt->siglength; J = wt->J; wt->length[0] = wt->length[J] = temp_len; wt->outlength = wt->length[J + 1] = (J + 1) * temp_len; - len_filt = wt->wave->filtlength; M = 1; for (iter = 1; iter < J; ++iter) { M = 2 * M; @@ -2196,13 +2200,9 @@ static void swt_direct(wt_object wt, double *inp) { lenacc -= temp_len; if (iter > 0) { M = 2 * M; - N = M * len_filt; - } - else { - N = len_filt; } - swt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len); + swt_per(wt, M, wt->params, temp_len, cA, temp_len, cD); for (i = 0; i < temp_len; ++i) { @@ -2218,7 +2218,7 @@ static void swt_direct(wt_object wt, double *inp) { } -void swt(wt_object wt, double *inp) { +void swt(wt_object wt, const double *inp) { if (!strcmp(wt->method, "swt") && !strcmp(wt->cmethod, "direct") ) { swt_direct(wt,inp); } @@ -2376,7 +2376,7 @@ void iswt(wt_object wt, double *swtop) { free(det2); } -static void modwt_per(wt_object wt, int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { +static void modwt_per(wt_object wt, int M, double *inp, double *cA, int len_cA, double *cD) { int l, i, t, len_avg; double s; double *filt; @@ -2410,16 +2410,15 @@ static void modwt_per(wt_object wt, int M, double *inp, int N, double *cA, int l free(filt); } -void modwt(wt_object wt, double *inp) { - int i, J, temp_len, iter, M, N; - int lenacc, len_filt; +void modwt(wt_object wt, const double *inp) { + int i, J, temp_len, iter, M; + int lenacc; double *cA, *cD; temp_len = wt->siglength; J = wt->J; wt->length[0] = wt->length[J] = temp_len; wt->outlength = wt->length[J + 1] = (J + 1) * temp_len; - len_filt = wt->wave->filtlength; M = 1; for (iter = 1; iter < J; ++iter) { M = 2 * M; @@ -2442,13 +2441,9 @@ void modwt(wt_object wt, double *inp) { lenacc -= temp_len; if (iter > 0) { M = 2 * M; - N = M * len_filt; - } - else { - N = len_filt; } - modwt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len); + modwt_per(wt, M, wt->params, cA, temp_len, cD); for (i = 0; i < temp_len; ++i) { @@ -2463,7 +2458,7 @@ void modwt(wt_object wt, double *inp) { } -static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, int len_cD, double *X) { +static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, double *X) { int len_avg, i, l, t; double s; double *filt; @@ -2497,14 +2492,12 @@ static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, i } void imodwt(wt_object wt, double *dwtop) { - int N, lf, iter, i, J, j, U; + int N, iter, i, J, j; int lenacc,M; double *X; N = wt->siglength; J = wt->J; - U = 2; - lf = wt->wave->lpr_len; lenacc = N; M = (int)pow(2.0, (double)J - 1.0); //M = 1; @@ -2518,7 +2511,7 @@ void imodwt(wt_object wt, double *dwtop) { if (iter > 0) { M = M / 2; } - imodwt_per(wt, M, dwtop, N, wt->params + lenacc, N, X); + imodwt_per(wt, M, dwtop, N, wt->params + lenacc, X); /* for (j = lf - 1; j < N; ++j) { dwtop[j - lf + 1] = X[j]; @@ -2536,7 +2529,7 @@ void imodwt(wt_object wt, double *dwtop) { free(X); } -void setDWTExtension(wt_object wt, char *extension) { +void setDWTExtension(wt_object wt, const char *extension) { if (!strcmp(extension, "sym")) { strcpy(wt->ext, "sym"); } @@ -2549,7 +2542,7 @@ void setDWTExtension(wt_object wt, char *extension) { } } -void setWTREEExtension(wtree_object wt, char *extension) { +void setWTREEExtension(wtree_object wt, const char *extension) { if (!strcmp(extension, "sym")) { strcpy(wt->ext, "sym"); } @@ -2562,7 +2555,7 @@ void setWTREEExtension(wtree_object wt, char *extension) { } } -void setDWPTExtension(wpt_object wt, char *extension) { +void setDWPTExtension(wpt_object wt, const char *extension) { if (!strcmp(extension, "sym")) { strcpy(wt->ext, "sym"); } @@ -2575,7 +2568,7 @@ void setDWPTExtension(wpt_object wt, char *extension) { } } -void setDWPTEntropy(wpt_object wt, char *entropy, double eparam) { +void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam) { if (!strcmp(entropy, "shannon")) { strcpy(wt->entropy, "shannon"); } @@ -2596,7 +2589,7 @@ void setDWPTEntropy(wpt_object wt, char *entropy, double eparam) { } } -void setWTConv(wt_object wt, char *cmethod) { +void setWTConv(wt_object wt, const char *cmethod) { if (!strcmp(cmethod, "fft") || !strcmp(cmethod, "FFT")) { strcpy(wt->cmethod, "fft"); } diff --git a/src/wavelib.h b/src/wavelib.h deleted file mode 100644 index 86387f4..0000000 --- a/src/wavelib.h +++ /dev/null @@ -1,230 +0,0 @@ -/* -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_ */ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a184644..1e5f51c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,32 +1,42 @@ add_executable(cwttest cwttest.c) -target_link_libraries(cwttest wavelib m) +target_link_libraries(cwttest wavelib) add_executable(dwttest dwttest.c) -target_link_libraries(dwttest wavelib m) +target_link_libraries(dwttest wavelib) add_executable(swttest swttest.c) -target_link_libraries(swttest wavelib m) +target_link_libraries(swttest wavelib) add_executable(modwttest modwttest.c) -target_link_libraries(modwttest wavelib m) +target_link_libraries(modwttest wavelib) add_executable(dwpttest dwpttest.c) -target_link_libraries(dwpttest wavelib m) +target_link_libraries(dwpttest wavelib) add_executable(wtreetest wtreetest.c) -target_link_libraries(wtreetest wavelib m) +target_link_libraries(wtreetest wavelib) add_executable(denoisetest denoisetest.c) -target_link_libraries(denoisetest wauxlib wavelib m) +target_link_libraries(denoisetest wauxlib wavelib) + +if(UNIX) + target_link_libraries(cwttest m) + target_link_libraries(dwttest m) + target_link_libraries(swttest m) + target_link_libraries(modwttest m) + target_link_libraries(dwpttest m) + target_link_libraries(wtreetest m) + target_link_libraries(denoisetest m) +endif() set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/test" - ) \ No newline at end of file + ) diff --git a/unitTests/wavelibBoostTests/tst_dwt.cpp b/unitTests/wavelibBoostTests/tst_dwt.cpp index 6fe71e6..ad51519 100644 --- a/unitTests/wavelibBoostTests/tst_dwt.cpp +++ b/unitTests/wavelibBoostTests/tst_dwt.cpp @@ -4,10 +4,12 @@ #include #include -#include -#include #include +#include +#include +#include + #include "wavelib.h" #include @@ -117,7 +119,6 @@ void ReconstructionTest() double *inp,*out; int N, i,J; double epsilon = 1e-15; - char *type = (char*) "dwt"; N = 79926; @@ -346,10 +347,10 @@ void DWPTReconstructionTest() } void CWTReconstructionTest() { - int i, j,N, J,subscale,a0,iter; + int i, N, J,subscale,a0; double *inp,*oup; double dt, dj,s0, pi,t; - double val, epsilon; + double epsilon; int it1,it2; cwt_object wt;