From b344fc7ac130f15f6130c969d43f65660adec059 Mon Sep 17 00:00:00 2001 From: Rafat Hussain Date: Mon, 15 Dec 2014 15:47:46 +0530 Subject: [PATCH] Commit : Code Added --- COPYRIGHT | 18 + header/wavelib.h | 129 ++ src/conv.c | 208 +++ src/conv.h | 57 + src/hsfft.c | 1860 ++++++++++++++++++++++++++ src/hsfft.h | 74 ++ src/real.c | 111 ++ src/real.h | 36 + src/wavefilt.c | 3225 ++++++++++++++++++++++++++++++++++++++++++++++ src/wavefilt.h | 22 + src/wavelib.c | 1274 ++++++++++++++++++ src/wavelib.h | 85 ++ src/wtmath.c | 239 ++++ src/wtmath.h | 31 + test/dwttest.c | 83 ++ test/modwttest.c | 86 ++ test/signal.txt | 1 + test/swttest.c | 87 ++ 18 files changed, 7626 insertions(+) create mode 100644 COPYRIGHT create mode 100644 header/wavelib.h create mode 100644 src/conv.c create mode 100644 src/conv.h create mode 100644 src/hsfft.c create mode 100644 src/hsfft.h create mode 100644 src/real.c create mode 100644 src/real.h create mode 100644 src/wavefilt.c create mode 100644 src/wavefilt.h create mode 100644 src/wavelib.c create mode 100644 src/wavelib.h create mode 100644 src/wtmath.c create mode 100644 src/wtmath.h create mode 100644 test/dwttest.c create mode 100644 test/modwttest.c create mode 100644 test/signal.txt create mode 100644 test/swttest.c diff --git a/COPYRIGHT b/COPYRIGHT new file mode 100644 index 0000000..19c19fb --- /dev/null +++ b/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. \ No newline at end of file diff --git a/header/wavelib.h b/header/wavelib.h new file mode 100644 index 0000000..be1581a --- /dev/null +++ b/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_ */ + + diff --git a/src/conv.c b/src/conv.c new file mode 100644 index 0000000..768319f --- /dev/null +++ b/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); +} diff --git a/src/conv.h b/src/conv.h new file mode 100644 index 0000000..40b2548 --- /dev/null +++ b/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_ */ diff --git a/src/hsfft.c b/src/hsfft.c new file mode 100644 index 0000000..bf4c650 --- /dev/null +++ b/src/hsfft.c @@ -0,0 +1,1860 @@ +/* + * hsfft.c + * + * Created on: Apr 14, 2013 + * Author: Rafat Hussain + */ + + +#include "hsfft.h" + +fft_object fft_init(int N, int sgn) { + fft_object obj = NULL; + // Change N/2 to N-1 for longvector case + + int twi_len,ct,out; + out = dividebyN(N); + + if (out == 1) { + obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (N-1)); + obj->lf = factors(N,obj->factors); + longvectorN(obj->twiddle,N,obj->factors,obj->lf); + twi_len = N; + obj->lt = 0; + } else { + int K,M; + K = (int) pow(2.0,ceil(log10(N)/log10(2.0))); + + if (K < 2 * N - 2) { + M = K * 2; + } else { + M = K; + } + obj = (fft_object) malloc (sizeof(struct fft_set) + sizeof(fft_data)* (M-1)); + obj->lf = factors(M,obj->factors); + longvectorN(obj->twiddle,M,obj->factors,obj->lf); + obj->lt = 1; + twi_len = M; + } + + + obj->N = N; + obj->sgn = sgn; + + if (sgn == -1) { + for(ct = 0; ct < twi_len;ct++) { + (obj->twiddle+ct)->im = - (obj->twiddle+ct)->im; + + } + } + + return obj; +} + + + +static void mixed_radix_dit_rec(fft_data *op,fft_data *ip,const fft_object obj, int sgn, int N,int l,int inc) { + + int radix,m,ll; + if (N > 1) { + radix = obj->factors[inc]; + //printf("%d \n",radix); + } + + if (N == 1) { + + op[0].re = ip[0].re; + op[0].im = ip[0].im; + + } else if (N == 2) { + fft_type tau1r,tau1i; + op[0].re = ip[0].re; + op[0].im = ip[0].im; + + op[1].re = ip[l].re; + op[1].im = ip[l].im; + + tau1r = op[0].re; + tau1i = op[0].im; + + + op[0].re = tau1r + op[1].re; + op[0].im = tau1i + op[1].im; + + op[1].re = tau1r - op[1].re; + op[1].im = tau1i - op[1].im; + + } else if (N == 3) { + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i; + op[0].re = ip[0].re; + op[0].im = ip[0].im; + + op[1].re = ip[l].re; + op[1].im = ip[l].im; + + op[2].re = ip[2*l].re; + op[2].im = ip[2*l].im; + + tau0r = op[1].re + op[2].re; + tau0i = op[1].im + op[2].im; + + tau1r = sgn * 0.86602540378 * (op[1].re - op[2].re); + tau1i = sgn * 0.86602540378 * (op[1].im - op[2].im); + + tau2r = op[0].re - tau0r * 0.5000000000; + tau2i = op[0].im - tau0i * 0.5000000000; + + op[0].re = tau0r + op[0].re ; + op[0].im = tau0i + op[0].im; + + op[1].re = tau2r + tau1i; + op[1].im = tau2i - tau1r; + + op[2].re = tau2r - tau1i; + op[2].im = tau2i + tau1r; + + return; + + + } else if (N == 4) { + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i; + op[0].re = ip[0].re; + op[0].im = ip[0].im; + + op[1].re = ip[l].re; + op[1].im = ip[l].im; + + op[2].re = ip[2*l].re; + op[2].im = ip[2*l].im; + + op[3].re = ip[3*l].re; + op[3].im = ip[3*l].im; + + tau0r = op[0].re + op[2].re; + tau0i = op[0].im + op[2].im; + + tau1r = op[0].re - op[2].re; + tau1i = op[0].im - op[2].im; + + tau2r = op[1].re + op[3].re; + tau2i = op[1].im + op[3].im; + + tau3r = sgn* (op[1].re - op[3].re); + tau3i = sgn* (op[1].im - op[3].im); + + op[0].re = tau0r + tau2r ; + op[0].im = tau0i + tau2i; + + op[1].re = tau1r + tau3i; + op[1].im = tau1i - tau3r; + + op[2].re = tau0r - tau2r; + op[2].im = tau0i - tau2i; + + op[3].re = tau1r - tau3i; + op[3].im = tau1i + tau3r; + + + + } else if (N == 5) { + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i,tau4r,tau4i,tau5r,tau5i,tau6r,tau6i; + fft_type c1,c2,s1,s2; + op[0].re = ip[0].re; + op[0].im = ip[0].im; + + op[1].re = ip[l].re; + op[1].im = ip[l].im; + + op[2].re = ip[2*l].re; + op[2].im = ip[2*l].im; + + op[3].re = ip[3*l].re; + op[3].im = ip[3*l].im; + + op[4].re = ip[4*l].re; + op[4].im = ip[4*l].im; + + c1 = 0.30901699437; + c2 = -0.80901699437; + s1 = 0.95105651629; + s2 = 0.58778525229; + + tau0r = op[1].re + op[4].re; + tau2r = op[1].re - op[4].re; + tau0i = op[1].im + op[4].im; + tau2i = op[1].im - op[4].im; + + tau1r = op[2].re + op[3].re; + tau3r = op[2].re - op[3].re; + tau1i = op[2].im + op[3].im; + tau3i = op[2].im - op[3].im; + + + tau4r = c1 * tau0r + c2 * tau1r; + tau4i = c1 * tau0i + c2 * tau1i; + + //tau5r = sgn * ( s1 * tau2r + s2 * tau3r); + //tau5i = sgn * ( s1 * tau2i + s2 * tau3i); + + if (sgn == 1) { + tau5r = s1 * tau2r + s2 * tau3r; + tau5i = s1 * tau2i + s2 * tau3i; + + } else { + tau5r = -s1 * tau2r - s2 * tau3r; + tau5i = -s1 * tau2i - s2 * tau3i; + + } + + tau6r = op[0].re + tau4r; + tau6i = op[0].im + tau4i; + + op[1].re = tau6r + tau5i; + op[1].im = tau6i - tau5r; + + op[4].re = tau6r - tau5i; + op[4].im = tau6i + tau5r; + + tau4r = c2 * tau0r + c1 * tau1r; + tau4i = c2 * tau0i + c1 * tau1i; + + //tau5r = sgn * ( s2 * tau2r - s1 * tau3r); + //tau5i = sgn * ( s2 * tau2i - s1 * tau3i); + + if (sgn == 1) { + tau5r = s2 * tau2r - s1 * tau3r; + tau5i = s2 * tau2i - s1 * tau3i; + + } else { + tau5r = -s2 * tau2r + s1 * tau3r; + tau5i = -s2 * tau2i + s1 * tau3i; + + } + + tau6r = op[0].re + tau4r; + tau6i = op[0].im + tau4i; + + op[2].re = tau6r + tau5i; + op[2].im = tau6i - tau5r; + + op[3].re = tau6r - tau5i; + op[3].im = tau6i + tau5r; + + + + op[0].re += tau0r + tau1r; + op[0].im += tau0i + tau1i; + + + } else if (N == 7) { + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i,tau4r,tau4i,tau5r,tau5i,tau6r,tau6i,tau7r,tau7i; + fft_type c1,c2,c3,s1,s2,s3; + op[0].re = ip[0].re; + op[0].im = ip[0].im; + + op[1].re = ip[l].re; + op[1].im = ip[l].im; + + op[2].re = ip[2*l].re; + op[2].im = ip[2*l].im; + + op[3].re = ip[3*l].re; + op[3].im = ip[3*l].im; + + op[4].re = ip[4*l].re; + op[4].im = ip[4*l].im; + + op[5].re = ip[5*l].re; + op[5].im = ip[5*l].im; + + op[6].re = ip[6*l].re; + op[6].im = ip[6*l].im; + + c1 = 0.62348980185; + c2 = -0.22252093395; + c3 = -0.9009688679; + s1 = 0.78183148246; + s2 = 0.97492791218; + s3 = 0.43388373911; + + tau0r = op[1].re + op[6].re; + tau3r = op[1].re - op[6].re; + + tau0i = op[1].im + op[6].im; + tau3i = op[1].im - op[6].im; + + tau1r = op[2].re + op[5].re; + tau4r = op[2].re - op[5].re; + + tau1i = op[2].im + op[5].im; + tau4i = op[2].im - op[5].im; + + tau2r = op[3].re + op[4].re; + tau5r = op[3].re - op[4].re; + + tau2i = op[3].im + op[4].im; + tau5i = op[3].im - op[4].im; + + + tau6r = op[0].re + c1 * tau0r + c2 * tau1r + c3 * tau2r; + tau6i = op[0].im + c1 * tau0i + c2 * tau1i + c3 * tau2i; + + //tau7r = sgn * ( -s1 * tau3r - s2 * tau4r - s3 * tau5r); + //tau7i = sgn * ( -s1 * tau3i - s2 * tau4i - s3 * tau5i); + + if (sgn == 1) { + tau7r = -s1 * tau3r - s2 * tau4r - s3 * tau5r; + tau7i = -s1 * tau3i - s2 * tau4i - s3 * tau5i; + + } else { + tau7r = s1 * tau3r + s2 * tau4r + s3 * tau5r; + tau7i = s1 * tau3i + s2 * tau4i + s3 * tau5i; + + } + + + op[1].re = tau6r - tau7i; + op[6].re = tau6r + tau7i; + + op[1].im = tau6i + tau7r; + op[6].im = tau6i - tau7r; + + tau6r = op[0].re + c2 * tau0r + c3 * tau1r + c1 * tau2r; + tau6i = op[0].im + c2 * tau0i + c3 * tau1i + c1 * tau2i; + + //tau7r = sgn * ( -s2 * tau3r + s3 * tau4r + s1 * tau5r); + //tau7i = sgn * ( -s2 * tau3i + s3 * tau4i + s1 * tau5i); + + if (sgn == 1) { + tau7r = -s2 * tau3r + s3 * tau4r + s1 * tau5r; + tau7i = -s2 * tau3i + s3 * tau4i + s1 * tau5i; + } else { + tau7r = s2 * tau3r - s3 * tau4r - s1 * tau5r; + tau7i = s2 * tau3i - s3 * tau4i - s1 * tau5i; + + } + + + op[2].re = tau6r - tau7i; + op[5].re = tau6r + tau7i; + op[2].im = tau6i + tau7r; + op[5].im = tau6i - tau7r; + + tau6r = op[0].re + c3 * tau0r + c1 * tau1r + c2 * tau2r; + tau6i = op[0].im + c3 * tau0i + c1 * tau1i + c2 * tau2i; + + //tau7r = sgn * ( -s3 * tau3r + s1 * tau4r - s2 * tau5r); + //tau7i = sgn * ( -s3 * tau3i + s1 * tau4i - s2 * tau5i); + + if (sgn == 1) { + tau7r = -s3 * tau3r + s1 * tau4r - s2 * tau5r; + tau7i = -s3 * tau3i + s1 * tau4i - s2 * tau5i; + + } else { + tau7r = s3 * tau3r - s1 * tau4r + s2 * tau5r; + tau7i = s3 * tau3i - s1 * tau4i + s2 * tau5i; + + } + + + op[3].re = tau6r - tau7i; + op[4].re = tau6r + tau7i; + op[3].im = tau6i + tau7r; + op[4].im = tau6i - tau7r; + + op[0].re += tau0r + tau1r + tau2r; + op[0].im += tau0i + tau1i + tau2i; + + + } else if (N == 8) { + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i,tau4r,tau4i,tau5r,tau5i,tau6r,tau6i,tau7r,tau7i,tau8r,tau8i,tau9r,tau9i; + fft_type c1,s1,temp1r,temp1i,temp2r,temp2i; + op[0].re = ip[0].re; + op[0].im = ip[0].im; + + op[1].re = ip[l].re; + op[1].im = ip[l].im; + + op[2].re = ip[2*l].re; + op[2].im = ip[2*l].im; + + op[3].re = ip[3*l].re; + op[3].im = ip[3*l].im; + + op[4].re = ip[4*l].re; + op[4].im = ip[4*l].im; + + op[5].re = ip[5*l].re; + op[5].im = ip[5*l].im; + + op[6].re = ip[6*l].re; + op[6].im = ip[6*l].im; + + op[7].re = ip[7*l].re; + op[7].im = ip[7*l].im; + + c1 = 0.70710678118654752440084436210485; + s1 = 0.70710678118654752440084436210485; + + tau0r = op[0].re + op[4].re; + tau4r = op[0].re - op[4].re; + + tau0i = op[0].im + op[4].im; + tau4i = op[0].im - op[4].im; + + tau1r = op[1].re + op[7].re; + tau5r = op[1].re - op[7].re; + + tau1i = op[1].im + op[7].im; + tau5i = op[1].im - op[7].im; + + tau2r = op[3].re + op[5].re; + tau6r = op[3].re - op[5].re; + + tau2i = op[3].im + op[5].im; + tau6i = op[3].im - op[5].im; + + tau3r = op[2].re + op[6].re; + tau7r = op[2].re - op[6].re; + + tau3i = op[2].im + op[6].im; + tau7i = op[2].im - op[6].im; + + op[0].re = tau0r + tau1r + tau2r + tau3r; + op[0].im = tau0i + tau1i + tau2i + tau3i; + + op[4].re = tau0r - tau1r - tau2r + tau3r; + op[4].im = tau0i - tau1i - tau2i + tau3i; + + temp1r = tau1r - tau2r; + temp1i = tau1i - tau2i; + + temp2r = tau5r + tau6r; + temp2i = tau5i + tau6i; + + tau8r = tau4r + c1 * temp1r; + tau8i = tau4i + c1 * temp1i; + + //tau9r = sgn * ( -s1 * temp2r - tau7r); + //tau9i = sgn * ( -s1 * temp2i - tau7i); + + if (sgn == 1) { + tau9r = -s1 * temp2r - tau7r; + tau9i = -s1 * temp2i - tau7i; + + } else { + tau9r = s1 * temp2r + tau7r; + tau9i = s1 * temp2i + tau7i; + + } + + + op[1].re = tau8r - tau9i; + op[1].im = tau8i + tau9r; + + op[7].re = tau8r + tau9i; + op[7].im = tau8i - tau9r; + + tau8r = tau0r - tau3r; + tau8i = tau0i - tau3i; + + //tau9r = sgn * ( -tau5r + tau6r); + //tau9i = sgn * ( -tau5i + tau6i); + + if (sgn == 1) { + tau9r = -tau5r + tau6r; + tau9i = -tau5i + tau6i; + + } else { + tau9r = tau5r - tau6r; + tau9i = tau5i - tau6i; + + } + + + op[2].re = tau8r - tau9i; + op[2].im = tau8i + tau9r; + + op[6].re = tau8r + tau9i; + op[6].im = tau8i - tau9r; + + tau8r = tau4r - c1 * temp1r; + tau8i = tau4i - c1 * temp1i; + + //tau9r = sgn * ( -s1 * temp2r + tau7r); + //tau9i = sgn * ( -s1 * temp2i + tau7i); + + if (sgn == 1) { + tau9r = -s1 * temp2r + tau7r; + tau9i = -s1 * temp2i + tau7i; + + } else { + tau9r = s1 * temp2r - tau7r; + tau9i = s1 * temp2i - tau7i; + + } + + + op[3].re = tau8r - tau9i; + op[3].im = tau8i + tau9r; + + op[5].re = tau8r + tau9i; + op[5].im = tau8i - tau9r; + + + + + } else if (radix == 2) { + int k,tkm1,ind; + fft_type wlr,wli; + fft_type tau1r,tau1i,tau2r,tau2i; + m = N/2; + ll = 2*l; + mixed_radix_dit_rec(op,ip,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+m,ip+l,obj,sgn,m,ll,inc+1); + + for (k = 0; k < m;k++) { + ind = m-1+k; + wlr = (obj->twiddle+ind)->re; + wli = (obj->twiddle+ind)->im; + + tkm1 = k+m; + + tau1r = op[k].re; + tau1i = op[k].im; + + tau2r = op[tkm1].re * wlr - op[tkm1].im * wli; + tau2i = op[tkm1].im * wlr + op[tkm1].re * wli; + + op[k].re = tau1r + tau2r; + op[k].im = tau1i + tau2i; + + op[tkm1].re = tau1r - tau2r; + op[tkm1].im = tau1i - tau2i; + + + + } + + } else if (radix == 3) { + int k,tkm1,tkm2,ind; + fft_type wlr,wli,wl2r,wl2i; + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i; + fft_type ar,ai,br,bi,cr,ci; + m = N/3; + ll = 3*l; + mixed_radix_dit_rec(op,ip,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+m,ip+l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+2*m,ip+2*l,obj,sgn,m,ll,inc+1); + //printf("%d \n",inc); + //mixed_radix3_dit_rec(op,ip,obj,sgn,ll,m); + + + for (k = 0; k < m; ++k) { + ind = m-1+2*k; + wlr = (obj->twiddle+ind)->re; + wli = (obj->twiddle+ind)->im; + ind++; + wl2r = (obj->twiddle+ind)->re; + wl2i = (obj->twiddle+ind)->im; + tkm1 = k + m; + tkm2 = tkm1 + m; + + ar = op[k].re; + ai = op[k].im; + + br = op[tkm1].re * wlr - op[tkm1].im * wli; + bi = op[tkm1].im * wlr + op[tkm1].re * wli; + + cr = op[tkm2].re * wl2r - op[tkm2].im * wl2i; + ci = op[tkm2].im * wl2r + op[tkm2].re * wl2i; + + tau0r = br + cr; + tau0i = bi + ci; + + tau1r = sgn * 0.86602540378 * (br - cr); + tau1i = sgn * 0.86602540378 * (bi - ci); + + tau2r = ar - tau0r * 0.5000000000; + tau2i = ai - tau0i * 0.5000000000; + + + op[k].re = ar + tau0r; + op[k].im = ai + tau0i; + + op[tkm1].re = tau2r + tau1i; + op[tkm1].im = tau2i - tau1r; + + op[tkm2].re = tau2r - tau1i; + op[tkm2].im = tau2i + tau1r; + + } + + } else if (radix == 4) { + int k,tkm1,tkm2,tkm3,ind; + fft_type wlr,wli,wl2r,wl2i,wl3r,wl3i; + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i; + fft_type ar,ai,br,bi,cr,ci,dr,di; + m = N/4; + ll = 4*l; + mixed_radix_dit_rec(op,ip,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+m,ip+l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+2*m,ip+2*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+3*m,ip+3*l,obj,sgn,m,ll,inc+1); + + //mixed_radix4_dit_rec(op,ip,obj,sgn,ll,m); + + tkm1 = m; + tkm2 = tkm1 + m; + tkm3 = tkm2 + m; + + ar = op[0].re; + ai = op[0].im; + + br = op[tkm1].re; + bi = op[tkm1].im; + + cr = op[tkm2].re; + ci = op[tkm2].im; + + dr = op[tkm3].re; + di = op[tkm3].im; + + tau0r = ar + cr; + tau0i = ai + ci; + + tau1r = ar - cr; + tau1i = ai - ci; + + tau2r = br + dr; + tau2i = bi + di; + + tau3r = sgn* (br - dr); + tau3i = sgn* (bi - di); + + op[0].re = tau0r + tau2r ; + op[0].im = tau0i + tau2i; + + op[tkm1].re = tau1r + tau3i; + op[tkm1].im = tau1i - tau3r; + + op[tkm2].re = tau0r - tau2r; + op[tkm2].im = tau0i - tau2i; + + op[tkm3].re = tau1r - tau3i; + op[tkm3].im = tau1i + tau3r; + + + for (k = 1; k < m; k++) { + ind = m-1+3*k; + wlr = (obj->twiddle+ind)->re; + wli = (obj->twiddle+ind)->im; + ind++; + wl2r = (obj->twiddle+ind)->re; + wl2i = (obj->twiddle+ind)->im; + ind++; + wl3r = (obj->twiddle+ind)->re; + wl3i = (obj->twiddle+ind)->im; + + tkm1 = k+m; + tkm2 = tkm1 + m; + tkm3 = tkm2 + m; + + ar = op[k].re; + ai = op[k].im; + + br = op[tkm1].re * wlr - op[tkm1].im * wli; + bi = op[tkm1].im * wlr + op[tkm1].re * wli; + + cr = op[tkm2].re * wl2r - op[tkm2].im * wl2i; + ci = op[tkm2].im * wl2r + op[tkm2].re * wl2i; + + dr = op[tkm3].re * wl3r - op[tkm3].im * wl3i; + di = op[tkm3].im * wl3r + op[tkm3].re * wl3i; + + tau0r = ar + cr; + tau0i = ai + ci; + + tau1r = ar - cr; + tau1i = ai - ci; + + tau2r = br + dr; + tau2i = bi + di; + + tau3r = sgn* (br - dr); + tau3i = sgn* (bi - di); + + op[k].re = tau0r + tau2r ; + op[k].im = tau0i + tau2i; + + op[tkm1].re = tau1r + tau3i; + op[tkm1].im = tau1i - tau3r; + + op[tkm2].re = tau0r - tau2r; + op[tkm2].im = tau0i - tau2i; + + op[tkm3].re = tau1r - tau3i; + op[tkm3].im = tau1i + tau3r; + + } + + } else if (radix == 5) { + int k,tkm1,tkm2,tkm3,tkm4,ind; + fft_type wlr,wli,wl2r,wl2i,wl3r,wl3i,wl4r,wl4i; + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i; + fft_type ar,ai,br,bi,cr,ci,dr,di,er,ei; + fft_type tau4r,tau4i,tau5r,tau5i,tau6r,tau6i; + fft_type c1,c2,s1,s2; + m = N/5; + ll = 5*l; + mixed_radix_dit_rec(op,ip,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+m,ip+l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+2*m,ip+2*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+3*m,ip+3*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+4*m,ip+4*l,obj,sgn,m,ll,inc+1); + //printf("%d \n",inc); + //mixed_radix3_dit_rec(op,ip,obj,sgn,ll,m); + + c1 = 0.30901699437; + c2 = -0.80901699437; + s1 = 0.95105651629; + s2 = 0.58778525229; + + tkm1 = m; + tkm2 = tkm1 + m; + tkm3 = tkm2 + m; + tkm4 = tkm3 + m; + + ar = op[0].re; + ai = op[0].im; + + br = op[tkm1].re; + bi = op[tkm1].im; + + cr = op[tkm2].re; + ci = op[tkm2].im; + + dr = op[tkm3].re; + di = op[tkm3].im; + + er = op[tkm4].re; + ei = op[tkm4].im; + + tau0r = br + er; + tau0i = bi + ei; + + tau1r = cr + dr; + tau1i = ci + di; + + tau2r = br - er; + tau2i = bi - ei; + + tau3r = cr - dr; + tau3i = ci - di; + + op[0].re = ar + tau0r + tau1r; + op[0].im = ai + tau0i + tau1i; + + tau4r = c1 * tau0r + c2 * tau1r; + tau4i = c1 * tau0i + c2 * tau1i; + + tau5r = sgn * ( s1 * tau2r + s2 * tau3r); + tau5i = sgn * ( s1 * tau2i + s2 * tau3i); + + tau6r = ar + tau4r; + tau6i = ai + tau4i; + + op[tkm1].re = tau6r + tau5i; + op[tkm1].im = tau6i - tau5r; + + op[tkm4].re = tau6r - tau5i; + op[tkm4].im = tau6i + tau5r; + + tau4r = c2 * tau0r + c1 * tau1r; + tau4i = c2 * tau0i + c1 * tau1i; + + tau5r = sgn * ( s2 * tau2r - s1 * tau3r); + tau5i = sgn * ( s2 * tau2i - s1 * tau3i); + + tau6r = ar + tau4r; + tau6i = ai + tau4i; + + op[tkm2].re = tau6r + tau5i; + op[tkm2].im = tau6i - tau5r; + + op[tkm3].re = tau6r - tau5i; + op[tkm3].im = tau6i + tau5r; + + for (k = 1; k < m; k++) { + ind = m-1+4*k; + wlr = (obj->twiddle+ind)->re; + wli = (obj->twiddle+ind)->im; + ind++; + wl2r = (obj->twiddle+ind)->re; + wl2i = (obj->twiddle+ind)->im; + ind++; + wl3r = (obj->twiddle+ind)->re; + wl3i = (obj->twiddle+ind)->im; + ind++; + wl4r = (obj->twiddle+ind)->re; + wl4i = (obj->twiddle+ind)->im; + + tkm1 = k + m; + tkm2 = tkm1 + m; + tkm3 = tkm2 + m; + tkm4 = tkm3 + m; + + ar = op[k].re; + ai = op[k].im; + + br = op[tkm1].re * wlr - op[tkm1].im * wli; + bi = op[tkm1].im * wlr + op[tkm1].re * wli; + + cr = op[tkm2].re * wl2r - op[tkm2].im * wl2i; + ci = op[tkm2].im * wl2r + op[tkm2].re * wl2i; + + dr = op[tkm3].re * wl3r - op[tkm3].im * wl3i; + di = op[tkm3].im * wl3r + op[tkm3].re * wl3i; + + er = op[tkm4].re * wl4r - op[tkm4].im * wl4i; + ei = op[tkm4].im * wl4r + op[tkm4].re * wl4i; + + tau0r = br + er; + tau0i = bi + ei; + + tau1r = cr + dr; + tau1i = ci + di; + + tau2r = br - er; + tau2i = bi - ei; + + tau3r = cr - dr; + tau3i = ci - di; + + op[k].re = ar + tau0r + tau1r; + op[k].im = ai + tau0i + tau1i; + + tau4r = c1 * tau0r + c2 * tau1r; + tau4i = c1 * tau0i + c2 * tau1i; + + //tau5r = sgn * ( s1 * tau2r + s2 * tau3r); + //tau5i = sgn * ( s1 * tau2i + s2 * tau3i); + + if (sgn == 1) { + tau5r = s1 * tau2r + s2 * tau3r; + tau5i = s1 * tau2i + s2 * tau3i; + + } else { + tau5r = -s1 * tau2r - s2 * tau3r; + tau5i = -s1 * tau2i - s2 * tau3i; + + } + + tau6r = ar + tau4r; + tau6i = ai + tau4i; + + op[tkm1].re = tau6r + tau5i; + op[tkm1].im = tau6i - tau5r; + + op[tkm4].re = tau6r - tau5i; + op[tkm4].im = tau6i + tau5r; + + tau4r = c2 * tau0r + c1 * tau1r; + tau4i = c2 * tau0i + c1 * tau1i; + + //tau5r = sgn * ( s2 * tau2r - s1 * tau3r); + //tau5i = sgn * ( s2 * tau2i - s1 * tau3i); + + if (sgn == 1) { + tau5r = s2 * tau2r - s1 * tau3r; + tau5i = s2 * tau2i - s1 * tau3i; + + } else { + tau5r = -s2 * tau2r + s1 * tau3r; + tau5i = -s2 * tau2i + s1 * tau3i; + + } + + tau6r = ar + tau4r; + tau6i = ai + tau4i; + + op[tkm2].re = tau6r + tau5i; + op[tkm2].im = tau6i - tau5r; + + op[tkm3].re = tau6r - tau5i; + op[tkm3].im = tau6i + tau5r; + + } + + } else if (radix == 7) { + int k,tkm1,tkm2,tkm3,tkm4,tkm5,tkm6,ind; + fft_type wlr,wli,wl2r,wl2i,wl3r,wl3i,wl4r,wl4i,wl5r,wl5i,wl6r,wl6i; + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i; + fft_type ar,ai,br,bi,cr,ci,dr,di,er,ei,fr,fi,gr,gi; + fft_type tau4r,tau4i,tau5r,tau5i,tau6r,tau6i,tau7r,tau7i; + fft_type c1,c2,c3,s1,s2,s3; + m = N/7; + ll = 7*l; + mixed_radix_dit_rec(op,ip,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+m,ip+l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+2*m,ip+2*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+3*m,ip+3*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+4*m,ip+4*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+5*m,ip+5*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+6*m,ip+6*l,obj,sgn,m,ll,inc+1); + //printf("%d \n",inc); + //mixed_radix3_dit_rec(op,ip,obj,sgn,ll,m); + + c1 = 0.62348980185; + c2 = -0.22252093395; + c3 = -0.9009688679; + s1 = 0.78183148246; + s2 = 0.97492791218; + s3 = 0.43388373911; + + tkm1 = m; + tkm2 = tkm1 + m; + tkm3 = tkm2 + m; + tkm4 = tkm3 + m; + tkm5 = tkm4 + m; + tkm6 = tkm5 + m; + + ar = op[0].re; + ai = op[0].im; + + br = op[tkm1].re; + bi = op[tkm1].im; + + cr = op[tkm2].re; + ci = op[tkm2].im; + + dr = op[tkm3].re; + di = op[tkm3].im; + + er = op[tkm4].re; + ei = op[tkm4].im; + + fr = op[tkm5].re; + fi = op[tkm5].im; + + gr = op[tkm6].re; + gi = op[tkm6].im; + + tau0r = br + gr; + tau3r = br - gr; + tau0i = bi + gi; + tau3i = bi - gi; + + tau1r = cr + fr; + tau4r = cr - fr; + tau1i = ci + fi; + tau4i = ci - fi; + + tau2r = dr + er; + tau5r = dr - er; + tau2i = di + ei; + tau5i = di - ei; + + op[0].re = ar + tau0r + tau1r + tau2r; + op[0].im = ai + tau0i + tau1i + tau2i; + + tau6r = ar + c1 * tau0r + c2 * tau1r + c3 * tau2r; + tau6i = ai + c1 * tau0i + c2 * tau1i + c3 * tau2i; + + //tau7r = sgn * ( -s1 * tau3r - s2 * tau4r - s3 * tau5r); + //tau7i = sgn * ( -s1 * tau3i - s2 * tau4i - s3 * tau5i); + + if (sgn == 1) { + tau7r = -s1 * tau3r - s2 * tau4r - s3 * tau5r; + tau7i = -s1 * tau3i - s2 * tau4i - s3 * tau5i; + + } else { + tau7r = s1 * tau3r + s2 * tau4r + s3 * tau5r; + tau7i = s1 * tau3i + s2 * tau4i + s3 * tau5i; + + } + + + op[tkm1].re = tau6r - tau7i; + op[tkm1].im = tau6i + tau7r; + + op[tkm6].re = tau6r + tau7i; + op[tkm6].im = tau6i - tau7r; + + tau6r = ar + c2 * tau0r + c3 * tau1r + c1 * tau2r; + tau6i = ai + c2 * tau0i + c3 * tau1i + c1 * tau2i; + + //tau7r = sgn * ( -s2 * tau3r + s3 * tau4r + s1 * tau5r); + //tau7i = sgn * ( -s2 * tau3i + s3 * tau4i + s1 * tau5i); + + if (sgn == 1) { + tau7r = -s2 * tau3r + s3 * tau4r + s1 * tau5r; + tau7i = -s2 * tau3i + s3 * tau4i + s1 * tau5i; + + } else { + tau7r = s2 * tau3r - s3 * tau4r - s1 * tau5r; + tau7i = s2 * tau3i - s3 * tau4i - s1 * tau5i; + + } + + + op[tkm2].re = tau6r - tau7i; + op[tkm2].im = tau6i + tau7r; + + op[tkm5].re = tau6r + tau7i; + op[tkm5].im = tau6i - tau7r; + + tau6r = ar + c3 * tau0r + c1 * tau1r + c2 * tau2r; + tau6i = ai + c3 * tau0i + c1 * tau1i + c2 * tau2i; + + //tau7r = sgn * ( -s3 * tau3r + s1 * tau4r - s2 * tau5r); + //tau7i = sgn * ( -s3 * tau3i + s1 * tau4i - s2 * tau5i); + + if (sgn == 1) { + tau7r = -s3 * tau3r + s1 * tau4r - s2 * tau5r; + tau7i = -s3 * tau3i + s1 * tau4i - s2 * tau5i; + + } else { + tau7r = s3 * tau3r - s1 * tau4r + s2 * tau5r; + tau7i = s3 * tau3i - s1 * tau4i + s2 * tau5i; + + } + + + op[tkm3].re = tau6r - tau7i; + op[tkm3].im = tau6i + tau7r; + + op[tkm4].re = tau6r + tau7i; + op[tkm4].im = tau6i - tau7r; + + + for (k = 1; k < m; k++) { + ind = m-1+6*k; + wlr = (obj->twiddle+ind)->re; + wli = (obj->twiddle+ind)->im; + ind++; + wl2r = (obj->twiddle+ind)->re; + wl2i = (obj->twiddle+ind)->im; + ind++; + wl3r = (obj->twiddle+ind)->re; + wl3i = (obj->twiddle+ind)->im; + ind++; + wl4r = (obj->twiddle+ind)->re; + wl4i = (obj->twiddle+ind)->im; + ind++; + wl5r = (obj->twiddle+ind)->re; + wl5i = (obj->twiddle+ind)->im; + ind++; + wl6r = (obj->twiddle+ind)->re; + wl6i = (obj->twiddle+ind)->im; + + tkm1 = k + m; + tkm2 = tkm1 + m; + tkm3 = tkm2 + m; + tkm4 = tkm3 + m; + tkm5 = tkm4 + m; + tkm6 = tkm5 + m; + + ar = op[k].re; + ai = op[k].im; + + br = op[tkm1].re * wlr - op[tkm1].im * wli; + bi = op[tkm1].im * wlr + op[tkm1].re * wli; + + cr = op[tkm2].re * wl2r - op[tkm2].im * wl2i; + ci = op[tkm2].im * wl2r + op[tkm2].re * wl2i; + + dr = op[tkm3].re * wl3r - op[tkm3].im * wl3i; + di = op[tkm3].im * wl3r + op[tkm3].re * wl3i; + + er = op[tkm4].re * wl4r - op[tkm4].im * wl4i; + ei = op[tkm4].im * wl4r + op[tkm4].re * wl4i; + + fr = op[tkm5].re * wl5r - op[tkm5].im * wl5i; + fi = op[tkm5].im * wl5r + op[tkm5].re * wl5i; + + gr = op[tkm6].re * wl6r - op[tkm6].im * wl6i; + gi = op[tkm6].im * wl6r + op[tkm6].re * wl6i; + + tau0r = br + gr; + tau3r = br - gr; + tau0i = bi + gi; + tau3i = bi - gi; + + tau1r = cr + fr; + tau4r = cr - fr; + tau1i = ci + fi; + tau4i = ci - fi; + + tau2r = dr + er; + tau5r = dr - er; + tau2i = di + ei; + tau5i = di - ei; + + op[k].re = ar + tau0r + tau1r + tau2r; + op[k].im = ai + tau0i + tau1i + tau2i; + + tau6r = ar + c1 * tau0r + c2 * tau1r + c3 * tau2r; + tau6i = ai + c1 * tau0i + c2 * tau1i + c3 * tau2i; + + //tau7r = sgn * ( -s1 * tau3r - s2 * tau4r - s3 * tau5r); + //tau7i = sgn * ( -s1 * tau3i - s2 * tau4i - s3 * tau5i); + + if (sgn == 1) { + tau7r = -s1 * tau3r - s2 * tau4r - s3 * tau5r; + tau7i = -s1 * tau3i - s2 * tau4i - s3 * tau5i; + + } else { + tau7r = s1 * tau3r + s2 * tau4r + s3 * tau5r; + tau7i = s1 * tau3i + s2 * tau4i + s3 * tau5i; + + } + + + op[tkm1].re = tau6r - tau7i; + op[tkm1].im = tau6i + tau7r; + + op[tkm6].re = tau6r + tau7i; + op[tkm6].im = tau6i - tau7r; + + tau6r = ar + c2 * tau0r + c3 * tau1r + c1 * tau2r; + tau6i = ai + c2 * tau0i + c3 * tau1i + c1 * tau2i; + + //tau7r = sgn * ( -s2 * tau3r + s3 * tau4r + s1 * tau5r); + //tau7i = sgn * ( -s2 * tau3i + s3 * tau4i + s1 * tau5i); + + if (sgn == 1) { + tau7r = -s2 * tau3r + s3 * tau4r + s1 * tau5r; + tau7i = -s2 * tau3i + s3 * tau4i + s1 * tau5i; + + } else { + tau7r = s2 * tau3r - s3 * tau4r - s1 * tau5r; + tau7i = s2 * tau3i - s3 * tau4i - s1 * tau5i; + + } + + + op[tkm2].re = tau6r - tau7i; + op[tkm2].im = tau6i + tau7r; + + op[tkm5].re = tau6r + tau7i; + op[tkm5].im = tau6i - tau7r; + + tau6r = ar + c3 * tau0r + c1 * tau1r + c2 * tau2r; + tau6i = ai + c3 * tau0i + c1 * tau1i + c2 * tau2i; + + //tau7r = sgn * ( -s3 * tau3r + s1 * tau4r - s2 * tau5r); + //tau7i = sgn * ( -s3 * tau3i + s1 * tau4i - s2 * tau5i); + + if (sgn == 1) { + tau7r = -s3 * tau3r + s1 * tau4r - s2 * tau5r; + tau7i = -s3 * tau3i + s1 * tau4i - s2 * tau5i; + + } else { + tau7r = s3 * tau3r - s1 * tau4r + s2 * tau5r; + tau7i = s3 * tau3i - s1 * tau4i + s2 * tau5i; + + } + + + op[tkm3].re = tau6r - tau7i; + op[tkm3].im = tau6i + tau7r; + + op[tkm4].re = tau6r + tau7i; + op[tkm4].im = tau6i - tau7r; + + } + + } else if (radix == 8) { + int k,tkm1,tkm2,tkm3,tkm4,tkm5,tkm6,tkm7,ind; + fft_type wlr,wli,wl2r,wl2i,wl3r,wl3i,wl4r,wl4i,wl5r,wl5i,wl6r,wl6i,wl7r,wl7i; + fft_type tau0r,tau0i,tau1r,tau1i,tau2r,tau2i,tau3r,tau3i; + fft_type ar,ai,br,bi,cr,ci,dr,di,er,ei,fr,fi,gr,gi,hr,hi; + fft_type tau4r,tau4i,tau5r,tau5i,tau6r,tau6i,tau7r,tau7i,tau8r,tau8i,tau9r,tau9i; + fft_type c1,s1,temp1r,temp1i,temp2r,temp2i; + m = N/8; + ll = 8*l; + mixed_radix_dit_rec(op,ip,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+m,ip+l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+2*m,ip+2*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+3*m,ip+3*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+4*m,ip+4*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+5*m,ip+5*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+6*m,ip+6*l,obj,sgn,m,ll,inc+1); + mixed_radix_dit_rec(op+7*m,ip+7*l,obj,sgn,m,ll,inc+1); + //printf("%d \n",inc); + //mixed_radix3_dit_rec(op,ip,obj,sgn,ll,m); + + c1 = 0.70710678118654752440084436210485; + s1 = 0.70710678118654752440084436210485; + + + for (k = 0; k < m; k++) { + ind = m-1+7*k; + wlr = (obj->twiddle+ind)->re; + wli = (obj->twiddle+ind)->im; + ind++; + wl2r = (obj->twiddle+ind)->re; + wl2i = (obj->twiddle+ind)->im; + ind++; + wl3r = (obj->twiddle+ind)->re; + wl3i = (obj->twiddle+ind)->im; + ind++; + wl4r = (obj->twiddle+ind)->re; + wl4i = (obj->twiddle+ind)->im; + ind++; + wl5r = (obj->twiddle+ind)->re; + wl5i = (obj->twiddle+ind)->im; + ind++; + wl6r = (obj->twiddle+ind)->re; + wl6i = (obj->twiddle+ind)->im; + ind++; + wl7r = (obj->twiddle+ind)->re; + wl7i = (obj->twiddle+ind)->im; + + tkm1 = k + m; + tkm2 = tkm1 + m; + tkm3 = tkm2 + m; + tkm4 = tkm3 + m; + tkm5 = tkm4 + m; + tkm6 = tkm5 + m; + tkm7 = tkm6 + m; + + ar = op[k].re; + ai = op[k].im; + + br = op[tkm1].re * wlr - op[tkm1].im * wli; + bi = op[tkm1].im * wlr + op[tkm1].re * wli; + + cr = op[tkm2].re * wl2r - op[tkm2].im * wl2i; + ci = op[tkm2].im * wl2r + op[tkm2].re * wl2i; + + dr = op[tkm3].re * wl3r - op[tkm3].im * wl3i; + di = op[tkm3].im * wl3r + op[tkm3].re * wl3i; + + er = op[tkm4].re * wl4r - op[tkm4].im * wl4i; + ei = op[tkm4].im * wl4r + op[tkm4].re * wl4i; + + fr = op[tkm5].re * wl5r - op[tkm5].im * wl5i; + fi = op[tkm5].im * wl5r + op[tkm5].re * wl5i; + + gr = op[tkm6].re * wl6r - op[tkm6].im * wl6i; + gi = op[tkm6].im * wl6r + op[tkm6].re * wl6i; + + hr = op[tkm7].re * wl7r - op[tkm7].im * wl7i; + hi = op[tkm7].im * wl7r + op[tkm7].re * wl7i; + + tau0r = ar + er; + tau4r = ar - er; + tau0i = ai + ei; + tau4i = ai - ei; + + tau1r = br + hr; + tau5r = br - hr; + tau1i = bi + hi; + tau5i = bi - hi; + + tau2r = dr + fr; + tau6r = dr - fr; + tau6i = di - fi; + tau2i = di + fi; + + tau3r = cr + gr; + tau7r = cr - gr; + tau7i = ci - gi; + tau3i = ci + gi; + + op[k].re = tau0r + tau1r + tau2r + tau3r; + op[k].im = tau0i + tau1i + tau2i + tau3i; + + op[tkm4].re = tau0r - tau1r - tau2r + tau3r; + op[tkm4].im = tau0i - tau1i - tau2i + tau3i; + + temp1r = tau1r - tau2r; + temp1i = tau1i - tau2i; + + temp2r = tau5r + tau6r; + temp2i = tau5i + tau6i; + + tau8r = tau4r + c1 * temp1r; + tau8i = tau4i + c1 * temp1i; + + //tau9r = sgn * ( -s1 * temp2r - tau7r); + //tau9i = sgn * ( -s1 * temp2i - tau7i); + + if (sgn == 1) { + tau9r = -s1 * temp2r - tau7r; + tau9i = -s1 * temp2i - tau7i; + + } else { + tau9r = s1 * temp2r + tau7r; + tau9i = s1 * temp2i + tau7i; + + } + + + op[tkm1].re = tau8r - tau9i; + op[tkm1].im = tau8i + tau9r; + + op[tkm7].re = tau8r + tau9i; + op[tkm7].im = tau8i - tau9r; + + tau8r = tau0r - tau3r; + tau8i = tau0i - tau3i; + + //tau9r = sgn * ( -tau5r + tau6r); + //tau9i = sgn * ( -tau5i + tau6i); + + if (sgn == 1) { + tau9r = -tau5r + tau6r; + tau9i = -tau5i + tau6i; + + } else { + tau9r = tau5r - tau6r; + tau9i = tau5i - tau6i; + + } + + + op[tkm2].re = tau8r - tau9i; + op[tkm2].im = tau8i + tau9r; + + op[tkm6].re = tau8r + tau9i; + op[tkm6].im = tau8i - tau9r; + + tau8r = tau4r - c1 * temp1r; + tau8i = tau4i - c1 * temp1i; + + //tau9r = sgn * ( -s1 * temp2r + tau7r); + //tau9i = sgn * ( -s1 * temp2i + tau7i); + + if (sgn == 1) { + tau9r = -s1 * temp2r + tau7r; + tau9i = -s1 * temp2i + tau7i; + + } else { + tau9r = s1 * temp2r - tau7r; + tau9i = s1 * temp2i - tau7i; + + } + + + op[tkm3].re = tau8r - tau9i; + op[tkm3].im = tau8i + tau9r; + + op[tkm5].re = tau8r + tau9i; + op[tkm5].im = tau8i - tau9r; + + } + + } else { + int k,i,ind; + int M,tkm,u,v,t,tt; + fft_type temp1r,temp1i,temp2r,temp2i; + fft_type* wlr = (fft_type*) malloc (sizeof(fft_type) * (radix-1)); + fft_type* wli = (fft_type*) malloc (sizeof(fft_type) * (radix-1)); + fft_type* taur = (fft_type*) malloc (sizeof(fft_type) * (radix-1)); + fft_type* taui = (fft_type*) malloc (sizeof(fft_type) * (radix-1)); + fft_type* c1 = (fft_type*) malloc (sizeof(fft_type) * (radix-1)); + fft_type* s1 = (fft_type*) malloc (sizeof(fft_type) * (radix-1)); + fft_type* yr = (fft_type*) malloc (sizeof(fft_type) * (radix)); + fft_type* yi = (fft_type*) malloc (sizeof(fft_type) * (radix)); + + m = N /radix; + ll = radix * l; + + for (i = 0; i < radix;++i) { + mixed_radix_dit_rec(op+i*m,ip+i*l,obj,sgn,m,ll,inc+1); + } + + M = (radix - 1 )/2; + + for (i = 1; i < M+1;++i) { + c1[i-1] = cos(i*PI2/radix); + s1[i-1] = sin(i*PI2/radix); + } + + for (i = 0; i < M;++i) { + s1[i+M] = -s1[M-1-i]; + c1[i+M] = c1[M-1-i]; + } + + for (k = 0; k < m;++k) { + ind = m-1+(radix-1)*k; + yr[0] = op[k].re; + yi[0] = op[k].im; + for(i =0; i < radix -1;++i) { + wlr[i] = (obj->twiddle+ind)->re; + wli[i] = (obj->twiddle+ind)->im; + tkm = k + (i+1)*m; + yr[i+1] = op[tkm].re * wlr[i] - op[tkm].im * wli[i]; + yi[i+1] = op[tkm].im * wlr[i] + op[tkm].re * wli[i]; + ind++; + } + + for (i = 0; i < M; ++i) { + taur[i] = yr[i+1] + yr[radix-1-i]; + taui[i+M] = yi[i+1] - yi[radix-1-i]; + taui[i] = yi[i+1] + yi[radix-1-i]; + taur[i+M] = yr[i+1] - yr[radix-1-i]; + } + + temp1r = yr[0]; + temp1i = yi[0]; + + for (i = 0; i < M; ++i) { + temp1r+= taur[i]; + temp1i+= taui[i]; + } + + op[k].re = temp1r; + op[k].im = temp1i; + + for (u = 0; u < M; u++) { + temp1r = yr[0]; + temp1i = yi[0]; + temp2r = 0.0; + temp2i = 0.0; + for (v = 0; v < M; v++) { + //int ind2 = (u+v)%M; + t = (u+1)*(v+1); + while(t >= radix) + t-=radix; + tt = t-1; + + temp1r+= c1[tt]*taur[v]; + temp1i+= c1[tt]*taui[v]; + temp2r-= s1[tt]*taur[v+M]; + temp2i-= s1[tt]*taui[v+M]; + } + temp2r = sgn * temp2r; + temp2i = sgn * temp2i; + + + op[k + (u+1)*m].re = temp1r - temp2i; + op[k + (u+1)*m].im = temp1i + temp2r; + + op[k + (radix-u-1)*m].re = temp1r + temp2i; + op[k + (radix-u-1)*m].im = temp1i - temp2r; + } + + + } + free(wlr); + free(wli); + free(taur); + free(taui); + free(c1); + free(s1); + free(yr); + free(yi); + + } + + +} + +static void bluestein_exp(fft_data *hl, fft_data *hlt, int len, int M) { + fft_type PI,theta,angle; + int l2,len2,i; + PI = 3.1415926535897932384626433832795; + theta = PI / len; + l2 = 0; + len2 = 2 * len; + + for (i = 0 ; i < len; ++i) { + angle = theta * l2; + hlt[i].re = cos(angle); + hlt[i].im = sin(angle); + hl[i].re = hlt[i].re; + hl[i].im = hlt[i].im; + l2+=2*i+1; + while (l2 > len2) { + l2-=len2; + } + + } + + for (i = len; i < M-len+1; i++) { + hl[i].re = 0.0; + hl[i].im = 0.0; + } + + for (i = M - len + 1; i < M; i++) { + hl[i].re = hlt[M-i].re; + hl[i].im = hlt[M-i].im; + } + +} + +static void bluestein_fft(fft_data *data, fft_data *oup,fft_object obj,int sgn, int N) { + + int K,M,ii,i; + int def_lt,def_N,def_sgn; + fft_type scale,temp; + fft_data* yn; + fft_data* hk; + fft_data* tempop; + fft_data* yno; + fft_data* hlt; + obj->lt = 0; + K = (int) pow(2.0,ceil((double) log10((double) N)/log10((double) 2.0))); + def_lt = 1; + def_sgn = obj->sgn; + def_N = obj->N; + + if (K < 2 * N - 2) { + M = K * 2; + } else { + M = K; + } + obj->N = M; + + yn = (fft_data*) malloc (sizeof(fft_data) * M); + hk = (fft_data*) malloc (sizeof(fft_data) * M); + tempop = (fft_data*) malloc (sizeof(fft_data) * M); + yno = (fft_data*) malloc (sizeof(fft_data) * M); + hlt = (fft_data*) malloc (sizeof(fft_data) * N); + //fft_data* twi = (fft_data*) malloc (sizeof(fft_data) * M); + + bluestein_exp(tempop,hlt,N,M); + scale = 1.0/M; + + for (ii = 0; ii < M; ++ii) { + tempop[ii].im *= scale; + tempop[ii].re *= scale; + } + + //fft_object obj = initialize_fft2(M,1); + fft_exec(obj,tempop,hk); + + if (sgn == 1) { + for (i = 0; i < N; i++) { + tempop[i].re = data[i].re * hlt[i].re + data[i].im * hlt[i].im; + tempop[i].im = -data[i].re * hlt[i].im + data[i].im * hlt[i].re; + } + } else { + for (i = 0; i < N; i++) { + tempop[i].re = data[i].re * hlt[i].re - data[i].im * hlt[i].im; + tempop[i].im = data[i].re * hlt[i].im + data[i].im * hlt[i].re; + } + } + + for (i = N;i < M;i++) { + tempop[i].re = 0.0; + tempop[i].im = 0.0; + } + + fft_exec(obj,tempop,yn); + + if (sgn == 1) { + for (i = 0; i < M; i++) { + temp = yn[i].re * hk[i].re - yn[i].im * hk[i].im; + yn[i].im = yn[i].re * hk[i].im + yn[i].im * hk[i].re; + yn[i].re = temp; + } + } else { + for (i = 0; i < M; i++) { + temp = yn[i].re * hk[i].re + yn[i].im * hk[i].im; + yn[i].im = -yn[i].re * hk[i].im + yn[i].im * hk[i].re; + yn[i].re = temp; + } + + } + + //IFFT + + for (ii = 0; ii < M; ++ii) { + (obj->twiddle+ii)->im = -(obj->twiddle+ii)->im; + } + + obj->sgn = -1*sgn; + + fft_exec(obj,yn,yno); + + if (sgn == 1) { + for (i = 0; i < N; i++) { + oup[i].re = yno[i].re * hlt[i].re + yno[i].im * hlt[i].im; + oup[i].im = -yno[i].re * hlt[i].im + yno[i].im * hlt[i].re; + } + } else { + for (i = 0; i < N; i++) { + oup[i].re = yno[i].re * hlt[i].re - yno[i].im * hlt[i].im; + oup[i].im = yno[i].re * hlt[i].im + yno[i].im * hlt[i].re; + } + + } + + obj->sgn = def_sgn; + obj->N = def_N; + obj->lt = def_lt; + for (ii = 0; ii < M; ++ii) { + (obj->twiddle+ii)->im = -(obj->twiddle+ii)->im; + } + + free(yn); + free(yno); + free(tempop); + free(hk); + free(hlt); + +} + + + +void fft_exec(fft_object obj,fft_data *inp,fft_data *oup) { + if (obj->lt == 0) { + //fftct_radix3_dit_rec(inp,oup,obj, obj->sgn, obj->N); + //fftct_mixed_rec(inp,oup,obj, obj->sgn, obj->N); + //printf("%f \n", 1.785); + int l,inc; + int nn,sgn1; + nn = obj->N; + sgn1 = obj->sgn; + l = 1; + inc = 0; + //radix3_dit_rec(oup,inp,obj,sgn1,nn,l); + mixed_radix_dit_rec(oup,inp,obj,sgn1,nn,l,inc); + } else if (obj->lt == 1){ + //printf("%f \n", 1.785); + int nn,sgn1; + nn = obj->N; + sgn1 = obj->sgn; + bluestein_fft(inp,oup,obj,sgn1,nn); + + } + +} + +int divideby(int M,int d) { + while (M%d == 0) { + M = M/d; + } + if (M == 1) { + return 1; + } + return 0; +} + +int dividebyN(int N) { + while (N%53 == 0) { + N = N/53; + } + while (N%47 == 0) { + N = N/47; + } + while (N%43 == 0) { + N = N/43; + } + while (N%41 == 0) { + N = N/41; + } + while (N%37 == 0) { + N = N/37; + } + while (N%31 == 0) { + N = N/31; + } + while (N%29 == 0) { + N = N/29; + } + while (N%23 == 0) { + N = N/23; + } + while (N%17 == 0) { + N = N/17; + } + while (N%13 == 0) { + N = N/13; + } + while (N%11 == 0) { + N = N/11; + } + while (N%8 == 0) { + N = N/8; + } + while (N%7 == 0) { + N = N/7; + } + while (N%5 == 0) { + N = N/5; + } + while (N%4 == 0) { + N = N/4; + } + while (N%3 == 0) { + N = N/3; + } + while (N%2 == 0) { + N = N/2; + } + if (N == 1) { + return 1; + } + return 0; + +} + +int factors(int M, int* arr) { + int i,N,num,mult,m1,m2; + i = 0; + N = M; + while (N%53 == 0) { + N = N/53; + arr[i] = 53; + i++; + } + while (N%47 == 0) { + N = N/47; + arr[i] = 47; + i++; + } + while (N%43 == 0) { + N = N/43; + arr[i] = 43; + i++; + } + while (N%41 == 0) { + N = N/41; + arr[i] = 41; + i++; + } + while (N%37 == 0) { + N = N/37; + arr[i] = 37; + i++; + } + while (N%31 == 0) { + N = N/31; + arr[i] = 31; + i++; + } + while (N%29 == 0) { + N = N/29; + arr[i] = 29; + i++; + } + while (N%23 == 0) { + N = N/23; + arr[i] = 23; + i++; + } + while (N%19 == 0) { + N = N/19; + arr[i] = 19; + i++; + } + while (N%17 == 0) { + N = N/17; + arr[i] = 17; + i++; + } + while (N%13 == 0) { + N = N/13; + arr[i] = 13; + i++; + } + while (N%11 == 0) { + N = N/11; + arr[i] = 11; + i++; + } + while (N%8 == 0) { + N = N/8; + arr[i] = 8; + i++; + } + while (N%7 == 0) { + N = N/7; + arr[i] = 7; + i++; + } + while (N%5 == 0) { + N = N/5; + arr[i] = 5; + i++; + } + while (N%4 == 0) { + N = N/4; + arr[i] = 4; + i++; + } + while (N%3 == 0) { + N = N/3; + arr[i] = 3; + i++; + } + while (N%2 == 0) { + N = N/2; + arr[i] = 2; + i++; + } + if (N > 31) { + num = 2; + + while (N > 1) { + mult = num*6; + m1 = mult-1; + m2 = mult+1; + while (N%m1 == 0 ) { + arr[i] = m1; + i++; + N = N / m1; + } + while (N%m2 == 0 ) { + arr[i] = m2; + i++; + N = N / m2; + } + num+=1; + + } + } + return i; + +} + + +void twiddle(fft_data *vec,int N, int radix) { + int K,KL; + fft_type theta,theta2; + theta = PI2/N; + KL = N/radix; + vec[0].re = 1.0; + vec[0].im = 0.0; + + for (K = 1; K < KL;K++) { + theta2 = theta * K; + vec[K].re = cos(theta2); + vec[K].im = -sin(theta2); + } + +} + +void longvectorN(fft_data *sig,int N, int *array, int tx) { + int L,i,Ls,ct,j,k; + fft_type theta; + L = 1; + ct = 0; + for (i = 0; i < tx; i++) { + L = L * array[tx-1-i]; + Ls = L / array[tx-1-i]; + theta = -1.0 * PI2/L; + for (j = 0; j < Ls;j++) { + for (k = 0; k < array[tx-1-i] -1 ;k++) { + sig[ct].re = cos((k+1)*j*theta); + sig[ct].im = sin((k+1)*j*theta); + ct++; + } + } + + } + +} + + + + +void free_fft(fft_object object) { + free(object); +} diff --git a/src/hsfft.h b/src/hsfft.h new file mode 100644 index 0000000..95d6b71 --- /dev/null +++ b/src/hsfft.h @@ -0,0 +1,74 @@ +/* + * hsfft.h + * + * Created on: Apr 14, 2013 + * Author: Rafat Hussain + */ + +#ifndef HSFFT_H_ +#define HSFFT_H_ + +#include +#include +#include + +#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_ */ diff --git a/src/real.c b/src/real.c new file mode 100644 index 0000000..c37dbc7 --- /dev/null +++ b/src/real.c @@ -0,0 +1,111 @@ +/* + * real.c + * + * Created on: Apr 20, 2013 + * Author: Rafat Hussain + */ +#include +#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); +} + diff --git a/src/real.h b/src/real.h new file mode 100644 index 0000000..42d30eb --- /dev/null +++ b/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_ */ diff --git a/src/wavefilt.c b/src/wavefilt.c new file mode 100644 index 0000000..71b39de --- /dev/null +++ b/src/wavefilt.c @@ -0,0 +1,3225 @@ +#include "wavefilt.h" + +int filtlength(char* name) { + if (!strcmp(name,"haar") || !strcmp(name,"db1")) { + return 2; + } + else if (!strcmp(name,"db2")){ + return 4; + } + else if (!strcmp(name,"db3")){ + return 6; + } + else if (!strcmp(name,"db4")){ + return 8; + } + else if (!strcmp(name,"db5")){ + return 10; + } + else if (!strcmp(name,"db6")){ + return 12; + } + else if (!strcmp(name,"db7")){ + return 14; + } + else if (!strcmp(name,"db8")){ + return 16; + } + else if (!strcmp(name,"db9")){ + return 18; + } + + else if (!strcmp(name,"db10")){ + return 20; + } + + else if (!strcmp(name,"db12")){ + return 24; + } + else if (!strcmp(name,"db13")){ + return 26; + } + + else if (!strcmp(name,"db11")){ + return 22; + } + + else if (!strcmp(name,"db14")){ + return 28; + } + else if (!strcmp(name,"db15")){ + return 30; + } + else if (!strcmp(name,"bior1.1")){ + return 2; + } + + else if (!strcmp(name,"bior1.3")){ + return 6; + } + + else if (!strcmp(name,"bior1.5")){ + return 10; + } + + else if (!strcmp(name,"bior2.2")){ + return 6; + } + + else if (!strcmp(name,"bior2.4")){ + return 10; + } + + else if (!strcmp(name,"bior2.6")){ + return 14; + } + else if (!strcmp(name,"bior2.8")){ + return 18; + } + + else if (!strcmp(name,"bior3.1")){ + return 8; + } + else if (!strcmp(name,"bior3.5")){ + return 12; + } + + else if (!strcmp(name,"bior3.7")){ + return 16; + } + else if (!strcmp(name,"bior3.9")){ + return 20; + } + else if (!strcmp(name,"bior4.4")){ + return 10; + } + else if (!strcmp(name,"bior5.5")){ + return 12; + } + else if (!strcmp(name,"bior6.8")){ + return 18; + } + + else if (!strcmp(name,"coif1")){ + return 6; + } + else if (!strcmp(name,"coif2")){ + return 12; + } + else if (!strcmp(name,"coif3")){ + return 18; + } + else if (!strcmp(name,"coif4")){ + return 24; + } + else if (!strcmp(name,"coif5")){ + return 30; + } + + else if (!strcmp(name,"sym2")){ + return 4; + } + + else if (!strcmp(name,"sym3")){ + return 6; + } + + else if (!strcmp(name,"sym4")){ + return 8; + } + + else if (!strcmp(name,"sym5")){ + return 10; + } + + else if (!strcmp(name,"sym6")){ + return 12; + } + + else if (!strcmp(name,"sym7")){ + return 14; + } + + else if (!strcmp(name,"sym8")){ + return 16; + } + + else if (!strcmp(name,"sym9")){ + return 18; + } + + else if (!strcmp(name,"sym10")){ + return 20; + } + else { + printf("\n Filter Not in Database \n"); + return -1; + } + + return 0; +} + +int filtcoef(char* name, double *lp1, double *hp1, double *lp2, double *hp2) { + int i; + if (!strcmp(name,"haar") || !strcmp(name,"db1")) { + double lp1_a[] = { 0.7071, 0.7071 }; + double hp1_a[] = { -0.7071, 0.7071 }; + double lp2_a[] = { 0.7071, 0.7071 }; + double hp2_a[] = { 0.7071, -0.7071 }; + for (i = 0; i < 2; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 2; + } + else if (!strcmp(name,"db2")){ + double lp1_a[] = { -0.12940952255092145, 0.22414386804185735, 0.83651630373746899, + 0.48296291314469025 }; + + double hp1_a[] = { -0.48296291314469025, 0.83651630373746899, -0.22414386804185735, + -0.12940952255092145 }; + + double lp2_a[] = { 0.48296291314469025, 0.83651630373746899, 0.22414386804185735, + -0.12940952255092145 }; + + double hp2_a[] = { -0.12940952255092145, -0.22414386804185735, 0.83651630373746899, + -0.48296291314469025 }; + for (i = 0; i < 4; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 4; + } + else if (!strcmp(name,"db3")){ + double lp1_a[] = { 0.035226291882100656, -0.085441273882241486, -0.13501102001039084, + 0.45987750211933132, 0.80689150931333875, 0.33267055295095688 }; + + double hp1_a[] = { -0.33267055295095688, 0.80689150931333875, -0.45987750211933132, + -0.13501102001039084, 0.085441273882241486, 0.035226291882100656 }; + + double lp2_a[] = { 0.33267055295095688, 0.80689150931333875, 0.45987750211933132, + -0.13501102001039084, -0.085441273882241486, 0.035226291882100656 }; + + double hp2_a[] = { 0.035226291882100656, 0.085441273882241486, -0.13501102001039084, + -0.45987750211933132, 0.80689150931333875, -0.33267055295095688 }; + for (i = 0; i < 6; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 6; + } + else if (!strcmp(name,"db4")){ + double lp1_a[] = { -0.010597401784997278, 0.032883011666982945, 0.030841381835986965, + -0.18703481171888114, -0.027983769416983849, 0.63088076792959036, + 0.71484657055254153, 0.23037781330885523 }; + + double hp1_a[] = { -0.23037781330885523, 0.71484657055254153, -0.63088076792959036, + -0.027983769416983849, 0.18703481171888114, 0.030841381835986965, + -0.032883011666982945, -0.010597401784997278 }; + + double lp2_a[] = { 0.23037781330885523, 0.71484657055254153, 0.63088076792959036, + -0.027983769416983849, -0.18703481171888114, 0.030841381835986965, + 0.032883011666982945, -0.010597401784997278 }; + + double hp2_a[] = { -0.010597401784997278, -0.032883011666982945, 0.030841381835986965, + 0.18703481171888114, -0.027983769416983849, -0.63088076792959036, + 0.71484657055254153, -0.23037781330885523 }; + for (i = 0; i < 8; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 8; + } + else if (!strcmp(name,"db5")){ + double lp1_a[] = { 0.0033357252850015492, -0.012580751999015526, -0.0062414902130117052, + 0.077571493840065148, -0.03224486958502952, -0.24229488706619015, + 0.13842814590110342, 0.72430852843857441, 0.60382926979747287, + 0.16010239797412501 }; + + double hp1_a[] = { -0.16010239797412501, 0.60382926979747287, -0.72430852843857441, + 0.13842814590110342, 0.24229488706619015, -0.03224486958502952, + -0.077571493840065148, -0.0062414902130117052, 0.012580751999015526, + 0.0033357252850015492 }; + + double lp2_a[] = { 0.16010239797412501, 0.60382926979747287, 0.72430852843857441, + 0.13842814590110342, -0.24229488706619015, -0.03224486958502952, + 0.077571493840065148, -0.0062414902130117052, -0.012580751999015526, + 0.0033357252850015492 }; + + double hp2_a[] = { 0.0033357252850015492, 0.012580751999015526, -0.0062414902130117052, + -0.077571493840065148, -0.03224486958502952, 0.24229488706619015, + 0.13842814590110342, -0.72430852843857441, 0.60382926979747287, + -0.16010239797412501 }; + for (i = 0; i < 10; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 10; + } + else if (!strcmp(name,"db6")){ + double lp1_a[] = { -0.0010773010849955799, + 0.0047772575110106514, + 0.0005538422009938016, + -0.031582039318031156, + 0.027522865530016288, + 0.097501605587079362, + -0.12976686756709563, + -0.22626469396516913, + 0.3152503517092432, + 0.75113390802157753, + 0.49462389039838539, + 0.11154074335008017 + }; + + double hp1_a[] = { -0.11154074335008017, + 0.49462389039838539, + -0.75113390802157753, + 0.3152503517092432, + 0.22626469396516913, + -0.12976686756709563, + -0.097501605587079362, + 0.027522865530016288, + 0.031582039318031156, + 0.0005538422009938016, + -0.0047772575110106514, + -0.0010773010849955799 + }; + + double lp2_a[] = { 0.11154074335008017, + 0.49462389039838539, + 0.75113390802157753, + 0.3152503517092432, + -0.22626469396516913, + -0.12976686756709563, + 0.097501605587079362, + 0.027522865530016288, + -0.031582039318031156, + 0.0005538422009938016, + 0.0047772575110106514, + -0.0010773010849955799 + }; + + double hp2_a[] = { -0.0010773010849955799, + -0.0047772575110106514, + 0.0005538422009938016, + 0.031582039318031156, + 0.027522865530016288, + -0.097501605587079362, + -0.12976686756709563, + 0.22626469396516913, + 0.3152503517092432, + -0.75113390802157753, + 0.49462389039838539, + -0.11154074335008017 + }; + for (i = 0; i < 12; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 12; + } + else if (!strcmp(name,"db7")){ + double lp1_a[] = { 0.00035371380000103988, + -0.0018016407039998328, + 0.00042957797300470274, + 0.012550998556013784, + -0.01657454163101562, + -0.038029936935034633, + 0.080612609151065898, + 0.071309219267050042, + -0.22403618499416572, + -0.14390600392910627, + 0.4697822874053586, + 0.72913209084655506, + 0.39653931948230575, + 0.077852054085062364 + }; + + double hp1_a[] = { -0.077852054085062364, + 0.39653931948230575, + -0.72913209084655506, + 0.4697822874053586, + 0.14390600392910627, + -0.22403618499416572, + -0.071309219267050042, + 0.080612609151065898, + 0.038029936935034633, + -0.01657454163101562, + -0.012550998556013784, + 0.0004295779730047027, + 0.0018016407039998328, + 0.00035371380000103988 + }; + + double lp2_a[] = { 0.077852054085062364, + 0.39653931948230575, + 0.72913209084655506, + 0.4697822874053586, + -0.14390600392910627, + -0.22403618499416572, + 0.071309219267050042, + 0.080612609151065898, + -0.038029936935034633, + -0.01657454163101562, + 0.012550998556013784, + 0.00042957797300470274, + -0.0018016407039998328, + 0.00035371380000103988 + }; + + double hp2_a[] = { 0.00035371380000103988, + 0.0018016407039998328, + 0.00042957797300470274, + -0.01255099855601378, + -0.01657454163101562, + 0.038029936935034633, + 0.080612609151065898, + -0.071309219267050042, + -0.22403618499416572, + 0.14390600392910627, + 0.4697822874053586, + -0.72913209084655506, + 0.39653931948230575, + -0.077852054085062364 + }; + for (i = 0; i < 14; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 14; + } + else if (!strcmp(name,"db8")){ + double lp1_a[] = { -0.00011747678400228192, + 0.00067544940599855677, + -0.00039174037299597711, + -0.0048703529930106603, + 0.0087460940470156547, + 0.013981027917015516, + -0.044088253931064719, + -0.017369301002022108, + 0.12874742662018601, + 0.00047248457399797254, + -0.28401554296242809, + -0.015829105256023893, + 0.58535468365486909, + 0.67563073629801285, + 0.31287159091446592, + 0.054415842243081609 + }; + + double hp1_a[] = { -0.054415842243081609, + 0.31287159091446592, + -0.67563073629801285, + 0.58535468365486909, + 0.015829105256023893, + -0.28401554296242809, + -0.00047248457399797254, + 0.12874742662018601, + 0.017369301002022108, + -0.044088253931064719, + -0.013981027917015516, + 0.0087460940470156547, + 0.0048703529930106603, + -0.00039174037299597711, + -0.00067544940599855677, + -0.00011747678400228192 + }; + + double lp2_a[] = { 0.054415842243081609, + 0.31287159091446592, + 0.67563073629801285, + 0.58535468365486909, + -0.015829105256023893, + -0.28401554296242809, + 0.00047248457399797254, + 0.12874742662018601, + -0.017369301002022108, + -0.044088253931064719, + 0.013981027917015516, + 0.0087460940470156547, + -0.0048703529930106603, + -0.00039174037299597711, + 0.00067544940599855677, + -0.00011747678400228192 + }; + + double hp2_a[] = { -0.00011747678400228192, + -0.00067544940599855677, + -0.00039174037299597711, + 0.0048703529930106603, + 0.0087460940470156547, + -0.013981027917015516, + -0.044088253931064719, + 0.017369301002022108, + 0.12874742662018601, + -0.00047248457399797254, + -0.28401554296242809, + 0.015829105256023893, + 0.58535468365486909, + -0.67563073629801285, + 0.31287159091446592, + -0.054415842243081609 + }; + for (i = 0; i < 16; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 16; + } + else if (!strcmp(name,"db9")){ + double lp1_a[] = { 3.9347319995026124e-05, + -0.00025196318899817888, + 0.00023038576399541288, + 0.0018476468829611268, + -0.0042815036819047227, + -0.004723204757894831, + 0.022361662123515244, + 0.00025094711499193845, + -0.067632829059523988, + 0.030725681478322865, + 0.14854074933476008, + -0.096840783220879037, + -0.29327378327258685, + 0.13319738582208895, + 0.65728807803663891, + 0.6048231236767786, + 0.24383467463766728, + 0.038077947363167282 + }; + + double hp1_a[] = { -0.038077947363167282, + 0.24383467463766728, + -0.6048231236767786, + 0.65728807803663891, + -0.13319738582208895, + -0.29327378327258685, + 0.096840783220879037, + 0.14854074933476008, + -0.030725681478322865, + -0.067632829059523988, + -0.00025094711499193845, + 0.022361662123515244, + 0.004723204757894831, + -0.0042815036819047227, + -0.0018476468829611268, + 0.00023038576399541288, + 0.00025196318899817888, + 3.9347319995026124e-05 + }; + + double lp2_a[] = { 0.038077947363167282, + 0.24383467463766728, + 0.6048231236767786, + 0.65728807803663891, + 0.13319738582208895, + -0.29327378327258685, + -0.096840783220879037, + 0.14854074933476008, + 0.030725681478322865, + -0.067632829059523988, + 0.00025094711499193845, + 0.022361662123515244, + -0.004723204757894831, + -0.0042815036819047227, + 0.0018476468829611268, + 0.00023038576399541288, + -0.00025196318899817888, + 3.9347319995026124e-05 + }; + + double hp2_a[] = { 3.9347319995026124e-05, + 0.00025196318899817888, + 0.00023038576399541288, + -0.0018476468829611268, + -0.0042815036819047227, + 0.004723204757894831, + 0.022361662123515244, + -0.00025094711499193845, + -0.067632829059523988, + -0.030725681478322865, + 0.14854074933476008, + 0.096840783220879037, + -0.29327378327258685, + -0.13319738582208895, + 0.65728807803663891, + -0.6048231236767786, + 0.24383467463766728, + -0.038077947363167282 + }; + for (i = 0; i < 18; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 18; + } + + else if (!strcmp(name,"db10")){ + double lp1_a[] = { -1.3264203002354869e-05, + 9.3588670001089845e-05, + -0.0001164668549943862, + -0.00068585669500468248, + 0.0019924052949908499, + 0.0013953517469940798, + -0.010733175482979604, + 0.0036065535669883944, + 0.033212674058933238, + -0.029457536821945671, + -0.071394147165860775, + 0.093057364603806592, + 0.12736934033574265, + -0.19594627437659665, + -0.24984642432648865, + 0.28117234366042648, + 0.68845903945259213, + 0.52720118893091983, + 0.18817680007762133, + 0.026670057900950818 + }; + + double hp1_a[] = { -0.026670057900950818, + 0.18817680007762133, + -0.52720118893091983, + 0.68845903945259213, + -0.28117234366042648, + -0.24984642432648865, + 0.19594627437659665, + 0.12736934033574265, + -0.093057364603806592, + -0.071394147165860775, + 0.029457536821945671, + 0.033212674058933238, + -0.0036065535669883944, + -0.010733175482979604, + -0.0013953517469940798, + 0.0019924052949908499, + 0.00068585669500468248, + -0.0001164668549943862, + -9.3588670001089845e-05, + -1.3264203002354869e-05 + }; + + double lp2_a[] = { 0.026670057900950818, + 0.18817680007762133, + 0.52720118893091983, + 0.68845903945259213, + 0.28117234366042648, + -0.24984642432648865, + -0.19594627437659665, + 0.12736934033574265, + 0.093057364603806592, + -0.071394147165860775, + -0.029457536821945671, + 0.033212674058933238, + 0.0036065535669883944, + -0.010733175482979604, + 0.0013953517469940798, + 0.0019924052949908499, + -0.00068585669500468248, + -0.0001164668549943862, + 9.3588670001089845e-05, + -1.3264203002354869e-05 + }; + + double hp2_a[] = { -1.3264203002354869e-05, + -9.3588670001089845e-05, + -0.0001164668549943862, + 0.00068585669500468248, + 0.0019924052949908499, + -0.0013953517469940798, + -0.010733175482979604, + -0.0036065535669883944, + 0.033212674058933238, + 0.029457536821945671, + -0.071394147165860775, + -0.093057364603806592, + 0.12736934033574265, + 0.19594627437659665, + -0.24984642432648865, + -0.28117234366042648, + 0.68845903945259213, + -0.52720118893091983, + 0.18817680007762133, + -0.026670057900950818 + }; + for (i = 0; i < 20; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 20; + } + + else if (!strcmp(name,"db12")){ + double lp1_a[] = { -1.5290717580684923e-06, + 1.2776952219379579e-05, + -2.4241545757030318e-05, + -8.8504109208203182e-05, + 0.00038865306282092672, + 6.5451282125215034e-06, + -0.0021795036186277044, + 0.0022486072409952287, + 0.0067114990087955486, + -0.012840825198299882, + -0.01221864906974642, + 0.041546277495087637, + 0.010849130255828966, + -0.09643212009649671, + 0.0053595696743599965, + 0.18247860592758275, + -0.023779257256064865, + -0.31617845375277914, + -0.044763885653777619, + 0.51588647842780067, + 0.65719872257929113, + 0.37735513521420411, + 0.10956627282118277, + 0.013112257957229239 + }; + + double hp1_a[] = { -0.013112257957229239, + 0.10956627282118277, + -0.37735513521420411, + 0.65719872257929113, + -0.51588647842780067, + -0.044763885653777619, + 0.31617845375277914, + -0.023779257256064865, + -0.18247860592758275, + 0.0053595696743599965, + 0.09643212009649671, + 0.010849130255828966, + -0.041546277495087637, + -0.01221864906974642, + 0.012840825198299882, + 0.0067114990087955486, + -0.0022486072409952287, + -0.0021795036186277044, + -6.5451282125215034e-06, + 0.00038865306282092672, + 8.8504109208203182e-05, + -2.4241545757030318e-05, + -1.2776952219379579e-05, + -1.5290717580684923e-06 + }; + + double lp2_a[] = { 0.013112257957229239, + 0.10956627282118277, + 0.37735513521420411, + 0.65719872257929113, + 0.51588647842780067, + -0.044763885653777619, + -0.31617845375277914, + -0.023779257256064865, + 0.18247860592758275, + 0.0053595696743599965, + -0.09643212009649671, + 0.010849130255828966, + 0.041546277495087637, + -0.01221864906974642, + -0.012840825198299882, + 0.0067114990087955486, + 0.0022486072409952287, + -0.0021795036186277044, + 6.5451282125215034e-06, + 0.00038865306282092672, + -8.8504109208203182e-05, + -2.4241545757030318e-05, + 1.2776952219379579e-05, + -1.5290717580684923e-06 + }; + + double hp2_a[] = { -1.5290717580684923e-06, + -1.2776952219379579e-05, + -2.4241545757030318e-05, + 8.8504109208203182e-05, + 0.00038865306282092672, + -6.5451282125215034e-06, + -0.0021795036186277044, + -0.0022486072409952287, + 0.0067114990087955486, + 0.012840825198299882, + -0.01221864906974642, + -0.041546277495087637, + 0.010849130255828966, + 0.09643212009649671, + 0.0053595696743599965, + -0.18247860592758275, + -0.023779257256064865, + 0.31617845375277914, + -0.044763885653777619, + -0.51588647842780067, + 0.65719872257929113, + -0.37735513521420411, + 0.10956627282118277, + -0.013112257957229239 + }; + for (i = 0; i < 24; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 24; + } + else if (!strcmp(name,"db13")){ + double lp1_a[] = { 5.2200350984547998e-07, + -4.7004164793608082e-06, + 1.0441930571407941e-05, + 3.0678537579324358e-05, + -0.00016512898855650571, + 4.9251525126285676e-05, + 0.00093232613086724904, + -0.0013156739118922766, + -0.002761911234656831, + 0.0072555894016171187, + 0.0039239414487955773, + -0.023831420710327809, + 0.0023799722540522269, + 0.056139477100276156, + -0.026488406475345658, + -0.10580761818792761, + 0.072948933656788742, + 0.17947607942935084, + -0.12457673075080665, + -0.31497290771138414, + 0.086985726179645007, + 0.58888957043121193, + 0.61105585115878114, + 0.31199632216043488, + 0.082861243872901946, + 0.0092021335389622788 + }; + + double hp1_a[] = { -0.0092021335389622788, + 0.082861243872901946, + -0.31199632216043488, + 0.61105585115878114, + -0.58888957043121193, + 0.086985726179645007, + 0.31497290771138414, + -0.12457673075080665, + -0.17947607942935084, + 0.072948933656788742, + 0.10580761818792761, + -0.026488406475345658, + -0.056139477100276156, + 0.0023799722540522269, + 0.023831420710327809, + 0.0039239414487955773, + -0.0072555894016171187, + -0.002761911234656831, + 0.0013156739118922766, + 0.00093232613086724904, + -4.9251525126285676e-05, + -0.00016512898855650571, + -3.0678537579324358e-05, + 1.0441930571407941e-05, + 4.7004164793608082e-06, + 5.2200350984547998e-07 + }; + + double lp2_a[] = { 0.0092021335389622788, + 0.082861243872901946, + 0.31199632216043488, + 0.61105585115878114, + 0.58888957043121193, + 0.086985726179645007, + -0.31497290771138414, + -0.12457673075080665, + 0.17947607942935084, + 0.072948933656788742, + -0.10580761818792761, + -0.026488406475345658, + 0.056139477100276156, + 0.0023799722540522269, + -0.023831420710327809, + 0.0039239414487955773, + 0.0072555894016171187, + -0.002761911234656831, + -0.0013156739118922766, + 0.00093232613086724904, + 4.9251525126285676e-05, + -0.00016512898855650571, + 3.0678537579324358e-05, + 1.0441930571407941e-05, + -4.7004164793608082e-06, + 5.2200350984547998e-07 + }; + + double hp2_a[] = { 5.2200350984547998e-07, + 4.7004164793608082e-06, + 1.0441930571407941e-05, + -3.0678537579324358e-05, + -0.00016512898855650571, + -4.9251525126285676e-05, + 0.00093232613086724904, + 0.0013156739118922766, + -0.002761911234656831, + -0.0072555894016171187, + 0.0039239414487955773, + 0.023831420710327809, + 0.0023799722540522269, + -0.056139477100276156, + -0.026488406475345658, + 0.10580761818792761, + 0.072948933656788742, + -0.17947607942935084, + -0.12457673075080665, + 0.31497290771138414, + 0.086985726179645007, + -0.58888957043121193, + 0.61105585115878114, + -0.31199632216043488, + 0.082861243872901946, + -0.0092021335389622788 + }; + for (i = 0; i < 26; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 26; + } + + else if (!strcmp(name,"db11")){ + double lp1_a[] = { 4.4942742772363519e-06, + -3.4634984186983789e-05, + 5.4439074699366381e-05, + 0.00024915252355281426, + -0.00089302325066623663, + -0.00030859285881515924, + 0.0049284176560587777, + -0.0033408588730145018, + -0.015364820906201324, + 0.020840904360180039, + 0.031335090219045313, + -0.066438785695020222, + -0.04647995511667613, + 0.14981201246638268, + 0.066043588196690886, + -0.27423084681792875, + -0.16227524502747828, + 0.41196436894789695, + 0.68568677491617847, + 0.44989976435603013, + 0.14406702115061959, + 0.018694297761470441 + }; + + double hp1_a[] = { -0.018694297761470441, + 0.14406702115061959, + -0.44989976435603013, + 0.68568677491617847, + -0.41196436894789695, + -0.16227524502747828, + 0.27423084681792875, + 0.066043588196690886, + -0.14981201246638268, + -0.04647995511667613, + 0.066438785695020222, + 0.031335090219045313, + -0.020840904360180039, + -0.015364820906201324, + 0.0033408588730145018, + 0.0049284176560587777, + 0.00030859285881515924, + -0.00089302325066623663, + -0.00024915252355281426, + 5.4439074699366381e-05, + 3.4634984186983789e-05, + 4.4942742772363519e-06 + }; + + double lp2_a[] = { 0.018694297761470441, + 0.14406702115061959, + 0.44989976435603013, + 0.68568677491617847, + 0.41196436894789695, + -0.16227524502747828, + -0.27423084681792875, + 0.066043588196690886, + 0.14981201246638268, + -0.04647995511667613, + -0.066438785695020222, + 0.031335090219045313, + 0.020840904360180039, + -0.015364820906201324, + -0.0033408588730145018, + 0.0049284176560587777, + -0.00030859285881515924, + -0.00089302325066623663, + 0.00024915252355281426, + 5.4439074699366381e-05, + -3.4634984186983789e-05, + 4.4942742772363519e-06 + }; + + double hp2_a[] = { 4.4942742772363519e-06, + 3.4634984186983789e-05, + 5.4439074699366381e-05, + -0.00024915252355281426, + -0.00089302325066623663, + 0.00030859285881515924, + 0.0049284176560587777, + 0.0033408588730145018, + -0.015364820906201324, + -0.020840904360180039, + 0.031335090219045313, + 0.066438785695020222, + -0.04647995511667613, + -0.14981201246638268, + 0.066043588196690886, + 0.27423084681792875, + -0.16227524502747828, + -0.41196436894789695, + 0.68568677491617847, + -0.44989976435603013, + 0.14406702115061959, + -0.018694297761470441 + }; + for (i = 0; i < 22; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 22; + } + + else if (!strcmp(name,"db14")){ + double lp1_a[] = { -1.7871399683109222e-07, + 1.7249946753674012e-06, + -4.3897049017804176e-06, + -1.0337209184568496e-05, + 6.875504252695734e-05, + -4.1777245770370672e-05, + -0.00038683194731287514, + 0.00070802115423540481, + 0.001061691085606874, + -0.003849638868019787, + -0.00074621898926387534, + 0.012789493266340071, + -0.0056150495303375755, + -0.030185351540353976, + 0.026981408307947971, + 0.05523712625925082, + -0.071548955503983505, + -0.086748411568110598, + 0.13998901658445695, + 0.13839521386479153, + -0.21803352999321651, + -0.27168855227867705, + 0.21867068775886594, + 0.63118784910471981, + 0.55430561794077093, + 0.25485026779256437, + 0.062364758849384874, + 0.0064611534600864905 + }; + + double hp1_a[] = { -0.0064611534600864905, + 0.062364758849384874, + -0.25485026779256437, + 0.55430561794077093, + -0.63118784910471981, + 0.21867068775886594, + 0.27168855227867705, + -0.21803352999321651, + -0.13839521386479153, + 0.13998901658445695, + 0.086748411568110598, + -0.071548955503983505, + -0.05523712625925082, + 0.026981408307947971, + 0.030185351540353976, + -0.0056150495303375755, + -0.012789493266340071, + -0.00074621898926387534, + 0.003849638868019787, + 0.001061691085606874, + -0.00070802115423540481, + -0.00038683194731287514, + 4.1777245770370672e-05, + 6.875504252695734e-05, + 1.0337209184568496e-05, + -4.3897049017804176e-06, + -1.7249946753674012e-06, + -1.7871399683109222e-07 + }; + + double lp2_a[] = { 0.0064611534600864905, + 0.062364758849384874, + 0.25485026779256437, + 0.55430561794077093, + 0.63118784910471981, + 0.21867068775886594, + -0.27168855227867705, + -0.21803352999321651, + 0.13839521386479153, + 0.13998901658445695, + -0.086748411568110598, + -0.071548955503983505, + 0.05523712625925082, + 0.026981408307947971, + -0.030185351540353976, + -0.0056150495303375755, + 0.012789493266340071, + -0.00074621898926387534, + -0.003849638868019787, + 0.001061691085606874, + 0.00070802115423540481, + -0.00038683194731287514, + -4.1777245770370672e-05, + 6.875504252695734e-05, + -1.0337209184568496e-05, + -4.3897049017804176e-06, + 1.7249946753674012e-06, + -1.7871399683109222e-07 + }; + + double hp2_a[] = { -1.7871399683109222e-07, + -1.7249946753674012e-06, + -4.3897049017804176e-06, + 1.0337209184568496e-05, + 6.875504252695734e-05, + 4.1777245770370672e-05, + -0.00038683194731287514, + -0.00070802115423540481, + 0.001061691085606874, + 0.003849638868019787, + -0.00074621898926387534, + -0.012789493266340071, + -0.0056150495303375755, + 0.030185351540353976, + 0.026981408307947971, + -0.05523712625925082, + -0.071548955503983505, + 0.086748411568110598, + 0.13998901658445695, + -0.13839521386479153, + -0.21803352999321651, + 0.27168855227867705, + 0.21867068775886594, + -0.63118784910471981, + 0.55430561794077093, + -0.25485026779256437, + 0.062364758849384874, + -0.0064611534600864905 + }; + for (i = 0; i < 28; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 28; + } + else if (!strcmp(name,"db15")){ + double lp1_a[] = { 6.1333599133037138e-08, + -6.3168823258794506e-07, + 1.8112704079399406e-06, + 3.3629871817363823e-06, + -2.8133296266037558e-05, + 2.579269915531323e-05, + 0.00015589648992055726, + -0.00035956524436229364, + -0.00037348235413726472, + 0.0019433239803823459, + -0.00024175649075894543, + -0.0064877345603061454, + 0.0051010003604228726, + 0.015083918027862582, + -0.020810050169636805, + -0.025767007328366939, + 0.054780550584559995, + 0.033877143923563204, + -0.11112093603713753, + -0.039666176555733602, + 0.19014671400708816, + 0.065282952848765688, + -0.28888259656686216, + -0.19320413960907623, + 0.33900253545462167, + 0.64581314035721027, + 0.49263177170797529, + 0.20602386398692688, + 0.046743394892750617, + 0.0045385373615773762 + }; + + double hp1_a[] = { -0.0045385373615773762, + 0.046743394892750617, + -0.20602386398692688, + 0.49263177170797529, + -0.64581314035721027, + 0.33900253545462167, + 0.19320413960907623, + -0.28888259656686216, + -0.065282952848765688, + 0.19014671400708816, + 0.039666176555733602, + -0.11112093603713753, + -0.033877143923563204, + 0.054780550584559995, + 0.025767007328366939, + -0.020810050169636805, + -0.015083918027862582, + 0.0051010003604228726, + 0.0064877345603061454, + -0.00024175649075894543, + -0.0019433239803823459, + -0.00037348235413726472, + 0.00035956524436229364, + 0.00015589648992055726, + -2.579269915531323e-05, + -2.8133296266037558e-05, + -3.3629871817363823e-06, + 1.8112704079399406e-06, + 6.3168823258794506e-07, + 6.1333599133037138e-08 + }; + + double lp2_a[] = { 0.0045385373615773762, + 0.046743394892750617, + 0.20602386398692688, + 0.49263177170797529, + 0.64581314035721027, + 0.33900253545462167, + -0.19320413960907623, + -0.28888259656686216, + 0.065282952848765688, + 0.19014671400708816, + -0.039666176555733602, + -0.11112093603713753, + 0.033877143923563204, + 0.054780550584559995, + -0.025767007328366939, + -0.020810050169636805, + 0.015083918027862582, + 0.0051010003604228726, + -0.0064877345603061454, + -0.00024175649075894543, + 0.0019433239803823459, + -0.00037348235413726472, + -0.00035956524436229364, + 0.00015589648992055726, + 2.579269915531323e-05, + -2.8133296266037558e-05, + 3.3629871817363823e-06, + 1.8112704079399406e-06, + -6.3168823258794506e-07, + 6.1333599133037138e-08 + }; + + double hp2_a[] = { 6.1333599133037138e-08, + 6.3168823258794506e-07, + 1.8112704079399406e-06, + -3.3629871817363823e-06, + -2.8133296266037558e-05, + -2.579269915531323e-05, + 0.00015589648992055726, + 0.00035956524436229364, + -0.00037348235413726472, + -0.0019433239803823459, + -0.00024175649075894543, + 0.0064877345603061454, + 0.0051010003604228726, + -0.015083918027862582, + -0.020810050169636805, + 0.025767007328366939, + 0.054780550584559995, + -0.033877143923563204, + -0.11112093603713753, + 0.039666176555733602, + 0.19014671400708816, + -0.065282952848765688, + -0.28888259656686216, + 0.19320413960907623, + 0.33900253545462167, + -0.64581314035721027, + 0.49263177170797529, + -0.20602386398692688, + 0.046743394892750617, + -0.0045385373615773762 + }; + for (i = 0; i < 30; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 30; + } + else if (!strcmp(name,"bior1.1")){ + double lp1_a[] = { 0.70710678118654757, + 0.70710678118654757 + }; + + double hp1_a[] = { -0.70710678118654757, + 0.70710678118654757 + }; + + double lp2_a[] = { 0.70710678118654757, + 0.70710678118654757 + }; + + double hp2_a[] = { 0.70710678118654757, + -0.70710678118654757 + }; + for (i = 0; i < 2; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 2; + } + + else if (!strcmp(name,"bior1.3")){ + double lp1_a[] = { -0.088388347648318447, + 0.088388347648318447, + 0.70710678118654757, + 0.70710678118654757, + 0.088388347648318447, + -0.088388347648318447, + }; + + double hp1_a[] = { 0.0, + 0.0, + -0.70710678118654757, + 0.70710678118654757, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.70710678118654757, + 0.70710678118654757, + 0.0, + 0.0 + }; + + double hp2_a[] = { -0.088388347648318447, + -0.088388347648318447, + 0.70710678118654757, + -0.70710678118654757, + 0.088388347648318447, + 0.088388347648318447 + }; + for (i = 0; i < 6; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 6; + } + + else if (!strcmp(name,"bior1.5")){ + double lp1_a[] = { 0.01657281518405971, + -0.01657281518405971, + -0.12153397801643787, + 0.12153397801643787, + 0.70710678118654757, + 0.70710678118654757, + 0.12153397801643787, + -0.12153397801643787, + -0.01657281518405971, + 0.01657281518405971 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + -0.70710678118654757, + 0.70710678118654757, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.70710678118654757, + 0.70710678118654757, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.01657281518405971, + 0.01657281518405971, + -0.12153397801643787, + -0.12153397801643787, + 0.70710678118654757, + -0.70710678118654757, + 0.12153397801643787, + 0.12153397801643787, + -0.01657281518405971, + -0.01657281518405971 + }; + for (i = 0; i < 10; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 10; + } + + else if (!strcmp(name,"bior2.2")){ + double lp1_a[] = { 0.0, + -0.17677669529663689, + 0.35355339059327379, + 1.0606601717798214, + 0.35355339059327379, + -0.17677669529663689 + }; + + double hp1_a[] = { 0.0, + 0.35355339059327379, + -0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.35355339059327379, + 0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.0, + 0.17677669529663689, + 0.35355339059327379, + -1.0606601717798214, + 0.35355339059327379, + 0.17677669529663689 + + }; + for (i = 0; i < 6; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 6; + } + + else if (!strcmp(name,"bior2.4")){ + double lp1_a[] = { 0.0, + 0.033145630368119419, + -0.066291260736238838, + -0.17677669529663689, + 0.4198446513295126, + 0.99436891104358249, + 0.4198446513295126, + -0.17677669529663689, + -0.066291260736238838, + 0.033145630368119419 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.35355339059327379, + -0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0, + 0.0, + 0.0 + + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.35355339059327379, + 0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0, + 0.0, + 0.0 + + }; + + double hp2_a[] = { 0.0, + -0.033145630368119419, + -0.066291260736238838, + 0.17677669529663689, + 0.4198446513295126, + -0.99436891104358249, + 0.4198446513295126, + 0.17677669529663689, + -0.066291260736238838, + -0.033145630368119419 + }; + for (i = 0; i < 10; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 10; + } + + else if (!strcmp(name,"bior2.6")){ + double lp1_a[] = { 0.0, + -0.0069053396600248784, + 0.013810679320049757, + 0.046956309688169176, + -0.10772329869638811, + -0.16987135563661201, + 0.44746600996961211, + 0.96674755240348298, + 0.44746600996961211, + -0.16987135563661201, + -0.10772329869638811, + 0.046956309688169176, + 0.013810679320049757, + -0.0069053396600248784 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.35355339059327379, + -0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.35355339059327379, + 0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.0, + 0.0069053396600248784, + 0.013810679320049757, + -0.046956309688169176, + -0.10772329869638811, + 0.16987135563661201, + 0.44746600996961211, + -0.96674755240348298, + 0.44746600996961211, + 0.16987135563661201, + -0.10772329869638811, + -0.046956309688169176, + 0.013810679320049757, + 0.0069053396600248784 + }; + for (i = 0; i < 14; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 14; + } + else if (!strcmp(name,"bior2.8")){ + double lp1_a[] = { 0.0, + 0.0015105430506304422, + -0.0030210861012608843, + -0.012947511862546647, + 0.028916109826354178, + 0.052998481890690945, + -0.13491307360773608, + -0.16382918343409025, + 0.46257144047591658, + 0.95164212189717856, + 0.46257144047591658, + -0.16382918343409025, + -0.13491307360773608, + 0.052998481890690945, + 0.028916109826354178, + -0.012947511862546647, + -0.0030210861012608843, + 0.0015105430506304422 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.35355339059327379, + -0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.35355339059327379, + 0.70710678118654757, + 0.35355339059327379, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.0, + -0.0015105430506304422, + -0.0030210861012608843, + 0.012947511862546647, + 0.028916109826354178, + -0.052998481890690945, + -0.13491307360773608, + 0.16382918343409025, + 0.46257144047591658, + -0.95164212189717856, + 0.46257144047591658, + 0.16382918343409025, + -0.13491307360773608, + -0.052998481890690945, + 0.028916109826354178, + 0.012947511862546647, + -0.0030210861012608843, + -0.0015105430506304422 + }; + for (i = 0; i < 18; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 18; + } + + else if (!strcmp(name,"bior3.1")){ + double lp1_a[] = { -0.35355339059327379, + 1.0606601717798214, + 1.0606601717798214, + -0.35355339059327379 + }; + + double hp1_a[] = { -0.17677669529663689, + 0.53033008588991071, + -0.53033008588991071, + 0.17677669529663689 + }; + + double lp2_a[] = { 0.17677669529663689, + 0.53033008588991071, + 0.53033008588991071, + 0.17677669529663689 + }; + + double hp2_a[] = { -0.35355339059327379, + -1.0606601717798214, + 1.0606601717798214, + 0.35355339059327379 + }; + for (i = 0; i < 4; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 4; + } + else if (!strcmp(name,"bior3.3")){ + double lp1_a[] = { 0.066291260736238838, + -0.19887378220871652, + -0.15467960838455727, + 0.99436891104358249, + 0.99436891104358249, + -0.15467960838455727, + -0.19887378220871652, + 0.066291260736238838 + }; + + double hp1_a[] = { 0.0, + 0.0, + -0.17677669529663689, + 0.53033008588991071, + -0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.17677669529663689, + 0.53033008588991071, + 0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.066291260736238838, + 0.19887378220871652, + -0.15467960838455727, + -0.99436891104358249, + 0.99436891104358249, + 0.15467960838455727, + -0.19887378220871652, + -0.066291260736238838 + }; + for (i = 0; i < 8; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 8; + } + else if (!strcmp(name,"bior3.5")){ + double lp1_a[] = { -0.013810679320049757, + 0.041432037960149271, + 0.052480581416189075, + -0.26792717880896527, + -0.071815532464258744, + 0.96674755240348298, + 0.96674755240348298, + -0.071815532464258744, + -0.26792717880896527, + 0.052480581416189075, + 0.041432037960149271, + -0.013810679320049757 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + -0.17677669529663689, + 0.53033008588991071, + -0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.17677669529663689, + 0.53033008588991071, + 0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double hp2_a[] = { -0.013810679320049757, + -0.041432037960149271, + 0.052480581416189075, + 0.26792717880896527, + -0.071815532464258744, + -0.96674755240348298, + 0.96674755240348298, + 0.071815532464258744, + -0.26792717880896527, + -0.052480581416189075, + 0.041432037960149271, + 0.013810679320049757 + }; + for (i = 0; i < 12; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 12; + } + + else if (!strcmp(name,"bior3.7")){ + double lp1_a[] = { 0.0030210861012608843, + -0.0090632583037826529, + -0.016831765421310641, + 0.074663985074019001, + 0.031332978707362888, + -0.301159125922835, + -0.026499240945345472, + 0.95164212189717856, + 0.95164212189717856, + -0.026499240945345472, + -0.301159125922835, + 0.031332978707362888, + 0.074663985074019001, + -0.016831765421310641, + -0.0090632583037826529, + 0.0030210861012608843 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + -0.17677669529663689, + 0.53033008588991071, + -0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.17677669529663689, + 0.53033008588991071, + 0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.0030210861012608843, + 0.0090632583037826529, + -0.016831765421310641, + -0.074663985074019001, + 0.031332978707362888, + 0.301159125922835, + -0.026499240945345472, + -0.95164212189717856, + 0.95164212189717856, + 0.026499240945345472, + -0.301159125922835, + -0.031332978707362888, + 0.074663985074019001, + 0.016831765421310641, + -0.0090632583037826529, + -0.0030210861012608843 + }; + for (i = 0; i < 16; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 16; + } + else if (!strcmp(name,"bior3.9")){ + double lp1_a[] = { -0.00067974437278369901, + 0.0020392331183510968, + 0.0050603192196119811, + -0.020618912641105536, + -0.014112787930175846, + 0.09913478249423216, + 0.012300136269419315, + -0.32019196836077857, + 0.0020500227115698858, + 0.94212570067820678, + 0.94212570067820678, + 0.0020500227115698858, + -0.32019196836077857, + 0.012300136269419315, + 0.09913478249423216, + -0.014112787930175846, + -0.020618912641105536, + 0.0050603192196119811, + 0.0020392331183510968, + -0.00067974437278369901 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + -0.17677669529663689, + 0.53033008588991071, + -0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.17677669529663689, + 0.53033008588991071, + 0.53033008588991071, + 0.17677669529663689, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double hp2_a[] = { -0.00067974437278369901, + -0.0020392331183510968, + 0.0050603192196119811, + 0.020618912641105536, + -0.014112787930175846, + -0.09913478249423216, + 0.012300136269419315, + 0.32019196836077857, + 0.0020500227115698858, + -0.94212570067820678, + 0.94212570067820678, + -0.0020500227115698858, + -0.32019196836077857, + -0.012300136269419315, + 0.09913478249423216, + 0.014112787930175846, + -0.020618912641105536, + -0.0050603192196119811, + 0.0020392331183510968, + 0.00067974437278369901 + }; + for (i = 0; i < 20; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 20; + } + else if (!strcmp(name,"bior4.4")){ + double lp1_a[] = { 0.0, + 0.03782845550726404, + -0.023849465019556843, + -0.11062440441843718, + 0.37740285561283066, + 0.85269867900889385, + 0.37740285561283066, + -0.11062440441843718, + -0.023849465019556843, + 0.03782845550726404 + }; + + double hp1_a[] = { 0.0, + -0.064538882628697058, + 0.040689417609164058, + 0.41809227322161724, + -0.7884856164055829, + 0.41809227322161724, + 0.040689417609164058, + -0.064538882628697058, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + -0.064538882628697058, + -0.040689417609164058, + 0.41809227322161724, + 0.7884856164055829, + 0.41809227322161724, + -0.040689417609164058, + -0.064538882628697058, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.0, + -0.03782845550726404, + -0.023849465019556843, + 0.11062440441843718, + 0.37740285561283066, + -0.85269867900889385, + 0.37740285561283066, + 0.11062440441843718, + -0.023849465019556843, + -0.03782845550726404 + }; + for (i = 0; i < 10; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 10; + } + else if (!strcmp(name,"bior5.5")){ + double lp1_a[] = { 0.0, + 0.0, + 0.03968708834740544, + 0.0079481086372403219, + -0.054463788468236907, + 0.34560528195603346, + 0.73666018142821055, + 0.34560528195603346, + -0.054463788468236907, + 0.0079481086372403219, + 0.03968708834740544, + 0.0 + }; + + double hp1_a[] = { -0.013456709459118716, + -0.0026949668801115071, + 0.13670658466432914, + -0.093504697400938863, + -0.47680326579848425, + 0.89950610974864842, + -0.47680326579848425, + -0.093504697400938863, + 0.13670658466432914, + -0.0026949668801115071, + -0.013456709459118716, + 0.0 + }; + + double lp2_a[] = { 0.013456709459118716, + -0.0026949668801115071, + -0.13670658466432914, + -0.093504697400938863, + 0.47680326579848425, + 0.89950610974864842, + 0.47680326579848425, + -0.093504697400938863, + -0.13670658466432914, + -0.0026949668801115071, + 0.013456709459118716, + 0.0 + }; + + double hp2_a[] = { 0.0, + 0.0, + 0.03968708834740544, + -0.0079481086372403219, + -0.054463788468236907, + -0.34560528195603346, + 0.73666018142821055, + -0.34560528195603346, + -0.054463788468236907, + -0.0079481086372403219, + 0.03968708834740544, + 0.0 + }; + for (i = 0; i < 12; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 12; + } + else if (!strcmp(name,"bior6.8")){ + double lp1_a[] = { 0.0, + 0.0019088317364812906, + -0.0019142861290887667, + -0.016990639867602342, + 0.01193456527972926, + 0.04973290349094079, + -0.077263173167204144, + -0.09405920349573646, + 0.42079628460982682, + 0.82592299745840225, + 0.42079628460982682, + -0.09405920349573646, + -0.077263173167204144, + 0.04973290349094079, + 0.01193456527972926, + -0.016990639867602342, + -0.0019142861290887667, + 0.0019088317364812906 + }; + + double hp1_a[] = { 0.0, + 0.0, + 0.0, + 0.014426282505624435, + -0.014467504896790148, + -0.078722001062628819, + 0.040367979030339923, + 0.41784910915027457, + -0.75890772945365415, + 0.41784910915027457, + 0.040367979030339923, + -0.078722001062628819, + -0.014467504896790148, + 0.014426282505624435, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double lp2_a[] = { 0.0, + 0.0, + 0.0, + 0.014426282505624435, + 0.014467504896790148, + -0.078722001062628819, + -0.040367979030339923, + 0.41784910915027457, + 0.75890772945365415, + 0.41784910915027457, + -0.040367979030339923, + -0.078722001062628819, + 0.014467504896790148, + 0.014426282505624435, + 0.0, + 0.0, + 0.0, + 0.0 + }; + + double hp2_a[] = { 0.0, + -0.0019088317364812906, + -0.0019142861290887667, + 0.016990639867602342, + 0.01193456527972926, + -0.04973290349094079, + -0.077263173167204144, + 0.09405920349573646, + 0.42079628460982682, + -0.82592299745840225, + 0.42079628460982682, + 0.09405920349573646, + -0.077263173167204144, + -0.04973290349094079, + 0.01193456527972926, + 0.016990639867602342, + -0.0019142861290887667, + -0.0019088317364812906 + }; + for (i = 0; i < 18; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 18; + } + + else if (!strcmp(name,"coif1")){ + double lp1_a[] = { -0.01565572813546454, + -0.072732619512853897, + 0.38486484686420286, + 0.85257202021225542, + 0.33789766245780922, + -0.072732619512853897 + }; + + double hp1_a[] = { 0.072732619512853897, + 0.33789766245780922, + -0.85257202021225542, + 0.38486484686420286, + 0.072732619512853897, + -0.01565572813546454 + }; + + double lp2_a[] = { -0.072732619512853897, + 0.33789766245780922, + 0.85257202021225542, + 0.38486484686420286, + -0.072732619512853897, + -0.01565572813546454 + }; + + double hp2_a[] = { -0.01565572813546454, + 0.072732619512853897, + 0.38486484686420286, + -0.85257202021225542, + 0.33789766245780922, + 0.072732619512853897 + }; + for (i = 0; i < 6; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 6; + } + else if (!strcmp(name,"coif2")){ + double lp1_a[] = { -0.00072054944536451221, + -0.0018232088707029932, + 0.0056114348193944995, + 0.023680171946334084, + -0.059434418646456898, + -0.076488599078306393, + 0.41700518442169254, + 0.81272363544554227, + 0.38611006682116222, + -0.067372554721963018, + -0.041464936781759151, + 0.016387336463522112 + }; + + double hp1_a[] = { -0.016387336463522112, + -0.041464936781759151, + 0.067372554721963018, + 0.38611006682116222, + -0.81272363544554227, + 0.41700518442169254, + 0.076488599078306393, + -0.059434418646456898, + -0.023680171946334084, + 0.0056114348193944995, + 0.0018232088707029932, + -0.00072054944536451221 + }; + + double lp2_a[] = { 0.016387336463522112, + -0.041464936781759151, + -0.067372554721963018, + 0.38611006682116222, + 0.81272363544554227, + 0.41700518442169254, + -0.076488599078306393, + -0.059434418646456898, + 0.023680171946334084, + 0.0056114348193944995, + -0.0018232088707029932, + -0.00072054944536451221 + }; + + double hp2_a[] = { -0.00072054944536451221, + 0.0018232088707029932, + 0.0056114348193944995, + -0.023680171946334084, + -0.059434418646456898, + 0.076488599078306393, + 0.41700518442169254, + -0.81272363544554227, + 0.38611006682116222, + 0.067372554721963018, + -0.041464936781759151, + -0.016387336463522112 + }; + for (i = 0; i < 12; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 12; + } + else if (!strcmp(name,"coif3")){ + double lp1_a[] = { -3.4599772836212559e-05, + -7.0983303138141252e-05, + 0.00046621696011288631, + 0.0011175187708906016, + -0.0025745176887502236, + -0.0090079761366615805, + 0.015880544863615904, + 0.034555027573061628, + -0.082301927106885983, + -0.071799821619312018, + 0.42848347637761874, + 0.79377722262562056, + 0.4051769024096169, + -0.061123390002672869, + -0.0657719112818555, + 0.023452696141836267, + 0.0077825964273254182, + -0.0037935128644910141 + }; + + double hp1_a[] = { 0.0037935128644910141, + 0.0077825964273254182, + -0.023452696141836267, + -0.0657719112818555, + 0.061123390002672869, + 0.4051769024096169, + -0.79377722262562056, + 0.42848347637761874, + 0.071799821619312018, + -0.082301927106885983, + -0.034555027573061628, + 0.015880544863615904, + 0.0090079761366615805, + -0.0025745176887502236, + -0.0011175187708906016, + 0.00046621696011288631, + 7.0983303138141252e-05, + -3.4599772836212559e-05 + }; + + double lp2_a[] = { -0.0037935128644910141, + 0.0077825964273254182, + 0.023452696141836267, + -0.0657719112818555, + -0.061123390002672869, + 0.4051769024096169, + 0.79377722262562056, + 0.42848347637761874, + -0.071799821619312018, + -0.082301927106885983, + 0.034555027573061628, + 0.015880544863615904, + -0.0090079761366615805, + -0.0025745176887502236, + 0.0011175187708906016, + 0.00046621696011288631, + -7.0983303138141252e-05, + -3.4599772836212559e-05 + }; + + double hp2_a[] = { -3.4599772836212559e-05, + 7.0983303138141252e-05, + 0.00046621696011288631, + -0.0011175187708906016, + -0.0025745176887502236, + 0.0090079761366615805, + 0.015880544863615904, + -0.034555027573061628, + -0.082301927106885983, + 0.071799821619312018, + 0.42848347637761874, + -0.79377722262562056, + 0.4051769024096169, + 0.061123390002672869, + -0.0657719112818555, + -0.023452696141836267, + 0.0077825964273254182, + 0.0037935128644910141 + }; + for (i = 0; i < 18; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 18; + } + else if (!strcmp(name,"coif4")){ + double lp1_a[] = { -1.7849850030882614e-06, + -3.2596802368833675e-06, + 3.1229875865345646e-05, + 6.2339034461007128e-05, + -0.00025997455248771324, + -0.00058902075624433831, + 0.0012665619292989445, + 0.0037514361572784571, + -0.0056582866866107199, + -0.015211731527946259, + 0.025082261844864097, + 0.039334427123337491, + -0.096220442033987982, + -0.066627474263425038, + 0.4343860564914685, + 0.78223893092049901, + 0.41530840703043026, + -0.056077313316754807, + -0.081266699680878754, + 0.026682300156053072, + 0.016068943964776348, + -0.0073461663276420935, + -0.0016294920126017326, + 0.00089231366858231456 + }; + + double hp1_a[] = { -0.00089231366858231456, + -0.0016294920126017326, + 0.0073461663276420935, + 0.016068943964776348, + -0.026682300156053072, + -0.081266699680878754, + 0.056077313316754807, + 0.41530840703043026, + -0.78223893092049901, + 0.4343860564914685, + 0.066627474263425038, + -0.096220442033987982, + -0.039334427123337491, + 0.025082261844864097, + 0.015211731527946259, + -0.0056582866866107199, + -0.0037514361572784571, + 0.0012665619292989445, + 0.00058902075624433831, + -0.00025997455248771324, + -6.2339034461007128e-05, + 3.1229875865345646e-05, + 3.2596802368833675e-06, + -1.7849850030882614e-06 + }; + + double lp2_a[] = { 0.00089231366858231456, + -0.0016294920126017326, + -0.0073461663276420935, + 0.016068943964776348, + 0.026682300156053072, + -0.081266699680878754, + -0.056077313316754807, + 0.41530840703043026, + 0.78223893092049901, + 0.4343860564914685, + -0.066627474263425038, + -0.096220442033987982, + 0.039334427123337491, + 0.025082261844864097, + -0.015211731527946259, + -0.0056582866866107199, + 0.0037514361572784571, + 0.0012665619292989445, + -0.00058902075624433831, + -0.00025997455248771324, + 6.2339034461007128e-05, + 3.1229875865345646e-05, + -3.2596802368833675e-06, + -1.7849850030882614e-06 + }; + + double hp2_a[] = { -1.7849850030882614e-06, + 3.2596802368833675e-06, + 3.1229875865345646e-05, + -6.2339034461007128e-05, + -0.00025997455248771324, + 0.00058902075624433831, + 0.0012665619292989445, + -0.0037514361572784571, + -0.0056582866866107199, + 0.015211731527946259, + 0.025082261844864097, + -0.039334427123337491, + -0.096220442033987982, + 0.066627474263425038, + 0.4343860564914685, + -0.78223893092049901, + 0.41530840703043026, + 0.056077313316754807, + -0.081266699680878754, + -0.026682300156053072, + 0.016068943964776348, + 0.0073461663276420935, + -0.0016294920126017326, + -0.00089231366858231456 + }; + for (i = 0; i < 24; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 24; + } + else if (!strcmp(name,"coif5")){ + double lp1_a[] = { -9.517657273819165e-08, + -1.6744288576823017e-07, + 2.0637618513646814e-06, + 3.7346551751414047e-06, + -2.1315026809955787e-05, + -4.1340432272512511e-05, + 0.00014054114970203437, + 0.00030225958181306315, + -0.00063813134304511142, + -0.0016628637020130838, + 0.0024333732126576722, + 0.0067641854480530832, + -0.0091642311624818458, + -0.019761778942572639, + 0.032683574267111833, + 0.041289208750181702, + -0.10557420870333893, + -0.062035963962903569, + 0.43799162617183712, + 0.77428960365295618, + 0.42156620669085149, + -0.052043163176243773, + -0.091920010559696244, + 0.02816802897093635, + 0.023408156785839195, + -0.010131117519849788, + -0.004159358781386048, + 0.0021782363581090178, + 0.00035858968789573785, + -0.00021208083980379827 + }; + + double hp1_a[] = { 0.00021208083980379827, + 0.00035858968789573785, + -0.0021782363581090178, + -0.004159358781386048, + 0.010131117519849788, + 0.023408156785839195, + -0.02816802897093635, + -0.091920010559696244, + 0.052043163176243773, + 0.42156620669085149, + -0.77428960365295618, + 0.43799162617183712, + 0.062035963962903569, + -0.10557420870333893, + -0.041289208750181702, + 0.032683574267111833, + 0.019761778942572639, + -0.0091642311624818458, + -0.0067641854480530832, + 0.0024333732126576722, + 0.0016628637020130838, + -0.00063813134304511142, + -0.00030225958181306315, + 0.00014054114970203437, + 4.1340432272512511e-05, + -2.1315026809955787e-05, + -3.7346551751414047e-06, + 2.0637618513646814e-06, + 1.6744288576823017e-07, + -9.517657273819165e-08 + }; + + double lp2_a[] = { -0.00021208083980379827, + 0.00035858968789573785, + 0.0021782363581090178, + -0.004159358781386048, + -0.010131117519849788, + 0.023408156785839195, + 0.02816802897093635, + -0.091920010559696244, + -0.052043163176243773, + 0.42156620669085149, + 0.77428960365295618, + 0.43799162617183712, + -0.062035963962903569, + -0.10557420870333893, + 0.041289208750181702, + 0.032683574267111833, + -0.019761778942572639, + -0.0091642311624818458, + 0.0067641854480530832, + 0.0024333732126576722, + -0.0016628637020130838, + -0.00063813134304511142, + 0.00030225958181306315, + 0.00014054114970203437, + -4.1340432272512511e-05, + -2.1315026809955787e-05, + 3.7346551751414047e-06, + 2.0637618513646814e-06, + -1.6744288576823017e-07, + -9.517657273819165e-08 + }; + + double hp2_a[] = { -9.517657273819165e-08, + 1.6744288576823017e-07, + 2.0637618513646814e-06, + -3.7346551751414047e-06, + -2.1315026809955787e-05, + 4.1340432272512511e-05, + 0.00014054114970203437, + -0.00030225958181306315, + -0.00063813134304511142, + 0.0016628637020130838, + 0.0024333732126576722, + -0.0067641854480530832, + -0.0091642311624818458, + 0.019761778942572639, + 0.032683574267111833, + -0.041289208750181702, + -0.10557420870333893, + 0.062035963962903569, + 0.43799162617183712, + -0.77428960365295618, + 0.42156620669085149, + 0.052043163176243773, + -0.091920010559696244, + -0.02816802897093635, + 0.023408156785839195, + 0.010131117519849788, + -0.004159358781386048, + -0.0021782363581090178, + 0.00035858968789573785, + 0.00021208083980379827 + }; + for (i = 0; i < 30; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 30; + } + + else if (!strcmp(name,"sym2")){ + double lp1_a[] = { -0.12940952255092145, + 0.22414386804185735, + 0.83651630373746899, + 0.48296291314469025 + }; + + double hp1_a[] = { -0.48296291314469025, + 0.83651630373746899, + -0.22414386804185735, + -0.12940952255092145 + }; + + double lp2_a[] = { 0.48296291314469025, + 0.83651630373746899, + 0.22414386804185735, + -0.12940952255092145 + }; + + double hp2_a[] = { -0.12940952255092145, + -0.22414386804185735, + 0.83651630373746899, + -0.48296291314469025 + + }; + for (i = 0; i < 4; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 4; + } + + else if (!strcmp(name,"sym3")){ + double lp1_a[] = { 0.035226291882100656, + -0.085441273882241486, + -0.13501102001039084, + 0.45987750211933132, + 0.80689150931333875, + 0.33267055295095688 + + }; + + double hp1_a[] = { -0.33267055295095688, + 0.80689150931333875, + -0.45987750211933132, + -0.13501102001039084, + 0.085441273882241486, + 0.035226291882100656 + }; + + double lp2_a[] = { 0.33267055295095688, + 0.80689150931333875, + 0.45987750211933132, + -0.13501102001039084, + -0.085441273882241486, + 0.035226291882100656 + }; + + double hp2_a[] = { 0.035226291882100656, + 0.085441273882241486, + -0.13501102001039084, + -0.45987750211933132, + 0.80689150931333875, + -0.33267055295095688 + }; + for (i = 0; i < 6; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 6; + } + + else if (!strcmp(name,"sym4")){ + double lp1_a[] = { -0.075765714789273325, + -0.02963552764599851, + 0.49761866763201545, + 0.80373875180591614, + 0.29785779560527736, + -0.099219543576847216, + -0.012603967262037833, + 0.032223100604042702 + }; + + double hp1_a[] = { -0.032223100604042702, + -0.012603967262037833, + 0.099219543576847216, + 0.29785779560527736, + -0.80373875180591614, + 0.49761866763201545, + 0.02963552764599851, + -0.075765714789273325 + }; + + double lp2_a[] = { 0.032223100604042702, + -0.012603967262037833, + -0.099219543576847216, + 0.29785779560527736, + 0.80373875180591614, + 0.49761866763201545, + -0.02963552764599851, + -0.075765714789273325 + }; + + double hp2_a[] = { -0.075765714789273325, + 0.02963552764599851, + 0.49761866763201545, + -0.80373875180591614, + 0.29785779560527736, + 0.099219543576847216, + -0.012603967262037833, + -0.032223100604042702 + }; + for (i = 0; i < 8; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 8; + } + + else if (!strcmp(name,"sym5")){ + double lp1_a[] = { 0.027333068345077982, + 0.029519490925774643, + -0.039134249302383094, + 0.1993975339773936, + 0.72340769040242059, + 0.63397896345821192, + 0.016602105764522319, + -0.17532808990845047, + -0.021101834024758855, + 0.019538882735286728 + }; + + double hp1_a[] = { -0.019538882735286728, + -0.021101834024758855, + 0.17532808990845047, + 0.016602105764522319, + -0.63397896345821192, + 0.72340769040242059, + -0.1993975339773936, + -0.039134249302383094, + -0.029519490925774643, + 0.027333068345077982 + }; + + double lp2_a[] = { 0.019538882735286728, + -0.021101834024758855, + -0.17532808990845047, + 0.016602105764522319, + 0.63397896345821192, + 0.72340769040242059, + 0.1993975339773936, + -0.039134249302383094, + 0.029519490925774643, + 0.027333068345077982 + + }; + + double hp2_a[] = { 0.027333068345077982, + -0.029519490925774643, + -0.039134249302383094, + -0.1993975339773936, + 0.72340769040242059, + -0.63397896345821192, + 0.016602105764522319, + 0.17532808990845047, + -0.021101834024758855, + -0.019538882735286728 + }; + for (i = 0; i < 10; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 10; + } + + else if (!strcmp(name,"sym6")){ + double lp1_a[] = { 0.015404109327027373, + 0.0034907120842174702, + -0.11799011114819057, + -0.048311742585632998, + 0.49105594192674662, + 0.787641141030194, + 0.3379294217276218, + -0.072637522786462516, + -0.021060292512300564, + 0.044724901770665779, + 0.0017677118642428036, + -0.007800708325034148 + }; + + double hp1_a[] = { 0.007800708325034148, + 0.0017677118642428036, + -0.044724901770665779, + -0.021060292512300564, + 0.072637522786462516, + 0.3379294217276218, + -0.787641141030194, + 0.49105594192674662, + 0.048311742585632998, + -0.11799011114819057, + -0.0034907120842174702, + 0.015404109327027373 + }; + + double lp2_a[] = { -0.007800708325034148, + 0.0017677118642428036, + 0.044724901770665779, + -0.021060292512300564, + -0.072637522786462516, + 0.3379294217276218, + 0.787641141030194, + 0.49105594192674662, + -0.048311742585632998, + -0.11799011114819057, + 0.0034907120842174702, + 0.015404109327027373 + }; + + double hp2_a[] = { 0.015404109327027373, + -0.0034907120842174702, + -0.11799011114819057, + 0.048311742585632998, + 0.49105594192674662, + -0.787641141030194, + 0.3379294217276218, + 0.072637522786462516, + -0.021060292512300564, + -0.044724901770665779, + 0.0017677118642428036, + 0.007800708325034148 + }; + for (i = 0; i < 12; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 12; + } + + else if (!strcmp(name,"sym7")){ + double lp1_a[] = { 0.0026818145682578781, + -0.0010473848886829163, + -0.01263630340325193, + 0.03051551316596357, + 0.067892693501372697, + -0.049552834937127255, + 0.017441255086855827, + 0.5361019170917628, + 0.76776431700316405, + 0.28862963175151463, + -0.14004724044296152, + -0.10780823770381774, + 0.0040102448715336634, + 0.010268176708511255 + }; + + double hp1_a[] = { -0.010268176708511255, + 0.0040102448715336634, + 0.10780823770381774, + -0.14004724044296152, + -0.28862963175151463, + 0.76776431700316405, + -0.5361019170917628, + 0.017441255086855827, + 0.049552834937127255, + 0.067892693501372697, + -0.03051551316596357, + -0.01263630340325193, + 0.0010473848886829163, + 0.0026818145682578781 + + }; + + double lp2_a[] = { 0.010268176708511255, + 0.0040102448715336634, + -0.10780823770381774, + -0.14004724044296152, + 0.28862963175151463, + 0.76776431700316405, + 0.5361019170917628, + 0.017441255086855827, + -0.049552834937127255, + 0.067892693501372697, + 0.03051551316596357, + -0.01263630340325193, + -0.0010473848886829163, + 0.0026818145682578781 + }; + + double hp2_a[] = { 0.0026818145682578781, + 0.0010473848886829163, + -0.01263630340325193, + -0.03051551316596357, + 0.067892693501372697, + 0.049552834937127255, + 0.017441255086855827, + -0.5361019170917628, + 0.76776431700316405, + -0.28862963175151463, + -0.14004724044296152, + 0.10780823770381774, + 0.0040102448715336634, + -0.010268176708511255 + }; + for (i = 0; i < 14; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 14; + } + + else if (!strcmp(name,"sym8")){ + double lp1_a[] = { -0.0033824159510061256, + -0.00054213233179114812, + 0.031695087811492981, + 0.0076074873249176054, + -0.14329423835080971, + -0.061273359067658524, + 0.48135965125837221, + 0.77718575170052351, + 0.3644418948353314, + -0.051945838107709037, + -0.027219029917056003, + 0.049137179673607506, + 0.0038087520138906151, + -0.014952258337048231, + -0.0003029205147213668, + 0.0018899503327594609 + }; + + double hp1_a[] = { -0.0018899503327594609, + -0.0003029205147213668, + 0.014952258337048231, + 0.0038087520138906151, + -0.049137179673607506, + -0.027219029917056003, + 0.051945838107709037, + 0.3644418948353314, + -0.77718575170052351, + 0.48135965125837221, + 0.061273359067658524, + -0.14329423835080971, + -0.0076074873249176054, + 0.031695087811492981, + 0.00054213233179114812, + -0.0033824159510061256 + }; + + double lp2_a[] = { 0.0018899503327594609, + -0.0003029205147213668, + -0.014952258337048231, + 0.0038087520138906151, + 0.049137179673607506, + -0.027219029917056003, + -0.051945838107709037, + 0.3644418948353314, + 0.77718575170052351, + 0.48135965125837221, + -0.061273359067658524, + -0.14329423835080971, + 0.0076074873249176054, + 0.031695087811492981, + -0.00054213233179114812, + -0.0033824159510061256 + }; + + double hp2_a[] = { -0.0033824159510061256, + 0.00054213233179114812, + 0.031695087811492981, + -0.0076074873249176054, + -0.14329423835080971, + 0.061273359067658524, + 0.48135965125837221, + -0.77718575170052351, + 0.3644418948353314, + 0.051945838107709037, + -0.027219029917056003, + -0.049137179673607506, + 0.0038087520138906151, + 0.014952258337048231, + -0.0003029205147213668, + -0.0018899503327594609 + }; + for (i = 0; i < 16; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 16; + } + + else if (!strcmp(name,"sym9")){ + double lp1_a[] = { 0.0014009155259146807, + 0.00061978088898558676, + -0.013271967781817119, + -0.01152821020767923, + 0.03022487885827568, + 0.00058346274612580684, + -0.054568958430834071, + 0.238760914607303, + 0.717897082764412, + 0.61733844914093583, + 0.035272488035271894, + -0.19155083129728512, + -0.018233770779395985, + 0.06207778930288603, + 0.0088592674934004842, + -0.010264064027633142, + -0.00047315449868008311, + 0.0010694900329086053 + }; + + double hp1_a[] = { -0.0010694900329086053, + -0.00047315449868008311, + 0.010264064027633142, + 0.0088592674934004842, + -0.06207778930288603, + -0.018233770779395985, + 0.19155083129728512, + 0.035272488035271894, + -0.61733844914093583, + 0.717897082764412, + -0.238760914607303, + -0.054568958430834071, + -0.00058346274612580684, + 0.03022487885827568, + 0.01152821020767923, + -0.013271967781817119, + -0.00061978088898558676, + 0.0014009155259146807 + }; + + double lp2_a[] = { 0.0010694900329086053, + -0.00047315449868008311, + -0.010264064027633142, + 0.0088592674934004842, + 0.06207778930288603, + -0.018233770779395985, + -0.19155083129728512, + 0.035272488035271894, + 0.61733844914093583, + 0.717897082764412, + 0.238760914607303, + -0.054568958430834071, + 0.00058346274612580684, + 0.03022487885827568, + -0.01152821020767923, + -0.013271967781817119, + 0.00061978088898558676, + 0.0014009155259146807 + }; + + double hp2_a[] = { 0.0014009155259146807, + -0.00061978088898558676, + -0.013271967781817119, + 0.01152821020767923, + 0.03022487885827568, + -0.00058346274612580684, + -0.054568958430834071, + -0.238760914607303, + 0.717897082764412, + -0.61733844914093583, + 0.035272488035271894, + 0.19155083129728512, + -0.018233770779395985, + -0.06207778930288603, + 0.0088592674934004842, + 0.010264064027633142, + -0.00047315449868008311, + -0.0010694900329086053 + }; + for (i = 0; i < 18; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 18; + } + + else if (!strcmp(name,"sym10")){ + double lp1_a[] = { 0.00077015980911449011, + 9.5632670722894754e-05, + -0.0086412992770224222, + -0.0014653825813050513, + 0.045927239231092203, + 0.011609893903711381, + -0.15949427888491757, + -0.070880535783243853, + 0.47169066693843925, + 0.7695100370211071, + 0.38382676106708546, + -0.035536740473817552, + -0.0319900568824278, + 0.049994972077376687, + 0.0057649120335819086, + -0.02035493981231129, + -0.00080435893201654491, + 0.0045931735853118284, + 5.7036083618494284e-05, + -0.00045932942100465878 + }; + + double hp1_a[] = { 0.00045932942100465878, + 5.7036083618494284e-05, + -0.0045931735853118284, + -0.00080435893201654491, + 0.02035493981231129, + 0.0057649120335819086, + -0.049994972077376687, + -0.0319900568824278, + 0.035536740473817552, + 0.38382676106708546, + -0.7695100370211071, + 0.47169066693843925, + 0.070880535783243853, + -0.15949427888491757, + -0.011609893903711381, + 0.045927239231092203, + 0.0014653825813050513, + -0.0086412992770224222, + -9.5632670722894754e-05, + 0.00077015980911449011 + }; + + double lp2_a[] = { -0.00045932942100465878, + 5.7036083618494284e-05, + 0.0045931735853118284, + -0.00080435893201654491, + -0.02035493981231129, + 0.0057649120335819086, + 0.049994972077376687, + -0.0319900568824278, + -0.035536740473817552, + 0.38382676106708546, + 0.7695100370211071, + 0.47169066693843925, + -0.070880535783243853, + -0.15949427888491757, + 0.011609893903711381, + 0.045927239231092203, + -0.0014653825813050513, + -0.0086412992770224222, + 9.5632670722894754e-05, + 0.00077015980911449011 + }; + + double hp2_a[] = { 0.00077015980911449011, + -9.5632670722894754e-05, + -0.0086412992770224222, + 0.0014653825813050513, + 0.045927239231092203, + -0.011609893903711381, + -0.15949427888491757, + 0.070880535783243853, + 0.47169066693843925, + -0.7695100370211071, + 0.38382676106708546, + 0.035536740473817552, + -0.0319900568824278, + -0.049994972077376687, + 0.0057649120335819086, + 0.02035493981231129, + -0.00080435893201654491, + -0.0045931735853118284, + 5.7036083618494284e-05, + 0.00045932942100465878 + }; + for (i = 0; i < 20; ++i) { + lp1[i] = lp1_a[i]; + hp1[i] = hp1_a[i]; + lp2[i] = lp2_a[i]; + hp2[i] = hp2_a[i]; + } + return 20; + } + + + else { + printf("\n Filter Not in Database \n"); + return -1; + } + + return 0; +} \ No newline at end of file diff --git a/src/wavefilt.h b/src/wavefilt.h new file mode 100644 index 0000000..24f16e5 --- /dev/null +++ b/src/wavefilt.h @@ -0,0 +1,22 @@ +#ifndef WAVEFILT_H_ +#define WAVEFILT_H_ + +#include +#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_ */ \ No newline at end of file diff --git a/src/wavelib.c b/src/wavelib.c new file mode 100644 index 0000000..cfd271f --- /dev/null +++ b/src/wavelib.c @@ -0,0 +1,1274 @@ +/* + Copyright (c) 2014, Rafat Hussain +*/ + +#include "wavelib.h" + +wave_object wave_init(char* wname) { + wave_object obj = NULL; + int retval; + retval = 0; + + if (wname != NULL) { + retval = filtlength(wname); + //obj->filtlength = retval; + //strcopy(obj->wname, wname); + } + + obj = (wave_object)malloc(sizeof(struct wave_set) + sizeof(double)* 4 * retval); + + obj->filtlength = retval; + obj->lpd_len = obj->hpd_len = obj->lpr_len = obj->hpr_len = obj->filtlength; + strcpy(obj->wname, wname); + if (wname != NULL) { + filtcoef(wname,obj->params,obj->params+retval,obj->params+2*retval,obj->params+3*retval); + } + obj->lpd = &obj->params[0]; + obj->hpd = &obj->params[retval]; + obj->lpr = &obj->params[2 * retval]; + obj->hpr = &obj->params[3 * retval]; + return obj; +} + +wt_object wt_init(wave_object wave,char* method, int siglength,int J) { + int size,i,MaxIter; + wt_object obj = NULL; + + size = wave->filtlength; + + if (J > 100) { + printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n"); + exit(-1); + } + + MaxIter = wmaxiter(siglength, size); + + if (J > MaxIter) { + printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter); + exit(-1); + } + + if (method == NULL) { + obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength + 2 * J * (size+1))); + obj->outlength = siglength + 2 * J * (size + 1); // Default + strcpy(obj->ext, "sym"); // Default + } + else if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) { + obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength + 2 * J * (size + 1))); + obj->outlength = siglength + 2 * J * (size + 1); // Default + strcpy(obj->ext, "sym"); // Default + } + else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) { + if (!testSWTlength(siglength, J)) { + printf("\n For SWT the signal length must be a multiple of 2^J. \n"); + exit(-1); + } + + obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1))); + obj->outlength = siglength * (J + 1); // Default + strcpy(obj->ext, "per"); // Default + } + else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) { + + if (!strstr(wave->wname,"db")) { + if (!strstr(wave->wname, "sym")) { + if (!strstr(wave->wname, "coif")) { + printf("\n MODWT is only implemented for orthogonal wavelet families - db, sym and coif \n"); + exit(-1); + } + } + } + + obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1))); + obj->outlength = siglength * (J + 1); // Default + strcpy(obj->ext, "per"); // Default + } + + obj->wave = wave; + obj->siglength = siglength; + obj->J = J; + obj->MaxIter = MaxIter; + strcpy(obj->method, method); + + if (siglength % 2 == 0) { + obj->even = 1; + } + else { + obj->even = 0; + } + + obj->cobj = NULL; + + strcpy(obj->cmethod, "direct"); // Default + obj->cfftset = 0; + obj->lenlength = J + 2; + obj->output = &obj->params[0]; + if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) { + for (i = 0; i < siglength + 2 * J * (size + 1); ++i) { + obj->params[i] = 0.0; + } + } + else if (!strcmp(method, "swt") || !strcmp(method, "SWT") || !strcmp(method, "modwt") || !strcmp(method, "MODWT")) { + for (i = 0; i < siglength * (J + 1); ++i) { + obj->params[i] = 0.0; + } + } + //wave_summary(obj->wave); + + return obj; +} + +static void wconv(wt_object wt, double *sig, int N, double *filt, int L, double *oup) { + if (!strcmp(wt->cmethod,"direct")) { + conv_direct(sig, N, filt, L, oup); + } + else if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) { + if (wt->cfftset == 0) { + wt->cobj = conv_init(N, L); + conv_fft(wt->cobj, sig, filt, oup); + free_conv(wt->cobj); + } + else { + conv_fft(wt->cobj, sig, filt, oup); + } + } + else { + printf("Convolution Only accepts two methods - direct and fft"); + exit(-1); + } +} + + +static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { + int l,l2,isodd,i,t,len_avg; + + len_avg = wt->wave->lpd_len; + l2 = len_avg / 2; + isodd = N % 2; + + for (i = 0; i < len_cA; ++i) { + t = 2 * i + l2; + cA[i] = 0.0; + cD[i] = 0.0; + for (l = 0; l < len_avg; ++l) { + if ((t - l) >= l2 && (t - l) < N) { + cA[i] += wt->wave->lpd[l] * inp[t - l]; + cD[i] += wt->wave->hpd[l] * inp[t - l]; + } + else if ((t - l) < l2 && (t - l) >= 0) { + cA[i] += wt->wave->lpd[l] * inp[t - l]; + cD[i] += wt->wave->hpd[l] * inp[t - l]; + } + else if ((t - l) < 0) { + cA[i] += wt->wave->lpd[l] * inp[t - l + N]; + cD[i] += wt->wave->hpd[l] * inp[t - l + N]; + } + else if ((t - l) >= N && isodd == 0) { + cA[i] += wt->wave->lpd[l] * inp[t - l - N]; + cD[i] += wt->wave->hpd[l] * inp[t - l - N]; + } + else if ((t - l) >= N && isodd == 1) { + if (t - l != N) { + cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)]; + cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)]; + } + else { + cA[i] += wt->wave->lpd[l] * inp[N - 1]; + cD[i] += wt->wave->hpd[l] * inp[N - 1]; + } + } + + } + } + + +} + +static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { + int i,l, t, len_avg; + + len_avg = wt->wave->lpd_len; + + for (i = 0; i < len_cA; ++i) { + t = 2 * i + 1; + cA[i] = 0.0; + cD[i] = 0.0; + for (l = 0; l < len_avg; ++l) { + if ((t - l) >= 0 && (t - l) < N) { + cA[i] += wt->wave->lpd[l] * inp[t - l]; + cD[i] += wt->wave->hpd[l] * inp[t - l]; + } + else if ((t - l) < 0) { + cA[i] += wt->wave->lpd[l] * inp[-t + l - 1]; + cD[i] += wt->wave->hpd[l] * inp[-t + l - 1]; + } + else if ((t - l) >= N) { + cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1]; + cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1]; + } + + } + } + + +} + +static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, int len_cA, double *cD, int len_cD) { + int len_avg,D,lf; + double *signal,*cA_undec; + len_avg = (wt->wave->lpd_len + wt->wave->hpd_len) / 2; + //len_sig = 2 * (int)ceil((double)len_sig / 2.0); + + D = 2; + + if (!strcmp(wt->ext, "per")) { + signal = (double*)malloc(sizeof(double)* (len_sig + len_avg + (len_sig % 2))); + + len_sig = per_ext(sig, len_sig, len_avg / 2, signal); + + cA_undec = (double*)malloc(sizeof(double)* (len_sig + len_avg + wt->wave->lpd_len - 1)); + + if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + wt->cobj = conv_init(len_sig + len_avg, wt->wave->lpd_len); + wt->cfftset = 1; + } + else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) { + printf("Decomposition Filters must have the same length."); + exit(-1); + } + + wconv(wt, signal, len_sig + len_avg, wt->wave->lpd, wt->wave->lpd_len, cA_undec); + + downsamp(cA_undec + len_avg, len_sig, D, cA); + + wconv(wt, signal, len_sig + len_avg, wt->wave->hpd, wt->wave->hpd_len, cA_undec); + + downsamp(cA_undec + len_avg, len_sig, D, cD); + } + else if (!strcmp(wt->ext, "sym")) { + //printf("\n YES %s \n", wt->ext); + lf = wt->wave->lpd_len;// lpd and hpd have the same length + + signal = (double*)malloc(sizeof(double)* (len_sig + 2 * (lf - 1))); + + len_sig = symm_ext(sig, len_sig, lf - 1, signal); + + cA_undec = (double*)malloc(sizeof(double)* (len_sig + 3 * (lf - 1))); + + if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + wt->cobj = conv_init(len_sig + 2 * (lf - 1), lf); + wt->cfftset = 1; + } + else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) { + printf("Decomposition Filters must have the same length."); + exit(-1); + } + + + wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->lpd, wt->wave->lpd_len, cA_undec); + + downsamp(cA_undec + lf, len_sig + lf - 2, D, cA); + + wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->hpd, wt->wave->hpd_len, cA_undec); + + downsamp(cA_undec + lf, len_sig + lf - 2, D, cD); + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } + + + if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + free_conv(wt->cobj); + wt->cfftset = 0; + } + + free(signal); + free(cA_undec); +} + +void dwt(wt_object wt,double *inp) { + int i,J,temp_len,iter,N,lp; + int len_cA; + double *orig,*orig2; + + temp_len = wt->siglength; + J = wt->J; + wt->length[J + 1] = temp_len; + wt->outlength = 0; + if ((temp_len % 2) == 0) { + wt->zpad = 0; + orig = (double*)malloc(sizeof(double)* temp_len); + orig2 = (double*)malloc(sizeof(double)* temp_len); + } + else { + wt->zpad = 1; + temp_len++; + orig = (double*)malloc(sizeof(double)* temp_len); + orig2 = (double*)malloc(sizeof(double)* temp_len); + } + + for (i = 0; i < wt->siglength; ++i) { + orig[i] = inp[i]; + } + + if (wt->zpad == 1) { + orig[temp_len - 1] = orig[temp_len - 2]; + } + + N = temp_len; + lp = wt->wave->lpd_len; + + if (!strcmp(wt->ext,"per")) { + i = J; + while (i > 0) { + N = (int)ceil((double)N / 2.0); + wt->length[i] = N; + wt->outlength += wt->length[i]; + i--; + } + wt->length[0] = wt->length[1]; + wt->outlength += wt->length[0]; + N = wt->outlength; + + for (iter = 0; iter < J; ++iter) { + len_cA = wt->length[J - iter]; + N -= len_cA; + if ( !strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT") ) { + dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + } + else { + dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + } + temp_len = wt->length[J - iter]; + if (iter == J - 1) { + for (i = 0; i < len_cA; ++i) { + wt->params[i] = orig2[i]; + } + } + else { + for (i = 0; i < len_cA; ++i) { + orig[i] = orig2[i]; + } + } + } + } + else if (!strcmp(wt->ext,"sym")) { + //printf("\n YES %s \n", wt->ext); + i = J; + while (i > 0) { + N = N + lp - 2; + N = (int) ceil((double)N / 2.0); + wt->length[i] = N; + wt->outlength += wt->length[i]; + i--; + } + wt->length[0] = wt->length[1]; + wt->outlength += wt->length[0]; + N = wt->outlength; + + for (iter = 0; iter < J; ++iter) { + len_cA = wt->length[J - iter]; + N -= len_cA; + if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) { + dwt1(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + } + else { + dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N, len_cA); + } + temp_len = wt->length[J - iter]; + + if (iter == J - 1) { + for (i = 0; i < len_cA; ++i) { + wt->params[i] = orig2[i]; + } + } + else { + for (i = 0; i < len_cA; ++i) { + orig[i] = orig2[i]; + } + } + } + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } + + free(orig); + free(orig2); +} + +static void idwt1(wt_object wt,double *temp, double *cA_up,double *cA, int len_cA,double *cD,int len_cD,double *X_lp,double *X_hp,double *X) { + int len_avg, N, U,N2,i; + + len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; + N = 2 * len_cD; + U = 2; + + upsamp2(cA, len_cA, U, cA_up); + + per_ext(cA_up, 2 * len_cA, len_avg / 2, temp); + + N2 = 2 * len_cA + len_avg; + + if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + wt->cobj = conv_init(N2, len_avg); + wt->cfftset = 1; + } + else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) { + printf("Decomposition Filters must have the same length."); + exit(-1); + } + + wconv(wt, temp, N2, wt->wave->lpr, len_avg, X_lp); + + upsamp2(cD, len_cD, U, cA_up); + + per_ext(cA_up, 2 * len_cD, len_avg / 2, temp); + + N2 = 2 * len_cD + len_avg; + + wconv(wt, temp, N2, wt->wave->hpr, len_avg, X_hp); + + + for (i = len_avg - 1; i < N + len_avg - 1; ++i) { + X[i - len_avg + 1] = X_lp[i] + X_hp[i]; + } + + if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + free_conv(wt->cobj); + wt->cfftset = 0; + } + +} + +static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { + int len_avg, N, i,l,m,n,t,v,l2; + + len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; + N = 2 * len_cA; + l2 = len_avg / 2; + m = -2; + n = -1; + + for (i = 0; i < len_cA + l2 - 1; ++i) { + m += 2; + n += 2; + X[m] = 0.0; + X[n] = 0.0; + for (l = 0; l < l2; ++l) { + t = 2 * l; + if ((i - l) >= 0 && (i - l) < len_cA) { + X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l]; + X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l]; + } + else if ((i - l) >= len_cA && (i-l) < len_cA + len_avg - 1) { + X[m] += wt->wave->lpr[t] * cA[i - l - len_cA] + wt->wave->hpr[t] * cD[i - l - len_cA]; + X[n] += wt->wave->lpr[t + 1] * cA[i - l - len_cA] + wt->wave->hpr[t + 1] * cD[i - l - len_cA]; + } + else if ((i - l) < 0 && (i-l) > -l2) { + X[m] += wt->wave->lpr[t] * cA[len_cA + i - l] + wt->wave->hpr[t] * cD[len_cA + i - l]; + X[n] += wt->wave->lpr[t + 1] * cA[len_cA + i - l] + wt->wave->hpr[t + 1] * cD[len_cA + i - l]; + } + } + } +} + +static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, int len_cD, double *X) { + int len_avg, N, i, l, m, n, t, v,l2; + + len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; + N = 2 * len_cA - 1; + m = -2; + n = -1; + l2 = (len_avg - 2) / 2; + + for (v = 0; v < len_cA; ++v) { + i = v; + m += 2; + n += 2; + X[m] = 0.0; + X[n] = 0.0; + for (l = 0; l < len_avg / 2; ++l) { + t = 2 * l; + if ((i - l) >= 0 && (i - l) < len_cA) { + X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l]; + X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l]; + } + } + } +} + + +void idwt(wt_object wt, double *dwtop) { + int J,U,i,lf,N,N2,iter,k; + int app_len, det_len; + double *cA_up, *X_lp, *X_hp,*out,*temp; + + J = wt->J; + U = 2; + app_len = wt->length[0]; + out = (double*)malloc(sizeof(double)* (wt->siglength + 1)); + if (!strcmp(wt->ext, "per") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + app_len = wt->length[0]; + det_len = wt->length[1]; + N = 2 * wt->length[J]; + lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; + + cA_up = (double*)malloc(sizeof(double)* N); + temp = (double*)malloc(sizeof(double)* (N + lf)); + X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); + X_hp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); + iter = app_len; + + for (i = 0; i < app_len; ++i) { + out[i] = wt->output[i]; + } + + for (i = 0; i < J; ++i) { + + 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); + for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) { + out[k - lf/2 + 1] = X_lp[k]; + } + */ + iter += det_len; + det_len = wt->length[i + 2]; + } + free(cA_up); + free(X_lp); + free(X_hp); + free(temp); + + } + else if (!strcmp(wt->ext, "per") && !strcmp(wt->cmethod, "direct")) { + app_len = wt->length[0]; + det_len = wt->length[1]; + N = 2 * wt->length[J]; + lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; + + X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); + iter = app_len; + + for (i = 0; i < app_len; ++i) { + out[i] = wt->output[i]; + } + + for (i = 0; i < J; ++i) { + + //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); + for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) { + out[k - lf/2 + 1] = X_lp[k]; + } + + iter += det_len; + det_len = wt->length[i + 2]; + } + + free(X_lp); + + } + else if (!strcmp(wt->ext, "sym") && !strcmp(wt->cmethod, "direct")) { + app_len = wt->length[0]; + det_len = wt->length[1]; + N = 2 * wt->length[J] - 1; + lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2; + + X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1)); + iter = app_len; + + for (i = 0; i < app_len; ++i) { + out[i] = wt->output[i]; + } + + for (i = 0; i < J; ++i) { + + //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); + for (k = lf-2; k < 2 * det_len; ++k) { + out[k - lf + 2] = X_lp[k]; + } + + iter += det_len; + det_len = wt->length[i + 2]; + } + + free(X_lp); + + } + else if (!strcmp(wt->ext, "sym") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + lf = wt->wave->lpd_len;// lpd and hpd have the same length + + N = 2 * wt->length[J] - 1; + cA_up = (double*)malloc(sizeof(double)* N); + X_lp = (double*)malloc(sizeof(double)* (N + lf - 1)); + X_hp = (double*)malloc(sizeof(double)* (N + lf - 1)); + + for (i = 0; i < app_len; ++i) { + out[i] = wt->output[i]; + } + + iter = app_len; + + for (i = 0; i < J; ++i) { + det_len = wt->length[i + 1]; + upsamp(out, det_len, U, cA_up); + N2 = 2 * wt->length[i + 1] - 1; + + if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + wt->cobj = conv_init(N2, lf); + wt->cfftset = 1; + } + else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) { + printf("Decomposition Filters must have the same length."); + exit(-1); + } + + wconv(wt, cA_up, N2, wt->wave->lpr, lf, X_lp); + + upsamp(wt->output + iter, det_len, U, cA_up); + + wconv(wt, cA_up, N2, wt->wave->hpr, lf, X_hp); + + for (k = lf - 2; k < N2 + 1; ++k) { + out[k - lf + 2] = X_lp[k] + X_hp[k]; + } + iter += det_len; + if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + free_conv(wt->cobj); + wt->cfftset = 0; + } + } + + free(cA_up); + free(X_lp); + free(X_hp); + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } + + for (i = 0; i < wt->siglength; ++i) { + dwtop[i] = out[i]; + } + + + free(out); + +} + +static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD, int len_cD) { + int l, l2, isodd, i, t, len_avg,j; + + len_avg = M * wt->wave->lpd_len; + l2 = len_avg / 2; + isodd = N % 2; + + for (i = 0; i < len_cA; ++i) { + t = i + l2; + cA[i] = 0.0; + cD[i] = 0.0; + l = -1; + for (j = 0; j < len_avg; j+=M) { + l++; + while (j >= len_cA) { + j -= len_cA; + } + if ((t - j) >= l2 && (t -j) < N) { + cA[i] += wt->wave->lpd[l] * inp[t - j]; + cD[i] += wt->wave->hpd[l] * inp[t - j]; + } + else if ((t - j) < l2 && (t - j) >= 0) { + cA[i] += wt->wave->lpd[l] * inp[t - j]; + cD[i] += wt->wave->hpd[l] * inp[t - j]; + } + else if ((t - j) < 0) { + cA[i] += wt->wave->lpd[l] * inp[t - j + N]; + cD[i] += wt->wave->hpd[l] * inp[t - j + N]; + } + else if ((t - j) >= N && isodd == 0) { + cA[i] += wt->wave->lpd[l] * inp[t - j - N]; + cD[i] += wt->wave->hpd[l] * inp[t - j - N]; + } + else if ((t - j) >= N && isodd == 1) { + if (t - l != N) { + cA[i] += wt->wave->lpd[l] * inp[t - j - (N + 1)]; + cD[i] += wt->wave->hpd[l] * inp[t - j - (N + 1)]; + } + else { + cA[i] += wt->wave->lpd[l] * inp[N - 1]; + cD[i] += wt->wave->hpd[l] * inp[N - 1]; + } + } + + } + } + + +} + +static void swt_fft(wt_object wt, double *inp) { + int i, J, temp_len, iter, M, N, len_filt; + int lenacc; + double *low_pass, *high_pass,*sig,*cA,*cD; + + temp_len = wt->siglength; + J = wt->J; + wt->length[0] = wt->length[J] = temp_len; + wt->outlength = wt->length[J+1] = (J + 1) * temp_len; + M = 1; + for (iter = 1; iter < J; ++iter) { + M = 2 * M; + wt->length[iter] = temp_len; + } + + len_filt = wt->wave->filtlength; + + low_pass = (double*)malloc(sizeof(double)* M * len_filt); + high_pass = (double*)malloc(sizeof(double)* M * len_filt); + sig = (double*)malloc(sizeof(double)* (M * len_filt + temp_len + (temp_len%2))); + cA = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1); + cD = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1); + + M = 1; + + for (i = 0; i < temp_len; ++i) { + wt->params[i] = inp[i]; + } + + lenacc = wt->outlength; + + for (iter = 0; iter < J; ++iter) { + lenacc -= temp_len; + if (iter > 0) { + M = 2 * M; + N = M * len_filt; + upsamp2(wt->wave->lpd, wt->wave->lpd_len, M, low_pass); + upsamp2(wt->wave->hpd, wt->wave->hpd_len, M, high_pass); + } + else { + N = len_filt; + for (i = 0; i < N; ++i) { + low_pass[i] = wt->wave->lpd[i]; + high_pass[i] = wt->wave->hpd[i]; + } + } + + //swt_per(wt,M, wt->params, temp_len, cA, temp_len, cD,temp_len); + + per_ext(wt->params, temp_len, N / 2, sig); + + if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + wt->cobj = conv_init(N + temp_len + (temp_len % 2), N); + wt->cfftset = 1; + } + else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) { + printf("Decomposition Filters must have the same length."); + exit(-1); + } + + + wconv(wt, sig, N + temp_len + (temp_len % 2), low_pass, N, cA); + + wconv(wt, sig, N + temp_len + (temp_len % 2), high_pass, N, cD); + + if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + free_conv(wt->cobj); + wt->cfftset = 0; + } + + for (i = 0; i < temp_len; ++i) { + wt->params[i] = cA[N + i]; + wt->params[lenacc + i] = cD[N + i]; + } + + + } + + free(low_pass); + free(high_pass); + free(sig); + free(cA); + free(cD); +} + +static void swt_direct(wt_object wt, double *inp) { + int i, J, temp_len, iter, M, N; + int lenacc, len_filt; + double *cA, *cD; + + temp_len = wt->siglength; + J = wt->J; + wt->length[0] = wt->length[J] = temp_len; + wt->outlength = wt->length[J + 1] = (J + 1) * temp_len; + len_filt = wt->wave->filtlength; + M = 1; + for (iter = 1; iter < J; ++iter) { + M = 2 * M; + wt->length[iter] = temp_len; + } + + + cA = (double*)malloc(sizeof(double)* temp_len); + cD = (double*)malloc(sizeof(double)* temp_len); + + M = 1; + + for (i = 0; i < temp_len; ++i) { + wt->params[i] = inp[i]; + } + + lenacc = wt->outlength; + + for (iter = 0; iter < J; ++iter) { + lenacc -= temp_len; + if (iter > 0) { + 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); + + + for (i = 0; i < temp_len; ++i) { + wt->params[i] = cA[i]; + wt->params[lenacc + i] = cD[i]; + } + + } + + free(cA); + free(cD); + +} + + +void swt(wt_object wt, double *inp) { + if (!strcmp(wt->method, "swt") && !strcmp(wt->cmethod, "direct") ) { + swt_direct(wt,inp); + } + else if (!strcmp(wt->method, "swt") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + swt_fft(wt, inp); + } + else { + printf("SWT Only accepts two methods - direct and fft"); + exit(-1); + } +} + +void iswt(wt_object wt, double *swtop) { + int N, lf, iter,i,J,index,value,count,len; + int index_shift,len0,U,N1,N2,index2; + double *low_pass, *high_pass; + double *appx1, *det1,*appx_sig,*det_sig,*cL0,*cH0,*tempx,*oup00L,*oup00H,*oup00,*oup01,*appx2,*det2; + + N = wt->siglength; + J = wt->J; + U = 2; + lf = wt->wave->lpr_len; + + appx_sig = (double*)malloc(sizeof(double)* N); + det_sig = (double*)malloc(sizeof(double)* N); + appx1 = (double*)malloc(sizeof(double)* N); + det1 = (double*)malloc(sizeof(double)* N); + appx2 = (double*)malloc(sizeof(double)* N); + det2 = (double*)malloc(sizeof(double)* N); + tempx = (double*)malloc(sizeof(double)* N); + cL0 = (double*)malloc(sizeof(double)* (N + (N%2) + lf)); + cH0 = (double*)malloc(sizeof(double)* (N + (N % 2) + lf)); + oup00L = (double*)malloc(sizeof(double)* (N + 2 * lf)); + oup00H = (double*)malloc(sizeof(double)* (N + 2 * lf)); + oup00 = (double*)malloc(sizeof(double)* N); + oup01 = (double*)malloc(sizeof(double)* N); + + + + for (iter = 0; iter < J; ++iter) { + for (i = 0; i < N; ++i) { + swtop[i] = 0.0; + } + if (iter == 0) { + for (i = 0; i < N; ++i) { + appx_sig[i] = wt->output[i]; + det_sig[i] = wt->output[N + i]; + } + } + else { + for (i = 0; i < N; ++i) { + det_sig[i] = wt->output[(iter + 1) * N + i]; + } + } + + value = (int)pow(2.0, (double) (J - 1 - iter)); + + for (count = 0; count < value; count++) { + len = 0; + for (index = count; index < N; index += value) { + appx1[len] = appx_sig[index]; + det1[len] = det_sig[index]; + len++; + } + + + //SHIFT 0 + len0 = 0; + + for (index_shift = 0; index_shift < len; index_shift += 2) { + appx2[len0] = appx1[index_shift]; + det2[len0] = det1[index_shift]; + len0++; + } + upsamp2(appx2, len0, U, tempx); + per_ext(tempx, 2 * len0, lf / 2, cL0); + + upsamp2(det2, len0, U, tempx); + per_ext(tempx, 2 * len0, lf / 2, cH0); + + N1 = 2 * len0 + lf; + + if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) { + wt->cobj = conv_init(N1, lf); + wt->cfftset = 1; + } + else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) { + printf("Decomposition Filters must have the same length."); + exit(-1); + } + + wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L); + + wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H); + + for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) { + oup00[i - lf + 1] = oup00L[i] + oup00H[i]; + } + + //SHIFT 1 + + len0 = 0; + + for (index_shift = 1; index_shift < len; index_shift += 2) { + appx2[len0] = appx1[index_shift]; + det2[len0] = det1[index_shift]; + len0++; + } + + upsamp2(appx2, len0, U, tempx); + per_ext(tempx, 2 * len0, lf / 2, cL0); + + upsamp2(det2, len0, U, tempx); + per_ext(tempx, 2 * len0, lf / 2, cH0); + + N1 = 2 * len0 + lf; + + wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L); + + wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H); + + for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) { + oup01[i - lf + 1] = oup00L[i] + oup00H[i]; + } + + circshift(oup01, 2*len0, -1); + + index2 = 0; + + for (index = count; index < N; index += value) { + swtop[index] = (oup00[index2] + oup01[index2]) / 2.0; + index2++; + } + + } + for (i = 0; i < N; ++i) { + appx_sig[i] = swtop[i]; + } + + } + + + + free(appx_sig); + free(det_sig); + free(appx1); + free(det1); + free(tempx); + free(cL0); + free(cH0); + free(oup00L); + free(oup00H); + free(oup00); + free(oup01); + free(appx2); + 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) { + int l, i, t, len_avg; + double s; + double *filt; + len_avg = wt->wave->lpd_len; + + filt = (double*)malloc(sizeof(double)* 2 * len_avg); + s = sqrt(2.0); + for (i = 0; i < len_avg; ++i) { + filt[i] = wt->wave->lpd[i] / s; + filt[len_avg + i] = wt->wave->hpd[i] / s; + } + + for (i = 0; i < len_cA; ++i) { + t = i; + cA[i] = filt[0] * inp[t]; + cD[i] = filt[len_avg] * inp[t]; + for (l = 1; l < len_avg; l++) { + t -= M; + while (t >= len_cA) { + t -= len_cA; + } + while (t < 0) { + t += len_cA; + } + + cA[i] += filt[l] * inp[t]; + cD[i] += filt[len_avg + l] * inp[t]; + + } + } + free(filt); +} + +void modwt(wt_object wt, double *inp) { + int i, J, temp_len, iter, M, N; + int lenacc, len_filt; + double *cA, *cD; + + temp_len = wt->siglength; + J = wt->J; + wt->length[0] = wt->length[J] = temp_len; + wt->outlength = wt->length[J + 1] = (J + 1) * temp_len; + len_filt = wt->wave->filtlength; + M = 1; + for (iter = 1; iter < J; ++iter) { + M = 2 * M; + wt->length[iter] = temp_len; + } + + + cA = (double*)malloc(sizeof(double)* temp_len); + cD = (double*)malloc(sizeof(double)* temp_len); + + M = 1; + + for (i = 0; i < temp_len; ++i) { + wt->params[i] = inp[i]; + } + + lenacc = wt->outlength; + + for (iter = 0; iter < J; ++iter) { + lenacc -= temp_len; + if (iter > 0) { + 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); + + + for (i = 0; i < temp_len; ++i) { + wt->params[i] = cA[i]; + wt->params[lenacc + i] = cD[i]; + } + + } + + free(cA); + free(cD); + +} + +static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, int len_cD, double *X) { + int len_avg, i, l, t, j; + double s; + double *filt; + len_avg = wt->wave->lpd_len; + + filt = (double*)malloc(sizeof(double)* 2 * len_avg); + s = sqrt(2.0); + for (i = 0; i < len_avg; ++i) { + filt[i] = wt->wave->lpd[i] / s; + filt[len_avg + i] = wt->wave->hpd[i] / s; + } + + + for (i = 0; i < len_cA; ++i) { + t = i; + X[i] = (filt[0] * cA[t]) + (filt[len_avg] * cD[t]); + for (l = 1; l < len_avg; l++) { + t += M; + while (t >= len_cA) { + t -= len_cA; + } + while (t < 0) { + t += len_cA; + } + + X[i] += (filt[l] * cA[t]) + (filt[len_avg + l] * cD[t]); + + } + } + free(filt); +} + +void imodwt(wt_object wt, double *dwtop) { + int N, lf, iter, i, J, index, value, j, len,U; + int lenacc,M; + double *X; + + N = wt->siglength; + J = wt->J; + U = 2; + lf = wt->wave->lpr_len; + lenacc = N; + M = (int)pow(2.0, (double)J - 1.0); + //M = 1; + X = (double*)malloc(sizeof(double)* N); + + for (i = 0; i < N; ++i) { + dwtop[i] = wt->output[i]; + } + + for (iter = 0; iter < J; ++iter) { + if (iter > 0) { + M = M / 2; + } + imodwt_per(wt, M, dwtop, N, wt->params + lenacc, N, X); + /* + for (j = lf - 1; j < N; ++j) { + dwtop[j - lf + 1] = X[j]; + } + for (j = 0; j < lf - 1; ++j) { + dwtop[N - lf + 1 + j] = X[j]; + } + */ + for (j = 0; j < N; ++j) { + dwtop[j] = X[j]; + } + + lenacc += N; + } + free(X); +} + +void setDWTExtension(wt_object wt, char *extension) { + if (!strcmp(extension, "sym")) { + strcpy(wt->ext, "sym"); + } + else if (!strcmp(extension, "per")) { + strcpy(wt->ext, "per"); + } + else { + printf("Signal extension can be either per or sym"); + exit(-1); + } +} + +void setWTConv(wt_object wt, char *cmethod) { + if (!strcmp(cmethod, "fft") || !strcmp(cmethod, "FFT")) { + strcpy(wt->cmethod, "fft"); + } + else if (!strcmp(cmethod, "direct")) { + strcpy(wt->cmethod, "direct"); + } + else { + printf("Convolution Only accepts two methods - direct and fft"); + exit(-1); + } +} + +void wave_summary(wave_object obj) { + int i,N; + N = obj->filtlength; + printf("\n"); + printf("Wavelet Name : %s \n",obj->wname); + printf("\n"); + printf("Wavelet Filters \n\n"); + printf("lpd : ["); + for (i = 0; i < N-1; ++i) { + printf("%g,", obj->lpd[i]); + } + printf("%g", obj->lpd[N-1]); + printf("] \n\n"); + printf("hpd : ["); + for (i = 0; i < N-1; ++i) { + printf("%g,", obj->hpd[i]); + } + printf("%g", obj->hpd[N - 1]); + printf("] \n\n"); + printf("lpr : ["); + for (i = 0; i < N-1; ++i) { + printf("%g,", obj->lpr[i]); + } + printf("%g", obj->lpr[N - 1]); + printf("] \n\n"); + printf("hpr : ["); + for (i = 0; i < N-1; ++i) { + printf("%g,", obj->hpr[i]); + } + printf("%g", obj->hpr[N - 1]); + printf("] \n\n"); +} + +void wt_summary(wt_object wt) { + int i; + int J,t; + J = wt->J; + wave_summary(wt->wave); + printf("\n"); + printf("Wavelet Transform : %s \n", wt->method); + printf("\n"); + printf("Signal Extension : %s \n", wt->ext); + printf("\n"); + printf("Convolutional Method : %s \n", wt->cmethod); + printf("\n"); + printf("Number of Decomposition Levels %d \n", wt->J); + printf("\n"); + printf("Length of Input Signal %d \n", wt->siglength); + printf("\n"); + printf("Length of WT Output Vector %d \n", wt->outlength); + printf("\n"); + printf("Wavelet Coefficients are contained in vector : %s \n", "output"); + printf("\n"); + printf("Approximation Coefficients \n"); + printf("Level %d Access : output[%d] Length : %d \n", 1, 0, wt->length[0]); + printf("\n"); + printf("Detail Coefficients \n"); + t = wt->length[0]; + for (i = 0; i < J; ++i) { + printf("Level %d Access : output[%d] Length : %d \n", i + 1,t,wt->length[i+1]); + t += wt->length[i+1]; + } + printf("\n"); + +} + +void wave_free(wave_object object) { + free(object); +} + +void wt_free(wt_object object) { + free(object); +} \ No newline at end of file diff --git a/src/wavelib.h b/src/wavelib.h new file mode 100644 index 0000000..8cd6a41 --- /dev/null +++ b/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_ */ + + diff --git a/src/wtmath.c b/src/wtmath.c new file mode 100644 index 0000000..42a81b5 --- /dev/null +++ b/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; +} \ No newline at end of file diff --git a/src/wtmath.h b/src/wtmath.h new file mode 100644 index 0000000..076c707 --- /dev/null +++ b/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_ */ \ No newline at end of file diff --git a/test/dwttest.c b/test/dwttest.c new file mode 100644 index 0000000..ce2a92e --- /dev/null +++ b/test/dwttest.c @@ -0,0 +1,83 @@ +#include +#include +#include +#include +#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; +} diff --git a/test/modwttest.c b/test/modwttest.c new file mode 100644 index 0000000..426a41c --- /dev/null +++ b/test/modwttest.c @@ -0,0 +1,86 @@ +#include +#include +#include +#include +#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; +} + diff --git a/test/signal.txt b/test/signal.txt new file mode 100644 index 0000000..9154ffd --- /dev/null +++ b/test/signal.txt @@ -0,0 +1 @@ + -18.3237 -18.2232 -18.0974 -17.9410 -17.7480 -17.5113 -17.2230 -16.8744 -16.4558 -15.9565 -15.3653 -14.6701 -13.8586 -12.9182 -11.8363 -10.6008 -9.2006 -7.6257 -5.8680 -3.9217 -1.7839 0.5452 3.0614 5.7562 8.6167 11.6252 14.7591 17.9909 21.2884 24.6155 27.9319 31.1947 34.3587 37.3775 40.2049 42.7957 13.2164 14.2125 15.0317 15.6595 16.0845 16.2990 16.2990 16.0845 15.6595 15.0317 14.2125 13.2164 12.0608 10.7654 9.3517 34.3587 31.1947 27.9319 24.6155 21.2884 17.9909 14.7591 11.6252 8.6167 5.7562 3.0614 0.5452 -1.7839 -3.9217 -5.8680 -7.6257 -9.2006 -10.6008 -11.8363 -12.9182 -13.8586 -14.6701 -15.3653 -15.9565 -16.4558 -16.8744 -17.2230 -17.5113 -17.7480 -17.9410 -18.0974 -18.2232 -18.3237 -18.4035 -18.0080 -17.8889 -17.7403 -17.5533 -17.3156 -17.0102 -16.6129 -16.0884 -15.3848 -14.4239 -13.0840 -11.1708 -8.3634 -4.1098 2.5833 13.6048 32.7934 28.0187 10.9660 1.0776 -4.9459 -8.7354 -11.1225 -12.4865 -12.8019 -11.2050 -3.3124 1.8995 -11.3573 -15.0684 -16.5028 -17.1937 -17.5831 -17.8288 -17.9968 -18.1185 -18.2103 -18.2818 -18.3388 -18.3849 -18.4229 -18.4545 -18.4810 -17.4642 -17.2104 -16.9033 -16.5317 -16.0822 -15.5384 -14.8804 -14.0844 -13.1214 -11.9563 -10.5467 -8.8414 -6.7782 -4.2822 -1.2624 2.3911 6.8111 12.1585 18.6280 26.4549 35.9241 35.9241 26.4549 18.6280 12.1585 6.8111 2.3911 -1.2624 -4.2822 -6.7782 -8.8414 -10.5467 -11.9563 -13.1214 -14.0844 -14.8804 -15.5384 -16.0822 -16.5317 -16.9033 -17.2104 -17.4642 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 6.3259 34.8066 34.6752 34.5285 34.3645 34.1812 33.9763 33.7474 33.4917 33.2058 32.8863 32.5294 32.1304 31.6846 31.1864 30.6296 30.0074 29.3121 28.5350 27.6667 26.6963 25.6118 24.3999 23.0456 21.5322 19.8408 17.9507 15.8385 13.4781 10.8403 7.8925 4.5982 0.9168 -3.1972 -7.7947 -12.9325 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.6741 -18.3237 \ No newline at end of file diff --git a/test/swttest.c b/test/swttest.c new file mode 100644 index 0000000..b6da495 --- /dev/null +++ b/test/swttest.c @@ -0,0 +1,87 @@ +#include +#include +#include +#include +#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; +}