diff --git a/auxiliary/waux.c b/auxiliary/waux.c index 2249b7f..285761e 100644 --- a/auxiliary/waux.c +++ b/auxiliary/waux.c @@ -73,7 +73,7 @@ void getDWTAppx(wt_object wt, double *appx,int N) { void getDWTDetail(wt_object wt, double *detail, int N, int level) { /* - returns Detail coefficents at the jth level where j = 1,2,.., J + returns Detail coefficents at the jth level where j = J,J-1,...,1 and Wavelet decomposition is stored as [A(J) D(J) D(J-1) ..... D(1)] in wt->output vector Use getDWTAppx() to get A(J) @@ -85,13 +85,13 @@ void getDWTDetail(wt_object wt, double *detail, int N, int level) { int i, iter, J; J = wt->J; - if (level > J) { - printf("The decomposition only has %d levels", J); + if (level > J || level < 1) { + printf("The decomposition only has 1,..,%d levels", J); } iter = wt->length[0]; - for (i = 1; i < level; ++i) { + for (i = 1; i < J-level; ++i) { iter += wt->length[i]; } @@ -100,4 +100,143 @@ void getDWTDetail(wt_object wt, double *detail, int N, int level) { } } +void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, int J,double *lpr, + double *hpr,int lf,int siglength,double *reccoeff) { + int i,j,k,det_len,N,l,m,n,v,t,l2; + double *out,*X_lp,*filt; + out = (double*)malloc(sizeof(double)* (siglength + 1)); + l2 = lf / 2; + m = -2; + n = -1; + if (!strcmp(ext, "per")) { + if (!strcmp((ctype), "appx")) { + det_len = length[0]; + } else { + det_len = length[J - level + 1]; + } + + N = 2 * length[J]; + + X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); + + for (i = 0; i < det_len; ++i) { + out[i] = coeff[i]; + } + + for (j = level; j > 0; --j) { + + //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); + + if (!strcmp((ctype), "det") && j == level) { + filt = hpr; + } else { + filt = lpr; + } + + //idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp); + + for (i = 0; i < lf + l2 - 1; ++i) { + m += 2; + n += 2; + X_lp[m] = 0.0; + X_lp[n] = 0.0; + for (l = 0; l < l2; ++l) { + t = 2 * l; + if ((i - l) >= 0 && (i - l) < det_len) { + X_lp[m] += filt[t] * out[i - l]; + X_lp[n] += filt[t + 1] * out[i - l]; + } + else if ((i - l) >= det_len && (i-l) < det_len + lf - 1) { + X_lp[m] += filt[t] * out[i - l - det_len]; + X_lp[n] += filt[t + 1] * out[i - l - det_len]; + } + else if ((i - l) < 0 && (i-l) > -l2) { + X_lp[m] += filt[t] * out[det_len + i - l] ; + X_lp[n] += filt[t + 1] * out[det_len + i - l]; + } + } + } + + for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) { + out[k - lf/2 + 1] = X_lp[k]; + } + + if (j != 1) { + det_len = length[J - j + 2]; + } + } + + free(X_lp); + + } + else if (!strcmp(ext, "sym")) { + if (!strcmp((ctype), "appx")) { + det_len = length[0]; + } else { + det_len = length[J - level + 1]; + } + + N = 2 * length[J] - 1; + + X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); + + + for (i = 0; i < det_len; ++i) { + out[i] = coeff[i]; + } + + for (j = level; j > 0; --j) { + + //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); + + if (!strcmp((ctype), "det") && j == level) { + filt = hpr; + } else { + filt = lpr; + } + //idwt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp); + + m = -2; + n = -1; + + for (v = 0; v < det_len; ++v) { + i = v; + m += 2; + n += 2; + X_lp[m] = 0.0; + X_lp[n] = 0.0; + for (l = 0; l < lf / 2; ++l) { + t = 2 * l; + if ((i - l) >= 0 && (i - l) < det_len) { + X_lp[m] += filt[t] * out[i - l]; + X_lp[n] += filt[t + 1] * out[i - l] + filt[t + 1] * out[i - l]; + } + } + } + + for (k = lf-2; k < 2 * det_len; ++k) { + out[k - lf + 2] = X_lp[k]; + } + + + if (j != 1) { + det_len = length[J - j + 2]; + } + } + + free(X_lp); + + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } + + for (i = 0; i < siglength; ++i) { + reccoeff[i] = out[i]; + } + + free(out); + +} diff --git a/src/wavelib.c b/src/wavelib.c index bdb1d7c..2936572 100644 --- a/src/wavelib.c +++ b/src/wavelib.c @@ -2663,12 +2663,12 @@ void wt_summary(wt_object wt) { 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("Level %d Access : output[%d] Length : %d \n", J, 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]); + printf("Level %d Access : output[%d] Length : %d \n", J - i,t,wt->length[i+1]); t += wt->length[i+1]; } printf("\n");