1 /*
2 Copyright (c) 2014, Rafat Hussain
3 */
4
5 #include <math.h>
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <string.h>
9
10 #include "cwt.h"
11 #include "../header/wavelib.h"
12 #include "wtmath.h"
13
wave_init(const char * wname)14 wave_object wave_init(const char* wname) {
15 wave_object obj = NULL;
16 int retval;
17 retval = 0;
18
19 if (wname != NULL) {
20 retval = filtlength(wname);
21 //obj->filtlength = retval;
22 //strcopy(obj->wname, wname);
23 }
24
25 obj = (wave_object)malloc(sizeof(struct wave_set) + sizeof(double)* 4 * retval);
26
27 obj->filtlength = retval;
28 obj->lpd_len = obj->hpd_len = obj->lpr_len = obj->hpr_len = obj->filtlength;
29 strcpy(obj->wname, wname);
30 if (wname != NULL) {
31 filtcoef(wname,obj->params,obj->params+retval,obj->params+2*retval,obj->params+3*retval);
32 }
33 obj->lpd = &obj->params[0];
34 obj->hpd = &obj->params[retval];
35 obj->lpr = &obj->params[2 * retval];
36 obj->hpr = &obj->params[3 * retval];
37 return obj;
38 }
39
wt_init(wave_object wave,const char * method,int siglength,int J)40 wt_object wt_init(wave_object wave,const char* method, int siglength,int J) {
41 int size,i,MaxIter;
42 wt_object obj = NULL;
43
44 size = wave->filtlength;
45
46 if (J > 100) {
47 printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
48 exit(-1);
49 }
50
51 MaxIter = wmaxiter(siglength, size);
52
53 if (J > MaxIter) {
54 printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n",MaxIter);
55 exit(-1);
56 }
57
58 if (method == NULL) {
59 obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength + 2 * J * (size+1)));
60 obj->outlength = siglength + 2 * J * (size + 1); // Default
61 strcpy(obj->ext, "sym"); // Default
62 }
63 else if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
64 obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength + 2 * J * (size + 1)));
65 obj->outlength = siglength + 2 * J * (size + 1); // Default
66 strcpy(obj->ext, "sym"); // Default
67 }
68 else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) {
69 if (!testSWTlength(siglength, J)) {
70 printf("\n For SWT the signal length must be a multiple of 2^J. \n");
71 exit(-1);
72 }
73
74 obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * (J + 1)));
75 obj->outlength = siglength * (J + 1); // Default
76 strcpy(obj->ext, "per"); // Default
77 }
78 else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
79
80 if (!strstr(wave->wname,"haar")) {
81 if (!strstr(wave->wname,"db")) {
82 if (!strstr(wave->wname, "sym")) {
83 if (!strstr(wave->wname, "coif")) {
84 printf("\n MODWT is only implemented for orthogonal wavelet families - db, sym and coif \n");
85 exit(-1);
86 }
87 }
88 }
89 }
90
91 obj = (wt_object)malloc(sizeof(struct wt_set) + sizeof(double)* (siglength * 2 * (J + 1)));
92 obj->outlength = siglength * (J + 1); // Default
93 strcpy(obj->ext, "per"); // Default
94 }
95
96 obj->wave = wave;
97 obj->siglength = siglength;
98 obj->modwtsiglength = siglength;
99 obj->J = J;
100 obj->MaxIter = MaxIter;
101 strcpy(obj->method, method);
102
103 if (siglength % 2 == 0) {
104 obj->even = 1;
105 }
106 else {
107 obj->even = 0;
108 }
109
110 obj->cobj = NULL;
111
112 strcpy(obj->cmethod, "direct"); // Default
113 obj->cfftset = 0;
114 obj->lenlength = J + 2;
115 obj->output = &obj->params[0];
116 if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
117 for (i = 0; i < siglength + 2 * J * (size + 1); ++i) {
118 obj->params[i] = 0.0;
119 }
120 }
121 else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) {
122 for (i = 0; i < siglength * (J + 1); ++i) {
123 obj->params[i] = 0.0;
124 }
125 }
126 else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
127 for (i = 0; i < siglength * 2 * (J + 1); ++i) {
128 obj->params[i] = 0.0;
129 }
130 }
131 //wave_summary(obj->wave);
132
133 return obj;
134 }
135
wtree_init(wave_object wave,int siglength,int J)136 wtree_object wtree_init(wave_object wave, int siglength,int J) {
137 int size,i,MaxIter,temp,temp2,elength,nodes;
138 wtree_object obj = NULL;
139
140 size = wave->filtlength;
141
142 if (J > 100) {
143 printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
144 exit(-1);
145 }
146
147
148 MaxIter = wmaxiter(siglength, size);
149 if (J > MaxIter) {
150 printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
151 exit(-1);
152 }
153 temp = 1;
154 elength = 0;
155 nodes = 0;
156 for(i = 0; i < J;++i) {
157 temp *= 2;
158 nodes += temp;
159 temp2 = (size - 2) * (temp - 1);
160 elength += temp2;
161 }
162
163 obj = (wtree_object)malloc(sizeof(struct wtree_set) + sizeof(double)* (siglength * (J + 1) + elength + nodes + J + 1));
164 obj->outlength = siglength * (J + 1) + elength;
165 strcpy(obj->ext, "sym");
166
167 obj->wave = wave;
168 obj->siglength = siglength;
169 obj->J = J;
170 obj->MaxIter = MaxIter;
171 strcpy(obj->method, "dwt");
172
173 if (siglength % 2 == 0) {
174 obj->even = 1;
175 }
176 else {
177 obj->even = 0;
178 }
179
180 obj->cobj = NULL;
181 obj->nodes = nodes;
182
183 obj->cfftset = 0;
184 obj->lenlength = J + 2;
185 obj->output = &obj->params[0];
186 obj->nodelength = (int*) &obj->params[siglength * (J + 1) + elength];
187 obj->coeflength = (int*)&obj->params[siglength * (J + 1) + elength + nodes];
188
189 for (i = 0; i < siglength * (J + 1) + elength + nodes + J + 1; ++i) {
190 obj->params[i] = 0.0;
191 }
192
193 //wave_summary(obj->wave);
194
195 return obj;
196 }
197
wpt_init(wave_object wave,int siglength,int J)198 wpt_object wpt_init(wave_object wave, int siglength, int J) {
199 int size, i, MaxIter, temp, nodes,elength,p2,N,lp;
200 wpt_object obj = NULL;
201
202 size = wave->filtlength;
203
204 if (J > 100) {
205 printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
206 exit(-1);
207 }
208
209
210 MaxIter = wmaxiter(siglength, size);
211 if (J > MaxIter) {
212 printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
213 exit(-1);
214 }
215 temp = 1;
216 nodes = 0;
217 for (i = 0; i < J; ++i) {
218 temp *= 2;
219 nodes += temp;
220 }
221
222 i = J;
223 p2 = 2;
224 N = siglength;
225 lp = size;
226 elength = 0;
227 while (i > 0) {
228 N = N + lp - 2;
229 N = (int)ceil((double)N / 2.0);
230 elength = p2 * N;
231 i--;
232 p2 *= 2;
233 }
234 //printf("elength %d", elength);
235
236 obj = (wpt_object)malloc(sizeof(struct wpt_set) + sizeof(double)* (elength + 4 * nodes + 2 * J + 6));
237 obj->outlength = siglength + 2 * (J + 1) * (size + 1);
238 strcpy(obj->ext, "sym");
239 strcpy(obj->entropy, "shannon");
240 obj->eparam = 0.0;
241
242 obj->wave = wave;
243 obj->siglength = siglength;
244 obj->J = J;
245 obj->MaxIter = MaxIter;
246
247 if (siglength % 2 == 0) {
248 obj->even = 1;
249 }
250 else {
251 obj->even = 0;
252 }
253
254 obj->cobj = NULL;
255 obj->nodes = nodes;
256
257 obj->lenlength = J + 2;
258 obj->output = &obj->params[0];
259 obj->costvalues = &obj->params[elength];
260 obj->basisvector = &obj->params[elength + nodes + 1];
261 obj->nodeindex = (int*)&obj->params[elength + 2*nodes + 2];
262 obj->numnodeslevel = (int*)&obj->params[elength + 4 * nodes + 4];
263 obj->coeflength = (int*)&obj->params[elength + 4 * nodes + J + 5];
264
265 for (i = 0; i < elength + 4 * nodes + 2 * J + 6; ++i) {
266 obj->params[i] = 0.0;
267 }
268
269 //wave_summary(obj->wave);
270
271 return obj;
272 }
273
cwt_init(const char * wave,double param,int siglength,double dt,int J)274 cwt_object cwt_init(const char* wave, double param,int siglength, double dt, int J) {
275 cwt_object obj = NULL;
276 int N, i,nj2,ibase2,mother;
277 double s0, dj;
278 double t1;
279 int m, odd;
280 const char *pdefault = "pow";
281
282 m = (int)param;
283 odd = 1;
284 if (2 * (m / 2) == m) {
285 odd = 0;
286 }
287
288 N = siglength;
289 nj2 = 2 * N * J;
290 obj = (cwt_object)malloc(sizeof(struct cwt_set) + sizeof(double)* (nj2 + 2 * J + N));
291
292 if (!strcmp(wave, "morlet") || !strcmp(wave, "morl")) {
293 s0 = 2 * dt;
294 dj = 0.4875;
295 mother = 0;
296 if (param < 0.0) {
297 printf("\n Morlet Wavelet Parameter should be >= 0 \n");
298 exit(-1);
299 }
300 if (param == 0) {
301 param = 6.0;
302 }
303 strcpy(obj->wave,"morlet");
304
305 }
306 else if (!strcmp(wave, "paul")) {
307 s0 = 2 * dt;
308 dj = 0.4875;
309 mother = 1;
310 if (param < 0 || param > 20) {
311 printf("\n Paul Wavelet Parameter should be > 0 and <= 20 \n");
312 exit(-1);
313 }
314 if (param == 0) {
315 param = 4.0;
316 }
317 strcpy(obj->wave,"paul");
318
319 }
320 else if (!strcmp(wave, "dgauss") || !strcmp(wave, "dog")) {
321 s0 = 2 * dt;
322 dj = 0.4875;
323 mother = 2;
324 if (param < 0 || odd == 1) {
325 printf("\n DOG Wavelet Parameter should be > 0 and even \n");
326 exit(-1);
327 }
328 if (param == 0) {
329 param = 2.0;
330 }
331 strcpy(obj->wave,"dog");
332 }
333
334 obj->pow = 2;
335 strcpy(obj->type, pdefault);
336
337 obj->s0 = s0;
338 obj->dj = dj;
339 obj->dt = dt;
340 obj->J = J;
341 obj->siglength = siglength;
342 obj->sflag = 0;
343 obj->pflag = 1;
344 obj->mother = mother;
345 obj->m = param;
346
347 t1 = 0.499999 + log((double)N) / log(2.0);
348 ibase2 = 1 + (int)t1;
349
350 obj->npad = (int)pow(2.0, (double)ibase2);
351
352 obj->output = (cplx_data*) &obj->params[0];
353 obj->scale = &obj->params[nj2];
354 obj->period = &obj->params[nj2 + J];
355 obj->coi = &obj->params[nj2 + 2*J];
356
357 for (i = 0; i < nj2 + 2 * J + N; ++i) {
358 obj->params[i] = 0.0;
359 }
360
361 return obj;
362 }
363
wt2_init(wave_object wave,const char * method,int rows,int cols,int J)364 wt2_object wt2_init(wave_object wave, const char* method, int rows, int cols, int J) {
365 int size, i, MaxIter, MaxRows, MaxCols,sumacc;
366 wt2_object obj = NULL;
367
368 size = wave->filtlength;
369
370 MaxRows = wmaxiter(rows, size);
371 MaxCols = wmaxiter(cols, size);
372
373 MaxIter = (MaxRows < MaxCols) ? MaxRows : MaxCols;
374
375 if (J > MaxIter) {
376 printf("\n Error - The Signal Can only be iterated %d times using this wavelet. Exiting\n", MaxIter);
377 exit(-1);
378 }
379
380 if (J == 1) {
381 sumacc = 4;
382 }
383 else if (J > 1) {
384 sumacc = J * 3 + 1;
385 }
386 else {
387 printf("Error : J should be >= 1 \n");
388 exit(-1);
389 }
390
391 if (method == NULL) {
392 obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
393 obj->outlength = 0; // Default
394 strcpy(obj->ext, "per");
395 }
396 else if (!strcmp(method, "dwt") || !strcmp(method, "DWT")) {
397 obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
398 obj->outlength = 0; // Default
399 strcpy(obj->ext, "per");
400 }
401 else if (!strcmp(method, "swt") || !strcmp(method, "SWT")) {
402 if (!testSWTlength(rows, J) || !testSWTlength(cols, J)) {
403 printf("\n For SWT data rows and columns must be a multiple of 2^J. \n");
404 exit(-1);
405 }
406
407 obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
408 obj->outlength = 0; // Default
409 strcpy(obj->ext, "per");
410 }
411 else if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
412 if (!strstr(wave->wname, "haar")) {
413 if (!strstr(wave->wname, "db")) {
414 if (!strstr(wave->wname, "sym")) {
415 if (!strstr(wave->wname, "coif")) {
416 printf("\n MODWT is only implemented for orthogonal wavelet families - db, sym and coif \n");
417 exit(-1);
418 }
419 }
420 }
421 }
422 obj = (wt2_object)malloc(sizeof(struct wt2_set) + sizeof(int)* (2 * J + sumacc));
423 obj->outlength = 0; // Default
424 strcpy(obj->ext, "per");
425 }
426
427
428 obj->wave = wave;
429 obj->rows = rows;
430 obj->cols = cols;
431 obj->J = J;
432 obj->MaxIter = MaxIter;
433 strcpy(obj->method, method);
434 obj->coeffaccesslength = sumacc;
435
436 obj->dimensions = &obj->params[0];
437 obj->coeffaccess = &obj->params[2 * J];
438 for (i = 0; i < (2 * J + sumacc); ++i) {
439 obj->params[i] = 0;
440 }
441
442 return obj;
443 }
444
445
wconv(wt_object wt,double * sig,int N,double * filt,int L,double * oup)446 static void wconv(wt_object wt, double *sig, int N, double *filt, int L, double *oup) {
447 if (!strcmp(wt->cmethod,"direct")) {
448 conv_direct(sig, N, filt, L, oup);
449 }
450 else if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) {
451 if (wt->cfftset == 0) {
452 wt->cobj = conv_init(N, L);
453 conv_fft(wt->cobj, sig, filt, oup);
454 free_conv(wt->cobj);
455 }
456 else {
457 conv_fft(wt->cobj, sig, filt, oup);
458 }
459 }
460 else {
461 printf("Convolution Only accepts two methods - direct and fft");
462 exit(-1);
463 }
464 }
465
466
dwt_per(wt_object wt,double * inp,int N,double * cA,int len_cA,double * cD)467 static void dwt_per(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
468
469 dwt_per_stride(inp,N,wt->wave->lpd,wt->wave->hpd,wt->wave->lpd_len,cA,len_cA,cD,1,1);
470
471 }
472
wtree_per(wtree_object wt,double * inp,int N,double * cA,int len_cA,double * cD)473 static void wtree_per(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
474 int l, l2, isodd, i, t, len_avg;
475
476 len_avg = wt->wave->lpd_len;
477 l2 = len_avg / 2;
478 isodd = N % 2;
479
480 for (i = 0; i < len_cA; ++i) {
481 t = 2 * i + l2;
482 cA[i] = 0.0;
483 cD[i] = 0.0;
484 for (l = 0; l < len_avg; ++l) {
485 if ((t - l) >= l2 && (t - l) < N) {
486 cA[i] += wt->wave->lpd[l] * inp[t - l];
487 cD[i] += wt->wave->hpd[l] * inp[t - l];
488 }
489 else if ((t - l) < l2 && (t - l) >= 0) {
490 cA[i] += wt->wave->lpd[l] * inp[t - l];
491 cD[i] += wt->wave->hpd[l] * inp[t - l];
492 }
493 else if ((t - l) < 0 && isodd == 0) {
494 cA[i] += wt->wave->lpd[l] * inp[t - l + N];
495 cD[i] += wt->wave->hpd[l] * inp[t - l + N];
496 }
497 else if ((t - l) < 0 && isodd == 1) {
498 if ((t - l) != -1) {
499 cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
500 cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
501 }
502 else {
503 cA[i] += wt->wave->lpd[l] * inp[N - 1];
504 cD[i] += wt->wave->hpd[l] * inp[N - 1];
505 }
506 }
507 else if ((t - l) >= N && isodd == 0) {
508 cA[i] += wt->wave->lpd[l] * inp[t - l - N];
509 cD[i] += wt->wave->hpd[l] * inp[t - l - N];
510 }
511 else if ((t - l) >= N && isodd == 1) {
512 if (t - l != N) {
513 cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
514 cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
515 }
516 else {
517 cA[i] += wt->wave->lpd[l] * inp[N - 1];
518 cD[i] += wt->wave->hpd[l] * inp[N - 1];
519 }
520 }
521
522 }
523 }
524
525
526
527 }
528
dwpt_per(wpt_object wt,double * inp,int N,double * cA,int len_cA,double * cD)529 static void dwpt_per(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
530 int l, l2, isodd, i, t, len_avg;
531
532 len_avg = wt->wave->lpd_len;
533 l2 = len_avg / 2;
534 isodd = N % 2;
535
536 for (i = 0; i < len_cA; ++i) {
537 t = 2 * i + l2;
538 cA[i] = 0.0;
539 cD[i] = 0.0;
540 for (l = 0; l < len_avg; ++l) {
541 if ((t - l) >= l2 && (t - l) < N) {
542 cA[i] += wt->wave->lpd[l] * inp[t - l];
543 cD[i] += wt->wave->hpd[l] * inp[t - l];
544 }
545 else if ((t - l) < l2 && (t - l) >= 0) {
546 cA[i] += wt->wave->lpd[l] * inp[t - l];
547 cD[i] += wt->wave->hpd[l] * inp[t - l];
548 }
549 else if ((t - l) < 0 && isodd == 0) {
550 cA[i] += wt->wave->lpd[l] * inp[t - l + N];
551 cD[i] += wt->wave->hpd[l] * inp[t - l + N];
552 }
553 else if ((t - l) < 0 && isodd == 1) {
554 if ((t - l) != -1) {
555 cA[i] += wt->wave->lpd[l] * inp[t - l + N + 1];
556 cD[i] += wt->wave->hpd[l] * inp[t - l + N + 1];
557 }
558 else {
559 cA[i] += wt->wave->lpd[l] * inp[N - 1];
560 cD[i] += wt->wave->hpd[l] * inp[N - 1];
561 }
562 }
563 else if ((t - l) >= N && isodd == 0) {
564 cA[i] += wt->wave->lpd[l] * inp[t - l - N];
565 cD[i] += wt->wave->hpd[l] * inp[t - l - N];
566 }
567 else if ((t - l) >= N && isodd == 1) {
568 if (t - l != N) {
569 cA[i] += wt->wave->lpd[l] * inp[t - l - (N + 1)];
570 cD[i] += wt->wave->hpd[l] * inp[t - l - (N + 1)];
571 }
572 else {
573 cA[i] += wt->wave->lpd[l] * inp[N - 1];
574 cD[i] += wt->wave->hpd[l] * inp[N - 1];
575 }
576 }
577
578 }
579 }
580
581
582
583 }
584
dwt_sym(wt_object wt,double * inp,int N,double * cA,int len_cA,double * cD)585 static void dwt_sym(wt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
586
587 dwt_sym_stride(inp,N,wt->wave->lpd,wt->wave->hpd,wt->wave->lpd_len,cA,len_cA,cD,1,1);
588 }
589
wtree_sym(wtree_object wt,double * inp,int N,double * cA,int len_cA,double * cD)590 static void wtree_sym(wtree_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
591 int i, l, t, len_avg;
592
593 len_avg = wt->wave->lpd_len;
594
595 for (i = 0; i < len_cA; ++i) {
596 t = 2 * i + 1;
597 cA[i] = 0.0;
598 cD[i] = 0.0;
599 for (l = 0; l < len_avg; ++l) {
600 if ((t - l) >= 0 && (t - l) < N) {
601 cA[i] += wt->wave->lpd[l] * inp[t - l];
602 cD[i] += wt->wave->hpd[l] * inp[t - l];
603 }
604 else if ((t - l) < 0) {
605 cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
606 cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
607 }
608 else if ((t - l) >= N) {
609 cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
610 cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
611 }
612
613 }
614 }
615
616
617 }
618
dwpt_sym(wpt_object wt,double * inp,int N,double * cA,int len_cA,double * cD)619 static void dwpt_sym(wpt_object wt, double *inp, int N, double *cA, int len_cA, double *cD) {
620 int i, l, t, len_avg;
621
622 len_avg = wt->wave->lpd_len;
623
624 for (i = 0; i < len_cA; ++i) {
625 t = 2 * i + 1;
626 cA[i] = 0.0;
627 cD[i] = 0.0;
628 for (l = 0; l < len_avg; ++l) {
629 if ((t - l) >= 0 && (t - l) < N) {
630 cA[i] += wt->wave->lpd[l] * inp[t - l];
631 cD[i] += wt->wave->hpd[l] * inp[t - l];
632 }
633 else if ((t - l) < 0) {
634 cA[i] += wt->wave->lpd[l] * inp[-t + l - 1];
635 cD[i] += wt->wave->hpd[l] * inp[-t + l - 1];
636 }
637 else if ((t - l) >= N) {
638 cA[i] += wt->wave->lpd[l] * inp[2 * N - t + l - 1];
639 cD[i] += wt->wave->hpd[l] * inp[2 * N - t + l - 1];
640 }
641
642 }
643 }
644
645
646 }
647
dwt1(wt_object wt,double * sig,int len_sig,double * cA,double * cD)648 static void dwt1(wt_object wt,double *sig,int len_sig, double *cA, double *cD) {
649 int len_avg,D,lf;
650 double *signal,*cA_undec;
651 len_avg = (wt->wave->lpd_len + wt->wave->hpd_len) / 2;
652 //len_sig = 2 * (int)ceil((double)len_sig / 2.0);
653
654 D = 2;
655
656 if (!strcmp(wt->ext, "per")) {
657 signal = (double*)malloc(sizeof(double)* (len_sig + len_avg + (len_sig % 2)));
658
659 len_sig = per_ext(sig, len_sig, len_avg / 2, signal);
660
661 cA_undec = (double*)malloc(sizeof(double)* (len_sig + len_avg + wt->wave->lpd_len - 1));
662
663 if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
664 wt->cobj = conv_init(len_sig + len_avg, wt->wave->lpd_len);
665 wt->cfftset = 1;
666 }
667 else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
668 printf("Decomposition Filters must have the same length.");
669 exit(-1);
670 }
671
672 wconv(wt, signal, len_sig + len_avg, wt->wave->lpd, wt->wave->lpd_len, cA_undec);
673
674 downsamp(cA_undec + len_avg, len_sig, D, cA);
675
676 wconv(wt, signal, len_sig + len_avg, wt->wave->hpd, wt->wave->hpd_len, cA_undec);
677
678 downsamp(cA_undec + len_avg, len_sig, D, cD);
679 }
680 else if (!strcmp(wt->ext, "sym")) {
681 //printf("\n YES %s \n", wt->ext);
682 lf = wt->wave->lpd_len;// lpd and hpd have the same length
683
684 signal = (double*)malloc(sizeof(double)* (len_sig + 2 * (lf - 1)));
685
686 len_sig = symm_ext(sig, len_sig, lf - 1, signal);
687
688 cA_undec = (double*)malloc(sizeof(double)* (len_sig + 3 * (lf - 1)));
689
690 if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
691 wt->cobj = conv_init(len_sig + 2 * (lf - 1), lf);
692 wt->cfftset = 1;
693 }
694 else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
695 printf("Decomposition Filters must have the same length.");
696 exit(-1);
697 }
698
699
700 wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->lpd, wt->wave->lpd_len, cA_undec);
701
702 downsamp(cA_undec + lf, len_sig + lf - 2, D, cA);
703
704 wconv(wt, signal, len_sig + 2 * (lf - 1), wt->wave->hpd, wt->wave->hpd_len, cA_undec);
705
706 downsamp(cA_undec + lf, len_sig + lf - 2, D, cD);
707 }
708 else {
709 printf("Signal extension can be either per or sym");
710 exit(-1);
711 }
712
713
714 if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
715 free_conv(wt->cobj);
716 wt->cfftset = 0;
717 }
718
719 free(signal);
720 free(cA_undec);
721 }
722
dwt(wt_object wt,const double * inp)723 void dwt(wt_object wt,const double *inp) {
724 int i,J,temp_len,iter,N,lp;
725 int len_cA;
726 double *orig,*orig2;
727
728 temp_len = wt->siglength;
729 J = wt->J;
730 wt->length[J + 1] = temp_len;
731 wt->outlength = 0;
732 wt->zpad = 0;
733 orig = (double*)malloc(sizeof(double)* temp_len);
734 orig2 = (double*)malloc(sizeof(double)* temp_len);
735 /*
736 if ((temp_len % 2) == 0) {
737 wt->zpad = 0;
738 orig = (double*)malloc(sizeof(double)* temp_len);
739 orig2 = (double*)malloc(sizeof(double)* temp_len);
740 }
741 else {
742 wt->zpad = 1;
743 temp_len++;
744 orig = (double*)malloc(sizeof(double)* temp_len);
745 orig2 = (double*)malloc(sizeof(double)* temp_len);
746 }
747 */
748
749 for (i = 0; i < wt->siglength; ++i) {
750 orig[i] = inp[i];
751 }
752
753 if (wt->zpad == 1) {
754 orig[temp_len - 1] = orig[temp_len - 2];
755 }
756
757 N = temp_len;
758 lp = wt->wave->lpd_len;
759
760 if (!strcmp(wt->ext,"per")) {
761 i = J;
762 while (i > 0) {
763 N = (int)ceil((double)N / 2.0);
764 wt->length[i] = N;
765 wt->outlength += wt->length[i];
766 i--;
767 }
768 wt->length[0] = wt->length[1];
769 wt->outlength += wt->length[0];
770 N = wt->outlength;
771
772 for (iter = 0; iter < J; ++iter) {
773 len_cA = wt->length[J - iter];
774 N -= len_cA;
775 if ( !strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT") ) {
776 dwt1(wt, orig, temp_len, orig2, wt->params + N);
777 }
778 else {
779 dwt_per(wt, orig, temp_len, orig2, len_cA, wt->params + N);
780 }
781 temp_len = wt->length[J - iter];
782 if (iter == J - 1) {
783 for (i = 0; i < len_cA; ++i) {
784 wt->params[i] = orig2[i];
785 }
786 }
787 else {
788 for (i = 0; i < len_cA; ++i) {
789 orig[i] = orig2[i];
790 }
791 }
792 }
793 }
794 else if (!strcmp(wt->ext,"sym")) {
795 //printf("\n YES %s \n", wt->ext);
796 i = J;
797 while (i > 0) {
798 N = N + lp - 2;
799 N = (int) ceil((double)N / 2.0);
800 wt->length[i] = N;
801 wt->outlength += wt->length[i];
802 i--;
803 }
804 wt->length[0] = wt->length[1];
805 wt->outlength += wt->length[0];
806 N = wt->outlength;
807
808 for (iter = 0; iter < J; ++iter) {
809 len_cA = wt->length[J - iter];
810 N -= len_cA;
811 if (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT")) {
812 dwt1(wt, orig, temp_len, orig2, wt->params + N);
813 }
814 else {
815 dwt_sym(wt, orig, temp_len, orig2, len_cA, wt->params + N);
816 }
817 temp_len = wt->length[J - iter];
818
819 if (iter == J - 1) {
820 for (i = 0; i < len_cA; ++i) {
821 wt->params[i] = orig2[i];
822 }
823 }
824 else {
825 for (i = 0; i < len_cA; ++i) {
826 orig[i] = orig2[i];
827 }
828 }
829 }
830 }
831 else {
832 printf("Signal extension can be either per or sym");
833 exit(-1);
834 }
835
836 free(orig);
837 free(orig2);
838 }
839
getDWTRecCoeff(double * coeff,int * length,const char * ctype,const char * ext,int level,int J,double * lpr,double * hpr,int lf,int siglength,double * reccoeff)840 static void getDWTRecCoeff(double *coeff, int *length, const char *ctype, const char *ext, int level, int J, double *lpr,
841 double *hpr, int lf, int siglength, double *reccoeff) {
842
843 int i, j, k, det_len, N, l, m, n, v, t, l2;
844 double *out, *X_lp, *filt;
845 out = (double*)malloc(sizeof(double)* (siglength + 1));
846 l2 = lf / 2;
847 m = -2;
848 n = -1;
849 if (!strcmp(ext, "per")) {
850 if (!strcmp((ctype), "appx")) {
851 det_len = length[0];
852 }
853 else {
854 det_len = length[J - level + 1];
855 }
856
857 N = 2 * length[J];
858
859 X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
860
861 for (i = 0; i < det_len; ++i) {
862 out[i] = coeff[i];
863 }
864
865 for (j = level; j > 0; --j) {
866
867 if (!strcmp((ctype), "det") && j == level) {
868 filt = hpr;
869 }
870 else {
871 filt = lpr;
872 }
873
874 m = -2;
875 n = -1;
876
877 for (i = 0; i < det_len + l2 - 1; ++i) {
878 m += 2;
879 n += 2;
880 X_lp[m] = 0.0;
881 X_lp[n] = 0.0;
882 for (l = 0; l < l2; ++l) {
883 t = 2 * l;
884 if ((i - l) >= 0 && (i - l) < det_len) {
885 X_lp[m] += filt[t] * out[i - l];
886 X_lp[n] += filt[t + 1] * out[i - l];
887 }
888 else if ((i - l) >= det_len && (i - l) < det_len + lf - 1) {
889 X_lp[m] += filt[t] * out[i - l - det_len];
890 X_lp[n] += filt[t + 1] * out[i - l - det_len];
891 }
892 else if ((i - l) < 0 && (i - l) > -l2) {
893 X_lp[m] += filt[t] * out[det_len + i - l];
894 X_lp[n] += filt[t + 1] * out[det_len + i - l];
895 }
896 }
897 }
898
899 for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
900 out[k - lf / 2 + 1] = X_lp[k];
901 }
902
903 if (j != 1) {
904 det_len = length[J - j + 2];
905 }
906 }
907
908 free(X_lp);
909
910 }
911 else if (!strcmp(ext, "sym")) {
912 if (!strcmp((ctype), "appx")) {
913 det_len = length[0];
914 }
915 else {
916 det_len = length[J - level + 1];
917 }
918
919 N = 2 * length[J] - 1;
920
921 X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
922
923 for (i = 0; i < det_len; ++i) {
924 out[i] = coeff[i];
925 }
926
927 for (j = level; j > 0; --j) {
928
929 if (!strcmp((ctype), "det") && j == level) {
930 filt = hpr;
931 }
932 else {
933 filt = lpr;
934 }
935
936 m = -2;
937 n = -1;
938
939 for (v = 0; v < det_len; ++v) {
940 i = v;
941 m += 2;
942 n += 2;
943 X_lp[m] = 0.0;
944 X_lp[n] = 0.0;
945 for (l = 0; l < lf / 2; ++l) {
946 t = 2 * l;
947 if ((i - l) >= 0 && (i - l) < det_len) {
948 X_lp[m] += filt[t] * out[i - l];
949 X_lp[n] += filt[t + 1] * out[i - l];
950 }
951 }
952 }
953
954 for (k = lf - 2; k < 2 * det_len; ++k) {
955 out[k - lf + 2] = X_lp[k];
956 }
957
958
959 if (j != 1) {
960 det_len = length[J - j + 2];
961 }
962 }
963
964 free(X_lp);
965
966 }
967 else {
968 printf("Signal extension can be either per or sym");
969 exit(-1);
970 }
971
972 for (i = 0; i < siglength; ++i) {
973 reccoeff[i] = out[i];
974 }
975
976 free(out);
977
978 }
979
980
getDWTmra(wt_object wt,double * wavecoeffs)981 double *getDWTmra(wt_object wt, double *wavecoeffs) {
982 int i, J, access,N;
983 double *mra;
984 J = wt->J;
985 mra = (double*)malloc(sizeof(double)* wt->siglength*(J + 1));
986 access = 0;
987
988
989 // Approximation MRA
990 getDWTRecCoeff(wt->output + access, wt->length, "appx", wt->ext, J, J, wt->wave->lpr, wt->wave->hpr, wt->wave->lpr_len, wt->siglength, mra);
991
992 // Details MRA
993 N = wt->siglength;
994 for (i = J; i > 0; --i) {
995 access += wt->length[J - i];
996 getDWTRecCoeff(wt->output + access, wt->length, "det", wt->ext, i, J, wt->wave->lpr, wt->wave->hpr, wt->wave->lpr_len, wt->siglength, mra+N);
997 N += wt->siglength;
998 }
999
1000 return mra;
1001 }
1002
wtree(wtree_object wt,const double * inp)1003 void wtree(wtree_object wt,const double *inp) {
1004 int i,J,temp_len,iter,N,lp,p2,k,N2,Np;
1005 int len_cA,t,t2,it1;
1006 double *orig;
1007
1008 temp_len = wt->siglength;
1009 J = wt->J;
1010 wt->length[J + 1] = temp_len;
1011 wt->outlength = 0;
1012 wt->zpad = 0;
1013 orig = (double*)malloc(sizeof(double)* temp_len);
1014 /*
1015 if ((temp_len % 2) == 0) {
1016 wt->zpad = 0;
1017 orig = (double*)malloc(sizeof(double)* temp_len);
1018 }
1019 else {
1020 wt->zpad = 1;
1021 temp_len++;
1022 orig = (double*)malloc(sizeof(double)* temp_len);
1023 }
1024 */
1025 for (i = 0; i < wt->siglength; ++i) {
1026 orig[i] = inp[i];
1027 }
1028
1029 if (wt->zpad == 1) {
1030 orig[temp_len - 1] = orig[temp_len - 2];
1031 }
1032
1033 N = temp_len;
1034 lp = wt->wave->lpd_len;
1035 p2 = 1;
1036
1037 if (!strcmp(wt->ext,"per")) {
1038 i = J;
1039 p2 = 2;
1040 while (i > 0) {
1041 N = (int)ceil((double)N / 2.0);
1042 wt->length[i] = N;
1043 wt->outlength += p2 * (wt->length[i]);
1044 i--;
1045 p2 *= 2;
1046 }
1047 wt->length[0] = wt->length[1];
1048
1049 N2 = N = wt->outlength;
1050 p2 = 1;
1051 for (iter = 0; iter < J; ++iter) {
1052 len_cA = wt->length[J - iter];
1053 N2 -= 2 * p2 * len_cA;
1054 N = N2;
1055 for(k = 0; k < p2;++k) {
1056 if (iter == 0) {
1057 wtree_per(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
1058 } else {
1059 wtree_per(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
1060 }
1061 N += 2 * len_cA;
1062 }
1063
1064 temp_len = wt->length[J - iter];
1065 p2 = 2 * p2;
1066 Np = N2;
1067 }
1068 }
1069 else if (!strcmp(wt->ext,"sym")) {
1070 //printf("\n YES %s \n", wt->ext);
1071 i = J;
1072 p2 = 2;
1073 while (i > 0) {
1074 N = N + lp - 2;
1075 N = (int) ceil((double)N / 2.0);
1076 wt->length[i] = N;
1077 wt->outlength += p2 * (wt->length[i]);
1078 i--;
1079 p2 *= 2;
1080 }
1081 wt->length[0] = wt->length[1];
1082
1083 N2 = N = wt->outlength;
1084 p2 = 1;
1085
1086 for (iter = 0; iter < J; ++iter) {
1087 len_cA = wt->length[J - iter];
1088 N2 -= 2 * p2 * len_cA;
1089 N = N2;
1090 for(k = 0; k < p2;++k) {
1091 if (iter == 0) {
1092 wtree_sym(wt, orig, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
1093 } else {
1094 wtree_sym(wt, wt->params + Np + k * temp_len, temp_len, wt->params + N, len_cA, wt->params + N + len_cA);
1095 }
1096 N += 2 * len_cA;
1097 }
1098
1099 temp_len = wt->length[J - iter];
1100 p2 = 2 * p2;
1101 Np = N2;
1102 }
1103
1104 }
1105 else {
1106 printf("Signal extension can be either per or sym");
1107 exit(-1);
1108 }
1109
1110 J = wt->J;
1111 t2 = wt->outlength - 2 * wt->length[J];
1112 p2 = 2;
1113 it1 = 0;
1114 for (i = 0; i < J; ++i) {
1115 t = t2;
1116 for (k = 0; k < p2; ++k) {
1117 wt->nodelength[it1] = t;
1118 it1++;
1119 t += wt->length[J - i];
1120 }
1121 p2 *= 2;
1122 t2 = t2 - p2 * wt->length[J - i - 1];
1123 }
1124
1125 wt->coeflength[0] = wt->siglength;
1126
1127 for (i = 1; i < J + 1; ++i) {
1128 wt->coeflength[i] = wt->length[J - i + 1];
1129 }
1130
1131 free(orig);
1132 }
1133
ipow2(int n)1134 static int ipow2(int n) {
1135 int p,i;
1136 p = 1;
1137
1138 for (i = 0; i < n; ++i) {
1139 p *= 2;
1140 }
1141
1142 return p;
1143 }
1144
dwpt(wpt_object wt,const double * inp)1145 void dwpt(wpt_object wt, const double *inp) {
1146 int i, J, temp_len, iter, N, lp, p2, k, N2, Np;
1147 int temp, elength, temp2,size,nodes,llb,n1,j;
1148 double eparam,v1,v2;
1149 int len_cA, t, t2, it1,it2;
1150 double *orig,*tree;
1151 int *nodelength;
1152
1153 temp_len = wt->siglength;
1154 J = wt->J;
1155 wt->length[J + 1] = temp_len;
1156 wt->outlength = 0;
1157 temp = 1;
1158 elength = 0;
1159 size = wt->wave->filtlength;
1160 nodes = wt->nodes;
1161 n1 = nodes + 1;
1162 for (i = 0; i < J; ++i) {
1163 temp *= 2;
1164 temp2 = (size - 2) * (temp - 1);
1165 elength += temp2;
1166 }
1167 eparam = wt->eparam;
1168 orig = (double*)malloc(sizeof(double)* temp_len);
1169 tree = (double*)malloc(sizeof(double)* (temp_len * (J + 1) + elength));
1170 nodelength = (int*)malloc(sizeof(int)* nodes);
1171
1172 for (i = 0; i < wt->siglength; ++i) {
1173 orig[i] = inp[i];
1174 }
1175
1176 for (i = 0; i < temp_len * (J + 1) + elength; ++i) {
1177 tree[i] = 0.0;
1178 }
1179
1180 for (i = 0; i < nodes + 1; ++i) {
1181 wt->basisvector[i] = 0.0;
1182 wt->costvalues[i] = 0.0;
1183 }
1184
1185 N = temp_len;
1186 lp = wt->wave->lpd_len;
1187 p2 = 1;
1188
1189 //set eparam value here
1190 wt->costvalues[0] = costfunc(orig, wt->siglength, wt->entropy, eparam);
1191 it2 = 1;
1192 if (!strcmp(wt->ext, "per")) {
1193 i = J;
1194 p2 = 2;
1195 while (i > 0) {
1196 N = (int)ceil((double)N / 2.0);
1197 wt->length[i] = N;
1198 wt->outlength += p2 * (wt->length[i]);
1199 i--;
1200 p2 *= 2;
1201 }
1202 wt->length[0] = wt->length[1];
1203
1204 N2 = N = wt->outlength;
1205 p2 = 1;
1206 for (iter = 0; iter < J; ++iter) {
1207 len_cA = wt->length[J - iter];
1208 N2 -= 2 * p2 * len_cA;
1209 N = N2;
1210 for (k = 0; k < p2; ++k) {
1211 if (iter == 0) {
1212 dwpt_per(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA);
1213 }
1214 else {
1215 dwpt_per(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA);
1216 }
1217 wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
1218 it2++;
1219 wt->costvalues[it2] = costfunc(tree + N +len_cA, len_cA, wt->entropy, eparam);
1220 it2++;
1221 N += 2 * len_cA;
1222 }
1223
1224 temp_len = wt->length[J - iter];
1225 p2 = 2 * p2;
1226 Np = N2;
1227 }
1228 }
1229 else if (!strcmp(wt->ext, "sym")) {
1230 //printf("\n YES %s \n", wt->ext);
1231 i = J;
1232 p2 = 2;
1233 while (i > 0) {
1234 N = N + lp - 2;
1235 N = (int)ceil((double)N / 2.0);
1236 wt->length[i] = N;
1237 wt->outlength += p2 * (wt->length[i]);
1238 i--;
1239 p2 *= 2;
1240 }
1241 wt->length[0] = wt->length[1];
1242
1243 N2 = N = wt->outlength;
1244 p2 = 1;
1245
1246 for (iter = 0; iter < J; ++iter) {
1247 len_cA = wt->length[J - iter];
1248 N2 -= 2 * p2 * len_cA;
1249 N = N2;
1250 for (k = 0; k < p2; ++k) {
1251 if (iter == 0) {
1252 dwpt_sym(wt, orig, temp_len, tree + N, len_cA, tree + N + len_cA);
1253 }
1254 else {
1255 dwpt_sym(wt, tree + Np + k * temp_len, temp_len, tree + N, len_cA, tree + N + len_cA);
1256 }
1257 wt->costvalues[it2] = costfunc(tree + N, len_cA, wt->entropy, eparam);
1258 it2++;
1259 wt->costvalues[it2] = costfunc(tree + N + len_cA, len_cA, wt->entropy, eparam);
1260 it2++;
1261 N += 2 * len_cA;
1262 }
1263
1264 temp_len = wt->length[J - iter];
1265 p2 = 2 * p2;
1266 Np = N2;
1267 }
1268
1269 }
1270 else {
1271 printf("Signal extension can be either per or sym");
1272 exit(-1);
1273 }
1274
1275 J = wt->J;
1276 t2 = wt->outlength - 2 * wt->length[J];
1277 p2 = 2;
1278 it1 = 0;
1279 for (i = 0; i < J; ++i) {
1280 t = t2;
1281 for (k = 0; k < p2; ++k) {
1282 nodelength[it1] = t;
1283 it1++;
1284 t += wt->length[J - i];
1285 }
1286 p2 *= 2;
1287 t2 = t2 - p2 * wt->length[J - i - 1];
1288 }
1289
1290
1291 J = wt->J;
1292 llb = 1;
1293 for (i = 0; i < J; ++i) {
1294 llb *= 2;
1295 }
1296
1297 for (i = n1 - llb; i < n1; ++i) {
1298 wt->basisvector[i] = 1;
1299 }
1300
1301 for (j = J - 1; j >= 0; --j) {
1302 for (k = ipow2(j) - 1; k < ipow2(j + 1) - 1; ++k) {
1303 v1 = wt->costvalues[k];
1304 v2 = wt->costvalues[2 * k + 1] + wt->costvalues[2 * k + 2];
1305 //printf(" %g %g", v1,v2);
1306 if (v1 <= v2) {
1307 wt->basisvector[k] = 1;
1308 }
1309 else {
1310 wt->costvalues[k] = v2;
1311 }
1312 }
1313 //printf("\n");
1314 }
1315
1316 for (k = 0; k < nodes / 2; ++k) {
1317 if (wt->basisvector[k] == 1 || wt->basisvector[k] == 2) {
1318 wt->basisvector[2 * k + 1] = 2;
1319 wt->basisvector[2 * k + 2] = 2;
1320 }
1321 }
1322
1323 for (k = 0; k < n1; ++k) {
1324 if (wt->basisvector[k] == 2) {
1325 wt->basisvector[k] = 0;
1326 }
1327 }
1328
1329 N2 = 0;
1330 it1 = n1;
1331 it2 = 0;
1332 wt->nodes = 0;
1333 wt->numnodeslevel[0] = 0;
1334 //printf("Start \n");
1335
1336 if (wt->basisvector[0] == 1) {
1337 wt->outlength = wt->siglength;
1338 for (i = 0; i < wt->siglength; ++i) {
1339 wt->output[i] = inp[i];
1340 }
1341 wt->nodes = 1;
1342 wt->nodeindex[0] = 0;
1343 wt->nodeindex[1] = 0;
1344 wt->numnodeslevel[0] = 1;
1345 }
1346 else {
1347 for (i = J; i > 0; --i) {
1348 llb = ipow2(i);
1349 it1 -= llb;
1350 wt->numnodeslevel[i] = 0;
1351 for (j = 0; j < llb; ++j) {
1352 if (wt->basisvector[it1 + j] == 1) {
1353 //printf("NODE %d %d %d \n", i, j, wt->length[J - i + 1]);
1354 wt->nodeindex[2 * wt->nodes] = i;
1355 wt->nodeindex[2 * wt->nodes + 1] = j;
1356 wt->nodes += 1;
1357 wt->numnodeslevel[i] += 1;
1358 for (k = 0; k < wt->length[J - i + 1]; ++k) {
1359 wt->output[it2 + k] = tree[nodelength[it1 - 1 + j] + k];// access tree
1360 }
1361 it2 += wt->length[J - i + 1];
1362 }
1363 }
1364 }
1365 wt->outlength = it2;
1366 }
1367
1368 wt->coeflength[0] = wt->siglength;
1369
1370 for (i = 1; i < J + 1; ++i) {
1371 wt->coeflength[i] = wt->length[J - i + 1];
1372 }
1373
1374 free(orig);
1375 free(tree);
1376 free(nodelength);
1377 }
1378
getWTREENodelength(wtree_object wt,int X)1379 int getWTREENodelength(wtree_object wt, int X) {
1380 int N;
1381 N = -1;
1382 /*
1383 X - Level. All Nodes at any level have the same length
1384 */
1385 if (X <= 0 || X > wt->J) {
1386 printf("X co-ordinate must be >= 1 and <= %d", wt->J);
1387 exit(-1);
1388 }
1389
1390 N = wt->length[wt->J -X + 1];
1391
1392 return N;
1393 }
1394
getDWPTNodelength(wpt_object wt,int X)1395 int getDWPTNodelength(wpt_object wt, int X) {
1396 int N;
1397 N = -1;
1398 /*
1399 X - Level. All Nodes at any level have the same length
1400 */
1401 if (X <= 0 || X > wt->J) {
1402 printf("X co-ordinate must be >= 1 and <= %d", wt->J);
1403 exit(-1);
1404 }
1405
1406 N = wt->length[wt->J - X + 1];
1407
1408 return N;
1409 }
1410
getWTREECoeffs(wtree_object wt,int X,int Y,double * coeffs,int N)1411 void getWTREECoeffs(wtree_object wt, int X,int Y,double *coeffs,int N) {
1412 int ymax,i,t,t2;
1413
1414 if (X <= 0 || X > wt->J) {
1415 printf("X co-ordinate must be >= 1 and <= %d", wt->J);
1416 exit(-1);
1417 }
1418 ymax = 1;
1419 for (i = 0; i < X; ++i) {
1420 ymax *= 2;
1421 }
1422
1423 ymax -= 1;
1424
1425 if (Y < 0 ||Y > ymax) {
1426 printf("Y co-ordinate must be >= 0 and <= %d", ymax);
1427 exit(-1);
1428 }
1429
1430 if (X == 1) {
1431 t = 0;
1432 }
1433 else {
1434 t = 0;
1435 t2 = 1;
1436 for (i = 0; i < X - 1; ++i) {
1437 t2 *= 2;
1438 t += t2;
1439 }
1440 }
1441
1442 t += Y;
1443 t2 = wt->nodelength[t];
1444 for (i = 0; i < N; ++i) {
1445 coeffs[i] = wt->output[t2+i];
1446 }
1447
1448 }
1449
getDWPTCoeffs(wpt_object wt,int X,int Y,double * coeffs,int N)1450 void getDWPTCoeffs(wpt_object wt, int X, int Y, double *coeffs, int N) {
1451 int ymax, i;
1452 int np,citer;
1453 int flag;
1454
1455 if (X <= 0 || X > wt->J) {
1456 printf("X co-ordinate must be >= 1 and <= %d", wt->J);
1457 exit(-1);
1458 }
1459 ymax = 1;
1460 for (i = 0; i < X; ++i) {
1461 ymax *= 2;
1462 }
1463
1464 ymax -= 1;
1465
1466 if (Y < 0 || Y > ymax) {
1467 printf("Y co-ordinate must be >= 0 and <= %d", ymax);
1468 exit(-1);
1469 }
1470
1471 np = 0;
1472 citer = 0;
1473
1474 for (i = wt->J; i > X; --i) {
1475 np += wt->numnodeslevel[i];
1476 citer += wt->numnodeslevel[i] * wt->coeflength[i];
1477 }
1478
1479 i = 0;
1480 flag = 0;
1481 for (i = 0; i < wt->numnodeslevel[X]; ++i) {
1482 if (wt->nodeindex[2 * np + 1] == Y) {
1483 flag = 1;
1484 break;
1485 }
1486 np++;
1487 citer += wt->coeflength[X];
1488 }
1489
1490 if (flag == 0) {
1491 printf("The Node is Not Part Of The Best Basis Tree Use wpt_summary function to list available nodes \n");
1492 exit(-1);
1493 }
1494
1495 for (i = 0; i < N; ++i) {
1496 coeffs[i] = wt->output[citer + i];
1497 }
1498
1499 }
1500
getCWTScaleLength(int N)1501 int getCWTScaleLength(int N) {
1502 int J;
1503 double temp,dj;
1504
1505 dj = 0.4875;
1506
1507 temp = (log((double)N / 2.0) / log(2.0)) / dj;
1508 J = (int)temp;
1509
1510 return J;
1511 }
1512
setCWTScales(cwt_object wt,double s0,double dj,const char * type,int power)1513 void setCWTScales(cwt_object wt, double s0, double dj,const char *type,int power) {
1514 int i;
1515 strcpy(wt->type,type);
1516 //s0*pow(2.0, (double)(j - 1)*dj);
1517 if (!strcmp(wt->type, "pow") || !strcmp(wt->type, "power")) {
1518 for (i = 0; i < wt->J; ++i) {
1519 wt->scale[i] = s0*pow((double) power, (double)(i)*dj);
1520 }
1521 wt->sflag = 1;
1522 wt->pow = power;
1523
1524 }
1525 else if (!strcmp(wt->type, "lin") || !strcmp(wt->type, "linear")) {
1526 for (i = 0; i < wt->J; ++i) {
1527 wt->scale[i] = s0 + (double)i * dj;
1528 }
1529 wt->sflag = 1;
1530 }
1531 else {
1532 printf("\n Type accepts only two values : pow and lin\n");
1533 exit(-1);
1534 }
1535 wt->s0 = s0;
1536 wt->dj = dj;
1537 }
1538
setCWTScaleVector(cwt_object wt,const double * scale,int J,double s0,double dj)1539 void setCWTScaleVector(cwt_object wt, const double *scale, int J,double s0,double dj) {
1540 int i;
1541
1542 if (J != wt->J) {
1543 printf("\n CWT object is only valid for %d scales\n", wt->J);
1544 exit(-1);
1545 }
1546
1547 for (i = 0; i < wt->J; ++i) {
1548 wt->scale[i] = scale[i];
1549 }
1550 wt->dj = dj;
1551 wt->s0 = s0;
1552 wt->sflag = 1;
1553 }
1554
setCWTPadding(cwt_object wt,int pad)1555 void setCWTPadding(cwt_object wt, int pad) {
1556 if (pad == 0) {
1557 wt->pflag = 0;
1558 }
1559 else {
1560 wt->pflag = 1;
1561 }
1562 }
1563
cwt(cwt_object wt,const double * inp)1564 void cwt(cwt_object wt, const double *inp) {
1565 int i, N, npad,nj2,j,j2;
1566 N = wt->siglength;
1567 if (wt->sflag == 0) {
1568 for (i = 0; i < wt->J; ++i) {
1569 wt->scale[i] = wt->s0*pow(2.0, (double)(i)*wt->dj);
1570 }
1571 wt->sflag = 1;
1572 }
1573
1574 if (wt->pflag == 0) {
1575 npad = N;
1576 }
1577 else {
1578 npad = wt->npad;
1579 }
1580
1581 nj2 = 2 * N * wt->J;
1582 j = wt->J;
1583 j2 = 2 * j;
1584
1585 wt->smean = 0.0;
1586
1587 for (i = 0; i < N; ++i) {
1588 wt->smean += inp[i];
1589 }
1590 wt->smean /= N;
1591
1592 cwavelet(inp, N, wt->dt, wt->mother, wt->m, wt->s0,wt->dj,wt->J,npad,wt->params, wt->params+nj2, wt->params+nj2+j, wt->params+nj2+j2);
1593
1594 }
1595
icwt(cwt_object wt,double * cwtop)1596 void icwt(cwt_object wt, double *cwtop) {
1597 double psi, cdel;
1598 int real,i,N,nj2;
1599
1600 N = wt->siglength;
1601 nj2 = N * 2 * wt->J;
1602
1603 psi0(wt->mother, wt->m, &psi, &real);
1604 cdel = cdelta(wt->mother, wt->m, psi);
1605
1606 //printf("\n PSI %g CDEL %g param %g mother %d \n", psi, cdel,wt->m,wt->mother);
1607 if ((!strcmp(wt->type, "pow") || !strcmp(wt->type, "power")) && wt->pow == 2) {
1608 icwavelet(wt->params, N, wt->params+nj2, wt->J, wt->dt, wt->dj, cdel, psi, cwtop);
1609 } else {
1610 printf("Inverse CWT is only available for power of 2.0 scales \n");
1611 exit(-1);
1612 }
1613 for(i = 0; i < N;++i) {
1614 cwtop[i] += wt->smean;
1615 }
1616
1617 }
1618
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)1619 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) {
1620 int len_avg, N, U,N2,i;
1621
1622 len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
1623 N = 2 * len_cD;
1624 U = 2;
1625
1626 upsamp2(cA, len_cA, U, cA_up);
1627
1628 per_ext(cA_up, 2 * len_cA, len_avg / 2, temp);
1629
1630 N2 = 2 * len_cA + len_avg;
1631
1632 if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
1633 wt->cobj = conv_init(N2, len_avg);
1634 wt->cfftset = 1;
1635 }
1636 else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) {
1637 printf("Decomposition Filters must have the same length.");
1638 exit(-1);
1639 }
1640
1641 wconv(wt, temp, N2, wt->wave->lpr, len_avg, X_lp);
1642
1643 upsamp2(cD, len_cD, U, cA_up);
1644
1645 per_ext(cA_up, 2 * len_cD, len_avg / 2, temp);
1646
1647 N2 = 2 * len_cD + len_avg;
1648
1649 wconv(wt, temp, N2, wt->wave->hpr, len_avg, X_hp);
1650
1651
1652 for (i = len_avg - 1; i < N + len_avg - 1; ++i) {
1653 X[i - len_avg + 1] = X_lp[i] + X_hp[i];
1654 }
1655
1656 if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
1657 free_conv(wt->cobj);
1658 wt->cfftset = 0;
1659 }
1660
1661 }
1662
idwt_per(wt_object wt,double * cA,int len_cA,double * cD,double * X)1663 static void idwt_per(wt_object wt, double *cA, int len_cA, double *cD, double *X) {
1664 idwt_per_stride(cA,len_cA,cD, wt->wave->lpr, wt->wave->hpr, wt->wave->lpr_len, X,1,1);
1665 }
1666
idwt_sym(wt_object wt,double * cA,int len_cA,double * cD,double * X)1667 static void idwt_sym(wt_object wt, double *cA, int len_cA, double *cD, double *X) {
1668 idwt_sym_stride(cA,len_cA,cD, wt->wave->lpr, wt->wave->hpr, wt->wave->lpr_len, X,1,1);
1669 }
1670
1671
idwt(wt_object wt,double * dwtop)1672 void idwt(wt_object wt, double *dwtop) {
1673 int J,U,i,lf,N,N2,iter,k;
1674 int app_len, det_len;
1675 double *cA_up, *X_lp, *X_hp,*out,*temp;
1676
1677 J = wt->J;
1678 U = 2;
1679 app_len = wt->length[0];
1680 out = (double*)malloc(sizeof(double)* (wt->siglength + 1));
1681 if (!strcmp(wt->ext, "per") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
1682 app_len = wt->length[0];
1683 det_len = wt->length[1];
1684 N = 2 * wt->length[J];
1685 lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
1686
1687 cA_up = (double*)malloc(sizeof(double)* N);
1688 temp = (double*)malloc(sizeof(double)* (N + lf));
1689 X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
1690 X_hp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
1691 iter = app_len;
1692
1693 for (i = 0; i < app_len; ++i) {
1694 out[i] = wt->output[i];
1695 }
1696
1697 for (i = 0; i < J; ++i) {
1698
1699 idwt1(wt,temp,cA_up,out,det_len,wt->output+iter,det_len,X_lp,X_hp,out);
1700 /*
1701 idwt_per(wt,out, det_len, wt->output + iter, det_len, X_lp);
1702 for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
1703 out[k - lf/2 + 1] = X_lp[k];
1704 }
1705 */
1706 iter += det_len;
1707 det_len = wt->length[i + 2];
1708 }
1709 free(cA_up);
1710 free(X_lp);
1711 free(X_hp);
1712 free(temp);
1713
1714 }
1715 else if (!strcmp(wt->ext, "per") && !strcmp(wt->cmethod, "direct")) {
1716 app_len = wt->length[0];
1717 det_len = wt->length[1];
1718 N = 2 * wt->length[J];
1719 lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
1720
1721 X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
1722 iter = app_len;
1723
1724 for (i = 0; i < app_len; ++i) {
1725 out[i] = wt->output[i];
1726 }
1727
1728 for (i = 0; i < J; ++i) {
1729
1730 //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
1731
1732 idwt_per(wt,out, det_len, wt->output + iter, X_lp);
1733 for (k = lf/2 - 1; k < 2 * det_len + lf/2 - 1; ++k) {
1734 out[k - lf/2 + 1] = X_lp[k];
1735 }
1736
1737 iter += det_len;
1738 det_len = wt->length[i + 2];
1739 }
1740
1741 free(X_lp);
1742
1743 }
1744 else if (!strcmp(wt->ext, "sym") && !strcmp(wt->cmethod, "direct")) {
1745 app_len = wt->length[0];
1746 det_len = wt->length[1];
1747 N = 2 * wt->length[J] - 1;
1748 lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
1749
1750 X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
1751 iter = app_len;
1752
1753 for (i = 0; i < app_len; ++i) {
1754 out[i] = wt->output[i];
1755 }
1756
1757 for (i = 0; i < J; ++i) {
1758
1759 //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
1760
1761 idwt_sym(wt, out, det_len, wt->output + iter, X_lp);
1762 for (k = lf-2; k < 2 * det_len; ++k) {
1763 out[k - lf + 2] = X_lp[k];
1764 }
1765
1766 iter += det_len;
1767 det_len = wt->length[i + 2];
1768 }
1769
1770 free(X_lp);
1771
1772 }
1773 else if (!strcmp(wt->ext, "sym") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
1774 lf = wt->wave->lpd_len;// lpd and hpd have the same length
1775
1776 N = 2 * wt->length[J] - 1;
1777 cA_up = (double*)malloc(sizeof(double)* N);
1778 X_lp = (double*)malloc(sizeof(double)* (N + lf - 1));
1779 X_hp = (double*)malloc(sizeof(double)* (N + lf - 1));
1780
1781 for (i = 0; i < app_len; ++i) {
1782 out[i] = wt->output[i];
1783 }
1784
1785 iter = app_len;
1786
1787 for (i = 0; i < J; ++i) {
1788 det_len = wt->length[i + 1];
1789 upsamp(out, det_len, U, cA_up);
1790 N2 = 2 * wt->length[i + 1] - 1;
1791
1792 if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
1793 wt->cobj = conv_init(N2, lf);
1794 wt->cfftset = 1;
1795 }
1796 else if (!(wt->wave->lpr_len == wt->wave->hpr_len)) {
1797 printf("Decomposition Filters must have the same length.");
1798 exit(-1);
1799 }
1800
1801 wconv(wt, cA_up, N2, wt->wave->lpr, lf, X_lp);
1802
1803 upsamp(wt->output + iter, det_len, U, cA_up);
1804
1805 wconv(wt, cA_up, N2, wt->wave->hpr, lf, X_hp);
1806
1807 for (k = lf - 2; k < N2 + 1; ++k) {
1808 out[k - lf + 2] = X_lp[k] + X_hp[k];
1809 }
1810 iter += det_len;
1811 if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
1812 free_conv(wt->cobj);
1813 wt->cfftset = 0;
1814 }
1815 }
1816
1817 free(cA_up);
1818 free(X_lp);
1819 free(X_hp);
1820 }
1821 else {
1822 printf("Signal extension can be either per or sym");
1823 exit(-1);
1824 }
1825
1826 for (i = 0; i < wt->siglength; ++i) {
1827 dwtop[i] = out[i];
1828 }
1829
1830
1831 free(out);
1832
1833 }
1834
idwpt_per(wpt_object wt,double * cA,int len_cA,double * cD,double * X)1835 static void idwpt_per(wpt_object wt, double *cA, int len_cA, double *cD, double *X) {
1836 int len_avg, i, l, m, n, t, l2;
1837
1838 len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
1839 l2 = len_avg / 2;
1840 m = -2;
1841 n = -1;
1842
1843 for (i = 0; i < len_cA + l2 - 1; ++i) {
1844 m += 2;
1845 n += 2;
1846 X[m] = 0.0;
1847 X[n] = 0.0;
1848 for (l = 0; l < l2; ++l) {
1849 t = 2 * l;
1850 if ((i - l) >= 0 && (i - l) < len_cA) {
1851 X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
1852 X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
1853 }
1854 else if ((i - l) >= len_cA && (i - l) < len_cA + len_avg - 1) {
1855 X[m] += wt->wave->lpr[t] * cA[i - l - len_cA] + wt->wave->hpr[t] * cD[i - l - len_cA];
1856 X[n] += wt->wave->lpr[t + 1] * cA[i - l - len_cA] + wt->wave->hpr[t + 1] * cD[i - l - len_cA];
1857 }
1858 else if ((i - l) < 0 && (i - l) > -l2) {
1859 X[m] += wt->wave->lpr[t] * cA[len_cA + i - l] + wt->wave->hpr[t] * cD[len_cA + i - l];
1860 X[n] += wt->wave->lpr[t + 1] * cA[len_cA + i - l] + wt->wave->hpr[t + 1] * cD[len_cA + i - l];
1861 }
1862 }
1863 }
1864 }
1865
idwpt_sym(wpt_object wt,double * cA,int len_cA,double * cD,double * X)1866 static void idwpt_sym(wpt_object wt, double *cA, int len_cA, double *cD, double *X) {
1867 int len_avg, i, l, m, n, t, v;
1868
1869 len_avg = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
1870 m = -2;
1871 n = -1;
1872
1873 for (v = 0; v < len_cA; ++v) {
1874 i = v;
1875 m += 2;
1876 n += 2;
1877 X[m] = 0.0;
1878 X[n] = 0.0;
1879 for (l = 0; l < len_avg / 2; ++l) {
1880 t = 2 * l;
1881 if ((i - l) >= 0 && (i - l) < len_cA) {
1882 X[m] += wt->wave->lpr[t] * cA[i - l] + wt->wave->hpr[t] * cD[i - l];
1883 X[n] += wt->wave->lpr[t + 1] * cA[i - l] + wt->wave->hpr[t + 1] * cD[i - l];
1884 }
1885 }
1886 }
1887 }
1888
idwpt(wpt_object wt,double * dwtop)1889 void idwpt(wpt_object wt, double *dwtop) {
1890 int J, i, lf, k,p,l;
1891 int app_len, det_len, index, n1, llb, index2, index3, index4,indexp,xlen;
1892 double *X_lp, *X, *out, *out2;
1893 int *prep,*ptemp;
1894
1895 J = wt->J;
1896 app_len = wt->length[0];
1897 p = ipow2(J);
1898 lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
1899 xlen = p * (app_len + 2 * lf);
1900
1901 X_lp = (double*)malloc(sizeof(double)* 2 * (wt->length[J] + lf));
1902 X = (double*)malloc(sizeof(double)* xlen);
1903 out = (double*)malloc(sizeof(double)* wt->length[J]);
1904 out2 = (double*)malloc(sizeof(double)* wt->length[J]);
1905 prep = (int*)malloc(sizeof(int)* p);
1906 ptemp = (int*)malloc(sizeof(int)* p);
1907 n1 = 1;
1908 llb = 1;
1909 index2 = xlen / p;
1910 indexp = 0;
1911 if (wt->basisvector[0] == 1) {
1912 for (i = 0; i < wt->siglength; ++i) {
1913 dwtop[i] = wt->output[i];
1914 }
1915
1916 }
1917 else {
1918 for (i = 0; i < J; ++i) {
1919 llb *= 2;
1920 n1 += llb;
1921 }
1922
1923 for (i = 0; i < xlen; ++i) {
1924 X[i] = 0.0;
1925 }
1926
1927 for (i = 0; i < llb; ++i) {
1928 prep[i] = (int)wt->basisvector[n1 - llb + i];
1929 ptemp[i] = 0;
1930 }
1931
1932 if (!strcmp(wt->ext, "per")) {
1933 app_len = wt->length[0];
1934 det_len = wt->length[1];
1935 index = 0;
1936
1937
1938 for (i = 0; i < J; ++i) {
1939 p = ipow2(J - i - 1);
1940 det_len = wt->length[i + 1];
1941 index2 *= 2;
1942 index3 = 0;
1943 index4 = 0;
1944 //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
1945 n1 -= llb;
1946 for (l = 0; l < llb; ++l) {
1947 if (ptemp[l] != 2) {
1948 prep[l] = (int)wt->basisvector[n1 + l];
1949 }
1950 else {
1951 prep[l] = ptemp[l];
1952 }
1953 ptemp[l] = 0;
1954 }
1955
1956
1957 for (l = 0; l < p; ++l) {
1958 if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) {
1959 for (k = 0; k < det_len; ++k) {
1960 out[k] = wt->output[index + k];
1961 out2[k] = wt->output[index + det_len + k];
1962 }
1963 idwpt_per(wt, out, det_len, out2, X_lp);
1964 for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
1965 X[index3 + k - lf / 2 + 1] = X_lp[k];
1966 }
1967 index += 2 * det_len;
1968 index3 += index2;
1969 index4 += 2 * indexp;
1970 ptemp[l] = 2;
1971 }
1972 else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) {
1973 index4 += indexp;
1974 for (k = 0; k < det_len; ++k) {
1975 out[k] = wt->output[index + k];
1976 out2[k] = X[index4 + k];
1977 }
1978 idwpt_per(wt, out, det_len, out2, X_lp);
1979 for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
1980 X[index3 + k - lf / 2 + 1] = X_lp[k];
1981 }
1982 index += det_len;
1983 index3 += index2;
1984 index4 += indexp;
1985 ptemp[l] = 2;
1986 }
1987 else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) {
1988 for (k = 0; k < det_len; ++k) {
1989 out[k] = X[index4 + k];
1990 out2[k] = wt->output[index + k];
1991 }
1992 idwpt_per(wt, out, det_len, out2, X_lp);
1993 for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
1994 X[index3 + k - lf / 2 + 1] = X_lp[k];
1995 }
1996 index += det_len;
1997 index3 += index2;
1998 index4 += 2 * indexp;
1999 ptemp[l] = 2;
2000 }
2001 else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) {
2002 for (k = 0; k < det_len; ++k) {
2003 out[k] = X[index4 + k];
2004 out2[k] = X[index4 + indexp + k];
2005 }
2006 idwpt_per(wt, out, det_len, out2, X_lp);
2007 for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
2008 X[index3 + k - lf / 2 + 1] = X_lp[k];
2009 }
2010 index4 += 2 * indexp;
2011 index3 += index2;
2012 ptemp[l] = 2;
2013 }
2014 else {
2015 index3 += index2;
2016 index4 += 2 * indexp;
2017 }
2018
2019 }
2020
2021
2022 /*
2023 idwt_per(wt, out, det_len, wt->output + iter, det_len, X_lp);
2024 for (k = lf / 2 - 1; k < 2 * det_len + lf / 2 - 1; ++k) {
2025 out[k - lf / 2 + 1] = X_lp[k];
2026 }
2027
2028 iter += det_len;
2029 det_len = wt->length[i + 2];
2030 */
2031 llb /= 2;
2032 indexp = index2;
2033 }
2034
2035 //free(X_lp);
2036
2037 }
2038 else if (!strcmp(wt->ext, "sym")) {
2039 app_len = wt->length[0];
2040 det_len = wt->length[1];
2041
2042 //X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
2043 index = 0;
2044
2045 for (i = 0; i < J; ++i) {
2046 p = ipow2(J - i - 1);
2047 det_len = wt->length[i + 1];
2048 index2 *= 2;
2049 index3 = 0;
2050 index4 = 0;
2051 //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
2052 n1 -= llb;
2053 for (l = 0; l < llb; ++l) {
2054 if (ptemp[l] != 2) {
2055 prep[l] = (int)wt->basisvector[n1 + l];
2056 }
2057 else {
2058 prep[l] = ptemp[l];
2059 }
2060 ptemp[l] = 0;
2061 }
2062
2063
2064 for (l = 0; l < p; ++l) {
2065 if (prep[2 * l] == 1 && prep[2 * l + 1] == 1) {
2066 for (k = 0; k < det_len; ++k) {
2067 out[k] = wt->output[index + k];
2068 out2[k] = wt->output[index + det_len + k];
2069 }
2070 idwpt_sym(wt, out, det_len, out2, X_lp);
2071 for (k = lf - 2; k < 2 * det_len; ++k) {
2072 X[index3 + k - lf + 2] = X_lp[k];
2073 }
2074 index += 2 * det_len;
2075 index3 += index2;
2076 index4 += 2 * indexp;
2077 ptemp[l] = 2;
2078 }
2079 else if (prep[2 * l] == 1 && prep[2 * l + 1] == 2) {
2080 index4 += indexp;
2081 for (k = 0; k < det_len; ++k) {
2082 out[k] = wt->output[index + k];
2083 out2[k] = X[index4 + k];
2084 }
2085 idwpt_sym(wt, out, det_len, out2, X_lp);
2086 for (k = lf - 2; k < 2 * det_len; ++k) {
2087 X[index3 + k - lf + 2] = X_lp[k];
2088 }
2089 index += det_len;
2090 index3 += index2;
2091 index4 += indexp;
2092 ptemp[l] = 2;
2093 }
2094 else if (prep[2 * l] == 2 && prep[2 * l + 1] == 1) {
2095 for (k = 0; k < det_len; ++k) {
2096 out[k] = X[index4 + k];
2097 out2[k] = wt->output[index + k];
2098 }
2099 idwpt_sym(wt, out, det_len, out2, X_lp);
2100 for (k = lf - 2; k < 2 * det_len; ++k) {
2101 X[index3 + k - lf + 2] = X_lp[k];
2102 }
2103 index += det_len;
2104 index3 += index2;
2105 index4 += 2 * indexp;
2106 ptemp[l] = 2;
2107 }
2108 else if (prep[2 * l] == 2 && prep[2 * l + 1] == 2) {
2109 for (k = 0; k < det_len; ++k) {
2110 out[k] = X[index4 + k];
2111 out2[k] = X[index4 + indexp + k];
2112 }
2113 idwpt_sym(wt, out, det_len, out2, X_lp);
2114 for (k = lf - 2; k < 2 * det_len; ++k) {
2115 X[index3 + k - lf + 2] = X_lp[k];
2116 }
2117 index4 += 2 * indexp;
2118 index3 += index2;
2119 ptemp[l] = 2;
2120 }
2121 else {
2122 index3 += index2;
2123 index4 += 2 * indexp;
2124 }
2125
2126 }
2127
2128 //idwt1(wt, temp, cA_up, out, det_len, wt->output + iter, det_len, X_lp, X_hp, out);
2129 /*
2130 idwpt_sym(wt, out, det_len, wt->output + iter, det_len, X_lp);
2131 for (k = lf - 2; k < 2 * det_len; ++k) {
2132 out[k - lf + 2] = X_lp[k];
2133 }
2134
2135 iter += det_len;
2136 det_len = wt->length[i + 2];
2137 */
2138 llb /= 2;
2139 indexp = index2;
2140 }
2141
2142 //free(X_lp);
2143
2144 }
2145 else {
2146 printf("Signal extension can be either per or sym");
2147 exit(-1);
2148 }
2149
2150 for (i = 0; i < wt->siglength; ++i) {
2151 //printf("%g ", X[i]);
2152 dwtop[i] = X[i];
2153 }
2154
2155 }
2156
2157
2158 free(out);
2159 free(X_lp);
2160 free(X);
2161 free(out2);
2162 free(prep);
2163 free(ptemp);
2164 }
2165
2166
swt_per(wt_object wt,int M,double * inp,int N,double * cA,int len_cA,double * cD)2167 static void swt_per(wt_object wt,int M, double *inp, int N, double *cA, int len_cA, double *cD) {
2168
2169 swt_per_stride(M,inp,N,wt->wave->lpd,wt->wave->hpd, wt->wave->lpd_len, cA,len_cA,cD,1,1);
2170 }
2171
swt_fft(wt_object wt,const double * inp)2172 static void swt_fft(wt_object wt, const double *inp) {
2173 int i, J, temp_len, iter, M, N, len_filt;
2174 int lenacc;
2175 double *low_pass, *high_pass,*sig,*cA,*cD;
2176
2177 temp_len = wt->siglength;
2178 J = wt->J;
2179 wt->length[0] = wt->length[J] = temp_len;
2180 wt->outlength = wt->length[J+1] = (J + 1) * temp_len;
2181 M = 1;
2182 for (iter = 1; iter < J; ++iter) {
2183 M = 2 * M;
2184 wt->length[iter] = temp_len;
2185 }
2186
2187 len_filt = wt->wave->filtlength;
2188
2189 low_pass = (double*)malloc(sizeof(double)* M * len_filt);
2190 high_pass = (double*)malloc(sizeof(double)* M * len_filt);
2191 sig = (double*)malloc(sizeof(double)* (M * len_filt + temp_len + (temp_len%2)));
2192 cA = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1);
2193 cD = (double*)malloc(sizeof(double)* (2 * M * len_filt + temp_len + (temp_len % 2)) - 1);
2194
2195 M = 1;
2196
2197 for (i = 0; i < temp_len; ++i) {
2198 wt->params[i] = inp[i];
2199 }
2200
2201 lenacc = wt->outlength;
2202
2203 for (iter = 0; iter < J; ++iter) {
2204 lenacc -= temp_len;
2205 if (iter > 0) {
2206 M = 2 * M;
2207 N = M * len_filt;
2208 upsamp2(wt->wave->lpd, wt->wave->lpd_len, M, low_pass);
2209 upsamp2(wt->wave->hpd, wt->wave->hpd_len, M, high_pass);
2210 }
2211 else {
2212 N = len_filt;
2213 for (i = 0; i < N; ++i) {
2214 low_pass[i] = wt->wave->lpd[i];
2215 high_pass[i] = wt->wave->hpd[i];
2216 }
2217 }
2218
2219 //swt_per(wt,M, wt->params, temp_len, cA, temp_len, cD,temp_len);
2220
2221 per_ext(wt->params, temp_len, N / 2, sig);
2222
2223 if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
2224 wt->cobj = conv_init(N + temp_len + (temp_len % 2), N);
2225 wt->cfftset = 1;
2226 }
2227 else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
2228 printf("Decomposition Filters must have the same length.");
2229 exit(-1);
2230 }
2231
2232
2233 wconv(wt, sig, N + temp_len + (temp_len % 2), low_pass, N, cA);
2234
2235 wconv(wt, sig, N + temp_len + (temp_len % 2), high_pass, N, cD);
2236
2237 if (wt->wave->lpd_len == wt->wave->hpd_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
2238 free_conv(wt->cobj);
2239 wt->cfftset = 0;
2240 }
2241
2242 for (i = 0; i < temp_len; ++i) {
2243 wt->params[i] = cA[N + i];
2244 wt->params[lenacc + i] = cD[N + i];
2245 }
2246
2247
2248 }
2249
2250 free(low_pass);
2251 free(high_pass);
2252 free(sig);
2253 free(cA);
2254 free(cD);
2255 }
2256
swt_direct(wt_object wt,const double * inp)2257 static void swt_direct(wt_object wt, const double *inp) {
2258 int i, J, temp_len, iter, M;
2259 int lenacc;
2260 double *cA, *cD;
2261
2262 temp_len = wt->siglength;
2263 J = wt->J;
2264 wt->length[0] = wt->length[J] = temp_len;
2265 wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
2266 M = 1;
2267 for (iter = 1; iter < J; ++iter) {
2268 M = 2 * M;
2269 wt->length[iter] = temp_len;
2270 }
2271
2272
2273 cA = (double*)malloc(sizeof(double)* temp_len);
2274 cD = (double*)malloc(sizeof(double)* temp_len);
2275
2276 M = 1;
2277
2278 for (i = 0; i < temp_len; ++i) {
2279 wt->params[i] = inp[i];
2280 }
2281
2282 lenacc = wt->outlength;
2283
2284 for (iter = 0; iter < J; ++iter) {
2285 lenacc -= temp_len;
2286 if (iter > 0) {
2287 M = 2 * M;
2288 }
2289
2290 swt_per(wt, M, wt->params, temp_len, cA, temp_len, cD);
2291
2292
2293 for (i = 0; i < temp_len; ++i) {
2294 wt->params[i] = cA[i];
2295 wt->params[lenacc + i] = cD[i];
2296 }
2297
2298 }
2299
2300 free(cA);
2301 free(cD);
2302
2303 }
2304
2305
swt(wt_object wt,const double * inp)2306 void swt(wt_object wt, const double *inp) {
2307 if (!strcmp(wt->method, "swt") && !strcmp(wt->cmethod, "direct") ) {
2308 swt_direct(wt,inp);
2309 }
2310 else if (!strcmp(wt->method, "swt") && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
2311 swt_fft(wt, inp);
2312 }
2313 else {
2314 printf("SWT Only accepts two methods - direct and fft");
2315 exit(-1);
2316 }
2317 }
2318
getSWTRecCoeff(double * coeff,int * length,const char * ctype,int level,int J,double * lpr,double * hpr,int lf,int siglength,double * swtop)2319 static void getSWTRecCoeff(double *coeff, int *length, const char *ctype, int level, int J, double *lpr,
2320 double *hpr, int lf, int siglength, double *swtop) {
2321 int N, iter, i, index, value, count, len;
2322 int index_shift, len0, U, N1, index2;
2323 double *appx1, *det1, *appx_sig, *det_sig, *cL0, *cH0, *tempx, *oup00L, *oup00H, *oup00, *oup01, *appx2, *det2;
2324
2325 N = siglength;
2326 U = 2;
2327
2328 appx_sig = (double*)malloc(sizeof(double)* N);
2329 det_sig = (double*)malloc(sizeof(double)* N);
2330 appx1 = (double*)malloc(sizeof(double)* N);
2331 det1 = (double*)malloc(sizeof(double)* N);
2332 appx2 = (double*)malloc(sizeof(double)* N);
2333 det2 = (double*)malloc(sizeof(double)* N);
2334 tempx = (double*)malloc(sizeof(double)* N);
2335 cL0 = (double*)malloc(sizeof(double)* (N + (N % 2) + lf));
2336 cH0 = (double*)malloc(sizeof(double)* (N + (N % 2) + lf));
2337 oup00L = (double*)malloc(sizeof(double)* (N + 2 * lf));
2338 oup00H = (double*)malloc(sizeof(double)* (N + 2 * lf));
2339 oup00 = (double*)malloc(sizeof(double)* N);
2340 oup01 = (double*)malloc(sizeof(double)* N);
2341
2342
2343
2344 for (iter = J-level; iter < J; ++iter) {
2345 for (i = 0; i < N; ++i) {
2346 swtop[i] = 0.0;
2347 }
2348 if (!strcmp((ctype), "appx") && (iter == (J-level))) {
2349 for (i = 0; i < N; ++i) {
2350 appx_sig[i] = coeff[i];
2351 det_sig[i] = 0.0;
2352 }
2353 }
2354 else if (!strcmp((ctype), "det") && (iter == (J-level))) {
2355 for (i = 0; i < N; ++i) {
2356 det_sig[i] = coeff[i];
2357 appx_sig[i] = 0.0;
2358 }
2359 }
2360 else {
2361 for (i = 0; i < N; ++i) {
2362 det_sig[i] = 0.0;
2363 }
2364 }
2365
2366 value = (int)pow(2.0, (double)(J - 1 - iter));
2367
2368 for (count = 0; count < value; count++) {
2369 len = 0;
2370 for (index = count; index < N; index += value) {
2371 appx1[len] = appx_sig[index];
2372 det1[len] = det_sig[index];
2373 len++;
2374 }
2375
2376
2377 //SHIFT 0
2378 len0 = 0;
2379
2380 for (index_shift = 0; index_shift < len; index_shift += 2) {
2381 appx2[len0] = appx1[index_shift];
2382 det2[len0] = det1[index_shift];
2383 len0++;
2384 }
2385 upsamp2(appx2, len0, U, tempx);
2386 per_ext(tempx, 2 * len0, lf / 2, cL0);
2387
2388 upsamp2(det2, len0, U, tempx);
2389 per_ext(tempx, 2 * len0, lf / 2, cH0);
2390
2391 N1 = 2 * len0 + lf;
2392
2393
2394
2395 conv_direct(cL0, N1, lpr, lf, oup00L);
2396
2397
2398 conv_direct(cH0, N1, hpr, lf, oup00H);
2399
2400 for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
2401 oup00[i - lf + 1] = oup00L[i] + oup00H[i];
2402 }
2403
2404 //SHIFT 1
2405
2406 len0 = 0;
2407
2408 for (index_shift = 1; index_shift < len; index_shift += 2) {
2409 appx2[len0] = appx1[index_shift];
2410 det2[len0] = det1[index_shift];
2411 len0++;
2412 }
2413
2414 upsamp2(appx2, len0, U, tempx);
2415 per_ext(tempx, 2 * len0, lf / 2, cL0);
2416
2417 upsamp2(det2, len0, U, tempx);
2418 per_ext(tempx, 2 * len0, lf / 2, cH0);
2419
2420 N1 = 2 * len0 + lf;
2421
2422 conv_direct(cL0, N1, lpr, lf, oup00L);
2423
2424 conv_direct(cH0, N1, hpr, lf, oup00H);
2425
2426 for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
2427 oup01[i - lf + 1] = oup00L[i] + oup00H[i];
2428 }
2429
2430 circshift(oup01, 2 * len0, -1);
2431
2432 index2 = 0;
2433
2434 for (index = count; index < N; index += value) {
2435 swtop[index] = (oup00[index2] + oup01[index2]) / 2.0;
2436 index2++;
2437 }
2438
2439 }
2440 for (i = 0; i < N; ++i) {
2441 appx_sig[i] = swtop[i];
2442 }
2443
2444 }
2445
2446
2447
2448 free(appx_sig);
2449 free(det_sig);
2450 free(appx1);
2451 free(det1);
2452 free(tempx);
2453 free(cL0);
2454 free(cH0);
2455 free(oup00L);
2456 free(oup00H);
2457 free(oup00);
2458 free(oup01);
2459 free(appx2);
2460 free(det2);
2461 }
2462
getSWTmra(wt_object wt,double * wavecoeffs)2463 double *getSWTmra(wt_object wt, double *wavecoeffs) {
2464 int i, J, access, N;
2465 double *mra;
2466 J = wt->J;
2467 mra = (double*)malloc(sizeof(double)* wt->siglength*(J + 1));
2468 access = 0;
2469
2470
2471 // Approximation MRA
2472 getSWTRecCoeff(wt->output + access, wt->length, "appx", J, J, wt->wave->lpr, wt->wave->hpr, wt->wave->lpr_len, wt->siglength, mra);
2473 // Details MRA
2474 N = wt->siglength;
2475
2476 for (i = J; i > 0; --i) {
2477 access += wt->length[J - i];
2478 getSWTRecCoeff(wt->output + access, wt->length, "det", i, J, wt->wave->lpr, wt->wave->hpr, wt->wave->lpr_len, wt->siglength, mra+N);
2479 N += wt->siglength;
2480 }
2481
2482 return mra;
2483 }
2484
iswt(wt_object wt,double * swtop)2485 void iswt(wt_object wt, double *swtop) {
2486 int N, lf, iter,i,J,index,value,count,len;
2487 int index_shift,len0,U,N1,index2;
2488 double *appx1, *det1,*appx_sig,*det_sig,*cL0,*cH0,*tempx,*oup00L,*oup00H,*oup00,*oup01,*appx2,*det2;
2489
2490 N = wt->siglength;
2491 J = wt->J;
2492 U = 2;
2493 lf = wt->wave->lpr_len;
2494
2495 appx_sig = (double*)malloc(sizeof(double)* N);
2496 det_sig = (double*)malloc(sizeof(double)* N);
2497 appx1 = (double*)malloc(sizeof(double)* N);
2498 det1 = (double*)malloc(sizeof(double)* N);
2499 appx2 = (double*)malloc(sizeof(double)* N);
2500 det2 = (double*)malloc(sizeof(double)* N);
2501 tempx = (double*)malloc(sizeof(double)* N);
2502 cL0 = (double*)malloc(sizeof(double)* (N + (N%2) + lf));
2503 cH0 = (double*)malloc(sizeof(double)* (N + (N % 2) + lf));
2504 oup00L = (double*)malloc(sizeof(double)* (N + 2 * lf));
2505 oup00H = (double*)malloc(sizeof(double)* (N + 2 * lf));
2506 oup00 = (double*)malloc(sizeof(double)* N);
2507 oup01 = (double*)malloc(sizeof(double)* N);
2508
2509
2510
2511 for (iter = 0; iter < J; ++iter) {
2512 for (i = 0; i < N; ++i) {
2513 swtop[i] = 0.0;
2514 }
2515 if (iter == 0) {
2516 for (i = 0; i < N; ++i) {
2517 appx_sig[i] = wt->output[i];
2518 det_sig[i] = wt->output[N + i];
2519 }
2520 }
2521 else {
2522 for (i = 0; i < N; ++i) {
2523 det_sig[i] = wt->output[(iter + 1) * N + i];
2524 }
2525 }
2526
2527 value = (int)pow(2.0, (double) (J - 1 - iter));
2528
2529 for (count = 0; count < value; count++) {
2530 len = 0;
2531 for (index = count; index < N; index += value) {
2532 appx1[len] = appx_sig[index];
2533 det1[len] = det_sig[index];
2534 len++;
2535 }
2536
2537
2538 //SHIFT 0
2539 len0 = 0;
2540
2541 for (index_shift = 0; index_shift < len; index_shift += 2) {
2542 appx2[len0] = appx1[index_shift];
2543 det2[len0] = det1[index_shift];
2544 len0++;
2545 }
2546 upsamp2(appx2, len0, U, tempx);
2547 per_ext(tempx, 2 * len0, lf / 2, cL0);
2548
2549 upsamp2(det2, len0, U, tempx);
2550 per_ext(tempx, 2 * len0, lf / 2, cH0);
2551
2552 N1 = 2 * len0 + lf;
2553
2554 if (wt->wave->lpr_len == wt->wave->hpr_len && (!strcmp(wt->cmethod, "fft") || !strcmp(wt->cmethod, "FFT"))) {
2555 wt->cobj = conv_init(N1, lf);
2556 wt->cfftset = 1;
2557 }
2558 else if (!(wt->wave->lpd_len == wt->wave->hpd_len)) {
2559 printf("Decomposition Filters must have the same length.");
2560 exit(-1);
2561 }
2562
2563 wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L);
2564
2565 wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H);
2566
2567 for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
2568 oup00[i - lf + 1] = oup00L[i] + oup00H[i];
2569 }
2570
2571 //SHIFT 1
2572
2573 len0 = 0;
2574
2575 for (index_shift = 1; index_shift < len; index_shift += 2) {
2576 appx2[len0] = appx1[index_shift];
2577 det2[len0] = det1[index_shift];
2578 len0++;
2579 }
2580
2581 upsamp2(appx2, len0, U, tempx);
2582 per_ext(tempx, 2 * len0, lf / 2, cL0);
2583
2584 upsamp2(det2, len0, U, tempx);
2585 per_ext(tempx, 2 * len0, lf / 2, cH0);
2586
2587 N1 = 2 * len0 + lf;
2588
2589 wconv(wt, cL0, N1, wt->wave->lpr, lf, oup00L);
2590
2591 wconv(wt, cH0, N1, wt->wave->hpr, lf, oup00H);
2592
2593 for (i = lf - 1; i < 2 * len0 + lf - 1; ++i) {
2594 oup01[i - lf + 1] = oup00L[i] + oup00H[i];
2595 }
2596
2597 circshift(oup01, 2*len0, -1);
2598
2599 index2 = 0;
2600
2601 for (index = count; index < N; index += value) {
2602 swtop[index] = (oup00[index2] + oup01[index2]) / 2.0;
2603 index2++;
2604 }
2605
2606 }
2607 for (i = 0; i < N; ++i) {
2608 appx_sig[i] = swtop[i];
2609 }
2610
2611 }
2612
2613
2614
2615 free(appx_sig);
2616 free(det_sig);
2617 free(appx1);
2618 free(det1);
2619 free(tempx);
2620 free(cL0);
2621 free(cH0);
2622 free(oup00L);
2623 free(oup00H);
2624 free(oup00);
2625 free(oup01);
2626 free(appx2);
2627 free(det2);
2628 }
2629
modwt_per(wt_object wt,int M,double * inp,double * cA,int len_cA,double * cD)2630 static void modwt_per(wt_object wt, int M, double *inp, double *cA, int len_cA, double *cD) {
2631 int l, i, t, len_avg;
2632 double s;
2633 double *filt;
2634 len_avg = wt->wave->lpd_len;
2635
2636 filt = (double*)malloc(sizeof(double)* 2 * len_avg);
2637 s = sqrt(2.0);
2638 for (i = 0; i < len_avg; ++i) {
2639 filt[i] = wt->wave->lpd[i] / s;
2640 filt[len_avg + i] = wt->wave->hpd[i] / s;
2641 }
2642
2643 for (i = 0; i < len_cA; ++i) {
2644 t = i;
2645 cA[i] = filt[0] * inp[t];
2646 cD[i] = filt[len_avg] * inp[t];
2647 for (l = 1; l < len_avg; l++) {
2648 t -= M;
2649 while (t >= len_cA) {
2650 t -= len_cA;
2651 }
2652 while (t < 0) {
2653 t += len_cA;
2654 }
2655
2656 cA[i] += filt[l] * inp[t];
2657 cD[i] += filt[len_avg + l] * inp[t];
2658
2659 }
2660 }
2661 free(filt);
2662 }
2663
modwt_direct(wt_object wt,const double * inp)2664 static void modwt_direct(wt_object wt, const double *inp) {
2665 int i, J, temp_len, iter, M;
2666 int lenacc;
2667 double *cA, *cD;
2668
2669 if (strcmp(wt->ext, "per")) {
2670 printf("MODWT direct method only uses periodic extension per. \n");
2671 printf(" Use MODWT fft method for symmetric extension sym \n");
2672 exit(-1);
2673 }
2674
2675 temp_len = wt->siglength;
2676 J = wt->J;
2677 wt->length[0] = wt->length[J] = temp_len;
2678 wt->outlength = wt->length[J + 1] = (J + 1) * temp_len;
2679 M = 1;
2680 for (iter = 1; iter < J; ++iter) {
2681 M = 2 * M;
2682 wt->length[iter] = temp_len;
2683 }
2684
2685
2686 cA = (double*)malloc(sizeof(double)* temp_len);
2687 cD = (double*)malloc(sizeof(double)* temp_len);
2688
2689 M = 1;
2690
2691 for (i = 0; i < temp_len; ++i) {
2692 wt->params[i] = inp[i];
2693 }
2694
2695 lenacc = wt->outlength;
2696
2697 for (iter = 0; iter < J; ++iter) {
2698 lenacc -= temp_len;
2699 if (iter > 0) {
2700 M = 2 * M;
2701 }
2702
2703 modwt_per(wt, M, wt->params, cA, temp_len, cD);
2704
2705
2706 for (i = 0; i < temp_len; ++i) {
2707 wt->params[i] = cA[i];
2708 wt->params[lenacc + i] = cD[i];
2709 }
2710
2711 }
2712
2713 free(cA);
2714 free(cD);
2715
2716 }
2717
modwt_fft(wt_object wt,const double * inp)2718 static void modwt_fft(wt_object wt, const double *inp) {
2719 int i, J, temp_len, iter, M,N, len_avg;
2720 int lenacc;
2721 double s,tmp1,tmp2;
2722 fft_data *cA, *cD, *low_pass,*high_pass,*sig;
2723 int *index;
2724 fft_object fft_fd = NULL;
2725 fft_object fft_bd = NULL;
2726
2727 temp_len = wt->siglength;
2728 len_avg = wt->wave->lpd_len;
2729 if (!strcmp(wt->ext, "sym")) {
2730 N = 2 * temp_len;
2731 } else if (!strcmp(wt->ext, "per")) {
2732 N = temp_len;
2733 }
2734 J = wt->J;
2735 wt->modwtsiglength = N;
2736 wt->length[0] = wt->length[J] = N;
2737 wt->outlength = wt->length[J + 1] = (J + 1) * N;
2738
2739 s = sqrt(2.0);
2740 for (iter = 1; iter < J; ++iter) {
2741 wt->length[iter] = N;
2742 }
2743
2744 fft_fd = fft_init(N, 1);
2745 fft_bd = fft_init(N, -1);
2746
2747 sig = (fft_data*)malloc(sizeof(fft_data)* N);
2748 cA = (fft_data*)malloc(sizeof(fft_data)* N);
2749 cD = (fft_data*)malloc(sizeof(fft_data)* N);
2750 low_pass = (fft_data*)malloc(sizeof(fft_data)* N);
2751 high_pass = (fft_data*)malloc(sizeof(fft_data)* N);
2752 index = (int*)malloc(sizeof(int)*N);
2753
2754
2755 // N-point FFT of low pass and high pass filters
2756
2757 // Low Pass Filter
2758
2759 for(i = 0; i < len_avg;++i) {
2760 sig[i].re = (fft_type) wt->wave->lpd[i] / s;
2761 sig[i].im = 0.0;
2762 }
2763 for(i = len_avg; i < N;++i) {
2764 sig[i].re = 0.0;
2765 sig[i].im = 0.0;
2766 }
2767
2768
2769 fft_exec(fft_fd, sig, low_pass);
2770
2771 // High Pass Filter
2772
2773 for (i = 0; i < len_avg; ++i) {
2774 sig[i].re = (fft_type)wt->wave->hpd[i] / s;
2775 sig[i].im = 0.0;
2776 }
2777 for (i = len_avg; i < N; ++i) {
2778 sig[i].re = 0.0;
2779 sig[i].im = 0.0;
2780 }
2781
2782 fft_exec(fft_fd, sig, high_pass);
2783
2784 // symmetric extension
2785 for (i = 0; i < temp_len; ++i) {
2786 sig[i].re = (fft_type)inp[i];
2787 sig[i].im = 0.0;
2788 }
2789 for (i = temp_len; i < N; ++i) {
2790 sig[i].re = (fft_type) inp[N-i-1];
2791 sig[i].im = 0.0;
2792 }
2793
2794 // FFT of data
2795
2796 fft_exec(fft_fd, sig, cA);
2797
2798 lenacc = wt->outlength;
2799
2800 M = 1;
2801
2802 for (iter = 0; iter < J; ++iter) {
2803 lenacc -= N;
2804
2805 for (i = 0; i < N; ++i) {
2806 index[i] = (M *i) % N;
2807 }
2808
2809 for (i = 0; i < N; ++i) {
2810 tmp1 = cA[i].re;
2811 tmp2 = cA[i].im;
2812 cA[i].re = low_pass[index[i]].re*tmp1 - low_pass[index[i]].im*tmp2;
2813 cA[i].im = low_pass[index[i]].re*tmp2 + low_pass[index[i]].im*tmp1;
2814
2815 cD[i].re = high_pass[index[i]].re*tmp1 - high_pass[index[i]].im*tmp2;
2816 cD[i].im = high_pass[index[i]].re*tmp2 + high_pass[index[i]].im*tmp1;
2817 }
2818
2819 fft_exec(fft_bd, cD, sig);
2820
2821 for (i = 0; i < N; ++i) {
2822 wt->params[lenacc + i] = sig[i].re/N;
2823 }
2824
2825 M *= 2;
2826 }
2827
2828 fft_exec(fft_bd, cA, sig);
2829
2830 for (i = 0; i < N; ++i) {
2831 wt->params[i] = sig[i].re/N;
2832 }
2833
2834 free(sig);
2835 free(cA);
2836 free(cD);
2837 free(low_pass);
2838 free(high_pass);
2839 free_fft(fft_fd);
2840 free_fft(fft_bd);
2841 }
2842
modwt(wt_object wt,const double * inp)2843 void modwt(wt_object wt, const double *inp) {
2844 if (!strcmp(wt->cmethod, "direct")) {
2845 modwt_direct(wt, inp);
2846 }
2847 else if (!strcmp(wt->cmethod, "fft")) {
2848 modwt_fft(wt, inp);
2849 }
2850 else {
2851 printf("Error- Available Choices for this method are - direct and fft \n");
2852 exit(-1);
2853 }
2854 }
2855
conj_complex(fft_data * x,int N)2856 static void conj_complex(fft_data *x, int N) {
2857 int i;
2858
2859 for (i = 0; i < N; ++i) {
2860 x[i].im *= (-1.0);
2861 }
2862 }
2863
getMODWTRecCoeff(fft_object fft_fd,fft_object fft_bd,fft_data * appx,fft_data * det,fft_data * cA,fft_data * cD,int * index,const char * ctype,int level,int J,fft_data * low_pass,fft_data * high_pass,int N)2864 static void getMODWTRecCoeff(fft_object fft_fd, fft_object fft_bd, fft_data* appx,fft_data *det, fft_data* cA,fft_data *cD, int *index, const char *ctype, int level,
2865 int J, fft_data *low_pass, fft_data *high_pass, int N) {
2866
2867 int iter,M,i;
2868 fft_type tmp1, tmp2;
2869
2870 M = (int)pow(2.0, (double)level - 1.0);
2871
2872 if (!strcmp((ctype), "appx")) {
2873 for (iter = 0; iter < level; ++iter) {
2874 fft_exec(fft_fd, appx, cA);
2875 fft_exec(fft_fd, det, cD);
2876
2877 for (i = 0; i < N; ++i) {
2878 index[i] = (M *i) % N;
2879 }
2880
2881 for (i = 0; i < N; ++i) {
2882 tmp1 = cA[i].re;
2883 tmp2 = cA[i].im;
2884 cA[i].re = low_pass[index[i]].re*tmp1 - low_pass[index[i]].im*tmp2 + high_pass[index[i]].re*cD[i].re - high_pass[index[i]].im*cD[i].im;
2885 cA[i].im = low_pass[index[i]].re*tmp2 + low_pass[index[i]].im*tmp1 + high_pass[index[i]].re*cD[i].im + high_pass[index[i]].im*cD[i].re;
2886 }
2887
2888 fft_exec(fft_bd, cA, appx);
2889
2890 for (i = 0; i < N; ++i) {
2891 appx[i].re /= N;
2892 appx[i].im /= N;
2893 }
2894
2895 M /= 2;
2896 }
2897 }
2898 else if (!strcmp(ctype, "det")) {
2899 for (iter = 0; iter < level; ++iter) {
2900 fft_exec(fft_fd, appx, cA);
2901 fft_exec(fft_fd, det, cD);
2902
2903 for (i = 0; i < N; ++i) {
2904 index[i] = (M *i) % N;
2905 }
2906
2907 for (i = 0; i < N; ++i) {
2908 tmp1 = cA[i].re;
2909 tmp2 = cA[i].im;
2910 cA[i].re = low_pass[index[i]].re*tmp1 - low_pass[index[i]].im*tmp2 + high_pass[index[i]].re*cD[i].re - high_pass[index[i]].im*cD[i].im;
2911 cA[i].im = low_pass[index[i]].re*tmp2 + low_pass[index[i]].im*tmp1 + high_pass[index[i]].re*cD[i].im + high_pass[index[i]].im*cD[i].re;
2912 }
2913
2914 fft_exec(fft_bd, cA, appx);
2915
2916 for (i = 0; i < N; ++i) {
2917 appx[i].re /= N;
2918 appx[i].im /= N;
2919 det[i].re = 0.0;
2920 det[i].im = 0.0;
2921 }
2922
2923 M /= 2;
2924 }
2925 }
2926 else {
2927 printf("ctype can only be one of appx or det \n");
2928 exit(-1);
2929 }
2930
2931
2932 }
2933
getMODWTmra(wt_object wt,double * wavecoeffs)2934 double* getMODWTmra(wt_object wt, double *wavecoeffs) {
2935 double *mra;
2936 int i, J, temp_len, iter, M, N, len_avg,lmra;
2937 int lenacc;
2938 double s;
2939 fft_data *cA, *cD, *low_pass, *high_pass, *sig,*ninp;
2940 int *index;
2941 fft_object fft_fd = NULL;
2942 fft_object fft_bd = NULL;
2943
2944 N = wt->modwtsiglength;
2945 len_avg = wt->wave->lpd_len;
2946 if (!strcmp(wt->ext, "sym")) {
2947 temp_len = N / 2;
2948 }
2949 else if (!strcmp(wt->ext, "per")) {
2950 temp_len = N;
2951 }
2952 J = wt->J;
2953
2954 s = sqrt(2.0);
2955 fft_fd = fft_init(N, 1);
2956 fft_bd = fft_init(N, -1);
2957
2958 sig = (fft_data*)malloc(sizeof(fft_data)* N);
2959 cA = (fft_data*)malloc(sizeof(fft_data)* N);
2960 cD = (fft_data*)malloc(sizeof(fft_data)* N);
2961 ninp = (fft_data*)malloc(sizeof(fft_data)* N);
2962 low_pass = (fft_data*)malloc(sizeof(fft_data)* N);
2963 high_pass = (fft_data*)malloc(sizeof(fft_data)* N);
2964 index = (int*)malloc(sizeof(int)*N);
2965 mra = (double*)malloc(sizeof(double)*temp_len*(J + 1));
2966
2967 // N-point FFT of low pass and high pass filters
2968
2969 // Low Pass Filter
2970
2971 for (i = 0; i < len_avg; ++i) {
2972 sig[i].re = (fft_type)wt->wave->lpd[i] / s;
2973 sig[i].im = 0.0;
2974 }
2975 for (i = len_avg; i < N; ++i) {
2976 sig[i].re = 0.0;
2977 sig[i].im = 0.0;
2978 }
2979
2980
2981 fft_exec(fft_fd, sig, low_pass);
2982
2983 // High Pass Filter
2984
2985 for (i = 0; i < len_avg; ++i) {
2986 sig[i].re = (fft_type)wt->wave->hpd[i] / s;
2987 sig[i].im = 0.0;
2988 }
2989 for (i = len_avg; i < N; ++i) {
2990 sig[i].re = 0.0;
2991 sig[i].im = 0.0;
2992 }
2993
2994 fft_exec(fft_fd, sig, high_pass);
2995
2996
2997 // Complex conjugate of the two filters
2998
2999 conj_complex(low_pass, N);
3000 conj_complex(high_pass, N);
3001
3002 M = (int)pow(2.0, (double)J - 1.0);
3003 lenacc = N;
3004
3005 //
3006 for (i = 0; i < N; ++i) {
3007 sig[i].re = (fft_type)wt->output[i];
3008 sig[i].im = 0.0;
3009 ninp[i].re = 0.0;
3010 ninp[i].im = 0.0;
3011 }
3012
3013 // Find Approximation MRA
3014
3015 getMODWTRecCoeff(fft_fd,fft_bd,sig,ninp,cA,cD,index,"appx",J,J,low_pass,high_pass,N);
3016
3017 for (i = 0; i < wt->siglength; ++i) {
3018 mra[i] = sig[i].re;
3019 }
3020 lmra = wt->siglength;
3021 // Find Details MRA
3022
3023 for (iter = 0; iter < J; ++iter) {
3024 for (i = 0; i < N; ++i) {
3025 sig[i].re = (fft_type)wt->output[lenacc+i];
3026 sig[i].im = 0.0;
3027 ninp[i].re = 0.0;
3028 ninp[i].im = 0.0;
3029 }
3030
3031 getMODWTRecCoeff(fft_fd, fft_bd, sig, ninp,cA, cD, index, "det", J-iter, J, low_pass, high_pass, N);
3032
3033 for (i = 0; i < wt->siglength; ++i) {
3034 mra[lmra+i] = sig[i].re;
3035 }
3036
3037 lenacc += N;
3038 lmra += wt->siglength;
3039 }
3040
3041 free(ninp);
3042 free(index);
3043 free(sig);
3044 free(cA);
3045 free(cD);
3046 free(low_pass);
3047 free(high_pass);
3048 free_fft(fft_fd);
3049 free_fft(fft_bd);
3050
3051 return mra;
3052 }
3053
imodwt_fft(wt_object wt,double * oup)3054 void imodwt_fft(wt_object wt, double *oup) {
3055 int i, J, temp_len, iter, M, N, len_avg;
3056 int lenacc;
3057 double s, tmp1, tmp2;
3058 fft_data *cA, *cD, *low_pass, *high_pass, *sig;
3059 int *index;
3060 fft_object fft_fd = NULL;
3061 fft_object fft_bd = NULL;
3062
3063 N = wt->modwtsiglength;
3064 len_avg = wt->wave->lpd_len;
3065 if (!strcmp(wt->ext, "sym")) {
3066 temp_len = N/2;
3067 }
3068 else if (!strcmp(wt->ext, "per")) {
3069 temp_len = N;
3070 }
3071 J = wt->J;
3072
3073 s = sqrt(2.0);
3074 fft_fd = fft_init(N, 1);
3075 fft_bd = fft_init(N, -1);
3076
3077 sig = (fft_data*)malloc(sizeof(fft_data)* N);
3078 cA = (fft_data*)malloc(sizeof(fft_data)* N);
3079 cD = (fft_data*)malloc(sizeof(fft_data)* N);
3080 low_pass = (fft_data*)malloc(sizeof(fft_data)* N);
3081 high_pass = (fft_data*)malloc(sizeof(fft_data)* N);
3082 index = (int*)malloc(sizeof(int)*N);
3083
3084
3085 // N-point FFT of low pass and high pass filters
3086
3087 // Low Pass Filter
3088
3089 for (i = 0; i < len_avg; ++i) {
3090 sig[i].re = (fft_type)wt->wave->lpd[i] / s;
3091 sig[i].im = 0.0;
3092 }
3093 for (i = len_avg; i < N; ++i) {
3094 sig[i].re = 0.0;
3095 sig[i].im = 0.0;
3096 }
3097
3098
3099 fft_exec(fft_fd, sig, low_pass);
3100
3101 // High Pass Filter
3102
3103 for (i = 0; i < len_avg; ++i) {
3104 sig[i].re = (fft_type)wt->wave->hpd[i] / s;
3105 sig[i].im = 0.0;
3106 }
3107 for (i = len_avg; i < N; ++i) {
3108 sig[i].re = 0.0;
3109 sig[i].im = 0.0;
3110 }
3111
3112 fft_exec(fft_fd, sig, high_pass);
3113
3114
3115 // Complex conjugate of the two filters
3116
3117 conj_complex(low_pass, N);
3118 conj_complex(high_pass, N);
3119
3120 M = (int)pow(2.0, (double)J - 1.0);
3121 lenacc = N;
3122
3123 //
3124 for (i = 0; i < N; ++i) {
3125 sig[i].re = (fft_type)wt->output[i];
3126 sig[i].im = 0.0;
3127 }
3128
3129 for (iter = 0; iter < J; ++iter) {
3130 fft_exec(fft_fd, sig, cA);
3131 for (i = 0; i < N; ++i) {
3132 sig[i].re = wt->output[lenacc+i];
3133 sig[i].im = 0.0;
3134 }
3135 fft_exec(fft_fd, sig, cD);
3136
3137 for (i = 0; i < N; ++i) {
3138 index[i] = (M *i) % N;
3139 }
3140
3141 for (i = 0; i < N; ++i) {
3142 tmp1 = cA[i].re;
3143 tmp2 = cA[i].im;
3144 cA[i].re = low_pass[index[i]].re*tmp1 - low_pass[index[i]].im*tmp2 + high_pass[index[i]].re*cD[i].re - high_pass[index[i]].im*cD[i].im;
3145 cA[i].im = low_pass[index[i]].re*tmp2 + low_pass[index[i]].im*tmp1 + high_pass[index[i]].re*cD[i].im + high_pass[index[i]].im*cD[i].re;
3146 }
3147
3148 fft_exec(fft_bd, cA, sig);
3149
3150 for (i = 0; i < N; ++i) {
3151 sig[i].re /= N;
3152 sig[i].im /= N;
3153 }
3154 M /= 2;
3155 lenacc += N;
3156 }
3157
3158 for (i = 0; i < wt->siglength; ++i) {
3159 oup[i] = sig[i].re;
3160 }
3161
3162 free(sig);
3163 free(cA);
3164 free(cD);
3165 free(low_pass);
3166 free(high_pass);
3167 free_fft(fft_fd);
3168 free_fft(fft_bd);
3169 }
3170
3171
imodwt_per(wt_object wt,int M,double * cA,int len_cA,double * cD,double * X)3172 static void imodwt_per(wt_object wt,int M, double *cA, int len_cA, double *cD, double *X) {
3173 int len_avg, i, l, t;
3174 double s;
3175 double *filt;
3176 len_avg = wt->wave->lpd_len;
3177
3178 filt = (double*)malloc(sizeof(double)* 2 * len_avg);
3179 s = sqrt(2.0);
3180 for (i = 0; i < len_avg; ++i) {
3181 filt[i] = wt->wave->lpd[i] / s;
3182 filt[len_avg + i] = wt->wave->hpd[i] / s;
3183 }
3184
3185
3186 for (i = 0; i < len_cA; ++i) {
3187 t = i;
3188 X[i] = (filt[0] * cA[t]) + (filt[len_avg] * cD[t]);
3189 for (l = 1; l < len_avg; l++) {
3190 t += M;
3191 while (t >= len_cA) {
3192 t -= len_cA;
3193 }
3194 while (t < 0) {
3195 t += len_cA;
3196 }
3197
3198 X[i] += (filt[l] * cA[t]) + (filt[len_avg + l] * cD[t]);
3199
3200 }
3201 }
3202 free(filt);
3203 }
3204
imodwt_direct(wt_object wt,double * dwtop)3205 static void imodwt_direct(wt_object wt, double *dwtop) {
3206 int N, iter, i, J, j;
3207 int lenacc,M;
3208 double *X;
3209
3210 N = wt->siglength;
3211 J = wt->J;
3212 lenacc = N;
3213 M = (int)pow(2.0, (double)J - 1.0);
3214 //M = 1;
3215 X = (double*)malloc(sizeof(double)* N);
3216
3217 for (i = 0; i < N; ++i) {
3218 dwtop[i] = wt->output[i];
3219 }
3220
3221 for (iter = 0; iter < J; ++iter) {
3222 if (iter > 0) {
3223 M = M / 2;
3224 }
3225 imodwt_per(wt, M, dwtop, N, wt->params + lenacc, X);
3226 /*
3227 for (j = lf - 1; j < N; ++j) {
3228 dwtop[j - lf + 1] = X[j];
3229 }
3230 for (j = 0; j < lf - 1; ++j) {
3231 dwtop[N - lf + 1 + j] = X[j];
3232 }
3233 */
3234 for (j = 0; j < N; ++j) {
3235 dwtop[j] = X[j];
3236 }
3237
3238 lenacc += N;
3239 }
3240 free(X);
3241 }
3242
imodwt(wt_object wt,double * oup)3243 void imodwt(wt_object wt, double *oup) {
3244 if (!strcmp(wt->cmethod, "direct")) {
3245 imodwt_direct(wt, oup);
3246 }
3247 else if (!strcmp(wt->cmethod, "fft")) {
3248 imodwt_fft(wt, oup);
3249 }
3250 else {
3251 printf("Error- Available Choices for this method are - direct and fft \n");
3252 exit(-1);
3253 }
3254 }
3255
3256
setDWTExtension(wt_object wt,const char * extension)3257 void setDWTExtension(wt_object wt, const char *extension) {
3258 if (!strcmp(extension, "sym")) {
3259 strcpy(wt->ext, "sym");
3260 }
3261 else if (!strcmp(extension, "per")) {
3262 strcpy(wt->ext, "per");
3263 }
3264 else {
3265 printf("Signal extension can be either per or sym");
3266 exit(-1);
3267 }
3268 }
3269
setWTREEExtension(wtree_object wt,const char * extension)3270 void setWTREEExtension(wtree_object wt, const char *extension) {
3271 if (!strcmp(extension, "sym")) {
3272 strcpy(wt->ext, "sym");
3273 }
3274 else if (!strcmp(extension, "per")) {
3275 strcpy(wt->ext, "per");
3276 }
3277 else {
3278 printf("Signal extension can be either per or sym");
3279 exit(-1);
3280 }
3281 }
3282
setDWPTExtension(wpt_object wt,const char * extension)3283 void setDWPTExtension(wpt_object wt, const char *extension) {
3284 if (!strcmp(extension, "sym")) {
3285 strcpy(wt->ext, "sym");
3286 }
3287 else if (!strcmp(extension, "per")) {
3288 strcpy(wt->ext, "per");
3289 }
3290 else {
3291 printf("Signal extension can be either per or sym");
3292 exit(-1);
3293 }
3294 }
3295
setDWT2Extension(wt2_object wt,const char * extension)3296 void setDWT2Extension(wt2_object wt, const char *extension) {
3297 if (!strcmp(wt->method, "dwt")) {
3298 if (!strcmp(extension, "sym")) {
3299 strcpy(wt->ext, "sym");
3300 }
3301 else if (!strcmp(extension, "per")) {
3302 strcpy(wt->ext, "per");
3303 }
3304 else {
3305 printf("Signal extension can be either per or sym");
3306 exit(-1);
3307 }
3308 }
3309 else if (!strcmp(wt->method, "swt") || !strcmp(wt->method, "modwt")) {
3310 if (!strcmp(extension, "per")) {
3311 strcpy(wt->ext, "per");
3312 }
3313 else {
3314 printf("Signal extension can only be per");
3315 exit(-1);
3316 }
3317 }
3318
3319 }
3320
setDWPTEntropy(wpt_object wt,const char * entropy,double eparam)3321 void setDWPTEntropy(wpt_object wt, const char *entropy, double eparam) {
3322 if (!strcmp(entropy, "shannon")) {
3323 strcpy(wt->entropy, "shannon");
3324 }
3325 else if (!strcmp(entropy, "threshold")) {
3326 strcpy(wt->entropy, "threshold");
3327 wt ->eparam = eparam;
3328 }
3329 else if (!strcmp(entropy, "norm")) {
3330 strcpy(wt->entropy, "norm");
3331 wt->eparam = eparam;
3332 }
3333 else if (!strcmp(entropy, "logenergy") || !strcmp(entropy, "log energy") || !strcmp(entropy, "energy")) {
3334 strcpy(wt->entropy, "logenergy");
3335 }
3336 else {
3337 printf("Entropy should be one of shannon, threshold, norm or logenergy");
3338 exit(-1);
3339 }
3340 }
3341
setWTConv(wt_object wt,const char * cmethod)3342 void setWTConv(wt_object wt, const char *cmethod) {
3343 if (!strcmp(cmethod, "fft") || !strcmp(cmethod, "FFT")) {
3344 strcpy(wt->cmethod, "fft");
3345 }
3346 else if (!strcmp(cmethod, "direct")) {
3347 strcpy(wt->cmethod, "direct");
3348 }
3349 else {
3350 printf("Convolution Only accepts two methods - direct and fft");
3351 exit(-1);
3352 }
3353 }
3354
dwt2(wt2_object wt,double * inp)3355 double* dwt2(wt2_object wt, double *inp) {
3356 double *wavecoeff;
3357 int i, J, iter, N, lp, rows_n, cols_n, rows_i, cols_i;
3358 int ir, ic, istride,ostride;
3359 int aLL, aLH, aHL, aHH, cdim,clen;
3360 double *orig, *lp_dn1,*hp_dn1;
3361 J = wt->J;
3362 wt->outlength = 0;
3363
3364 rows_n = wt->rows;
3365 cols_n = wt->cols;
3366 lp = wt->wave->lpd_len;
3367 clen = J * 3;
3368 if (!strcmp(wt->ext, "per")) {
3369 i = 2 * J;
3370 while (i > 0) {
3371 rows_n = (int)ceil((double)rows_n / 2.0);
3372 cols_n = (int)ceil((double)cols_n / 2.0);
3373 wt->dimensions[i - 1] = cols_n;
3374 wt->dimensions[i - 2] = rows_n;
3375 wt->outlength += (rows_n * cols_n) * 3;
3376 i = i - 2;
3377 }
3378 wt->outlength += (rows_n * cols_n);
3379 N = wt->outlength;
3380 wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
3381
3382 orig = inp;
3383 ir = wt->rows;
3384 ic = wt->cols;
3385 cols_i = wt->dimensions[2 * J - 1];
3386
3387 lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3388 hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3389
3390 for (iter = 0; iter < J; ++iter) {
3391 rows_i = wt->dimensions[2*J - 2*iter - 2];
3392 cols_i = wt->dimensions[2*J - 2*iter - 1];
3393 istride = 1;
3394 ostride = 1;
3395 cdim = rows_i * cols_i;
3396 // Row filtering and column subsampling
3397 for (i = 0; i < ir; ++i) {
3398 dwt_per_stride(orig+i*ic, ic, wt->wave->lpd,wt->wave->hpd,lp, lp_dn1+i*cols_i, cols_i, hp_dn1+i*cols_i, istride, ostride);
3399 }
3400
3401 // Column Filtering and Row subsampling
3402 aHH = N - cdim;
3403 wt->coeffaccess[clen] = aHH;
3404 aHL = aHH - cdim;
3405 wt->coeffaccess[clen-1] = aHL;
3406 aLH = aHL - cdim;
3407 wt->coeffaccess[clen-2] = aLH;
3408 aLL = aLH - cdim;
3409
3410 N -= 3 * cdim;
3411 ic = cols_i;
3412 istride = ic;
3413 ostride = ic;
3414
3415 for (i = 0; i < ic; ++i) {
3416 dwt_per_stride(lp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff+aLL+i, rows_i, wavecoeff+aLH+i, istride, ostride);
3417 }
3418
3419
3420 for (i = 0; i < ic; ++i) {
3421 dwt_per_stride(hp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff+aHL+i, rows_i, wavecoeff+aHH+i, istride, ostride);
3422 }
3423
3424 ir = rows_i;
3425 orig = wavecoeff+aLL;
3426 clen -= 3;
3427
3428 }
3429 wt->coeffaccess[0] = 0;
3430 free(lp_dn1);
3431 free(hp_dn1);
3432 }
3433 else if (!strcmp(wt->ext, "sym")) {
3434 i = 2 * J;
3435 while (i > 0) {
3436 rows_n += lp - 2;
3437 cols_n += lp - 2;
3438 rows_n = (int)ceil((double)rows_n / 2.0);
3439 cols_n = (int)ceil((double)cols_n / 2.0);
3440 wt->dimensions[i - 1] = cols_n;
3441 wt->dimensions[i - 2] = rows_n;
3442 wt->outlength += (rows_n * cols_n) * 3;
3443 i = i - 2;
3444 }
3445 wt->outlength += (rows_n * cols_n);
3446 N = wt->outlength;
3447 wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
3448
3449 orig = inp;
3450 ir = wt->rows;
3451 ic = wt->cols;
3452 cols_i = wt->dimensions[2 * J - 1];
3453
3454 lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3455 hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3456
3457 for (iter = 0; iter < J; ++iter) {
3458 rows_i = wt->dimensions[2 * J - 2 * iter - 2];
3459 cols_i = wt->dimensions[2 * J - 2 * iter - 1];
3460 istride = 1;
3461 ostride = 1;
3462 cdim = rows_i * cols_i;
3463 // Row filtering and column subsampling
3464 for (i = 0; i < ir; ++i) {
3465 dwt_sym_stride(orig + i*ic, ic, wt->wave->lpd, wt->wave->hpd, lp, lp_dn1 + i*cols_i, cols_i, hp_dn1 + i*cols_i, istride, ostride);
3466 }
3467
3468 // Column Filtering and Row subsampling
3469 aHH = N - cdim;
3470 wt->coeffaccess[clen] = aHH;
3471 aHL = aHH - cdim;
3472 wt->coeffaccess[clen - 1] = aHL;
3473 aLH = aHL - cdim;
3474 wt->coeffaccess[clen - 2] = aLH;
3475 aLL = aLH - cdim;
3476 N -= 3 * cdim;
3477 ic = cols_i;
3478 istride = ic;
3479 ostride = ic;
3480
3481 for (i = 0; i < ic; ++i) {
3482 dwt_sym_stride(lp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aLL + i, rows_i, wavecoeff + aLH + i, istride, ostride);
3483 }
3484
3485 for (i = 0; i < ic; ++i) {
3486 dwt_sym_stride(hp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aHL + i, rows_i, wavecoeff + aHH + i, istride, ostride);
3487 }
3488
3489 ir = rows_i;
3490 orig = wavecoeff + aLL;
3491 clen -= 3;
3492
3493 }
3494 wt->coeffaccess[0] = 0;
3495 free(lp_dn1);
3496 free(hp_dn1);
3497 }
3498
3499 return wavecoeff;
3500 }
3501
idwt2(wt2_object wt,double * wavecoeff,double * oup)3502 void idwt2(wt2_object wt, double *wavecoeff, double *oup) {
3503 int i, k, rows, cols, N, ir,ic,lf,dim1,dim2;
3504 int istride, ostride, iter, J;
3505 int aLL, aLH, aHL, aHH;
3506 double *cL, *cH, *X_lp,*orig;
3507
3508 rows = wt->rows;
3509 cols = wt->cols;
3510 J = wt->J;
3511 double *out;
3512
3513
3514 if (!strcmp(wt->ext, "per")) {
3515 N = rows > cols ? 2 * rows : 2 * cols;
3516 lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
3517
3518 i = J;
3519 dim1 = wt->dimensions[0];
3520 dim2 = wt->dimensions[1];
3521 k = 0;
3522 while (i > 0) {
3523 k += 1;
3524 dim1 *= 2;
3525 dim2 *= 2;
3526 i--;
3527 }
3528
3529
3530
3531 X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
3532 cL = (double*)calloc(dim1*dim2, sizeof(double));
3533 cH = (double*)calloc(dim1*dim2, sizeof(double));
3534 out = (double*)calloc(dim1*dim2, sizeof(double));
3535 aLL = wt->coeffaccess[0];
3536 orig = wavecoeff + aLL;
3537 for (iter = 0; iter < J; ++iter) {
3538 ir = wt->dimensions[2 * iter];
3539 ic = wt->dimensions[2 * iter + 1];
3540 istride = ic;
3541 ostride = 1;
3542 aLH = wt->coeffaccess[iter*3 + 1];
3543 aHL = wt->coeffaccess[iter*3 + 2];
3544 aHH = wt->coeffaccess[iter*3 + 3];
3545 for (i = 0; i < ic; ++i) {
3546 idwt_per_stride(orig+i, ir, wavecoeff+aLH+i, wt->wave->lpr, wt->wave->hpr, lf, X_lp,istride,ostride);
3547
3548 for (k = lf / 2 - 1; k < 2 * ir + lf / 2 - 1; ++k) {
3549 cL[(k - lf / 2 + 1)*ic + i] = X_lp[k];
3550 }
3551
3552 idwt_per_stride(wavecoeff + aHL+i, ir, wavecoeff + aHH+i, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
3553
3554 for (k = lf / 2 - 1; k < 2 * ir + lf / 2 - 1; ++k) {
3555 cH[(k - lf / 2 + 1)*ic + i] = X_lp[k];
3556 }
3557 }
3558
3559 ir *= 2;
3560 istride = 1;
3561 ostride = 1;
3562
3563 for (i = 0; i < ir; ++i) {
3564 idwt_per_stride(cL+i*ic, ic, cH+i*ic, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
3565
3566 for (k = lf / 2 - 1; k < 2 * ic + lf / 2 - 1; ++k) {
3567 out[(k - lf / 2 + 1) + i*ic*2] = X_lp[k];
3568 }
3569 }
3570 ic *= 2;
3571 if (iter == J - 1) {
3572 for (i = 0; i < wt->rows; ++i) {
3573 for (k = 0; k < wt->cols; ++k) {
3574 oup[k + i*wt->cols] = out[k + i*ic];
3575 }
3576 }
3577 }
3578 else {
3579 for (i = 0; i < wt->dimensions[2 * (iter+1)]; ++i) {
3580 for (k = 0; k < wt->dimensions[2 * (iter + 1)+1]; ++k) {
3581 oup[k + i*wt->dimensions[2 * (iter + 1) + 1]] = out[k + i*ic];
3582 }
3583 }
3584 }
3585
3586
3587 orig = oup;
3588 }
3589 free(X_lp);
3590 free(cL);
3591 free(cH);
3592 }
3593 else if (!strcmp(wt->ext, "sym")) {
3594 N = rows > cols ? 2 * rows - 1 : 2 * cols - 1;
3595 lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
3596
3597 i = J;
3598 dim1 = wt->dimensions[0];
3599 dim2 = wt->dimensions[1];
3600 k = 0;
3601 while (i > 0) {
3602 k += 1;
3603 dim1 *= 2;
3604 dim2 *= 2;
3605 i--;
3606 }
3607
3608
3609
3610 X_lp = (double*)malloc(sizeof(double)* (N + 2 * lf - 1));
3611 cL = (double*)calloc(dim1*dim2, sizeof(double));
3612 cH = (double*)calloc(dim1*dim2, sizeof(double));
3613 out = (double*)calloc(dim1*dim2, sizeof(double));
3614 aLL = wt->coeffaccess[0];
3615 orig = wavecoeff + aLL;
3616 for (iter = 0; iter < J; ++iter) {
3617 ir = wt->dimensions[2 * iter];
3618 ic = wt->dimensions[2 * iter + 1];
3619 istride = ic;
3620 ostride = 1;
3621 aLH = wt->coeffaccess[iter * 3 + 1];
3622 aHL = wt->coeffaccess[iter * 3 + 2];
3623 aHH = wt->coeffaccess[iter * 3 + 3];
3624 for (i = 0; i < ic; ++i) {
3625 idwt_sym_stride(orig + i, ir, wavecoeff + aLH + i, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
3626
3627 for (k = lf - 2; k < 2 * ir; ++k) {
3628 cL[(k - lf + 2)*ic + i] = X_lp[k];
3629 }
3630
3631 idwt_sym_stride(wavecoeff + aHL + i, ir, wavecoeff + aHH + i, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
3632
3633 for (k = lf - 2; k < 2 * ir; ++k) {
3634 cH[(k - lf + 2)*ic + i] = X_lp[k];
3635 }
3636 }
3637
3638 ir *= 2;
3639 istride = 1;
3640 ostride = 1;
3641
3642 for (i = 0; i < ir; ++i) {
3643 idwt_sym_stride(cL + i*ic, ic, cH + i*ic, wt->wave->lpr, wt->wave->hpr, lf, X_lp, istride, ostride);
3644
3645 for (k = lf - 2; k < 2 * ic; ++k) {
3646 out[(k - lf + 2) + i*ic * 2] = X_lp[k];
3647 }
3648 }
3649 ic *= 2;
3650 if (iter == J - 1) {
3651 for (i = 0; i < wt->rows; ++i) {
3652 for (k = 0; k < wt->cols; ++k) {
3653 oup[k + i*wt->cols] = out[k + i*ic];
3654 }
3655 }
3656 }
3657 else {
3658 for (i = 0; i < wt->dimensions[2 * (iter + 1)]; ++i) {
3659 for (k = 0; k < wt->dimensions[2 * (iter + 1) + 1]; ++k) {
3660 oup[k + i*wt->dimensions[2 * (iter + 1) + 1]] = out[k + i*ic];
3661 }
3662 }
3663 }
3664
3665
3666 orig = oup;
3667 }
3668 free(X_lp);
3669 free(cL);
3670 free(cH);
3671 }
3672
3673 free(out);
3674
3675 }
3676
swt2(wt2_object wt,double * inp)3677 double* swt2(wt2_object wt, double *inp) {
3678 double *wavecoeff;
3679 int i, J, iter, M, N, lp, rows_n, cols_n, rows_i, cols_i;
3680 int ir, ic, istride, ostride;
3681 int aLL, aLH, aHL, aHH, cdim, clen;
3682 double *orig, *lp_dn1, *hp_dn1;
3683
3684 J = wt->J;
3685 M = 1;
3686 wt->outlength = 0;
3687
3688 rows_n = wt->rows;
3689 cols_n = wt->cols;
3690 lp = wt->wave->lpd_len;
3691 clen = J * 3;
3692
3693 i = 2 * J;
3694 while (i > 0) {
3695 wt->dimensions[i - 1] = cols_n;
3696 wt->dimensions[i - 2] = rows_n;
3697 wt->outlength += (rows_n * cols_n) * 3;
3698 i = i - 2;
3699 }
3700 wt->outlength += (rows_n * cols_n);
3701 N = wt->outlength;
3702 wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
3703
3704 orig = inp;
3705 ir = wt->rows;
3706 ic = wt->cols;
3707 cols_i = wt->dimensions[2 * J - 1];
3708
3709 lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3710 hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3711
3712 for (iter = 0; iter < J; ++iter) {
3713 if (iter > 0) {
3714 M = 2 * M;
3715 }
3716 rows_i = wt->dimensions[2 * J - 2 * iter - 2];
3717 cols_i = wt->dimensions[2 * J - 2 * iter - 1];
3718 istride = 1;
3719 ostride = 1;
3720 cdim = rows_i * cols_i;
3721 // Row filtering and column subsampling
3722 for (i = 0; i < ir; ++i) {
3723 swt_per_stride(M,orig + i*ic, ic, wt->wave->lpd, wt->wave->hpd, lp, lp_dn1 + i*cols_i, cols_i, hp_dn1 + i*cols_i, istride, ostride);
3724 }
3725 // Column Filtering and Row subsampling
3726 aHH = N - cdim;
3727 wt->coeffaccess[clen] = aHH;
3728 aHL = aHH - cdim;
3729 wt->coeffaccess[clen - 1] = aHL;
3730 aLH = aHL - cdim;
3731 wt->coeffaccess[clen - 2] = aLH;
3732 aLL = aLH - cdim;
3733
3734 N -= 3 * cdim;
3735 ic = cols_i;
3736 istride = ic;
3737 ostride = ic;
3738 for (i = 0; i < ic; ++i) {
3739 swt_per_stride(M,lp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aLL + i, rows_i, wavecoeff + aLH + i, istride, ostride);
3740 }
3741
3742 for (i = 0; i < ic; ++i) {
3743 swt_per_stride(M,hp_dn1 + i, ir, wt->wave->lpd, wt->wave->hpd, lp, wavecoeff + aHL + i, rows_i, wavecoeff + aHH + i, istride, ostride);
3744 }
3745
3746 ir = rows_i;
3747 orig = wavecoeff + aLL;
3748 clen -= 3;
3749 }
3750 wt->coeffaccess[0] = 0;
3751 free(lp_dn1);
3752 free(hp_dn1);
3753
3754 return wavecoeff;
3755 }
3756
iswt2(wt2_object wt,double * wavecoeffs,double * oup)3757 void iswt2(wt2_object wt, double *wavecoeffs, double *oup) {
3758 int i, k, iter, it2, it3, J, M, rows, cols, lf, ir,ic,k1,i1;
3759 double *A, *H, *V, *D,*oup1,*oup2;
3760 int aLL, aLH, aHL, aHH,shift;
3761 J = wt->J;
3762 rows = wt->rows;
3763 cols = wt->cols;
3764 lf = wt->wave->lpd_len;
3765 A = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
3766 H = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
3767 V = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
3768 D = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
3769 oup1 = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
3770 oup2 = (double*)calloc((rows + lf)*(cols + lf), sizeof(double));
3771
3772 aLL = wt->coeffaccess[0];
3773
3774 for (i = 0; i < rows; ++i) {
3775 for (k = 0; k < cols; ++k) {
3776 oup[i*cols + k] = wavecoeffs[aLL + i*cols + k];
3777 }
3778 }
3779
3780 for (iter = J; iter > 0; iter--) {
3781 aLH = wt->coeffaccess[(J - iter) * 3 + 1];
3782 aHL = wt->coeffaccess[(J - iter) * 3 + 2];
3783 aHH = wt->coeffaccess[(J - iter) * 3 + 3];
3784 M = (int)pow(2.0, (double)iter - 1);
3785
3786 for (it2 = 0; it2 < M; ++it2) {
3787 ir = 0;
3788 ic = 0;
3789 it3 = 0;
3790 // oup1
3791 for (i = it2; i < rows; i += 2 * M) {
3792 ic = 0;
3793 for (k = it2; k < cols; k += 2 * M) {
3794 A[it3] = oup[i*cols + k];
3795 H[it3] = wavecoeffs[aLH + i*cols + k];
3796 V[it3] = wavecoeffs[aHL + i*cols + k];
3797 D[it3] = wavecoeffs[aHH + i*cols + k];
3798 it3++;
3799 ic++;
3800 }
3801 ir++;
3802 }
3803 shift = 0;
3804 idwt2_shift(shift, ir, ic, wt->wave->lpr, wt->wave->hpr, wt->wave->lpd_len, A, H, V, D, oup1);
3805 //oup2
3806 ir = 0;
3807 ic = 0;
3808 it3 = 0;
3809 for (i = it2 + M; i < rows; i += 2 * M) {
3810 ic = 0;
3811 for (k = it2 + M; k < cols; k += 2 * M) {
3812 A[it3] = oup[i*cols + k];
3813 H[it3] = wavecoeffs[aLH + i*cols + k];
3814 V[it3] = wavecoeffs[aHL + i*cols + k];
3815 D[it3] = wavecoeffs[aHH + i*cols + k];
3816 it3++;
3817 ic++;
3818 }
3819 ir++;
3820 }
3821 shift = -1;
3822 idwt2_shift(shift, ir, ic, wt->wave->lpr, wt->wave->hpr, wt->wave->lpd_len, A, H, V, D, oup2);
3823 // Shift oup1 and oup2. Then add them to get A.
3824 i1 = 0;
3825 for (i = it2; i < rows; i += M) {
3826 k1 = 0;
3827 for (k = it2; k < cols; k += M) {
3828 oup[i*cols+k] = 0.5*(oup1[i1*2*ic + k1] + oup2[i1*2*ic+k1]);
3829 k1++;
3830 }
3831 i1++;
3832 }
3833 }
3834
3835 }
3836
3837 free(A);
3838 free(H);
3839 free(V);
3840 free(D);
3841 free(oup1);
3842 free(oup2);
3843 }
3844
modwt2(wt2_object wt,double * inp)3845 double* modwt2(wt2_object wt, double *inp) {
3846 double *wavecoeff;
3847 int i, J, iter, M, N, lp, rows_n, cols_n, rows_i, cols_i;
3848 int ir, ic, istride, ostride;
3849 int aLL, aLH, aHL, aHH, cdim, clen;
3850 double *orig, *lp_dn1, *hp_dn1,*filt;
3851 double s;
3852
3853 J = wt->J;
3854 M = 1;
3855 wt->outlength = 0;
3856
3857 rows_n = wt->rows;
3858 cols_n = wt->cols;
3859 lp = wt->wave->lpd_len;
3860 clen = J * 3;
3861
3862 i = 2 * J;
3863 while (i > 0) {
3864 wt->dimensions[i - 1] = cols_n;
3865 wt->dimensions[i - 2] = rows_n;
3866 wt->outlength += (rows_n * cols_n) * 3;
3867 i = i - 2;
3868 }
3869 wt->outlength += (rows_n * cols_n);
3870 N = wt->outlength;
3871 wavecoeff = (double*)calloc(wt->outlength, sizeof(double));
3872 filt = (double*)malloc(sizeof(double)* 2 * lp);
3873 s = sqrt(2.0);
3874 for (i = 0; i < lp; ++i) {
3875 filt[i] = wt->wave->lpd[i] / s;
3876 filt[lp + i] = wt->wave->hpd[i] / s;
3877 }
3878
3879 orig = inp;
3880 ir = wt->rows;
3881 ic = wt->cols;
3882 cols_i = wt->dimensions[2 * J - 1];
3883
3884 lp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3885 hp_dn1 = (double*)malloc(sizeof(double)*ir*cols_i);
3886
3887 for (iter = 0; iter < J; ++iter) {
3888 if (iter > 0) {
3889 M = 2 * M;
3890 }
3891 rows_i = wt->dimensions[2 * J - 2 * iter - 2];
3892 cols_i = wt->dimensions[2 * J - 2 * iter - 1];
3893 istride = 1;
3894 ostride = 1;
3895 cdim = rows_i * cols_i;
3896 // Row filtering and column subsampling
3897 for (i = 0; i < ir; ++i) {
3898 modwt_per_stride(M, orig + i*ic, ic, filt, lp, lp_dn1 + i*cols_i, cols_i, hp_dn1 + i*cols_i, istride, ostride);
3899 }
3900 // Column Filtering and Row subsampling
3901 aHH = N - cdim;
3902 wt->coeffaccess[clen] = aHH;
3903 aHL = aHH - cdim;
3904 wt->coeffaccess[clen - 1] = aHL;
3905 aLH = aHL - cdim;
3906 wt->coeffaccess[clen - 2] = aLH;
3907 aLL = aLH - cdim;
3908 N -= 3 * cdim;
3909 ic = cols_i;
3910 istride = ic;
3911 ostride = ic;
3912 for (i = 0; i < ic; ++i) {
3913 modwt_per_stride(M, lp_dn1 + i, ir, filt, lp, wavecoeff + aLL + i, rows_i, wavecoeff + aLH + i, istride, ostride);
3914 }
3915
3916
3917 for (i = 0; i < ic; ++i) {
3918 modwt_per_stride(M, hp_dn1 + i, ir, filt, lp, wavecoeff + aHL + i, rows_i, wavecoeff + aHH + i, istride, ostride);
3919 }
3920
3921
3922 ir = rows_i;
3923 orig = wavecoeff + aLL;
3924 clen -= 3;
3925 }
3926 wt->coeffaccess[0] = 0;
3927 free(lp_dn1);
3928 free(hp_dn1);
3929 free(filt);
3930 return wavecoeff;
3931 }
3932
imodwt2(wt2_object wt,double * wavecoeff,double * oup)3933 void imodwt2(wt2_object wt, double *wavecoeff, double *oup) {
3934 int i, rows, cols, M, N, ir, ic, lf;
3935 int istride, ostride, iter, J;
3936 int aLL, aLH, aHL, aHH;
3937 double *cL, *cH, *orig,*filt;
3938 double s;
3939
3940 rows = wt->rows;
3941 cols = wt->cols;
3942 J = wt->J;
3943
3944
3945 M = (int)pow(2.0, (double)J - 1.0);
3946 N = rows > cols ? rows : cols;
3947 lf = (wt->wave->lpr_len + wt->wave->hpr_len) / 2;
3948
3949 filt = (double*)malloc(sizeof(double)* 2 * lf);
3950 s = sqrt(2.0);
3951 for (i = 0; i < lf; ++i) {
3952 filt[i] = wt->wave->lpd[i] / s;
3953 filt[lf + i] = wt->wave->hpd[i] / s;
3954 }
3955
3956
3957 cL = (double*)calloc(rows*cols, sizeof(double));
3958 cH = (double*)calloc(rows*cols, sizeof(double));
3959 aLL = wt->coeffaccess[0];
3960 orig = wavecoeff + aLL;
3961 for (iter = 0; iter < J; ++iter) {
3962 if (iter > 0) {
3963 M = M / 2;
3964 }
3965 ir = wt->dimensions[2 * iter];
3966 ic = wt->dimensions[2 * iter + 1];
3967 istride = ic;
3968 ostride = ic;
3969 aLH = wt->coeffaccess[iter * 3 + 1];
3970 aHL = wt->coeffaccess[iter * 3 + 2];
3971 aHH = wt->coeffaccess[iter * 3 + 3];
3972 for (i = 0; i < ic; ++i) {
3973 imodwt_per_stride(M,orig + i, ir, wavecoeff + aLH + i, filt, lf, cL + i, istride, ostride);
3974
3975 imodwt_per_stride(M,wavecoeff + aHL + i, ir, wavecoeff + aHH + i, filt, lf, cH + i, istride, ostride);
3976
3977 }
3978
3979 istride = 1;
3980 ostride = 1;
3981
3982 for (i = 0; i < ir; ++i) {
3983 imodwt_per_stride(M,cL + i*ic, ic, cH + i*ic, filt, lf, oup+i*ic, istride, ostride);
3984
3985 }
3986
3987 orig = oup;
3988 }
3989
3990 free(cL);
3991 free(cH);
3992 free(filt);
3993 }
3994
getWT2Coeffs(wt2_object wt,double * wcoeffs,int level,char * type,int * rows,int * cols)3995 double* getWT2Coeffs(wt2_object wt,double* wcoeffs, int level,char *type, int *rows, int *cols) {
3996 int J,iter,t;
3997 double *ptr;
3998 J = wt->J;
3999 // Error Check
4000
4001 if (level > J || level < 1) {
4002 printf("Error : The data is decomposed into %d levels so the acceptable values of level are between 1 and %d", J, J);
4003 exit(-1);
4004 }
4005
4006 if (!strcmp(type, "A") && level != J) {
4007 printf("Approximation Coefficients are only available for level %d", J);
4008 exit(-1);
4009 }
4010
4011 if (!strcmp(type, "A")) {
4012 t = 0;
4013 iter = t;
4014 }
4015 else if (!strcmp(type, "H")) {
4016 t = 1;
4017 iter = t;
4018 }
4019 else if (!strcmp(type, "V")) {
4020 t = 2;
4021 iter = t;
4022 }
4023 else if (!strcmp(type, "D")) {
4024 t = 3;
4025 iter = t;
4026 }
4027 else {
4028 printf("Only four types of coefficients are accessible A, H, V and D \n");
4029 exit(-1);
4030 }
4031
4032 iter += (J - level) * 3;
4033
4034 ptr = wcoeffs+wt->coeffaccess[iter];
4035 *rows = wt->dimensions[2 * (J - level)];
4036 *cols = wt->dimensions[2 * (J - level)+1];
4037
4038 return ptr;
4039 }
4040
dispWT2Coeffs(double * A,int row,int col)4041 void dispWT2Coeffs(double *A, int row, int col) {
4042 int i, j;
4043 printf("\n MATRIX Order : %d X %d \n \n", row, col);
4044
4045 for (i = 0; i < row; i++) {
4046 printf("R%d: ", i);
4047 for (j = 0; j < col; j++) {
4048 printf("%g ", A[i*col + j]);
4049 }
4050 printf(":R%d \n", i);
4051 }
4052 }
4053
wave_summary(wave_object obj)4054 void wave_summary(wave_object obj) {
4055 int i,N;
4056 N = obj->filtlength;
4057 printf("\n");
4058 printf("Wavelet Name : %s \n",obj->wname);
4059 printf("\n");
4060 printf("Wavelet Filters \n\n");
4061 printf("lpd : [");
4062 for (i = 0; i < N-1; ++i) {
4063 printf("%g,", obj->lpd[i]);
4064 }
4065 printf("%g", obj->lpd[N-1]);
4066 printf("] \n\n");
4067 printf("hpd : [");
4068 for (i = 0; i < N-1; ++i) {
4069 printf("%g,", obj->hpd[i]);
4070 }
4071 printf("%g", obj->hpd[N - 1]);
4072 printf("] \n\n");
4073 printf("lpr : [");
4074 for (i = 0; i < N-1; ++i) {
4075 printf("%g,", obj->lpr[i]);
4076 }
4077 printf("%g", obj->lpr[N - 1]);
4078 printf("] \n\n");
4079 printf("hpr : [");
4080 for (i = 0; i < N-1; ++i) {
4081 printf("%g,", obj->hpr[i]);
4082 }
4083 printf("%g", obj->hpr[N - 1]);
4084 printf("] \n\n");
4085 }
4086
wt_summary(wt_object wt)4087 void wt_summary(wt_object wt) {
4088 int i;
4089 int J,t;
4090 J = wt->J;
4091 wave_summary(wt->wave);
4092 printf("\n");
4093 printf("Wavelet Transform : %s \n", wt->method);
4094 printf("\n");
4095 printf("Signal Extension : %s \n", wt->ext);
4096 printf("\n");
4097 printf("Convolutional Method : %s \n", wt->cmethod);
4098 printf("\n");
4099 printf("Number of Decomposition Levels %d \n", wt->J);
4100 printf("\n");
4101 printf("Length of Input Signal %d \n", wt->siglength);
4102 printf("\n");
4103 printf("Length of WT Output Vector %d \n", wt->outlength);
4104 printf("\n");
4105 printf("Wavelet Coefficients are contained in vector : %s \n", "output");
4106 printf("\n");
4107 printf("Approximation Coefficients \n");
4108 printf("Level %d Access : output[%d] Length : %d \n", J, 0, wt->length[0]);
4109 printf("\n");
4110 printf("Detail Coefficients \n");
4111 t = wt->length[0];
4112 for (i = 0; i < J; ++i) {
4113 printf("Level %d Access : output[%d] Length : %d \n", J - i,t,wt->length[i+1]);
4114 t += wt->length[i+1];
4115 }
4116 printf("\n");
4117
4118 }
4119
wtree_summary(wtree_object wt)4120 void wtree_summary(wtree_object wt) {
4121 int i,k,p2;
4122 int J,t;
4123 J = wt->J;
4124 wave_summary(wt->wave);
4125 printf("\n");
4126 printf("Wavelet Transform : %s \n", wt->method);
4127 printf("\n");
4128 printf("Signal Extension : %s \n", wt->ext);
4129 printf("\n");
4130 printf("Number of Decomposition Levels %d \n", wt->J);
4131 printf("\n");
4132 printf("Length of Input Signal %d \n", wt->siglength);
4133 printf("\n");
4134 printf("Length of WT Output Vector %d \n", wt->outlength);
4135 printf("\n");
4136 printf("Wavelet Coefficients are contained in vector : %s \n", "output");
4137 printf("\n");
4138 printf("Coefficients Access \n");
4139 t = 0;
4140 p2 = 2;
4141 for (i = 0; i < J; ++i) {
4142 for (k = 0; k < p2; ++k) {
4143 printf("Node %d %d Access : output[%d] Length : %d \n", i + 1, k, wt->nodelength[t], wt->length[J - i]);
4144 t++;
4145 }
4146 p2 *= 2;
4147 }
4148 printf("\n");
4149
4150 }
4151
wpt_summary(wpt_object wt)4152 void wpt_summary(wpt_object wt) {
4153 int i, k, p2;
4154 int J, it1,it2;
4155 J = wt->J;
4156 wave_summary(wt->wave);
4157 printf("\n");
4158 printf("Signal Extension : %s \n", wt->ext);
4159 printf("\n");
4160 printf("Entropy : %s \n", wt->entropy);
4161 printf("\n");
4162 printf("Number of Decomposition Levels %d \n", wt->J);
4163 printf("\n");
4164 printf("Number of Active Nodes %d \n", wt->nodes);
4165 printf("\n");
4166 printf("Length of Input Signal %d \n", wt->siglength);
4167 printf("\n");
4168 printf("Length of WT Output Vector %d \n", wt->outlength);
4169 printf("\n");
4170 printf("Wavelet Coefficients are contained in vector : %s \n", "output");
4171 printf("\n");
4172 printf("Coefficients Access \n");
4173 it1 = 1;
4174 it2 = 0;
4175 for (i = 0; i < J; ++i) {
4176 it1 += ipow2(i + 1);
4177 }
4178 for (i = J; i > 0; --i) {
4179 p2 = ipow2(i);
4180 it1 -= p2;
4181 for (k = 0; k < p2; ++k) {
4182 if (wt->basisvector[it1 + k] == 1) {
4183 printf("Node %d %d Access : output[%d] Length : %d \n", i, k, it2, wt->length[J - i + 1]);
4184 it2 += wt->length[J - i + 1];
4185 }
4186 }
4187 }
4188
4189 printf("\n");
4190
4191 }
4192
cwt_summary(cwt_object wt)4193 void cwt_summary(cwt_object wt) {
4194
4195 printf("\n");
4196 printf("Wavelet : %s Parameter %lf \n", wt->wave,wt->m);
4197 printf("\n");
4198 printf("Length of Input Signal : %d \n", wt->siglength);
4199 printf("\n");
4200 printf("Sampling Rate : %g \n", wt->dt);
4201 printf("\n");
4202 printf("Total Number of Scales : %d \n", wt->J);
4203 printf("\n");
4204 printf("Smallest Scale (s0) : %lf \n", wt->s0);
4205 printf("\n");
4206 printf("Separation Between Scales (dj) %lf \n", wt->dj);
4207 printf("\n");
4208 printf("Scale Type %s \n", wt->type);
4209 printf("\n");
4210 printf("Complex CWT Output Vector is of size %d * %d stored in Row Major format \n",wt->J,wt->siglength);
4211 printf("\n");
4212 printf("The ith real value can be accessed using wt->output[i].re and imaginary value by wt->output[i].im \n");
4213 printf("\n");
4214
4215 }
4216
wt2_summary(wt2_object wt)4217 void wt2_summary(wt2_object wt) {
4218 int i;
4219 int J, t,rows,cols,vsize;
4220 J = wt->J;
4221 wave_summary(wt->wave);
4222 printf("\n");
4223 printf("Wavelet Transform : %s \n", wt->method);
4224 printf("\n");
4225 printf("Signal Extension : %s \n", wt->ext);
4226 printf("\n");
4227 printf("Number of Decomposition Levels %d \n", wt->J);
4228 printf("\n");
4229 printf("Input Signal Rows %d \n", wt->rows);
4230 printf("\n");
4231 printf("Input Signal Cols %d \n", wt->cols);
4232 printf("\n");
4233 printf("Length of Wavelet Coefficients Vector %d \n", wt->outlength);
4234 printf("\n");
4235 t = 0;
4236 for (i = J; i > 0; --i) {
4237 rows = wt->dimensions[2 * (J - i)];
4238 cols = wt->dimensions[2 * (J - i) + 1];
4239 vsize = rows*cols;
4240 printf("Level %d Decomposition Rows :%d Columns:%d Vector Size (Rows*Cols):%d \n", i, rows, cols,vsize);
4241 printf("Access Row values stored at wt->dimensions[%d]\n", 2 * (J - i));
4242 printf("Access Column values stored at wt->dimensions[%d]\n\n", 2 * (J - i) + 1);
4243
4244 if (i == J) {
4245 printf("Approximation Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n", t,wt->coeffaccess[t],vsize);
4246 }
4247
4248 t += 1;
4249 printf("Horizontal Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n", t, wt->coeffaccess[t], vsize);
4250 t += 1;
4251 printf("Vertical Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n", t, wt->coeffaccess[t], vsize);
4252 t += 1;
4253 printf("Diagonal Coefficients access at wt->coeffaccess[%d]=%d, Vector size:%d \n\n", t, wt->coeffaccess[t], vsize);
4254 }
4255
4256 }
4257
wave_free(wave_object object)4258 void wave_free(wave_object object) {
4259 free(object);
4260 }
4261
wt_free(wt_object object)4262 void wt_free(wt_object object) {
4263 free(object);
4264 }
4265
wtree_free(wtree_object object)4266 void wtree_free(wtree_object object) {
4267 free(object);
4268 }
4269
wpt_free(wpt_object object)4270 void wpt_free(wpt_object object) {
4271 free(object);
4272 }
4273
cwt_free(cwt_object object)4274 void cwt_free(cwt_object object) {
4275 free(object);
4276 }
4277
wt2_free(wt2_object wt)4278 void wt2_free(wt2_object wt) {
4279 free(wt);
4280 }
4281