From a26e271971f98ed9490b0f722ed5ac1a47c71a3f Mon Sep 17 00:00:00 2001 From: Rafat Hussain Date: Sat, 23 Sep 2017 17:19:28 +0530 Subject: [PATCH] commit : visushrink/sureshrink added for dwt/swt --- auxiliary/denoise.c | 139 +++++++++++++++++++++++++++++--------------- auxiliary/denoise.h | 4 +- header/denoise.h | 4 +- test/denoisetest.c | 9 +-- 4 files changed, 102 insertions(+), 54 deletions(-) diff --git a/auxiliary/denoise.c b/auxiliary/denoise.c index 7f5e1eb..ae9a7a5 100644 --- a/auxiliary/denoise.c +++ b/auxiliary/denoise.c @@ -1,12 +1,12 @@ #include "denoise.h" -void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,double *denoised) { - int filt_len,iter,i,dlen,dwt_len,sgn, MaxIter; +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; wave_object wave; wt_object wt; - double *dout; + double *dout,*lnoise; wave = wave_init(wname); @@ -29,8 +29,11 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch modwt(wt,signal); } else { printf("Acceptable WT methods are - dwt,swt and modwt\n"); + exit(-1); } + lnoise = (double*)malloc(sizeof(double) * J); + //Set sigma iter = wt->length[0]; @@ -38,36 +41,61 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch dout = (double*)malloc(sizeof(double) * dlen); - 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); - dwt_len = wt->outlength; - - td = sqrt(2.0 * log(dwt_len)) * sigma / 0.6745; - - - if(!strcmp(thresh,"hard")) { - for(i = wt->length[0]; i < dwt_len;++i) { - if (fabs(wt->output[i]) < td) { - wt->output[i] = 0; - } + if(!strcmp(level,"first")) { + for (i = 1; i < J; ++i) { + iter += wt->length[i]; } - } else if(!strcmp(thresh,"soft")) { - for(i = wt->length[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; + + 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) { + lnoise[it] = sigma; + } + } else if(!strcmp(level,"all")){ + for(it = 0; it < J;++it) { + dlen = wt->length[it+1]; + for(i = 0; i < dlen;++i) { + dout[i] = fabs(wt->output[iter+i]); + } + sigma = median(dout,dlen) / 0.6745; + lnoise[it] = sigma; + iter += dlen; + } + + } else { + printf("Acceptable Noise estimation level values are - first and all \n"); + exit(-1); + } + + dwt_len = wt->outlength; + iter = wt->length[0]; + for(it = 0; it < J;++it) { + sigma = lnoise[it]; + dlen = wt->length[it+1]; + td = sqrt(2.0 * log(dwt_len)) * sigma; + + if(!strcmp(thresh,"hard")) { + for(i = 0; i < dlen;++i) { + if (fabs(wt->output[iter+i]) < td) { + wt->output[iter+i] = 0; } } + } else if(!strcmp(thresh,"soft")) { + for(i = 0; i < dlen;++i) { + if (fabs(wt->output[iter + i]) < td) { + wt->output[iter+i] = 0; + } else { + sgn = wt->output[iter+i] >= 0 ? 1 : -1; + tmp = sgn * (fabs(wt->output[iter+i]) - td); + wt->output[iter+i] = tmp; + } + } + } + + iter += wt->length[it+1]; } if(!strcmp(method,"dwt")) { @@ -79,16 +107,17 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch } free(dout); + free(lnoise); wave_free(wave); wt_free(wt); } -void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,double *denoised) { +void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised) { 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,*dsum; + double *dout,*risk,*dsum,*lnoise; wave = wave_init(wname); @@ -108,10 +137,9 @@ void sureshrink(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"); + printf("Acceptable WT methods are - dwt and swt\n"); + exit(-1); } len = wt->length[0]; @@ -120,22 +148,43 @@ 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); + lnoise = (double*)malloc(sizeof(double) * J); iter = wt->length[0]; - for (i = 1; i < J; ++i) { - iter += wt->length[i]; + if(!strcmp(level,"first")) { + 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) { + lnoise[it] = sigma; + } + } else if(!strcmp(level,"all")){ + for(it = 0; it < J;++it) { + dlen = wt->length[it+1]; + for(i = 0; i < dlen;++i) { + dout[i] = fabs(wt->output[iter+i]); + } + sigma = median(dout,dlen) / 0.6745; + lnoise[it] = sigma; + iter += dlen; + } + + } else { + printf("Acceptable Noise estimation level values are - first and all \n"); + exit(-1); } - 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]; - + sigma = lnoise[it]; if ( sigma < 0.00000001) { td = 0; @@ -202,13 +251,11 @@ void sureshrink(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); free(dsum); free(risk); + free(lnoise); wave_free(wave); wt_free(wt); } diff --git a/auxiliary/denoise.h b/auxiliary/denoise.h index 004a84d..8ce8cd3 100644 --- a/auxiliary/denoise.h +++ b/auxiliary/denoise.h @@ -12,9 +12,9 @@ extern "C" { #endif //depends on J -void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,double *denoised); +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,double *denoised); +void sureshrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); #ifdef __cplusplus diff --git a/header/denoise.h b/header/denoise.h index 113c7e0..2c4ccba 100644 --- a/header/denoise.h +++ b/header/denoise.h @@ -12,9 +12,9 @@ extern "C" { //depends on J -void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,double *denoised); +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,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); diff --git a/test/denoisetest.c b/test/denoisetest.c index 3a15f3c..c800b2a 100644 --- a/test/denoisetest.c +++ b/test/denoisetest.c @@ -52,10 +52,11 @@ int main() { FILE *ifp,*ofp; double temp[2400]; - char *wname = "sym6"; + char *wname = "db5"; char *method = "dwt"; char *ext = "sym"; char *thresh = "soft"; + char *level = "all"; ifp = fopen("pieceregular1024.txt", "r"); i = 0; @@ -72,7 +73,7 @@ int main() { fclose(ifp); N = i; - J = 6; + J = 4; sig = (double*)malloc(sizeof(double)* N); inp = (double*)malloc(sizeof(double)* N); @@ -99,8 +100,8 @@ int main() { for(i = 0; i < N;++i) { inp[i] = temp[i]; } - //visushrink(inp,N,J,wname,method,ext,thresh,oup); - sureshrink(inp,N,J,wname,method,ext,thresh,oup); + 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");