From e16859af59b0e311a0e5d70b39a756810f9d1d77 Mon Sep 17 00:00:00 2001 From: Rafat Hussain Date: Wed, 20 Mar 2019 18:27:57 +0530 Subject: [PATCH] MODWT branch added. Test for memory leaks --- auxiliary/denoise.c | 132 ++++++++++- header/wauxlib.h | 4 + header/wavelib.h | 1 + src/wavelib.c | 298 +++++++++++++++++++++++- test/CMakeLists.txt | 7 +- test/denoisetest.c | 21 +- test/modwtdenoisetest.c | 123 ++++++++++ unitTests/wavelibBoostTests/tst_dwt.cpp | 96 ++++++++ 8 files changed, 665 insertions(+), 17 deletions(-) mode change 100644 => 100755 header/wavelib.h create mode 100755 test/modwtdenoisetest.c diff --git a/auxiliary/denoise.c b/auxiliary/denoise.c index 5354252..cb87011 100644 --- a/auxiliary/denoise.c +++ b/auxiliary/denoise.c @@ -19,9 +19,10 @@ denoise_object denoise_init(int length, int J,const char* wname) { //Set Default Values strcpy(obj->dmethod,"sureshrink"); strcpy(obj->ext,"sym"); - strcpy(obj->level,"first"); + strcpy(obj->level,"all"); strcpy(obj->thresh,"soft"); strcpy(obj->wmethod,"dwt"); + strcpy(obj->cmethod,"direct"); return obj; } @@ -280,11 +281,130 @@ void sureshrink(double *signal,int N,int J,const char *wname,const char *method, wt_free(wt); } +void modwtshrink(double *signal, int N, int J, const char *wname, const char *cmethod, const char *ext, const char *thresh, double *denoised) { + int filt_len, iter, i, dlen, sgn, MaxIter, it; + double sigma, td, tmp, M, llen; + wave_object wave; + wt_object wt; + double *dout, *lnoise; + + wave = wave_init(wname); + + filt_len = wave->filtlength; + + MaxIter = (int)(log((double)N / ((double)filt_len - 1.0)) / log(2.0)); + + if (J > MaxIter) { + printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter); + exit(-1); + } + + wt = wt_init(wave, "modwt", N, J); + + if (!strcmp(ext, "sym") && !strcmp(cmethod,"fft")) { + setWTConv(wt, "fft"); + setDWTExtension(wt, "sym"); + } + else if (!strcmp(ext, "sym") && !strcmp(cmethod, "direct")) { + printf("Symmetric Extension is not available for direct method"); + exit(-1); + } + else if (!strcmp(ext, "per") && !strcmp(cmethod, "direct")) { + setWTConv(wt, "direct"); + setDWTExtension(wt, "per"); + } + else if (!strcmp(ext, "per") && !strcmp(cmethod, "fft")) { + setWTConv(wt, "fft"); + setDWTExtension(wt, "per"); + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } + + modwt(wt, signal); + + lnoise = (double*)malloc(sizeof(double)* J); + + //Set sigma + + iter = wt->length[0]; + dlen = wt->length[J]; + dout = (double*)malloc(sizeof(double)* dlen); + + for (it = 0; it < J; ++it) { + dlen = wt->length[it + 1]; + for (i = 0; i < dlen; ++i) { + dout[i] = fabs(wt->output[iter + i]); + } + + sigma = sqrt(2.0) * median(dout, dlen) / 0.6745; + lnoise[it] = sigma; + iter += dlen; + } + + M = pow(2.0,J); + llen = log((double)wt->modwtsiglength); + // Thresholding + + iter = wt->length[0]; + for (it = 0; it < J; ++it) { + sigma = lnoise[it]; + dlen = wt->length[it + 1]; + td = sqrt(2.0 * llen / M) * sigma; + + if (!strcmp(thresh, "hard")) { + for (i = 0; i < dlen; ++i) { + if (fabs(wt->output[iter + i]) < td) { + wt->output[iter + i] = 0; + } + } + } + else if (!strcmp(thresh, "soft")) { + for (i = 0; i < dlen; ++i) { + if (fabs(wt->output[iter + i]) < td) { + wt->output[iter + i] = 0; + } + else { + sgn = wt->output[iter + i] >= 0 ? 1 : -1; + tmp = sgn * (fabs(wt->output[iter + i]) - td); + wt->output[iter + i] = tmp; + } + } + } + + iter += wt->length[it + 1]; + M /= 2.0; + } + + imodwt(wt, denoised); + + free(dout); + free(lnoise); + wave_free(wave); + wt_free(wt); +} + + void denoise(denoise_object obj, double *signal,double *denoised) { if(!strcmp(obj->dmethod,"sureshrink")) { + if (!strcmp(obj->wmethod, "modwt")) { + printf("sureshrink method only works with swt and dwt. Please use setDenoiseWTMethod to set the correct method\n"); + exit(-1); + } sureshrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised); } else if(!strcmp(obj->dmethod,"visushrink")) { + if (!strcmp(obj->wmethod, "modwt")) { + printf("visushrink method only works with swt and dwt. Please use setDenoiseWTMethod to set the correct method\n"); + exit(-1); + } visushrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);; + } else if(!strcmp(obj->dmethod,"modwtshrink")) { + if (strcmp(obj->wmethod, "modwt")) { + printf("modwtshrink method only works with modwt. Please use setDenoiseWTMethod to set the correct method\n"); + exit(-1); + } + modwtshrink(signal,obj->N,obj->J,obj->wname,obj->cmethod,obj->ext,obj->thresh,denoised);; } else { printf("Acceptable Denoising methods are - sureshrink and visushrink\n"); exit(-1); @@ -297,9 +417,12 @@ void setDenoiseMethod(denoise_object obj, const char *dmethod) { } else if (!strcmp(dmethod, "visushrink")) { strcpy(obj->dmethod, "visushrink"); + } + else if (!strcmp(dmethod, "modwtshrink")) { + strcpy(obj->dmethod, "modwtshrink"); } else { - printf("Acceptable Denoising methods are - sureshrink and visushrink\n"); + printf("Acceptable Denoising methods are - sureshrink, visushrink and modwtshrink\n"); exit(-1); } } @@ -310,9 +433,12 @@ void setDenoiseWTMethod(denoise_object obj, const char *wmethod) { } else if (!strcmp(wmethod, "swt")) { strcpy(obj->wmethod, "swt"); + } + else if (!strcmp(wmethod, "modwt")) { + strcpy(obj->wmethod, "modwt"); } else { - printf("Wavelet decomposition method can be either dwt or swt"); + printf("Wavelet decomposition method can be one of dwt, modwt or swt.\n"); exit(-1); } } diff --git a/header/wauxlib.h b/header/wauxlib.h index 53e4647..65c7e1e 100644 --- a/header/wauxlib.h +++ b/header/wauxlib.h @@ -19,6 +19,8 @@ struct denoise_set{ int J; // Levels of Wavelet decomposition char wname[10]; //Wavelet name char wmethod[10]; //Wavelet decomposition method - dwt or swt + char cmethod[10]; //Cnvolution Method - direct or fft . Available only for modwt. + // SWT and DWT only use direct method. char ext[10]; // Signal Extension - sym or per char thresh[10]; // thresholding - soft or hard char level[10]; // Noise Estimation level - first or all @@ -30,6 +32,8 @@ void visushrink(double *signal,int N,int J,const char *wname,const char *method, void sureshrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised); +void modwtshrink(double *signal, int N, int J, const char *wname, const char *cmethod, const char *ext, const char *thresh, double *denoised); + void denoise(denoise_object obj, double *signal,double *denoised); void setDenoiseMethod(denoise_object obj, const char *dmethod); diff --git a/header/wavelib.h b/header/wavelib.h old mode 100644 new mode 100755 index e36096b..9fd18d5 --- a/header/wavelib.h +++ b/header/wavelib.h @@ -90,6 +90,7 @@ struct wt_set{ conv_object cobj; char method[10]; int siglength;// Length of the original signal. + int modwtsiglength; // Modified signal length for MODWT int outlength;// Length of the output DWT vector int lenlength;// Length of the Output Dimension Vector "length" int J; // Number of decomposition Levels diff --git a/src/wavelib.c b/src/wavelib.c index aaf028c..5615b2f 100644 --- a/src/wavelib.c +++ b/src/wavelib.c @@ -88,13 +88,14 @@ wt_object wt_init(wave_object wave,const char* method, int siglength,int J) { } } - obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1))); + obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * 2 * (J + 1))); obj->outlength = siglength * (J + 1); // Default strcpy(obj->ext, "per"); // Default } obj->wave = wave; obj->siglength = siglength; + obj->modwtsiglength = siglength; obj->J = J; obj->MaxIter = MaxIter; strcpy(obj->method, method); @@ -117,11 +118,16 @@ wt_object wt_init(wave_object wave,const char* method, int siglength,int J) { obj->params[i] = 0.0; } } - else if (!strcmp(method, "swt") || !strcmp(method, "SWT") || !strcmp(method, "modwt") || !strcmp(method, "MODWT")) { + else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) { for (i = 0; i < siglength * (J + 1); ++i) { obj->params[i] = 0.0; } } + else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) { + for (i = 0; i < siglength * 2 * (J + 1); ++i) { + obj->params[i] = 0.0; + } + } //wave_summary(obj->wave); return obj; @@ -2410,11 +2416,17 @@ static void modwt_per(wt_object wt, int M, double *inp, double *cA, int len_cA, free(filt); } -void modwt(wt_object wt, const double *inp) { +static void modwt_direct(wt_object wt, const double *inp) { int i, J, temp_len, iter, M; int lenacc; double *cA, *cD; + if (strcmp(wt->ext, "per")) { + printf("MODWT direct method only uses periodic extension per. \n"); + printf(" Use MODWT fft method for symmetric extension sym \n"); + exit(-1); + } + temp_len = wt->siglength; J = wt->J; wt->length[0] = wt->length[J] = temp_len; @@ -2458,6 +2470,270 @@ void modwt(wt_object wt, const double *inp) { } +static void modwt_fft(wt_object wt, const double *inp) { + int i, J, temp_len, iter, M,N, len_avg; + int lenacc; + double s,tmp1,tmp2; + fft_data *cA, *cD, *low_pass,*high_pass,*sig; + int *index; + fft_object fft_fd = NULL; + fft_object fft_bd = NULL; + + temp_len = wt->siglength; + len_avg = wt->wave->lpd_len; + if (!strcmp(wt->ext, "sym")) { + N = 2 * temp_len; + } else if (!strcmp(wt->ext, "per")) { + N = temp_len; + } + J = wt->J; + wt->modwtsiglength = N; + wt->length[0] = wt->length[J] = N; + wt->outlength = wt->length[J + 1] = (J + 1) * N; + + s = sqrt(2.0); + for (iter = 1; iter < J; ++iter) { + wt->length[iter] = N; + } + + fft_fd = fft_init(N, 1); + fft_bd = fft_init(N, -1); + + sig = (fft_data*)malloc(sizeof(fft_data)* N); + cA = (fft_data*)malloc(sizeof(fft_data)* N); + cD = (fft_data*)malloc(sizeof(fft_data)* N); + low_pass = (fft_data*)malloc(sizeof(fft_data)* N); + high_pass = (fft_data*)malloc(sizeof(fft_data)* N); + index = (int*)malloc(sizeof(int)*N); + + + // N-point FFT of low pass and high pass filters + + // Low Pass Filter + + for(i = 0; i < len_avg;++i) { + sig[i].re = (fft_type) wt->wave->lpd[i] / s; + sig[i].im = 0.0; + } + for(i = len_avg; i < N;++i) { + sig[i].re = 0.0; + sig[i].im = 0.0; + } + + + fft_exec(fft_fd, sig, low_pass); + + // High Pass Filter + + for (i = 0; i < len_avg; ++i) { + sig[i].re = (fft_type)wt->wave->hpd[i] / s; + sig[i].im = 0.0; + } + for (i = len_avg; i < N; ++i) { + sig[i].re = 0.0; + sig[i].im = 0.0; + } + + fft_exec(fft_fd, sig, high_pass); + + // symmetric extension + for (i = 0; i < temp_len; ++i) { + sig[i].re = (fft_type)inp[i]; + sig[i].im = 0.0; + } + for (i = temp_len; i < N; ++i) { + sig[i].re = (fft_type) inp[N-i-1]; + sig[i].im = 0.0; + } + + // FFT of data + + fft_exec(fft_fd, sig, cA); + + lenacc = wt->outlength; + + M = 1; + + for (iter = 0; iter < J; ++iter) { + lenacc -= N; + + for (i = 0; i < N; ++i) { + index[i] = (M *i) % N; + } + + for (i = 0; i < N; ++i) { + tmp1 = cA[i].re; + tmp2 = cA[i].im; + cA[i].re = low_pass[index[i]].re*tmp1 - low_pass[index[i]].im*tmp2; + cA[i].im = low_pass[index[i]].re*tmp2 + low_pass[index[i]].im*tmp1; + + cD[i].re = high_pass[index[i]].re*tmp1 - high_pass[index[i]].im*tmp2; + cD[i].im = high_pass[index[i]].re*tmp2 + high_pass[index[i]].im*tmp1; + } + + fft_exec(fft_bd, cD, sig); + + for (i = 0; i < N; ++i) { + wt->params[lenacc + i] = sig[i].re/N; + } + + M *= 2; + } + + fft_exec(fft_bd, cA, sig); + + for (i = 0; i < N; ++i) { + wt->params[i] = sig[i].re/N; + } + + free(sig); + free(cA); + free(cD); + free(low_pass); + free(high_pass); + free_fft(fft_fd); + free_fft(fft_bd); +} + +void modwt(wt_object wt, const double *inp) { + if (!strcmp(wt->cmethod, "direct")) { + modwt_direct(wt, inp); + } + else if (!strcmp(wt->cmethod, "fft")) { + modwt_fft(wt, inp); + } + else { + printf("Error- Available Choices for this method are - direct and fft \n"); + exit(-1); + } +} + +static void conj_complex(fft_data *x, int N) { + int i; + + for (i = 0; i < N; ++i) { + x[i].im *= (-1.0); + } +} + +void imodwt_fft(wt_object wt, double *oup) { + int i, J, temp_len, iter, M, N, len_avg; + int lenacc; + double s, tmp1, tmp2; + fft_data *cA, *cD, *low_pass, *high_pass, *sig; + int *index; + fft_object fft_fd = NULL; + fft_object fft_bd = NULL; + + N = wt->modwtsiglength; + len_avg = wt->wave->lpd_len; + if (!strcmp(wt->ext, "sym")) { + temp_len = N/2; + } + else if (!strcmp(wt->ext, "per")) { + temp_len = N; + } + J = wt->J; + + s = sqrt(2.0); + fft_fd = fft_init(N, 1); + fft_bd = fft_init(N, -1); + + sig = (fft_data*)malloc(sizeof(fft_data)* N); + cA = (fft_data*)malloc(sizeof(fft_data)* N); + cD = (fft_data*)malloc(sizeof(fft_data)* N); + low_pass = (fft_data*)malloc(sizeof(fft_data)* N); + high_pass = (fft_data*)malloc(sizeof(fft_data)* N); + index = (int*)malloc(sizeof(int)*N); + + + // N-point FFT of low pass and high pass filters + + // Low Pass Filter + + for (i = 0; i < len_avg; ++i) { + sig[i].re = (fft_type)wt->wave->lpd[i] / s; + sig[i].im = 0.0; + } + for (i = len_avg; i < N; ++i) { + sig[i].re = 0.0; + sig[i].im = 0.0; + } + + + fft_exec(fft_fd, sig, low_pass); + + // High Pass Filter + + for (i = 0; i < len_avg; ++i) { + sig[i].re = (fft_type)wt->wave->hpd[i] / s; + sig[i].im = 0.0; + } + for (i = len_avg; i < N; ++i) { + sig[i].re = 0.0; + sig[i].im = 0.0; + } + + fft_exec(fft_fd, sig, high_pass); + + + // Complex conjugate of the two filters + + conj_complex(low_pass, N); + conj_complex(high_pass, N); + + M = (int)pow(2.0, (double)J - 1.0); + lenacc = N; + + // + for (i = 0; i < N; ++i) { + sig[i].re = (fft_type)wt->output[i]; + sig[i].im = 0.0; + } + + for (iter = 0; iter < J; ++iter) { + fft_exec(fft_fd, sig, cA); + for (i = 0; i < N; ++i) { + sig[i].re = wt->output[lenacc+i]; + sig[i].im = 0.0; + } + fft_exec(fft_fd, sig, cD); + + for (i = 0; i < N; ++i) { + index[i] = (M *i) % N; + } + + for (i = 0; i < N; ++i) { + tmp1 = cA[i].re; + tmp2 = cA[i].im; + cA[i].re = low_pass[index[i]].re*tmp1 - low_pass[index[i]].im*tmp2 + high_pass[index[i]].re*cD[i].re - high_pass[index[i]].im*cD[i].im; + cA[i].im = low_pass[index[i]].re*tmp2 + low_pass[index[i]].im*tmp1 + high_pass[index[i]].re*cD[i].im + high_pass[index[i]].im*cD[i].re; + } + + fft_exec(fft_bd, cA, sig); + + for (i = 0; i < N; ++i) { + sig[i].re /= N; + sig[i].im /= N; + } + M /= 2; + lenacc += N; + } + + for (i = 0; i < wt->siglength; ++i) { + oup[i] = sig[i].re; + } + + free(sig); + free(cA); + free(cD); + free(low_pass); + free(high_pass); + free_fft(fft_fd); + free_fft(fft_bd); +} + + 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; @@ -2491,7 +2767,7 @@ static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, d free(filt); } -void imodwt(wt_object wt, double *dwtop) { +static void imodwt_direct(wt_object wt, double *dwtop) { int N, iter, i, J, j; int lenacc,M; double *X; @@ -2529,6 +2805,20 @@ void imodwt(wt_object wt, double *dwtop) { free(X); } +void imodwt(wt_object wt, double *oup) { + if (!strcmp(wt->cmethod, "direct")) { + imodwt_direct(wt, oup); + } + else if (!strcmp(wt->cmethod, "fft")) { + imodwt_fft(wt, oup); + } + else { + printf("Error- Available Choices for this method are - direct and fft \n"); + exit(-1); + } +} + + void setDWTExtension(wt_object wt, const char *extension) { if (!strcmp(extension, "sym")) { strcpy(wt->ext, "sym"); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1e5f51c..ac9a1f8 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -26,6 +26,10 @@ add_executable(denoisetest denoisetest.c) target_link_libraries(denoisetest wauxlib wavelib) +add_executable(modwtdenoisetest modwtdenoisetest.c) + +target_link_libraries(modwtdenoisetest wauxlib wavelib) + if(UNIX) target_link_libraries(cwttest m) target_link_libraries(dwttest m) @@ -34,9 +38,10 @@ if(UNIX) target_link_libraries(dwpttest m) target_link_libraries(wtreetest m) target_link_libraries(denoisetest m) + target_link_libraries(modwtdenoisetest m) endif() -set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest +set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest modwtdenoisetest PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/test" ) diff --git a/test/denoisetest.c b/test/denoisetest.c index 1f17d28..8b7e6ec 100644 --- a/test/denoisetest.c +++ b/test/denoisetest.c @@ -56,10 +56,12 @@ int main() { double temp[2400]; char *wname = "db5"; - char *method = "dwt"; - char *ext = "sym"; - char *thresh = "soft"; - char *level = "all"; + char *method = "dwt";// Available - dwt, swt and modwt. modwt works only with modwtshrink. The other two methods work with + // visushrink and sureshrink + char *ext = "sym";// sym and per work with dwt. swt and modwt only use per extension when called through denoise. + // You can use sy extension if you directly call modwtshrink with cmethod set to fft. See modwtdenoisetest.c file + char *thresh = "soft";// soft or hard + char *level = "all"; // noise estimation at "first" or "all" levels. modwt only has the option of "all" ifp = fopen("pieceregular1024.txt", "r"); i = 0; @@ -104,19 +106,20 @@ int main() { 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 + setDenoiseMethod(obj,"visushrink");// sureshrink is also the default. The other option with dwt and swt is visushrink. + // modwt works only with modwtshrink method + setDenoiseWTMethod(obj,method);// Default is dwt. the other options are swt and modwt setDenoiseWTExtension(obj,ext);// Default is sym. the other option is per setDenoiseParameters(obj,thresh,level);// Default for thresh is soft. Other option is hard - // Default for level is first. The other option is all + // Default for level is all. The other option is first denoise(obj,inp,oup); // Alternative to denoise_object - // Just use visushrink and sureshrink functions + // Just use visushrink, modwtshrink and sureshrink functions //visushrink(inp,N,J,wname,method,ext,thresh,level,oup); //sureshrink(inp,N,J,wname,method,ext,thresh,level,oup); - + // modwtshrink(sig,N,J,wname,cmethod,ext,thresh,oup); See modwtdenoisetest.c //ofp = fopen("denoiseds.txt", "w"); diff --git a/test/modwtdenoisetest.c b/test/modwtdenoisetest.c new file mode 100755 index 0000000..62f9728 --- /dev/null +++ b/test/modwtdenoisetest.c @@ -0,0 +1,123 @@ +#include +#include +#include +#include + +#include "../header/wauxlib.h" + +static double rmse(int N,double *x,double *y) { + double rms; + int i; + + rms = 0.0; + + for(i = 0; i < N;++i) { + rms += (x[i] - y[i]) * (x[i] - y[i]); + } + + rms = sqrt(rms/(double)N); + + return rms; +} + +static double corrcoef(int N,double *x,double *y) { + double cc,xm,ym,tx,ty,num,den1,den2; + int i; + xm = ym = 0.0; + for(i = 0; i < N;++i) { + xm += x[i]; + ym += y[i]; + } + + xm = xm/N; + ym = ym / N; + num = den1 = den2 = 0.0; + + for(i = 0; i < N;++i) { + tx = x[i] - xm; + ty = y[i] - ym; + num += (tx*ty); + den1 += (tx*tx); + den2 += (ty*ty); + } + + cc = num / sqrt(den1*den2); + + return cc; +} + +int main() { + // gcc -Wall -I../header -L../Bin modwtdenoisetest.c -o modwtdenoise -lwauxlib -lwavelib -lm + /* + modwtshrink can also be called from the denoise object. See denoisetest.c for more information + */ + double *sig,*inp,*oup; + int i,N,J; + FILE *ifp; + + double temp[2400]; + + char *wname = "db5"; + char *ext = "per";// The other option sym is only available with "fft" cmethod + char *thresh = "soft"; + char *cmethod = "direct";// The other option is "fft" + + ifp = fopen("pieceregular1024.txt", "r"); + i = 0; + if (!ifp) { + printf("Cannot Open File"); + exit(100); + } + + while (!feof(ifp)) { + fscanf(ifp, "%lf \n", &temp[i]); + i++; + } + + fclose(ifp); + + N = i; + J = 4; + + sig = (double*)malloc(sizeof(double)* N); + inp = (double*)malloc(sizeof(double)* N); + oup = (double*)malloc(sizeof(double)* N); + + for(i = 0; i < N;++i) { + sig[i] = temp[i]; + } + + ifp = fopen("PieceRegular10.txt", "r"); + i = 0; + if (!ifp) { + printf("Cannot Open File"); + exit(100); + } + + while (!feof(ifp)) { + fscanf(ifp, "%lf \n", &temp[i]); + i++; + } + + fclose(ifp); + + for(i = 0; i < N;++i) { + inp[i] = temp[i]; + } + + modwtshrink(sig,N,J,wname,cmethod,ext,thresh,oup); + + + printf("Signal - Noisy Signal Stats \n"); + printf("RMSE %g\n",rmse(N,sig,inp)); + printf("Corr Coeff %g\n",corrcoef(N,sig,inp)); + + printf("Signal - DeNoised Signal Stats \n"); + printf("RMSE %g\n",rmse(N,sig,oup)); + printf("Corr Coeff %g\n",corrcoef(N,sig,oup)); + + free(sig); + free(inp); + free(oup); + return 0; +} diff --git a/unitTests/wavelibBoostTests/tst_dwt.cpp b/unitTests/wavelibBoostTests/tst_dwt.cpp index ad51519..7d028a5 100644 --- a/unitTests/wavelibBoostTests/tst_dwt.cpp +++ b/unitTests/wavelibBoostTests/tst_dwt.cpp @@ -229,6 +229,99 @@ void ReconstructionTest() free(inp); } +void MODWTReconstructionTest() +{ + + wave_object obj; + wt_object wt; + double *inp,*out; + int N, i,J; + double epsilon = 1e-15; + double err; + + 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 waveletNames; + + for (unsigned int j = 0; j < 15; j++) + { + waveletNames.push_back(std::string("db") + patch::to_string(j + 1)); + } + for (unsigned int j = 0; j < 5; j++) + { + waveletNames.push_back(std::string("coif") + patch::to_string(j + 1)); + } + for (unsigned int j = 1; j < 10; j++) + { + waveletNames.push_back(std::string("sym") + patch::to_string(j + 1)); + } + + + 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*) "modwt", N, J);// Initialize the wavelet transform object + + if (direct_fft == 0) + setWTConv(wt, (char*) "direct"); + else + setWTConv(wt, (char*) "fft"); + + if (sym_per == 0) + setDWTExtension(wt, (char*) "per");// Options are "per" and "sym". Symmetric is the default option + else if (sym_per == 1 && direct_fft == 1) + setDWTExtension(wt, (char*) "sym"); + else break; + + modwt(wt, inp);// Perform DWT + + imodwt(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)); + err = RMS_Error(out, inp, wt->siglength); + //printf("%d %d %g \n",direct_fft,sym_per,err); + if (err > 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() { @@ -649,6 +742,9 @@ int main() { printf("Running DWT ReconstructionTests ... "); ReconstructionTest(); printf("DONE \n"); + printf("Running MODWT ReconstructionTests ... "); + MODWTReconstructionTest(); + printf("DONE \n"); printf("Running DWPT ReconstructionTests ... "); DWPTReconstructionTest(); printf("DONE \n");