commit : DWPT and Wavelet Tree added

README.md

@@ -1,9 +1,9 @@
 wavelib
 =======
 
-C Implementation of Wavelet Transform (DWT,SWT and MODWT)
+C Implementation of Wavelet Transform (DWT,SWT and MODWT) and Packet Transform ( Full Tree Decomposition and Best Basis DPWT).
 
-Methods Implemented
+Discrete Wavelet Transform Methods Implemented
 
 DWT/IDWT A decimated Discrete Wavelet Transform implementation using implicit signal extension and up/downsampling so it is a fast implementation. A FFT based implementation is optional but will not be usually needed. Both periodic and symmetric options are available.
 
@@ -11,6 +11,12 @@ SWT/ISWT Stationary Wavelet Transform. It works only for signal lengths that are
 
 MODWT/IMODWT Maximal Overlap Discrete Wavelet Transform is another undecimated transform. It is implemented for signals of any length but only orthogonal wavelets (Daubechies, Symlets and Coiflets) can be deployed. This implementation is based on the method laid out in "Wavelet Methods For Wavelet Analysis" by Donald Percival and Andrew Walden. 
 
+Discrete Wavelet Packet Transform Methods Implemented
+
+WTREE A Fully Decimated Wavelet Tree Decomposition. This is a highly redundant transform and retains all coefficients at each node. This is not recommended for compression and denoising applications.
+
+DWPT/IDWPT Is a derivative of WTREE method which retains coefficients based on entropy methods. This is a non-redundant transform and output length is of the same order as the input.
+
 Documentation Available at - https://github.com/rafat/wavelib/wiki
 
 Live Demo (Emscripten) - http://rafat.github.io/wavelib/

header/wavelib.h

@@ -83,21 +83,81 @@ struct wt_set{
 	char ext[10];// Type of Extension used - "per" or "sym"
 	char cmethod[10]; // Convolution Method - "direct" or "FFT"
 
-	int N; // 
-	int cfftset; 
+	int N; //
+	int cfftset;
 	int zpad;
 	int length[102];
 	double *output;
 	double params[0];
 };
 
-void dwt11(wt_object wt, double *Vin, int M, double *Wout,
-	double *Vout);
+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);
@@ -108,16 +168,38 @@ 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
 }
@@ -125,5 +207,3 @@ void wt_free(wt_object object);
 
 
 #endif /* WAVELIB_H_ */
-
-

src/wavelib.h

@@ -42,18 +42,81 @@ struct wt_set{
 	char ext[10];// Type of Extension used - "per" or "sym"
 	char cmethod[10]; // Convolution Method - "direct" or "FFT"
 
-	int N; // 
-	int cfftset; 
+	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);
@@ -64,16 +127,38 @@ 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
 }
@@ -81,5 +166,3 @@ void wt_free(wt_object object);
 
 
 #endif /* WAVELIB_H_ */
-
-

src/wtmath.c

@@ -236,4 +236,95 @@ int wmaxiter(int sig_len, int filt_len) {
 	lev = (int)temp;
 
 	return lev;
-}
+}
+
+static double entropy_s(double *x,int N) {
+  int i;
+  double val,x2;
+
+  val = 0.0;
+
+  for(i = 0; i < N; ++i) {
+    if (x[i] != 0) {
+      x2 = x[i] * x[i];
+      val -= x2 * log(x2);
+    }
+  }
+  return val;
+}
+
+static double entropy_t(double *x,int N, double t) {
+  int i;
+  double val,x2;
+  if (t < 0) {
+    printf("Threshold value must be >= 0");
+    exit(1);
+  }
+  val = 0.0;
+
+  for(i = 0; i < N; ++i) {
+    x2 = fabs(x[i]);
+    if (x2 > t) {
+      val += 1;
+    }
+
+  }
+
+  return val;
+
+}
+
+static double entropy_n(double *x,int N,double p) {
+  int i;
+  double val,x2;
+  if (p < 1) {
+    printf("Norm power value must be >= 1");
+    exit(1);
+  }
+  val = 0.0;
+  for(i = 0; i < N; ++i) {
+    x2 = fabs(x[i]);
+    val += pow(x2,(double)p);
+
+  }
+
+  return val;
+}
+
+static double entropy_l(double *x,int N) {
+  int i;
+  double val,x2;
+
+  val = 0.0;
+
+  for(i = 0; i < N; ++i) {
+    if (x[i] != 0) {
+      x2 = x[i] * x[i];
+      val += log(x2);
+    }
+  }
+  return val;
+}
+
+double costfunc(double *x, int N ,char *entropy,double p) {
+	double val;
+
+	if (!strcmp(entropy, "shannon")) {
+		val = entropy_s(x, N);
+	}
+	else if (!strcmp(entropy, "threshold")) {
+		val = entropy_t(x, N,p);
+	}
+	else if (!strcmp(entropy, "norm")) {
+		val = entropy_n(x, N,p);
+	}
+	else if (!strcmp(entropy, "logenergy") || !strcmp(entropy, "log energy") || !strcmp(entropy, "energy")) {
+		val = entropy_l(x, N);
+	}
+	else {
+		printf("Entropy must be one of shannon, threshold, norm or energy");
+		exit(-1);
+	}
+
+	return val;
+}

src/wtmath.h

@@ -23,9 +23,11 @@ int testSWTlength(int N, int J);
 
 int wmaxiter(int sig_len, int filt_len);
 
+double costfunc(double *x, int N, char *entropy, double p);
+
 #ifdef __cplusplus
 }
 #endif
 
 
-#endif /* WAVELIB_H_ */
+#endif /* WAVELIB_H_ */

test/dwpttest.c

@@ -0,0 +1,61 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+int main() {
+	int i, J, N;
+	wave_object obj;
+	wpt_object wt;
+	double *inp, *oup, *diff;
+
+	char *name = "db4";
+	obj = wave_init(name);// Initialize the wavelet
+	N = 788 + 23;
+	inp = (double*)malloc(sizeof(double)* N);
+	oup = (double*)malloc(sizeof(double)* N);
+	diff = (double*)malloc(sizeof(double)* N);
+	for (i = 1; i < N + 1; ++i) {
+		//inp[i - 1] = -0.25*i*i*i + 25 * i *i + 10 * i;
+		inp[i - 1] = i;
+	}
+	J = 4;
+
+	wt = wpt_init(obj, N, J);// Initialize the wavelet transform Tree object
+	setDWPTExtension(wt, "per");// Options are "per" and "sym". Symmetric is the default option
+	setDWPTEntropy(wt, "logenergy", 0);
+
+	dwpt(wt, inp); // Discrete Wavelet Packet Transform
+
+	idwpt(wt, oup); // Inverse Discrete Wavelet Packet Transform
+
+	for (i = 0; i < N; ++i) {
+		diff[i] = (inp[i] - oup[i])/inp[i];
+	}
+
+	wpt_summary(wt); // Tree Summary
+
+	printf("\n MAX %g \n", absmax(diff, wt->siglength)); // If Reconstruction succeeded then the output should be a small value.
+
+	free(inp);
+	free(oup);
+	free(diff);
+	wave_free(obj);
+	wpt_free(wt);
+	return 0;
+}

test/wtreetest.c

@@ -0,0 +1,47 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+int main() {
+	int i, J, N, len;
+	int X, Y;
+	wave_object obj;
+	wtree_object wt;
+	double *inp, *oup;
+
+	char *name = "db3";
+	obj = wave_init(name);// Initialize the wavelet
+	N = 147;
+	inp = (double*)malloc(sizeof(double)* N);
+	for (i = 1; i < N + 1; ++i) {
+		inp[i - 1] = -0.25*i*i*i + 25 * i *i + 10 * i;
+	}
+	J = 3;
+
+	wt = wtree_init(obj, N, J);// Initialize the wavelet transform object
+	setWTREEExtension(wt, "sym");// Options are "per" and "sym". Symmetric is the default option
+
+	wtree(wt, inp);
+	wtree_summary(wt);
+	X = 3;
+	Y = 5;
+	len = getWTREENodelength(wt, X);
+	printf("\n %d", len);
+	printf("\n");
+	oup = (double*)malloc(sizeof(double)* len);
+
+	printf("Node [%d %d] Coefficients : \n",X,Y);
+	getWTREECoeffs(wt, X, Y, oup, len);
+	for (i = 0; i < len; ++i) {
+		printf("%g ", oup[i]);
+	}
+	printf("\n");
+
+	free(inp);
+	free(oup);
+	wave_free(obj);
+	wtree_free(wt);
+	return 0;
+}