mirror of
https://github.com/simon987/wavelib.git
synced 2025-11-04 02:16:53 +00:00
490 lines
12 KiB
C
490 lines
12 KiB
C
#include <math.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
#include "waux.h"
|
|
#include "../header/wauxlib.h"
|
|
|
|
denoise_object denoise_init(int length, int J,const 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,"all");
|
|
strcpy(obj->thresh,"soft");
|
|
strcpy(obj->wmethod,"dwt");
|
|
strcpy(obj->cmethod,"direct");
|
|
|
|
return obj;
|
|
}
|
|
|
|
void visushrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const 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,*lnoise;
|
|
|
|
wave = wave_init(wname);
|
|
|
|
filt_len = wave->filtlength;
|
|
|
|
MaxIter = (int) (log((double)N / ((double)filt_len - 1.0)) / log(2.0));
|
|
|
|
if (J > MaxIter) {
|
|
printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter);
|
|
exit(-1);
|
|
}
|
|
|
|
wt = wt_init(wave,method,N,J);
|
|
if(!strcmp(method,"dwt")) {
|
|
setDWTExtension(wt,ext);
|
|
dwt(wt,signal);
|
|
} else if(!strcmp(method,"swt")) {
|
|
swt(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];
|
|
dlen = wt->length[J];
|
|
|
|
dout = (double*)malloc(sizeof(double) * dlen);
|
|
|
|
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);
|
|
}
|
|
|
|
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")) {
|
|
idwt(wt,denoised);
|
|
} else if(!strcmp(method,"swt")) {
|
|
iswt(wt,denoised);
|
|
}
|
|
|
|
free(dout);
|
|
free(lnoise);
|
|
wave_free(wave);
|
|
wt_free(wt);
|
|
}
|
|
|
|
void sureshrink(double *signal,int N,int J,const char *wname,const char *method,const char *ext,const char *thresh,const 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,*lnoise;
|
|
|
|
wave = wave_init(wname);
|
|
|
|
filt_len = wave->filtlength;
|
|
|
|
MaxIter = (int) (log((double)N / ((double)filt_len - 1.0)) / log(2.0));
|
|
// Depends on J
|
|
if (J > MaxIter) {
|
|
printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter);
|
|
exit(-1);
|
|
}
|
|
|
|
wt = wt_init(wave,method,N,J);
|
|
|
|
if(!strcmp(method,"dwt")) {
|
|
setDWTExtension(wt,ext);
|
|
dwt(wt,signal);
|
|
} else if(!strcmp(method,"swt")) {
|
|
swt(wt,signal);
|
|
} else {
|
|
printf("Acceptable WT methods are - dwt and swt\n");
|
|
exit(-1);
|
|
}
|
|
|
|
len = wt->length[0];
|
|
dlen = wt->length[J];
|
|
|
|
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];
|
|
|
|
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(it = 0; it < J;++it) {
|
|
dwt_len = wt->length[it+1];
|
|
sigma = lnoise[it];
|
|
|
|
if ( sigma < 0.00000001) {
|
|
td = 0;
|
|
} else {
|
|
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] /(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);
|
|
|
|
if (te < ct) {
|
|
td = tv;
|
|
} else {
|
|
x_sum = 0.0;
|
|
|
|
for(i = 0; i < dwt_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 - 2 * ((double)i + 1) +dsum[i] +
|
|
dout[i]*((double)dwt_len - 1 -(double) i))/(double)dwt_len;
|
|
}
|
|
min_index = minindex(risk,dwt_len);
|
|
thr = sqrt(dout[min_index]);
|
|
td = thr < tv ? thr : tv;
|
|
}
|
|
}
|
|
|
|
td = td * sigma;
|
|
|
|
if(!strcmp(thresh,"hard")) {
|
|
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 = 0; i < dwt_len;++i) {
|
|
if (fabs(wt->output[len + i]) < td) {
|
|
wt->output[len+i] = 0;
|
|
} else {
|
|
sgn = wt->output[len+i] >= 0 ? 1 : -1;
|
|
temp = sgn * (fabs(wt->output[len+i]) - td);
|
|
wt->output[len+i] = temp;
|
|
}
|
|
}
|
|
}
|
|
|
|
len += wt->length[it+1];
|
|
}
|
|
|
|
if(!strcmp(method,"dwt")) {
|
|
idwt(wt,denoised);
|
|
} else if(!strcmp(method,"swt")) {
|
|
iswt(wt,denoised);
|
|
}
|
|
free(dout);
|
|
free(dsum);
|
|
free(risk);
|
|
free(lnoise);
|
|
wave_free(wave);
|
|
wt_free(wt);
|
|
}
|
|
|
|
void modwtshrink(double *signal, int N, int J, const char *wname, const char *cmethod, const char *ext, const char *thresh, double *denoised) {
|
|
int filt_len, iter, i, dlen, sgn, MaxIter, it;
|
|
double sigma, td, tmp, M, llen;
|
|
wave_object wave;
|
|
wt_object wt;
|
|
double *dout, *lnoise;
|
|
|
|
wave = wave_init(wname);
|
|
|
|
filt_len = wave->filtlength;
|
|
|
|
MaxIter = (int)(log((double)N / ((double)filt_len - 1.0)) / log(2.0));
|
|
|
|
if (J > MaxIter) {
|
|
printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
|
|
exit(-1);
|
|
}
|
|
|
|
wt = wt_init(wave, "modwt", N, J);
|
|
|
|
if (!strcmp(ext, "sym") && !strcmp(cmethod,"fft")) {
|
|
setWTConv(wt, "fft");
|
|
setDWTExtension(wt, "sym");
|
|
}
|
|
else if (!strcmp(ext, "sym") && !strcmp(cmethod, "direct")) {
|
|
printf("Symmetric Extension is not available for direct method");
|
|
exit(-1);
|
|
}
|
|
else if (!strcmp(ext, "per") && !strcmp(cmethod, "direct")) {
|
|
setWTConv(wt, "direct");
|
|
setDWTExtension(wt, "per");
|
|
}
|
|
else if (!strcmp(ext, "per") && !strcmp(cmethod, "fft")) {
|
|
setWTConv(wt, "fft");
|
|
setDWTExtension(wt, "per");
|
|
}
|
|
else {
|
|
printf("Signal extension can be either per or sym");
|
|
exit(-1);
|
|
}
|
|
|
|
modwt(wt, signal);
|
|
|
|
lnoise = (double*)malloc(sizeof(double)* J);
|
|
|
|
//Set sigma
|
|
|
|
iter = wt->length[0];
|
|
dlen = wt->length[J];
|
|
dout = (double*)malloc(sizeof(double)* dlen);
|
|
|
|
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 = sqrt(2.0) * median(dout, dlen) / 0.6745;
|
|
lnoise[it] = sigma;
|
|
iter += dlen;
|
|
}
|
|
|
|
M = pow(2.0,J);
|
|
llen = log((double)wt->modwtsiglength);
|
|
// Thresholding
|
|
|
|
iter = wt->length[0];
|
|
for (it = 0; it < J; ++it) {
|
|
sigma = lnoise[it];
|
|
dlen = wt->length[it + 1];
|
|
td = sqrt(2.0 * llen / M) * 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];
|
|
M /= 2.0;
|
|
}
|
|
|
|
imodwt(wt, denoised);
|
|
|
|
free(dout);
|
|
free(lnoise);
|
|
wave_free(wave);
|
|
wt_free(wt);
|
|
}
|
|
|
|
|
|
void denoise(denoise_object obj, double *signal,double *denoised) {
|
|
if(!strcmp(obj->dmethod,"sureshrink")) {
|
|
if (!strcmp(obj->wmethod, "modwt")) {
|
|
printf("sureshrink method only works with swt and dwt. Please use setDenoiseWTMethod to set the correct method\n");
|
|
exit(-1);
|
|
}
|
|
sureshrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);
|
|
} else if(!strcmp(obj->dmethod,"visushrink")) {
|
|
if (!strcmp(obj->wmethod, "modwt")) {
|
|
printf("visushrink method only works with swt and dwt. Please use setDenoiseWTMethod to set the correct method\n");
|
|
exit(-1);
|
|
}
|
|
visushrink(signal,obj->N,obj->J,obj->wname,obj->wmethod,obj->ext,obj->thresh,obj->level,denoised);;
|
|
} else if(!strcmp(obj->dmethod,"modwtshrink")) {
|
|
if (strcmp(obj->wmethod, "modwt")) {
|
|
printf("modwtshrink method only works with modwt. Please use setDenoiseWTMethod to set the correct method\n");
|
|
exit(-1);
|
|
}
|
|
modwtshrink(signal,obj->N,obj->J,obj->wname,obj->cmethod,obj->ext,obj->thresh,denoised);;
|
|
} else {
|
|
printf("Acceptable Denoising methods are - sureshrink and visushrink\n");
|
|
exit(-1);
|
|
}
|
|
}
|
|
|
|
void setDenoiseMethod(denoise_object obj, const char *dmethod) {
|
|
if (!strcmp(dmethod, "sureshrink")) {
|
|
strcpy(obj->dmethod, "sureshrink");
|
|
}
|
|
else if (!strcmp(dmethod, "visushrink")) {
|
|
strcpy(obj->dmethod, "visushrink");
|
|
}
|
|
else if (!strcmp(dmethod, "modwtshrink")) {
|
|
strcpy(obj->dmethod, "modwtshrink");
|
|
}
|
|
else {
|
|
printf("Acceptable Denoising methods are - sureshrink, visushrink and modwtshrink\n");
|
|
exit(-1);
|
|
}
|
|
}
|
|
|
|
void setDenoiseWTMethod(denoise_object obj, const char *wmethod) {
|
|
if (!strcmp(wmethod, "dwt")) {
|
|
strcpy(obj->wmethod, "dwt");
|
|
}
|
|
else if (!strcmp(wmethod, "swt")) {
|
|
strcpy(obj->wmethod, "swt");
|
|
}
|
|
else if (!strcmp(wmethod, "modwt")) {
|
|
strcpy(obj->wmethod, "modwt");
|
|
}
|
|
else {
|
|
printf("Wavelet decomposition method can be one of dwt, modwt or swt.\n");
|
|
exit(-1);
|
|
}
|
|
}
|
|
|
|
void setDenoiseWTExtension(denoise_object obj, const 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, const char *thresh,const 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);
|
|
}
|