diff --git a/README.md b/README.md index 025a18b..cb891f2 100644 --- a/README.md +++ b/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/ diff --git a/header/wavelib.h b/header/wavelib.h index be1581a..1e80b30 100644 --- a/header/wavelib.h +++ b/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_ */ - - diff --git a/src/wavelib.c b/src/wavelib.c index b7709a3..dc4f358 100644 --- a/src/wavelib.c +++ b/src/wavelib.c @@ -118,6 +118,145 @@ wt_object wt_init(wave_object wave,char* method, int siglength,int J) { 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); @@ -192,6 +331,119 @@ static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, do } + +} + +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) { @@ -221,6 +473,64 @@ static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, do } +} + +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) { @@ -312,15 +622,15 @@ void dwt(wt_object wt,double *inp) { 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); + 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); + wt->zpad = 1; + temp_len++; + orig = (double*)malloc(sizeof(double)* temp_len); + orig2 = (double*)malloc(sizeof(double)* temp_len); } */ @@ -415,6 +725,504 @@ void dwt(wt_object wt,double *inp) { 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; @@ -460,10 +1268,9 @@ static void idwt1(wt_object wt,double *temp, double *cA_up,double *cA, int len_c } static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { - int len_avg, N, i,l,m,n,t,v,l2; + int len_avg,i,l,m,n,t,l2; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; - N = 2 * len_cA; l2 = len_avg / 2; m = -2; n = -1; @@ -492,13 +1299,11 @@ static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_c } static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { - int len_avg, N, i, l, m, n, t, v,l2; + int len_avg, i, l, m, n, t, v; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; - N = 2 * len_cA - 1; m = -2; n = -1; - l2 = (len_avg - 2) / 2; for (v = 0; v < len_cA; ++v) { i = v; @@ -680,6 +1485,331 @@ void idwt(wt_object wt, double *dwtop) { } +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; + 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]; + 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]; + index2 *= 2; + index3 = 0; + index4 = 0; + //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); + n1 -= llb; + for (l = 0; l < llb; ++l) { + if (ptemp[l] != 2) { + 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) { + out[k] = wt->output[index + k]; + out2[k] = wt->output[index + det_len + k]; + } + idwpt_per(wt, out, det_len, out2, det_len, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index += 2 * det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) { + index4 += indexp; + for (k = 0; k < det_len; ++k) { + out[k] = wt->output[index + k]; + out2[k] = X[index4 + k]; + } + idwpt_per(wt, out, det_len, out2, det_len, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = wt->output[index + k]; + } + idwpt_per(wt, out, det_len, out2, det_len, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = X[index4 + indexp + k]; + } + idwpt_per(wt, out, det_len, out2, det_len, X_lp); + for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { + X[index3 + k - lf / 2 + 1] = X_lp[k]; + } + index4 += 2 * indexp; + index3 += index2; + ptemp[l] = 2; + } + else { + 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]; + } + + 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]; + index2 *= 2; + index3 = 0; + index4 = 0; + //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); + n1 -= llb; + for (l = 0; l < llb; ++l) { + if (ptemp[l] != 2) { + 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) { + out[k] = wt->output[index + k]; + out2[k] = wt->output[index + det_len + k]; + } + idwpt_sym(wt, out, det_len, out2, det_len, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index += 2 * det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) { + index4 += indexp; + for (k = 0; k < det_len; ++k) { + out[k] = wt->output[index + k]; + out2[k] = X[index4 + k]; + } + idwpt_sym(wt, out, det_len, out2, det_len, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = wt->output[index + k]; + } + idwpt_sym(wt, out, det_len, out2, det_len, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index += det_len; + index3 += index2; + index4 += 2 * indexp; + ptemp[l] = 2; + } + else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) { + for (k = 0; k < det_len; ++k) { + out[k] = X[index4 + k]; + out2[k] = X[index4 + indexp + k]; + } + idwpt_sym(wt, out, det_len, out2, det_len, X_lp); + for (k = lf - 2; k < 2 * det_len; ++k) { + X[index3 + k - lf + 2] = X_lp[k]; + } + index4 += 2 * indexp; + index3 += index2; + ptemp[l] = 2; + } + else { + 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]; + } + + 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; @@ -884,8 +2014,7 @@ void swt(wt_object wt, double *inp) { void iswt(wt_object wt, double *swtop) { int N, lf, iter,i,J,index,value,count,len; - int index_shift,len0,U,N1,N2,index2; - double *low_pass, *high_pass; + 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; @@ -1116,7 +2245,7 @@ void modwt(wt_object wt, double *inp) { } 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, j; + int len_avg, i, l, t; double s; double *filt; len_avg = wt->wave->lpd_len; @@ -1149,7 +2278,7 @@ static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, i } void imodwt(wt_object wt, double *dwtop) { - int N, lf, iter, i, J, index, value, j, len,U; + int N, lf, iter, i, J, j, U; int lenacc,M; double *X; @@ -1201,6 +2330,53 @@ void setDWTExtension(wt_object wt, char *extension) { } } +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"); @@ -1280,6 +2456,79 @@ void wt_summary(wt_object wt) { } +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); } @@ -1287,3 +2536,11 @@ void wave_free(wave_object 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); +} diff --git a/src/wavelib.h b/src/wavelib.h index 8cd6a41..6507d86 100644 --- a/src/wavelib.h +++ b/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_ */ - - diff --git a/src/wtmath.c b/src/wtmath.c index 42a81b5..dc1700e 100644 --- a/src/wtmath.c +++ b/src/wtmath.c @@ -236,4 +236,95 @@ int wmaxiter(int sig_len, int filt_len) { lev = (int)temp; return lev; -} \ No newline at end of file +} + +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; +} diff --git a/src/wtmath.h b/src/wtmath.h index 076c707..c7c7a52 100644 --- a/src/wtmath.h +++ b/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_ */ \ No newline at end of file +#endif /* WAVELIB_H_ */ diff --git a/test/dwpttest.c b/test/dwpttest.c new file mode 100644 index 0000000..c2d0547 --- /dev/null +++ b/test/dwpttest.c @@ -0,0 +1,61 @@ +#include +#include +#include +#include +#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; +} diff --git a/test/wtreetest.c b/test/wtreetest.c new file mode 100644 index 0000000..da3b75e --- /dev/null +++ b/test/wtreetest.c @@ -0,0 +1,47 @@ +#include +#include +#include +#include +#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; +} diff --git a/wavelib-doc.pdf b/wavelib-doc.pdf new file mode 100644 index 0000000..0a98f87 Binary files /dev/null and b/wavelib-doc.pdf differ