commit : visualshrink code added

CMakeLists.txt

@@ -45,5 +45,6 @@ if(BUILD_UT)
 endif()
 
 add_subdirectory(src)
+add_subdirectory(denoise)
 
 install(DIRECTORY ${WAVELIB_SRC_ROOT}/include/ DESTINATION include FILES_MATCHING PATTERN "*.h")

denoise/denoise.c

@@ -1,5 +1,6 @@
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
 #include "../header/wavelib.h"
 
 int compare_double(const void* a, const void* b)
@@ -42,13 +43,106 @@ double mad(double *x, int N) {
 	return sigma;
 }
 
-void visushrink(wave_object wave,double *signal,int N) {
-	int J,filt_len;
-	double sigma;
+void getDWTAppx(wt_object wt, double *appx,int N) {
+	/*
+	Wavelet decomposition is stored as
+	[A(J) D(J) D(J-1) ..... D(1)] in wt->output vector
+
+	Length of A(J) , N = wt->length[0]
+	*/
+	int i;
+
+	for (i = 0; i < N; ++i) {
+		appx[i] = wt->output[i];
+	}
+}
+
+void getDWTDetail(wt_object wt, double *detail, int N, int level) {
+	/*
+	returns Detail coefficents at the jth level where j = 1,2,.., J
+	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)
+	Level 1 : Length of D(J), ie N, is stored in wt->length[1]
+	Level 2 :Length of D(J-1), ie N, is stored in wt->length[2]
+	....
+	Level J : Length of D(1), ie N, is stored in wt->length[J]
+	*/
+	int i, iter, J;
+	J = wt->J;
+
+	if (level > J) {
+		printf("The decomposition only has %d levels", J);
+	}
+
+	iter = wt->length[0];
+
+	for (i = 1; i < level; ++i) {
+		iter += wt->length[i];
+	}
+
+	for (i = 0; i < N; ++i) {
+		detail[i] = wt->output[i + iter];
+	}
+}
+
+void visushrink(double *signal,int N,char *wname,char *method,char *ext,char *thresh,double *denoised) {
+	int J,filt_len,iter,i,dlen,dwt_len,sgn;
+	double sigma,td,tmp;
+	wave_object wave;
+	wt_object wt;
+
+	wave = wave_init(wname);
 	
 	filt_len = wave->filtlength;
 	
 	J = (int) (log((double)N / ((double)filt_len - 1.0)) / log(2.0));
+
+	if (J > 50) {
+		J = 50;
+	}
+
+	wt = wt_init(wave,method,N,J);
+	setDWTExtension(wt,ext);
+
+	dwt(wt,signal);
+
+	//Set sigma
+
+	iter = wt->length[0];
+	dlen = wt->length[J];
+
+	for (i = 1; i < J; ++i) {
+		iter += wt->length[i];
+	}
+
+	sigma = mad(wt->output+iter,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) {
+			if (fabs(wt->output[i]) < td) {
+				wt->output[i] = 0;
+			}
+		}
+	} else if(!strcmp(thresh,"soft")) {
+		for(i = 0; i < dwt_len;++i) {
+				if (fabs(wt->output[i]) < td) {
+					wt->output[i] = 0;
+				} else {
+					sgn = wt->output[i] >= 0 ? 1 : -1;
+					tmp = sgn * (fabs(wt->output[i]) - td);
+					wt->output[i] = tmp;
+				}
+			}
+	}
+
+	idwt(wt,denoised);
+
+	wave_free(wave);
+	wt_free(wt);
 }
  
 int main(void)