From bd2314727e85c4b65637637219b508d5c8e778bb Mon Sep 17 00:00:00 2001
From: Rafat Hussain <rafat.hsn@gmail.com>
Date: Fri, 22 Sep 2017 17:01:13 +0530
Subject: [PATCH] commit : sureshrink corrected

---
 auxiliary/denoise.c | 39 ++++++++++++++++++++++++---------------
 1 file changed, 24 insertions(+), 15 deletions(-)

diff --git a/auxiliary/denoise.c b/auxiliary/denoise.c
index 5c9829c..7f5e1eb 100644
--- a/auxiliary/denoise.c
+++ b/auxiliary/denoise.c
@@ -84,11 +84,11 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 }
 
 void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,double *denoised) {
-	int filt_len,i,it,len,dlen,dwt_len,min_index,sgn, MaxIter;
+	int filt_len,i,it,len,dlen,dwt_len,min_index,sgn, MaxIter,iter;
 	double sigma,norm,td,tv,te,ct,thr,temp,x_sum;
 	wave_object wave;
 	wt_object wt;
-	double *dout,*risk;
+	double *dout,*risk,*dsum;
 
 	wave = wave_init(wname);
 
@@ -119,15 +119,23 @@ void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 
 	dout = (double*)malloc(sizeof(double) * dlen);
 	risk = (double*)malloc(sizeof(double) * dlen);
+	dsum = (double*)malloc(sizeof(double) * dlen);
+
+	iter = wt->length[0];
+
+	for (i = 1; i < J; ++i) {
+		iter += wt->length[i];
+	}
+
+	for(i = 0; i < dlen;++i) {
+		dout[i] = fabs(wt->output[iter+i]);
+	}
+
+	sigma = median(dout,dlen) / 0.6745;
 
 	for(it = 0; it < J;++it) {
 		dwt_len = wt->length[it+1];
 
-		for(i = 0; i < dwt_len;++i) {
-			dout[i] = fabs(wt->output[len+i]);
-		}
-
-		sigma = median(dout,dwt_len);
 
 		if ( sigma < 0.00000001) {
 			td = 0;
@@ -135,7 +143,7 @@ void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 			tv = sqrt(2.0 * log(dwt_len));
 			norm = 0.0;
 			for(i = 0; i < dwt_len;++i) {
-				norm += (wt->output[len+i] *wt->output[len+i]);
+				norm += (wt->output[len+i] *wt->output[len+i] /(sigma*sigma));
 			}
 			te =(norm - (double) dwt_len)/(double) dwt_len;
 			ct = pow(log((double) dwt_len)/log(2.0),1.5)/sqrt((double) dwt_len);
@@ -146,37 +154,37 @@ void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 				x_sum = 0.0;
 
 				for(i = 0; i < dwt_len;++i) {
-					dout[i] = wt->output[len+i];
+					dout[i] = fabs(wt->output[len+i]/sigma);
 				}
 
 				qsort(dout, dwt_len, sizeof(double), compare_double);
 				for(i = 0; i < dwt_len;++i) {
 					dout[i] = (dout[i]*dout[i]);
 					x_sum += dout[i];
+					dsum[i] = x_sum;
 				}
 
 				for(i = 0;i < dwt_len;++i) {
-					risk[i] = ((double)dwt_len + 1 - 2 * ((double)i + 1) +x_sum +
+					risk[i] = ((double)dwt_len - 2 * ((double)i + 1) +dsum[i] +
 							dout[i]*((double)dwt_len - 1 -(double) i))/(double)dwt_len;
 				}
-
+				printf(" \n");
 				min_index = minindex(risk,dwt_len);
 				thr = sqrt(dout[min_index]);
-
 				td = thr < tv ? thr : tv;
 			}
 		}
 
-		td = td * sigma / 0.6745;
+		td = td * sigma;
 
 		if(!strcmp(thresh,"hard")) {
-			for(i = wt->length[0]; i < dwt_len;++i) {
+			for(i = 0; i < dwt_len;++i) {
 				if (fabs(wt->output[len+i]) < td) {
 					wt->output[len+i] = 0;
 				}
 			}
 		} else if(!strcmp(thresh,"soft")) {
-			for(i = wt->length[0]; i < dwt_len;++i) {
+			for(i = 0; i < dwt_len;++i) {
 					if (fabs(wt->output[len + i]) < td) {
 						wt->output[len+i] = 0;
 					} else {
@@ -199,6 +207,7 @@ void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 	}
 
 	free(dout);
+	free(dsum);
 	free(risk);
 	wave_free(wave);
 	wt_free(wt);