MODWT branch added. Test for memory leaks

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);
 	}
 }

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);

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

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");

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"
         )

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");
 
 

test/modwtdenoisetest.c

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

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<std::string > 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");