commit : visushrink updated

auxiliary/denoise.c

@@ -43,22 +43,23 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 	}
 
 	for(i = 0; i < dlen;++i) {
-		dout[i] = wt->output[iter+i];
+		dout[i] = fabs(wt->output[iter+i]);
 	}
 
-	sigma = mad(dout,dlen);
+	sigma = median(dout,dlen);
 	dwt_len = wt->outlength;
 
 	td = sqrt(2.0 * log(dwt_len)) * sigma / 0.6745;
 
+
 	if(!strcmp(thresh,"hard")) {
-		for(i = 0; i < dwt_len;++i) {
+		for(i = wt->length[0]; i < dwt_len;++i) {
 			if (fabs(wt->output[i]) < td) {
 				wt->output[i] = 0;
 			}
 		}
 	} else if(!strcmp(thresh,"soft")) {
-		for(i = 0; i < dwt_len;++i) {
+		for(i = wt->length[0]; i < dwt_len;++i) {
 				if (fabs(wt->output[i]) < td) {
 					wt->output[i] = 0;
 				} else {
@@ -123,10 +124,10 @@ void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 		dwt_len = wt->length[it+1];
 
 		for(i = 0; i < dwt_len;++i) {
-			dout[i] = wt->output[len+i];
+			dout[i] = fabs(wt->output[len+i]);
 		}
 
-		sigma = mad(dout,dwt_len);
+		sigma = median(dout,dwt_len);
 
 		if ( sigma < 0.00000001) {
 			td = 0;
@@ -169,13 +170,13 @@ void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
 		td = td * sigma / 0.6745;
 
 		if(!strcmp(thresh,"hard")) {
-			for(i = 0; i < dwt_len;++i) {
+			for(i = wt->length[0]; i < dwt_len;++i) {
 				if (fabs(wt->output[len+i]) < td) {
 					wt->output[len+i] = 0;
 				}
 			}
 		} else if(!strcmp(thresh,"soft")) {
-			for(i = 0; i < dwt_len;++i) {
+			for(i = wt->length[0]; i < dwt_len;++i) {
 					if (fabs(wt->output[len + i]) < td) {
 						wt->output[len+i] = 0;
 					} else {

test/denoisetest.c

@@ -4,19 +4,60 @@
 #include <math.h>
 #include "../header/denoise.h"
 
+static double rmse(int N,double *x,double *y) {
+	double rms;
+	int i;
+
+	rms = 0.0;
+
+	for(i = 0; i < N;++i) {
+		rms += (x[i] - y[i]) * (x[i] - y[i]);
+	}
+
+	rms = sqrt(rms/(double)N);
+
+	return rms;
+}
+
+static double corrcoef(int N,double *x,double *y) {
+	double cc,xm,ym,tx,ty,num,den1,den2;
+	int i;
+	xm = ym = 0.0;
+	for(i = 0; i < N;++i) {
+		xm += x[i];
+		ym += y[i];
+	}
+
+	xm = xm/N;
+	ym = ym / N;
+	num = den1 = den2 = 0.0;
+
+	for(i = 0; i < N;++i) {
+		tx = x[i] - xm;
+		ty = y[i] - ym;
+		num += (tx*ty);
+		den1 += (tx*tx);
+		den2 += (ty*ty);
+	}
+
+	cc = num / sqrt(den1*den2);
+
+	return cc;
+}
+
 int main() {
 	// gcc -Wall -I../header -L../Bin denoisetest.c -o denoise -lwauxlib -lwavelib -lm
-	double *inp,*oup;
+	double *sig,*inp,*oup;
 	int i,N,J;
 	FILE *ifp,*ofp;
 	double temp[2400];
 
-	char *wname = "sym8";
+	char *wname = "sym6";
-	char *method = "swt";
+	char *method = "dwt";
 	char *ext = "sym";
 	char *thresh = "soft";
 
-	ifp = fopen("noisybumps.txt", "r");
+	ifp = fopen("pieceregular1024.txt", "r");
 	i = 0;
 	if (!ifp) {
 	    printf("Cannot Open File");
@@ -33,24 +74,49 @@ int main() {
 	N = i;
 	J = 6;
 
+	sig = (double*)malloc(sizeof(double)* N);
 	inp = (double*)malloc(sizeof(double)* N);
 	oup = (double*)malloc(sizeof(double)* N);
 
 	for(i = 0; i < N;++i) {
-		inp[i] = temp[i];
+		sig[i] = temp[i];
 	}
 
-	//visushrink(inp,N,wname,method,ext,thresh,oup);
+	ifp = fopen("PieceRegular10.txt", "r");
-	sureshrink(inp,N,J,wname,method,ext,thresh,oup);
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
 
-	ofp = fopen("denoiseds.txt", "w");
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+
+	fclose(ifp);
 
 	for(i = 0; i < N;++i) {
-		fprintf(ofp,"%g \n",oup[i]);
+		inp[i] = temp[i];
+	}
+	//visushrink(inp,N,J,wname,method,ext,thresh,oup);
+	sureshrink(inp,N,J,wname,method,ext,thresh,oup);
+
+	//ofp = fopen("denoiseds.txt", "w");
+
+	for(i = 0; i < N;++i) {
+		//fprintf(ofp,"%g \n",oup[i]);
 	}
 
-	fclose(ofp);
+	//fclose(ofp);
 
+	printf("RMSE %g\n",rmse(N,sig,inp));
+	printf("Corr Coeff %g\n",corrcoef(N,sig,inp));
+
+	printf("RMSE %g\n",rmse(N,sig,oup));
+		printf("Corr Coeff %g\n",corrcoef(N,sig,oup));
+
+	free(sig);
 	free(inp);
 	free(oup);
 	return 0;