Merge branch 'aitap-const-correct'

Const-correctness, unused parameter cleanup, build system improvements
This commit is contained in:
Rafat Hussain 2018-12-14 16:38:13 +05:30
commit 924d3a3c9b
20 changed files with 402 additions and 739 deletions

View File

@ -4,9 +4,7 @@ set(SOURCE_FILES denoise.c
waux.c waux.c
) )
set(HEADER_FILES denoise.h set(HEADER_FILES waux.h)
waux.h
)
add_library(wauxlib STATIC ${SOURCE_FILES} ${HEADER_FILES}) add_library(wauxlib STATIC ${SOURCE_FILES} ${HEADER_FILES})

View File

@ -1,7 +1,12 @@
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "denoise.h" #include "waux.h"
#include "wauxlib.h"
denoise_object denoise_init(int length, int J,char* wname) { denoise_object denoise_init(int length, int J,const char* wname) {
denoise_object obj = NULL; denoise_object obj = NULL;
obj = (denoise_object)malloc(sizeof(struct denoise_set) +sizeof(double)); obj = (denoise_object)malloc(sizeof(struct denoise_set) +sizeof(double));
@ -21,7 +26,7 @@ denoise_object denoise_init(int length, int J,char* wname) {
return obj; return obj;
} }
void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised) { 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; int filt_len,iter,i,dlen,dwt_len,sgn, MaxIter,it;
double sigma,td,tmp; double sigma,td,tmp;
wave_object wave; wave_object wave;
@ -128,7 +133,7 @@ void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,ch
wt_free(wt); wt_free(wt);
} }
void sureshrink(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,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; 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; double sigma,norm,td,tv,te,ct,thr,temp,x_sum;
wave_object wave; wave_object wave;
@ -286,7 +291,7 @@ void denoise(denoise_object obj, double *signal,double *denoised) {
} }
} }
void setDenoiseMethod(denoise_object obj, char *dmethod) { void setDenoiseMethod(denoise_object obj, const char *dmethod) {
if (!strcmp(dmethod, "sureshrink")) { if (!strcmp(dmethod, "sureshrink")) {
strcpy(obj->dmethod, "sureshrink"); strcpy(obj->dmethod, "sureshrink");
} }
@ -299,7 +304,7 @@ void setDenoiseMethod(denoise_object obj, char *dmethod) {
} }
} }
void setDenoiseWTMethod(denoise_object obj, char *wmethod) { void setDenoiseWTMethod(denoise_object obj, const char *wmethod) {
if (!strcmp(wmethod, "dwt")) { if (!strcmp(wmethod, "dwt")) {
strcpy(obj->wmethod, "dwt"); strcpy(obj->wmethod, "dwt");
} }
@ -312,7 +317,7 @@ void setDenoiseWTMethod(denoise_object obj, char *wmethod) {
} }
} }
void setDenoiseWTExtension(denoise_object obj, char *extension) { void setDenoiseWTExtension(denoise_object obj, const char *extension) {
if (!strcmp(extension, "sym")) { if (!strcmp(extension, "sym")) {
strcpy(obj->ext, "sym"); strcpy(obj->ext, "sym");
} }
@ -325,7 +330,7 @@ void setDenoiseWTExtension(denoise_object obj, char *extension) {
} }
} }
void setDenoiseParameters(denoise_object obj, char *thresh,char *level) { void setDenoiseParameters(denoise_object obj, const char *thresh,const char *level) {
//Set thresholding //Set thresholding
if (!strcmp(thresh, "soft")) { if (!strcmp(thresh, "soft")) {

View File

@ -1,52 +0,0 @@
/*
Copyright (c) 2017, Rafat Hussain
*/
#ifndef DENOISE_H_
#define DENOISE_H_
#include "waux.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);
#ifdef __cplusplus
}
#endif
#endif /* DENOISE_H_ */

View File

@ -1,3 +1,4 @@
#include "wauxlib.h"
#include "waux.h" #include "waux.h"
int compare_double(const void* a, const void* b) int compare_double(const void* a, const void* b)
@ -11,7 +12,7 @@ int compare_double(const void* a, const void* b)
} }
double mean(double* vec, int N) { double mean(const double* vec, int N) {
int i; int i;
double m; double m;
m = 0.0; m = 0.0;
@ -23,7 +24,7 @@ double mean(double* vec, int N) {
return m; return m;
} }
double var(double* vec, int N) { double var(const double* vec, int N) {
double v,temp,m; double v,temp,m;
int i; int i;
v = 0.0; v = 0.0;
@ -69,7 +70,7 @@ double mad(double *x, int N) {
return sigma; return sigma;
} }
int minindex(double *arr, int N) { int minindex(const double *arr, int N) {
double min; double min;
int index,i; int index,i;
@ -130,7 +131,7 @@ 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, void getDWTRecCoeff(double *coeff,int *length,const char *ctype,const char *ext, int level, int J,double *lpr,
double *hpr,int lf,int siglength,double *reccoeff) { double *hpr,int lf,int siglength,double *reccoeff) {
int i,j,k,det_len,N,l,m,n,v,t,l2; int i,j,k,det_len,N,l,m,n,v,t,l2;
@ -274,7 +275,7 @@ void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level,
} }
void autocovar(double* vec,int N, double* acov,int M) { void autocovar(const double* vec,int N, double* acov,int M) {
double m,temp1,temp2; double m,temp1,temp2;
int i,t; int i,t;
m = mean(vec,N); m = mean(vec,N);
@ -301,7 +302,7 @@ void autocovar(double* vec,int N, double* acov,int M) {
} }
void autocorr(double* vec,int N,double* acorr, int M) { void autocorr(const double* vec,int N,double* acorr, int M) {
double var; double var;
int i; int i;
if (M > N) { if (M > N) {

View File

@ -21,26 +21,21 @@ extern "C" {
int compare_double(const void* a, const void* b); int compare_double(const void* a, const void* b);
double mean(double* vec, int N); double mean(const double* vec, int N);
double var(double* vec, int N); double var(const double* vec, int N);
double median(double *x, int N); double median(double *x, int N);
double mad(double *x, int N); int minindex(const double *arr, int N);
int minindex(double *arr, int N);
void getDWTAppx(wt_object wt, double *appx,int N); void getDWTAppx(wt_object wt, double *appx,int N);
void getDWTDetail(wt_object wt, double *detail, int N, int level); 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, void autocovar(const double* vec,int N,double* acov, int M);
double *hpr,int lf,int siglength,double *reccoeff);
void autocovar(double* vec,int N,double* acov, int M); void autocorr(const double* vec,int N,double* acorr, int M);
void autocorr(double* vec,int N,double* acorr, int M);
#ifdef __cplusplus #ifdef __cplusplus

View File

@ -12,7 +12,7 @@ extern "C" {
typedef struct denoise_set* denoise_object; typedef struct denoise_set* denoise_object;
denoise_object denoise_init(int length, int J,char* wname); denoise_object denoise_init(int length, int J,const char* wname);
struct denoise_set{ struct denoise_set{
int N; //signal length int N; //signal length
@ -26,23 +26,23 @@ struct denoise_set{
//double params[0]; //double params[0];
}; };
void visushrink(double *signal,int N,int J,char *wname,char *method,char *ext,char *thresh,char *level,double *denoised); 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);
void sureshrink(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,const char *wname,const char *method,const char *ext,const char *thresh,const char *level,double *denoised);
void denoise(denoise_object obj, double *signal,double *denoised); void denoise(denoise_object obj, double *signal,double *denoised);
void setDenoiseMethod(denoise_object obj, char *dmethod); void setDenoiseMethod(denoise_object obj, const char *dmethod);
void setDenoiseWTMethod(denoise_object obj, char *wmethod); void setDenoiseWTMethod(denoise_object obj, const char *wmethod);
void setDenoiseWTExtension(denoise_object obj, char *extension); void setDenoiseWTExtension(denoise_object obj, const char *extension);
void setDenoiseParameters(denoise_object obj, char *thresh,char *level); void setDenoiseParameters(denoise_object obj, const char *thresh,const char *level);
void denoise_free(denoise_object object); void denoise_free(denoise_object object);
void getDWTRecCoeff(double *coeff,int *length,char *ctype,char *ext, int level, int J,double *lpr, void getDWTRecCoeff(double *coeff,int *length,const char *ctype,const char *ext, int level, int J,double *lpr,
double *hpr,int lf,int siglength,double *reccoeff); double *hpr,int lf,int siglength,double *reccoeff);
double mad(double *x, int N); double mad(double *x, int N);

View File

@ -26,7 +26,7 @@ typedef struct cplx_t {
typedef struct wave_set* wave_object; typedef struct wave_set* wave_object;
wave_object wave_init(char* wname); wave_object wave_init(const char* wname);
struct wave_set{ struct wave_set{
char wname[50]; char wname[50];
@ -83,7 +83,7 @@ struct conv_set{
typedef struct wt_set* wt_object; typedef struct wt_set* wt_object;
wt_object wt_init(wave_object wave,char* method, int siglength, int J); wt_object wt_init(wave_object wave,const char* method, int siglength, int J);
struct wt_set{ struct wt_set{
wave_object wave; wave_object wave;
@ -165,7 +165,7 @@ struct wpt_set{
typedef struct cwt_set* cwt_object; typedef struct cwt_set* cwt_object;
cwt_object cwt_init(char* wave, double param, int siglength,double dt, int J); cwt_object cwt_init(const char* wave, double param, int siglength,double dt, int J);
struct cwt_set{ struct cwt_set{
char wave[10];// Wavelet - morl/morlet,paul,dog/dgauss char wave[10];// Wavelet - morl/morlet,paul,dog/dgauss
@ -191,33 +191,33 @@ struct cwt_set{
}; };
void dwt(wt_object wt, double *inp); void dwt(wt_object wt, const double *inp);
void idwt(wt_object wt, double *dwtop); void idwt(wt_object wt, double *dwtop);
void wtree(wtree_object wt, double *inp); void wtree(wtree_object wt, const double *inp);
void dwpt(wpt_object wt, double *inp); void dwpt(wpt_object wt, const double *inp);
void idwpt(wpt_object wt, double *dwtop); void idwpt(wpt_object wt, double *dwtop);
void swt(wt_object wt, double *inp); void swt(wt_object wt, const double *inp);
void iswt(wt_object wt, double *swtop); void iswt(wt_object wt, double *swtop);
void modwt(wt_object wt, double *inp); void modwt(wt_object wt, const double *inp);
void imodwt(wt_object wt, double *dwtop); void imodwt(wt_object wt, double *dwtop);
void setDWTExtension(wt_object wt, char *extension); void setDWTExtension(wt_object wt, const char *extension);
void setWTREEExtension(wtree_object wt, char *extension); void setWTREEExtension(wtree_object wt, const char *extension);
void setDWPTExtension(wpt_object wt, char *extension); void setDWPTExtension(wpt_object wt, const char *extension);
void setDWPTEntropy(wpt_object wt, char *entropy, double eparam); void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam);
void setWTConv(wt_object wt, char *cmethod); void setWTConv(wt_object wt, const char *cmethod);
int getWTREENodelength(wtree_object wt, int X); int getWTREENodelength(wtree_object wt, int X);
@ -227,13 +227,13 @@ int getDWPTNodelength(wpt_object wt, int X);
void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N); void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N);
void setCWTScales(cwt_object wt, double s0, double dj, char *type, int power); void setCWTScales(cwt_object wt, double s0, double dj, const char *type, int power);
void setCWTScaleVector(cwt_object wt, double *scale, int J, double s0, double dj); void setCWTScaleVector(cwt_object wt, const double *scale, int J, double s0, double dj);
void setCWTPadding(cwt_object wt, int pad); void setCWTPadding(cwt_object wt, int pad);
void cwt(cwt_object wt, double *inp); void cwt(cwt_object wt, const double *inp);
void icwt(cwt_object wt, double *cwtop); void icwt(cwt_object wt, double *cwtop);

View File

@ -20,7 +20,6 @@ set(HEADER_FILES conv.h
real.h real.h
wavefilt.h wavefilt.h
wavefunc.h wavefunc.h
wavelib.h
wtmath.h wtmath.h
) )
@ -28,5 +27,5 @@ add_library(wavelib STATIC ${SOURCE_FILES} ${HEADER_FILES})
set_property(TARGET wavelib PROPERTY FOLDER "lib") set_property(TARGET wavelib PROPERTY FOLDER "lib")
target_include_directories(wavelib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) target_include_directories(wavelib PUBLIC ${CMAKE_SOURCE_DIR}/header)

View File

@ -17,18 +17,8 @@ extern "C" {
#define MIN(a,b) (((a)<(b))?(a):(b)) #define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b)) #define MAX(a,b) (((a)>(b))?(a):(b))
typedef struct conv_set* conv_object;
conv_object conv_init(int N, int L); conv_object conv_init(int N, int L);
struct conv_set{
fft_real_object fobj;
fft_real_object iobj;
int ilen1;
int ilen2;
int clen;
};
int factorf(int M); int factorf(int M);
int findnext(int M); int findnext(int M);

View File

@ -8,37 +8,8 @@ C. Torrence and G. Compo, and is available at URL: http://atoc.colorado.edu/rese
#include "cwt.h" #include "cwt.h"
static int factorial2(int N) {
int factorial,i;
factorial = 1;
for (i = 1; i <= N;++i) {
factorial *= i;
}
return factorial;
}
static double factorial3(int N) {
int i;
double factorial;
factorial = 1;
for (i = 1; i <= N; ++i) {
factorial *= i;
}
return factorial;
}
double factorial(int N) { double factorial(int N) {
if (N > 40) { static const double fact[41] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000,
printf("This program is only valid for N <= 40 \n");
return -1.0;
}
double fact[41] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000,
20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000, 51090942171709440000.0, 1124000727777607680000.0, 20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000, 51090942171709440000.0, 1124000727777607680000.0,
25852016738884976640000.0, 620448401733239439360000.0, 15511210043330985984000000.0, 403291461126605635584000000.0, 10888869450418352160768000000.0, 25852016738884976640000.0, 620448401733239439360000.0, 15511210043330985984000000.0, 403291461126605635584000000.0, 10888869450418352160768000000.0,
304888344611713860501504000000.0, 8841761993739701954543616000000.0, 265252859812191058636308480000000.0, 8222838654177922817725562880000000.0, 304888344611713860501504000000.0, 8841761993739701954543616000000.0, 265252859812191058636308480000000.0, 8222838654177922817725562880000000.0,
@ -46,6 +17,11 @@ double factorial(int N) {
371993326789901217467999448150835200000000.0, 13763753091226345046315979581580902400000000.0, 523022617466601111760007224100074291200000000.0, 371993326789901217467999448150835200000000.0, 13763753091226345046315979581580902400000000.0, 523022617466601111760007224100074291200000000.0,
20397882081197443358640281739902897356800000000.0, 815915283247897734345611269596115894272000000000.0 }; 20397882081197443358640281739902897356800000000.0, 815915283247897734345611269596115894272000000000.0 };
if (N > 40 || N < 0) {
printf("This program is only valid for 0 <= N <= 40 \n");
return -1.0;
}
return fact[N]; return fact[N];
} }
@ -142,7 +118,7 @@ static void wave_function(int nk, double dt,int mother, double param,double scal
} }
} }
void cwavelet(double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad, void cwavelet(const double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad,
double *wave, double *scale, double *period, double *coi) { double *wave, double *scale, double *period, double *coi) {
int i, j, k, iter; int i, j, k, iter;
@ -153,6 +129,8 @@ void cwavelet(double *y, int N, double dt, int mother, double param, double s0,
fft_object obj, iobj; fft_object obj, iobj;
fft_data *ypad, *yfft,*daughter; fft_data *ypad, *yfft,*daughter;
(void)s0; (void)dj; /* yes, we need these parameters unused */
pi = 4.0 * atan(1.0); pi = 4.0 * atan(1.0);
if (npad < N) { if (npad < N) {

View File

@ -7,7 +7,7 @@
extern "C" { extern "C" {
#endif #endif
void cwavelet(double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad, void cwavelet(const double *y, int N, double dt, int mother, double param, double s0, double dj, int jtot, int npad,
double *wave, double *scale, double *period, double *coi); double *wave, double *scale, double *period, double *coi);
void psi0(int mother, double param, double *val, int *real); void psi0(int mother, double param, double *val, int *real);

View File

@ -18,7 +18,7 @@ fft_object fft_init(int N, int sgn) {
if (out == 1) { if (out == 1) {
obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (N-1)); obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (N-1));
obj->lf = factors(N,obj->factors); obj->lf = factors(N,obj->factors);
longvectorN(obj->twiddle,N,obj->factors,obj->lf); longvectorN(obj->twiddle,obj->factors,obj->lf);
twi_len = N; twi_len = N;
obj->lt = 0; obj->lt = 0;
} else { } else {
@ -32,7 +32,7 @@ fft_object fft_init(int N, int sgn) {
} }
obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (M-1)); obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (M-1));
obj->lf = factors(M,obj->factors); obj->lf = factors(M,obj->factors);
longvectorN(obj->twiddle,M,obj->factors,obj->lf); longvectorN(obj->twiddle,obj->factors,obj->lf);
obj->lt = 1; obj->lt = 1;
twi_len = M; twi_len = M;
} }
@ -1831,7 +1831,7 @@ void twiddle(fft_data *vec,int N, int radix) {
} }
void longvectorN(fft_data *sig,int N, int *array, int tx) { void longvectorN(fft_data *sig, int *array, int tx) {
int L,i,Ls,ct,j,k; int L,i,Ls,ct,j,k;
fft_type theta; fft_type theta;
L = 1; L = 1;

View File

@ -12,6 +12,8 @@
#include <math.h> #include <math.h>
#include <string.h> #include <string.h>
#include "wavelib.h"
#ifdef __cplusplus #ifdef __cplusplus
extern "C" { extern "C" {
#endif #endif
@ -22,11 +24,6 @@ extern "C" {
#define fft_type double #define fft_type double
#endif #endif
typedef struct fft_t {
fft_type re;
fft_type im;
} fft_data;
/* /*
#define SADD(a,b) ((a)+(b)) #define SADD(a,b) ((a)+(b))
@ -35,19 +32,8 @@ typedef struct fft_t {
#define SMUL(a,b) ((a)*(b)) #define SMUL(a,b) ((a)*(b))
*/ */
typedef struct fft_set* fft_object;
fft_object fft_init(int N, int sgn); fft_object fft_init(int N, int sgn);
struct fft_set{
int N;
int sgn;
int factors[64];
int lf;
int lt;
fft_data twiddle[1];
};
void fft_exec(fft_object obj,fft_data *inp,fft_data *oup); void fft_exec(fft_object obj,fft_data *inp,fft_data *oup);
int divideby(int M,int d); int divideby(int M,int d);
@ -60,7 +46,7 @@ int factors(int M, int* arr);
void twiddle(fft_data *sig,int N, int radix); void twiddle(fft_data *sig,int N, int radix);
void longvectorN(fft_data *sig,int N, int *array, int M); void longvectorN(fft_data *sig, int *array, int M);
void free_fft(fft_object object); void free_fft(fft_object object);

View File

@ -9,11 +9,9 @@
fft_real_object fft_real_init(int N, int sgn) { fft_real_object fft_real_init(int N, int sgn) {
fft_real_object obj = NULL; fft_real_object obj = NULL;
fft_type PI, theta; fft_type theta;
int k; int k;
PI = 3.1415926535897932384626433832795;
obj = (fft_real_object) malloc (sizeof(struct fft_real_set) + sizeof(fft_data)* (N/2)); obj = (fft_real_object) malloc (sizeof(struct fft_real_set) + sizeof(fft_data)* (N/2));
obj->cobj = fft_init(N/2,sgn); obj->cobj = fft_init(N/2,sgn);
@ -79,10 +77,9 @@ void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup) {
fft_data* cinp; fft_data* cinp;
fft_data* coup; fft_data* coup;
int i,N2,N; int i,N2;
fft_type temp1,temp2; fft_type temp1,temp2;
N2 = obj->cobj->N; N2 = obj->cobj->N;
N = N2*2;
cinp = (fft_data*) malloc (sizeof(fft_data) * N2); cinp = (fft_data*) malloc (sizeof(fft_data) * N2);
coup = (fft_data*) malloc (sizeof(fft_data) * N2); coup = (fft_data*) malloc (sizeof(fft_data) * N2);

View File

@ -14,15 +14,8 @@
extern "C" { extern "C" {
#endif #endif
typedef struct fft_real_set* fft_real_object;
fft_real_object fft_real_init(int N, int sgn); fft_real_object fft_real_init(int N, int sgn);
struct fft_real_set{
fft_object cobj;
fft_data twiddle2[1];
};
void fft_r2c_exec(fft_real_object obj,fft_type *inp,fft_data *oup); void fft_r2c_exec(fft_real_object obj,fft_type *inp,fft_data *oup);
void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup); void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup);

View File

@ -3311,7 +3311,6 @@ void copy_reverse(const double *in, int N,double *out)
void qmf_wrev(const double *in, int N, double *out) void qmf_wrev(const double *in, int N, double *out)
{ {
int count = 0;
double *sigOutTemp; double *sigOutTemp;
sigOutTemp = (double*)malloc(N*sizeof(double)); sigOutTemp = (double*)malloc(N*sizeof(double));

View File

@ -2,9 +2,16 @@
Copyright (c) 2014, Rafat Hussain Copyright (c) 2014, Rafat Hussain
*/ */
#include "wavelib.h" #include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
wave_object wave_init(char* wname) { #include "cwt.h"
#include "wavelib.h"
#include "wtmath.h"
wave_object wave_init(const char* wname) {
wave_object obj = NULL; wave_object obj = NULL;
int retval; int retval;
retval = 0; retval = 0;
@ -30,7 +37,7 @@ wave_object wave_init(char* wname) {
return obj; return obj;
} }
wt_object wt_init(wave_object wave,char* method, int siglength,int J) { wt_object wt_init(wave_object wave,const char* method, int siglength,int J) {
int size,i,MaxIter; int size,i,MaxIter;
wt_object obj = NULL; wt_object obj = NULL;
@ -258,13 +265,13 @@ wpt_object wpt_init(wave_object wave, int siglength, int J) {
return obj; return obj;
} }
cwt_object cwt_init(char* wave, double param,int siglength, double dt, int J) { cwt_object cwt_init(const char* wave, double param,int siglength, double dt, int J) {
cwt_object obj = NULL; cwt_object obj = NULL;
int N, i,nj2,ibase2,mother; int N, i,nj2,ibase2,mother;
double s0, dj; double s0, dj;
double t1; double t1;
int m, odd; int m, odd;
char *pdefault = "pow"; const char *pdefault = "pow";
m = (int)param; m = (int)param;
odd = 1; odd = 1;
@ -370,7 +377,7 @@ static void wconv(wt_object wt, double *sig, int N, double *filt, int L, double
} }
static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
int l,l2,isodd,i,t,len_avg; int l,l2,isodd,i,t,len_avg;
len_avg = wt->wave->lpd_len; len_avg = wt->wave->lpd_len;
@ -426,7 +433,7 @@ static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, do
} }
static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
int l, l2, isodd, i, t, len_avg; int l, l2, isodd, i, t, len_avg;
len_avg = wt->wave->lpd_len; len_avg = wt->wave->lpd_len;
@ -482,7 +489,7 @@ static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_c
} }
static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
int l, l2, isodd, i, t, len_avg; int l, l2, isodd, i, t, len_avg;
len_avg = wt->wave->lpd_len; len_avg = wt->wave->lpd_len;
@ -538,7 +545,7 @@ static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA,
} }
static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
int i,l, t, len_avg; int i,l, t, len_avg;
len_avg = wt->wave->lpd_len; len_avg = wt->wave->lpd_len;
@ -567,7 +574,7 @@ static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, do
} }
static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
int i, l, t, len_avg; int i, l, t, len_avg;
len_avg = wt->wave->lpd_len; len_avg = wt->wave->lpd_len;
@ -596,7 +603,7 @@ static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_c
} }
static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
int i, l, t, len_avg; int i, l, t, len_avg;
len_avg = wt->wave->lpd_len; len_avg = wt->wave->lpd_len;
@ -625,7 +632,7 @@ static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA,
} }
static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, int len_cA, double *cD, int len_cD) { static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, double *cD) {
int len_avg,D,lf; int len_avg,D,lf;
double *signal,*cA_undec; double *signal,*cA_undec;
len_avg = (wt->wave->lpd_len + wt->wave->hpd_len) / 2; len_avg = (wt->wave->lpd_len + wt->wave->hpd_len) / 2;
@ -700,7 +707,7 @@ static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, int len_cA, d
free(cA_undec); free(cA_undec);
} }
void dwt(wt_object wt,double *inp) { void dwt(wt_object wt,const double *inp) {
int i,J,temp_len,iter,N,lp; int i,J,temp_len,iter,N,lp;
int len_cA; int len_cA;
double *orig,*orig2; double *orig,*orig2;
@ -753,10 +760,10 @@ void dwt(wt_object wt,double *inp) {
len_cA = wt->length[J - iter]; len_cA = wt->length[J - iter];
N -= len_cA; N -= len_cA;
if ( !strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT") ) { if ( !strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT") ) {
dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); dwt1(wt, orig, temp_len, orig2, wt->params + N);
} }
else { else {
dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N);
} }
temp_len = wt->length[J - iter]; temp_len = wt->length[J - iter];
if (iter == J - 1) { if (iter == J - 1) {
@ -789,10 +796,10 @@ void dwt(wt_object wt,double *inp) {
len_cA = wt->length[J - iter]; len_cA = wt->length[J - iter];
N -= len_cA; N -= len_cA;
if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) { if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) {
dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); dwt1(wt, orig, temp_len, orig2, wt->params + N);
} }
else { else {
dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N);
} }
temp_len = wt->length[J - iter]; temp_len = wt->length[J - iter];
@ -817,7 +824,7 @@ void dwt(wt_object wt,double *inp) {
free(orig2); free(orig2);
} }
void wtree(wtree_object wt,double *inp) { void wtree(wtree_object wt,const double *inp) {
int i,J,temp_len,iter,N,lp,p2,k,N2,Np; int i,J,temp_len,iter,N,lp,p2,k,N2,Np;
int len_cA,t,t2,it1; int len_cA,t,t2,it1;
double *orig; double *orig;
@ -871,9 +878,9 @@ void wtree(wtree_object wt,double *inp) {
N = N2; N = N2;
for(k = 0; k < p2;++k) { for(k = 0; k < p2;++k) {
if (iter == 0) { if (iter == 0) {
wtree_per(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); wtree_per(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
} else { } else {
wtree_per(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); wtree_per(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
} }
N += 2 * len_cA; N += 2 * len_cA;
} }
@ -906,9 +913,9 @@ void wtree(wtree_object wt,double *inp) {
N = N2; N = N2;
for(k = 0; k < p2;++k) { for(k = 0; k < p2;++k) {
if (iter == 0) { if (iter == 0) {
wtree_sym(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); wtree_sym(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
} else { } else {
wtree_sym(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA, len_cA); wtree_sym(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
} }
N += 2 * len_cA; N += 2 * len_cA;
} }
@ -959,7 +966,7 @@ static int ipow2(int n) {
return p; return p;
} }
void dwpt(wpt_object wt, double *inp) { void dwpt(wpt_object wt, const double *inp) {
int i, J, temp_len, iter, N, lp, p2, k, N2, Np; int i, J, temp_len, iter, N, lp, p2, k, N2, Np;
int temp, elength, temp2,size,nodes,llb,n1,j; int temp, elength, temp2,size,nodes,llb,n1,j;
double eparam,v1,v2; double eparam,v1,v2;
@ -1026,10 +1033,10 @@ void dwpt(wpt_object wt, double *inp) {
N = N2; N = N2;
for (k = 0; k < p2; ++k) { for (k = 0; k < p2; ++k) {
if (iter == 0) { if (iter == 0) {
dwpt_per(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); dwpt_per(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA);
} }
else { else {
dwpt_per(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); dwpt_per(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA);
} }
wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam); wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
it2++; it2++;
@ -1066,10 +1073,10 @@ void dwpt(wpt_object wt, double *inp) {
N = N2; N = N2;
for (k = 0; k < p2; ++k) { for (k = 0; k < p2; ++k) {
if (iter == 0) { if (iter == 0) {
dwpt_sym(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); dwpt_sym(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA);
} }
else { else {
dwpt_sym(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA, len_cA); dwpt_sym(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA);
} }
wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam); wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
it2++; it2++;
@ -1327,7 +1334,7 @@ int getCWTScaleLength(int N) {
return J; return J;
} }
void setCWTScales(cwt_object wt, double s0, double dj,char *type,int power) { void setCWTScales(cwt_object wt, double s0, double dj,const char *type,int power) {
int i; int i;
strcpy(wt->type,type); strcpy(wt->type,type);
//s0*pow(2.0, (double)(j - 1)*dj); //s0*pow(2.0, (double)(j - 1)*dj);
@ -1353,7 +1360,7 @@ void setCWTScales(cwt_object wt, double s0, double dj,char *type,int power) {
wt->dj = dj; wt->dj = dj;
} }
void setCWTScaleVector(cwt_object wt, double *scale, int J,double s0,double dj) { void setCWTScaleVector(cwt_object wt, const double *scale, int J,double s0,double dj) {
int i; int i;
if (J != wt->J) { if (J != wt->J) {
@ -1378,7 +1385,7 @@ void setCWTPadding(cwt_object wt, int pad) {
} }
} }
void cwt(cwt_object wt, double *inp) { void cwt(cwt_object wt, const double *inp) {
int i, N, npad,nj2,j,j2; int i, N, npad,nj2,j,j2;
N = wt->siglength; N = wt->siglength;
if (wt->sflag == 0) { if (wt->sflag == 0) {
@ -1406,7 +1413,7 @@ void cwt(cwt_object wt, double *inp) {
} }
wt->smean /= N; wt->smean /= N;
cwavelet(inp, N, wt->dt, wt->mother, wt->m, wt->s0, wt->dj, wt->J,npad,wt->params, wt->params+nj2, wt->params+nj2+j, wt->params+nj2+j2); cwavelet(inp, N, wt->dt, wt->mother, wt->m, wt->s0,wt->dj,wt->J,npad,wt->params, wt->params+nj2, wt->params+nj2+j, wt->params+nj2+j2);
} }
@ -1477,7 +1484,7 @@ static void idwt1(wt_object wt,double *temp, double *cA_up,double *cA, int len_c
} }
static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, double *X) {
int len_avg,i,l,m,n,t,l2; int len_avg,i,l,m,n,t,l2;
len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
@ -1508,7 +1515,7 @@ static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_c
} }
} }
static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, double *X) {
int len_avg, i, l, m, n, t, v; int len_avg, i, l, m, n, t, v;
len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
@ -1592,7 +1599,7 @@ void idwt(wt_object wt, double *dwtop) {
//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp); idwt_per(wt,out, det_len, wt->output + iter, X_lp);
for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) { for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
out[k - lf/2 + 1] = X_lp[k]; out[k - lf/2 + 1] = X_lp[k];
} }
@ -1621,7 +1628,7 @@ void idwt(wt_object wt, double *dwtop) {
//idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out); //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
idwt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp); idwt_sym(wt, out, det_len, wt->output + iter, X_lp);
for (k = lf-2; k < 2 * det_len; ++k) { for (k = lf-2; k < 2 * det_len; ++k) {
out[k - lf + 2] = X_lp[k]; out[k - lf + 2] = X_lp[k];
} }
@ -1695,7 +1702,7 @@ void idwt(wt_object wt, double *dwtop) {
} }
static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, double *X) {
int len_avg, i, l, m, n, t, l2; int len_avg, i, l, m, n, t, l2;
len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
@ -1726,7 +1733,7 @@ static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, int len
} }
} }
static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, double *X) {
int len_avg, i, l, m, n, t, v; int len_avg, i, l, m, n, t, v;
len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
@ -1750,13 +1757,12 @@ static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, int len
} }
void idwpt(wpt_object wt, double *dwtop) { void idwpt(wpt_object wt, double *dwtop) {
int J, U, i, lf, N, k,p,l; int J, i, lf, k,p,l;
int app_len, det_len, index, n1, llb, index2, index3, index4,indexp,xlen; int app_len, det_len, index, n1, llb, index2, index3, index4,indexp,xlen;
double *X_lp, *X, *out, *out2; double *X_lp, *X, *out, *out2;
int *prep,*ptemp; int *prep,*ptemp;
J = wt->J; J = wt->J;
U = 2;
app_len = wt->length[0]; app_len = wt->length[0];
p = ipow2(J); p = ipow2(J);
lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
@ -1824,7 +1830,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = wt->output[index + k]; out[k] = wt->output[index + k];
out2[k] = wt->output[index + det_len + k]; out2[k] = wt->output[index + det_len + k];
} }
idwpt_per(wt, out, det_len, out2, det_len, X_lp); idwpt_per(wt, out, det_len, out2, X_lp);
for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
X[index3 + k - lf / 2 + 1] = X_lp[k]; X[index3 + k - lf / 2 + 1] = X_lp[k];
} }
@ -1839,7 +1845,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = wt->output[index + k]; out[k] = wt->output[index + k];
out2[k] = X[index4 + k]; out2[k] = X[index4 + k];
} }
idwpt_per(wt, out, det_len, out2, det_len, X_lp); idwpt_per(wt, out, det_len, out2, X_lp);
for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
X[index3 + k - lf / 2 + 1] = X_lp[k]; X[index3 + k - lf / 2 + 1] = X_lp[k];
} }
@ -1853,7 +1859,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = X[index4 + k]; out[k] = X[index4 + k];
out2[k] = wt->output[index + k]; out2[k] = wt->output[index + k];
} }
idwpt_per(wt, out, det_len, out2, det_len, X_lp); idwpt_per(wt, out, det_len, out2, X_lp);
for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
X[index3 + k - lf / 2 + 1] = X_lp[k]; X[index3 + k - lf / 2 + 1] = X_lp[k];
} }
@ -1867,7 +1873,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = X[index4 + k]; out[k] = X[index4 + k];
out2[k] = X[index4 + indexp + k]; out2[k] = X[index4 + indexp + k];
} }
idwpt_per(wt, out, det_len, out2, det_len, X_lp); idwpt_per(wt, out, det_len, out2, X_lp);
for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) { for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
X[index3 + k - lf / 2 + 1] = X_lp[k]; X[index3 + k - lf / 2 + 1] = X_lp[k];
} }
@ -1902,7 +1908,6 @@ void idwpt(wpt_object wt, double *dwtop) {
else if (!strcmp(wt->ext, "sym")) { else if (!strcmp(wt->ext, "sym")) {
app_len = wt->length[0]; app_len = wt->length[0];
det_len = wt->length[1]; det_len = wt->length[1];
N = 2 * wt->length[J] - 1;
//X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); //X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
index = 0; index = 0;
@ -1932,7 +1937,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = wt->output[index + k]; out[k] = wt->output[index + k];
out2[k] = wt->output[index + det_len + k]; out2[k] = wt->output[index + det_len + k];
} }
idwpt_sym(wt, out, det_len, out2, det_len, X_lp); idwpt_sym(wt, out, det_len, out2, X_lp);
for (k = lf - 2; k < 2 * det_len; ++k) { for (k = lf - 2; k < 2 * det_len; ++k) {
X[index3 + k - lf + 2] = X_lp[k]; X[index3 + k - lf + 2] = X_lp[k];
} }
@ -1947,7 +1952,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = wt->output[index + k]; out[k] = wt->output[index + k];
out2[k] = X[index4 + k]; out2[k] = X[index4 + k];
} }
idwpt_sym(wt, out, det_len, out2, det_len, X_lp); idwpt_sym(wt, out, det_len, out2, X_lp);
for (k = lf - 2; k < 2 * det_len; ++k) { for (k = lf - 2; k < 2 * det_len; ++k) {
X[index3 + k - lf + 2] = X_lp[k]; X[index3 + k - lf + 2] = X_lp[k];
} }
@ -1961,7 +1966,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = X[index4 + k]; out[k] = X[index4 + k];
out2[k] = wt->output[index + k]; out2[k] = wt->output[index + k];
} }
idwpt_sym(wt, out, det_len, out2, det_len, X_lp); idwpt_sym(wt, out, det_len, out2, X_lp);
for (k = lf - 2; k < 2 * det_len; ++k) { for (k = lf - 2; k < 2 * det_len; ++k) {
X[index3 + k - lf + 2] = X_lp[k]; X[index3 + k - lf + 2] = X_lp[k];
} }
@ -1975,7 +1980,7 @@ void idwpt(wpt_object wt, double *dwtop) {
out[k] = X[index4 + k]; out[k] = X[index4 + k];
out2[k] = X[index4 + indexp + k]; out2[k] = X[index4 + indexp + k];
} }
idwpt_sym(wt, out, det_len, out2, det_len, X_lp); idwpt_sym(wt, out, det_len, out2, X_lp);
for (k = lf - 2; k < 2 * det_len; ++k) { for (k = lf - 2; k < 2 * det_len; ++k) {
X[index3 + k - lf + 2] = X_lp[k]; X[index3 + k - lf + 2] = X_lp[k];
} }
@ -2029,7 +2034,7 @@ void idwpt(wpt_object wt, double *dwtop) {
} }
static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD) {
int l, l2, isodd, i, t, len_avg,j; int l, l2, isodd, i, t, len_avg,j;
len_avg = M * wt->wave->lpd_len; len_avg = M * wt->wave->lpd_len;
@ -2079,7 +2084,7 @@ static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_
} }
static void swt_fft(wt_object wt, double *inp) { static void swt_fft(wt_object wt, const double *inp) {
int i, J, temp_len, iter, M, N, len_filt; int i, J, temp_len, iter, M, N, len_filt;
int lenacc; int lenacc;
double *low_pass, *high_pass,*sig,*cA,*cD; double *low_pass, *high_pass,*sig,*cA,*cD;
@ -2164,16 +2169,15 @@ static void swt_fft(wt_object wt, double *inp) {
free(cD); free(cD);
} }
static void swt_direct(wt_object wt, double *inp) { static void swt_direct(wt_object wt, const double *inp) {
int i, J, temp_len, iter, M, N; int i, J, temp_len, iter, M;
int lenacc, len_filt; int lenacc;
double *cA, *cD; double *cA, *cD;
temp_len = wt->siglength; temp_len = wt->siglength;
J = wt->J; J = wt->J;
wt->length[0] = wt->length[J] = temp_len; wt->length[0] = wt->length[J] = temp_len;
wt->outlength = wt->length[J + 1] = (J + 1) * temp_len; wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
len_filt = wt->wave->filtlength;
M = 1; M = 1;
for (iter = 1; iter < J; ++iter) { for (iter = 1; iter < J; ++iter) {
M = 2 * M; M = 2 * M;
@ -2196,13 +2200,9 @@ static void swt_direct(wt_object wt, double *inp) {
lenacc -= temp_len; lenacc -= temp_len;
if (iter > 0) { if (iter > 0) {
M = 2 * M; M = 2 * M;
N = M * len_filt;
}
else {
N = len_filt;
} }
swt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len); swt_per(wt, M, wt->params, temp_len, cA, temp_len, cD);
for (i = 0; i < temp_len; ++i) { for (i = 0; i < temp_len; ++i) {
@ -2218,7 +2218,7 @@ static void swt_direct(wt_object wt, double *inp) {
} }
void swt(wt_object wt, double *inp) { void swt(wt_object wt, const double *inp) {
if (!strcmp(wt->method, "swt") && !strcmp(wt->cmethod, "direct") ) { if (!strcmp(wt->method, "swt") && !strcmp(wt->cmethod, "direct") ) {
swt_direct(wt,inp); swt_direct(wt,inp);
} }
@ -2376,7 +2376,7 @@ void iswt(wt_object wt, double *swtop) {
free(det2); free(det2);
} }
static void modwt_per(wt_object wt, int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { static void modwt_per(wt_object wt, int M, double *inp, double *cA, int len_cA, double *cD) {
int l, i, t, len_avg; int l, i, t, len_avg;
double s; double s;
double *filt; double *filt;
@ -2410,16 +2410,15 @@ static void modwt_per(wt_object wt, int M, double *inp, int N, double *cA, int l
free(filt); free(filt);
} }
void modwt(wt_object wt, double *inp) { void modwt(wt_object wt, const double *inp) {
int i, J, temp_len, iter, M, N; int i, J, temp_len, iter, M;
int lenacc, len_filt; int lenacc;
double *cA, *cD; double *cA, *cD;
temp_len = wt->siglength; temp_len = wt->siglength;
J = wt->J; J = wt->J;
wt->length[0] = wt->length[J] = temp_len; wt->length[0] = wt->length[J] = temp_len;
wt->outlength = wt->length[J + 1] = (J + 1) * temp_len; wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
len_filt = wt->wave->filtlength;
M = 1; M = 1;
for (iter = 1; iter < J; ++iter) { for (iter = 1; iter < J; ++iter) {
M = 2 * M; M = 2 * M;
@ -2442,13 +2441,9 @@ void modwt(wt_object wt, double *inp) {
lenacc -= temp_len; lenacc -= temp_len;
if (iter > 0) { if (iter > 0) {
M = 2 * M; M = 2 * M;
N = M * len_filt;
}
else {
N = len_filt;
} }
modwt_per(wt, M, wt->params, temp_len, cA, temp_len, cD, temp_len); modwt_per(wt, M, wt->params, cA, temp_len, cD);
for (i = 0; i < temp_len; ++i) { for (i = 0; i < temp_len; ++i) {
@ -2463,7 +2458,7 @@ void modwt(wt_object wt, double *inp) {
} }
static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, int len_cD, double *X) { static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, double *X) {
int len_avg, i, l, t; int len_avg, i, l, t;
double s; double s;
double *filt; double *filt;
@ -2497,14 +2492,12 @@ static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, i
} }
void imodwt(wt_object wt, double *dwtop) { void imodwt(wt_object wt, double *dwtop) {
int N, lf, iter, i, J, j, U; int N, iter, i, J, j;
int lenacc,M; int lenacc,M;
double *X; double *X;
N = wt->siglength; N = wt->siglength;
J = wt->J; J = wt->J;
U = 2;
lf = wt->wave->lpr_len;
lenacc = N; lenacc = N;
M = (int)pow(2.0, (double)J - 1.0); M = (int)pow(2.0, (double)J - 1.0);
//M = 1; //M = 1;
@ -2518,7 +2511,7 @@ void imodwt(wt_object wt, double *dwtop) {
if (iter > 0) { if (iter > 0) {
M = M / 2; M = M / 2;
} }
imodwt_per(wt, M, dwtop, N, wt->params + lenacc, N, X); imodwt_per(wt, M, dwtop, N, wt->params + lenacc, X);
/* /*
for (j = lf - 1; j < N; ++j) { for (j = lf - 1; j < N; ++j) {
dwtop[j - lf + 1] = X[j]; dwtop[j - lf + 1] = X[j];
@ -2536,7 +2529,7 @@ void imodwt(wt_object wt, double *dwtop) {
free(X); free(X);
} }
void setDWTExtension(wt_object wt, char *extension) { void setDWTExtension(wt_object wt, const char *extension) {
if (!strcmp(extension, "sym")) { if (!strcmp(extension, "sym")) {
strcpy(wt->ext, "sym"); strcpy(wt->ext, "sym");
} }
@ -2549,7 +2542,7 @@ void setDWTExtension(wt_object wt, char *extension) {
} }
} }
void setWTREEExtension(wtree_object wt, char *extension) { void setWTREEExtension(wtree_object wt, const char *extension) {
if (!strcmp(extension, "sym")) { if (!strcmp(extension, "sym")) {
strcpy(wt->ext, "sym"); strcpy(wt->ext, "sym");
} }
@ -2562,7 +2555,7 @@ void setWTREEExtension(wtree_object wt, char *extension) {
} }
} }
void setDWPTExtension(wpt_object wt, char *extension) { void setDWPTExtension(wpt_object wt, const char *extension) {
if (!strcmp(extension, "sym")) { if (!strcmp(extension, "sym")) {
strcpy(wt->ext, "sym"); strcpy(wt->ext, "sym");
} }
@ -2575,7 +2568,7 @@ void setDWPTExtension(wpt_object wt, char *extension) {
} }
} }
void setDWPTEntropy(wpt_object wt, char *entropy, double eparam) { void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam) {
if (!strcmp(entropy, "shannon")) { if (!strcmp(entropy, "shannon")) {
strcpy(wt->entropy, "shannon"); strcpy(wt->entropy, "shannon");
} }
@ -2596,7 +2589,7 @@ void setDWPTEntropy(wpt_object wt, char *entropy, double eparam) {
} }
} }
void setWTConv(wt_object wt, char *cmethod) { void setWTConv(wt_object wt, const char *cmethod) {
if (!strcmp(cmethod, "fft") || !strcmp(cmethod, "FFT")) { if (!strcmp(cmethod, "fft") || !strcmp(cmethod, "FFT")) {
strcpy(wt->cmethod, "fft"); strcpy(wt->cmethod, "fft");
} }

View File

@ -1,230 +0,0 @@
/*
Copyright (c) 2014, Rafat Hussain
*/
#ifndef WAVELIB_H_
#define WAVELIB_H_
#include "wtmath.h"
#include "cwt.h"
#ifdef __cplusplus
extern "C" {
#endif
#if defined(_MSC_VER)
#pragma warning(disable : 4200)
#pragma warning(disable : 4996)
#endif
#ifndef cplx_type
#define cplx_type double
#endif
typedef struct cplx_t {
cplx_type re;
cplx_type im;
} cplx_data;
typedef struct wave_set* wave_object;
wave_object wave_init(char* wname);
struct wave_set{
char wname[50];
int filtlength;// When all filters are of the same length. [Matlab uses zero-padding to make all filters of the same length]
int lpd_len;// Default filtlength = lpd_len = lpr_len = hpd_len = hpr_len
int hpd_len;
int lpr_len;
int hpr_len;
double *lpd;
double *hpd;
double *lpr;
double *hpr;
double params[0];
};
typedef struct wt_set* wt_object;
wt_object wt_init(wave_object wave,char* method, int siglength, int J);
struct wt_set{
wave_object wave;
conv_object cobj;
char method[10];
int siglength;// Length of the original signal.
int outlength;// Length of the output DWT vector
int lenlength;// Length of the Output Dimension Vector "length"
int J; // Number of decomposition Levels
int MaxIter;// Maximum Iterations J <= MaxIter
int even;// even = 1 if signal is of even length. even = 0 otherwise
char ext[10];// Type of Extension used - "per" or "sym"
char cmethod[10]; // Convolution Method - "direct" or "FFT"
int N; //
int cfftset;
int zpad;
int length[102];
double *output;
double params[0];
};
typedef struct wtree_set* wtree_object;
wtree_object wtree_init(wave_object wave, int siglength, int J);
struct wtree_set{
wave_object wave;
conv_object cobj;
char method[10];
int siglength;// Length of the original signal.
int outlength;// Length of the output DWT vector
int lenlength;// Length of the Output Dimension Vector "length"
int J; // Number of decomposition Levels
int MaxIter;// Maximum Iterations J <= MaxIter
int even;// even = 1 if signal is of even length. even = 0 otherwise
char ext[10];// Type of Extension used - "per" or "sym"
int N; //
int nodes;
int cfftset;
int zpad;
int length[102];
double *output;
int *nodelength;
int *coeflength;
double params[0];
};
typedef struct wpt_set* wpt_object;
wpt_object wpt_init(wave_object wave, int siglength, int J);
struct wpt_set{
wave_object wave;
conv_object cobj;
int siglength;// Length of the original signal.
int outlength;// Length of the output DWT vector
int lenlength;// Length of the Output Dimension Vector "length"
int J; // Number of decomposition Levels
int MaxIter;// Maximum Iterations J <= MaxIter
int even;// even = 1 if signal is of even length. even = 0 otherwise
char ext[10];// Type of Extension used - "per" or "sym"
char entropy[20];
double eparam;
int N; //
int nodes;
int length[102];
double *output;
double *costvalues;
double *basisvector;
int *nodeindex;
int *numnodeslevel;
int *coeflength;
double params[0];
};
typedef struct cwt_set* cwt_object;
cwt_object cwt_init(char* wave, double param, int siglength,double dt, int J);
struct cwt_set{
char wave[10];// Wavelet - morl/morlet,paul,dog/dgauss
int siglength;// Length of Input Data
int J;// Total Number of Scales
double s0;// Smallest scale. It depends on the sampling rate. s0 <= 2 * dt for most wavelets
double dt;// Sampling Rate
double dj;// Separation between scales. eg., scale = s0 * 2 ^ ( [0:N-1] *dj ) or scale = s0 *[0:N-1] * dj
char type[10];// Scale Type - Power or Linear
int pow;// Base of Power in case type = pow. Typical value is pow = 2
int sflag;
int pflag;
int npad;
int mother;
double m;// Wavelet parameter param
double smean;// Input Signal mean
cplx_data *output;
double *scale;
double *period;
double *coi;
double params[0];
};
void dwt(wt_object wt, double *inp);
void idwt(wt_object wt, double *dwtop);
void wtree(wtree_object wt, double *inp);
void dwpt(wpt_object wt, double *inp);
void idwpt(wpt_object wt, double *dwtop);
void swt(wt_object wt, double *inp);
void iswt(wt_object wt, double *swtop);
void modwt(wt_object wt, double *inp);
void imodwt(wt_object wt, double *dwtop);
void setDWTExtension(wt_object wt, char *extension);
void setWTREEExtension(wtree_object wt, char *extension);
void setDWPTExtension(wpt_object wt, char *extension);
void setDWPTEntropy(wpt_object wt, char *entropy, double eparam);
void setWTConv(wt_object wt, char *cmethod);
int getWTREENodelength(wtree_object wt, int X);
void getWTREECoeffs(wtree_object wt, int X, int Y, double *coeffs, int N);
int getDWPTNodelength(wpt_object wt, int X);
void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N);
void setCWTScales(cwt_object wt, double s0, double dj, char *type, int power);
void setCWTScaleVector(cwt_object wt, double *scale, int J, double s0, double dj);
void setCWTPadding(cwt_object wt, int pad);
void cwt(cwt_object wt, double *inp);
void icwt(cwt_object wt, double *cwtop);
int getCWTScaleLength(int N);
void wave_summary(wave_object obj);
void wt_summary(wt_object wt);
void wtree_summary(wtree_object wt);
void wpt_summary(wpt_object wt);
void cwt_summary(cwt_object wt);
void wave_free(wave_object object);
void wt_free(wt_object object);
void wtree_free(wtree_object object);
void wpt_free(wpt_object object);
void cwt_free(cwt_object object);
#ifdef __cplusplus
}
#endif
#endif /* WAVELIB_H_ */

View File

@ -1,30 +1,40 @@
add_executable(cwttest cwttest.c) add_executable(cwttest cwttest.c)
target_link_libraries(cwttest wavelib m) target_link_libraries(cwttest wavelib)
add_executable(dwttest dwttest.c) add_executable(dwttest dwttest.c)
target_link_libraries(dwttest wavelib m) target_link_libraries(dwttest wavelib)
add_executable(swttest swttest.c) add_executable(swttest swttest.c)
target_link_libraries(swttest wavelib m) target_link_libraries(swttest wavelib)
add_executable(modwttest modwttest.c) add_executable(modwttest modwttest.c)
target_link_libraries(modwttest wavelib m) target_link_libraries(modwttest wavelib)
add_executable(dwpttest dwpttest.c) add_executable(dwpttest dwpttest.c)
target_link_libraries(dwpttest wavelib m) target_link_libraries(dwpttest wavelib)
add_executable(wtreetest wtreetest.c) add_executable(wtreetest wtreetest.c)
target_link_libraries(wtreetest wavelib m) target_link_libraries(wtreetest wavelib)
add_executable(denoisetest denoisetest.c) add_executable(denoisetest denoisetest.c)
target_link_libraries(denoisetest wauxlib wavelib m) target_link_libraries(denoisetest wauxlib wavelib)
if(UNIX)
target_link_libraries(cwttest m)
target_link_libraries(dwttest m)
target_link_libraries(swttest m)
target_link_libraries(modwttest m)
target_link_libraries(dwpttest m)
target_link_libraries(wtreetest m)
target_link_libraries(denoisetest m)
endif()
set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest set_target_properties(cwttest dwttest swttest modwttest dwpttest wtreetest denoisetest
PROPERTIES PROPERTIES

View File

@ -4,10 +4,12 @@
#include <sstream> #include <sstream>
#include <iostream> #include <iostream>
#include <cstdlib>
#include <string>
#include <cmath> #include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "wavelib.h" #include "wavelib.h"
#include<vector> #include<vector>
@ -117,7 +119,6 @@ void ReconstructionTest()
double *inp,*out; double *inp,*out;
int N, i,J; int N, i,J;
double epsilon = 1e-15; double epsilon = 1e-15;
char *type = (char*) "dwt";
N = 79926; N = 79926;
@ -346,10 +347,10 @@ void DWPTReconstructionTest()
} }
void CWTReconstructionTest() { void CWTReconstructionTest() {
int i, j,N, J,subscale,a0,iter; int i, N, J,subscale,a0;
double *inp,*oup; double *inp,*oup;
double dt, dj,s0, pi,t; double dt, dj,s0, pi,t;
double val, epsilon; double epsilon;
int it1,it2; int it1,it2;
cwt_object wt; cwt_object wt;