Commit : Code Added

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_ */