simon987/wavelib
1
0
Fork 0
mirror of https://github.com/simon987/wavelib.git synced 2025-04-10 14:06:46 +00:00
Code Issues Projects Releases Wiki Activity

commit : sureshrink corrected

Browse Source
This commit is contained in:
Rafat Hussain 2017-09-22 17:01:13 +05:30
parent 310b7759f1
commit bd2314727e
1 changed files with 24 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
auxiliary/denoise.c
View File

@ -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);