From f906ae8d669adfe8a0154310ce8ec26284188381 Mon Sep 17 00:00:00 2001 From: Rafat Hussain Date: Sun, 24 Sep 2017 18:24:45 +0530 Subject: [PATCH] commit : denoise object added --- auxiliary/denoise.c | 106 ++++++++++++++++++++++++++++++++++++++++++-- auxiliary/denoise.h | 29 +++++++++++- auxiliary/waux.c | 78 ++++++++++++++++++++++++++++++++ auxiliary/waux.h | 9 +++- header/denoise.h | 30 ------------- header/wauxlib.h | 56 +++++++++++++++++++++++ test/denoisetest.c | 26 +++++++---- 7 files changed, 290 insertions(+), 44 deletions(-) delete mode 100644 header/denoise.h create mode 100644 header/wauxlib.h diff --git a/auxiliary/denoise.c b/auxiliary/denoise.c index ae9a7a5..677094c 100644 --- a/auxiliary/denoise.c +++ b/auxiliary/denoise.c @@ -1,6 +1,26 @@ #include "denoise.h" +denoise_object denoise_init(int length, int J,char* wname) { + denoise_object obj = NULL; + + obj = (denoise_object)malloc(sizeof(struct denoise_set) +sizeof(double)); + + obj->N = length; + obj->J = J; + + strcpy(obj->wname,wname); + + //Set Default Values + strcpy(obj->dmethod,"sureshrink"); + strcpy(obj->ext,"sym"); + strcpy(obj->level,"first"); + strcpy(obj->thresh,"soft"); + strcpy(obj->wmethod,"dwt"); + + return obj; +} + void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised) { int filt_len,iter,i,dlen,dwt_len,sgn, MaxIter,it; double sigma,td,tmp; @@ -25,8 +45,6 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch dwt(wt,signal); } else if(!strcmp(method,"swt")) { swt(wt,signal); - } else if(!strcmp(method,"modwt")) { - modwt(wt,signal); } else { printf("Acceptable WT methods are - dwt,swt and modwt\n"); exit(-1); @@ -102,8 +120,6 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch idwt(wt,denoised); } else if(!strcmp(method,"swt")) { iswt(wt,denoised); - } else if(!strcmp(method,"modwt")) { - imodwt(wt,denoised); } free(dout); @@ -259,3 +275,85 @@ void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch wave_free(wave); wt_free(wt); } + +void denoise(denoise_object obj, double *signal,double *denoised) { + if(!strcmp(obj->dmethod,"sureshrink")) { + sureshrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised); + } else if(!strcmp(obj->dmethod,"visushrink")) { + visushrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);; + } else { + printf("Acceptable Denoising methods are - sureshrink and visushrink\n"); + exit(-1); + } +} + +void setDenoiseMethod(denoise_object obj, char *dmethod) { + if (!strcmp(dmethod, "sureshrink")) { + strcpy(obj->dmethod, "sureshrink"); + } + else if (!strcmp(dmethod, "visushrink")) { + strcpy(obj->dmethod, "visushrink"); + } + else { + printf("Acceptable Denoising methods are - sureshrink and visushrink\n"); + exit(-1); + } +} + +void setDenoiseWTMethod(denoise_object obj, char *wmethod) { + if (!strcmp(wmethod, "dwt")) { + strcpy(obj->wmethod, "dwt"); + } + else if (!strcmp(wmethod, "swt")) { + strcpy(obj->wmethod, "swt"); + } + else { + printf("Wavelet decomposition method can be either dwt or swt"); + exit(-1); + } +} + +void setDenoiseWTExtension(denoise_object obj, char *extension) { + if (!strcmp(extension, "sym")) { + strcpy(obj->ext, "sym"); + } + else if (!strcmp(extension, "per")) { + strcpy(obj->ext, "per"); + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } +} + +void setDenoiseParameters(denoise_object obj, char *thresh,char *level) { + + //Set thresholding + if (!strcmp(thresh, "soft")) { + strcpy(obj->thresh, "soft"); + } + else if (!strcmp(thresh, "hard")) { + strcpy(obj->thresh, "hard"); + } + else { + printf("Thresholding Method - soft or hard"); + exit(-1); + } + + // Set Noise estimation at the first level or at all levels + + if (!strcmp(level, "first")) { + strcpy(obj->level, "first"); + } + else if (!strcmp(level, "all")) { + strcpy(obj->level, "all"); + } + else { + printf("Noise Estimation at level - first or all"); + exit(-1); + } +} + +void denoise_free(denoise_object object) { + free(object); +} diff --git a/auxiliary/denoise.h b/auxiliary/denoise.h index 8ce8cd3..11e0e99 100644 --- a/auxiliary/denoise.h +++ b/auxiliary/denoise.h @@ -11,11 +11,38 @@ Copyright (c) 2017, Rafat Hussain extern "C" { #endif -//depends on J +typedef struct denoise_set* denoise_object; + +denoise_object denoise_init(int length, int J,char* wname); + +struct denoise_set{ + int N; //signal length + int J; // Levels of Wavelet decomposition + char wname[10]; //Wavelet name + char wmethod[10]; //Wavelet decomposition method - dwt or swt + char ext[10]; // Signal Extension - sym or per + char thresh[10]; // thresholding - soft or hard + char level[10]; // Noise Estimation level - first or all + char dmethod[20]; //Denoising Method -sureshrink or visushrink + //double params[0]; +}; + void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); +void denoise(denoise_object obj, double *signal,double *denoised); + +void setDenoiseMethod(denoise_object obj, char *dmethod); + +void setDenoiseWTMethod(denoise_object obj, char *wmethod); + +void setDenoiseWTExtension(denoise_object obj, char *extension); + +void setDenoiseParameters(denoise_object obj, char *thresh,char *level); + +void denoise_free(denoise_object object); + #ifdef __cplusplus } diff --git a/auxiliary/waux.c b/auxiliary/waux.c index b45454c..0626d8b 100644 --- a/auxiliary/waux.c +++ b/auxiliary/waux.c @@ -11,6 +11,35 @@ int compare_double(const void* a, const void* b) } +double mean(double* vec, int N) { + int i; + double m; + m = 0.0; + + for (i = 0; i < N; ++i) { + m+= vec[i]; + } + m = m / N; + return m; +} + +double var(double* vec, int N) { + double v,temp,m; + int i; + v = 0.0; + m = mean(vec,N); + + for (i = 0; i < N; ++i) { + temp = vec[i] - m; + v+= temp*temp; + } + + v = v / N; + + return v; + +} + double median(double *x, int N) { double sigma; @@ -242,3 +271,52 @@ void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, free(out); } + + +void autocovar(double* vec,int N, double* acov,int M) { + double m,temp1,temp2; + int i,t; + m = mean(vec,N); + + if ( M > N) { + M = N-1; + printf("\n Lag is greater than the length N of the input vector. It is automatically set to length N - 1.\n"); + printf("\n The Output Vector only contains N calculated values."); + } else if ( M < 0) { + M = 0; + } + + for(i = 0; i < M;i++) { + acov[i] = 0.0; + for (t = 0; t < N-i;t++) { + temp1 = vec[t] - m; + temp2 = vec[t+i] - m; + acov[i]+= temp1*temp2; + } + acov[i] = acov[i] / N; + + } + + +} + +void autocorr(double* vec,int N,double* acorr, int M) { + double var; + int i; + if (M > N) { + M = N - 1; + printf("\n Lag is greater than the length N of the input vector. It is automatically set to length N - 1.\n"); + printf("\n The Output Vector only contains N calculated values."); + } + else if (M < 0) { + M = 0; + } + autocovar(vec,N,acorr,M); + var = acorr[0]; + acorr[0] = 1.0; + + for(i = 1; i < M; i++) { + acorr[i] = acorr[i]/var; + } + +} diff --git a/auxiliary/waux.h b/auxiliary/waux.h index 093b040..d0a3dbb 100644 --- a/auxiliary/waux.h +++ b/auxiliary/waux.h @@ -19,9 +19,12 @@ extern "C" { #endif - int compare_double(const void* a, const void* b); +double mean(double* vec, int N); + +double var(double* vec, int N); + double median(double *x, int N); double mad(double *x, int N); @@ -35,6 +38,10 @@ 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); +void autocovar(double* vec,int N,double* acov, int M); + +void autocorr(double* vec,int N,double* acorr, int M); + #ifdef __cplusplus } diff --git a/header/denoise.h b/header/denoise.h deleted file mode 100644 index 2c4ccba..0000000 --- a/header/denoise.h +++ /dev/null @@ -1,30 +0,0 @@ -/* -Copyright (c) 2017, Rafat Hussain -*/ -#ifndef DENOISE_H_ -#define DENOISE_H_ - -#include "wavelib.h" - -#ifdef __cplusplus -extern "C" { -#endif - -//depends on J - -void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); - -void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); - -void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, int J,double *lpr, - double *hpr,int lf,int siglength,double *reccoeff); - -double mad(double *x, int N); - - -#ifdef __cplusplus -} -#endif - - -#endif /* DENOISE_H_ */ diff --git a/header/wauxlib.h b/header/wauxlib.h new file mode 100644 index 0000000..85ca443 --- /dev/null +++ b/header/wauxlib.h @@ -0,0 +1,56 @@ +/* +Copyright (c) 2017, Rafat Hussain +*/ +#ifndef WAUXLIB_H_ +#define WAUXLIB_H_ + +#include "wavelib.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct denoise_set* denoise_object; + +denoise_object denoise_init(int length, int J,char* wname); + +struct denoise_set{ + int N; //signal length + int J; // Levels of Wavelet decomposition + char wname[10]; //Wavelet name + char wmethod[10]; //Wavelet decomposition method - dwt or swt + char ext[10]; // Signal Extension - sym or per + char thresh[10]; // thresholding - soft or hard + char level[10]; // Noise Estimation level - first or all + char dmethod[20]; //Denoising Method -sureshrink or visushrink + //double params[0]; +}; + +void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); + +void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); + +void denoise(denoise_object obj, double *signal,double *denoised); + +void setDenoiseMethod(denoise_object obj, char *dmethod); + +void setDenoiseWTMethod(denoise_object obj, char *wmethod); + +void setDenoiseWTExtension(denoise_object obj, char *extension); + +void setDenoiseParameters(denoise_object obj, char *thresh,char *level); + +void denoise_free(denoise_object object); + +void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, int J,double *lpr, + double *hpr,int lf,int siglength,double *reccoeff); + +double mad(double *x, int N); + + +#ifdef __cplusplus +} +#endif + + +#endif /* WAUXLIB_H_ */ diff --git a/test/denoisetest.c b/test/denoisetest.c index c800b2a..8460240 100644 --- a/test/denoisetest.c +++ b/test/denoisetest.c @@ -2,7 +2,8 @@ #include #include #include -#include "../header/denoise.h" + +#include "../header/wauxlib.h" static double rmse(int N,double *x,double *y) { double rms; @@ -49,7 +50,9 @@ int main() { // gcc -Wall -I../header -L../Bin denoisetest.c -o denoise -lwauxlib -lwavelib -lm double *sig,*inp,*oup; int i,N,J; - FILE *ifp,*ofp; + FILE *ifp; + + denoise_object obj; double temp[2400]; char *wname = "db5"; @@ -100,16 +103,23 @@ int main() { for(i = 0; i < N;++i) { inp[i] = temp[i]; } - visushrink(inp,N,J,wname,method,ext,thresh,level,oup); + obj = denoise_init(N,J,wname); + setDenoiseMethod(obj,"visushrink");// sureshrink is also the default. The other option is visushrink + setDenoiseWTMethod(obj,method);// Default is dwt. the other option is swt + setDenoiseWTExtension(obj,ext);// Default is sym. the other option is per + setDenoiseParameters(obj,thresh,level);// Default for thresh is soft. Other option is hard + // Default for level is first. The other option is all + + denoise(obj,inp,oup); + + // Alternative to denoise_object + // Just use visushrink and sureshrink functions + //visushrink(inp,N,J,wname,method,ext,thresh,level,oup); //sureshrink(inp,N,J,wname,method,ext,thresh,level,oup); //ofp = fopen("denoiseds.txt", "w"); - for(i = 0; i < N;++i) { - //fprintf(ofp,"%g \n",oup[i]); - } - //fclose(ofp); printf("RMSE %g\n",rmse(N,sig,inp)); printf("Corr Coeff %g\n",corrcoef(N,sig,inp)); @@ -119,6 +129,6 @@ int main() { free(sig); free(inp); - free(oup); + denoise_free(obj); return 0; }