Unit Tests updated : wavecoeffs memory Leak plugged

Merge MODWT with the master branch

.gitignore

@@ -1,6 +1,7 @@
 #Folders Ignore
 Bin/
 Testing/
+.vscode/
 
 #cmake ignore
 CMakeCache.txt

README.md

@@ -1,3 +1,5 @@
+[![Build Status](https://travis-ci.org/rafat/wavelib.svg?branch=master)](https://travis-ci.org/rafat/wavelib)
+
 wavelib
 =======
 

header/wavelib.h

@@ -191,11 +191,33 @@ struct cwt_set{
 	double params[0];
 };
 
+typedef struct wt2_set* wt2_object;
+
+wt2_object wt2_init(wave_object wave, const char* method, int rows, int cols, int J);
+
+struct wt2_set{
+	wave_object wave;
+	char method[10];
+	int rows;// Matrix Number of rows
+	int cols; // Matrix Number of columns
+	int outlength;// Length of the output DWT vector
+	int J; // Number of decomposition Levels
+	int MaxIter;// Maximum Iterations J <= MaxIter
+	char ext[10];// Type of Extension used - "per" or "sym"
+	int coeffaccesslength;
+
+	int N; //
+	int *dimensions;
+	int *coeffaccess;
+	int params[0];
+};
 
 void dwt(wt_object wt, const double *inp);
 
 void idwt(wt_object wt, double *dwtop);
 
+double *getDWTmra(wt_object wt, double *wavecoeffs);
+
 void wtree(wtree_object wt, const double *inp);
 
 void dwpt(wpt_object wt, const double *inp);
@@ -206,16 +228,22 @@ void swt(wt_object wt, const double *inp);
 
 void iswt(wt_object wt, double *swtop);
 
+double *getSWTmra(wt_object wt, double *wavecoeffs);
+
 void modwt(wt_object wt, const double *inp);
 
 void imodwt(wt_object wt, double *dwtop);
 
+double* getMODWTmra(wt_object wt, double *wavecoeffs);
+
 void setDWTExtension(wt_object wt, const char *extension);
 
 void setWTREEExtension(wtree_object wt, const char *extension);
 
 void setDWPTExtension(wpt_object wt, const char *extension);
 
+void setDWT2Extension(wt2_object wt, const char *extension);
+
 void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam);
 
 void setWTConv(wt_object wt, const char *cmethod);
@@ -240,6 +268,22 @@ void icwt(cwt_object wt, double *cwtop);
 
 int getCWTScaleLength(int N);
 
+double* dwt2(wt2_object wt, double *inp);
+
+void idwt2(wt2_object wt,double *wavecoeff, double *oup);
+
+double* swt2(wt2_object wt, double *inp);
+
+void iswt2(wt2_object wt, double *wavecoeffs, double *oup);
+
+double* modwt2(wt2_object wt, double *inp);
+
+void imodwt2(wt2_object wt, double *wavecoeff, double *oup);
+
+double* getWT2Coeffs(wt2_object wt,double* wcoeffs, int level,char *type, int *rows, int *cols);
+
+void dispWT2Coeffs(double *A, int row, int col);
+
 void wave_summary(wave_object obj);
 
 void wt_summary(wt_object wt);
@@ -248,7 +292,9 @@ void wtree_summary(wtree_object wt);
 
 void wpt_summary(wpt_object wt);
 
-void cwt_summary(cwt_object wt);;
+void cwt_summary(cwt_object wt);
+
+void wt2_summary(wt2_object wt);
 
 void wave_free(wave_object object);
 
@@ -260,6 +306,8 @@ void wpt_free(wpt_object object);
 
 void cwt_free(cwt_object object);
 
+void wt2_free(wt2_object wt);
+
 
 #ifdef __cplusplus
 }

src/wtmath.c

@@ -1,8 +1,355 @@
 /*
-Copyright (c) 2014, Rafat Hussain
+Copyright (c) 2018, Rafat Hussain
 */
 #include "wtmath.h"
 
+void dwt_per_stride(double *inp, int N, double *lpd,double*hpd,int lpd_len,double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int l, l2, isodd, i, t, len_avg,is,os;
+
+	len_avg = lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + l2;
+		os = i *ostride;
+		cA[os] = 0.0;
+		cD[os] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= l2 && (t - l) < N) {
+				is = (t - l) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < l2 && (t - l) >= 0) {
+				is = (t - l) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < 0 && isodd == 0) {
+				is = (t - l + N) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < 0 && isodd == 1) {
+				if ((t - l) != -1) {
+					is = (t - l + N + 1) * istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+				else {
+					is = (N - 1) * istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+			}
+			else if ((t - l) >= N && isodd == 0) {
+				is = (t - l - N) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) >= N && isodd == 1) {
+				is = (t - l - (N + 1)) * istride;
+				if (t - l != N) {
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+				else {
+					is = (N - 1) * istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+			}
+
+		}
+		
+	}
+
+}
+
+void dwt_sym_stride(double *inp, int N, double *lpd, double*hpd, int lpd_len, double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int i, l, t, len_avg;
+	int is, os;
+	len_avg = lpd_len;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + 1;
+		os = i *ostride;
+		cA[os] = 0.0;
+		cD[os] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= 0 && (t - l) < N) {
+				is = (t - l) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) < 0) {
+				is = (-t + l - 1) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - l) >= N) {
+				is = (2 * N - t + l - 1) * istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+
+		}
+	}
+
+
+}
+
+void modwt_per_stride(int M, double *inp, int N, double *filt, int lpd_len, double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int l, i, t, len_avg;
+	int is, os;
+	len_avg = lpd_len;
+
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i;
+		os = i *ostride;
+		is = t *istride;
+		cA[os] = filt[0] * inp[is];
+		cD[os] = filt[len_avg] * inp[is];
+		for (l = 1; l < len_avg; l++) {
+			t -= M;
+			while (t >= len_cA) {
+				t -= len_cA;
+			}
+			while (t < 0) {
+				t += len_cA;
+			}
+			os = i * ostride;
+			is = t * istride;
+			cA[os] += filt[l] * inp[is];
+			cD[os] += filt[len_avg + l] * inp[is];
+
+		}
+	}
+	
+}
+
+void swt_per_stride(int M, double *inp, int N, double *lpd, double*hpd, int lpd_len, double *cA, int len_cA, double *cD, int istride, int ostride) {
+	int l, l2, isodd, i, t, len_avg, j;
+	int is, os;
+	len_avg = M * lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i + l2;
+		os = i *ostride;
+		cA[os] = 0.0;
+		cD[os] = 0.0;
+		l = -1;
+		for (j = 0; j < len_avg; j += M) {
+			l++;
+			while (j >= len_cA) {
+				j -= len_cA;
+			}
+			if ((t - j) >= l2 && (t - j) < N) {
+				is = (t - j)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) < l2 && (t - j) >= 0) {
+				is = (t - j)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) < 0) {
+				is = (t - j + N)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) >= N && isodd == 0) {
+				is = (t - j - N)*istride;
+				cA[os] += lpd[l] * inp[is];
+				cD[os] += hpd[l] * inp[is];
+			}
+			else if ((t - j) >= N && isodd == 1) {
+				if (t - l != N) {
+					is = (t - j - (N + 1))*istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[is];
+				}
+				else {
+					is = (N - 1)*istride;
+					cA[os] += lpd[l] * inp[is];
+					cD[os] += hpd[l] * inp[N - 1];
+				}
+			}
+
+		}
+	}
+
+}
+
+void idwt_per_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, int lpr_len, double *X, int istride, int ostride) {
+	int len_avg, i, l, m, n, t, l2;
+	int is, ms, ns;
+
+	len_avg = lpr_len;
+	l2 = len_avg / 2;
+	m = -2;
+	n = -1;
+
+	for (i = 0; i < len_cA + l2 - 1; ++i) {
+		m += 2;
+		n += 2;
+		ms = m * ostride;
+		ns = n * ostride;
+		X[ms] = 0.0;
+		X[ns] = 0.0;
+		for (l = 0; l < l2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				is = (i - l) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+			else if ((i - l) >= len_cA && (i - l) < len_cA + len_avg - 1) {
+				is = (i - l - len_cA) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+			else if ((i - l) < 0 && (i - l) > -l2) {
+				is = (len_cA + i - l) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+		}
+	}
+}
+
+void idwt_sym_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, int lpr_len, double *X, int istride, int ostride) {
+	int len_avg, i, l, m, n, t, v;
+	int ms, ns, is;
+	len_avg = lpr_len;
+	m = -2;
+	n = -1;
+
+	for (v = 0; v < len_cA; ++v) {
+		i = v;
+		m += 2;
+		n += 2;
+		ms = m * ostride;
+		ns = n * ostride;
+		X[ms] = 0.0;
+		X[ns] = 0.0;
+		for (l = 0; l < len_avg / 2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				is = (i - l) * istride;
+				X[ms] += lpr[t] * cA[is] + hpr[t] * cD[is];
+				X[ns] += lpr[t + 1] * cA[is] + hpr[t + 1] * cD[is];
+			}
+		}
+	}
+}
+
+void imodwt_per_stride(int M, double *cA, int len_cA, double *cD, double *filt,int lf,double *X,int istride, int ostride) {
+	int len_avg, i, l, t;
+	int is, os;
+	
+	len_avg = lf;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i;
+		os = i * ostride;
+		is = t *istride;
+		X[os] = (filt[0] * cA[is]) + (filt[len_avg] * cD[is]);
+		for (l = 1; l < len_avg; l++) {
+			t += M;
+			while (t >= len_cA) {
+				t -= len_cA;
+			}
+			while (t < 0) {
+				t += len_cA;
+			}
+			is = t *istride;
+			X[os] += (filt[l] * cA[is]) + (filt[len_avg + l] * cD[is]);
+
+		}
+	}
+
+}
+
+void idwt2_shift(int shift, int rows, int cols, double *lpr, double *hpr, int lf, double *A,double *H, double *V,double *D, double *oup) {
+	int i, k, N, ir, ic, J, dim1, dim2;
+	int istride, ostride;
+	double *cL, *cH, *X_lp;
+
+
+	N = rows > cols ? 2 * rows : 2 * cols;
+
+	J = 1;
+	dim1 = 2 * rows;
+	dim2 = 2 * cols;
+
+	X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+	cL = (double*)calloc(dim1*dim2, sizeof(double));
+	cH = (double*)calloc(dim1*dim2, sizeof(double));
+
+	ir = rows;
+	ic = cols;
+	istride = ic;
+	ostride = 1;
+	for (i = 0; i < ic; ++i) {
+		idwt_per_stride(A+i, ir, H+i, lpr, hpr, lf, X_lp, istride, ostride);
+
+		for (k = lf / 2 - 1; k < 2 * ir + lf / 2 - 1; ++k) {
+			cL[(k - lf / 2 + 1)*ic + i] = X_lp[k];
+		}
+
+		idwt_per_stride(V+i, ir, D+i, lpr, hpr, lf, X_lp, istride, ostride);
+
+		for (k = lf / 2 - 1; k < 2 * ir + lf / 2 - 1; ++k) {
+			cH[(k - lf / 2 + 1)*ic + i] = X_lp[k];
+		}
+	}
+	ir *= 2;
+	istride = 1;
+	ostride = 1;
+
+	for (i = 0; i < ir; ++i) {
+		idwt_per_stride(cL + i*ic, ic, cH + i*ic, lpr, hpr, lf, X_lp, istride, ostride);
+
+		for (k = lf / 2 - 1; k < 2 * ic + lf / 2 - 1; ++k) {
+			oup[(k - lf / 2 + 1) + i*ic * 2] = X_lp[k];
+		}
+	}
+	ic *= 2;
+
+
+	if (shift == -1) {
+		//Save the last column
+		for (i = 0; i < ir; ++i) {
+			cL[i] = oup[(i + 1)*ic - 1];
+		}
+		// Save the last row
+		memcpy(cH, oup + (ir - 1)*ic, sizeof(double)*ic);
+		for (i = ir - 1; i > 0; --i) {
+			memcpy(oup + i*ic + 1, oup + (i - 1)*ic, sizeof(double)*(ic - 1));
+		}
+		oup[0] = cL[ir - 1];
+		for (i = 1; i < ir; ++i) {
+			oup[i*ic] = cL[i - 1];
+		}
+
+		for (i = 1; i < ic; ++i) {
+			oup[i] = cH[i - 1];
+		}
+	}
+
+
+	free(X_lp);
+	free(cL);
+	free(cH);
+
+}
+
 int upsamp(double *x, int lenx, int M, double *y) {
 	int N, i, j, k;
 

src/wtmath.h

@@ -10,6 +10,30 @@ Copyright (c) 2014, Rafat Hussain
 extern "C" {
 #endif
 
+void dwt_per_stride(double *inp, int N, double *lpd,double*hpd,int lpd_len,
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void dwt_sym_stride(double *inp, int N, double *lpd, double*hpd, int lpd_len,
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void modwt_per_stride(int M, double *inp, int N, double *filt, int lpd_len, 
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void swt_per_stride(int M, double *inp, int N, double *lpd, double*hpd, int lpd_len, 
+double *cA, int len_cA, double *cD, int istride, int ostride);
+
+void idwt_per_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, 
+int lpr_len, double *X, int istride, int ostride);
+
+void idwt_sym_stride(double *cA, int len_cA, double *cD, double *lpr, double *hpr, 
+int lpr_len, double *X, int istride, int ostride);
+
+void imodwt_per_stride(int M, double *cA, int len_cA, double *cD, double *filt,
+int lf,double *X,int istride, int ostride);
+
+void idwt2_shift(int shift, int rows, int cols, double *lpr, double *hpr, int lf, 
+double *A,double *H, double *V,double *D, double *oup);
+
 int upsamp(double *x, int lenx, int M, double *y);
 
 int upsamp2(double *x, int lenx, int M, double *y);

test/CMakeLists.txt

@@ -30,6 +30,18 @@ add_executable(modwtdenoisetest modwtdenoisetest.c)
 
 target_link_libraries(modwtdenoisetest wauxlib wavelib)
 
+add_executable(dwt2test dwt2test.c)
+
+target_link_libraries(dwt2test wavelib)
+
+add_executable(swt2test swt2test.c)
+
+target_link_libraries(swt2test wavelib)
+
+add_executable(modwt2test modwt2test.c)
+
+target_link_libraries(modwt2test wavelib)
+
 if(UNIX)
 	target_link_libraries(cwttest m)
 	target_link_libraries(dwttest m)
@@ -39,9 +51,12 @@ if(UNIX)
 	target_link_libraries(wtreetest m)
 	target_link_libraries(denoisetest m)
 	target_link_libraries(modwtdenoisetest m)
+	target_link_libraries(dwt2test m)
+	target_link_libraries(swt2test m)
+	target_link_libraries(modwt2test m)
 endif()
 
-set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest modwtdenoisetest
+set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest modwtdenoisetest dwt2test swt2test modwt2test
         PROPERTIES
         RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/test"
         )

test/dwt2test.c

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

test/modwt2test.c

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

test/swt2test.c

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

unitTests/wavelibBoostTests/tst_dwt.cpp

@@ -111,7 +111,15 @@ double REL_Error(double *data, double *rec, int N) {
     return sqrt(sum1)/sqrt(sum2);
 }
 
-void ReconstructionTest()
+double generate_rnd() {
+	double rnd;
+
+	rnd = (double) (rand() % 100 + 1);
+
+	return rnd;
+}
+
+void DWTReconstructionTest()
 {
 
 	wave_object obj;
@@ -229,6 +237,125 @@ void ReconstructionTest()
     free(inp);
 }
 
+void DWT2ReconstructionTest()
+{
+    wave_object obj;
+	wt2_object wt;
+	int i, k, J, N, rows, cols;
+	double *inp, *wavecoeffs,*out;
+    double epsilon;
+
+    rows = 1024;
+    cols = 1000;
+
+    N = rows*cols;
+
+    inp = (double*)malloc(sizeof(double)* N);
+	out = (double*)malloc(sizeof(double)* N);
+
+    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));
+    }
+    
+    waveletNames.push_back("bior1.1");
+    waveletNames.push_back("bior1.3");
+    waveletNames.push_back("bior1.5");
+    waveletNames.push_back("bior2.2");
+    waveletNames.push_back("bior2.4");
+    waveletNames.push_back("bior2.6");
+    waveletNames.push_back("bior2.8");
+    waveletNames.push_back("bior3.1");
+    waveletNames.push_back("bior3.3");
+    waveletNames.push_back("bior3.5");
+    waveletNames.push_back("bior3.7");
+    waveletNames.push_back("bior3.9");
+    waveletNames.push_back("bior4.4");
+    waveletNames.push_back("bior5.5");
+    waveletNames.push_back("bior6.8");
+    
+    waveletNames.push_back("rbior1.1");
+    waveletNames.push_back("rbior1.3");
+    waveletNames.push_back("rbior1.5");
+    waveletNames.push_back("rbior2.2");
+    waveletNames.push_back("rbior2.4");
+    waveletNames.push_back("rbior2.6");
+    waveletNames.push_back("rbior2.8");
+    waveletNames.push_back("rbior3.1");
+    waveletNames.push_back("rbior3.3");
+    waveletNames.push_back("rbior3.5");
+    waveletNames.push_back("rbior3.7");
+    waveletNames.push_back("rbior3.9");
+    waveletNames.push_back("rbior4.4");
+    waveletNames.push_back("rbior5.5");
+    waveletNames.push_back("rbior6.8");
+
+    for (i = 0; i < rows; ++i) {
+		for (k = 0; k < cols; ++k) {
+			//inp[i*cols + k] = i*cols + k;
+			inp[i*cols + k] = generate_rnd();
+			out[i*cols + k] = 0.0;
+		}
+	}
+
+    for (unsigned int direct_fft = 0; direct_fft < 1; 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 = wt2_init(obj,(char*) "dwt", rows,cols, J);// Initialize the wavelet transform object
+                    if (sym_per == 0)
+                        setDWT2Extension(wt, (char*) "sym");// Options are "per" and "sym". Symmetric is the default option
+                    else
+                        setDWT2Extension(wt, (char*) "per");
+
+                    wavecoeffs = dwt2(wt, inp);// Perform DWT
+
+                    idwt2(wt, wavecoeffs, out);// Perform IDWT (if needed)
+                    // Test Reconstruction
+
+                    if (direct_fft == 0)
+                        epsilon = 1e-8;
+                    else
+                        epsilon = 1e-10;
+                    //BOOST_CHECK_SMALL(RMS_Error(out, inp, wt->siglength), epsilon); // If Reconstruction succeeded then the output should be a small value.
+                    
+					//printf("%g ",RMS_Error(out, inp, wt->siglength));
+					if (RMS_Error(out, inp, N) > epsilon) {
+						printf("\n ERROR : DWT2 Reconstruction Unit Test Failed. Exiting. \n");
+						exit(-1);
+					}
+                    wt2_free(wt);
+                    free(wavecoeffs);
+                }
+                wave_free(obj);
+                delete[] name;
+            }
+        }
+    }
+
+    free(inp);
+    free(out);
+}
+
 void MODWTReconstructionTest()
 {
 
@@ -322,6 +449,333 @@ void MODWTReconstructionTest()
     free(inp);
 }
 
+void MODWT2ReconstructionTest()
+{
+    wave_object obj;
+	wt2_object wt;
+	int i, k, J, N, rows, cols;
+	double *inp, *wavecoeffs,*out;
+    double epsilon;
+
+    rows = 1024;
+    cols = 1000;
+
+    N = rows*cols;
+
+    inp = (double*)malloc(sizeof(double)* N);
+	out = (double*)malloc(sizeof(double)* N);
+
+    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 (i = 0; i < rows; ++i) {
+		for (k = 0; k < cols; ++k) {
+			//inp[i*cols + k] = i*cols + k;
+			inp[i*cols + k] = generate_rnd();
+			out[i*cols + k] = 0.0;
+		}
+	}
+
+    for (unsigned int direct_fft = 0; direct_fft < 1; direct_fft++)
+    {
+        for (unsigned int sym_per = 0; sym_per < 1; 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 = wt2_init(obj,(char*) "modwt", rows,cols, J);// Initialize the wavelet transform object
+                    if (sym_per == 0)
+                        setDWT2Extension(wt, (char*) "per");// Options are "per"
+
+                    wavecoeffs = modwt2(wt, inp);// Perform DWT
+
+                    imodwt2(wt, wavecoeffs, out);// Perform IDWT (if needed)
+                    // Test Reconstruction
+
+                    if (direct_fft == 0)
+                        epsilon = 1e-8;
+                    else
+                        epsilon = 1e-10;
+                    //BOOST_CHECK_SMALL(RMS_Error(out, inp, wt->siglength), epsilon); // If Reconstruction succeeded then the output should be a small value.
+                    
+					//printf("%g ",RMS_Error(out, inp, wt->siglength));
+					if (RMS_Error(out, inp, N) > epsilon) {
+						printf("\n ERROR : MODWT2 Reconstruction Unit Test Failed. Exiting. \n");
+						exit(-1);
+					}
+                    wt2_free(wt);
+                    free(wavecoeffs);
+                }
+                wave_free(obj);
+                delete[] name;
+            }
+        }
+    }
+
+    free(inp);
+    free(out);
+}
+
+void SWTReconstructionTest()
+{
+
+    wave_object obj;
+    wt_object wt;
+    double *inp,*out;
+    int N, i,J;
+    double epsilon = 1e-15;
+    double err;
+
+    N = 4000;
+    
+    //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));
+    }
+    
+    waveletNames.push_back("bior1.1");
+    waveletNames.push_back("bior1.3");
+    waveletNames.push_back("bior1.5");
+    waveletNames.push_back("bior2.2");
+    waveletNames.push_back("bior2.4");
+    waveletNames.push_back("bior2.6");
+    waveletNames.push_back("bior2.8");
+    waveletNames.push_back("bior3.1");
+    waveletNames.push_back("bior3.3");
+    waveletNames.push_back("bior3.5");
+    waveletNames.push_back("bior3.7");
+    waveletNames.push_back("bior3.9");
+    waveletNames.push_back("bior4.4");
+    waveletNames.push_back("bior5.5");
+    waveletNames.push_back("bior6.8");
+    
+    waveletNames.push_back("rbior1.1");
+    waveletNames.push_back("rbior1.3");
+    waveletNames.push_back("rbior1.5");
+    waveletNames.push_back("rbior2.2");
+    waveletNames.push_back("rbior2.4");
+    waveletNames.push_back("rbior2.6");
+    waveletNames.push_back("rbior2.8");
+    waveletNames.push_back("rbior3.1");
+    waveletNames.push_back("rbior3.3");
+    waveletNames.push_back("rbior3.5");
+    waveletNames.push_back("rbior3.7");
+    waveletNames.push_back("rbior3.9");
+    waveletNames.push_back("rbior4.4");
+    waveletNames.push_back("rbior5.5");
+    waveletNames.push_back("rbior6.8");
+
+    for (unsigned int direct_fft = 0; direct_fft < 2; direct_fft++)
+    {
+        for (unsigned int sym_per = 0; sym_per < 1; 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*) "swt", 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;
+
+                    swt(wt, inp);// Perform DWT
+
+                    iswt(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 : SWT Reconstruction Unit Test Failed. Exiting. \n");
+                        exit(-1);
+                    }
+                    wt_free(wt);
+                }
+                wave_free(obj);
+                delete[] name;
+            }
+        }
+    }
+    
+    free(out);
+    free(inp);
+}
+
+void SWT2ReconstructionTest()
+{
+    wave_object obj;
+	wt2_object wt;
+	int i, k, J, N, rows, cols;
+	double *inp, *wavecoeffs,*out;
+    double epsilon;
+
+    rows = 1024;
+    cols = 1000;
+
+    N = rows*cols;
+
+    inp = (double*)malloc(sizeof(double)* N);
+	out = (double*)malloc(sizeof(double)* N);
+
+    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));
+    }
+    
+    waveletNames.push_back("bior1.1");
+    waveletNames.push_back("bior1.3");
+    waveletNames.push_back("bior1.5");
+    waveletNames.push_back("bior2.2");
+    waveletNames.push_back("bior2.4");
+    waveletNames.push_back("bior2.6");
+    waveletNames.push_back("bior2.8");
+    waveletNames.push_back("bior3.1");
+    waveletNames.push_back("bior3.3");
+    waveletNames.push_back("bior3.5");
+    waveletNames.push_back("bior3.7");
+    waveletNames.push_back("bior3.9");
+    waveletNames.push_back("bior4.4");
+    waveletNames.push_back("bior5.5");
+    waveletNames.push_back("bior6.8");
+    
+    waveletNames.push_back("rbior1.1");
+    waveletNames.push_back("rbior1.3");
+    waveletNames.push_back("rbior1.5");
+    waveletNames.push_back("rbior2.2");
+    waveletNames.push_back("rbior2.4");
+    waveletNames.push_back("rbior2.6");
+    waveletNames.push_back("rbior2.8");
+    waveletNames.push_back("rbior3.1");
+    waveletNames.push_back("rbior3.3");
+    waveletNames.push_back("rbior3.5");
+    waveletNames.push_back("rbior3.7");
+    waveletNames.push_back("rbior3.9");
+    waveletNames.push_back("rbior4.4");
+    waveletNames.push_back("rbior5.5");
+    waveletNames.push_back("rbior6.8");
+
+    for (i = 0; i < rows; ++i) {
+		for (k = 0; k < cols; ++k) {
+			//inp[i*cols + k] = i*cols + k;
+			inp[i*cols + k] = generate_rnd();
+			out[i*cols + k] = 0.0;
+		}
+	}
+
+    for (unsigned int direct_fft = 0; direct_fft < 1; direct_fft++)
+    {
+        for (unsigned int sym_per = 0; sym_per < 1; 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 = wt2_init(obj,(char*) "swt", rows,cols, J);// Initialize the wavelet transform object
+                    if (sym_per == 0)
+                        setDWT2Extension(wt, (char*) "per");// Options are "per"
+
+                    wavecoeffs = swt2(wt, inp);// Perform DWT
+
+                    iswt2(wt, wavecoeffs, out);// Perform IDWT (if needed)
+                    // Test Reconstruction
+
+                    if (direct_fft == 0)
+                        epsilon = 1e-8;
+                    else
+                        epsilon = 1e-10;
+                    //BOOST_CHECK_SMALL(RMS_Error(out, inp, wt->siglength), epsilon); // If Reconstruction succeeded then the output should be a small value.
+                    
+					//printf("%g ",RMS_Error(out, inp, wt->siglength));
+					if (RMS_Error(out, inp, N) > epsilon) {
+						printf("\n ERROR : SWT2 Reconstruction Unit Test Failed. Exiting. \n");
+						exit(-1);
+					}
+                    wt2_free(wt);
+                    free(wavecoeffs);
+                }
+                wave_free(obj);
+                delete[] name;
+            }
+        }
+    }
+
+    free(inp);
+    free(out);
+}
+
 void DWPTReconstructionTest()
 {
 
@@ -740,10 +1194,13 @@ int main() {
 	RBiorCoefTests();
 	printf("DONE \n");
 	printf("Running DWT ReconstructionTests ... ");
-	ReconstructionTest();
+	DWTReconstructionTest();
 	printf("DONE \n");
     printf("Running MODWT ReconstructionTests ... ");
     MODWTReconstructionTest();
+    printf("DONE \n");
+    printf("Running SWT ReconstructionTests ... ");
+    SWTReconstructionTest();
     printf("DONE \n");
 	printf("Running DWPT ReconstructionTests ... ");
 	DWPTReconstructionTest();
@@ -751,6 +1208,15 @@ int main() {
 	printf("Running CWT ReconstructionTests ... ");
 	CWTReconstructionTest();
 	printf("DONE \n");
+    printf("Running DWT2 ReconstructionTests ... ");
+	DWT2ReconstructionTest();
+	printf("DONE \n");
+    printf("Running MODWT2 ReconstructionTests ... ");
+	MODWT2ReconstructionTest();
+	printf("DONE \n");
+    printf("Running SWT2 ReconstructionTests ... ");
+	SWT2ReconstructionTest();
+	printf("DONE \n");
 	printf("\n\nUnit Tests Successful\n\n");
 	return 0;