diff --git a/header/wavelib.h b/header/wavelib.h
index 4e4f6e9..b4bb293 100755
--- a/header/wavelib.h
+++ b/header/wavelib.h
@@ -191,6 +191,26 @@ 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);
 
@@ -222,6 +242,8 @@ 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);
@@ -246,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);
@@ -254,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);
 
@@ -266,6 +306,8 @@ void wpt_free(wpt_object object);
 
 void cwt_free(cwt_object object);
 
+void wt2_free(wt2_object wt);
+
 
 #ifdef __cplusplus
 }
diff --git a/src/wavelib.c b/src/wavelib.c
index 84c48e7..47eedba 100644
--- a/src/wavelib.c
+++ b/src/wavelib.c
@@ -361,6 +361,87 @@ cwt_object cwt_init(const char* wave, double param,int siglength, double dt, int
 	return obj;
 }
 
+wt2_object wt2_init(wave_object wave, const char* method, int rows, int cols, int J) {
+	int size, i, MaxIter, MaxRows, MaxCols,sumacc;
+	wt2_object obj = NULL;
+
+	size = wave->filtlength;
+
+	MaxRows = wmaxiter(rows, size);
+	MaxCols = wmaxiter(cols, size);
+
+	MaxIter = (MaxRows < MaxCols) ? MaxRows : MaxCols;
+
+	if (J > MaxIter) {
+		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
+		exit(-1);
+	}
+
+	if (J == 1) {
+		sumacc = 4;
+	}
+	else if (J > 1) {
+		sumacc = J * 3 + 1;
+	}
+	else {
+		printf("Error : J should be >= 1 \n");
+		exit(-1);
+	}
+
+	if (method == NULL) {
+		obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
+		obj->outlength = 0; // Default
+		strcpy(obj->ext, "per");
+	}
+	else if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
+		obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
+		obj->outlength = 0; // Default
+		strcpy(obj->ext, "per");
+	}
+	else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) {
+		if (!testSWTlength(rows, J) || !testSWTlength(cols, J)) {
+			printf("\n For SWT data rows and columns must be a multiple of 2^J. \n");
+			exit(-1);
+		}
+
+		obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
+		obj->outlength = 0; // Default
+		strcpy(obj->ext, "per");
+	}
+	else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
+		if (!strstr(wave->wname, "haar")) {
+			if (!strstr(wave->wname, "db")) {
+				if (!strstr(wave->wname, "sym")) {
+					if (!strstr(wave->wname, "coif")) {
+						printf("\n MODWT is only implemented for orthogonal wavelet families - db, sym and coif \n");
+						exit(-1);
+					}
+				}
+			}
+		}
+		obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
+		obj->outlength = 0; // Default
+		strcpy(obj->ext, "per");
+	}
+
+
+	obj->wave = wave;
+	obj->rows = rows;
+	obj->cols = cols;
+	obj->J = J;
+	obj->MaxIter = MaxIter;
+	strcpy(obj->method, method);
+	obj->coeffaccesslength = sumacc;
+
+	obj->dimensions = &obj->params[0];
+	obj->coeffaccess = &obj->params[2 * J];
+	for (i = 0; i < (2 * J + sumacc); ++i) {
+		obj->params[i] = 0;
+	}
+
+	return obj;
+}
+
 
 static void wconv(wt_object wt, double *sig, int N, double *filt, int L, double *oup) {
 	if (!strcmp(wt->cmethod,"direct")) {
@@ -3212,6 +3293,31 @@ void setDWPTExtension(wpt_object wt, const char *extension) {
 	}
 }
 
+void setDWT2Extension(wt2_object wt, const char *extension) {
+	if (!strcmp(wt->method, "dwt")) {
+		if (!strcmp(extension, "sym")) {
+			strcpy(wt->ext, "sym");
+		}
+		else if (!strcmp(extension, "per")) {
+			strcpy(wt->ext, "per");
+		}
+		else {
+			printf("Signal extension can be either per or sym");
+			exit(-1);
+		}
+	}
+	else if (!strcmp(wt->method, "swt") || !strcmp(wt->method, "modwt")) {
+		if (!strcmp(extension, "per")) {
+			strcpy(wt->ext, "per");
+		}
+		else {
+			printf("Signal extension can only be per");
+			exit(-1);
+		}
+	}
+	
+}
+
 void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam) {
 	if (!strcmp(entropy, "shannon")) {
 		strcpy(wt->entropy, "shannon");
@@ -3246,6 +3352,705 @@ void setWTConv(wt_object wt, const char *cmethod) {
 	}
 }
 
+double* dwt2(wt2_object wt, double *inp) {
+	double *wavecoeff;
+	int i, J, iter, N, lp, rows_n, cols_n, rows_i, cols_i;
+	int ir, ic, istride,ostride;
+	int aLL, aLH, aHL, aHH, cdim,clen;
+	double *orig, *lp_dn1,*hp_dn1;
+	J = wt->J;
+	wt->outlength = 0;
+
+	rows_n = wt->rows;
+	cols_n = wt->cols;
+	lp = wt->wave->lpd_len;
+	clen = J * 3;
+	if (!strcmp(wt->ext, "per")) {
+		i = 2 * J;
+		while (i > 0) {
+			rows_n = (int)ceil((double)rows_n / 2.0);
+			cols_n = (int)ceil((double)cols_n / 2.0);
+			wt->dimensions[i - 1] = cols_n;
+			wt->dimensions[i - 2] = rows_n;
+			wt->outlength += (rows_n * cols_n) * 3;
+			i = i - 2;
+		}
+		wt->outlength += (rows_n * cols_n);
+		N = wt->outlength;
+		wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
+
+		orig = inp;
+		ir = wt->rows;
+		ic = wt->cols;
+		cols_i = wt->dimensions[2 * J - 1];
+
+		lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+		hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+
+		for (iter = 0; iter < J; ++iter) {
+			rows_i = wt->dimensions[2*J - 2*iter - 2];
+			cols_i = wt->dimensions[2*J - 2*iter - 1];
+			istride = 1;
+			ostride = 1;
+			cdim = rows_i * cols_i;
+			// Row filtering and column subsampling
+			for (i = 0; i < ir; ++i) {
+				dwt_per_stride(orig+i*ic, ic, wt->wave->lpd,wt->wave->hpd,lp, lp_dn1+i*cols_i, cols_i, hp_dn1+i*cols_i, istride, ostride);
+			}
+
+			// Column Filtering and Row subsampling
+			aHH = N - cdim;
+			wt->coeffaccess[clen] = aHH;
+			aHL = aHH - cdim;
+			wt->coeffaccess[clen-1] = aHL;
+			aLH = aHL - cdim;
+			wt->coeffaccess[clen-2] = aLH;
+			aLL = aLH - cdim;
+			
+			N -= 3 * cdim;
+			ic = cols_i;
+			istride = ic;
+			ostride = ic;
+			
+			for (i = 0; i < ic; ++i) {
+				dwt_per_stride(lp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff+aLL+i, rows_i, wavecoeff+aLH+i, istride, ostride);
+			}
+
+
+			for (i = 0; i < ic; ++i) {
+				dwt_per_stride(hp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff+aHL+i, rows_i, wavecoeff+aHH+i, istride, ostride);
+			}
+
+			ir = rows_i;
+			orig = wavecoeff+aLL;
+			clen -= 3;
+
+		}
+		wt->coeffaccess[0] = 0;
+		free(lp_dn1);
+		free(hp_dn1);
+	}
+	else if (!strcmp(wt->ext, "sym")) {
+		i = 2 * J;
+		while (i > 0) {
+			rows_n += lp - 2;
+			cols_n += lp - 2;
+			rows_n = (int)ceil((double)rows_n / 2.0);
+			cols_n = (int)ceil((double)cols_n / 2.0);
+			wt->dimensions[i - 1] = cols_n;
+			wt->dimensions[i - 2] = rows_n;
+			wt->outlength += (rows_n * cols_n) * 3;
+			i = i - 2;
+		}
+		wt->outlength += (rows_n * cols_n);
+		N = wt->outlength;
+		wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
+
+		orig = inp;
+		ir = wt->rows;
+		ic = wt->cols;
+		cols_i = wt->dimensions[2 * J - 1];
+
+		lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+		hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+
+		for (iter = 0; iter < J; ++iter) {
+			rows_i = wt->dimensions[2 * J - 2 * iter - 2];
+			cols_i = wt->dimensions[2 * J - 2 * iter - 1];
+			istride = 1;
+			ostride = 1;
+			cdim = rows_i * cols_i;
+			// Row filtering and column subsampling
+			for (i = 0; i < ir; ++i) {
+				dwt_sym_stride(orig + i*ic, ic, wt->wave->lpd, wt->wave->hpd, lp, lp_dn1 + i*cols_i, cols_i, hp_dn1 + i*cols_i, istride, ostride);
+			}
+
+			// Column Filtering and Row subsampling
+			aHH = N - cdim;
+			wt->coeffaccess[clen] = aHH;
+			aHL = aHH - cdim;
+			wt->coeffaccess[clen - 1] = aHL;
+			aLH = aHL - cdim;
+			wt->coeffaccess[clen - 2] = aLH;
+			aLL = aLH - cdim;
+			N -= 3 * cdim;
+			ic = cols_i;
+			istride = ic;
+			ostride = ic;
+
+			for (i = 0; i < ic; ++i) {
+				dwt_sym_stride(lp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aLL + i, rows_i, wavecoeff + aLH + i, istride, ostride);
+			}
+
+			for (i = 0; i < ic; ++i) {
+				dwt_sym_stride(hp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aHL + i, rows_i, wavecoeff + aHH + i, istride, ostride);
+			}
+
+			ir = rows_i;
+			orig = wavecoeff + aLL;
+			clen -= 3;
+
+		}
+		wt->coeffaccess[0] = 0;
+		free(lp_dn1);
+		free(hp_dn1);
+	}
+
+	return wavecoeff;
+}
+
+void idwt2(wt2_object wt, double *wavecoeff, double *oup) {
+	int i, k, rows, cols, N, ir,ic,lf,dim1,dim2;
+	int istride, ostride, iter, J;
+	int aLL, aLH, aHL, aHH;
+	double *cL, *cH, *X_lp,*orig;
+
+	rows = wt->rows;
+	cols = wt->cols;
+	J = wt->J;
+	double *out;
+	
+
+	if (!strcmp(wt->ext, "per")) {
+		N = rows > cols ? 2 * rows : 2 * cols;
+		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+
+		i = J;
+		dim1 = wt->dimensions[0];
+		dim2 = wt->dimensions[1];
+		k = 0;
+		while (i > 0) {
+			k += 1;
+			dim1 *= 2;
+			dim2 *= 2;
+			i--;
+		}
+
+		
+
+		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+		cL = (double*)calloc(dim1*dim2, sizeof(double));
+		cH = (double*)calloc(dim1*dim2, sizeof(double));
+		out = (double*)calloc(dim1*dim2, sizeof(double));
+		aLL = wt->coeffaccess[0];
+		orig = wavecoeff + aLL;
+		for (iter = 0; iter < J; ++iter) {
+			ir = wt->dimensions[2 * iter];
+			ic = wt->dimensions[2 * iter + 1];
+			istride = ic;
+			ostride = 1;
+			aLH = wt->coeffaccess[iter*3 + 1];
+			aHL = wt->coeffaccess[iter*3 + 2];
+			aHH = wt->coeffaccess[iter*3 + 3];
+			for (i = 0; i < ic; ++i) {
+				idwt_per_stride(orig+i, ir, wavecoeff+aLH+i, wt->wave->lpr, wt->wave->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(wavecoeff + aHL+i, ir, wavecoeff + aHH+i, wt->wave->lpr, wt->wave->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, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
+
+				for (k = lf / 2 - 1; k < 2 * ic + lf / 2 - 1; ++k) {
+					out[(k - lf / 2 + 1) + i*ic*2] = X_lp[k];
+				}
+			}
+			ic *= 2;
+			if (iter == J - 1) {
+				for (i = 0; i < wt->rows; ++i) {
+					for (k = 0; k < wt->cols; ++k) {
+						oup[k + i*wt->cols] = out[k + i*ic];
+					}
+				}
+			}
+			else {
+				for (i = 0; i < wt->dimensions[2 * (iter+1)]; ++i) {
+					for (k = 0; k < wt->dimensions[2 * (iter + 1)+1]; ++k) {
+						oup[k + i*wt->dimensions[2 * (iter + 1) + 1]] = out[k + i*ic];
+					}
+				}
+			}
+			
+
+			orig = oup;
+		}
+		free(X_lp);
+		free(cL);
+		free(cH);
+	}
+	else if (!strcmp(wt->ext, "sym")) {
+		N = rows > cols ? 2 * rows - 1 : 2 * cols - 1;
+		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+
+		i = J;
+		dim1 = wt->dimensions[0];
+		dim2 = wt->dimensions[1];
+		k = 0;
+		while (i > 0) {
+			k += 1;
+			dim1 *= 2;
+			dim2 *= 2;
+			i--;
+		}
+
+
+
+		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+		cL = (double*)calloc(dim1*dim2, sizeof(double));
+		cH = (double*)calloc(dim1*dim2, sizeof(double));
+		out = (double*)calloc(dim1*dim2, sizeof(double));
+		aLL = wt->coeffaccess[0];
+		orig = wavecoeff + aLL;
+		for (iter = 0; iter < J; ++iter) {
+			ir = wt->dimensions[2 * iter];
+			ic = wt->dimensions[2 * iter + 1];
+			istride = ic;
+			ostride = 1;
+			aLH = wt->coeffaccess[iter * 3 + 1];
+			aHL = wt->coeffaccess[iter * 3 + 2];
+			aHH = wt->coeffaccess[iter * 3 + 3];
+			for (i = 0; i < ic; ++i) {
+				idwt_sym_stride(orig + i, ir, wavecoeff + aLH + i, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
+
+				for (k = lf - 2; k < 2 * ir; ++k) {
+					cL[(k - lf + 2)*ic + i] = X_lp[k];
+				}
+
+				idwt_per_stride(wavecoeff + aHL + i, ir, wavecoeff + aHH + i, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
+
+				for (k = lf - 2; k < 2 * ir; ++k) {
+					cH[(k - lf + 2)*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, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
+
+				for (k = lf - 2; k < 2 * ic; ++k) {
+					out[(k - lf + 2) + i*ic * 2] = X_lp[k];
+				}
+			}
+			ic *= 2;
+			if (iter == J - 1) {
+				for (i = 0; i < wt->rows; ++i) {
+					for (k = 0; k < wt->cols; ++k) {
+						oup[k + i*wt->cols] = out[k + i*ic];
+					}
+				}
+			}
+			else {
+				for (i = 0; i < wt->dimensions[2 * (iter + 1)]; ++i) {
+					for (k = 0; k < wt->dimensions[2 * (iter + 1) + 1]; ++k) {
+						oup[k + i*wt->dimensions[2 * (iter + 1) + 1]] = out[k + i*ic];
+					}
+				}
+			}
+
+
+			orig = oup;
+		}
+		free(X_lp);
+		free(cL);
+		free(cH);
+	}
+
+	free(out);
+	
+}
+
+double* swt2(wt2_object wt, double *inp) {
+	double *wavecoeff;
+	int i, J, iter, M, N, lp, rows_n, cols_n, rows_i, cols_i;
+	int ir, ic, istride, ostride;
+	int aLL, aLH, aHL, aHH, cdim, clen;
+	double *orig, *lp_dn1, *hp_dn1;
+
+	J = wt->J;
+	M = 1;
+	wt->outlength = 0;
+
+	rows_n = wt->rows;
+	cols_n = wt->cols;
+	lp = wt->wave->lpd_len;
+	clen = J * 3;
+
+	i = 2 * J;
+	while (i > 0) {
+		wt->dimensions[i - 1] = cols_n;
+		wt->dimensions[i - 2] = rows_n;
+		wt->outlength += (rows_n * cols_n) * 3;
+		i = i - 2;
+	}
+	wt->outlength += (rows_n * cols_n);
+	N = wt->outlength;
+	wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
+
+	orig = inp;
+	ir = wt->rows;
+	ic = wt->cols;
+	cols_i = wt->dimensions[2 * J - 1];
+
+	lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+	hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+
+	for (iter = 0; iter < J; ++iter) {
+		if (iter > 0) {
+			M = 2 * M;
+		}
+		rows_i = wt->dimensions[2 * J - 2 * iter - 2];
+		cols_i = wt->dimensions[2 * J - 2 * iter - 1];
+		istride = 1;
+		ostride = 1;
+		cdim = rows_i * cols_i;
+		// Row filtering and column subsampling
+		for (i = 0; i < ir; ++i) {
+			swt_per_stride(M,orig + i*ic, ic, wt->wave->lpd, wt->wave->hpd, lp, lp_dn1 + i*cols_i, cols_i, hp_dn1 + i*cols_i, istride, ostride);
+		}
+		// Column Filtering and Row subsampling
+		aHH = N - cdim;
+		wt->coeffaccess[clen] = aHH;
+		aHL = aHH - cdim;
+		wt->coeffaccess[clen - 1] = aHL;
+		aLH = aHL - cdim;
+		wt->coeffaccess[clen - 2] = aLH;
+		aLL = aLH - cdim;
+
+		N -= 3 * cdim;
+		ic = cols_i;
+		istride = ic;
+		ostride = ic;
+		for (i = 0; i < ic; ++i) {
+			swt_per_stride(M,lp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aLL + i, rows_i, wavecoeff + aLH + i, istride, ostride);
+		}
+
+		for (i = 0; i < ic; ++i) {
+			swt_per_stride(M,hp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aHL + i, rows_i, wavecoeff + aHH + i, istride, ostride);
+		}
+
+		ir = rows_i;
+		orig = wavecoeff + aLL;
+		clen -= 3;
+	}
+	wt->coeffaccess[0] = 0;
+	free(lp_dn1);
+	free(hp_dn1);
+
+	return wavecoeff;
+}
+
+void iswt2(wt2_object wt, double *wavecoeffs, double *oup) {
+	int i, k, iter, it2, it3, J, M, rows, cols, lf, ir,ic,k1,i1;
+	double *A, *H, *V, *D,*oup1,*oup2;
+	int aLL, aLH, aHL, aHH,shift;
+	J = wt->J;
+	rows = wt->rows;
+	cols = wt->cols;
+	lf = wt->wave->lpd_len;
+	A = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
+	H = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
+	V = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
+	D = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
+	oup1 = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
+	oup2 = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
+
+	aLL = wt->coeffaccess[0];
+
+	for (i = 0; i < rows; ++i) {
+		for (k = 0; k < cols; ++k) {
+			oup[i*cols + k] = wavecoeffs[aLL + i*cols + k];
+		}
+	}
+
+	for (iter = J; iter > 0; iter--) {
+		aLH = wt->coeffaccess[(J - iter) * 3 + 1];
+		aHL = wt->coeffaccess[(J - iter) * 3 + 2];
+		aHH = wt->coeffaccess[(J - iter) * 3 + 3];
+		M = (int)pow(2.0, (double)iter - 1);
+
+		for (it2 = 0; it2 < M; ++it2) {
+			ir = 0;
+			ic = 0;
+			it3 = 0;
+			// oup1
+			for (i = it2; i < rows; i += 2 * M) {
+				ic = 0;
+				for (k = it2; k < cols; k += 2 * M) {
+					A[it3] = oup[i*cols + k];
+					H[it3] = wavecoeffs[aLH + i*cols + k];
+					V[it3] = wavecoeffs[aHL + i*cols + k];
+					D[it3] = wavecoeffs[aHH + i*cols + k];
+					it3++;
+					ic++;
+				}
+				ir++;
+			}
+			shift = 0;
+			idwt2_shift(shift, ir, ic, wt->wave->lpr, wt->wave->hpr, wt->wave->lpd_len, A, H, V, D, oup1);
+			//oup2
+			ir = 0;
+			ic = 0;
+			it3 = 0;
+			for (i = it2 + M; i < rows; i += 2 * M) {
+				ic = 0;
+				for (k = it2 + M; k < cols; k += 2 * M) {
+					A[it3] = oup[i*cols + k];
+					H[it3] = wavecoeffs[aLH + i*cols + k];
+					V[it3] = wavecoeffs[aHL + i*cols + k];
+					D[it3] = wavecoeffs[aHH + i*cols + k];
+					it3++;
+					ic++;
+				}
+				ir++;
+			}
+			shift = -1;
+			idwt2_shift(shift, ir, ic, wt->wave->lpr, wt->wave->hpr, wt->wave->lpd_len, A, H, V, D, oup2);
+			// Shift oup1 and oup2. Then add them to get A.
+			i1 = 0;
+			for (i = it2; i < rows; i += M) {
+				k1 = 0;
+				for (k = it2; k < cols; k += M) {
+					oup[i*cols+k] = 0.5*(oup1[i1*2*ic + k1] + oup2[i1*2*ic+k1]);
+					k1++;
+				}
+				i1++;
+			}
+		}
+
+	}
+
+	free(A);
+	free(H);
+	free(V);
+	free(D);
+	free(oup1);
+	free(oup2);
+}
+
+double* modwt2(wt2_object wt, double *inp) {
+	double *wavecoeff;
+	int i, J, iter, M, N, lp, rows_n, cols_n, rows_i, cols_i;
+	int ir, ic, istride, ostride;
+	int aLL, aLH, aHL, aHH, cdim, clen;
+	double *orig, *lp_dn1, *hp_dn1,*filt;
+	double s;
+
+	J = wt->J;
+	M = 1;
+	wt->outlength = 0;
+
+	rows_n = wt->rows;
+	cols_n = wt->cols;
+	lp = wt->wave->lpd_len;
+	clen = J * 3;
+
+	i = 2 * J;
+	while (i > 0) {
+		wt->dimensions[i - 1] = cols_n;
+		wt->dimensions[i - 2] = rows_n;
+		wt->outlength += (rows_n * cols_n) * 3;
+		i = i - 2;
+	}
+	wt->outlength += (rows_n * cols_n);
+	N = wt->outlength;
+	wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
+	filt = (double*)malloc(sizeof(double)* 2 * lp);
+	s = sqrt(2.0);
+	for (i = 0; i < lp; ++i) {
+		filt[i] = wt->wave->lpd[i] / s;
+		filt[lp + i] = wt->wave->hpd[i] / s;
+	}
+
+	orig = inp;
+	ir = wt->rows;
+	ic = wt->cols;
+	cols_i = wt->dimensions[2 * J - 1];
+
+	lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+	hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
+
+	for (iter = 0; iter < J; ++iter) {
+		if (iter > 0) {
+			M = 2 * M;
+		}
+		rows_i = wt->dimensions[2 * J - 2 * iter - 2];
+		cols_i = wt->dimensions[2 * J - 2 * iter - 1];
+		istride = 1;
+		ostride = 1;
+		cdim = rows_i * cols_i;
+		// Row filtering and column subsampling
+		for (i = 0; i < ir; ++i) {
+			modwt_per_stride(M, orig + i*ic, ic, filt, lp, lp_dn1 + i*cols_i, cols_i, hp_dn1 + i*cols_i, istride, ostride);
+		}
+		// Column Filtering and Row subsampling
+		aHH = N - cdim;
+		wt->coeffaccess[clen] = aHH;
+		aHL = aHH - cdim;
+		wt->coeffaccess[clen - 1] = aHL;
+		aLH = aHL - cdim;
+		wt->coeffaccess[clen - 2] = aLH;
+		aLL = aLH - cdim;
+		N -= 3 * cdim;
+		ic = cols_i;
+		istride = ic;
+		ostride = ic;
+		for (i = 0; i < ic; ++i) {
+			modwt_per_stride(M, lp_dn1 + i, ir, filt, lp, wavecoeff + aLL + i, rows_i, wavecoeff + aLH + i, istride, ostride);
+		}
+
+
+		for (i = 0; i < ic; ++i) {
+			modwt_per_stride(M, hp_dn1 + i, ir, filt, lp, wavecoeff + aHL + i, rows_i, wavecoeff + aHH + i, istride, ostride);
+		}
+
+
+		ir = rows_i;
+		orig = wavecoeff + aLL;
+		clen -= 3;
+	}
+	wt->coeffaccess[0] = 0;
+	free(lp_dn1);
+	free(hp_dn1);
+	free(filt);
+	return wavecoeff;
+}
+
+void imodwt2(wt2_object wt, double *wavecoeff, double *oup) {
+	int i, rows, cols, M, N, ir, ic, lf;
+	int istride, ostride, iter, J;
+	int aLL, aLH, aHL, aHH;
+	double *cL, *cH, *orig,*filt;
+	double s;
+
+	rows = wt->rows;
+	cols = wt->cols;
+	J = wt->J;
+
+
+	M = (int)pow(2.0, (double)J - 1.0);
+	N = rows > cols ? rows : cols;
+	lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+
+	filt = (double*)malloc(sizeof(double)* 2 * lf);
+	s = sqrt(2.0);
+	for (i = 0; i < lf; ++i) {
+		filt[i] = wt->wave->lpd[i] / s;
+		filt[lf + i] = wt->wave->hpd[i] / s;
+	}
+
+
+	cL = (double*)calloc(rows*cols, sizeof(double));
+	cH = (double*)calloc(rows*cols, sizeof(double));
+	aLL = wt->coeffaccess[0];
+	orig = wavecoeff + aLL;
+	for (iter = 0; iter < J; ++iter) {
+		if (iter > 0) {
+			M = M / 2;
+		}
+		ir = wt->dimensions[2 * iter];
+		ic = wt->dimensions[2 * iter + 1];
+		istride = ic;
+		ostride = ic;
+		aLH = wt->coeffaccess[iter * 3 + 1];
+		aHL = wt->coeffaccess[iter * 3 + 2];
+		aHH = wt->coeffaccess[iter * 3 + 3];
+		for (i = 0; i < ic; ++i) {
+			imodwt_per_stride(M,orig + i, ir, wavecoeff + aLH + i, filt, lf, cL + i, istride, ostride);
+
+			imodwt_per_stride(M,wavecoeff + aHL + i, ir, wavecoeff + aHH + i, filt, lf, cH + i, istride, ostride);
+
+		}
+
+		istride = 1;
+		ostride = 1;
+
+		for (i = 0; i < ir; ++i) {
+			imodwt_per_stride(M,cL + i*ic, ic, cH + i*ic, filt, lf, oup+i*ic, istride, ostride);
+
+		}
+
+		orig = oup;
+	}
+
+	free(cL);
+	free(cH);
+	free(filt);
+}
+
+double* getWT2Coeffs(wt2_object wt,double* wcoeffs, int level,char *type, int *rows, int *cols) {
+	int J,iter,t;
+	double *ptr;
+	J = wt->J;
+	// Error Check
+
+	if (level > J || level < 1) {
+		printf("Error : The data is decomposed into %d levels so the acceptable values of level are between 1 and %d", J, J);
+		exit(-1);
+	}
+
+	if (!strcmp(type, "A") && level != J) {
+		printf("Approximation Coefficients are only available for level %d", J);
+		exit(-1);
+	}
+
+	if (!strcmp(type, "A")) {
+		t = 0;
+		iter = t;
+	}
+	else if (!strcmp(type, "H")) {
+		t = 1;
+		iter = t;
+	}
+	else if (!strcmp(type, "V")) {
+		t = 2;
+		iter = t;
+	}
+	else if (!strcmp(type, "D")) {
+		t = 3;
+		iter = t;
+	}
+	else {
+		printf("Only four types of coefficients are accessible A, H, V and D \n");
+		exit(-1);
+	}
+
+	iter += (J - level) * 3;
+
+	ptr = wcoeffs+wt->coeffaccess[iter];
+	*rows = wt->dimensions[2 * (J - level)];
+	*cols = wt->dimensions[2 * (J - level)+1];
+
+	return ptr;
+}
+
+void dispWT2Coeffs(double *A, int row, int col) {
+	int i, j;
+	printf("\n MATRIX Order : %d X %d \n \n", row, col);
+
+	for (i = 0; i < row; i++) {
+		printf("R%d: ", i);
+		for (j = 0; j < col; j++) {
+			printf("%g ", A[i*col + j]);
+		}
+		printf(":R%d \n", i);
+	}
+}
+
 void wave_summary(wave_object obj) {
 	int i,N;
 	N = obj->filtlength;
@@ -3409,6 +4214,47 @@ void cwt_summary(cwt_object wt) {
 
 }
 
+void wt2_summary(wt2_object wt) {
+	int i;
+	int J, t,rows,cols,vsize;
+	J = wt->J;
+	wave_summary(wt->wave);
+	printf("\n");
+	printf("Wavelet Transform : %s \n", wt->method);
+	printf("\n");
+	printf("Signal Extension : %s \n", wt->ext);
+	printf("\n");
+	printf("Number of Decomposition Levels %d \n", wt->J);
+	printf("\n");
+	printf("Input Signal Rows %d \n", wt->rows);
+	printf("\n");
+	printf("Input Signal Cols %d \n", wt->cols);
+	printf("\n");
+	printf("Length of Wavelet Coefficients Vector %d \n", wt->outlength);
+	printf("\n");
+	t = 0;
+	for (i = J; i > 0; --i) {
+		rows = wt->dimensions[2 * (J - i)];
+		cols = wt->dimensions[2 * (J - i) + 1];
+		vsize = rows*cols;
+		printf("Level %d Decomposition Rows :%d Columns:%d Vector Size (Rows*Cols):%d \n", i, rows, cols,vsize);
+		printf("Access Row values stored at wt->dimensions[%d]\n", 2 * (J - i));
+		printf("Access Column values stored at wt->dimensions[%d]\n\n", 2 * (J - i) + 1);
+
+		if (i == J) {
+			printf("Approximation Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n", t,wt->coeffaccess[t],vsize);
+		}
+
+		t += 1;
+		printf("Horizontal Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n", t, wt->coeffaccess[t], vsize);
+		t += 1;
+		printf("Vertical Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n", t, wt->coeffaccess[t], vsize);
+		t += 1;
+		printf("Detail Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n\n", t, wt->coeffaccess[t], vsize);
+	}
+
+}
+
 void wave_free(wave_object object) {
 	free(object);
 }
@@ -3428,3 +4274,7 @@ void wpt_free(wpt_object object) {
 void cwt_free(cwt_object object) {
 	free(object);
 }
+
+void wt2_free(wt2_object wt) {
+	free(wt);
+}
diff --git a/src/wtmath.c b/src/wtmath.c
index 6685704..575df0c 100644
--- a/src/wtmath.c
+++ b/src/wtmath.c
@@ -276,6 +276,79 @@ void imodwt_per_stride(int M, double *cA, int len_cA, double *cD, double *filt,i
 
 }
 
+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;
diff --git a/src/wtmath.h b/src/wtmath.h
index 6a38eb9..224d41e 100644
--- a/src/wtmath.h
+++ b/src/wtmath.h
@@ -31,6 +31,9 @@ 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);