Merge remote-tracking branch 'origin' into gh-pages

COPYRIGHT

@@ -0,0 +1,18 @@
+Copyright (c) 2014, Rafat Hussain
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+ 1.  Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+ 2.  Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer 
+ 	in the documentation and/or other materials provided with the distribution.
+
+ 3.  The name of the author may not be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
+BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+ OF SUCH DAMAGE.

README.md

@@ -0,0 +1,20 @@
+wavelib
+=======
+
+C Implementation of Wavelet Transform (DWT,SWT and MODWT)
+
+Methods Implemented
+
+DWT/IDWT A decimated Discrete Wavelet Transform implementation using implicit signal extension and up/downsampling so it is a fast implementation. A FFT based implementation is optional but will not be usually needed. Both periodic and symmetric options are available.
+
+SWT/ISWT Stationary Wavelet Transform. It works only for signal lengths that are multiples of 2^J where J is the number of decomposition levels. For signals of other lengths see MODWT implementation.
+
+MODWT/IMODWT Maximal Overlap Discrete Wavelet Transform is another undecimated transform. It is implemented for signals of any length but only orthogonal wavelets (Daubechies, Symlets and Coiflets) can be deployed. This implementation is based on the method laid out in "Wavelet Methods For Wavelet Analysis" by Donald Percival and Andrew Walden. 
+
+Documentation Available at - https://github.com/rafat/wavelib/wiki
+
+Live Demo (Emscripten) - http://rafat.github.io/wavelib/
+
+License - BSD 3-Clause
+
+Contace - rafat.hsn@gmail.com

header/wavelib.h

@@ -0,0 +1,129 @@
+#ifndef WAVELIB_H_
+#define WAVELIB_H_
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#ifndef fft_type
+#define fft_type double
+#endif
+
+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 fft_t {
+  fft_type re;
+  fft_type im;
+} fft_data;
+
+typedef struct fft_set* fft_object;
+
+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];
+};
+
+typedef struct fft_real_set* fft_real_object;
+
+fft_real_object fft_real_init(int N, int sgn);
+
+struct fft_real_set{
+	fft_object cobj;
+	fft_data twiddle2[1];
+};
+
+typedef struct conv_set* conv_object;
+
+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;
+};
+
+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];
+};
+
+void dwt11(wt_object wt, double *Vin, int M, double *Wout,
+	double *Vout);
+
+void dwt(wt_object wt, double *inp);
+
+void idwt(wt_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 setWTConv(wt_object wt, char *cmethod);
+
+void wave_summary(wave_object obj);
+
+void wt_summary(wt_object wt);
+
+void wave_free(wave_object object);
+
+void wt_free(wt_object object);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */
+
+

src/conv.c

@@ -0,0 +1,208 @@
+/*
+ * conv.c
+ *
+ *  Created on: May 1, 2013
+ *      Author: Rafat Hussain
+ */
+
+#include "conv.h"
+
+int factorf(int M) {
+	int N;
+	N = M;
+	while (N%7 == 0){
+			N = N/7;
+		}
+	while (N%3 == 0){
+			N = N/3;
+		}
+	while (N%5 == 0){
+			N = N/5;
+		}
+	while (N%2 == 0){
+			N = N/2;
+		}
+
+	return N;
+}
+
+
+int findnext(int M) {
+	int N;
+	N = M;
+
+	while (factorf(N) != 1) {
+		++N;
+	}
+
+	return N;
+
+}
+
+int findnexte(int M) {
+	int N;
+	N = M;
+
+	while (factorf(N) != 1 || N%2 != 0) {
+		++N;
+	}
+
+	return N;
+
+}
+
+
+conv_object conv_init(int N, int L) {
+	
+	conv_object obj = NULL;
+	int conv_len;
+	conv_len = N + L - 1;
+		
+	obj = (conv_object) malloc (sizeof(struct conv_set));
+	 
+	//obj->clen = npow2(conv_len);
+	//obj->clen = conv_len;
+	obj->clen = findnexte(conv_len);
+	obj->ilen1 = N;
+	obj->ilen2 = L;
+	
+	obj->fobj = fft_real_init(obj->clen,1);
+	obj->iobj = fft_real_init(obj->clen,-1);
+	
+	return obj;
+	
+}
+
+void conv_directx(fft_type *inp1,int N, fft_type *inp2, int L,fft_type *oup){
+	int M,k,n;
+	
+	M = N + L - 1;
+	
+	for (k = 0; k < M;++k) {
+		oup[k] = 0.0;
+		for ( n = 0; n < N; ++n) {
+			if ( (k-n) >= 0 && (k-n) < L ) {
+				oup[k]+= inp1[n] * inp2[k-n];
+			}
+		}
+		
+	}
+		
+}
+
+void conv_direct(fft_type *inp1,int N, fft_type *inp2, int L,fft_type *oup) {
+
+	int M,k,m,i;
+	fft_type t1,tmin;
+
+	M = N + L -1;
+	i = 0;
+
+	if (N >= L) {
+
+		for (k = 0; k < L; k++) {
+			oup[k] = 0.0;
+			for (m = 0; m <= k;m++) {
+				oup[k]+= inp1[m] * inp2[k-m];
+			}
+		}
+
+		for (k = L; k < M; k++) {
+			oup[k] = 0.0;
+			i++;
+			t1 = L + i;
+			tmin = MIN(t1,N);
+			for (m = i; m < tmin;m++) {
+				oup[k]+= inp1[m] * inp2[k-m];
+			}
+		}
+
+
+	} else {
+		for (k = 0; k < N; k++) {
+			oup[k] = 0.0;
+			for (m = 0; m <= k;m++) {
+				oup[k]+= inp2[m] * inp1[k-m];
+			}
+		}
+
+		for (k = N; k < M; k++) {
+			oup[k] = 0.0;
+			i++;
+			t1 = N + i;
+			tmin = MIN(t1,L);
+			for (m = i; m < tmin;m++) {
+				oup[k]+= inp2[m] * inp1[k-m];
+			}
+		}
+
+	}
+
+
+}
+
+
+void conv_fft(const conv_object obj,fft_type *inp1,fft_type *inp2,fft_type *oup) {
+	int i,N,L1,L2,ls;
+	fft_type* a;
+	fft_type* b;
+	fft_data* c;
+	fft_data* ao;
+	fft_data* bo;
+	fft_type* co;
+	
+	N = obj->clen;
+	L1 = obj->ilen1;
+	L2 = obj->ilen2;
+	ls = L1 + L2 - 1;
+	
+	a = (fft_type*) malloc (sizeof(fft_data) * N);
+	b = (fft_type*) malloc (sizeof(fft_data) * N);
+	c = (fft_data*) malloc (sizeof(fft_data) * N);
+	ao = (fft_data*) malloc (sizeof(fft_data) * N);
+	bo = (fft_data*) malloc (sizeof(fft_data) * N);
+	co = (fft_type*) malloc (sizeof(fft_data) * N);
+	
+	for (i = 0; i < N;i++) {
+		if (i < L1) {
+			a[i] = inp1[i];
+		} else {
+			a[i] = 0.0;
+		}
+		
+		if (i < L2) {
+			b[i] = inp2[i];
+		} else {
+			b[i] = 0.0;
+		}
+	
+	}
+	
+	fft_r2c_exec(obj->fobj,a,ao);
+	fft_r2c_exec(obj->fobj,b,bo);
+	
+	for (i = 0; i < N;i++) {
+		c[i].re = ao[i].re * bo[i].re - ao[i].im * bo[i].im;
+		c[i].im = ao[i].im * bo[i].re + ao[i].re * bo[i].im;
+	}
+	
+	fft_c2r_exec(obj->iobj,c,co);
+	
+	for (i = 0; i < ls;i++) {
+		oup[i] = co[i]/N;
+	}
+	
+	free(a);
+	free(b);
+	free(c);
+	free(ao);
+	free(bo);
+	free(co);
+	
+	
+}
+
+
+void free_conv(conv_object object) {
+	free(object);
+}

src/conv.h

@@ -0,0 +1,57 @@
+/*
+ * conv.h
+ *
+ *  Created on: May 1, 2013
+ *      Author: Rafat Hussain
+ */
+
+#ifndef CONV_H_
+#define CONV_H_
+
+#include "real.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#define MIN(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);
+
+struct conv_set{
+	fft_real_object fobj;
+	fft_real_object iobj;
+	int ilen1;
+	int ilen2;
+	int clen;
+};
+
+int factorf(int M);
+
+int findnext(int M);
+
+int findnexte(int M);
+
+void conv_direct(fft_type *inp1,int N, fft_type *inp2, int L,fft_type *oup);
+
+void conv_directx(fft_type *inp1,int N, fft_type *inp2, int L,fft_type *oup);
+
+//void conv_fft(const conv_object obj,fft_type *inp1,fft_type *inp2,fft_type *oup);
+
+//void conv_fft(const conv_object obj,fft_type *inp1,fft_type *inp2,fft_type *oup);
+
+void conv_fft(const conv_object obj,fft_type *inp1,fft_type *inp2,fft_type *oup);
+
+//void free_conv(conv_object object);
+
+void free_conv(conv_object object);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* CONV_H_ */

src/hsfft.h

@@ -0,0 +1,74 @@
+/*
+ * hsfft.h
+ *
+ *  Created on: Apr 14, 2013
+ *      Author: Rafat Hussain
+ */
+
+#ifndef HSFFT_H_
+#define HSFFT_H_
+
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#define PI2 6.28318530717958647692528676655900577
+
+#ifndef fft_type
+#define fft_type double
+#endif
+
+
+typedef struct fft_t {
+  fft_type re;
+  fft_type im;
+} fft_data;
+/*
+#define SADD(a,b) ((a)+(b))
+
+#define SSUB(a,b) ((a)+(b))
+
+#define SMUL(a,b) ((a)*(b))
+*/
+
+typedef struct fft_set* fft_object;
+
+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);
+
+int divideby(int M,int d);
+
+int dividebyN(int N);
+
+//void arrrev(int M, int* arr);
+
+int factors(int M, int* arr);
+
+void twiddle(fft_data *sig,int N, int radix);
+
+void longvectorN(fft_data *sig,int N, int *array, int M);
+
+void free_fft(fft_object object);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+
+
+#endif /* HSFFT_H_ */

src/real.c

@@ -0,0 +1,111 @@
+/*
+ * real.c
+ *
+ *  Created on: Apr 20, 2013
+ *      Author: Rafat Hussain
+ */
+#include <stdio.h>
+#include "real.h"
+
+fft_real_object fft_real_init(int N, int sgn) {
+	fft_real_object obj = NULL;
+	fft_type PI, theta;
+	int k;
+	
+	PI = 3.1415926535897932384626433832795;
+	
+	obj = (fft_real_object) malloc (sizeof(struct fft_real_set) + sizeof(fft_data)* (N/2));
+	
+	obj->cobj = fft_init(N/2,sgn);
+	
+	for (k = 0; k < N/2;++k) {
+		theta = PI2*k/N;
+		obj->twiddle2[k].re = cos(theta);
+		obj->twiddle2[k].im = sin(theta);
+		
+	}
+
+
+	return obj;
+	
+	
+}
+
+void fft_r2c_exec(fft_real_object obj,fft_type *inp,fft_data *oup) {
+	fft_data* cinp;
+	fft_data* coup;
+	int i,N2,N;
+	fft_type temp1,temp2;
+	N2 = obj->cobj->N;
+	N = N2*2;
+	
+	cinp = (fft_data*) malloc (sizeof(fft_data) * N2);
+	coup = (fft_data*) malloc (sizeof(fft_data) * N2);
+	
+	for (i = 0; i < N2; ++i) {
+		cinp[i].re = inp[2*i];
+		cinp[i].im = inp[2*i+1];
+	}
+	
+	fft_exec(obj->cobj,cinp,coup);
+	
+	oup[0].re = coup[0].re + coup[0].im;
+	oup[0].im = 0.0;
+	
+	for (i = 1; i < N2; ++i) {
+		temp1 = coup[i].im + coup[N2-i].im ;
+		temp2 = coup[N2-i].re - coup[i].re ;
+		oup[i].re = (coup[i].re + coup[N2-i].re + (temp1 * obj->twiddle2[i].re) + (temp2 * obj->twiddle2[i].im)) / 2.0;
+		oup[i].im = (coup[i].im - coup[N2-i].im + (temp2 * obj->twiddle2[i].re) - (temp1 * obj->twiddle2[i].im)) / 2.0;
+	}
+	
+	
+	
+	oup[N2].re = coup[0].re - coup[0].im;
+	oup[N2].im = 0.0;
+	
+	for (i = 1; i < N2;++i) {
+		oup[N-i].re = oup[i].re; 
+		oup[N-i].im = -oup[i].im; 
+	}
+	
+		
+	free(cinp);
+	free(coup);
+
+}
+
+void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup) {
+	
+	fft_data* cinp;
+	fft_data* coup;
+	int i,N2,N;
+	fft_type temp1,temp2;
+	N2 = obj->cobj->N;
+	N = N2*2;
+	
+	cinp = (fft_data*) malloc (sizeof(fft_data) * N2);
+	coup = (fft_data*) malloc (sizeof(fft_data) * N2);
+	
+	for (i = 0; i < N2; ++i) {
+		temp1 = -inp[i].im - inp[N2-i].im ;
+		temp2 = -inp[N2-i].re + inp[i].re ;
+		cinp[i].re = inp[i].re + inp[N2-i].re + (temp1 * obj->twiddle2[i].re) - (temp2 * obj->twiddle2[i].im);
+		cinp[i].im = inp[i].im - inp[N2-i].im + (temp2 * obj->twiddle2[i].re) + (temp1 * obj->twiddle2[i].im);
+	}
+	
+	fft_exec(obj->cobj,cinp,coup);
+	for (i = 0; i < N2; ++i) {
+		oup[2*i] = coup[i].re;
+		oup[2*i+1] = coup[i].im;
+	}
+	free(cinp);
+	free(coup);
+
+
+}
+
+void free_real_fft(fft_real_object object) {
+	free(object);
+}
+

src/real.h

@@ -0,0 +1,36 @@
+/*
+ * real.h
+ *
+ *  Created on: Apr 20, 2013
+ *      Author: Rafat Hussain
+ */
+
+#ifndef REAL_H_
+#define REAL_H_
+
+#include "hsfft.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct fft_real_set* fft_real_object;
+
+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_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup);
+
+void free_real_fft(fft_real_object object);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* REAL_H_ */

src/wavefilt.h

@@ -0,0 +1,22 @@
+#ifndef WAVEFILT_H_
+#define WAVEFILT_H_
+
+#include <stdio.h>
+#include "conv.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+int filtlength(char* name);
+
+int filtcoef(char* name, double *lp1, double *hp1, double *lp2, double *hp2);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVEFILT_H_ */

src/wavelib.h

@@ -0,0 +1,85 @@
+#ifndef WAVELIB_H_
+#define WAVELIB_H_
+
+#include "wtmath.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+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];
+};
+
+void dwt(wt_object wt, double *inp);
+
+void idwt(wt_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 setWTConv(wt_object wt, char *cmethod);
+
+void wave_summary(wave_object obj);
+
+void wt_summary(wt_object wt);
+
+void wave_free(wave_object object);
+
+void wt_free(wt_object object);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */
+
+

src/wtmath.c

@@ -0,0 +1,239 @@
+#include "wtmath.h"
+
+int upsamp(double *x, int lenx, int M, double *y) {
+	int N, i, j, k;
+
+	if (M < 0) {
+		return -1;
+	}
+
+	if (M == 0) {
+		for (i = 0; i < lenx; ++i) {
+			y[i] = x[i];
+		}
+		return lenx;
+	}
+
+	N = M * (lenx - 1) + 1;
+	j = 1;
+	k = 0;
+
+	for (i = 0; i < N; ++i) {
+		j--;
+		y[i] = 0.0;
+		if (j == 0) {
+			y[i] = x[k];
+			k++;
+			j = M;
+		}
+	}
+
+	return N;
+}
+
+int upsamp2(double *x, int lenx, int M, double *y) {
+	int N, i, j, k;
+	// upsamp2 returns even numbered output. Last value is set to zero
+	if (M < 0) {
+		return -1;
+	}
+
+	if (M == 0) {
+		for (i = 0; i < lenx; ++i) {
+			y[i] = x[i];
+		}
+		return lenx;
+	}
+
+	N = M * lenx;
+	j = 1;
+	k = 0;
+
+	for (i = 0; i < N; ++i) {
+		j--;
+		y[i] = 0.0;
+		if (j == 0) {
+			y[i] = x[k];
+			k++;
+			j = M;
+		}
+	}
+
+	return N;
+}
+
+int downsamp(double *x, int lenx, int M, double *y) {
+	int N, i;
+
+	if (M < 0) {
+		return -1;
+	}
+	if (M == 0) {
+		for (i = 0; i < lenx; ++i) {
+			y[i] = x[i];
+		}
+		return lenx;
+	}
+
+	N = (lenx - 1) / M + 1;
+
+	for (i = 0; i < N; ++i) {
+		y[i] = x[i*M];
+	}
+
+	return N;
+}
+/*
+int per_ext(double *sig, int len, int a,double *oup) {
+	int i,len2;
+	// oup is of length len + (len%2) + 2 * a
+	for (i = 0; i < len; ++i) {
+		oup[a + i] = sig[i];
+	}
+	len2 = len;
+	if ((len % 2) != 0) {
+		len2 = len + 1;
+		oup[a + len] = sig[len - 1];
+	}
+	for (i = 0; i < a; ++i) {
+		oup[a-1-i] = sig[len - 1 - i];
+		oup[len2 + a + i] = sig[i];
+	}
+
+	return len2;
+
+}
+*/
+
+int per_ext(double *sig, int len, int a, double *oup) {
+	int i, len2;
+	double temp1;
+	double temp2;
+	for (i = 0; i < len; ++i) {
+		oup[a + i] = sig[i];
+	}
+	len2 = len;
+	if ((len % 2) != 0) {
+		len2 = len + 1;
+		oup[a + len] = sig[len - 1];
+	}
+	for (i = 0; i < a; ++i) {
+		temp1 = oup[a + i];
+		temp2 = oup[a + len2 - 1 - i];
+		oup[a - 1 - i] = temp2;
+		oup[len2 + a + i] = temp1;
+	}
+	return len2;
+}
+/*
+int symm_ext(double *sig, int len, int a, double *oup) {
+	int i, len2;
+	// oup is of length len + 2 * a
+	for (i = 0; i < len; ++i) {
+		oup[a + i] = sig[i];
+	}
+	len2 = len;
+	for (i = 0; i < a; ++i) {
+		oup[a - 1 - i] = sig[i];
+		oup[len2 + a + i] = sig[len - 1 - i];
+	}
+
+	return len2;
+
+}
+*/
+
+int symm_ext(double *sig, int len, int a, double *oup) {
+	int i, len2;
+	double temp1;
+	double temp2;
+	// oup is of length len + 2 * a
+	for (i = 0; i < len; ++i) {
+		oup[a + i] = sig[i];
+	}
+	len2 = len;
+	for (i = 0; i < a; ++i) {
+		temp1 = oup[a + i];
+		temp2 = oup[a + len2 - 1 - i];
+		oup[a - 1 - i] = temp1;
+		oup[len2 + a + i] = temp2;
+	}
+
+	return len2;
+
+}
+
+static int isign(int N) {
+	int M;
+	if (N >= 0) {
+		M = 1;
+	}
+	else {
+		M = -1;
+	}
+
+	return M;
+}
+
+static int iabs(int N) {
+	if (N >= 0) {
+		return N;
+	}
+	else {
+		return -N;
+	}
+}
+
+void circshift(double *array, int N, int L) {
+	int i;
+	double *temp;
+	if (iabs(L) > N) {
+		L = isign(L) * (iabs(L) % N);
+	}
+	if (L < 0) {
+		L = (N + L) % N;
+	}
+
+	temp = (double*)malloc(sizeof(double) * L);
+
+	for (i = 0; i < L; ++i) {
+		temp[i] = array[i];
+	}
+
+	for (i = 0; i < N - L; ++i) {
+		array[i] = array[i + L];
+	}
+
+	for (i = 0; i < L; ++i) {
+		array[N - L + i] = temp[i];
+	}
+
+	free(temp);
+}
+
+int testSWTlength(int N, int J) {
+	int ret,div,i;
+	ret = 1;
+
+	div = 1;
+	for (i = 0; i < J; ++i) {
+		div *= 2;
+	}
+
+	if (N % div) {
+		ret = 0;
+	}
+
+	return ret;
+
+}
+
+int wmaxiter(int sig_len, int filt_len) {
+	int lev;
+	double temp;
+
+	temp = log((double)sig_len / ((double)filt_len - 1.0)) / log(2.0);
+	lev = (int)temp;
+
+	return lev;
+}

src/wtmath.h

@@ -0,0 +1,31 @@
+#ifndef WTMATH_H_
+#define WTMATH_H_
+
+#include "wavefilt.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int upsamp(double *x, int lenx, int M, double *y);
+
+int upsamp2(double *x, int lenx, int M, double *y);
+
+int downsamp(double *x, int lenx, int M, double *y);
+
+int per_ext(double *sig, int len, int a,double *oup);
+
+int symm_ext(double *sig, int len, int a,double *oup);
+
+void circshift(double *array, int N, int L);
+
+int testSWTlength(int N, int J);
+
+int wmaxiter(int sig_len, int filt_len);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif /* WAVELIB_H_ */

test/dwttest.c

@@ -0,0 +1,83 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+int main() {
+	wave_object obj;
+	wt_object wt;
+	double *inp,*out,*diff;
+	int N, i,J;
+
+	FILE *ifp;
+	double temp[1200];
+
+	char *name = "db4";
+	obj = wave_init(name);// Initialize the wavelet
+
+	ifp = fopen("signal.txt", "r");
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+	N = 256;
+
+	inp = (double*)malloc(sizeof(double)* N);
+	out = (double*)malloc(sizeof(double)* N);
+	diff = (double*)malloc(sizeof(double)* N);
+	//wmean = mean(temp, N);
+
+	for (i = 0; i < N; ++i) {
+		inp[i] = temp[i];
+		//printf("%g \n",inp[i]);
+	}
+	J = 3;
+
+	wt = wt_init(obj, "dwt", N, J);// Initialize the wavelet transform object
+	setDWTExtension(wt, "sym");// Options are "per" and "sym". Symmetric is the default option
+	setWTConv(wt, "direct");
+	
+	dwt(wt, inp);// Perform DWT
+	//DWT output can be accessed using wt->output vector. Use wt_summary to find out how to extract appx and detail coefficients
+	
+	for (i = 0; i < wt->outlength; ++i) {
+		printf("%g ",wt->output[i]);
+	}
+	
+	idwt(wt, out);// Perform IDWT (if needed)
+	// Test Reconstruction
+	for (i = 0; i < wt->siglength; ++i) {
+		diff[i] = out[i] - inp[i];
+	}
+	
+	printf("\n MAX %g \n", absmax(diff, wt->siglength)); // If Reconstruction succeeded then the output should be a small value.
+	
+	wt_summary(wt);// Prints the full summary.
+	wave_free(obj);
+	wt_free(wt);
+
+	free(inp);
+	free(out);
+	free(diff);
+	return 0;
+}

test/modwttest.c

@@ -0,0 +1,86 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+int main() {
+	wave_object obj;
+	wt_object wt;
+	double *inp, *out, *diff;
+	int N, i, J;
+
+	FILE *ifp;
+	double temp[1200];
+
+	char *name = "db4";
+	obj = wave_init(name);
+	wave_summary(obj);
+
+	ifp = fopen("signal.txt", "r");
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+	N = 177;
+
+	inp = (double*)malloc(sizeof(double)* N);
+	out = (double*)malloc(sizeof(double)* N);
+	diff = (double*)malloc(sizeof(double)* N);
+	//wmean = mean(temp, N);
+
+	for (i = 0; i < N; ++i) {
+		inp[i] = temp[i];
+		//printf("%g \n",inp[i]);
+	}
+	J = 2;
+
+	wt = wt_init(obj, "modwt", N, J);// Initialize the wavelet transform object
+	
+	modwt(wt, inp);// Perform MODWT
+	//MODWT output can be accessed using wt->output vector. Use wt_summary to find out how to extract appx and detail coefficients
+	
+	for (i = 0; i < wt->outlength; ++i) {
+		printf("%g ",wt->output[i]);
+	}
+
+	imodwt(wt, out);// Perform ISWT (if needed)
+	// Test Reconstruction
+
+
+	for (i = 0; i < wt->siglength; ++i) {
+		diff[i] = out[i] - inp[i];
+	}
+
+	printf("\n MAX %g \n", absmax(diff, wt->siglength));// If Reconstruction succeeded then the output should be a small value.
+	
+	wt_summary(wt);// Prints the full summary.
+
+	wave_free(obj);
+	wt_free(wt);
+
+	free(inp);
+	free(out);
+	free(diff);
+	return 0;
+}
+

test/signal.txt

@@ -0,0 +1 @@
+  -18.3237

test/swttest.c

@@ -0,0 +1,87 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "../header/wavelib.h"
+
+double absmax(double *array, int N) {
+	double max;
+	int i;
+
+	max = 0.0;
+	for (i = 0; i < N; ++i) {
+		if (fabs(array[i]) >= max) {
+			max = fabs(array[i]);
+		}
+	}
+
+	return max;
+}
+
+int main() {
+	wave_object obj;
+	wt_object wt;
+	double *inp, *out, *diff;
+	int N, i, J;
+
+	FILE *ifp;
+	double temp[1200];
+
+	char *name = "bior3.5";
+	obj = wave_init(name);// Initialize the wavelet
+
+	ifp = fopen("signal.txt", "r");
+	i = 0;
+	if (!ifp) {
+		printf("Cannot Open File");
+		exit(100);
+	}
+	while (!feof(ifp)) {
+		fscanf(ifp, "%lf \n", &temp[i]);
+		i++;
+	}
+	N = 256;
+
+	inp = (double*)malloc(sizeof(double)* N);
+	out = (double*)malloc(sizeof(double)* N);
+	diff = (double*)malloc(sizeof(double)* N);
+	//wmean = mean(temp, N);
+
+	for (i = 0; i < N; ++i) {
+		inp[i] = temp[i];
+		//printf("%g \n",inp[i]);
+	}
+	J = 1;
+
+	wt = wt_init(obj, "swt", N, J);// Initialize the wavelet transform object
+	setWTConv(wt, "direct");
+
+	
+	swt(wt, inp);// Perform SWT
+	//SWT output can be accessed using wt->output vector. Use wt_summary to find out how to extract appx and detail coefficients
+	
+	for (i = 0; i < wt->outlength; ++i) {
+		printf("%g ",wt->output[i]);
+	}
+
+	iswt(wt, out);// Perform ISWT (if needed)
+	// Test Reconstruction
+	
+
+	for (i = 0; i < wt->siglength; ++i) {
+		diff[i] = out[i] - inp[i];
+	}
+
+	printf("\n MAX %g \n", absmax(diff, wt->siglength));// If Reconstruction succeeded then the output should be a small value.
+	
+	wt_summary(wt);// Prints the full summary.
+
+
+	wave_free(obj);
+	wt_free(wt);
+
+	free(inp);
+	free(out);
+	free(diff);
+	return 0;
+}