diff --git a/README.md b/README.md
index cb891f2..844ec85 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 wavelib
 =======
 
-C Implementation of Wavelet Transform (DWT,SWT and MODWT) and Packet Transform ( Full Tree Decomposition and Best Basis DPWT).
+C Implementation of Wavelet Transform (DWT,SWT and MODWT) and Packet Transform ( Full Tree Decomposition and Best Basis DWPT).
 
 Discrete Wavelet Transform Methods Implemented
 
diff --git a/header/wavelib.h b/header/wavelib.h
index 62a9f67..2399d8c 100644
--- a/header/wavelib.h
+++ b/header/wavelib.h
@@ -14,6 +14,16 @@ extern "C" {
 #define fft_type double
 #endif
 
+#ifndef cplx_type
+#define cplx_type double
+#endif
+
+
+typedef struct cplx_t {
+	cplx_type re;
+	cplx_type im;
+} cplx_data;
+
 typedef struct wave_set* wave_object;
 
 wave_object wave_init(char* wname);
@@ -153,6 +163,34 @@ struct wpt_set{
 };
 
 
+typedef struct cwt_set* cwt_object;
+
+cwt_object cwt_init(char* wave, double param, int siglength,double dt, int J);
+
+struct cwt_set{
+	char wave[10];// Wavelet - morl/morlet,paul,dog/dgauss
+	int siglength;// Length of Input Data
+	int J;// Total Number of Scales
+	double s0;// Smallest scale. It depends on the sampling rate. s0 <= 2 * dt for most wavelets
+	double dt;// Sampling Rate
+	double dj;// Separation between scales. eg., scale = s0 * 2 ^ ( [0:N-1] *dj ) or scale = s0 *[0:N-1] * dj
+	char type[10];// Scale Type - Power or Linear
+	int pow;// Base of Power in case type = pow. Typical value is pow = 2
+	int sflag;
+	int pflag;
+	int npad;
+	int mother;
+	double m;// Wavelet parameter param
+	double smean;// Input Signal mean
+
+	cplx_data *output;
+	double *scale;
+	double *period;
+	double *coi;
+	double params[0];
+};
+
+
 void dwt(wt_object wt, double *inp);
 
 void idwt(wt_object wt, double *dwtop);
@@ -189,6 +227,18 @@ int getDWPTNodelength(wpt_object wt, int X);
 
 void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N);
 
+void setCWTScales(cwt_object wt, double s0, double dj, char *type, int power);
+
+void setCWTScaleVector(cwt_object wt, double *scale, int J, double s0, double dj);
+
+void setCWTPadding(cwt_object wt, int pad);
+
+void cwt(cwt_object wt, double *inp);
+
+void icwt(cwt_object wt, double *cwtop);
+
+int getCWTScaleLength(int N);
+
 void wave_summary(wave_object obj);
 
 void wt_summary(wt_object wt);
@@ -197,6 +247,8 @@ void wtree_summary(wtree_object wt);
 
 void wpt_summary(wpt_object wt);
 
+void cwt_summary(cwt_object wt);;
+
 void wave_free(wave_object object);
 
 void wt_free(wt_object object);
@@ -205,6 +257,8 @@ void wtree_free(wtree_object object);
 
 void wpt_free(wpt_object object);
 
+void cwt_free(cwt_object object);
+
 
 #ifdef __cplusplus
 }
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 58f085e..758f5c3 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,17 +3,23 @@
 include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 
 set(SOURCE_FILES    conv.c
+					cwt.c
+					cwtmath.c
 					hsfft.c
 					real.c
 					wavefilt.c
+					wavefunc.c
 					wavelib.c
 					wtmath.c
                     )
 
 set(HEADER_FILES    conv.h
+					cwt.h
+					cwtmath.h
 					hsfft.h
 					real.h
 					wavefilt.h
+					wavefunc.h
 					wavelib.h
 					wtmath.h
                     )
diff --git a/src/cwt.c b/src/cwt.c
new file mode 100755
index 0000000..6e9efb1
--- /dev/null
+++ b/src/cwt.c
@@ -0,0 +1,420 @@
+/*
+  Copyright (c) 2015, Rafat Hussain
+*/
+/*
+This code is a C translation ( with some modifications) of Wavelet Software provided by
+C. Torrence and G. Compo, and is available at URL: http://atoc.colorado.edu/research/wavelets/''.
+*/
+
+#include "cwt.h"
+
+static int factorial2(int N) {
+	int factorial,i;
+
+	factorial = 1;
+
+	for (i = 1; i <= N;++i) {
+		factorial *= i;
+	}
+
+	return factorial;
+}
+
+static double factorial3(int N) {
+	int i;
+	double factorial;
+
+	factorial = 1;
+
+	for (i = 1; i <= N; ++i) {
+		factorial *= i;
+	}
+
+	return factorial;
+}
+
+double factorial(int N) {
+	if (N > 40) {
+		printf("This program is only valid for N <= 40 \n");
+		return -1.0;
+	}
+	double fact[41] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000,
+		20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000, 51090942171709440000.0, 1124000727777607680000.0,
+		25852016738884976640000.0, 620448401733239439360000.0, 15511210043330985984000000.0, 403291461126605635584000000.0, 10888869450418352160768000000.0,
+		304888344611713860501504000000.0, 8841761993739701954543616000000.0, 265252859812191058636308480000000.0, 8222838654177922817725562880000000.0,
+		263130836933693530167218012160000000.0, 8683317618811886495518194401280000000.0, 295232799039604140847618609643520000000.0, 10333147966386144929666651337523200000000.0,
+		371993326789901217467999448150835200000000.0, 13763753091226345046315979581580902400000000.0, 523022617466601111760007224100074291200000000.0,
+		20397882081197443358640281739902897356800000000.0, 815915283247897734345611269596115894272000000000.0 };
+
+	return fact[N];
+
+}
+static void wave_function(int nk, double dt,int mother, double param,double scale1, double *kwave, double pi,double *period1,
+	double *coi1, fft_data *daughter) {
+
+	double norm, expnt, fourier_factor;
+	int k, m;
+	double temp;
+	int sign,re;
+
+
+	if (mother == 0) {
+		//MORLET
+		if (param < 0.0) {
+			param = 6.0;
+		}
+		norm = sqrt(2.0*pi*scale1 / dt)*pow(pi,-0.25);
+
+		for (k = 1; k <= nk / 2 + 1; ++k) {
+			temp = (scale1*kwave[k-1] - param);
+			expnt = -0.5 * temp * temp;
+			daughter[k - 1].re = norm * exp(expnt);
+			daughter[k - 1].im = 0.0;
+		}
+		for (k = nk / 2 + 2; k <= nk; ++k) {
+			daughter[k - 1].re = daughter[k - 1].im = 0.0;
+		}
+		fourier_factor = (4.0*pi) / (param + sqrt(2.0 + param*param));
+		*period1 = scale1*fourier_factor;
+		*coi1 = fourier_factor / sqrt(2.0);
+	}
+	else if (mother == 1) {
+		// PAUL
+		if (param < 0.0) {
+			param = 4.0;
+		}
+		m = (int)param;
+		norm = sqrt(2.0*pi*scale1 / dt)*(pow(2.0,(double)m) / sqrt((double)(m*factorial(2 * m - 1))));
+		for (k = 1; k <= nk / 2 + 1; ++k) {
+			temp = scale1 * kwave[k - 1];
+			expnt = - temp;
+			daughter[k - 1].re = norm * pow(temp,(double)m) * exp(expnt);
+			daughter[k - 1].im = 0.0;
+		}
+		for (k = nk / 2 + 2; k <= nk; ++k) {
+			daughter[k - 1].re = daughter[k - 1].im = 0.0;
+		}
+		fourier_factor = (4.0*pi) / (2.0 * m + 1.0);
+		*period1 = scale1*fourier_factor;
+		*coi1 = fourier_factor * sqrt(2.0);
+	}
+	else if (mother == 2) {
+		if (param < 0.0) {
+			param = 2.0;
+		}
+		m = (int)param;
+
+		if (m % 2 == 0) {
+			re = 1;
+		}
+		else {
+			re = 0;
+		}
+
+		if (m % 4 == 0 || m % 4 == 1) {
+			sign = -1;
+		}
+		else {
+			sign = 1;
+		}
+
+
+		norm = sqrt(2.0*pi*scale1 / dt)*sqrt(1.0 / gamma(m + 0.50));
+		norm *= sign;
+
+		if (re == 1) {
+			for (k = 1; k <= nk; ++k) {
+				temp = scale1 * kwave[k - 1];
+				daughter[k - 1].re = norm*pow(temp,(double)m)*exp(-0.50*pow(temp,2.0));
+				daughter[k - 1].im = 0.0;
+			}
+		}
+		else if (re == 0) {
+			for (k = 1; k <= nk; ++k) {
+				temp = scale1 * kwave[k - 1];
+				daughter[k - 1].re = 0.0;
+				daughter[k - 1].im = norm*pow(temp, (double)m)*exp(-0.50*pow(temp, 2.0));
+			}
+		}
+		fourier_factor = (2.0*pi) * sqrt(2.0 / (2.0 * m + 1.0));
+		*period1 = scale1*fourier_factor;
+		*coi1 = fourier_factor / sqrt(2.0);
+	}
+}
+
+void cwavelet(double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad,
+	double *wave, double *scale, double *period, double *coi) {
+
+	int i, j, k, iter;
+	double ymean, freq1, pi, period1, coi1;
+	double tmp1, tmp2;
+	double scale1;
+	double *kwave;
+	fft_object obj, iobj;
+	fft_data *ypad, *yfft,*daughter;
+
+	pi = 4.0 * atan(1.0);
+
+	if (npad < N) {
+		printf("npad must be >= N \n");
+		exit(-1);
+	}
+
+	obj = fft_init(npad, 1);
+	iobj = fft_init(npad, -1);
+
+	ypad = (fft_data*)malloc(sizeof(fft_data)* npad);
+	yfft = (fft_data*)malloc(sizeof(fft_data)* npad);
+	daughter = (fft_data*)malloc(sizeof(fft_data)* npad);
+	kwave = (double*)malloc(sizeof(double)* npad);
+
+	ymean = 0.0;
+
+	for (i = 0; i < N; ++i) {
+		ymean += y[i];
+	}
+
+	ymean /= N;
+
+	for (i = 0; i < N; ++i) {
+		ypad[i].re = y[i] - ymean;
+		ypad[i].im = 0.0;
+	}
+
+	for (i = N; i < npad; ++i) {
+		ypad[i].re = ypad[i].im = 0.0;
+	}
+
+
+	// Find FFT of the input y (ypad)
+
+	fft_exec(obj, ypad, yfft);
+
+	for (i = 0; i < npad; ++i) {
+		yfft[i].re /= (double) npad;
+		yfft[i].im /= (double) npad;
+	}
+
+
+	//Construct the wavenumber array
+
+	freq1 = 2.0*pi / ((double)npad*dt);
+	kwave[0] = 0.0;
+
+	for (i = 1; i < npad / 2 + 1; ++i) {
+		kwave[i] = i * freq1;
+	}
+
+	for (i = npad / 2 + 1; i < npad; ++i) {
+		kwave[i] = -kwave[npad - i ];
+	}
+
+	
+	// Main loop
+
+	for (j = 1; j <= jtot; ++j) {
+		scale1 = scale[j - 1];// = s0*pow(2.0, (double)(j - 1)*dj);
+		wave_function(npad, dt, mother, param, scale1, kwave, pi,&period1,&coi1, daughter);
+		period[j - 1] = period1;
+		for (k = 0; k < npad; ++k) {
+			tmp1 = daughter[k].re * yfft[k].re - daughter[k].im * yfft[k].im;
+			tmp2 = daughter[k].re * yfft[k].im + daughter[k].im * yfft[k].re;
+			daughter[k].re = tmp1;
+			daughter[k].im = tmp2;
+		}
+		fft_exec(iobj, daughter, ypad);
+		iter = 2 * (j - 1) * N;
+		for (i = 0; i < N; ++i) {
+			wave[iter + 2 * i] = ypad[i].re;
+			wave[iter + 2 * i + 1] = ypad[i].im;
+		}
+	}
+	
+
+	for (i = 1; i <= (N + 1) / 2; ++i) {
+		coi[i - 1] = coi1 * dt * ((double)i - 1.0);
+		coi[N - i] = coi[i - 1];
+	}
+	
+
+	
+	free(kwave);
+	free(ypad);
+	free(yfft);
+	free(daughter);
+
+	free_fft(obj);
+	free_fft(iobj);
+
+}
+
+void psi0(int mother, double param,double *val,int *real) {
+	double pi,coeff;
+	int m,sign;
+
+	m = (int)param;
+	pi = 4.0 * atan(1.0);
+
+	if (mother == 0) {
+		// Morlet
+		*val = 1.0 / sqrt(sqrt(pi));
+		*real = 1;
+	}
+	else if (mother == 1) {
+		//Paul
+		if (m % 2 == 0) {
+			*real = 1;
+		}
+		else {
+			*real = 0;
+		}
+
+		if (m % 4 == 0 || m % 4 == 1) {
+			sign = 1;
+		}
+		else {
+			sign = -1;
+		}
+		*val = sign * pow(2.0, (double)m) * factorial(m) / (sqrt(pi * factorial(2 * m)));
+
+	}
+	else if (mother == 2) {
+		// D.O.G
+		*real = 1;
+
+		if (m % 2 == 0) {
+			if (m % 4 == 0) {
+				sign = -1;
+			}
+			else {
+				sign = 1;
+			}
+			coeff = sign * pow(2.0, (double)m / 2) / gamma(0.5);
+			*val = coeff * gamma(((double)m + 1.0) / 2.0) / sqrt(gamma(m + 0.50));
+		}
+		else {
+			*val = 0;
+		}
+	}
+}
+
+static int maxabs(double *array,int N) {
+	double maxval,temp;
+	int i,index;
+	maxval = 0.0;
+	index = -1;
+
+	for (i = 0; i < N; ++i) {
+		temp = fabs(array[i]);
+		if (temp >= maxval) {
+			maxval = temp;
+			index = i;
+		}
+	}
+
+	return index;
+}
+
+
+double cdelta(int mother, double param, double psi0 ) {
+	int N,i,j,iter;
+	double *delta, *scale,*period,*wave,*coi,*mval;
+	double den,cdel;
+	double subscale,dt,dj,s0;
+	int jtot;
+	int maxarr;
+
+	subscale = 8.0;
+	dt = 0.25;
+	if (mother == 0) {
+		N = 16;
+		s0 = dt/4;
+	}
+	else if (mother == 1) {
+		N = 16;
+		s0 = dt / 4.0;
+	}
+	else if (mother == 2)
+	{
+		s0 = dt/8.0;
+		N = 256;
+		if (param == 2.0) {
+			subscale = 16.0;
+			s0 = dt / 16.0;
+			N = 2048;
+		}
+	}
+
+	dj = 1.0 / subscale;
+	jtot = 16 * (int) subscale;
+
+	delta = (double*)malloc(sizeof(double)* N);
+	wave = (double*)malloc(sizeof(double)* 2 * N * jtot);
+	coi = (double*)malloc(sizeof(double)* N);
+	scale = (double*)malloc(sizeof(double)* jtot);
+	period = (double*)malloc(sizeof(double)* jtot);
+	mval = (double*)malloc(sizeof(double)* N);
+
+
+	delta[0] = 1;
+
+	for (i = 1; i < N; ++i) {
+		delta[i] = 0;
+	}
+
+	for (i = 0; i < jtot; ++i) {
+		scale[i] = s0*pow(2.0, (double)(i)*dj);
+	}
+
+	cwavelet(delta, N, dt, mother, param, s0, dj, jtot, N, wave, scale, period, coi);
+
+	for (i = 0; i < N; ++i) {
+		mval[i] = 0;
+	}
+
+	for (j = 0; j < jtot; ++j) {
+		iter = 2 * j * N;
+		den = sqrt(scale[j]);
+		for (i = 0; i < N; ++i) {
+			mval[i] += wave[iter + 2 * i]/den;
+		}
+	}
+
+
+	maxarr = maxabs(mval, N);
+
+	cdel = sqrt(dt) * dj * mval[maxarr] / psi0;
+
+	free(delta);
+	free(wave);
+
+	free(scale);
+	free(period);
+	free(coi);
+	free(mval);
+
+	return cdel;
+}
+
+void icwavelet(double *wave, int N, double *scale,int jtot,double dt,double dj,double cdelta,double psi0,double *oup) {
+	int i, j,iter;
+	double den, coeff;
+
+	coeff = sqrt(dt) * dj / (cdelta *psi0);
+
+	for (i = 0; i < N; ++i) {
+		oup[i] = 0.0;
+	}
+
+	for (j = 0; j < jtot; ++j) {
+		iter = 2 * j * N;
+		den = sqrt(scale[j]);
+		for (i = 0; i < N; ++i) {
+			oup[i] += wave[iter + 2 * i] / den;
+		}
+	}
+
+	for (i = 0; i < N; ++i) {
+		oup[i] *= coeff;
+	}
+}
diff --git a/src/cwt.h b/src/cwt.h
new file mode 100755
index 0000000..ed9ea9a
--- /dev/null
+++ b/src/cwt.h
@@ -0,0 +1,29 @@
+#ifndef CWT_H_
+#define CWT_H_
+
+#include "wavefunc.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void cwavelet(double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad,
+		double *wave, double *scale, double *period, double *coi);
+
+void psi0(int mother, double param, double *val, int *real);
+
+double factorial(int N);
+
+double cdelta(int mother, double param, double psi0);
+
+void icwavelet(double *wave, int N, double *scale, int jtot, double dt, double dj, double cdelta, double psi0, double *oup);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */
+
+
diff --git a/src/cwtmath.c b/src/cwtmath.c
new file mode 100755
index 0000000..490c95f
--- /dev/null
+++ b/src/cwtmath.c
@@ -0,0 +1,296 @@
+#include "cwtmath.h"
+
+static void nsfft_fd(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *w) {
+	int M,N,i,j,L;
+	double delta,den,theta,tempr,tempi,plb;
+	double *temp1,*temp2;
+
+	N = obj->N;
+	L = N/2;
+	//w = (double*)malloc(sizeof(double)*N);
+	
+	M = divideby(N, 2);
+	
+	if (M == 0) {
+		printf("The Non-Standard FFT Length must be a power of 2");
+		exit(1);
+	}
+	
+	temp1 = (double*)malloc(sizeof(double)*L);
+	temp2 = (double*)malloc(sizeof(double)*L);
+	
+	delta = (ub - lb)/ N;
+	j = -N;
+	den = 2 * (ub-lb);
+	
+	for(i = 0; i < N;++i) {
+		w[i] = (double)j/den;
+		j += 2;
+	}
+	
+	fft_exec(obj,inp,oup);
+	
+
+	for (i = 0; i < L; ++i) {
+		temp1[i] = oup[i].re;
+		temp2[i] = oup[i].im;
+	}
+
+	for (i = 0; i < N - L; ++i) {
+		oup[i].re = oup[i + L].re;
+		oup[i].im = oup[i + L].im;
+	}
+
+	for (i = 0; i < L; ++i) {
+		oup[N - L + i].re = temp1[i];
+		oup[N - L + i].im = temp2[i];
+	}
+	
+	plb = PI2 * lb;
+	
+	for(i = 0; i < N;++i) {
+		tempr = oup[i].re;
+		tempi = oup[i].im;
+		theta = w[i] * plb;
+		
+		oup[i].re = delta * (tempr*cos(theta) + tempi*sin(theta));
+		oup[i].im = delta * (tempi*cos(theta) - tempr*sin(theta));
+	}
+
+	
+	//free(w);
+	free(temp1);
+	free(temp2);
+}
+
+static void nsfft_bk(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *t) {
+	int M,N,i,j,L;
+	double *w;
+	double delta,den,plb,theta;
+	double *temp1,*temp2;
+	fft_data *inpt;
+
+	N = obj->N;
+	L = N/2;
+	
+	M = divideby(N, 2);
+	
+	if (M == 0) {
+		printf("The Non-Standard FFT Length must be a power of 2");
+		exit(1);
+	}
+	
+	temp1 = (double*)malloc(sizeof(double)*L);
+	temp2 = (double*)malloc(sizeof(double)*L);
+	w = (double*)malloc(sizeof(double)*N);
+	inpt = (fft_data*) malloc (sizeof(fft_data) * N);
+	
+	delta = (ub - lb)/ N;
+	j = -N;
+	den = 2 * (ub-lb);
+	
+	for(i = 0; i < N;++i) {
+		w[i] = (double)j/den;
+		j += 2;
+	}
+	
+	plb = PI2 * lb;
+	
+	for(i = 0; i < N;++i) {
+		theta = w[i] * plb;
+		
+		inpt[i].re = (inp[i].re*cos(theta) - inp[i].im*sin(theta))/delta;
+		inpt[i].im = (inp[i].im*cos(theta) + inp[i].re*sin(theta))/delta;
+	}
+	
+	for (i = 0; i < L; ++i) {
+		temp1[i] = inpt[i].re;
+		temp2[i] = inpt[i].im;
+	}
+
+	for (i = 0; i < N - L; ++i) {
+		inpt[i].re = inpt[i + L].re;
+		inpt[i].im = inpt[i + L].im;
+	}
+
+	for (i = 0; i < L; ++i) {
+		inpt[N - L + i].re = temp1[i];
+		inpt[N - L + i].im = temp2[i];
+	}
+	
+	fft_exec(obj,inpt,oup);
+	
+	for(i = 0; i < N;++i) {
+		t[i] = lb + i*delta;
+	}
+	
+	free(w);
+	free(temp1);
+	free(temp2);
+	free(inpt);
+}
+
+void nsfft_exec(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *w) {
+	if (obj->sgn == 1) {
+		nsfft_fd(obj,inp,oup,lb,ub,w);
+	} else if (obj->sgn == -1) {
+		nsfft_bk(obj,inp,oup,lb,ub,w);
+	}
+}
+
+static double fix(double x) {
+	// Rounds to the integer nearest to zero 
+	if (x >= 0.) {
+		return floor(x);
+	} else {
+		return ceil(x);
+	}
+}
+
+int nint(double N) {
+	int i;
+
+	i = (int)(N + 0.49999);
+
+	return i;
+}
+
+double gamma(double x) {
+	/*
+	 * This C program code is based on  W J Cody's fortran code.
+	 * http://www.netlib.org/specfun/gamma
+	 * 
+	 * References:
+   "An Overview of Software Development for Special Functions",
+	W. J. Cody, Lecture Notes in Mathematics, 506,
+	Numerical Analysis Dundee, 1975, G. A. Watson (ed.),
+	Springer Verlag, Berlin, 1976.
+
+   Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
+   */ 
+   
+	// numerator and denominator coefficients for 1 <= x <= 2
+	
+	double y,oup,fact,sum,y2,yi,z,nsum,dsum;
+	int swi,n,i;
+	
+	double spi = 0.9189385332046727417803297;
+    double pi  = 3.1415926535897932384626434;
+	double xmax = 171.624e+0;
+	double xinf = 1.79e308;
+	double eps = 2.22e-16;
+	double xninf = 1.79e-308;
+	
+	double num[8] = { -1.71618513886549492533811e+0,
+                        2.47656508055759199108314e+1,
+                       -3.79804256470945635097577e+2,
+                        6.29331155312818442661052e+2,
+                        8.66966202790413211295064e+2,
+                       -3.14512729688483675254357e+4,
+                       -3.61444134186911729807069e+4,
+                        6.64561438202405440627855e+4 };
+                       
+    double den[8] = { -3.08402300119738975254353e+1,
+                        3.15350626979604161529144e+2,
+                       -1.01515636749021914166146e+3,
+                       -3.10777167157231109440444e+3,
+                        2.25381184209801510330112e+4,
+                        4.75584627752788110767815e+3,
+                       -1.34659959864969306392456e+5,
+                       -1.15132259675553483497211e+5 }; 
+                       
+    // Coefficients for Hart's Minimax approximation x >= 12       
+    
+    
+	double c[7] = { -1.910444077728e-03,
+                        8.4171387781295e-04,
+                       -5.952379913043012e-04,
+                        7.93650793500350248e-04,
+                       -2.777777777777681622553e-03,
+                        8.333333333333333331554247e-02,
+                        5.7083835261e-03 };            
+    
+    y = x;             
+    swi = 0; 
+    fact = 1.0;
+    n = 0;      
+    
+    
+    if ( y < 0.) {
+		// Negative x
+		y = -x;
+		yi = fix(y);
+		oup = y - yi;
+		
+		if (oup != 0.0) {
+			if (yi != fix(yi * .5) * 2.) {
+				swi = 1;
+			}
+			fact = -pi / sin(pi * oup);
+			y += 1.;
+		} else {
+			return xinf;
+		}
+	}      
+	
+	if (y < eps) {
+		if (y >= xninf) {
+			oup = 1.0/y;
+		} else {
+			return xinf;
+		}
+		
+	} else if (y < 12.) {
+		yi = y;
+		if ( y < 1.) {
+			z = y;
+			y += 1.;
+		} else {
+			n = ( int ) y - 1;
+            y -= ( double ) n;
+            z = y - 1.0;
+		}
+		nsum = 0.;
+		dsum = 1.;
+		for (i = 0; i < 8; ++i) {
+			nsum = (nsum + num[i]) * z;
+			dsum = dsum * z + den[i];
+		}
+		oup = nsum / dsum + 1.;
+		
+		if (yi < y) {
+	   
+			oup /= yi;
+		} else if (yi > y) {
+	    
+			for (i = 0; i < n; ++i) {
+				oup *= y;
+				y += 1.;
+			}
+		
+		}      
+	
+	} else {
+		if (y <= xmax) {
+			y2 = y * y;
+			sum = c[6];
+			for (i = 0; i < 6; ++i) {
+				sum = sum / y2 + c[i];
+			}
+			sum = sum / y - y + spi;
+			sum += (y - .5) * log(y);
+			oup = exp(sum);
+		} else {
+			return(xinf);
+		}
+	}
+	
+	if (swi) {
+		oup = -oup;
+	}
+    if (fact != 1.) {
+		oup = fact / oup; 
+	}       
+	
+	return oup;
+}
diff --git a/src/cwtmath.h b/src/cwtmath.h
new file mode 100755
index 0000000..f989394
--- /dev/null
+++ b/src/cwtmath.h
@@ -0,0 +1,22 @@
+#ifndef CWTMATH_H_
+#define CWTMATH_H_
+
+#include "wtmath.h"
+#include "hsfft.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void nsfft_exec(fft_object obj, fft_data *inp, fft_data *oup,double lb,double ub,double *w);// lb -lower bound, ub - upper bound, w - time or frequency grid (Size N)
+
+double gamma(double x);
+
+int nint(double N);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */
diff --git a/src/wavefunc.c b/src/wavefunc.c
new file mode 100755
index 0000000..6fa2752
--- /dev/null
+++ b/src/wavefunc.c
@@ -0,0 +1,210 @@
+#include "wavefunc.h"
+
+void meyer(int N,double lb,double ub,double *phi,double *psi,double *tgrid) {
+	int M,i;
+	double *w;
+	double delta,j;
+	double theta,x,x2,x3,x4,v,cs,sn;
+	double wf;
+	fft_data *phiw,*psiw,*oup;
+	fft_object obj;
+	
+	M = divideby(N, 2);
+	
+	if (M == 0) {
+		printf("Size of Wavelet must be a power of 2");
+		exit(1);
+	}
+	if (lb >= ub) {
+		printf("upper bound must be greater than lower bound");
+		exit(1);
+	}
+	
+	obj = fft_init(N,-1);
+	w = (double*)malloc(sizeof(double)*N);
+	phiw = (fft_data*) malloc (sizeof(fft_data) * N);
+	psiw = (fft_data*) malloc (sizeof(fft_data) * N);
+	oup = (fft_data*) malloc (sizeof(fft_data) * N);
+	
+	delta = 2 * (ub-lb) / PI2;
+	
+	j = (double) N;
+	j *= -1.0;
+	
+	for(i = 0; i < N;++i) {
+		w[i] = j / delta;
+		j += 2.0; 
+		psiw[i].re = psiw[i].im = 0.0;
+		phiw[i].re = phiw[i].im = 0.0;
+	}
+	
+	
+	for(i = 0; i < N;++i) {
+		wf = fabs(w[i]);
+		if (wf <= PI2/3.0) {
+			phiw[i].re = 1.0;
+		}
+		if (wf > PI2/3.0 && wf <= 2 * PI2 / 3.0) {
+			x = (3 * wf / PI2) - 1.0;
+			x2 = x*x;
+			x3 = x2 * x;
+			x4 = x3 *x;
+			v = x4 *(35 - 84*x + 70*x2 - 20*x3);
+			theta = v * PI2 / 4.0;
+			cs = cos(theta);
+			sn = sin(theta);
+			
+			phiw[i].re = cs;
+			psiw[i].re = cos(w[i]/2.0) * sn;
+			psiw[i].im = sin(w[i]/2.0) * sn;
+		}
+		if (wf > 2.0 * PI2/3.0 && wf <= 4 * PI2 / 3.0) {
+			x = (1.5 * wf / PI2) - 1.0;
+			x2 = x*x;
+			x3 = x2 * x;
+			x4 = x3 *x;
+			v = x4 *(35 - 84*x + 70*x2 - 20*x3);
+			theta = v * PI2 / 4.0;
+			cs = cos(theta);
+			
+			psiw[i].re = cos(w[i]/2.0) * cs;
+			psiw[i].im = sin(w[i]/2.0) * cs;
+		}
+	}
+	
+	
+	nsfft_exec(obj,phiw,oup,lb,ub,tgrid);
+	
+	
+	for(i = 0; i < N;++i) {
+		phi[i] = oup[i].re/N;
+	}
+	
+	nsfft_exec(obj,psiw,oup,lb,ub,tgrid);
+	
+	
+	for(i = 0; i < N;++i) {
+		psi[i] = oup[i].re/N;
+	}
+	
+	
+	
+	free(oup);
+	free(phiw);
+	free(psiw);
+	free(w);
+}
+
+void gauss(int N,int p,double lb,double ub,double *psi,double *t) {
+  double delta,num,den,t2,t4;
+	int i;
+	
+	if (lb >= ub) {
+		printf("upper bound must be greater than lower bound");
+		exit(1);
+	}
+	
+	t[0] = lb;
+	t[N-1] = ub;
+	delta = (ub - lb) / (N-1);
+	for(i = 1; i < N-1;++i) {
+		t[i] = lb + delta * i;
+	}
+	
+	den = sqrt(gamma(p+0.5));
+	
+	if ((p+1)%2 == 0) {
+		num = 1.0;
+	} else {
+		num = -1.0;
+	}
+	
+	num /= den;
+	
+	//printf("\n%g\n",num);
+	
+	if (p == 1) {
+		for(i = 0; i < N;++i) {
+			psi[i] = -t[i] * exp(- t[i] * t[i]/2.0) * num; 
+		}
+	} else if (p == 2) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = (-1.0 + t2) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 3) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = t[i] * (3.0 - t2) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 4) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = (t2 * t2 - 6.0 * t2 + 3.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 5) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = t[i] * (-t2 * t2 + 10.0 * t2 - 15.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 6) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = (t2 * t2 * t2 - 15.0 * t2 * t2 + 45.0 * t2 - 15.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 7) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			psi[i] = t[i] * (-t2 * t2 * t2 + 21.0 * t2 * t2 - 105.0 * t2 + 105.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 8) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			t4 = t2 * t2;
+			psi[i] = (t4 * t4 - 28.0 * t4 * t2 + 210.0 * t4 - 420.0 * t2 + 105.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 9) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			t4 = t2 * t2;
+			psi[i] = t[i] * (- t4 * t4 + 36.0 * t4 * t2 - 378.0 * t4 + 1260.0 * t2 - 945.0) * exp(- t2/2.0) * num; 
+		}
+	} else if (p == 10) {
+		for(i = 0; i < N;++i) {
+			t2 = t[i] * t[i];
+			t4 = t2 * t2;
+			psi[i] = (t4 * t4 * t2 - 45.0 * t4 * t4 + 630.0 * t4 * t2 - 3150.0 * t4 + 4725.0 * t2 - 945.0) * exp(- t2/2.0) * num; 
+		}
+	} else {
+	  printf("\n The Gaussian Derivative Wavelet is only available for Derivatives 1 to 10");
+	  exit(1);
+	}
+
+	
+	
+}
+
+void mexhat(int N,double lb,double ub,double *psi,double *t) {
+  gauss(N,2,lb,ub,psi,t);
+}
+
+void morlet(int N,double lb,double ub,double *psi,double *t) {
+       int i;
+       double delta;
+
+  	if (lb >= ub) {
+		printf("upper bound must be greater than lower bound");
+		exit(1);
+	}
+	
+	t[0] = lb;
+	t[N-1] = ub;
+	delta = (ub - lb) / (N-1);
+	for(i = 1; i < N-1;++i) {
+		t[i] = lb + delta * i;
+	}
+
+	for(i = 0; i < N;++i) {
+	  psi[i] = exp(- t[i] * t[i] / 2.0) * cos(5 * t[i]);
+	}
+}
diff --git a/src/wavefunc.h b/src/wavefunc.h
new file mode 100755
index 0000000..2aa6ead
--- /dev/null
+++ b/src/wavefunc.h
@@ -0,0 +1,24 @@
+#ifndef WAVEFUNC_H_
+#define WAVEFUNC_H_
+
+#include "cwtmath.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void meyer(int N,double lb,double ub,double *phi,double *psi,double *tgrid);
+
+void gauss(int N,int p,double lb,double ub,double *psi,double *t);
+
+void mexhat(int N,double lb,double ub,double *psi,double *t);
+
+void morlet(int N,double lb,double ub,double *psi,double *t);  
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVEFUNC_H_ */
diff --git a/src/wavelib.c b/src/wavelib.c
index c2d7666..0ce0023 100644
--- a/src/wavelib.c
+++ b/src/wavelib.c
@@ -1,1594 +1,1797 @@
-/*
-  Copyright (c) 2014, Rafat Hussain
-*/
-
-#include "wavelib.h"
-
-wave_object wave_init(char* wname) {
-	wave_object obj = NULL;
-	int retval;
-	retval = 0;
-
-	if (wname != NULL) {
-		retval = filtlength(wname);
-		//obj->filtlength = retval;
-		//strcopy(obj->wname, wname);
-	}
-
-    obj = (wave_object)malloc(sizeof(struct wave_set) + sizeof(double) * 4 * retval);
-
-	obj->filtlength = retval;
-	obj->lpd_len = obj->hpd_len = obj->lpr_len = obj->hpr_len = obj->filtlength;
-	strcpy(obj->wname, wname);
-	if (wname != NULL) {
-		filtcoef(wname,obj->params,obj->params+retval,obj->params+2*retval,obj->params+3*retval);
-	}
-	obj->lpd = &obj->params[0];
-	obj->hpd = &obj->params[retval];
-	obj->lpr = &obj->params[2 * retval];
-	obj->hpr = &obj->params[3 * retval];
-	return obj;
-}
-
-wt_object wt_init(wave_object wave,char* method, int siglength,int J) {
-	int size,i,MaxIter;
-	wt_object obj = NULL;
-
-	size = wave->filtlength;
-
-	if (J > 100) {
-		printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
-		exit(-1);
-	}
-
-	MaxIter = wmaxiter(siglength, size);
-
-	if (J > MaxIter) {
-		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter);
-		exit(-1);
-	}
-
-	if (method == NULL) {
-        obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength + 2 * J * (size + 1)));
-		obj->outlength = siglength + 2 * J * (size + 1); // Default
-		strcpy(obj->ext, "sym"); // Default
-	}
-	else if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
-        obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength + 2 * J * (size + 1)));
-		obj->outlength = siglength + 2 * J * (size + 1); // Default
-		strcpy(obj->ext, "sym"); // Default
-	}
-	else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) {
-		if (!testSWTlength(siglength, J)) {
-			printf("\n For SWT the signal length must be a multiple of 2^J. \n");
-			exit(-1);
-		}
-
-        obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1)));
-		obj->outlength = siglength * (J + 1); // Default
-		strcpy(obj->ext, "per"); // Default
-	}
-	else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
-
-		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 = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1)));
-		obj->outlength = siglength * (J + 1); // Default
-		strcpy(obj->ext, "per"); // Default
-	}
-
-	obj->wave = wave;
-	obj->siglength = siglength;
-	obj->J = J;
-	obj->MaxIter = MaxIter;
-	strcpy(obj->method, method);
-
-	if (siglength % 2 == 0) {
-		obj->even = 1;
-	}
-	else {
-		obj->even = 0;
-	}
-
-	obj->cobj = NULL;
-
-	strcpy(obj->cmethod, "direct"); // Default
-	obj->cfftset = 0;
-	obj->lenlength = J + 2;
-	obj->output = &obj->params[0];
-	if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
-		for (i = 0; i < siglength + 2 * J * (size + 1); ++i) {
-			obj->params[i] = 0.0;
-		}
-	}
-	else if (!strcmp(method, "swt") || !strcmp(method, "SWT") || !strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
-		for (i = 0; i < siglength * (J + 1); ++i) {
-			obj->params[i] = 0.0;
-		}
-	}
-	//wave_summary(obj->wave);
-
-	return obj;
-}
-
-wtree_object wtree_init(wave_object wave, int siglength,int J) {
-    int size,i,MaxIter,temp,temp2,elength,nodes;
-	wtree_object obj = NULL;
-
-	size = wave->filtlength;
-
-	if (J > 100) {
-		printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
-		exit(-1);
-	}
-
-
-	MaxIter = wmaxiter(siglength, size);
-	if (J > MaxIter) {
-		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
-		exit(-1);
-	}
-	temp = 1;
-	elength = 0;
-	nodes = 0;
-	for(i = 0; i < J;++i) {
-	  temp *= 2;
-	  nodes += temp;
-	  temp2 = (size - 2) * (temp - 1);
-	  elength += temp2;
-	}
-
-    obj = (wtree_object)malloc(sizeof(struct wtree_set) + sizeof(double)* (siglength * (J + 1) + elength + nodes + J + 1));
-	obj->outlength = siglength * (J + 1) + elength;
-	strcpy(obj->ext, "sym");
-
-	obj->wave = wave;
-	obj->siglength = siglength;
-	obj->J = J;
-	obj->MaxIter = MaxIter;
-	strcpy(obj->method, "dwt");
-
-	if (siglength % 2 == 0) {
-		obj->even = 1;
-	}
-	else {
-		obj->even = 0;
-	}
-
-	obj->cobj = NULL;
-	obj->nodes = nodes;
-
-	obj->cfftset = 0;
-	obj->lenlength = J + 2;
-	obj->output = &obj->params[0];
-	obj->nodelength = (int*) &obj->params[siglength * (J + 1) + elength];
-	obj->coeflength = (int*)&obj->params[siglength * (J + 1) + elength + nodes];
-
-	for (i = 0; i < siglength * (J + 1) + elength + nodes + J + 1; ++i) {
-	       obj->params[i] = 0.0;
-	}
-
-	//wave_summary(obj->wave);
-
-	return obj;
-}
-
-wpt_object wpt_init(wave_object wave, int siglength, int J) {
-	int size, i, MaxIter, temp, nodes,elength,p2,N,lp;
-	wpt_object obj = NULL;
-
-	size = wave->filtlength;
-
-	if (J > 100) {
-		printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
-		exit(-1);
-	}
-
-
-	MaxIter = wmaxiter(siglength, size);
-	if (J > MaxIter) {
-		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
-		exit(-1);
-	}
-	temp = 1;
-	nodes = 0;
-	for (i = 0; i < J; ++i) {
-		temp *= 2;
-		nodes += temp;
-	}
-
-	i = J;
-	p2 = 2;
-	N = siglength;
-	lp = size;
-	elength = 0;
-	while (i > 0) {
-		N = N + lp - 2;
-		N = (int)ceil((double)N / 2.0);
-		elength = p2 * N;
-		i--;
-		p2 *= 2;
-	}
-	//printf("elength %d", elength);
-
-    obj = (wpt_object)malloc(sizeof(struct wpt_set) + sizeof(double)* (elength + 4 * nodes + 2 * J + 6));
-	obj->outlength = siglength + 2 * (J + 1) * (size + 1);
-	strcpy(obj->ext, "sym");
-	strcpy(obj->entropy, "shannon");
-	obj->eparam = 0.0;
-
-	obj->wave = wave;
-	obj->siglength = siglength;
-	obj->J = J;
-	obj->MaxIter = MaxIter;
-
-	if (siglength % 2 == 0) {
-		obj->even = 1;
-	}
-	else {
-		obj->even = 0;
-	}
-
-	obj->cobj = NULL;
-	obj->nodes = nodes;
-
-	obj->lenlength = J + 2;
-	obj->output = &obj->params[0];
-	obj->costvalues = &obj->params[elength];
-	obj->basisvector = &obj->params[elength + nodes + 1];
-	obj->nodeindex = (int*)&obj->params[elength + 2*nodes + 2];
-	obj->numnodeslevel = (int*)&obj->params[elength + 4 * nodes + 4];
-	obj->coeflength = (int*)&obj->params[elength + 4 * nodes + J + 5];
-
-	for (i = 0; i < elength + 4 * nodes + 2 * J + 6; ++i) {
-		obj->params[i] = 0.0;
-	}
-
-	//wave_summary(obj->wave);
-
-	return obj;
-}
-
-
-static void wconv(wt_object wt, double *sig, int N, double *filt, int L, double *oup) {
-	if (!strcmp(wt->cmethod,"direct")) {
-		conv_direct(sig, N, filt, L, oup);
-	}
-	else if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) {
-		if (wt->cfftset == 0) {
-			wt->cobj = conv_init(N, L);
-			conv_fft(wt->cobj, sig, filt, oup);
-			free_conv(wt->cobj);
-		}
-		else {
-			conv_fft(wt->cobj, sig, filt, oup);
-		}
-	}
-	else {
-		printf("Convolution Only accepts two methods - direct and fft");
-		exit(-1);
-	}
-}
-
-
-static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int l,l2,isodd,i,t,len_avg;
-
-	len_avg = wt->wave->lpd_len;
-	l2 = len_avg / 2;
-	isodd = N % 2;
-
-	for (i = 0; i < len_cA; ++i) {
-		t = 2 * i + l2;
-		cA[i] = 0.0;
-		cD[i] = 0.0;
-		for (l = 0; l < len_avg; ++l) {
-			if ((t - l) >= l2 && (t - l) < N) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < l2 && (t - l) >= 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < 0 && isodd == 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l + N];
-				cD[i] += wt->wave->hpd[l] * inp[t - l + N];
-			}
-			else if ((t - l) < 0 && isodd == 1) {
-				if ((t - l) != -1) {
-					cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
-					cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
-				}
-				else {
-					cA[i] += wt->wave->lpd[l] * inp[N - 1];
-					cD[i] += wt->wave->hpd[l] * inp[N - 1];
-				}
-			}
-			else if ((t - l) >= N && isodd == 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l - N];
-				cD[i] += wt->wave->hpd[l] * inp[t - l - N];
-			}
-			else if ((t - l) >= N && isodd == 1) {
-				if (t - l != N) {
-					cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
-					cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
-				}
-				else {
-					cA[i] += wt->wave->lpd[l] * inp[N - 1];
-					cD[i] += wt->wave->hpd[l] * inp[N - 1];
-				}
-			}
-
-		}
-	}
-
-
-
-}
-
-static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int l, l2, isodd, i, t, len_avg;
-
-	len_avg = wt->wave->lpd_len;
-	l2 = len_avg / 2;
-	isodd = N % 2;
-
-	for (i = 0; i < len_cA; ++i) {
-		t = 2 * i + l2;
-		cA[i] = 0.0;
-		cD[i] = 0.0;
-		for (l = 0; l < len_avg; ++l) {
-			if ((t - l) >= l2 && (t - l) < N) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < l2 && (t - l) >= 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < 0 && isodd == 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l + N];
-				cD[i] += wt->wave->hpd[l] * inp[t - l + N];
-			}
-			else if ((t - l) < 0 && isodd == 1) {
-				if ((t - l) != -1) {
-					cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
-					cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
-				}
-				else {
-					cA[i] += wt->wave->lpd[l] * inp[N - 1];
-					cD[i] += wt->wave->hpd[l] * inp[N - 1];
-				}
-			}
-			else if ((t - l) >= N && isodd == 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l - N];
-				cD[i] += wt->wave->hpd[l] * inp[t - l - N];
-			}
-			else if ((t - l) >= N && isodd == 1) {
-				if (t - l != N) {
-					cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
-					cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
-				}
-				else {
-					cA[i] += wt->wave->lpd[l] * inp[N - 1];
-					cD[i] += wt->wave->hpd[l] * inp[N - 1];
-				}
-			}
-
-		}
-	}
-
-
-
-}
-
-static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int l, l2, isodd, i, t, len_avg;
-
-	len_avg = wt->wave->lpd_len;
-	l2 = len_avg / 2;
-	isodd = N % 2;
-
-	for (i = 0; i < len_cA; ++i) {
-		t = 2 * i + l2;
-		cA[i] = 0.0;
-		cD[i] = 0.0;
-		for (l = 0; l < len_avg; ++l) {
-			if ((t - l) >= l2 && (t - l) < N) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < l2 && (t - l) >= 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < 0 && isodd == 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l + N];
-				cD[i] += wt->wave->hpd[l] * inp[t - l + N];
-			}
-			else if ((t - l) < 0 && isodd == 1) {
-				if ((t - l) != -1) {
-					cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
-					cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
-				}
-				else {
-					cA[i] += wt->wave->lpd[l] * inp[N - 1];
-					cD[i] += wt->wave->hpd[l] * inp[N - 1];
-				}
-			}
-			else if ((t - l) >= N && isodd == 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l - N];
-				cD[i] += wt->wave->hpd[l] * inp[t - l - N];
-			}
-			else if ((t - l) >= N && isodd == 1) {
-				if (t - l != N) {
-					cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
-					cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
-				}
-				else {
-					cA[i] += wt->wave->lpd[l] * inp[N - 1];
-					cD[i] += wt->wave->hpd[l] * inp[N - 1];
-				}
-			}
-
-		}
-	}
-
-
-
-}
-
-static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int i,l, t, len_avg;
-
-	len_avg = wt->wave->lpd_len;
-
-	for (i = 0; i < len_cA; ++i) {
-		t = 2 * i + 1;
-		cA[i] = 0.0;
-		cD[i] = 0.0;
-		for (l = 0; l < len_avg; ++l) {
-			if ((t - l) >= 0 && (t - l) < N) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < 0) {
-				cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
-				cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
-			}
-			else if ((t - l) >= N) {
-				cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
-				cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
-			}
-
-		}
-	}
-
-
-}
-
-static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int i, l, t, len_avg;
-
-	len_avg = wt->wave->lpd_len;
-
-	for (i = 0; i < len_cA; ++i) {
-		t = 2 * i + 1;
-		cA[i] = 0.0;
-		cD[i] = 0.0;
-		for (l = 0; l < len_avg; ++l) {
-			if ((t - l) >= 0 && (t - l) < N) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < 0) {
-				cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
-				cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
-			}
-			else if ((t - l) >= N) {
-				cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
-				cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
-			}
-
-		}
-	}
-
-
-}
-
-static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int i, l, t, len_avg;
-
-	len_avg = wt->wave->lpd_len;
-
-	for (i = 0; i < len_cA; ++i) {
-		t = 2 * i + 1;
-		cA[i] = 0.0;
-		cD[i] = 0.0;
-		for (l = 0; l < len_avg; ++l) {
-			if ((t - l) >= 0 && (t - l) < N) {
-				cA[i] += wt->wave->lpd[l] * inp[t - l];
-				cD[i] += wt->wave->hpd[l] * inp[t - l];
-			}
-			else if ((t - l) < 0) {
-				cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
-				cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
-			}
-			else if ((t - l) >= N) {
-				cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
-				cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
-			}
-
-		}
-	}
-
-
-}
-
-static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, int len_cA, double *cD, int len_cD) {
-	int len_avg,D,lf;
-	double *signal,*cA_undec;
-	len_avg = (wt->wave->lpd_len + wt->wave->hpd_len) / 2;
-	//len_sig = 2 * (int)ceil((double)len_sig / 2.0);
-
-	D = 2;
-
-	if (!strcmp(wt->ext, "per")) {
-		signal = (double*)malloc(sizeof(double)* (len_sig + len_avg + (len_sig % 2)));
-
-		len_sig = per_ext(sig, len_sig, len_avg / 2, signal);
-
-		cA_undec = (double*)malloc(sizeof(double)* (len_sig + len_avg + wt->wave->lpd_len - 1));
-
-		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-			wt->cobj = conv_init(len_sig + len_avg, wt->wave->lpd_len);
-			wt->cfftset = 1;
-		}
-		else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
-			printf("Decomposition Filters must have the same length.");
-			exit(-1);
-		}
-
-		wconv(wt, signal, len_sig + len_avg, wt->wave->lpd, wt->wave->lpd_len, cA_undec);
-
-		downsamp(cA_undec + len_avg, len_sig, D, cA);
-
-		wconv(wt, signal, len_sig + len_avg, wt->wave->hpd, wt->wave->hpd_len, cA_undec);
-
-		downsamp(cA_undec + len_avg, len_sig, D, cD);
-	}
-	else if (!strcmp(wt->ext, "sym")) {
-		//printf("\n YES %s \n", wt->ext);
-		lf = wt->wave->lpd_len;// lpd and hpd have the same length
-
-		signal = (double*)malloc(sizeof(double)* (len_sig + 2 * (lf - 1)));
-
-		len_sig = symm_ext(sig, len_sig, lf - 1, signal);
-
-		cA_undec = (double*)malloc(sizeof(double)* (len_sig + 3 * (lf - 1)));
-
-		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-			wt->cobj = conv_init(len_sig + 2 * (lf - 1), lf);
-			wt->cfftset = 1;
-		}
-		else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
-			printf("Decomposition Filters must have the same length.");
-			exit(-1);
-		}
-
-
-		wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->lpd, wt->wave->lpd_len, cA_undec);
-
-		downsamp(cA_undec + lf, len_sig + lf - 2, D, cA);
-
-		wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->hpd, wt->wave->hpd_len, cA_undec);
-
-		downsamp(cA_undec + lf, len_sig + lf - 2, D, cD);
-	}
-	else {
-		printf("Signal extension can be either per or sym");
-		exit(-1);
-	}
-
-
-	if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-		free_conv(wt->cobj);
-		wt->cfftset = 0;
-	}
-
-	free(signal);
-	free(cA_undec);
-}
-
-void dwt(wt_object wt,double *inp) {
-	int i,J,temp_len,iter,N,lp;
-	int len_cA;
-	double *orig,*orig2;
-
-	temp_len = wt->siglength;
-	J = wt->J;
-	wt->length[J + 1] = temp_len;
-	wt->outlength = 0;
-	wt->zpad = 0;
-	orig = (double*)malloc(sizeof(double)* temp_len);
-	orig2 = (double*)malloc(sizeof(double)* temp_len);
-	/*
-	if ((temp_len % 2) == 0) {
-	wt->zpad = 0;
-	orig = (double*)malloc(sizeof(double)* temp_len);
-	orig2 = (double*)malloc(sizeof(double)* temp_len);
-	}
-	else {
-	wt->zpad = 1;
-	temp_len++;
-	orig = (double*)malloc(sizeof(double)* temp_len);
-	orig2 = (double*)malloc(sizeof(double)* temp_len);
-	}
-	*/
-
-	for (i = 0; i < wt->siglength; ++i) {
-		orig[i] = inp[i];
-	}
-
-	if (wt->zpad == 1) {
-		orig[temp_len - 1] = orig[temp_len - 2];
-	}
-
-	N = temp_len;
-	lp = wt->wave->lpd_len;
-
-	if (!strcmp(wt->ext,"per")) {
-		i = J;
-		while (i > 0) {
-			N = (int)ceil((double)N / 2.0);
-			wt->length[i] = N;
-			wt->outlength += wt->length[i];
-			i--;
-		}
-		wt->length[0] = wt->length[1];
-		wt->outlength += wt->length[0];
-		N = wt->outlength;
-
-		for (iter = 0; iter < J; ++iter) {
-			len_cA = wt->length[J - iter];
-			N -= len_cA;
-			if ( !strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT") ) {
-				dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
-			}
-			else {
-				dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
-			}
-			temp_len = wt->length[J - iter];
-			if (iter == J - 1) {
-				for (i = 0; i < len_cA; ++i) {
-					wt->params[i] = orig2[i];
-				}
-			}
-			else {
-				for (i = 0; i < len_cA; ++i) {
-					orig[i] = orig2[i];
-				}
-			}
-		}
-	}
-	else if (!strcmp(wt->ext,"sym")) {
-		//printf("\n YES %s \n", wt->ext);
-		i = J;
-		while (i > 0) {
-			N = N + lp - 2;
-			N = (int) ceil((double)N / 2.0);
-			wt->length[i] = N;
-			wt->outlength += wt->length[i];
-			i--;
-		}
-		wt->length[0] = wt->length[1];
-		wt->outlength += wt->length[0];
-		N = wt->outlength;
-
-		for (iter = 0; iter < J; ++iter) {
-			len_cA = wt->length[J - iter];
-			N -= len_cA;
-			if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) {
-				dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
-			}
-			else {
-				dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
-			}
-			temp_len = wt->length[J - iter];
-
-			if (iter == J - 1) {
-				for (i = 0; i < len_cA; ++i) {
-					wt->params[i] = orig2[i];
-				}
-			}
-			else {
-				for (i = 0; i < len_cA; ++i) {
-					orig[i] = orig2[i];
-				}
-			}
-		}
-	}
-	else {
-		printf("Signal extension can be either per or sym");
-		exit(-1);
-	}
-
-	free(orig);
-	free(orig2);
-}
-
-void wtree(wtree_object wt,double *inp) {
-	int i,J,temp_len,iter,N,lp,p2,k,N2,Np;
-	int len_cA,t,t2,it1;
-	double *orig;
-
-	temp_len = wt->siglength;
-	J = wt->J;
-	wt->length[J + 1] = temp_len;
-	wt->outlength = 0;
-	wt->zpad = 0;
-	orig = (double*)malloc(sizeof(double)* temp_len);
-	/*
-	if ((temp_len % 2) == 0) {
-		wt->zpad = 0;
-		orig = (double*)malloc(sizeof(double)* temp_len);
-	}
-	else {
-		wt->zpad = 1;
-		temp_len++;
-		orig = (double*)malloc(sizeof(double)* temp_len);
-	}
-	*/
-	for (i = 0; i < wt->siglength; ++i) {
-		orig[i] = inp[i];
-	}
-
-	if (wt->zpad == 1) {
-		orig[temp_len - 1] = orig[temp_len - 2];
-	}
-
-	N = temp_len;
-	lp = wt->wave->lpd_len;
-        p2 = 1;
-
-	if (!strcmp(wt->ext,"per")) {
-		i = J;
-                p2 = 2;
-		while (i > 0) {
-			N = (int)ceil((double)N / 2.0);
-			wt->length[i] = N;
-			wt->outlength += p2 * (wt->length[i]);
-			i--;
-                        p2 *= 2;
-		}
-		wt->length[0] = wt->length[1];
-
-		N2 = N = wt->outlength;
-                p2 = 1;
-		for (iter = 0; iter < J; ++iter) {
-			len_cA = wt->length[J - iter];
-                        N2 -= 2 * p2 * len_cA;
-                        N = N2;
-                        for(k = 0; k < p2;++k) {
-                            if (iter == 0) {
-                               wtree_per(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
-                            } else {
-                                wtree_per(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
-                            }
-                            N += 2 * len_cA;
-                        }
-
-			temp_len = wt->length[J - iter];
-			p2 = 2 * p2;
-            Np = N2;
-		}
-	}
-	else if (!strcmp(wt->ext,"sym")) {
-		//printf("\n YES %s \n", wt->ext);
-		i = J;
-                p2 = 2;
-		while (i > 0) {
-			N = N + lp - 2;
-			N = (int) ceil((double)N / 2.0);
-			wt->length[i] = N;
-			wt->outlength += p2 * (wt->length[i]);
-			i--;
-                        p2 *= 2;
-		}
-		wt->length[0] = wt->length[1];
-
-		N2 = N = wt->outlength;
-                p2 = 1;
-
-                for (iter = 0; iter < J; ++iter) {
-			len_cA = wt->length[J - iter];
-                        N2 -= 2 * p2 * len_cA;
-                        N = N2;
-                        for(k = 0; k < p2;++k) {
-                            if (iter == 0) {
-                                wtree_sym(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
-                            } else {
-                                wtree_sym(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
-                            }
-                            N += 2 * len_cA;
-                        }
-
-			temp_len = wt->length[J - iter];
-			p2 = 2 * p2;
-            Np = N2;
-		}
-
-	}
-	else {
-		printf("Signal extension can be either per or sym");
-		exit(-1);
-	}
-
-	J = wt->J;
-	t2 = wt->outlength - 2 * wt->length[J];
-	p2 = 2;
-	it1 = 0;
-	for (i = 0; i < J; ++i) {
-		t = t2;
-		for (k = 0; k < p2; ++k) {
-			wt->nodelength[it1] = t;
-			it1++;
-			t += wt->length[J - i];
-		}
-		p2 *= 2;
-		t2 = t2 - p2 * wt->length[J - i - 1];
-	}
-
-	wt->coeflength[0] = wt->siglength;
-
-	for (i = 1; i < J + 1; ++i) {
-		wt->coeflength[i] = wt->length[J - i + 1];
-	}
-
-	free(orig);
-}
-
-static int ipow2(int n) {
-	int p,i;
-	p = 1;
-
-	for (i = 0; i < n; ++i) {
-		p *= 2;
-	}
-
-	return p;
-}
-
-void dwpt(wpt_object wt, double *inp) {
-	int i, J, temp_len, iter, N, lp, p2, k, N2, Np;
-	int temp, elength, temp2,size,nodes,llb,n1,j;
-	double eparam,v1,v2;
-	int len_cA, t, t2, it1,it2;
-	double *orig,*tree;
-	int *nodelength;
-
-	temp_len = wt->siglength;
-	J = wt->J;
-	wt->length[J + 1] = temp_len;
-	wt->outlength = 0;
-	temp = 1;
-	elength = 0;
-	size = wt->wave->filtlength;
-	nodes = wt->nodes;
-	n1 = nodes + 1;
-	for (i = 0; i < J; ++i) {
-		temp *= 2;
-		temp2 = (size - 2) * (temp - 1);
-		elength += temp2;
-	}
-	eparam = wt->eparam;
-	orig = (double*)malloc(sizeof(double)* temp_len);
-	tree = (double*)malloc(sizeof(double)* (temp_len * (J + 1) + elength));
-	nodelength = (int*)malloc(sizeof(int)* nodes);
-
-	for (i = 0; i < wt->siglength; ++i) {
-		orig[i] = inp[i];
-	}
-
-	for (i = 0; i < temp_len * (J + 1) + elength; ++i) {
-		tree[i] = 0.0;
-	}
-
-	for (i = 0; i < nodes + 1; ++i) {
-		wt->basisvector[i] = 0.0;
-		wt->costvalues[i] = 0.0;
-	}
-
-	N = temp_len;
-	lp = wt->wave->lpd_len;
-	p2 = 1;
-
-	//set eparam value here
-	wt->costvalues[0] = costfunc(orig, wt->siglength, wt->entropy, eparam);
-	it2 = 1;
-	if (!strcmp(wt->ext, "per")) {
-		i = J;
-		p2 = 2;
-		while (i > 0) {
-			N = (int)ceil((double)N / 2.0);
-			wt->length[i] = N;
-			wt->outlength += p2 * (wt->length[i]);
-			i--;
-			p2 *= 2;
-		}
-		wt->length[0] = wt->length[1];
-
-		N2 = N = wt->outlength;
-		p2 = 1;
-		for (iter = 0; iter < J; ++iter) {
-			len_cA = wt->length[J - iter];
-			N2 -= 2 * p2 * len_cA;
-			N = N2;
-			for (k = 0; k < p2; ++k) {
-				if (iter == 0) {
-					dwpt_per(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
-				}
-				else {
-					dwpt_per(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
-				}
-				wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
-				it2++;
-				wt->costvalues[it2] = costfunc(tree + N +len_cA, len_cA, wt->entropy, eparam);
-				it2++;
-				N += 2 * len_cA;
-			}
-
-			temp_len = wt->length[J - iter];
-			p2 = 2 * p2;
-			Np = N2;
-		}
-	}
-	else if (!strcmp(wt->ext, "sym")) {
-		//printf("\n YES %s \n", wt->ext);
-		i = J;
-		p2 = 2;
-		while (i > 0) {
-			N = N + lp - 2;
-			N = (int)ceil((double)N / 2.0);
-			wt->length[i] = N;
-			wt->outlength += p2 * (wt->length[i]);
-			i--;
-			p2 *= 2;
-		}
-		wt->length[0] = wt->length[1];
-
-		N2 = N = wt->outlength;
-		p2 = 1;
-
-		for (iter = 0; iter < J; ++iter) {
-			len_cA = wt->length[J - iter];
-			N2 -= 2 * p2 * len_cA;
-			N = N2;
-			for (k = 0; k < p2; ++k) {
-				if (iter == 0) {
-					dwpt_sym(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
-				}
-				else {
-					dwpt_sym(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
-				}
-				wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
-				it2++;
-				wt->costvalues[it2] = costfunc(tree + N + len_cA, len_cA, wt->entropy, eparam);
-				it2++;
-				N += 2 * len_cA;
-			}
-
-			temp_len = wt->length[J - iter];
-			p2 = 2 * p2;
-			Np = N2;
-		}
-
-	}
-	else {
-		printf("Signal extension can be either per or sym");
-		exit(-1);
-	}
-
-	J = wt->J;
-	t2 = wt->outlength - 2 * wt->length[J];
-	p2 = 2;
-	it1 = 0;
-	for (i = 0; i < J; ++i) {
-		t = t2;
-		for (k = 0; k < p2; ++k) {
-			nodelength[it1] = t;
-			it1++;
-			t += wt->length[J - i];
-		}
-		p2 *= 2;
-		t2 = t2 - p2 * wt->length[J - i - 1];
-	}
-
-
-	J = wt->J;
-	llb = 1;
-	for (i = 0; i < J; ++i) {
-		llb *= 2;
-	}
-
-	for (i = n1 - llb; i < n1; ++i) {
-		wt->basisvector[i] = 1;
-	}
-
-	for (j = J - 1; j >= 0; --j) {
-		for (k = ipow2(j) - 1; k < ipow2(j + 1) - 1; ++k) {
-			v1 = wt->costvalues[k];
-			v2 = wt->costvalues[2 * k + 1] + wt->costvalues[2 * k + 2];
-			//printf(" %g %g", v1,v2);
-			if (v1 <= v2) {
-				wt->basisvector[k] = 1;
-			}
-			else {
-				wt->costvalues[k] = v2;
-			}
-		}
-		//printf("\n");
-	}
-
-	for (k = 0; k < nodes / 2; ++k) {
-		if (wt->basisvector[k] == 1 || wt->basisvector[k] == 2) {
-			wt->basisvector[2 * k + 1] = 2;
-			wt->basisvector[2 * k + 2] = 2;
-		}
-	}
-
-	for (k = 0; k < n1; ++k) {
-		if (wt->basisvector[k] == 2) {
-			wt->basisvector[k] = 0;
-		}
-	}
-
-	N2 = 0;
-	it1 = n1;
-	it2 = 0;
-	wt->nodes = 0;
-	wt->numnodeslevel[0] = 0;
-	//printf("Start \n");
-
-	if (wt->basisvector[0] == 1) {
-		wt->outlength = wt->siglength;
-		for (i = 0; i < wt->siglength; ++i) {
-			wt->output[i] = inp[i];
-		}
-		wt->nodes = 1;
-		wt->nodeindex[0] = 0;
-		wt->nodeindex[1] = 0;
-		wt->numnodeslevel[0] = 1;
-	}
-	else {
-		for (i = J; i > 0; --i) {
-			llb = ipow2(i);
-			it1 -= llb;
-			wt->numnodeslevel[i] = 0;
-			for (j = 0; j < llb; ++j) {
-				if (wt->basisvector[it1 + j] == 1) {
-					//printf("NODE %d %d %d \n", i, j, wt->length[J - i + 1]);
-					wt->nodeindex[2 * wt->nodes] = i;
-					wt->nodeindex[2 * wt->nodes + 1] = j;
-					wt->nodes += 1;
-					wt->numnodeslevel[i] += 1;
-					for (k = 0; k < wt->length[J - i + 1]; ++k) {
-						wt->output[it2 + k] = tree[nodelength[it1 - 1 + j] + k];// access tree
-					}
-					it2 += wt->length[J - i + 1];
-				}
-			}
-		}
-		wt->outlength = it2;
-	}
-
-	wt->coeflength[0] = wt->siglength;
-
-	for (i = 1; i < J + 1; ++i) {
-		wt->coeflength[i] = wt->length[J - i + 1];
-	}
-
-	free(orig);
-	free(tree);
-	free(nodelength);
-}
-
-int getWTREENodelength(wtree_object wt, int X) {
-	int N;
-	N = -1;
-	/*
-	X - Level. All Nodes at any level have the same length
-	*/
-	if (X <= 0 || X > wt->J) {
-		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
-		exit(-1);
-	}
-
-	N = wt->length[wt->J -X + 1];
-
-	return N;
-}
-
-int getDWPTNodelength(wpt_object wt, int X) {
-	int N;
-	N = -1;
-	/*
-	X - Level. All Nodes at any level have the same length
-	*/
-	if (X <= 0 || X > wt->J) {
-		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
-		exit(-1);
-	}
-
-	N = wt->length[wt->J - X + 1];
-
-	return N;
-}
-
-void getWTREECoeffs(wtree_object wt, int X,int Y,double *coeffs,int N) {
-	int ymax,i,t,t2;
-
-	if (X <= 0 || X > wt->J) {
-		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
-		exit(-1);
-	}
-	ymax = 1;
-	for (i = 0; i < X; ++i) {
-		ymax *= 2;
-	}
-
-	ymax -= 1;
-
-	if (Y < 0 ||Y > ymax) {
-		printf("Y co-ordinate must be >= 0 and <= %d", ymax);
-		exit(-1);
-	}
-
-	if (X == 1) {
-		t = 0;
-	}
-	else {
-		t = 0;
-		t2 = 1;
-		for (i = 0; i < X - 1; ++i) {
-			t2 *= 2;
-			t += t2;
-		}
-	}
-
-	t += Y;
-	t2 = wt->nodelength[t];
-	for (i = 0; i < N; ++i) {
-		coeffs[i] = wt->output[t2+i];
-	}
-
-}
-
-void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N) {
-	int ymax, i, t, t2;
-	int np,nds,citer;
-	int flag;
-
-	if (X <= 0 || X > wt->J) {
-		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
-		exit(-1);
-	}
-	ymax = 1;
-	for (i = 0; i < X; ++i) {
-		ymax *= 2;
-	}
-
-	ymax -= 1;
-
-	if (Y < 0 || Y > ymax) {
-		printf("Y co-ordinate must be >= 0 and <= %d", ymax);
-		exit(-1);
-	}
-
-	np = 0;
-	citer = 0;
-
-	for (i = wt->J; i > X; --i) {
-		np += wt->numnodeslevel[i];
-		citer += wt->numnodeslevel[i] * wt->coeflength[i];
-	}
-
-	i = 0;
-	flag = 0;
-	for (i = 0; i < wt->numnodeslevel[X]; ++i) {
-		if (wt->nodeindex[2 * np + 1] == Y) {
-			flag = 1;
-			break;
-		}
-		np++;
-		citer += wt->coeflength[X];
-	}
-
-	if (flag == 0) {
-		printf("The Node is Not Part Of The Best Basis Tree Use wpt_summary function to list available nodes \n");
-		exit(-1);
-	}
-
-	for (i = 0; i < N; ++i) {
-		coeffs[i] = wt->output[citer + i];
-	}
-
-}
-
-static void idwt1(wt_object wt,double *temp, double *cA_up,double *cA, int len_cA,double *cD,int len_cD,double *X_lp,double *X_hp,double *X) {
-	int len_avg, N, U,N2,i;
-
-	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-	N = 2 * len_cD;
-	U = 2;
-
-	upsamp2(cA, len_cA, U, cA_up);
-
-	per_ext(cA_up, 2 * len_cA, len_avg / 2, temp);
-
-	N2 = 2 * len_cA + len_avg;
-
-	if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-		wt->cobj = conv_init(N2, len_avg);
-		wt->cfftset = 1;
-	}
-	else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) {
-		printf("Decomposition Filters must have the same length.");
-		exit(-1);
-	}
-
-	wconv(wt, temp, N2, wt->wave->lpr, len_avg, X_lp);
-
-	upsamp2(cD, len_cD, U, cA_up);
-
-	per_ext(cA_up, 2 * len_cD, len_avg / 2, temp);
-
-	N2 = 2 * len_cD + len_avg;
-
-	wconv(wt, temp, N2, wt->wave->hpr, len_avg, X_hp);
-
-
-	for (i = len_avg - 1; i < N + len_avg - 1; ++i) {
-		X[i - len_avg + 1] = X_lp[i] + X_hp[i];
-	}
-
-	if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-		free_conv(wt->cobj);
-		wt->cfftset = 0;
-	}
-
-}
-
-static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_cD,  double *X) {
-	int len_avg,i,l,m,n,t,l2;
-
-	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-	l2 = len_avg / 2;
-	m = -2;
-	n = -1;
-
-	for (i = 0; i < len_cA + l2 - 1; ++i) {
-		m += 2;
-		n += 2;
-		X[m] = 0.0;
-		X[n] = 0.0;
-		for (l = 0; l < l2; ++l) {
-			t = 2 * l;
-			if ((i - l) >= 0 && (i - l) < len_cA) {
-				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
-				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
-			}
-			else if ((i - l) >= len_cA && (i-l) < len_cA + len_avg - 1) {
-				X[m] += wt->wave->lpr[t] * cA[i - l - len_cA] + wt->wave->hpr[t] * cD[i - l - len_cA];
-				X[n] += wt->wave->lpr[t + 1] * cA[i - l - len_cA] + wt->wave->hpr[t + 1] * cD[i - l - len_cA];
-			}
-			else if ((i - l) < 0 && (i-l) > -l2) {
-				X[m] += wt->wave->lpr[t] * cA[len_cA + i - l] + wt->wave->hpr[t] * cD[len_cA + i - l];
-				X[n] += wt->wave->lpr[t + 1] * cA[len_cA + i - l] + wt->wave->hpr[t + 1] * cD[len_cA + i - l];
-			}
-		}
-	}
-}
-
-static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) {
-	int len_avg, i, l, m, n, t, v;
-
-	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-	m = -2;
-	n = -1;
-
-	for (v = 0; v < len_cA; ++v) {
-		i = v;
-		m += 2;
-		n += 2;
-		X[m] = 0.0;
-		X[n] = 0.0;
-		for (l = 0; l < len_avg / 2; ++l) {
-			t = 2 * l;
-			if ((i - l) >= 0 && (i - l) < len_cA) {
-				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
-				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
-			}
-		}
-	}
-}
-
-
-void idwt(wt_object wt, double *dwtop) {
-	int J,U,i,lf,N,N2,iter,k;
-	int app_len, det_len;
-	double *cA_up, *X_lp, *X_hp,*out,*temp;
-
-	J = wt->J;
-	U = 2;
-	app_len = wt->length[0];
-	out = (double*)malloc(sizeof(double)* (wt->siglength + 1));
-	if (!strcmp(wt->ext, "per") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-		app_len = wt->length[0];
-		det_len = wt->length[1];
-		N = 2 * wt->length[J];
-		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-
-		cA_up = (double*)malloc(sizeof(double)* N);
-		temp = (double*)malloc(sizeof(double)* (N + lf));
-		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
-		X_hp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
-		iter = app_len;
-
-		for (i = 0; i < app_len; ++i) {
-			out[i] = wt->output[i];
-		}
-
-		for (i = 0; i < J; ++i) {
-
-			idwt1(wt,temp,cA_up,out,det_len,wt->output+iter,det_len,X_lp,X_hp,out);
-			/*
-			idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp);
-			for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
-				out[k - lf/2 + 1] = X_lp[k];
-			}
-			*/
-			iter += det_len;
-			det_len = wt->length[i + 2];
-		}
-		free(cA_up);
-		free(X_lp);
-		free(X_hp);
-		free(temp);
-
-	}
-	else if (!strcmp(wt->ext, "per") && !strcmp(wt->cmethod, "direct")) {
-		app_len = wt->length[0];
-		det_len = wt->length[1];
-		N = 2 * wt->length[J];
-		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-
-		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
-		iter = app_len;
-
-		for (i = 0; i < app_len; ++i) {
-			out[i] = wt->output[i];
-		}
-
-		for (i = 0; i < J; ++i) {
-
-			//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
-
-			idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp);
-			for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
-			out[k - lf/2 + 1] = X_lp[k];
-			}
-
-			iter += det_len;
-			det_len = wt->length[i + 2];
-		}
-
-		free(X_lp);
-
-	}
-	else if (!strcmp(wt->ext, "sym") && !strcmp(wt->cmethod, "direct")) {
-		app_len = wt->length[0];
-		det_len = wt->length[1];
-		N = 2 * wt->length[J] - 1;
-		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-
-		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
-		iter = app_len;
-
-		for (i = 0; i < app_len; ++i) {
-			out[i] = wt->output[i];
-		}
-
-		for (i = 0; i < J; ++i) {
-
-			//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
-
-			idwt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp);
-			for (k = lf-2; k < 2 * det_len; ++k) {
-				out[k - lf + 2] = X_lp[k];
-			}
-
-			iter += det_len;
-			det_len = wt->length[i + 2];
-		}
-
-		free(X_lp);
-
-	}
-	else if (!strcmp(wt->ext, "sym") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-		lf = wt->wave->lpd_len;// lpd and hpd have the same length
-
-		N = 2 * wt->length[J] - 1;
-		cA_up = (double*)malloc(sizeof(double)* N);
-		X_lp = (double*)malloc(sizeof(double)* (N + lf - 1));
-		X_hp = (double*)malloc(sizeof(double)* (N + lf - 1));
-
-		for (i = 0; i < app_len; ++i) {
-			out[i] = wt->output[i];
-		}
-
-		iter = app_len;
-
-		for (i = 0; i < J; ++i) {
-			det_len = wt->length[i + 1];
-			upsamp(out, det_len, U, cA_up);
-			N2 = 2 * wt->length[i + 1] - 1;
-
-			if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-				wt->cobj = conv_init(N2, lf);
-				wt->cfftset = 1;
-			}
-			else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) {
-				printf("Decomposition Filters must have the same length.");
-				exit(-1);
-			}
-
-			wconv(wt, cA_up, N2, wt->wave->lpr, lf, X_lp);
-
-			upsamp(wt->output + iter, det_len, U, cA_up);
-
-			wconv(wt, cA_up, N2, wt->wave->hpr, lf, X_hp);
-
-			for (k = lf - 2; k < N2 + 1; ++k) {
-				out[k - lf + 2] = X_lp[k] + X_hp[k];
-			}
-			iter += det_len;
-			if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-				free_conv(wt->cobj);
-				wt->cfftset = 0;
-			}
-		}
-
-		free(cA_up);
-		free(X_lp);
-		free(X_hp);
-	}
-	else {
-		printf("Signal extension can be either per or sym");
-		exit(-1);
-	}
-
-	for (i = 0; i < wt->siglength; ++i) {
-		dwtop[i] = out[i];
-	}
-
-
-	free(out);
-
-}
-
-static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) {
-	int len_avg, i, l, m, n, t, l2;
-
-	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-	l2 = len_avg / 2;
-	m = -2;
-	n = -1;
-
-	for (i = 0; i < len_cA + l2 - 1; ++i) {
-		m += 2;
-		n += 2;
-		X[m] = 0.0;
-		X[n] = 0.0;
-		for (l = 0; l < l2; ++l) {
-			t = 2 * l;
-			if ((i - l) >= 0 && (i - l) < len_cA) {
-				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
-				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
-			}
-			else if ((i - l) >= len_cA && (i - l) < len_cA + len_avg - 1) {
-				X[m] += wt->wave->lpr[t] * cA[i - l - len_cA] + wt->wave->hpr[t] * cD[i - l - len_cA];
-				X[n] += wt->wave->lpr[t + 1] * cA[i - l - len_cA] + wt->wave->hpr[t + 1] * cD[i - l - len_cA];
-			}
-			else if ((i - l) < 0 && (i - l) > -l2) {
-				X[m] += wt->wave->lpr[t] * cA[len_cA + i - l] + wt->wave->hpr[t] * cD[len_cA + i - l];
-				X[n] += wt->wave->lpr[t + 1] * cA[len_cA + i - l] + wt->wave->hpr[t + 1] * cD[len_cA + i - l];
-			}
-		}
-	}
-}
-
-static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) {
-	int len_avg, i, l, m, n, t, v;
-
-	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-	m = -2;
-	n = -1;
-
-	for (v = 0; v < len_cA; ++v) {
-		i = v;
-		m += 2;
-		n += 2;
-		X[m] = 0.0;
-		X[n] = 0.0;
-		for (l = 0; l < len_avg / 2; ++l) {
-			t = 2 * l;
-			if ((i - l) >= 0 && (i - l) < len_cA) {
-				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
-				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
-			}
-		}
-	}
-}
-
-void idwpt(wpt_object wt, double *dwtop) {
-	int J, U, i, lf, N, k,p,l;
-	int app_len, det_len, index, n1, llb, index2, index3, index4,indexp,xlen;
-	double *X_lp, *X,  *out, *out2;
-	int *prep,*ptemp;
-
-	J = wt->J;
-	U = 2;
-	app_len = wt->length[0];
-	p = ipow2(J);
-	lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
-	xlen = p * (app_len + 2 * lf);
-
-	X_lp = (double*)malloc(sizeof(double)* 2 * (wt->length[J] + lf));
-	X = (double*)malloc(sizeof(double)* xlen);
-	out = (double*)malloc(sizeof(double)* wt->length[J]);
-	out2 = (double*)malloc(sizeof(double)* wt->length[J]);
-	prep = (int*)malloc(sizeof(int)* p);
-	ptemp = (int*)malloc(sizeof(int)* p);
-	n1 = 1;
-	llb = 1;
-	index2 = xlen / p;
-	indexp = 0;
-	
+/*
+  Copyright (c) 2014, Rafat Hussain
+*/
+
+#include "wavelib.h"
+
+wave_object wave_init(char* wname) {
+	wave_object obj = NULL;
+	int retval;
+	retval = 0;
+
+	if (wname != NULL) {
+		retval = filtlength(wname);
+		//obj->filtlength = retval;
+		//strcopy(obj->wname, wname);
+	}
+
+	obj = (wave_object)malloc(sizeof(struct wave_set) + sizeof(double)* 4 * retval);
+
+	obj->filtlength = retval;
+	obj->lpd_len = obj->hpd_len = obj->lpr_len = obj->hpr_len = obj->filtlength;
+	strcpy(obj->wname, wname);
+	if (wname != NULL) {
+		filtcoef(wname,obj->params,obj->params+retval,obj->params+2*retval,obj->params+3*retval);
+	}
+	obj->lpd = &obj->params[0];
+	obj->hpd = &obj->params[retval];
+	obj->lpr = &obj->params[2 * retval];
+	obj->hpr = &obj->params[3 * retval];
+	return obj;
+}
+
+wt_object wt_init(wave_object wave,char* method, int siglength,int J) {
+	int size,i,MaxIter;
+	wt_object obj = NULL;
+
+	size = wave->filtlength;
+
+	if (J > 100) {
+		printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
+		exit(-1);
+	}
+
+	MaxIter = wmaxiter(siglength, size);
+
+	if (J > MaxIter) {
+		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter);
+		exit(-1);
+	}
+
+	if (method == NULL) {
+		obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength +  2 * J * (size+1)));
+		obj->outlength = siglength + 2 * J * (size + 1); // Default
+		strcpy(obj->ext, "sym"); // Default
+	}
+	else if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
+		obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength + 2 * J * (size + 1)));
+		obj->outlength = siglength + 2 * J * (size + 1); // Default
+		strcpy(obj->ext, "sym"); // Default
+	}
+	else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) {
+		if (!testSWTlength(siglength, J)) {
+			printf("\n For SWT the signal length must be a multiple of 2^J. \n");
+			exit(-1);
+		}
+
+		obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1)));
+		obj->outlength = siglength * (J + 1); // Default
+		strcpy(obj->ext, "per"); // Default
+	}
+	else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
+
+		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 = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1)));
+		obj->outlength = siglength * (J + 1); // Default
+		strcpy(obj->ext, "per"); // Default
+	}
+
+	obj->wave = wave;
+	obj->siglength = siglength;
+	obj->J = J;
+	obj->MaxIter = MaxIter;
+	strcpy(obj->method, method);
+
+	if (siglength % 2 == 0) {
+		obj->even = 1;
+	}
+	else {
+		obj->even = 0;
+	}
+
+	obj->cobj = NULL;
+
+	strcpy(obj->cmethod, "direct"); // Default
+	obj->cfftset = 0;
+	obj->lenlength = J + 2;
+	obj->output = &obj->params[0];
+	if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
+		for (i = 0; i < siglength + 2 * J * (size + 1); ++i) {
+			obj->params[i] = 0.0;
+		}
+	}
+	else if (!strcmp(method, "swt") || !strcmp(method, "SWT") || !strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
+		for (i = 0; i < siglength * (J + 1); ++i) {
+			obj->params[i] = 0.0;
+		}
+	}
+	//wave_summary(obj->wave);
+
+	return obj;
+}
+
+wtree_object wtree_init(wave_object wave, int siglength,int J) {
+    int size,i,MaxIter,temp,temp2,elength,nodes;
+	wtree_object obj = NULL;
+
+	size = wave->filtlength;
+
+	if (J > 100) {
+		printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
+		exit(-1);
+	}
+
+
+	MaxIter = wmaxiter(siglength, size);
+	if (J > MaxIter) {
+		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
+		exit(-1);
+	}
+	temp = 1;
+	elength = 0;
+	nodes = 0;
+	for(i = 0; i < J;++i) {
+	  temp *= 2;
+	  nodes += temp;
+	  temp2 = (size - 2) * (temp - 1);
+	  elength += temp2;
+	}
+
+	obj = (wtree_object)malloc(sizeof(struct wtree_set) + sizeof(double)* (siglength * (J + 1) + elength + nodes + J + 1));
+	obj->outlength = siglength * (J + 1) + elength;
+	strcpy(obj->ext, "sym");
+
+	obj->wave = wave;
+	obj->siglength = siglength;
+	obj->J = J;
+	obj->MaxIter = MaxIter;
+	strcpy(obj->method, "dwt");
+
+	if (siglength % 2 == 0) {
+		obj->even = 1;
+	}
+	else {
+		obj->even = 0;
+	}
+
+	obj->cobj = NULL;
+	obj->nodes = nodes;
+
+	obj->cfftset = 0;
+	obj->lenlength = J + 2;
+	obj->output = &obj->params[0];
+	obj->nodelength = (int*) &obj->params[siglength * (J + 1) + elength];
+	obj->coeflength = (int*)&obj->params[siglength * (J + 1) + elength + nodes];
+
+	for (i = 0; i < siglength * (J + 1) + elength + nodes + J + 1; ++i) {
+	       obj->params[i] = 0.0;
+	}
+
+	//wave_summary(obj->wave);
+
+	return obj;
+}
+
+wpt_object wpt_init(wave_object wave, int siglength, int J) {
+	int size, i, MaxIter, temp, nodes,elength,p2,N,lp;
+	wpt_object obj = NULL;
+
+	size = wave->filtlength;
+
+	if (J > 100) {
+		printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
+		exit(-1);
+	}
+
+
+	MaxIter = wmaxiter(siglength, size);
+	if (J > MaxIter) {
+		printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
+		exit(-1);
+	}
+	temp = 1;
+	nodes = 0;
+	for (i = 0; i < J; ++i) {
+		temp *= 2;
+		nodes += temp;
+	}
+
+	i = J;
+	p2 = 2;
+	N = siglength;
+	lp = size;
+	elength = 0;
+	while (i > 0) {
+		N = N + lp - 2;
+		N = (int)ceil((double)N / 2.0);
+		elength = p2 * N;
+		i--;
+		p2 *= 2;
+	}
+	//printf("elength %d", elength);
+
+	obj = (wpt_object)malloc(sizeof(struct wpt_set) + sizeof(double)* (elength + 4 * nodes + 2 * J + 6));
+	obj->outlength = siglength + 2 * (J + 1) * (size + 1);
+	strcpy(obj->ext, "sym");
+	strcpy(obj->entropy, "shannon");
+	obj->eparam = 0.0;
+
+	obj->wave = wave;
+	obj->siglength = siglength;
+	obj->J = J;
+	obj->MaxIter = MaxIter;
+
+	if (siglength % 2 == 0) {
+		obj->even = 1;
+	}
+	else {
+		obj->even = 0;
+	}
+
+	obj->cobj = NULL;
+	obj->nodes = nodes;
+
+	obj->lenlength = J + 2;
+	obj->output = &obj->params[0];
+	obj->costvalues = &obj->params[elength];
+	obj->basisvector = &obj->params[elength + nodes + 1];
+	obj->nodeindex = (int*)&obj->params[elength + 2*nodes + 2];
+	obj->numnodeslevel = (int*)&obj->params[elength + 4 * nodes + 4];
+	obj->coeflength = (int*)&obj->params[elength + 4 * nodes + J + 5];
+
+	for (i = 0; i < elength + 4 * nodes + 2 * J + 6; ++i) {
+		obj->params[i] = 0.0;
+	}
+
+	//wave_summary(obj->wave);
+
+	return obj;
+}
+
+cwt_object cwt_init(char* wave, double param,int siglength, double dt, int J) {
+	cwt_object obj = NULL;
+	int N, i,nj2,ibase2,mother;
+	double s0, dj;
+	double t1;
+	int m, odd;
+	char *pdefault = "pow";
+
+	m = (int)param;
+	odd = 1;
+	if (2 * (m / 2) == m) {
+		odd = 0;
+	}
+
+	N = siglength;
+	nj2 = 2 * N * J;
+	obj = (cwt_object)malloc(sizeof(struct cwt_set) + sizeof(double)* (nj2 + 2 * J + N));
+
+	if (!strcmp(wave, "morlet") || !strcmp(wave, "morl")) {
+		s0 = 2 * dt;
+		dj = 0.4875;
+		mother = 0;
+		if (param < 0.0) {
+			printf("\n Morlet Wavelet Parameter should be >= 0 \n");
+			exit(-1);
+		}
+		if (param == 0) {
+			param = 6.0;
+		}
+		strcpy(obj->wave,"morlet");
+		
+	}
+	else if (!strcmp(wave, "paul")) {
+		s0 = 2 * dt;
+		dj = 0.4875;
+		mother = 1;
+		if (param < 0 || param > 20) {
+			printf("\n Paul Wavelet Parameter should be > 0 and <= 20 \n");
+			exit(-1);
+		}
+		if (param == 0) {
+			param = 4.0;
+		}
+		strcpy(obj->wave,"paul");
+	
+	}
+	else if (!strcmp(wave, "dgauss") || !strcmp(wave, "dog")) {
+		s0 = 2 * dt;
+		dj = 0.4875;
+		mother = 2;
+		if (param < 0 || odd == 1) {
+			printf("\n DOG Wavelet Parameter should be > 0 and even \n");
+			exit(-1);
+		}
+		if (param == 0) {
+			param = 2.0;
+		}
+		strcpy(obj->wave,"dog");
+	}
+
+	obj->pow = 2;
+	strcpy(obj->type, pdefault);
+
+	obj->s0 = s0;
+	obj->dj = dj;
+	obj->dt = dt;
+	obj->J = J;
+	obj->siglength = siglength;
+	obj->sflag = 0;
+	obj->pflag = 1;
+	obj->mother = mother;
+	obj->m = param;
+
+	t1 = 0.499999 + log((double)N) / log(2.0);
+	ibase2 = 1 + (int)t1;
+
+	obj->npad = (int)pow(2.0, (double)ibase2);
+
+	obj->output = (cplx_data*) &obj->params[0];
+	obj->scale = &obj->params[nj2];
+	obj->period = &obj->params[nj2 + J];
+	obj->coi = &obj->params[nj2 + 2*J];
+
+	for (i = 0; i < nj2 + 2 * J + N; ++i) {
+		obj->params[i] = 0.0;
+	}
+
+	return obj;
+}
+
+
+static void wconv(wt_object wt, double *sig, int N, double *filt, int L, double *oup) {
+	if (!strcmp(wt->cmethod,"direct")) {
+		conv_direct(sig, N, filt, L, oup);
+	}
+	else if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) {
+		if (wt->cfftset == 0) {
+			wt->cobj = conv_init(N, L);
+			conv_fft(wt->cobj, sig, filt, oup);
+			free_conv(wt->cobj);
+		}
+		else {
+			conv_fft(wt->cobj, sig, filt, oup);
+		}
+	}
+	else {
+		printf("Convolution Only accepts two methods - direct and fft");
+		exit(-1);
+	}
+}
+
+
+static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int l,l2,isodd,i,t,len_avg;
+
+	len_avg = wt->wave->lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + l2;
+		cA[i] = 0.0;
+		cD[i] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= l2 && (t - l) < N) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < l2 && (t - l) >= 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < 0 && isodd == 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l + N];
+				cD[i] += wt->wave->hpd[l] * inp[t - l + N];
+			}
+			else if ((t - l) < 0 && isodd == 1) {
+				if ((t - l) != -1) {
+					cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
+					cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
+				}
+				else {
+					cA[i] += wt->wave->lpd[l] * inp[N - 1];
+					cD[i] += wt->wave->hpd[l] * inp[N - 1];
+				}
+			}
+			else if ((t - l) >= N && isodd == 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l - N];
+				cD[i] += wt->wave->hpd[l] * inp[t - l - N];
+			}
+			else if ((t - l) >= N && isodd == 1) {
+				if (t - l != N) {
+					cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
+					cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
+				}
+				else {
+					cA[i] += wt->wave->lpd[l] * inp[N - 1];
+					cD[i] += wt->wave->hpd[l] * inp[N - 1];
+				}
+			}
+
+		}
+	}
+
+
+
+}
+
+static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int l, l2, isodd, i, t, len_avg;
+
+	len_avg = wt->wave->lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + l2;
+		cA[i] = 0.0;
+		cD[i] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= l2 && (t - l) < N) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < l2 && (t - l) >= 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < 0 && isodd == 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l + N];
+				cD[i] += wt->wave->hpd[l] * inp[t - l + N];
+			}
+			else if ((t - l) < 0 && isodd == 1) {
+				if ((t - l) != -1) {
+					cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
+					cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
+				}
+				else {
+					cA[i] += wt->wave->lpd[l] * inp[N - 1];
+					cD[i] += wt->wave->hpd[l] * inp[N - 1];
+				}
+			}
+			else if ((t - l) >= N && isodd == 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l - N];
+				cD[i] += wt->wave->hpd[l] * inp[t - l - N];
+			}
+			else if ((t - l) >= N && isodd == 1) {
+				if (t - l != N) {
+					cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
+					cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
+				}
+				else {
+					cA[i] += wt->wave->lpd[l] * inp[N - 1];
+					cD[i] += wt->wave->hpd[l] * inp[N - 1];
+				}
+			}
+
+		}
+	}
+
+
+
+}
+
+static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int l, l2, isodd, i, t, len_avg;
+
+	len_avg = wt->wave->lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + l2;
+		cA[i] = 0.0;
+		cD[i] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= l2 && (t - l) < N) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < l2 && (t - l) >= 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < 0 && isodd == 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l + N];
+				cD[i] += wt->wave->hpd[l] * inp[t - l + N];
+			}
+			else if ((t - l) < 0 && isodd == 1) {
+				if ((t - l) != -1) {
+					cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
+					cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
+				}
+				else {
+					cA[i] += wt->wave->lpd[l] * inp[N - 1];
+					cD[i] += wt->wave->hpd[l] * inp[N - 1];
+				}
+			}
+			else if ((t - l) >= N && isodd == 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l - N];
+				cD[i] += wt->wave->hpd[l] * inp[t - l - N];
+			}
+			else if ((t - l) >= N && isodd == 1) {
+				if (t - l != N) {
+					cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
+					cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
+				}
+				else {
+					cA[i] += wt->wave->lpd[l] * inp[N - 1];
+					cD[i] += wt->wave->hpd[l] * inp[N - 1];
+				}
+			}
+
+		}
+	}
+
+
+
+}
+
+static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int i,l, t, len_avg;
+
+	len_avg = wt->wave->lpd_len;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + 1;
+		cA[i] = 0.0;
+		cD[i] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= 0 && (t - l) < N) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < 0) {
+				cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
+				cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
+			}
+			else if ((t - l) >= N) {
+				cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
+				cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
+			}
+
+		}
+	}
+
+
+}
+
+static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int i, l, t, len_avg;
+
+	len_avg = wt->wave->lpd_len;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + 1;
+		cA[i] = 0.0;
+		cD[i] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= 0 && (t - l) < N) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < 0) {
+				cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
+				cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
+			}
+			else if ((t - l) >= N) {
+				cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
+				cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
+			}
+
+		}
+	}
+
+
+}
+
+static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int i, l, t, len_avg;
+
+	len_avg = wt->wave->lpd_len;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = 2 * i + 1;
+		cA[i] = 0.0;
+		cD[i] = 0.0;
+		for (l = 0; l < len_avg; ++l) {
+			if ((t - l) >= 0 && (t - l) < N) {
+				cA[i] += wt->wave->lpd[l] * inp[t - l];
+				cD[i] += wt->wave->hpd[l] * inp[t - l];
+			}
+			else if ((t - l) < 0) {
+				cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
+				cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
+			}
+			else if ((t - l) >= N) {
+				cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
+				cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
+			}
+
+		}
+	}
+
+
+}
+
+static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, int len_cA, double *cD, int len_cD) {
+	int len_avg,D,lf;
+	double *signal,*cA_undec;
+	len_avg = (wt->wave->lpd_len + wt->wave->hpd_len) / 2;
+	//len_sig = 2 * (int)ceil((double)len_sig / 2.0);
+
+	D = 2;
+
+	if (!strcmp(wt->ext, "per")) {
+		signal = (double*)malloc(sizeof(double)* (len_sig + len_avg + (len_sig % 2)));
+
+		len_sig = per_ext(sig, len_sig, len_avg / 2, signal);
+
+		cA_undec = (double*)malloc(sizeof(double)* (len_sig + len_avg + wt->wave->lpd_len - 1));
+
+		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+			wt->cobj = conv_init(len_sig + len_avg, wt->wave->lpd_len);
+			wt->cfftset = 1;
+		}
+		else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
+			printf("Decomposition Filters must have the same length.");
+			exit(-1);
+		}
+
+		wconv(wt, signal, len_sig + len_avg, wt->wave->lpd, wt->wave->lpd_len, cA_undec);
+
+		downsamp(cA_undec + len_avg, len_sig, D, cA);
+
+		wconv(wt, signal, len_sig + len_avg, wt->wave->hpd, wt->wave->hpd_len, cA_undec);
+
+		downsamp(cA_undec + len_avg, len_sig, D, cD);
+	}
+	else if (!strcmp(wt->ext, "sym")) {
+		//printf("\n YES %s \n", wt->ext);
+		lf = wt->wave->lpd_len;// lpd and hpd have the same length
+
+		signal = (double*)malloc(sizeof(double)* (len_sig + 2 * (lf - 1)));
+
+		len_sig = symm_ext(sig, len_sig, lf - 1, signal);
+
+		cA_undec = (double*)malloc(sizeof(double)* (len_sig + 3 * (lf - 1)));
+
+		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+			wt->cobj = conv_init(len_sig + 2 * (lf - 1), lf);
+			wt->cfftset = 1;
+		}
+		else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
+			printf("Decomposition Filters must have the same length.");
+			exit(-1);
+		}
+
+
+		wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->lpd, wt->wave->lpd_len, cA_undec);
+
+		downsamp(cA_undec + lf, len_sig + lf - 2, D, cA);
+
+		wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->hpd, wt->wave->hpd_len, cA_undec);
+
+		downsamp(cA_undec + lf, len_sig + lf - 2, D, cD);
+	}
+	else {
+		printf("Signal extension can be either per or sym");
+		exit(-1);
+	}
+
+
+	if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+		free_conv(wt->cobj);
+		wt->cfftset = 0;
+	}
+
+	free(signal);
+	free(cA_undec);
+}
+
+void dwt(wt_object wt,double *inp) {
+	int i,J,temp_len,iter,N,lp;
+	int len_cA;
+	double *orig,*orig2;
+
+	temp_len = wt->siglength;
+	J = wt->J;
+	wt->length[J + 1] = temp_len;
+	wt->outlength = 0;
+	wt->zpad = 0;
+	orig = (double*)malloc(sizeof(double)* temp_len);
+	orig2 = (double*)malloc(sizeof(double)* temp_len);
+	/*
+	if ((temp_len % 2) == 0) {
+	wt->zpad = 0;
+	orig = (double*)malloc(sizeof(double)* temp_len);
+	orig2 = (double*)malloc(sizeof(double)* temp_len);
+	}
+	else {
+	wt->zpad = 1;
+	temp_len++;
+	orig = (double*)malloc(sizeof(double)* temp_len);
+	orig2 = (double*)malloc(sizeof(double)* temp_len);
+	}
+	*/
+
+	for (i = 0; i < wt->siglength; ++i) {
+		orig[i] = inp[i];
+	}
+
+	if (wt->zpad == 1) {
+		orig[temp_len - 1] = orig[temp_len - 2];
+	}
+
+	N = temp_len;
+	lp = wt->wave->lpd_len;
+
+	if (!strcmp(wt->ext,"per")) {
+		i = J;
+		while (i > 0) {
+			N = (int)ceil((double)N / 2.0);
+			wt->length[i] = N;
+			wt->outlength += wt->length[i];
+			i--;
+		}
+		wt->length[0] = wt->length[1];
+		wt->outlength += wt->length[0];
+		N = wt->outlength;
+
+		for (iter = 0; iter < J; ++iter) {
+			len_cA = wt->length[J - iter];
+			N -= len_cA;
+			if ( !strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT") ) {
+				dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
+			}
+			else {
+				dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
+			}
+			temp_len = wt->length[J - iter];
+			if (iter == J - 1) {
+				for (i = 0; i < len_cA; ++i) {
+					wt->params[i] = orig2[i];
+				}
+			}
+			else {
+				for (i = 0; i < len_cA; ++i) {
+					orig[i] = orig2[i];
+				}
+			}
+		}
+	}
+	else if (!strcmp(wt->ext,"sym")) {
+		//printf("\n YES %s \n", wt->ext);
+		i = J;
+		while (i > 0) {
+			N = N + lp - 2;
+			N = (int) ceil((double)N / 2.0);
+			wt->length[i] = N;
+			wt->outlength += wt->length[i];
+			i--;
+		}
+		wt->length[0] = wt->length[1];
+		wt->outlength += wt->length[0];
+		N = wt->outlength;
+
+		for (iter = 0; iter < J; ++iter) {
+			len_cA = wt->length[J - iter];
+			N -= len_cA;
+			if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) {
+				dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
+			}
+			else {
+				dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA);
+			}
+			temp_len = wt->length[J - iter];
+
+			if (iter == J - 1) {
+				for (i = 0; i < len_cA; ++i) {
+					wt->params[i] = orig2[i];
+				}
+			}
+			else {
+				for (i = 0; i < len_cA; ++i) {
+					orig[i] = orig2[i];
+				}
+			}
+		}
+	}
+	else {
+		printf("Signal extension can be either per or sym");
+		exit(-1);
+	}
+
+	free(orig);
+	free(orig2);
+}
+
+void wtree(wtree_object wt,double *inp) {
+	int i,J,temp_len,iter,N,lp,p2,k,N2,Np;
+	int len_cA,t,t2,it1;
+	double *orig;
+
+	temp_len = wt->siglength;
+	J = wt->J;
+	wt->length[J + 1] = temp_len;
+	wt->outlength = 0;
+	wt->zpad = 0;
+	orig = (double*)malloc(sizeof(double)* temp_len);
+	/*
+	if ((temp_len % 2) == 0) {
+		wt->zpad = 0;
+		orig = (double*)malloc(sizeof(double)* temp_len);
+	}
+	else {
+		wt->zpad = 1;
+		temp_len++;
+		orig = (double*)malloc(sizeof(double)* temp_len);
+	}
+	*/
+	for (i = 0; i < wt->siglength; ++i) {
+		orig[i] = inp[i];
+	}
+
+	if (wt->zpad == 1) {
+		orig[temp_len - 1] = orig[temp_len - 2];
+	}
+
+	N = temp_len;
+	lp = wt->wave->lpd_len;
+        p2 = 1;
+
+	if (!strcmp(wt->ext,"per")) {
+		i = J;
+                p2 = 2;
+		while (i > 0) {
+			N = (int)ceil((double)N / 2.0);
+			wt->length[i] = N;
+			wt->outlength += p2 * (wt->length[i]);
+			i--;
+                        p2 *= 2;
+		}
+		wt->length[0] = wt->length[1];
+
+		N2 = N = wt->outlength;
+                p2 = 1;
+		for (iter = 0; iter < J; ++iter) {
+			len_cA = wt->length[J - iter];
+                        N2 -= 2 * p2 * len_cA;
+                        N = N2;
+                        for(k = 0; k < p2;++k) {
+                            if (iter == 0) {
+                               wtree_per(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
+                            } else {
+                                wtree_per(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
+                            }
+                            N += 2 * len_cA;
+                        }
+
+			temp_len = wt->length[J - iter];
+			p2 = 2 * p2;
+            Np = N2;
+		}
+	}
+	else if (!strcmp(wt->ext,"sym")) {
+		//printf("\n YES %s \n", wt->ext);
+		i = J;
+                p2 = 2;
+		while (i > 0) {
+			N = N + lp - 2;
+			N = (int) ceil((double)N / 2.0);
+			wt->length[i] = N;
+			wt->outlength += p2 * (wt->length[i]);
+			i--;
+                        p2 *= 2;
+		}
+		wt->length[0] = wt->length[1];
+
+		N2 = N = wt->outlength;
+                p2 = 1;
+
+                for (iter = 0; iter < J; ++iter) {
+			len_cA = wt->length[J - iter];
+                        N2 -= 2 * p2 * len_cA;
+                        N = N2;
+                        for(k = 0; k < p2;++k) {
+                            if (iter == 0) {
+                                wtree_sym(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
+                            } else {
+                                wtree_sym(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA);
+                            }
+                            N += 2 * len_cA;
+                        }
+
+			temp_len = wt->length[J - iter];
+			p2 = 2 * p2;
+            Np = N2;
+		}
+
+	}
+	else {
+		printf("Signal extension can be either per or sym");
+		exit(-1);
+	}
+
+	J = wt->J;
+	t2 = wt->outlength - 2 * wt->length[J];
+	p2 = 2;
+	it1 = 0;
+	for (i = 0; i < J; ++i) {
+		t = t2;
+		for (k = 0; k < p2; ++k) {
+			wt->nodelength[it1] = t;
+			it1++;
+			t += wt->length[J - i];
+		}
+		p2 *= 2;
+		t2 = t2 - p2 * wt->length[J - i - 1];
+	}
+
+	wt->coeflength[0] = wt->siglength;
+
+	for (i = 1; i < J + 1; ++i) {
+		wt->coeflength[i] = wt->length[J - i + 1];
+	}
+
+	free(orig);
+}
+
+static int ipow2(int n) {
+	int p,i;
+	p = 1;
+
+	for (i = 0; i < n; ++i) {
+		p *= 2;
+	}
+
+	return p;
+}
+
+void dwpt(wpt_object wt, double *inp) {
+	int i, J, temp_len, iter, N, lp, p2, k, N2, Np;
+	int temp, elength, temp2,size,nodes,llb,n1,j;
+	double eparam,v1,v2;
+	int len_cA, t, t2, it1,it2;
+	double *orig,*tree;
+	int *nodelength;
+
+	temp_len = wt->siglength;
+	J = wt->J;
+	wt->length[J + 1] = temp_len;
+	wt->outlength = 0;
+	temp = 1;
+	elength = 0;
+	size = wt->wave->filtlength;
+	nodes = wt->nodes;
+	n1 = nodes + 1;
+	for (i = 0; i < J; ++i) {
+		temp *= 2;
+		temp2 = (size - 2) * (temp - 1);
+		elength += temp2;
+	}
+	eparam = wt->eparam;
+	orig = (double*)malloc(sizeof(double)* temp_len);
+	tree = (double*)malloc(sizeof(double)* (temp_len * (J + 1) + elength));
+	nodelength = (int*)malloc(sizeof(int)* nodes);
+
+	for (i = 0; i < wt->siglength; ++i) {
+		orig[i] = inp[i];
+	}
+
+	for (i = 0; i < temp_len * (J + 1) + elength; ++i) {
+		tree[i] = 0.0;
+	}
+
+	for (i = 0; i < nodes + 1; ++i) {
+		wt->basisvector[i] = 0.0;
+		wt->costvalues[i] = 0.0;
+	}
+
+	N = temp_len;
+	lp = wt->wave->lpd_len;
+	p2 = 1;
+
+	//set eparam value here
+	wt->costvalues[0] = costfunc(orig, wt->siglength, wt->entropy, eparam);
+	it2 = 1;
+	if (!strcmp(wt->ext, "per")) {
+		i = J;
+		p2 = 2;
+		while (i > 0) {
+			N = (int)ceil((double)N / 2.0);
+			wt->length[i] = N;
+			wt->outlength += p2 * (wt->length[i]);
+			i--;
+			p2 *= 2;
+		}
+		wt->length[0] = wt->length[1];
+
+		N2 = N = wt->outlength;
+		p2 = 1;
+		for (iter = 0; iter < J; ++iter) {
+			len_cA = wt->length[J - iter];
+			N2 -= 2 * p2 * len_cA;
+			N = N2;
+			for (k = 0; k < p2; ++k) {
+				if (iter == 0) {
+					dwpt_per(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
+				}
+				else {
+					dwpt_per(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
+				}
+				wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
+				it2++;
+				wt->costvalues[it2] = costfunc(tree + N +len_cA, len_cA, wt->entropy, eparam);
+				it2++;
+				N += 2 * len_cA;
+			}
+
+			temp_len = wt->length[J - iter];
+			p2 = 2 * p2;
+			Np = N2;
+		}
+	}
+	else if (!strcmp(wt->ext, "sym")) {
+		//printf("\n YES %s \n", wt->ext);
+		i = J;
+		p2 = 2;
+		while (i > 0) {
+			N = N + lp - 2;
+			N = (int)ceil((double)N / 2.0);
+			wt->length[i] = N;
+			wt->outlength += p2 * (wt->length[i]);
+			i--;
+			p2 *= 2;
+		}
+		wt->length[0] = wt->length[1];
+
+		N2 = N = wt->outlength;
+		p2 = 1;
+
+		for (iter = 0; iter < J; ++iter) {
+			len_cA = wt->length[J - iter];
+			N2 -= 2 * p2 * len_cA;
+			N = N2;
+			for (k = 0; k < p2; ++k) {
+				if (iter == 0) {
+					dwpt_sym(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
+				}
+				else {
+					dwpt_sym(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA);
+				}
+				wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
+				it2++;
+				wt->costvalues[it2] = costfunc(tree + N + len_cA, len_cA, wt->entropy, eparam);
+				it2++;
+				N += 2 * len_cA;
+			}
+
+			temp_len = wt->length[J - iter];
+			p2 = 2 * p2;
+			Np = N2;
+		}
+
+	}
+	else {
+		printf("Signal extension can be either per or sym");
+		exit(-1);
+	}
+
+	J = wt->J;
+	t2 = wt->outlength - 2 * wt->length[J];
+	p2 = 2;
+	it1 = 0;
+	for (i = 0; i < J; ++i) {
+		t = t2;
+		for (k = 0; k < p2; ++k) {
+			nodelength[it1] = t;
+			it1++;
+			t += wt->length[J - i];
+		}
+		p2 *= 2;
+		t2 = t2 - p2 * wt->length[J - i - 1];
+	}
+
+
+	J = wt->J;
+	llb = 1;
+	for (i = 0; i < J; ++i) {
+		llb *= 2;
+	}
+
+	for (i = n1 - llb; i < n1; ++i) {
+		wt->basisvector[i] = 1;
+	}
+
+	for (j = J - 1; j >= 0; --j) {
+		for (k = ipow2(j) - 1; k < ipow2(j + 1) - 1; ++k) {
+			v1 = wt->costvalues[k];
+			v2 = wt->costvalues[2 * k + 1] + wt->costvalues[2 * k + 2];
+			//printf(" %g %g", v1,v2);
+			if (v1 <= v2) {
+				wt->basisvector[k] = 1;
+			}
+			else {
+				wt->costvalues[k] = v2;
+			}
+		}
+		//printf("\n");
+	}
+
+	for (k = 0; k < nodes / 2; ++k) {
+		if (wt->basisvector[k] == 1 || wt->basisvector[k] == 2) {
+			wt->basisvector[2 * k + 1] = 2;
+			wt->basisvector[2 * k + 2] = 2;
+		}
+	}
+
+	for (k = 0; k < n1; ++k) {
+		if (wt->basisvector[k] == 2) {
+			wt->basisvector[k] = 0;
+		}
+	}
+
+	N2 = 0;
+	it1 = n1;
+	it2 = 0;
+	wt->nodes = 0;
+	wt->numnodeslevel[0] = 0;
+	//printf("Start \n");
+
+	if (wt->basisvector[0] == 1) {
+		wt->outlength = wt->siglength;
+		for (i = 0; i < wt->siglength; ++i) {
+			wt->output[i] = inp[i];
+		}
+		wt->nodes = 1;
+		wt->nodeindex[0] = 0;
+		wt->nodeindex[1] = 0;
+		wt->numnodeslevel[0] = 1;
+	}
+	else {
+		for (i = J; i > 0; --i) {
+			llb = ipow2(i);
+			it1 -= llb;
+			wt->numnodeslevel[i] = 0;
+			for (j = 0; j < llb; ++j) {
+				if (wt->basisvector[it1 + j] == 1) {
+					//printf("NODE %d %d %d \n", i, j, wt->length[J - i + 1]);
+					wt->nodeindex[2 * wt->nodes] = i;
+					wt->nodeindex[2 * wt->nodes + 1] = j;
+					wt->nodes += 1;
+					wt->numnodeslevel[i] += 1;
+					for (k = 0; k < wt->length[J - i + 1]; ++k) {
+						wt->output[it2 + k] = tree[nodelength[it1 - 1 + j] + k];// access tree
+					}
+					it2 += wt->length[J - i + 1];
+				}
+			}
+		}
+		wt->outlength = it2;
+	}
+
+	wt->coeflength[0] = wt->siglength;
+
+	for (i = 1; i < J + 1; ++i) {
+		wt->coeflength[i] = wt->length[J - i + 1];
+	}
+
+	free(orig);
+	free(tree);
+	free(nodelength);
+}
+
+int getWTREENodelength(wtree_object wt, int X) {
+	int N;
+	N = -1;
+	/*
+	X - Level. All Nodes at any level have the same length
+	*/
+	if (X <= 0 || X > wt->J) {
+		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
+		exit(-1);
+	}
+
+	N = wt->length[wt->J -X + 1];
+
+	return N;
+}
+
+int getDWPTNodelength(wpt_object wt, int X) {
+	int N;
+	N = -1;
+	/*
+	X - Level. All Nodes at any level have the same length
+	*/
+	if (X <= 0 || X > wt->J) {
+		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
+		exit(-1);
+	}
+
+	N = wt->length[wt->J - X + 1];
+
+	return N;
+}
+
+void getWTREECoeffs(wtree_object wt, int X,int Y,double *coeffs,int N) {
+	int ymax,i,t,t2;
+
+	if (X <= 0 || X > wt->J) {
+		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
+		exit(-1);
+	}
+	ymax = 1;
+	for (i = 0; i < X; ++i) {
+		ymax *= 2;
+	}
+
+	ymax -= 1;
+
+	if (Y < 0 ||Y > ymax) {
+		printf("Y co-ordinate must be >= 0 and <= %d", ymax);
+		exit(-1);
+	}
+
+	if (X == 1) {
+		t = 0;
+	}
+	else {
+		t = 0;
+		t2 = 1;
+		for (i = 0; i < X - 1; ++i) {
+			t2 *= 2;
+			t += t2;
+		}
+	}
+
+	t += Y;
+	t2 = wt->nodelength[t];
+	for (i = 0; i < N; ++i) {
+		coeffs[i] = wt->output[t2+i];
+	}
+
+}
+
+void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N) {
+	int ymax, i;
+	int np,citer;
+	int flag;
+
+	if (X <= 0 || X > wt->J) {
+		printf("X co-ordinate must be >= 1 and <= %d", wt->J);
+		exit(-1);
+	}
+	ymax = 1;
+	for (i = 0; i < X; ++i) {
+		ymax *= 2;
+	}
+
+	ymax -= 1;
+
+	if (Y < 0 || Y > ymax) {
+		printf("Y co-ordinate must be >= 0 and <= %d", ymax);
+		exit(-1);
+	}
+
+	np = 0;
+	citer = 0;
+
+	for (i = wt->J; i > X; --i) {
+		np += wt->numnodeslevel[i];
+		citer += wt->numnodeslevel[i] * wt->coeflength[i];
+	}
+
+	i = 0;
+	flag = 0;
+	for (i = 0; i < wt->numnodeslevel[X]; ++i) {
+		if (wt->nodeindex[2 * np + 1] == Y) {
+			flag = 1;
+			break;
+		}
+		np++;
+		citer += wt->coeflength[X];
+	}
+
+	if (flag == 0) {
+		printf("The Node is Not Part Of The Best Basis Tree Use wpt_summary function to list available nodes \n");
+		exit(-1);
+	}
+
+	for (i = 0; i < N; ++i) {
+		coeffs[i] = wt->output[citer + i];
+	}
+
+}
+
+int getCWTScaleLength(int N) {
+	int J;
+	double temp,dj;
+
+	dj = 0.4875;
+
+	temp = (log((double)N / 2.0) / log(2.0)) / dj;
+	J = (int)temp;
+
+	return J;
+}
+
+void setCWTScales(cwt_object wt, double s0, double dj,char *type,int power) {
+	int i;
+	//s0*pow(2.0, (double)(j - 1)*dj);
+	if (!strcmp(wt->type, "pow") || !strcmp(wt->type, "power")) {
+		for (i = 0; i < wt->J; ++i) {
+			wt->scale[i] = s0*pow((double) power, (double)(i)*dj);
+		}
+		wt->sflag = 1;
+		
+	}
+	else if (!strcmp(wt->type, "lin") || !strcmp(wt->type, "linear")) {
+		for (i = 0; i < wt->J; ++i) {
+			wt->scale[i] = s0 + (double)i * dj;
+		}
+		wt->sflag = 1;
+	}
+	else {
+		printf("\n Type accepts only two values : pow and lin\n");
+		exit(-1);
+	}
+	wt->s0 = s0;
+	wt->dj = dj;
+}
+
+void setCWTScaleVector(cwt_object wt, double *scale, int J,double s0,double dj) {
+	int i;
+
+	if (J != wt->J) {
+		printf("\n CWT object is only valid for %d scales\n", wt->J);
+		exit(-1);
+	}
+
+	for (i = 0; i < wt->J; ++i) {
+		wt->scale[i] = scale[i];
+	}
+	wt->dj = dj;
+	wt->s0 = s0;
+	wt->sflag = 1;
+}
+
+void setCWTPadding(cwt_object wt, int pad) {
+	if (pad == 0) {
+		wt->pflag = 0;
+	}
+	else {
+		wt->pflag = 1;
+	}
+}
+
+void cwt(cwt_object wt, double *inp) {
+	int i, N, npad,nj2,j,j2;
+	N = wt->siglength;
+	if (wt->sflag == 0) {
+		for (i = 0; i < wt->J; ++i) {
+			wt->scale[i] = wt->s0*pow(2.0, (double)(i)*wt->dj);
+		}
+		wt->sflag = 1;
+	}
+
+	if (wt->pflag == 0) {
+		npad = N;
+	}
+	else {
+		npad = wt->npad;
+	}
+
+	nj2 = 2 * N * wt->J;
+	j = wt->J;
+	j2 = 2 * j;
+
+	wt->smean = 0.0;
+
+	for (i = 0; i < N; ++i) {
+		wt->smean += inp[i];
+	}
+	wt->smean /= N;
+
+	cwavelet(inp, N, wt->dt, wt->mother, wt->m, wt->s0, wt->dj, wt->J,npad,wt->params, wt->params+nj2, wt->params+nj2+j, wt->params+nj2+j2);
+
+}
+
+void icwt(cwt_object wt, double *cwtop) {
+	double psi, cdel;
+	int real,i,N,nj2;
+
+	N = wt->siglength;
+	nj2 = N * 2 * wt->J;
+
+	psi0(wt->mother, wt->m, &psi, &real);
+	cdel = cdelta(wt->mother, wt->m, psi);
+
+	//printf("\n PSI %g CDEL %g param %g mother %d \n", psi, cdel,wt->m,wt->mother);
+
+	icwavelet(wt->params, N, wt->params+nj2, wt->J, wt->dt, wt->dj, cdel, psi, cwtop);
+	
+	for(i = 0; i < N;++i) {
+		cwtop[i] += wt->smean;
+	}
+	
+}
+
+static void idwt1(wt_object wt,double *temp, double *cA_up,double *cA, int len_cA,double *cD,int len_cD,double *X_lp,double *X_hp,double *X) {
+	int len_avg, N, U,N2,i;
+
+	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+	N = 2 * len_cD;
+	U = 2;
+
+	upsamp2(cA, len_cA, U, cA_up);
+
+	per_ext(cA_up, 2 * len_cA, len_avg / 2, temp);
+
+	N2 = 2 * len_cA + len_avg;
+
+	if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+		wt->cobj = conv_init(N2, len_avg);
+		wt->cfftset = 1;
+	}
+	else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) {
+		printf("Decomposition Filters must have the same length.");
+		exit(-1);
+	}
+
+	wconv(wt, temp, N2, wt->wave->lpr, len_avg, X_lp);
+
+	upsamp2(cD, len_cD, U, cA_up);
+
+	per_ext(cA_up, 2 * len_cD, len_avg / 2, temp);
+
+	N2 = 2 * len_cD + len_avg;
+
+	wconv(wt, temp, N2, wt->wave->hpr, len_avg, X_hp);
+
+
+	for (i = len_avg - 1; i < N + len_avg - 1; ++i) {
+		X[i - len_avg + 1] = X_lp[i] + X_hp[i];
+	}
+
+	if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+		free_conv(wt->cobj);
+		wt->cfftset = 0;
+	}
+
+}
+
+static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_cD,  double *X) {
+	int len_avg,i,l,m,n,t,l2;
+
+	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+	l2 = len_avg / 2;
+	m = -2;
+	n = -1;
+
+	for (i = 0; i < len_cA + l2 - 1; ++i) {
+		m += 2;
+		n += 2;
+		X[m] = 0.0;
+		X[n] = 0.0;
+		for (l = 0; l < l2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
+				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
+			}
+			else if ((i - l) >= len_cA && (i-l) < len_cA + len_avg - 1) {
+				X[m] += wt->wave->lpr[t] * cA[i - l - len_cA] + wt->wave->hpr[t] * cD[i - l - len_cA];
+				X[n] += wt->wave->lpr[t + 1] * cA[i - l - len_cA] + wt->wave->hpr[t + 1] * cD[i - l - len_cA];
+			}
+			else if ((i - l) < 0 && (i-l) > -l2) {
+				X[m] += wt->wave->lpr[t] * cA[len_cA + i - l] + wt->wave->hpr[t] * cD[len_cA + i - l];
+				X[n] += wt->wave->lpr[t + 1] * cA[len_cA + i - l] + wt->wave->hpr[t + 1] * cD[len_cA + i - l];
+			}
+		}
+	}
+}
+
+static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) {
+	int len_avg, i, l, m, n, t, v;
+
+	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+	m = -2;
+	n = -1;
+
+	for (v = 0; v < len_cA; ++v) {
+		i = v;
+		m += 2;
+		n += 2;
+		X[m] = 0.0;
+		X[n] = 0.0;
+		for (l = 0; l < len_avg / 2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
+				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
+			}
+		}
+	}
+}
+
+
+void idwt(wt_object wt, double *dwtop) {
+	int J,U,i,lf,N,N2,iter,k;
+	int app_len, det_len;
+	double *cA_up, *X_lp, *X_hp,*out,*temp;
+
+	J = wt->J;
+	U = 2;
+	app_len = wt->length[0];
+	out = (double*)malloc(sizeof(double)* (wt->siglength + 1));
+	if (!strcmp(wt->ext, "per") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+		app_len = wt->length[0];
+		det_len = wt->length[1];
+		N = 2 * wt->length[J];
+		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+
+		cA_up = (double*)malloc(sizeof(double)* N);
+		temp = (double*)malloc(sizeof(double)* (N + lf));
+		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+		X_hp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+		iter = app_len;
+
+		for (i = 0; i < app_len; ++i) {
+			out[i] = wt->output[i];
+		}
+
+		for (i = 0; i < J; ++i) {
+
+			idwt1(wt,temp,cA_up,out,det_len,wt->output+iter,det_len,X_lp,X_hp,out);
+			/*
+			idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp);
+			for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
+				out[k - lf/2 + 1] = X_lp[k];
+			}
+			*/
+			iter += det_len;
+			det_len = wt->length[i + 2];
+		}
+		free(cA_up);
+		free(X_lp);
+		free(X_hp);
+		free(temp);
+
+	}
+	else if (!strcmp(wt->ext, "per") && !strcmp(wt->cmethod, "direct")) {
+		app_len = wt->length[0];
+		det_len = wt->length[1];
+		N = 2 * wt->length[J];
+		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+
+		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+		iter = app_len;
+
+		for (i = 0; i < app_len; ++i) {
+			out[i] = wt->output[i];
+		}
+
+		for (i = 0; i < J; ++i) {
+
+			//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
+
+			idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp);
+			for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
+			out[k - lf/2 + 1] = X_lp[k];
+			}
+
+			iter += det_len;
+			det_len = wt->length[i + 2];
+		}
+
+		free(X_lp);
+
+	}
+	else if (!strcmp(wt->ext, "sym") && !strcmp(wt->cmethod, "direct")) {
+		app_len = wt->length[0];
+		det_len = wt->length[1];
+		N = 2 * wt->length[J] - 1;
+		lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+
+		X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
+		iter = app_len;
+
+		for (i = 0; i < app_len; ++i) {
+			out[i] = wt->output[i];
+		}
+
+		for (i = 0; i < J; ++i) {
+
+			//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
+
+			idwt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp);
+			for (k = lf-2; k < 2 * det_len; ++k) {
+				out[k - lf + 2] = X_lp[k];
+			}
+
+			iter += det_len;
+			det_len = wt->length[i + 2];
+		}
+
+		free(X_lp);
+
+	}
+	else if (!strcmp(wt->ext, "sym") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+		lf = wt->wave->lpd_len;// lpd and hpd have the same length
+
+		N = 2 * wt->length[J] - 1;
+		cA_up = (double*)malloc(sizeof(double)* N);
+		X_lp = (double*)malloc(sizeof(double)* (N + lf - 1));
+		X_hp = (double*)malloc(sizeof(double)* (N + lf - 1));
+
+		for (i = 0; i < app_len; ++i) {
+			out[i] = wt->output[i];
+		}
+
+		iter = app_len;
+
+		for (i = 0; i < J; ++i) {
+			det_len = wt->length[i + 1];
+			upsamp(out, det_len, U, cA_up);
+			N2 = 2 * wt->length[i + 1] - 1;
+
+			if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+				wt->cobj = conv_init(N2, lf);
+				wt->cfftset = 1;
+			}
+			else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) {
+				printf("Decomposition Filters must have the same length.");
+				exit(-1);
+			}
+
+			wconv(wt, cA_up, N2, wt->wave->lpr, lf, X_lp);
+
+			upsamp(wt->output + iter, det_len, U, cA_up);
+
+			wconv(wt, cA_up, N2, wt->wave->hpr, lf, X_hp);
+
+			for (k = lf - 2; k < N2 + 1; ++k) {
+				out[k - lf + 2] = X_lp[k] + X_hp[k];
+			}
+			iter += det_len;
+			if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+				free_conv(wt->cobj);
+				wt->cfftset = 0;
+			}
+		}
+
+		free(cA_up);
+		free(X_lp);
+		free(X_hp);
+	}
+	else {
+		printf("Signal extension can be either per or sym");
+		exit(-1);
+	}
+
+	for (i = 0; i < wt->siglength; ++i) {
+		dwtop[i] = out[i];
+	}
+
+
+	free(out);
+
+}
+
+static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) {
+	int len_avg, i, l, m, n, t, l2;
+
+	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+	l2 = len_avg / 2;
+	m = -2;
+	n = -1;
+
+	for (i = 0; i < len_cA + l2 - 1; ++i) {
+		m += 2;
+		n += 2;
+		X[m] = 0.0;
+		X[n] = 0.0;
+		for (l = 0; l < l2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
+				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
+			}
+			else if ((i - l) >= len_cA && (i - l) < len_cA + len_avg - 1) {
+				X[m] += wt->wave->lpr[t] * cA[i - l - len_cA] + wt->wave->hpr[t] * cD[i - l - len_cA];
+				X[n] += wt->wave->lpr[t + 1] * cA[i - l - len_cA] + wt->wave->hpr[t + 1] * cD[i - l - len_cA];
+			}
+			else if ((i - l) < 0 && (i - l) > -l2) {
+				X[m] += wt->wave->lpr[t] * cA[len_cA + i - l] + wt->wave->hpr[t] * cD[len_cA + i - l];
+				X[n] += wt->wave->lpr[t + 1] * cA[len_cA + i - l] + wt->wave->hpr[t + 1] * cD[len_cA + i - l];
+			}
+		}
+	}
+}
+
+static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) {
+	int len_avg, i, l, m, n, t, v;
+
+	len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+	m = -2;
+	n = -1;
+
+	for (v = 0; v < len_cA; ++v) {
+		i = v;
+		m += 2;
+		n += 2;
+		X[m] = 0.0;
+		X[n] = 0.0;
+		for (l = 0; l < len_avg / 2; ++l) {
+			t = 2 * l;
+			if ((i - l) >= 0 && (i - l) < len_cA) {
+				X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
+				X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
+			}
+		}
+	}
+}
+
+void idwpt(wpt_object wt, double *dwtop) {
+	int J, U, i, lf, N, k,p,l;
+	int app_len, det_len, index, n1, llb, index2, index3, index4,indexp,xlen;
+	double *X_lp, *X,  *out, *out2;
+	int *prep,*ptemp;
+
+	J = wt->J;
+	U = 2;
+	app_len = wt->length[0];
+	p = ipow2(J);
+	lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
+	xlen = p * (app_len + 2 * lf);
+
+	X_lp = (double*)malloc(sizeof(double)* 2 * (wt->length[J] + lf));
+	X = (double*)malloc(sizeof(double)* xlen);
+	out = (double*)malloc(sizeof(double)* wt->length[J]);
+	out2 = (double*)malloc(sizeof(double)* wt->length[J]);
+	prep = (int*)malloc(sizeof(int)* p);
+	ptemp = (int*)malloc(sizeof(int)* p);
+	n1 = 1;
+	llb = 1;
+	index2 = xlen / p;
+	indexp = 0;
 	if (wt->basisvector[0] == 1) {
 		for (i = 0; i < wt->siglength; ++i) {
 			dwtop[i] = wt->output[i];
 		}
-		
-	} else {
+
+	}
+	else {
 		for (i = 0; i < J; ++i) {
 			llb *= 2;
 			n1 += llb;
 		}
-	
+
 		for (i = 0; i < xlen; ++i) {
 			X[i] = 0.0;
 		}
-	
+
 		for (i = 0; i < llb; ++i) {
-			prep[i] = (int) wt->basisvector[n1 - llb + i];
+			prep[i] = (int)wt->basisvector[n1 - llb + i];
 			ptemp[i] = 0;
 		}
-	
+
 		if (!strcmp(wt->ext, "per")) {
 			app_len = wt->length[0];
 			det_len = wt->length[1];
 			index = 0;
-	
-	
+
+
 			for (i = 0; i < J; ++i) {
 				p = ipow2(J - i - 1);
 				det_len = wt->length[i + 1];
@@ -1599,15 +1802,15 @@ void idwpt(wpt_object wt, double *dwtop) {
 				n1 -= llb;
 				for (l = 0; l < llb; ++l) {
 					if (ptemp[l] != 2) {
-						prep[l] = (int) wt->basisvector[n1 + l];
+						prep[l] = (int)wt->basisvector[n1 + l];
 					}
 					else {
 						prep[l] = ptemp[l];
 					}
 					ptemp[l] = 0;
 				}
-	
-	
+
+
 				for (l = 0; l < p; ++l) {
 					if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) {
 						for (k = 0; k < det_len; ++k) {
@@ -1669,34 +1872,34 @@ void idwpt(wpt_object wt, double *dwtop) {
 						index3 += index2;
 						index4 += 2 * indexp;
 					}
-	
+
 				}
-	
-	
+
+
 				/*
 				idwt_per(wt, out, det_len, wt->output + iter, det_len, X_lp);
 				for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
-					out[k - lf / 2 + 1] = X_lp[k];
+				out[k - lf / 2 + 1] = X_lp[k];
 				}
-	
+
 				iter += det_len;
 				det_len = wt->length[i + 2];
 				*/
 				llb /= 2;
 				indexp = index2;
 			}
-	
+
 			//free(X_lp);
-	
+
 		}
 		else if (!strcmp(wt->ext, "sym")) {
 			app_len = wt->length[0];
 			det_len = wt->length[1];
 			N = 2 * wt->length[J] - 1;
-	
+
 			//X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
 			index = 0;
-	
+
 			for (i = 0; i < J; ++i) {
 				p = ipow2(J - i - 1);
 				det_len = wt->length[i + 1];
@@ -1707,15 +1910,15 @@ void idwpt(wpt_object wt, double *dwtop) {
 				n1 -= llb;
 				for (l = 0; l < llb; ++l) {
 					if (ptemp[l] != 2) {
-						prep[l] = (int) wt->basisvector[n1 + l];
+						prep[l] = (int)wt->basisvector[n1 + l];
 					}
 					else {
 						prep[l] = ptemp[l];
 					}
 					ptemp[l] = 0;
 				}
-	
-	
+
+
 				for (l = 0; l < p; ++l) {
 					if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) {
 						for (k = 0; k < det_len; ++k) {
@@ -1777,779 +1980,807 @@ void idwpt(wpt_object wt, double *dwtop) {
 						index3 += index2;
 						index4 += 2 * indexp;
 					}
-	
+
 				}
-	
+
 				//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
 				/*
 				idwpt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp);
 				for (k = lf - 2; k < 2 * det_len; ++k) {
-					out[k - lf + 2] = X_lp[k];
+				out[k - lf + 2] = X_lp[k];
 				}
-	
+
 				iter += det_len;
 				det_len = wt->length[i + 2];
 				*/
 				llb /= 2;
 				indexp = index2;
 			}
-	
+
 			//free(X_lp);
-	
+
 		}
 		else {
 			printf("Signal extension can be either per or sym");
 			exit(-1);
 		}
-	
+
 		for (i = 0; i < wt->siglength; ++i) {
 			//printf("%g ", X[i]);
 			dwtop[i] = X[i];
 		}
-	
-	}
 
-
-	free(out);
-	free(X_lp);
-	free(X);
-	free(out2);
-	free(prep);
-	free(ptemp);
-}
-
-
-static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int l, l2, isodd, i, t, len_avg,j;
-
-	len_avg = M * wt->wave->lpd_len;
-	l2 = len_avg / 2;
-	isodd = N % 2;
-
-	for (i = 0; i < len_cA; ++i) {
-		t = i + l2;
-		cA[i] = 0.0;
-		cD[i] = 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) {
-				cA[i] += wt->wave->lpd[l] * inp[t - j];
-				cD[i] += wt->wave->hpd[l] * inp[t - j];
-			}
-			else if ((t - j) < l2 && (t - j) >= 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - j];
-				cD[i] += wt->wave->hpd[l] * inp[t - j];
-			}
-			else if ((t - j) < 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - j + N];
-				cD[i] += wt->wave->hpd[l] * inp[t - j + N];
-			}
-			else if ((t - j) >= N && isodd == 0) {
-				cA[i] += wt->wave->lpd[l] * inp[t - j - N];
-				cD[i] += wt->wave->hpd[l] * inp[t - j - N];
-			}
-			else if ((t - j) >= N && isodd == 1) {
-				if (t - l != N) {
-					cA[i] += wt->wave->lpd[l] * inp[t - j - (N + 1)];
-					cD[i] += wt->wave->hpd[l] * inp[t - j - (N + 1)];
-				}
-				else {
-					cA[i] += wt->wave->lpd[l] * inp[N - 1];
-					cD[i] += wt->wave->hpd[l] * inp[N - 1];
-				}
-			}
-
-		}
-	}
-
-
-}
-
-static void swt_fft(wt_object wt, double *inp) {
-	int i, J, temp_len, iter, M, N, len_filt;
-	int lenacc;
-	double *low_pass, *high_pass,*sig,*cA,*cD;
-
-	temp_len = wt->siglength;
-	J = wt->J;
-	wt->length[0] = wt->length[J] = temp_len;
-	wt->outlength = wt->length[J+1] = (J + 1) * temp_len;
-	M = 1;
-	for (iter = 1; iter < J; ++iter) {
-		M = 2 * M;
-		wt->length[iter] = temp_len;
-	}
-
-	len_filt = wt->wave->filtlength;
-
-	low_pass = (double*)malloc(sizeof(double)* M * len_filt);
-	high_pass = (double*)malloc(sizeof(double)* M * len_filt);
-	sig = (double*)malloc(sizeof(double)* (M * len_filt + temp_len + (temp_len%2)));
-	cA = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1);
-	cD = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1);
-
-	M = 1;
-
-	for (i = 0; i < temp_len; ++i) {
-		wt->params[i] = inp[i];
-	}
-
-	lenacc = wt->outlength;
-
-	for (iter = 0; iter < J; ++iter) {
-		lenacc -= temp_len;
-		if (iter > 0) {
-			M = 2 * M;
-			N = M * len_filt;
-			upsamp2(wt->wave->lpd, wt->wave->lpd_len, M, low_pass);
-			upsamp2(wt->wave->hpd, wt->wave->hpd_len, M, high_pass);
-		}
-		else {
-			N = len_filt;
-			for (i = 0; i < N; ++i) {
-				low_pass[i] = wt->wave->lpd[i];
-				high_pass[i] = wt->wave->hpd[i];
-			}
-		}
-
-		//swt_per(wt,M, wt->params, temp_len, cA, temp_len, cD,temp_len);
-
-		per_ext(wt->params, temp_len, N / 2, sig);
-
-		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-			wt->cobj = conv_init(N + temp_len + (temp_len % 2), N);
-			wt->cfftset = 1;
-		}
-		else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
-			printf("Decomposition Filters must have the same length.");
-			exit(-1);
-		}
-
-
-		wconv(wt, sig, N + temp_len + (temp_len % 2), low_pass, N, cA);
-
-		wconv(wt, sig, N + temp_len + (temp_len % 2), high_pass, N, cD);
-
-		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-			free_conv(wt->cobj);
-			wt->cfftset = 0;
-		}
-
-		for (i = 0; i < temp_len; ++i) {
-			wt->params[i] = cA[N + i];
-			wt->params[lenacc + i] = cD[N + i];
-		}
-
-
-	}
-
-	free(low_pass);
-	free(high_pass);
-	free(sig);
-	free(cA);
-	free(cD);
-}
-
-static void swt_direct(wt_object wt, double *inp) {
-	int i, J, temp_len, iter, M, N;
-	int lenacc, len_filt;
-	double  *cA, *cD;
-
-	temp_len = wt->siglength;
-	J = wt->J;
-	wt->length[0] = wt->length[J] = temp_len;
-	wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
-	len_filt = wt->wave->filtlength;
-	M = 1;
-	for (iter = 1; iter < J; ++iter) {
-		M = 2 * M;
-		wt->length[iter] = temp_len;
-	}
-
-
-	cA = (double*)malloc(sizeof(double)* temp_len);
-	cD = (double*)malloc(sizeof(double)* temp_len);
-
-	M = 1;
-
-	for (i = 0; i < temp_len; ++i) {
-		wt->params[i] = inp[i];
-	}
-
-	lenacc = wt->outlength;
-
-	for (iter = 0; iter < J; ++iter) {
-		lenacc -= temp_len;
-		if (iter > 0) {
-			M = 2 * M;
-			N = M * len_filt;
-		}
-		else {
-			N = len_filt;
-		}
-
-		swt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len);
-
-
-		for (i = 0; i < temp_len; ++i) {
-			wt->params[i] = cA[i];
-			wt->params[lenacc + i] = cD[i];
-		}
-
-	}
-
-	free(cA);
-	free(cD);
-
-}
-
-
-void swt(wt_object wt, double *inp) {
-	if (!strcmp(wt->method, "swt") && !strcmp(wt->cmethod, "direct") ) {
-		swt_direct(wt,inp);
-	}
-	else if (!strcmp(wt->method, "swt") &&  (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-		swt_fft(wt, inp);
-	}
-	else {
-		printf("SWT Only accepts two methods - direct and fft");
-		exit(-1);
-	}
-}
-
-void iswt(wt_object wt, double *swtop) {
-	int N, lf, iter,i,J,index,value,count,len;
-	int index_shift,len0,U,N1,index2;
-	double *appx1, *det1,*appx_sig,*det_sig,*cL0,*cH0,*tempx,*oup00L,*oup00H,*oup00,*oup01,*appx2,*det2;
-
-	N = wt->siglength;
-	J = wt->J;
-	U = 2;
-	lf = wt->wave->lpr_len;
-
-	appx_sig = (double*)malloc(sizeof(double)* N);
-	det_sig = (double*)malloc(sizeof(double)* N);
-	appx1 = (double*)malloc(sizeof(double)* N);
-	det1 = (double*)malloc(sizeof(double)* N);
-	appx2 = (double*)malloc(sizeof(double)* N);
-	det2 = (double*)malloc(sizeof(double)* N);
-	tempx = (double*)malloc(sizeof(double)* N);
-	cL0 = (double*)malloc(sizeof(double)* (N + (N%2) + lf));
-	cH0 = (double*)malloc(sizeof(double)* (N + (N % 2) + lf));
-	oup00L = (double*)malloc(sizeof(double)* (N + 2 * lf));
-	oup00H = (double*)malloc(sizeof(double)* (N + 2 * lf));
-	oup00 = (double*)malloc(sizeof(double)* N);
-	oup01 = (double*)malloc(sizeof(double)* N);
-
-
-
-	for (iter = 0; iter < J; ++iter) {
-		for (i = 0; i < N; ++i) {
-			swtop[i] = 0.0;
-		}
-		if (iter == 0) {
-			for (i = 0; i < N; ++i) {
-				appx_sig[i] = wt->output[i];
-				det_sig[i] = wt->output[N + i];
-			}
-		}
-		else {
-			for (i = 0; i < N; ++i) {
-				det_sig[i] = wt->output[(iter + 1) * N + i];
-			}
-		}
-
-		value = (int)pow(2.0, (double) (J - 1 - iter));
-
-		for (count = 0; count < value; count++) {
-			len = 0;
-			for (index = count; index < N; index += value) {
-				appx1[len] = appx_sig[index];
-				det1[len] = det_sig[index];
-				len++;
-			}
-
-
-			//SHIFT 0
-			len0 = 0;
-
-			for (index_shift = 0; index_shift < len; index_shift += 2) {
-				appx2[len0] = appx1[index_shift];
-				det2[len0] = det1[index_shift];
-				len0++;
-			}
-			upsamp2(appx2, len0, U, tempx);
-			per_ext(tempx, 2 * len0, lf / 2, cL0);
-
-			upsamp2(det2, len0, U, tempx);
-			per_ext(tempx, 2 * len0, lf / 2, cH0);
-
-			N1 = 2 * len0 + lf;
-
-			if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
-				wt->cobj = conv_init(N1, lf);
-				wt->cfftset = 1;
-			}
-			else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
-				printf("Decomposition Filters must have the same length.");
-				exit(-1);
-			}
-
-			wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L);
-
-			wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H);
-
-			for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
-				oup00[i - lf + 1] = oup00L[i] + oup00H[i];
-			}
-
-			//SHIFT 1
-
-			len0 = 0;
-
-			for (index_shift = 1; index_shift < len; index_shift += 2) {
-				appx2[len0] = appx1[index_shift];
-				det2[len0] = det1[index_shift];
-				len0++;
-			}
-
-			upsamp2(appx2, len0, U, tempx);
-			per_ext(tempx, 2 * len0, lf / 2, cL0);
-
-			upsamp2(det2, len0, U, tempx);
-			per_ext(tempx, 2 * len0, lf / 2, cH0);
-
-			N1 = 2 * len0 + lf;
-
-			wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L);
-
-			wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H);
-
-			for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
-				oup01[i - lf + 1] = oup00L[i] + oup00H[i];
-			}
-
-			circshift(oup01, 2*len0, -1);
-
-			index2 = 0;
-
-			for (index = count; index < N; index += value) {
-				swtop[index] = (oup00[index2] + oup01[index2]) / 2.0;
-				index2++;
-			}
-
-		}
-		for (i = 0; i < N; ++i) {
-			appx_sig[i] = swtop[i];
-		}
-
-	}
-
-
-
-	free(appx_sig);
-	free(det_sig);
-	free(appx1);
-	free(det1);
-	free(tempx);
-	free(cL0);
-	free(cH0);
-	free(oup00L);
-	free(oup00H);
-	free(oup00);
-	free(oup01);
-	free(appx2);
-	free(det2);
-}
-
-static void modwt_per(wt_object wt, int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
-	int l, i, t, len_avg;
-	double s;
-	double *filt;
-	len_avg = wt->wave->lpd_len;
-
-	filt = (double*)malloc(sizeof(double)* 2 * len_avg);
-	s = sqrt(2.0);
-	for (i = 0; i < len_avg; ++i) {
-		filt[i] = wt->wave->lpd[i] / s;
-		filt[len_avg + i] = wt->wave->hpd[i] / s;
-	}
-
-	for (i = 0; i < len_cA; ++i) {
-		t = i;
-		cA[i] = filt[0] * inp[t];
-		cD[i] = filt[len_avg] * inp[t];
-		for (l = 1; l < len_avg; l++) {
-			t -= M;
-			while (t >= len_cA) {
-				t -= len_cA;
-			}
-			while (t < 0) {
-				t += len_cA;
-			}
-
-			cA[i] += filt[l] * inp[t];
-			cD[i] += filt[len_avg + l] * inp[t];
-
-		}
-	}
-	free(filt);
-}
-
-void modwt(wt_object wt, double *inp) {
-	int i, J, temp_len, iter, M, N;
-	int lenacc, len_filt;
-	double  *cA, *cD;
-
-	temp_len = wt->siglength;
-	J = wt->J;
-	wt->length[0] = wt->length[J] = temp_len;
-	wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
-	len_filt = wt->wave->filtlength;
-	M = 1;
-	for (iter = 1; iter < J; ++iter) {
-		M = 2 * M;
-		wt->length[iter] = temp_len;
-	}
-
-
-	cA = (double*)malloc(sizeof(double)* temp_len);
-	cD = (double*)malloc(sizeof(double)* temp_len);
-
-	M = 1;
-
-	for (i = 0; i < temp_len; ++i) {
-		wt->params[i] = inp[i];
-	}
-
-	lenacc = wt->outlength;
-
-	for (iter = 0; iter < J; ++iter) {
-		lenacc -= temp_len;
-		if (iter > 0) {
-			M = 2 * M;
-			N = M * len_filt;
-		}
-		else {
-			N = len_filt;
-		}
-
-		modwt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len);
-
-
-		for (i = 0; i < temp_len; ++i) {
-			wt->params[i] = cA[i];
-			wt->params[lenacc + i] = cD[i];
-		}
-
-	}
-
-	free(cA);
-	free(cD);
-
-}
-
-static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, int len_cD, double *X) {
-	int len_avg, i, l, t;
-	double s;
-	double *filt;
-	len_avg = wt->wave->lpd_len;
-
-	filt = (double*)malloc(sizeof(double)* 2 * len_avg);
-	s = sqrt(2.0);
-	for (i = 0; i < len_avg; ++i) {
-		filt[i] = wt->wave->lpd[i] / s;
-		filt[len_avg + i] = wt->wave->hpd[i] / s;
-	}
-
-
-	for (i = 0; i < len_cA; ++i) {
-		t = i;
-		X[i] = (filt[0] * cA[t]) + (filt[len_avg] * cD[t]);
-		for (l = 1; l < len_avg; l++) {
-			t += M;
-			while (t >= len_cA) {
-				t -= len_cA;
-			}
-			while (t < 0) {
-				t += len_cA;
-			}
-
-			X[i] += (filt[l] * cA[t]) + (filt[len_avg + l] * cD[t]);
-
-		}
-	}
-	free(filt);
-}
-
-void imodwt(wt_object wt, double *dwtop) {
-	int N, lf, iter, i, J, j, U;
-	int lenacc,M;
-	double *X;
-
-	N = wt->siglength;
-	J = wt->J;
-	U = 2;
-	lf = wt->wave->lpr_len;
-	lenacc = N;
-	M = (int)pow(2.0, (double)J - 1.0);
-	//M = 1;
-	X = (double*)malloc(sizeof(double)* N);
-
-	for (i = 0; i < N; ++i) {
-		dwtop[i] = wt->output[i];
-	}
-
-	for (iter = 0; iter < J; ++iter) {
-		if (iter > 0) {
-			M = M / 2;
-		}
-		imodwt_per(wt, M, dwtop, N, wt->params + lenacc, N, X);
-		/*
-		for (j = lf - 1; j < N; ++j) {
-			dwtop[j - lf + 1] = X[j];
-		}
-		for (j = 0; j < lf - 1; ++j) {
-			dwtop[N - lf + 1 + j] = X[j];
-		}
-		*/
-		for (j = 0; j < N; ++j) {
-			dwtop[j] = X[j];
-		}
-
-		lenacc += N;
-	}
-	free(X);
-}
-
-void setDWTExtension(wt_object wt, char *extension) {
-	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);
-	}
-}
-
-void setWTREEExtension(wtree_object wt, char *extension) {
-	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);
-	}
-}
-
-void setDWPTExtension(wpt_object wt, char *extension) {
-	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);
-	}
-}
-
-void setDWPTEntropy(wpt_object wt, char *entropy, double eparam) {
-	if (!strcmp(entropy, "shannon")) {
-		strcpy(wt->entropy, "shannon");
-	}
-	else if (!strcmp(entropy, "threshold")) {
-		strcpy(wt->entropy, "threshold");
-		wt ->eparam = eparam;
-	}
-	else if (!strcmp(entropy, "norm")) {
-		strcpy(wt->entropy, "norm");
-		wt->eparam = eparam;
-	}
-	else if (!strcmp(entropy, "logenergy") || !strcmp(entropy, "log energy") || !strcmp(entropy, "energy")) {
-		strcpy(wt->entropy, "logenergy");
-	}
-	else {
-		printf("Entropy should be one of shannon, threshold, norm or logenergy");
-		exit(-1);
-	}
-}
-
-void setWTConv(wt_object wt, char *cmethod) {
-	if (!strcmp(cmethod, "fft") || !strcmp(cmethod, "FFT")) {
-		strcpy(wt->cmethod, "fft");
-	}
-	else if (!strcmp(cmethod, "direct")) {
-		strcpy(wt->cmethod, "direct");
-	}
-	else {
-		printf("Convolution Only accepts two methods - direct and fft");
-		exit(-1);
-	}
-}
-
-void wave_summary(wave_object obj) {
-	int i,N;
-	N = obj->filtlength;
-	printf("\n");
-	printf("Wavelet Name : %s \n",obj->wname);
-	printf("\n");
-	printf("Wavelet Filters \n\n");
-	printf("lpd : [");
-	for (i = 0; i < N-1; ++i) {
-		printf("%g,", obj->lpd[i]);
-	}
-	printf("%g", obj->lpd[N-1]);
-	printf("] \n\n");
-	printf("hpd : [");
-	for (i = 0; i < N-1; ++i) {
-		printf("%g,", obj->hpd[i]);
-	}
-	printf("%g", obj->hpd[N - 1]);
-	printf("] \n\n");
-	printf("lpr : [");
-	for (i = 0; i < N-1; ++i) {
-		printf("%g,", obj->lpr[i]);
-	}
-	printf("%g", obj->lpr[N - 1]);
-	printf("] \n\n");
-	printf("hpr : [");
-	for (i = 0; i < N-1; ++i) {
-		printf("%g,", obj->hpr[i]);
-	}
-	printf("%g", obj->hpr[N - 1]);
-	printf("] \n\n");
-}
-
-void wt_summary(wt_object wt) {
-	int i;
-	int J,t;
-	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("Convolutional Method : %s \n", wt->cmethod);
-	printf("\n");
-	printf("Number of Decomposition Levels %d \n", wt->J);
-	printf("\n");
-	printf("Length of Input Signal %d \n", wt->siglength);
-	printf("\n");
-	printf("Length of WT Output Vector %d \n", wt->outlength);
-	printf("\n");
-	printf("Wavelet Coefficients are contained in vector : %s \n", "output");
-	printf("\n");
-	printf("Approximation Coefficients \n");
-	printf("Level %d Access : output[%d] Length : %d \n", 1, 0, wt->length[0]);
-	printf("\n");
-	printf("Detail Coefficients \n");
-	t = wt->length[0];
-	for (i = 0; i < J; ++i) {
-		printf("Level %d Access : output[%d] Length : %d \n", i + 1,t,wt->length[i+1]);
-		t += wt->length[i+1];
-	}
-	printf("\n");
-
-}
-
-void wtree_summary(wtree_object wt) {
-	int i,k,p2;
-	int J,t;
-	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("Length of Input Signal %d \n", wt->siglength);
-	printf("\n");
-	printf("Length of WT Output Vector %d \n", wt->outlength);
-	printf("\n");
-	printf("Wavelet Coefficients are contained in vector : %s \n", "output");
-	printf("\n");
-	printf("Coefficients Access \n");
-	t = 0;
-	p2 = 2;
-	for (i = 0; i < J; ++i) {
-		for (k = 0; k < p2; ++k) {
-			printf("Node %d %d Access : output[%d] Length : %d \n", i + 1, k, wt->nodelength[t], wt->length[J - i]);
-			t++;
-		}
-		p2 *= 2;
-	}
-	printf("\n");
-
-}
-
-void wpt_summary(wpt_object wt) {
-	int i, k, p2;
-	int J, it1,it2;
-	J = wt->J;
-	wave_summary(wt->wave);
-	printf("\n");
-	printf("Signal Extension : %s \n", wt->ext);
-	printf("\n");
-	printf("Entropy : %s \n", wt->entropy);
-	printf("\n");
-	printf("Number of Decomposition Levels %d \n", wt->J);
-	printf("\n");
-	printf("Number of Active Nodes %d \n", wt->nodes);
-	printf("\n");
-	printf("Length of Input Signal %d \n", wt->siglength);
-	printf("\n");
-	printf("Length of WT Output Vector %d \n", wt->outlength);
-	printf("\n");
-	printf("Wavelet Coefficients are contained in vector : %s \n", "output");
-	printf("\n");
-	printf("Coefficients Access \n");
-	it1 = 1;
-	it2 = 0;
-	for (i = 0; i < J; ++i) {
-		it1 += ipow2(i + 1);
-	}
-	for (i = J; i > 0; --i) {
-		p2 = ipow2(i);
-		it1 -= p2;
-		for (k = 0; k < p2; ++k) {
-			if (wt->basisvector[it1 + k] == 1) {
-				printf("Node %d %d Access : output[%d] Length : %d \n", i, k, it2, wt->length[J - i + 1]);
-				it2 += wt->length[J - i + 1];
-			}
-		}
-	}
-
-	printf("\n");
-
-}
-
-void wave_free(wave_object object) {
-	free(object);
-}
-
-void wt_free(wt_object object) {
-	free(object);
-}
-
-void wtree_free(wtree_object object) {
-	free(object);
-}
-
-void wpt_free(wpt_object object) {
-	free(object);
-}
+	}
+
+
+	free(out);
+	free(X_lp);
+	free(X);
+	free(out2);
+	free(prep);
+	free(ptemp);
+}
+
+
+static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int l, l2, isodd, i, t, len_avg,j;
+
+	len_avg = M * wt->wave->lpd_len;
+	l2 = len_avg / 2;
+	isodd = N % 2;
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i + l2;
+		cA[i] = 0.0;
+		cD[i] = 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) {
+				cA[i] += wt->wave->lpd[l] * inp[t - j];
+				cD[i] += wt->wave->hpd[l] * inp[t - j];
+			}
+			else if ((t - j) < l2 && (t - j) >= 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - j];
+				cD[i] += wt->wave->hpd[l] * inp[t - j];
+			}
+			else if ((t - j) < 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - j + N];
+				cD[i] += wt->wave->hpd[l] * inp[t - j + N];
+			}
+			else if ((t - j) >= N && isodd == 0) {
+				cA[i] += wt->wave->lpd[l] * inp[t - j - N];
+				cD[i] += wt->wave->hpd[l] * inp[t - j - N];
+			}
+			else if ((t - j) >= N && isodd == 1) {
+				if (t - l != N) {
+					cA[i] += wt->wave->lpd[l] * inp[t - j - (N + 1)];
+					cD[i] += wt->wave->hpd[l] * inp[t - j - (N + 1)];
+				}
+				else {
+					cA[i] += wt->wave->lpd[l] * inp[N - 1];
+					cD[i] += wt->wave->hpd[l] * inp[N - 1];
+				}
+			}
+
+		}
+	}
+
+
+}
+
+static void swt_fft(wt_object wt, double *inp) {
+	int i, J, temp_len, iter, M, N, len_filt;
+	int lenacc;
+	double *low_pass, *high_pass,*sig,*cA,*cD;
+
+	temp_len = wt->siglength;
+	J = wt->J;
+	wt->length[0] = wt->length[J] = temp_len;
+	wt->outlength = wt->length[J+1] = (J + 1) * temp_len;
+	M = 1;
+	for (iter = 1; iter < J; ++iter) {
+		M = 2 * M;
+		wt->length[iter] = temp_len;
+	}
+
+	len_filt = wt->wave->filtlength;
+
+	low_pass = (double*)malloc(sizeof(double)* M * len_filt);
+	high_pass = (double*)malloc(sizeof(double)* M * len_filt);
+	sig = (double*)malloc(sizeof(double)* (M * len_filt + temp_len + (temp_len%2)));
+	cA = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1);
+	cD = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1);
+
+	M = 1;
+
+	for (i = 0; i < temp_len; ++i) {
+		wt->params[i] = inp[i];
+	}
+
+	lenacc = wt->outlength;
+
+	for (iter = 0; iter < J; ++iter) {
+		lenacc -= temp_len;
+		if (iter > 0) {
+			M = 2 * M;
+			N = M * len_filt;
+			upsamp2(wt->wave->lpd, wt->wave->lpd_len, M, low_pass);
+			upsamp2(wt->wave->hpd, wt->wave->hpd_len, M, high_pass);
+		}
+		else {
+			N = len_filt;
+			for (i = 0; i < N; ++i) {
+				low_pass[i] = wt->wave->lpd[i];
+				high_pass[i] = wt->wave->hpd[i];
+			}
+		}
+
+		//swt_per(wt,M, wt->params, temp_len, cA, temp_len, cD,temp_len);
+
+		per_ext(wt->params, temp_len, N / 2, sig);
+
+		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+			wt->cobj = conv_init(N + temp_len + (temp_len % 2), N);
+			wt->cfftset = 1;
+		}
+		else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
+			printf("Decomposition Filters must have the same length.");
+			exit(-1);
+		}
+
+
+		wconv(wt, sig, N + temp_len + (temp_len % 2), low_pass, N, cA);
+
+		wconv(wt, sig, N + temp_len + (temp_len % 2), high_pass, N, cD);
+
+		if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+			free_conv(wt->cobj);
+			wt->cfftset = 0;
+		}
+
+		for (i = 0; i < temp_len; ++i) {
+			wt->params[i] = cA[N + i];
+			wt->params[lenacc + i] = cD[N + i];
+		}
+
+
+	}
+
+	free(low_pass);
+	free(high_pass);
+	free(sig);
+	free(cA);
+	free(cD);
+}
+
+static void swt_direct(wt_object wt, double *inp) {
+	int i, J, temp_len, iter, M, N;
+	int lenacc, len_filt;
+	double  *cA, *cD;
+
+	temp_len = wt->siglength;
+	J = wt->J;
+	wt->length[0] = wt->length[J] = temp_len;
+	wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
+	len_filt = wt->wave->filtlength;
+	M = 1;
+	for (iter = 1; iter < J; ++iter) {
+		M = 2 * M;
+		wt->length[iter] = temp_len;
+	}
+
+
+	cA = (double*)malloc(sizeof(double)* temp_len);
+	cD = (double*)malloc(sizeof(double)* temp_len);
+
+	M = 1;
+
+	for (i = 0; i < temp_len; ++i) {
+		wt->params[i] = inp[i];
+	}
+
+	lenacc = wt->outlength;
+
+	for (iter = 0; iter < J; ++iter) {
+		lenacc -= temp_len;
+		if (iter > 0) {
+			M = 2 * M;
+			N = M * len_filt;
+		}
+		else {
+			N = len_filt;
+		}
+
+		swt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len);
+
+
+		for (i = 0; i < temp_len; ++i) {
+			wt->params[i] = cA[i];
+			wt->params[lenacc + i] = cD[i];
+		}
+
+	}
+
+	free(cA);
+	free(cD);
+
+}
+
+
+void swt(wt_object wt, double *inp) {
+	if (!strcmp(wt->method, "swt") && !strcmp(wt->cmethod, "direct") ) {
+		swt_direct(wt,inp);
+	}
+	else if (!strcmp(wt->method, "swt") &&  (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+		swt_fft(wt, inp);
+	}
+	else {
+		printf("SWT Only accepts two methods - direct and fft");
+		exit(-1);
+	}
+}
+
+void iswt(wt_object wt, double *swtop) {
+	int N, lf, iter,i,J,index,value,count,len;
+	int index_shift,len0,U,N1,index2;
+	double *appx1, *det1,*appx_sig,*det_sig,*cL0,*cH0,*tempx,*oup00L,*oup00H,*oup00,*oup01,*appx2,*det2;
+
+	N = wt->siglength;
+	J = wt->J;
+	U = 2;
+	lf = wt->wave->lpr_len;
+
+	appx_sig = (double*)malloc(sizeof(double)* N);
+	det_sig = (double*)malloc(sizeof(double)* N);
+	appx1 = (double*)malloc(sizeof(double)* N);
+	det1 = (double*)malloc(sizeof(double)* N);
+	appx2 = (double*)malloc(sizeof(double)* N);
+	det2 = (double*)malloc(sizeof(double)* N);
+	tempx = (double*)malloc(sizeof(double)* N);
+	cL0 = (double*)malloc(sizeof(double)* (N + (N%2) + lf));
+	cH0 = (double*)malloc(sizeof(double)* (N + (N % 2) + lf));
+	oup00L = (double*)malloc(sizeof(double)* (N + 2 * lf));
+	oup00H = (double*)malloc(sizeof(double)* (N + 2 * lf));
+	oup00 = (double*)malloc(sizeof(double)* N);
+	oup01 = (double*)malloc(sizeof(double)* N);
+
+
+
+	for (iter = 0; iter < J; ++iter) {
+		for (i = 0; i < N; ++i) {
+			swtop[i] = 0.0;
+		}
+		if (iter == 0) {
+			for (i = 0; i < N; ++i) {
+				appx_sig[i] = wt->output[i];
+				det_sig[i] = wt->output[N + i];
+			}
+		}
+		else {
+			for (i = 0; i < N; ++i) {
+				det_sig[i] = wt->output[(iter + 1) * N + i];
+			}
+		}
+
+		value = (int)pow(2.0, (double) (J - 1 - iter));
+
+		for (count = 0; count < value; count++) {
+			len = 0;
+			for (index = count; index < N; index += value) {
+				appx1[len] = appx_sig[index];
+				det1[len] = det_sig[index];
+				len++;
+			}
+
+
+			//SHIFT 0
+			len0 = 0;
+
+			for (index_shift = 0; index_shift < len; index_shift += 2) {
+				appx2[len0] = appx1[index_shift];
+				det2[len0] = det1[index_shift];
+				len0++;
+			}
+			upsamp2(appx2, len0, U, tempx);
+			per_ext(tempx, 2 * len0, lf / 2, cL0);
+
+			upsamp2(det2, len0, U, tempx);
+			per_ext(tempx, 2 * len0, lf / 2, cH0);
+
+			N1 = 2 * len0 + lf;
+
+			if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
+				wt->cobj = conv_init(N1, lf);
+				wt->cfftset = 1;
+			}
+			else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
+				printf("Decomposition Filters must have the same length.");
+				exit(-1);
+			}
+
+			wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L);
+
+			wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H);
+
+			for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
+				oup00[i - lf + 1] = oup00L[i] + oup00H[i];
+			}
+
+			//SHIFT 1
+
+			len0 = 0;
+
+			for (index_shift = 1; index_shift < len; index_shift += 2) {
+				appx2[len0] = appx1[index_shift];
+				det2[len0] = det1[index_shift];
+				len0++;
+			}
+
+			upsamp2(appx2, len0, U, tempx);
+			per_ext(tempx, 2 * len0, lf / 2, cL0);
+
+			upsamp2(det2, len0, U, tempx);
+			per_ext(tempx, 2 * len0, lf / 2, cH0);
+
+			N1 = 2 * len0 + lf;
+
+			wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L);
+
+			wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H);
+
+			for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
+				oup01[i - lf + 1] = oup00L[i] + oup00H[i];
+			}
+
+			circshift(oup01, 2*len0, -1);
+
+			index2 = 0;
+
+			for (index = count; index < N; index += value) {
+				swtop[index] = (oup00[index2] + oup01[index2]) / 2.0;
+				index2++;
+			}
+
+		}
+		for (i = 0; i < N; ++i) {
+			appx_sig[i] = swtop[i];
+		}
+
+	}
+
+
+
+	free(appx_sig);
+	free(det_sig);
+	free(appx1);
+	free(det1);
+	free(tempx);
+	free(cL0);
+	free(cH0);
+	free(oup00L);
+	free(oup00H);
+	free(oup00);
+	free(oup01);
+	free(appx2);
+	free(det2);
+}
+
+static void modwt_per(wt_object wt, int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) {
+	int l, i, t, len_avg;
+	double s;
+	double *filt;
+	len_avg = wt->wave->lpd_len;
+
+	filt = (double*)malloc(sizeof(double)* 2 * len_avg);
+	s = sqrt(2.0);
+	for (i = 0; i < len_avg; ++i) {
+		filt[i] = wt->wave->lpd[i] / s;
+		filt[len_avg + i] = wt->wave->hpd[i] / s;
+	}
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i;
+		cA[i] = filt[0] * inp[t];
+		cD[i] = filt[len_avg] * inp[t];
+		for (l = 1; l < len_avg; l++) {
+			t -= M;
+			while (t >= len_cA) {
+				t -= len_cA;
+			}
+			while (t < 0) {
+				t += len_cA;
+			}
+
+			cA[i] += filt[l] * inp[t];
+			cD[i] += filt[len_avg + l] * inp[t];
+
+		}
+	}
+	free(filt);
+}
+
+void modwt(wt_object wt, double *inp) {
+	int i, J, temp_len, iter, M, N;
+	int lenacc, len_filt;
+	double  *cA, *cD;
+
+	temp_len = wt->siglength;
+	J = wt->J;
+	wt->length[0] = wt->length[J] = temp_len;
+	wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
+	len_filt = wt->wave->filtlength;
+	M = 1;
+	for (iter = 1; iter < J; ++iter) {
+		M = 2 * M;
+		wt->length[iter] = temp_len;
+	}
+
+
+	cA = (double*)malloc(sizeof(double)* temp_len);
+	cD = (double*)malloc(sizeof(double)* temp_len);
+
+	M = 1;
+
+	for (i = 0; i < temp_len; ++i) {
+		wt->params[i] = inp[i];
+	}
+
+	lenacc = wt->outlength;
+
+	for (iter = 0; iter < J; ++iter) {
+		lenacc -= temp_len;
+		if (iter > 0) {
+			M = 2 * M;
+			N = M * len_filt;
+		}
+		else {
+			N = len_filt;
+		}
+
+		modwt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len);
+
+
+		for (i = 0; i < temp_len; ++i) {
+			wt->params[i] = cA[i];
+			wt->params[lenacc + i] = cD[i];
+		}
+
+	}
+
+	free(cA);
+	free(cD);
+
+}
+
+static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, int len_cD, double *X) {
+	int len_avg, i, l, t;
+	double s;
+	double *filt;
+	len_avg = wt->wave->lpd_len;
+
+	filt = (double*)malloc(sizeof(double)* 2 * len_avg);
+	s = sqrt(2.0);
+	for (i = 0; i < len_avg; ++i) {
+		filt[i] = wt->wave->lpd[i] / s;
+		filt[len_avg + i] = wt->wave->hpd[i] / s;
+	}
+
+
+	for (i = 0; i < len_cA; ++i) {
+		t = i;
+		X[i] = (filt[0] * cA[t]) + (filt[len_avg] * cD[t]);
+		for (l = 1; l < len_avg; l++) {
+			t += M;
+			while (t >= len_cA) {
+				t -= len_cA;
+			}
+			while (t < 0) {
+				t += len_cA;
+			}
+
+			X[i] += (filt[l] * cA[t]) + (filt[len_avg + l] * cD[t]);
+
+		}
+	}
+	free(filt);
+}
+
+void imodwt(wt_object wt, double *dwtop) {
+	int N, lf, iter, i, J, j, U;
+	int lenacc,M;
+	double *X;
+
+	N = wt->siglength;
+	J = wt->J;
+	U = 2;
+	lf = wt->wave->lpr_len;
+	lenacc = N;
+	M = (int)pow(2.0, (double)J - 1.0);
+	//M = 1;
+	X = (double*)malloc(sizeof(double)* N);
+
+	for (i = 0; i < N; ++i) {
+		dwtop[i] = wt->output[i];
+	}
+
+	for (iter = 0; iter < J; ++iter) {
+		if (iter > 0) {
+			M = M / 2;
+		}
+		imodwt_per(wt, M, dwtop, N, wt->params + lenacc, N, X);
+		/*
+		for (j = lf - 1; j < N; ++j) {
+			dwtop[j - lf + 1] = X[j];
+		}
+		for (j = 0; j < lf - 1; ++j) {
+			dwtop[N - lf + 1 + j] = X[j];
+		}
+		*/
+		for (j = 0; j < N; ++j) {
+			dwtop[j] = X[j];
+		}
+
+		lenacc += N;
+	}
+	free(X);
+}
+
+void setDWTExtension(wt_object wt, char *extension) {
+	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);
+	}
+}
+
+void setWTREEExtension(wtree_object wt, char *extension) {
+	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);
+	}
+}
+
+void setDWPTExtension(wpt_object wt, char *extension) {
+	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);
+	}
+}
+
+void setDWPTEntropy(wpt_object wt, char *entropy, double eparam) {
+	if (!strcmp(entropy, "shannon")) {
+		strcpy(wt->entropy, "shannon");
+	}
+	else if (!strcmp(entropy, "threshold")) {
+		strcpy(wt->entropy, "threshold");
+		wt ->eparam = eparam;
+	}
+	else if (!strcmp(entropy, "norm")) {
+		strcpy(wt->entropy, "norm");
+		wt->eparam = eparam;
+	}
+	else if (!strcmp(entropy, "logenergy") || !strcmp(entropy, "log energy") || !strcmp(entropy, "energy")) {
+		strcpy(wt->entropy, "logenergy");
+	}
+	else {
+		printf("Entropy should be one of shannon, threshold, norm or logenergy");
+		exit(-1);
+	}
+}
+
+void setWTConv(wt_object wt, char *cmethod) {
+	if (!strcmp(cmethod, "fft") || !strcmp(cmethod, "FFT")) {
+		strcpy(wt->cmethod, "fft");
+	}
+	else if (!strcmp(cmethod, "direct")) {
+		strcpy(wt->cmethod, "direct");
+	}
+	else {
+		printf("Convolution Only accepts two methods - direct and fft");
+		exit(-1);
+	}
+}
+
+void wave_summary(wave_object obj) {
+	int i,N;
+	N = obj->filtlength;
+	printf("\n");
+	printf("Wavelet Name : %s \n",obj->wname);
+	printf("\n");
+	printf("Wavelet Filters \n\n");
+	printf("lpd : [");
+	for (i = 0; i < N-1; ++i) {
+		printf("%g,", obj->lpd[i]);
+	}
+	printf("%g", obj->lpd[N-1]);
+	printf("] \n\n");
+	printf("hpd : [");
+	for (i = 0; i < N-1; ++i) {
+		printf("%g,", obj->hpd[i]);
+	}
+	printf("%g", obj->hpd[N - 1]);
+	printf("] \n\n");
+	printf("lpr : [");
+	for (i = 0; i < N-1; ++i) {
+		printf("%g,", obj->lpr[i]);
+	}
+	printf("%g", obj->lpr[N - 1]);
+	printf("] \n\n");
+	printf("hpr : [");
+	for (i = 0; i < N-1; ++i) {
+		printf("%g,", obj->hpr[i]);
+	}
+	printf("%g", obj->hpr[N - 1]);
+	printf("] \n\n");
+}
+
+void wt_summary(wt_object wt) {
+	int i;
+	int J,t;
+	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("Convolutional Method : %s \n", wt->cmethod);
+	printf("\n");
+	printf("Number of Decomposition Levels %d \n", wt->J);
+	printf("\n");
+	printf("Length of Input Signal %d \n", wt->siglength);
+	printf("\n");
+	printf("Length of WT Output Vector %d \n", wt->outlength);
+	printf("\n");
+	printf("Wavelet Coefficients are contained in vector : %s \n", "output");
+	printf("\n");
+	printf("Approximation Coefficients \n");
+	printf("Level %d Access : output[%d] Length : %d \n", 1, 0, wt->length[0]);
+	printf("\n");
+	printf("Detail Coefficients \n");
+	t = wt->length[0];
+	for (i = 0; i < J; ++i) {
+		printf("Level %d Access : output[%d] Length : %d \n", i + 1,t,wt->length[i+1]);
+		t += wt->length[i+1];
+	}
+	printf("\n");
+
+}
+
+void wtree_summary(wtree_object wt) {
+	int i,k,p2;
+	int J,t;
+	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("Length of Input Signal %d \n", wt->siglength);
+	printf("\n");
+	printf("Length of WT Output Vector %d \n", wt->outlength);
+	printf("\n");
+	printf("Wavelet Coefficients are contained in vector : %s \n", "output");
+	printf("\n");
+	printf("Coefficients Access \n");
+	t = 0;
+	p2 = 2;
+	for (i = 0; i < J; ++i) {
+		for (k = 0; k < p2; ++k) {
+			printf("Node %d %d Access : output[%d] Length : %d \n", i + 1, k, wt->nodelength[t], wt->length[J - i]);
+			t++;
+		}
+		p2 *= 2;
+	}
+	printf("\n");
+
+}
+
+void wpt_summary(wpt_object wt) {
+	int i, k, p2;
+	int J, it1,it2;
+	J = wt->J;
+	wave_summary(wt->wave);
+	printf("\n");
+	printf("Signal Extension : %s \n", wt->ext);
+	printf("\n");
+	printf("Entropy : %s \n", wt->entropy);
+	printf("\n");
+	printf("Number of Decomposition Levels %d \n", wt->J);
+	printf("\n");
+	printf("Number of Active Nodes %d \n", wt->nodes);
+	printf("\n");
+	printf("Length of Input Signal %d \n", wt->siglength);
+	printf("\n");
+	printf("Length of WT Output Vector %d \n", wt->outlength);
+	printf("\n");
+	printf("Wavelet Coefficients are contained in vector : %s \n", "output");
+	printf("\n");
+	printf("Coefficients Access \n");
+	it1 = 1;
+	it2 = 0;
+	for (i = 0; i < J; ++i) {
+		it1 += ipow2(i + 1);
+	}
+	for (i = J; i > 0; --i) {
+		p2 = ipow2(i);
+		it1 -= p2;
+		for (k = 0; k < p2; ++k) {
+			if (wt->basisvector[it1 + k] == 1) {
+				printf("Node %d %d Access : output[%d] Length : %d \n", i, k, it2, wt->length[J - i + 1]);
+				it2 += wt->length[J - i + 1];
+			}
+		}
+	}
+
+	printf("\n");
+
+}
+
+void cwt_summary(cwt_object wt) {
+
+	printf("\n");
+	printf("Wavelet : %s Parameter %lf \n", wt->wave,wt->m);
+	printf("\n");
+	printf("Length of Input Signal : %d \n", wt->siglength);
+	printf("\n");
+	printf("Sampling Rate : %g \n", wt->dt);
+	printf("\n");
+	printf("Total Number of Scales : %d \n", wt->J);
+	printf("\n");
+	printf("Smallest Scale (s0) : %lf \n", wt->s0);
+	printf("\n");
+	printf("Separation Between Scales (dj) %lf \n", wt->dj);
+	printf("\n");
+	printf("Scale Type %s \n", wt->type);
+	printf("\n");
+	printf("Complex CWT Output Vector is of size %d * %d stored in Row Major format \n",wt->J,wt->siglength);
+	printf("\n");
+	printf("The ith real value can be accessed using wt->output[i].re and imaginary value by wt->output[i].im \n");
+	printf("\n");
+
+}
+
+void wave_free(wave_object object) {
+	free(object);
+}
+
+void wt_free(wt_object object) {
+	free(object);
+}
+
+void wtree_free(wtree_object object) {
+	free(object);
+}
+
+void wpt_free(wpt_object object) {
+	free(object);
+}
+
+void cwt_free(cwt_object object) {
+	free(object);
+}
diff --git a/src/wavelib.h b/src/wavelib.h
index e4e7152..86387f4 100644
--- a/src/wavelib.h
+++ b/src/wavelib.h
@@ -1,176 +1,230 @@
 /*
 Copyright (c) 2014, Rafat Hussain
-*/
-#ifndef WAVELIB_H_
-#define WAVELIB_H_
-
-#include "wtmath.h"
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#if defined(_MSC_VER)
-#pragma warning(disable : 4200)
-#pragma warning(disable : 4996)
-#endif
-
-typedef struct wave_set* wave_object;
-
-wave_object wave_init(char* wname);
-
-struct wave_set{
-	char wname[50];
-	int filtlength;// When all filters are of the same length. [Matlab uses zero-padding to make all filters of the same length]
-	int lpd_len;// Default filtlength = lpd_len = lpr_len = hpd_len = hpr_len
-	int hpd_len;
-	int lpr_len;
-	int hpr_len;
-	double *lpd;
-	double *hpd;
-	double *lpr;
-	double *hpr;
-    double params[0];
-};
-
-typedef struct wt_set* wt_object;
-
-wt_object wt_init(wave_object wave,char* method, int siglength, int J);
-
-struct wt_set{
-	wave_object wave;
-	conv_object cobj;
-	char method[10];
-	int siglength;// Length of the original signal.
-	int outlength;// Length of the output DWT vector
-	int lenlength;// Length of the Output Dimension Vector "length"
-	int J; // Number of decomposition Levels
-	int MaxIter;// Maximum Iterations J <= MaxIter
-	int even;// even = 1 if signal is of even length. even = 0 otherwise
-	char ext[10];// Type of Extension used - "per" or "sym"
-	char cmethod[10]; // Convolution Method - "direct" or "FFT"
-
-	int N; //
-	int cfftset;
-	int zpad;
-	int length[102];
-	double *output;
-    double params[0];
-};
-
-typedef struct wtree_set* wtree_object;
-
-wtree_object wtree_init(wave_object wave, int siglength, int J);
-
-struct wtree_set{
-	wave_object wave;
-	conv_object cobj;
-	char method[10];
-	int siglength;// Length of the original signal.
-	int outlength;// Length of the output DWT vector
-	int lenlength;// Length of the Output Dimension Vector "length"
-	int J; // Number of decomposition Levels
-	int MaxIter;// Maximum Iterations J <= MaxIter
-	int even;// even = 1 if signal is of even length. even = 0 otherwise
-	char ext[10];// Type of Extension used - "per" or "sym"
-
-	int N; //
-	int nodes;
-	int cfftset;
-	int zpad;
-	int length[102];
-	double *output;
-	int *nodelength;
-	int *coeflength;
-    double params[0];
-};
-
-typedef struct wpt_set* wpt_object;
-
-wpt_object wpt_init(wave_object wave, int siglength, int J);
-
-struct wpt_set{
-	wave_object wave;
-	conv_object cobj;
-	int siglength;// Length of the original signal.
-	int outlength;// Length of the output DWT vector
-	int lenlength;// Length of the Output Dimension Vector "length"
-	int J; // Number of decomposition Levels
-	int MaxIter;// Maximum Iterations J <= MaxIter
-	int even;// even = 1 if signal is of even length. even = 0 otherwise
-	char ext[10];// Type of Extension used - "per" or "sym"
-	char entropy[20];
-	double eparam;
-
-	int N; //
-	int nodes;
-	int length[102];
-	double *output;
-	double *costvalues;
-	double *basisvector;
-	int *nodeindex;
-	int *numnodeslevel;
-	int *coeflength;
-    double params[0];
-};
-
-
-void dwt(wt_object wt, double *inp);
-
-void idwt(wt_object wt, double *dwtop);
-
-void wtree(wtree_object wt, double *inp);
-
-void dwpt(wpt_object wt, double *inp);
-
-void idwpt(wpt_object wt, double *dwtop);
-
-void swt(wt_object wt, double *inp);
-
-void iswt(wt_object wt, double *swtop);
-
-void modwt(wt_object wt, double *inp);
-
-void imodwt(wt_object wt, double *dwtop);
-
-void setDWTExtension(wt_object wt, char *extension);
-
-void setWTREEExtension(wtree_object wt, char *extension);
-
-void setDWPTExtension(wpt_object wt, char *extension);
-
-void setDWPTEntropy(wpt_object wt, char *entropy, double eparam);
-
-void setWTConv(wt_object wt, char *cmethod);
-
-int getWTREENodelength(wtree_object wt, int X);
-
-void getWTREECoeffs(wtree_object wt, int X, int Y, double *coeffs, int N);
-
-int getDWPTNodelength(wpt_object wt, int X);
-
-void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N);
-
-void wave_summary(wave_object obj);
-
-void wt_summary(wt_object wt);
-
-void wtree_summary(wtree_object wt);
-
-void wpt_summary(wpt_object wt);
-
-void wave_free(wave_object object);
-
-void wt_free(wt_object object);
-
-void wtree_free(wtree_object object);
-
-void wpt_free(wpt_object object);
-
-
-#ifdef __cplusplus
-}
-#endif
-
-
-#endif /* WAVELIB_H_ */
+*/
+#ifndef WAVELIB_H_
+#define WAVELIB_H_
+
+#include "wtmath.h"
+#include "cwt.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#if defined(_MSC_VER)
+#pragma warning(disable : 4200)
+#pragma warning(disable : 4996)
+#endif
+
+#ifndef cplx_type
+#define cplx_type double
+#endif
+
+
+typedef struct cplx_t {
+	cplx_type re;
+	cplx_type im;
+} cplx_data;
+
+typedef struct wave_set* wave_object;
+
+wave_object wave_init(char* wname);
+
+struct wave_set{
+	char wname[50];
+	int filtlength;// When all filters are of the same length. [Matlab uses zero-padding to make all filters of the same length]
+	int lpd_len;// Default filtlength = lpd_len = lpr_len = hpd_len = hpr_len
+	int hpd_len;
+	int lpr_len;
+	int hpr_len;
+	double *lpd;
+	double *hpd;
+	double *lpr;
+	double *hpr;
+	double params[0];
+};
+
+typedef struct wt_set* wt_object;
+
+wt_object wt_init(wave_object wave,char* method, int siglength, int J);
+
+struct wt_set{
+	wave_object wave;
+	conv_object cobj;
+	char method[10];
+	int siglength;// Length of the original signal.
+	int outlength;// Length of the output DWT vector
+	int lenlength;// Length of the Output Dimension Vector "length"
+	int J; // Number of decomposition Levels
+	int MaxIter;// Maximum Iterations J <= MaxIter
+	int even;// even = 1 if signal is of even length. even = 0 otherwise
+	char ext[10];// Type of Extension used - "per" or "sym"
+	char cmethod[10]; // Convolution Method - "direct" or "FFT"
+
+	int N; //
+	int cfftset;
+	int zpad;
+	int length[102];
+	double *output;
+	double params[0];
+};
+
+typedef struct wtree_set* wtree_object;
+
+wtree_object wtree_init(wave_object wave, int siglength, int J);
+
+struct wtree_set{
+	wave_object wave;
+	conv_object cobj;
+	char method[10];
+	int siglength;// Length of the original signal.
+	int outlength;// Length of the output DWT vector
+	int lenlength;// Length of the Output Dimension Vector "length"
+	int J; // Number of decomposition Levels
+	int MaxIter;// Maximum Iterations J <= MaxIter
+	int even;// even = 1 if signal is of even length. even = 0 otherwise
+	char ext[10];// Type of Extension used - "per" or "sym"
+
+	int N; //
+	int nodes;
+	int cfftset;
+	int zpad;
+	int length[102];
+	double *output;
+	int *nodelength;
+	int *coeflength;
+	double params[0];
+};
+
+typedef struct wpt_set* wpt_object;
+
+wpt_object wpt_init(wave_object wave, int siglength, int J);
+
+struct wpt_set{
+	wave_object wave;
+	conv_object cobj;
+	int siglength;// Length of the original signal.
+	int outlength;// Length of the output DWT vector
+	int lenlength;// Length of the Output Dimension Vector "length"
+	int J; // Number of decomposition Levels
+	int MaxIter;// Maximum Iterations J <= MaxIter
+	int even;// even = 1 if signal is of even length. even = 0 otherwise
+	char ext[10];// Type of Extension used - "per" or "sym"
+	char entropy[20];
+	double eparam;
+
+	int N; //
+	int nodes;
+	int length[102];
+	double *output;
+	double *costvalues;
+	double *basisvector;
+	int *nodeindex;
+	int *numnodeslevel;
+	int *coeflength;
+	double params[0];
+};
+
+typedef struct cwt_set* cwt_object;
+
+cwt_object cwt_init(char* wave, double param, int siglength,double dt, int J);
+
+struct cwt_set{
+	char wave[10];// Wavelet - morl/morlet,paul,dog/dgauss
+	int siglength;// Length of Input Data
+	int J;// Total Number of Scales
+	double s0;// Smallest scale. It depends on the sampling rate. s0 <= 2 * dt for most wavelets
+	double dt;// Sampling Rate
+	double dj;// Separation between scales. eg., scale = s0 * 2 ^ ( [0:N-1] *dj ) or scale = s0 *[0:N-1] * dj
+	char type[10];// Scale Type - Power or Linear
+	int pow;// Base of Power in case type = pow. Typical value is pow = 2
+	int sflag;
+	int pflag;
+	int npad;
+	int mother;
+	double m;// Wavelet parameter param
+	double smean;// Input Signal mean
+
+	cplx_data *output;
+	double *scale;
+	double *period;
+	double *coi;
+	double params[0];
+};
+
+
+void dwt(wt_object wt, double *inp);
+
+void idwt(wt_object wt, double *dwtop);
+
+void wtree(wtree_object wt, double *inp);
+
+void dwpt(wpt_object wt, double *inp);
+
+void idwpt(wpt_object wt, double *dwtop);
+
+void swt(wt_object wt, double *inp);
+
+void iswt(wt_object wt, double *swtop);
+
+void modwt(wt_object wt, double *inp);
+
+void imodwt(wt_object wt, double *dwtop);
+
+void setDWTExtension(wt_object wt, char *extension);
+
+void setWTREEExtension(wtree_object wt, char *extension);
+
+void setDWPTExtension(wpt_object wt, char *extension);
+
+void setDWPTEntropy(wpt_object wt, char *entropy, double eparam);
+
+void setWTConv(wt_object wt, char *cmethod);
+
+int getWTREENodelength(wtree_object wt, int X);
+
+void getWTREECoeffs(wtree_object wt, int X, int Y, double *coeffs, int N);
+
+int getDWPTNodelength(wpt_object wt, int X);
+
+void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N);
+
+void setCWTScales(cwt_object wt, double s0, double dj, char *type, int power);
+
+void setCWTScaleVector(cwt_object wt, double *scale, int J, double s0, double dj);
+
+void setCWTPadding(cwt_object wt, int pad);
+
+void cwt(cwt_object wt, double *inp);
+
+void icwt(cwt_object wt, double *cwtop);
+
+int getCWTScaleLength(int N);
+
+void wave_summary(wave_object obj);
+
+void wt_summary(wt_object wt);
+
+void wtree_summary(wtree_object wt);
+
+void wpt_summary(wpt_object wt);
+
+void cwt_summary(cwt_object wt);
+
+void wave_free(wave_object object);
+
+void wt_free(wt_object object);
+
+void wtree_free(wtree_object object);
+
+void wpt_free(wpt_object object);
+
+void cwt_free(cwt_object object);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */
diff --git a/test/cwttest.c b/test/cwttest.c
new file mode 100644
index 0000000..2da4ddc
--- /dev/null
+++ b/test/cwttest.c
@@ -0,0 +1,109 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+int main() {
+	int i, N, J,subscale,a0,iter,nd,k;
+	double *inp,*oup;
+	double dt, dj,s0, param,mn;
+	double td,tn,den, num, recon_mean, recon_var;
+	cwt_object wt;
+
+	FILE *ifp;
+	double temp[1200];
+
+	char *wave = "morlet";// Set Morlet wavelet. Other options "paul" and "dog"
+	char *type = "pow";
+
+	N = 504;
+	param = 6.0;
+	subscale = 4;
+	dt = 0.25;
+	s0 = dt;
+	dj = 1.0 / (double)subscale;
+	J = 11 * subscale; // Total Number of scales
+	a0 = 2;//power
+
+	ifp = fopen("sst_nino3.dat", "r");
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+
+	fclose(ifp);
+
+	wt = cwt_init(wave, param, N,dt, J);
+
+	inp = (double*)malloc(sizeof(double)* N);
+	oup = (double*)malloc(sizeof(double)* N);
+
+	for (i = 0; i < N; ++i) {
+		inp[i] = temp[i] ;
+	}
+
+	setCWTScales(wt, s0, dj, type, a0);
+
+	cwt(wt, inp);
+
+	printf("\n MEAN %g \n", wt->smean);
+	
+	mn = 0.0;
+
+	for (i = 0; i < N; ++i) {
+		mn += sqrt(wt->output[i].re * wt->output[i].re + wt->output[i].im * wt->output[i].im);
+	}
+	
+	cwt_summary(wt);
+
+	printf("\n abs mean %g \n", mn / N);
+	
+	printf("\n\n");
+	printf("Let CWT w = w(j, n/2 - 1) where n = %d\n\n", N);
+	nd = N/2 - 1;
+	
+	printf("%-15s%-15s%-15s%-15s \n","j","Scale","Period","ABS(w)^2");
+	for(k = 0; k < wt->J;++k) {
+		iter = nd + k * N;
+		printf("%-15d%-15lf%-15lf%-15lf \n",k,wt->scale[k],wt->period[k],
+		wt->output[iter].re * wt->output[iter].re + wt->output[iter].im * wt->output[iter].im);
+	}
+
+	icwt(wt, oup);
+
+	num = den = recon_var = recon_mean = 0.0;
+	printf("\n\n");
+	printf("Signal Reconstruction\n");
+	printf("%-15s%-15s%-15s \n","i","Input(i)","Output(i)");
+
+	for (i = N - 10; i < N; ++i) {
+		printf("%-15d%-15lf%-15lf \n", i,inp[i] , oup[i]);
+	}
+
+	for (i = 0; i < N; ++i) {
+		//printf("%g %g \n", oup[i] ,inp[i] - wt->smean);
+		td = inp[i] ;
+		tn = oup[i] - td;
+		num += (tn * tn);
+		den += (td * td);
+		recon_mean += oup[i];
+	}
+
+	recon_var = sqrt(num / N);
+	recon_mean /= N;
+
+	printf("\nRMS Error %g \n", sqrt(num) / sqrt(den));
+	printf("\nVariance %g \n", recon_var);
+	printf("\nMean %g \n", recon_mean);
+
+	free(inp);
+	free(oup);
+	cwt_free(wt);
+	return 0;
+}
diff --git a/test/sst_nino3.dat b/test/sst_nino3.dat
new file mode 100755
index 0000000..a60f6c5
--- /dev/null
+++ b/test/sst_nino3.dat
@@ -0,0 +1,504 @@
+-0.15
+-0.30
+-0.14
+-0.41
+-0.46
+-0.66
+-0.50
+-0.80
+-0.95
+-0.72
+-0.31
+-0.71
+-1.04
+-0.77
+-0.86
+-0.84
+-0.41
+-0.49
+-0.48
+-0.72
+-1.21
+-0.80
+ 0.16
+ 0.46
+ 0.40
+ 1.00
+ 2.17
+ 2.50
+ 2.34
+ 0.80
+ 0.14
+-0.06
+-0.34
+-0.71
+-0.34
+-0.73
+-0.48
+-0.11
+ 0.22
+ 0.51
+ 0.51
+ 0.25
+-0.10
+-0.33
+-0.42
+-0.23
+-0.53
+-0.44
+-0.30
+ 0.15
+ 0.09
+ 0.19
+-0.06
+ 0.25
+ 0.30
+ 0.81
+ 0.26
+ 0.10
+ 0.34
+ 1.01
+-0.31
+-0.90
+-0.73
+-0.92
+-0.73
+-0.31
+-0.03
+ 0.12
+ 0.37
+ 0.82
+ 1.22
+ 1.83
+ 1.60
+ 0.34
+-0.72
+-0.87
+-0.85
+-0.40
+-0.39
+-0.65
+ 0.07
+ 0.67
+ 0.39
+ 0.03
+-0.17
+-0.76
+-0.87
+-1.36
+-1.10
+-0.99
+-0.78
+-0.93
+-0.87
+-0.44
+-0.34
+-0.50
+-0.39
+-0.04
+ 0.42
+ 0.62
+ 0.17
+ 0.23
+ 1.03
+ 1.54
+ 1.09
+ 0.01
+ 0.12
+-0.27
+-0.47
+-0.41
+-0.37
+-0.36
+-0.39
+ 0.43
+ 1.05
+ 1.58
+ 1.25
+ 0.86
+ 0.60
+ 0.21
+ 0.19
+-0.23
+-0.29
+ 0.18
+ 0.12
+ 0.71
+ 1.42
+ 1.59
+ 0.93
+-0.25
+-0.66
+-0.95
+-0.47
+ 0.06
+ 0.70
+ 0.81
+ 0.78
+ 1.43
+ 1.22
+ 1.05
+ 0.44
+-0.35
+-0.67
+-0.84
+-0.66
+-0.45
+-0.12
+-0.20
+-0.16
+-0.47
+-0.52
+-0.79
+-0.80
+-0.62
+-0.86
+-1.29
+-1.04
+-1.05
+-0.75
+-0.81
+-0.90
+-0.25
+ 0.62
+ 1.22
+ 0.96
+ 0.21
+-0.11
+-0.25
+-0.24
+-0.43
+ 0.23
+ 0.67
+ 0.78
+ 0.41
+ 0.98
+ 1.28
+ 1.45
+ 1.02
+ 0.03
+-0.59
+-1.34
+-0.99
+-1.49
+-1.74
+-1.33
+-0.55
+-0.51
+-0.36
+-0.99
+ 0.32
+ 1.04
+ 1.41
+ 0.99
+ 0.66
+ 0.50
+ 0.22
+ 0.71
+-0.16
+ 0.38
+ 0.00
+-1.11
+-1.04
+ 0.05
+-0.64
+-0.34
+-0.50
+-1.85
+-0.94
+-0.78
+ 0.29
+ 0.27
+ 0.69
+-0.06
+-0.83
+-0.80
+-1.02
+-0.96
+-0.09
+ 0.62
+ 0.87
+ 1.03
+ 0.70
+-0.10
+-0.31
+ 0.04
+-0.46
+ 0.04
+ 0.24
+-0.08
+-0.28
+ 0.06
+ 0.05
+-0.31
+ 0.11
+ 0.27
+ 0.26
+ 0.04
+ 0.12
+ 1.11
+ 1.53
+ 1.23
+ 0.17
+-0.18
+-0.56
+ 0.05
+ 0.41
+ 0.22
+ 0.04
+-0.19
+-0.46
+-0.65
+-1.06
+-0.54
+ 0.14
+ 0.25
+-0.21
+-0.73
+-0.43
+ 0.48
+ 0.26
+ 0.05
+ 0.11
+-0.27
+-0.08
+-0.10
+ 0.29
+-0.15
+-0.28
+-0.55
+-0.44
+-1.40
+-0.55
+-0.69
+ 0.58
+ 0.37
+ 0.42
+ 1.83
+ 1.23
+ 0.65
+ 0.41
+ 1.03
+ 0.64
+-0.07
+ 0.98
+ 0.36
+-0.30
+-1.33
+-1.39
+-0.94
+ 0.34
+-0.00
+-0.15
+ 0.06
+ 0.39
+ 0.36
+-0.49
+-0.53
+ 0.35
+ 0.07
+-0.24
+ 0.20
+-0.22
+-0.68
+-0.44
+ 0.02
+-0.22
+-0.30
+-0.59
+ 0.10
+-0.02
+-0.27
+-0.60
+-0.48
+-0.37
+-0.53
+-1.35
+-1.22
+-0.99
+-0.34
+-0.79
+-0.24
+ 0.02
+ 0.69
+ 0.78
+ 0.17
+-0.17
+-0.29
+-0.27
+ 0.31
+ 0.44
+ 0.38
+ 0.24
+-0.13
+-0.89
+-0.76
+-0.71
+-0.37
+-0.59
+-0.63
+-1.47
+-0.40
+-0.18
+-0.37
+-0.43
+-0.06
+ 0.61
+ 1.33
+ 1.19
+ 1.13
+ 0.31
+ 0.14
+ 0.03
+ 0.21
+ 0.15
+-0.22
+-0.02
+ 0.03
+-0.17
+ 0.12
+-0.35
+-0.06
+ 0.38
+-0.45
+-0.32
+-0.33
+-0.49
+-0.14
+-0.56
+-0.18
+ 0.46
+ 1.09
+ 1.04
+ 0.23
+-0.99
+-0.59
+-0.92
+-0.28
+ 0.52
+ 1.31
+ 1.45
+ 0.61
+-0.11
+-0.18
+-0.39
+-0.39
+-0.36
+-0.50
+-0.81
+-1.10
+-0.29
+ 0.57
+ 0.68
+ 0.78
+ 0.78
+ 0.63
+ 0.98
+ 0.49
+-0.42
+-1.34
+-1.20
+-1.18
+-0.65
+-0.42
+-0.97
+-0.28
+ 0.77
+ 1.77
+ 2.22
+ 1.05
+-0.67
+-0.99
+-1.52
+-1.17
+-0.22
+-0.04
+-0.45
+-0.46
+-0.75
+-0.70
+-1.38
+-1.15
+-0.01
+ 0.97
+ 1.10
+ 0.68
+-0.02
+-0.04
+ 0.47
+ 0.30
+-0.55
+-0.51
+-0.09
+-0.01
+ 0.34
+ 0.61
+ 0.58
+ 0.33
+ 0.38
+ 0.10
+ 0.18
+-0.30
+-0.06
+-0.28
+ 0.12
+ 0.58
+ 0.89
+ 0.93
+ 2.39
+ 2.44
+ 1.92
+ 0.64
+-0.24
+ 0.27
+-0.13
+-0.16
+-0.54
+-0.13
+-0.37
+-0.78
+-0.22
+ 0.03
+ 0.25
+ 0.31
+ 1.03
+ 1.10
+ 1.05
+ 1.11
+ 1.28
+ 0.57
+-0.55
+-1.16
+-0.99
+-0.38
+ 0.01
+-0.29
+ 0.09
+ 0.46
+ 0.57
+ 0.24
+ 0.39
+ 0.49
+ 0.86
+ 0.51
+ 0.95
+ 1.25
+ 1.33
+-0.00
+ 0.34
+ 0.66
+ 1.11
+ 0.34
+ 0.48
+ 0.56
+ 0.39
+-0.17
+ 1.04
+ 0.77
+ 0.12
+-0.35
+-0.22
+ 0.08
+-0.08
+-0.18
+-0.06
diff --git a/unitTests/wavelibBoostTests/tst_dwt.cpp b/unitTests/wavelibBoostTests/tst_dwt.cpp
index c0de885..6fe71e6 100644
--- a/unitTests/wavelibBoostTests/tst_dwt.cpp
+++ b/unitTests/wavelibBoostTests/tst_dwt.cpp
@@ -98,6 +98,17 @@ double RMS_Error(double *data, double *rec, int N) {
     return sqrt(sum/((double)N-1));
 }
 
+double REL_Error(double *data, double *rec, int N) {
+    int i;
+    double sum1 = 0;
+    double sum2 = 0;
+    for (i = 0; i < N; ++i) {
+        sum1 += (data[i] - rec[i])*(data[i] - rec[i]);
+        sum2 += data[i] * data[i];
+    }
+    return sqrt(sum1)/sqrt(sum2);
+}
+
 void ReconstructionTest()
 {
 
@@ -334,6 +345,77 @@ void DWPTReconstructionTest()
     free(inp);
 }
 
+void CWTReconstructionTest() {
+	int i, j,N, J,subscale,a0,iter;
+	double *inp,*oup;
+	double dt, dj,s0, pi,t;
+	double val, epsilon;
+	int it1,it2;
+	cwt_object wt;
+	
+
+	char *wave[3];
+	wave[0] = (char*) "morl";
+	wave[1] =(char*) "paul";
+	wave[2] = (char*) "dog";
+	double param[30] = {4.5,5,5.5,6,6.5,8,10,13,17,20,
+		4,5,7,8,10,12,13,14,17,20,2,4,6,8,10,12,14,16,18,20};
+	char *type = (char*) "pow";
+	
+	epsilon = 0.01;
+	N = 2048;
+	dt = 0.000125;
+	subscale = 20;
+	dj = 1.0 / (double) subscale;
+	s0 = dt/32;	
+	J = 32 * subscale;
+	a0 = 2;//power
+	
+
+	inp = (double*)malloc(sizeof(double)* N);
+	oup = (double*)malloc(sizeof(double)* N);
+	
+	pi = 4.0 * atan(1.0);
+
+
+	for (i = 0; i < N; ++i) {
+		t = dt * i;
+		inp[i] = sin(2 * pi * 500 * t) + sin(2 * pi * 1000 * t) + 0.1 * sin(2 * pi * 8 * t);
+		if (i == 1200 || i ==1232) {
+			inp[i] += 5.0;
+		}
+	}
+	
+	for(it1 = 0; it1 < 3;++it1) {
+		for(it2 = 0; it2 < 10;++it2) {
+	
+			wt = cwt_init(wave[it1], param[it1*10+it2], N,dt, J);
+
+			setCWTScales(wt, s0, dj, type, a0);
+
+			cwt(wt, inp);
+
+
+			icwt(wt, oup);
+
+
+			//printf("\nWavelet : %s Parameter %g Error %g \n", wave[it1],param[it1*10+it2],REL_Error(inp,oup, wt->siglength));
+			if (REL_Error(inp,oup, wt->siglength) > epsilon) {
+				printf("\n ERROR : DWPT Reconstruction Unit Test Failed. Exiting. \n");
+				exit(-1);
+			}
+	
+			cwt_free(wt);
+		}
+	}
+	
+
+
+	free(inp);
+	free(oup);
+	
+}
+
 
 void DBCoefTests()
 {
@@ -569,6 +651,9 @@ int main() {
 	printf("Running DWPT ReconstructionTests ... ");
 	DWPTReconstructionTest();
 	printf("DONE \n");
+	printf("Running CWT ReconstructionTests ... ");
+	CWTReconstructionTest();
+	printf("DONE \n");
 	printf("\n\nUnit Tests Successful\n\n");
 	return 0;
 	
diff --git a/wavelib-doc.pdf b/wavelib-doc.pdf
index 0a98f87..2470781 100644
Binary files a/wavelib-doc.pdf and b/wavelib-doc.pdf differ