1 /*
2 * -------------------------------------------------------------------------
3 * swt.c -- 1-D and 2-D discrete stationary wavelet transform
4 * SWT - Scilab wavelet toolbox
5 * Copyright (C) 2005-2006 Roger Liu
6 * Copyright (C) 20010-2012 Holger Nahrstaedt
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 * -------------------------------------------------------------------------
22 */
23 #include "swtlib.h"
24
swt_conv(double * sigIn,int sigInLength,double * approx,int approxLength,double * detail,int detailLength,double * filterLow,double * filterHi,int filterLength)25 void swt_conv(double *sigIn, int sigInLength,
26 double *approx, int approxLength,
27 double *detail, int detailLength,
28 double *filterLow, double *filterHi,
29 int filterLength)
30 {
31 i_conv(sigIn,sigInLength,approx,sigInLength,filterLow,filterLength);
32 i_conv(sigIn,sigInLength,detail,sigInLength,filterHi,filterLength);
33 return;
34 }
35
swt_out1(double * sigIn,int sigInLength,double * sigOutMatrix,int rowLength,int colLength,double * filterLow,double * filterHi,int filterLength,int step)36 void swt_out1 (double *sigIn, int sigInLength,
37 double *sigOutMatrix, int rowLength, int colLength,
38 double *filterLow, double *filterHi, int filterLength, int step)
39 {
40 int row, count, *filterLen;
41 double *approxTemp,*approxOut,*sigOutTemp;
42 double *filterLowTemp, *filterHiTemp, *filterLowM, *filterHiM;
43
44 filterLen = malloc(step*sizeof(int));
45 approxTemp = malloc(colLength * sizeof(double));
46 approxOut = malloc(colLength*sizeof(double));
47 sigOutTemp = malloc(rowLength*colLength*sizeof(double));
48
49 filterLen[0] = filterLength;
50 if (step>1)
51 {
52 for(count=1;count<step;count++)
53 filterLen[count] = filterLen[count-1] * 2;
54 }
55
56 filterLowTemp = malloc((filterLen[step-1]+1)*sizeof(double));
57 filterHiTemp = malloc((filterLen[step-1]+1)*sizeof(double));
58 filterLowM = malloc((filterLen[step-1]+1)*sizeof(double));
59 filterHiM = malloc((filterLen[step-1]+1)*sizeof(double));
60 for(count=0;count<(filterLen[step-1]+1);count++)
61 {
62 filterLowM[count] = 0;
63 filterHiM[count] = 0;
64 }
65
66 verbatim_copy(filterLow,filterLength,filterLowTemp,filterLen[0]);
67 verbatim_copy(filterHi,filterLength,filterHiTemp,filterLen[0]);
68 verbatim_copy(sigIn,sigInLength,approxTemp,sigInLength);
69 for (row=0;row<step;row++)
70 {
71 swt_conv(approxTemp,sigInLength,
72 approxOut,sigInLength,
73 sigOutTemp+row*colLength,sigInLength,
74 filterLowTemp,filterHiTemp,filterLen[row]);
75 verbatim_copy(approxOut,sigInLength,approxTemp,sigInLength);
76 if (row<(step-1))
77 {
78 dyadup_1D_feed_even (filterLowTemp, filterLen[row],
79 filterLowM, filterLen[row]*2+1);
80 dyadup_1D_feed_even (filterHiTemp, filterLen[row],
81 filterHiM, filterLen[row]*2+1);
82 verbatim_copy(filterLowM+1,filterLen[row]*2,filterLowTemp,filterLen[row]*2);
83 verbatim_copy(filterHiM+1,filterLen[row]*2,filterHiTemp,filterLen[row]*2);
84 }
85 }
86 verbatim_copy(approxTemp,sigInLength,sigOutTemp+step*colLength,sigInLength);
87 free(approxTemp);
88 free(approxOut);
89 free(filterLowM);
90 free(filterHiM);
91 free(filterLowTemp);
92 free(filterHiTemp);
93 free(filterLen);
94 matrix_tran (sigOutTemp, rowLength, colLength,
95 sigOutMatrix, colLength, rowLength);
96 free(sigOutTemp);
97 return;
98 }
99
swt_out2(double * sigIn,int sigInLength,double * approxMatrix,double * detailMatrix,int rowLength,int colLength,double * filterLow,double * filterHi,int filterLength,int step)100 void swt_out2 (double *sigIn, int sigInLength,
101 double *approxMatrix, double *detailMatrix,
102 int rowLength, int colLength, double *filterLow,
103 double *filterHi, int filterLength, int step)
104 {
105 int row, count, *filterLen;
106 double *approxMatrixTemp, *detailMatrixTemp, *approxTemp;
107 double *filterLowTemp, *filterHiTemp, *filterLowM, *filterHiM;
108
109
110 filterLen = malloc(step*sizeof(int));
111 approxTemp = malloc(colLength * sizeof(double));
112 approxMatrixTemp = malloc(rowLength*colLength*sizeof(double));
113 detailMatrixTemp = malloc(rowLength*colLength*sizeof(double));
114
115 filterLen[0] = filterLength;
116 if (step>1)
117 {
118 for(count=1;count<step;count++)
119 filterLen[count] = filterLen[count-1] * 2;
120 }
121
122 filterLowTemp = malloc((filterLen[step-1])*sizeof(double));
123 filterHiTemp = malloc((filterLen[step-1])*sizeof(double));
124 filterLowM = malloc((filterLen[step-1])*sizeof(double));
125 filterHiM = malloc((filterLen[step-1])*sizeof(double));
126 for(count=0;count<(filterLen[step-1]);count++)
127 {
128 filterLowM[count] = 0;
129 filterHiM[count] = 0;
130 }
131
132 verbatim_copy(filterLow,filterLength,filterLowTemp,filterLen[0]);
133 verbatim_copy(filterHi,filterLength,filterHiTemp,filterLen[0]);
134 verbatim_copy(sigIn,sigInLength,approxTemp,sigInLength);
135 for (row=0;row<step;row++)
136 {
137 swt_conv(approxTemp,sigInLength,
138 approxMatrixTemp+row*colLength,sigInLength,
139 detailMatrixTemp+row*colLength,sigInLength,
140 filterLowTemp,filterHiTemp,filterLen[row]);
141 verbatim_copy(approxMatrixTemp+row*colLength,sigInLength,approxTemp,sigInLength);
142 if (row<(step-1))
143 {
144 for(count=0;count<filterLen[row];count++)
145 {
146 filterLowM[2*count] = filterLowTemp[count];
147 filterLowM[2*count+1] = 0;
148 filterHiM[2*count] = filterHiTemp[count];
149 filterHiM[2*count+1] = 0;
150 }
151
152 verbatim_copy(filterLowM,filterLen[row]*2,filterLowTemp,filterLen[row]*2);
153 verbatim_copy(filterHiM,filterLen[row]*2,filterHiTemp,filterLen[row]*2);
154 }
155 }
156 free(approxTemp);
157 matrix_tran (approxMatrixTemp, rowLength, colLength,
158 approxMatrix, colLength, rowLength);
159 matrix_tran (detailMatrixTemp, rowLength, colLength,
160 detailMatrix, colLength, rowLength);
161 free(approxMatrixTemp);
162 free(detailMatrixTemp);
163 free(filterLowM);
164 free(filterHiM);
165 free(filterLowTemp);
166 free(filterHiTemp);
167 free(filterLen);
168 return;
169 }
170
swt2_output4(double * matrixIn,int matrixInRow,int matrixInCol,double * matrixOutApprox,double * matrixOutColDetail,double * matrixOutRowDetail,double * matrixOutDetail,int matrixOutRow,int matrixOutCol,double * filterLow,double * filterHi,int filterLength,int step)171 void swt2_output4(double *matrixIn, int matrixInRow, int matrixInCol,
172 double *matrixOutApprox, double *matrixOutColDetail,
173 double *matrixOutRowDetail, double *matrixOutDetail,
174 int matrixOutRow, int matrixOutCol,
175 double *filterLow, double *filterHi,
176 int filterLength, int step)
177 {
178 int row, col, count;
179 int filterLen, twt;
180 double *matrixOutApproxPre, *matrixOutDetailPre;
181 double *filterLowM, *filterHiM;
182 double *pMatrixOutApproxPre, *pMatrixOutDetailPre;
183
184 //matrixInPre = malloc(matrixInRow*matrixInCol*sizeof(double));
185 //matrix_tran (matrixIn, matrixInCol, matrixInRow,
186 //matrixInPre, matrixInRow, matrixInCol);
187
188 matrixOutApproxPre = malloc (matrixInRow *
189 matrixOutCol * sizeof (double));
190 matrixOutDetailPre = malloc (matrixInRow *
191 matrixOutCol * sizeof (double));
192
193 swt_exp2(step,&twt);
194 filterLen = filterLength*twt;
195 filterLowM = malloc(filterLen*sizeof(double));
196 filterHiM = malloc(filterLen*sizeof(double));
197 for(count=0;count<filterLen;count++)
198 {
199 if (count%twt==0)
200 {
201 filterLowM[count] = filterLow[count/twt];
202 filterHiM[count] = filterHi[count/twt];
203 }
204 else
205 {
206 filterLowM[count] = 0;
207 filterHiM[count] = 0;
208 }
209 }
210
211
212 for (row=0;row<matrixInRow;row++)
213 {
214 swt_conv(matrixIn+row*matrixInCol,matrixInCol,
215 matrixOutApproxPre+row*matrixInCol,matrixInCol,
216 matrixOutDetailPre+row*matrixInCol,matrixInCol,
217 filterLowM,filterHiM,filterLen);
218
219 }
220
221 //free (matrixInPre);
222 pMatrixOutApproxPre = malloc (matrixInRow *
223 matrixOutCol * sizeof (double));
224 matrix_tran (matrixOutApproxPre, matrixInRow, matrixInCol,
225 pMatrixOutApproxPre, matrixInRow, matrixInCol);
226 free (matrixOutApproxPre);
227
228 pMatrixOutDetailPre = malloc (matrixInRow *
229 matrixOutCol * sizeof (double));
230 matrix_tran (matrixOutDetailPre, matrixInRow, matrixInCol,
231 pMatrixOutDetailPre, matrixInRow, matrixInCol);
232
233 free (matrixOutDetailPre);
234
235 for (col = 0; col < matrixInCol; col++)
236 {
237 swt_conv(pMatrixOutApproxPre+col*matrixInRow,matrixInRow,
238 matrixOutApprox+col*matrixInRow,matrixInRow,
239 matrixOutColDetail+col*matrixInRow,matrixInRow,
240 filterLowM,filterHiM,filterLen);
241 }
242 free (pMatrixOutApproxPre);
243
244 for (col = 0; col < matrixInCol; col++)
245 {
246 swt_conv(pMatrixOutDetailPre+col*matrixInRow,matrixInRow,
247 matrixOutRowDetail+col*matrixInRow,matrixInRow,
248 matrixOutDetail+col*matrixInRow,matrixInRow,
249 filterLowM,filterHiM,filterLen);
250 }
251 free (pMatrixOutDetailPre);
252 return;
253 }
254
swt2_output4_step(double * matrixIn,int matrixInRow,int matrixInCol,double * matrixOutApprox,double * matrixOutColDetail,double * matrixOutRowDetail,double * matrixOutDetail,int matrixOutRow,int matrixOutCol,double * filterLow,double * filterHi,int filterLength,int step)255 void swt2_output4_step(double *matrixIn, int matrixInRow, int matrixInCol,
256 double *matrixOutApprox, double *matrixOutColDetail,
257 double *matrixOutRowDetail, double *matrixOutDetail,
258 int matrixOutRow, int matrixOutCol,
259 double *filterLow, double *filterHi,
260 int filterLength, int step)
261 {
262 int count;
263 double *matrixApproxTemp,*matrixApproxOut;
264
265 matrixApproxTemp = malloc(matrixInRow*matrixInCol*sizeof(double));
266 matrixApproxOut = malloc(matrixInRow*matrixInCol*sizeof(double));
267 matrix_tran (matrixIn, matrixInCol, matrixInRow,
268 matrixApproxTemp, matrixInRow, matrixInCol);
269
270 for(count=0;count<step;count++)
271 {
272 swt2_output4(matrixApproxTemp, matrixInRow, matrixInCol,
273 matrixApproxOut, matrixOutColDetail+count*matrixInRow*matrixInCol,
274 matrixOutRowDetail+count*matrixInRow*matrixInCol,
275 matrixOutDetail+count*matrixInRow*matrixInCol,
276 matrixOutRow, matrixOutCol, filterLow, filterHi,
277 filterLength, count);
278 verbatim_copy(matrixApproxOut,matrixInRow*matrixInCol,
279 matrixOutApprox+count*matrixInRow*matrixInCol,
280 matrixInRow*matrixInCol);
281 matrix_tran (matrixApproxOut, matrixInCol, matrixInRow,
282 matrixApproxTemp, matrixInRow, matrixInCol);
283 }
284 free(matrixApproxOut);
285 free(matrixApproxTemp);
286
287 return;
288 }
289
swt2_output1_step(double * matrixIn,int matrixInRow,int matrixInCol,double * matrixOut,int matrixOutRow,int matrixOutCol,double * filterLow,double * filterHi,int filterLength,int step)290 void swt2_output1_step(double *matrixIn, int matrixInRow,
291 int matrixInCol, double *matrixOut,
292 int matrixOutRow, int matrixOutCol,
293 double *filterLow, double *filterHi,
294 int filterLength, int step)
295 {
296 int count;
297 double *matrixApproxTemp,*matrixApproxOut;
298
299 matrixApproxTemp = malloc(matrixInRow*matrixInCol*sizeof(double));
300 matrixApproxOut = malloc(matrixInRow*matrixInCol*sizeof(double));
301 matrix_tran (matrixIn, matrixInCol, matrixInRow,
302 matrixApproxTemp, matrixInRow, matrixInCol);
303
304 for(count=0;count<step;count++)
305 {
306 swt2_output4(matrixApproxTemp, matrixInRow, matrixInCol,
307 matrixApproxOut, matrixOut+count*matrixInRow*matrixInCol,
308 matrixOut+step*matrixInRow*matrixInCol+count*matrixInRow*matrixInCol,
309 matrixOut+2*step*matrixInRow*matrixInCol+count*matrixInRow*matrixInCol,
310 matrixInRow, matrixOutCol, filterLow, filterHi,
311 filterLength, count);
312
313 matrix_tran (matrixApproxOut, matrixInCol, matrixInRow,
314 matrixApproxTemp, matrixInRow, matrixInCol);
315 }
316 verbatim_copy(matrixApproxOut,matrixInRow*matrixInCol,
317 matrixOut+3*step*matrixInRow*matrixInCol,
318 matrixInRow*matrixInCol);
319 free(matrixApproxOut);
320 free(matrixApproxTemp);
321
322 return;
323 }
324
iswt_conv(double * approx,double * detail,int sigInLength,double * sigOut,int sigOutLength,double * filterLow,double * filterHi,int filterLength)325 void iswt_conv (double *approx, double *detail, int sigInLength,
326 double *sigOut, int sigOutLength, double *filterLow,
327 double *filterHi, int filterLength)
328 {
329 int length1, length2, count;
330 double *approxTempOdd, *detailTempOdd;
331 double *approxTempEven, *detailTempEven;
332 double *approxOutOdd, *detailOutOdd;
333 double *approxOutEven, *detailOutEven;
334 double *approxT, *detailT;
335 double *sigOdd, *sigEven;
336
337 length1 = (int)floor(sigInLength/2.0);
338 approxTempOdd = malloc(length1*sizeof(double));
339 detailTempOdd = malloc(length1*sizeof(double));
340 approxTempEven = malloc(length1*sizeof(double));
341 detailTempEven = malloc(length1*sizeof(double));
342
343 dyaddown_1D_keep_odd (approx, sigInLength, approxTempOdd, length1);
344 dyaddown_1D_keep_even (approx, sigInLength, approxTempEven, length1);
345 dyaddown_1D_keep_odd (detail, sigInLength, detailTempOdd, length1);
346 dyaddown_1D_keep_even (detail, sigInLength, detailTempEven, length1);
347
348 length2 = 2 * length1;
349 approxOutOdd = malloc(length2*sizeof(double));
350 approxOutEven = malloc(length2*sizeof(double));
351 detailOutOdd = malloc(length2*sizeof(double));
352 detailOutEven = malloc(length2*sizeof(double));
353
354 for (count=0;count<length1;count++)
355 {
356 approxOutOdd[2*count] = approxTempOdd[count];
357 detailOutOdd[2*count] = detailTempOdd[count];
358 approxOutOdd[2*count+1] = 0;
359 detailOutOdd[2*count+1] = 0;
360 approxOutEven[2*count] = 0;
361 detailOutEven[2*count] = 0;
362 approxOutEven[2*count+1] = approxTempEven[count];
363 detailOutEven[2*count+1] = detailTempEven[count];
364 }
365
366 free(approxTempOdd);
367 free(detailTempOdd);
368 free(approxTempEven);
369 free(detailTempEven);
370
371 approxT = malloc(length2*sizeof(double));
372 detailT = malloc(length2*sizeof(double));
373 sigOdd = malloc(length2*sizeof(double));
374 sigEven = malloc(length2*sizeof(double));
375
376 i_conv(approxOutOdd,length2,approxT,length2,filterLow,filterLength);
377 i_conv(detailOutOdd,length2,detailT,length2,filterHi,filterLength);
378 for(count=0;count<length2;count++)
379 sigOdd[count] = approxT[count] + detailT[count];
380 free(approxOutOdd);
381 free(detailOutOdd);
382
383 i_conv(approxOutEven,length2,approxT,length2,filterLow,filterLength);
384 i_conv(detailOutEven,length2,detailT,length2,filterHi,filterLength);
385 for(count=0;count<length2;count++)
386 sigEven[count] = approxT[count] + detailT[count];
387 free(approxOutEven);
388 free(detailOutEven);
389 free(approxT);
390 free(detailT);
391
392 for(count=(sigOutLength-filterLength-1);count<sigOutLength;count++)
393 sigOut[count+1+filterLength-sigOutLength] = (sigOdd[count] + sigEven[count]) / 2.0;
394 for(count=0;count<(sigOutLength - filterLength - 1);count++)
395 sigOut[count+filterLength+1] = (sigOdd[count]+sigEven[count])/2.0;
396 free(sigOdd);
397 free(sigEven);
398
399 return;
400 }
401
iswt_conv_step(double * approx,double * detail,int sigInLength,double * sigOut,int sigOutLength,double * filterLow,double * filterHi,int filterLength,int level)402 void iswt_conv_step (double *approx, double *detail, int sigInLength,
403 double *sigOut, int sigOutLength, double *filterLow,
404 double *filterHi, int filterLength, int level)
405 {
406 int length1, length2, count, twt, mu, su, startpoint;
407 double *approxTempOdd, *detailTempOdd;
408 double *approxTempEven, *detailTempEven;
409 double *approxOutOdd, *detailOutOdd;
410 double *approxOutEven, *detailOutEven;
411 double *approxT, *detailT;
412 double *sigOdd, *sigEven;
413 double *filterLowTemp, *filterHiTemp;
414
415 length1 = (int)floor(sigInLength/2.0);
416 approxTempOdd = malloc(length1*sizeof(double));
417 detailTempOdd = malloc(length1*sizeof(double));
418 approxTempEven = malloc(length1*sizeof(double));
419 detailTempEven = malloc(length1*sizeof(double));
420
421 dyaddown_1D_keep_odd (approx, sigInLength, approxTempOdd, length1);
422 dyaddown_1D_keep_even (approx, sigInLength, approxTempEven, length1);
423 dyaddown_1D_keep_odd (detail, sigInLength, detailTempOdd, length1);
424 dyaddown_1D_keep_even (detail, sigInLength, detailTempEven, length1);
425
426 length2 = 2 * length1;
427 approxOutOdd = malloc(length2*sizeof(double));
428 approxOutEven = malloc(length2*sizeof(double));
429 detailOutOdd = malloc(length2*sizeof(double));
430 detailOutEven = malloc(length2*sizeof(double));
431
432 for (count=0;count<length1;count++)
433 {
434 approxOutOdd[2*count] = approxTempOdd[count];
435 detailOutOdd[2*count] = detailTempOdd[count];
436 approxOutOdd[2*count+1] = 0;
437 detailOutOdd[2*count+1] = 0;
438 approxOutEven[2*count] = 0;
439 detailOutEven[2*count] = 0;
440 approxOutEven[2*count+1] = approxTempEven[count];
441 detailOutEven[2*count+1] = detailTempEven[count];
442 }
443
444 free(approxTempOdd);
445 free(detailTempOdd);
446 free(approxTempEven);
447 free(detailTempEven);
448
449 swt_exp2(level-1,&twt);
450 if (level==1)
451 {
452 mu = 1;
453 su = 0;
454 }
455 else
456 {
457 mu = twt;
458 su = twt - 1;
459 }
460
461 filterLowTemp = malloc(mu*filterLength*sizeof(double));
462 filterHiTemp = malloc(mu*filterLength*sizeof(double));
463
464 for(count=0;count<mu*filterLength;count++)
465 {
466 filterLowTemp[count] = 0;
467 filterHiTemp[count] = 0;
468 }
469
470 for(count=0;count<filterLength;count++)
471 {
472 filterLowTemp[mu*count] = filterLow[count];
473 filterHiTemp[mu*count] = filterHi[count];
474 }
475
476 approxT = malloc(length2*sizeof(double));
477 detailT = malloc(length2*sizeof(double));
478 sigOdd = malloc(length2*sizeof(double));
479 sigEven = malloc(length2*sizeof(double));
480
481 i_conv(approxOutOdd,length2,approxT,length2,filterLowTemp,mu*filterLength);
482 i_conv(detailOutOdd,length2,detailT,length2,filterHiTemp,mu*filterLength);
483 for(count=0;count<length2;count++)
484 sigOdd[count] = approxT[count] + detailT[count];
485 free(approxOutOdd);
486 free(detailOutOdd);
487
488 i_conv(approxOutEven,length2,approxT,length2,filterLowTemp,mu*filterLength);
489 i_conv(detailOutEven,length2,detailT,length2,filterHiTemp,mu*filterLength);
490 for(count=0;count<length2;count++)
491 sigEven[count] = approxT[count] + detailT[count];
492 free(approxOutEven);
493 free(detailOutEven);
494 free(approxT);
495 free(detailT);
496 free(filterLowTemp);
497 free(filterHiTemp);
498
499
500 startpoint = sigInLength - filterLength * mu - su - 1;
501
502 for(count=(startpoint);count<sigOutLength;count++)
503 sigOut[count-startpoint] = (sigOdd[count] + sigEven[count]) / 2.0;
504 for(count=0;count<startpoint;count++)
505 sigOut[count+1+filterLength*mu+su] = (sigOdd[count]+sigEven[count])/2.0;
506 free(sigOdd);
507 free(sigEven);
508
509 return;
510 }
511
iswt_input1(double * matrixIn,int rowLength,int colLength,double * sigOut,int sigOutLength,double * filterLow,double * filterHi,int filterLength)512 void iswt_input1 (double *matrixIn, int rowLength, int colLength,
513 double *sigOut, int sigOutLength, double *filterLow,
514 double *filterHi, int filterLength)
515 {
516 int step, count;
517 double *approxTemp, *approxOut, *matrixInTemp;
518
519 matrixInTemp = malloc(rowLength*colLength*sizeof(double));
520 matrix_tran(matrixIn, colLength, rowLength,
521 matrixInTemp, colLength, rowLength);
522
523 step = rowLength - 1;
524 approxTemp = malloc(colLength*sizeof(double));
525 approxOut = malloc(colLength*sizeof(double));
526
527
528 verbatim_copy(matrixInTemp+step*colLength,colLength,approxTemp,colLength);
529
530 for(count=0;count<step;count++)
531 {
532 iswt_conv_step(approxTemp,matrixInTemp+(step-count-1)*colLength,colLength,
533 approxOut,colLength,filterLow,filterHi,filterLength,step-count);
534 verbatim_copy(approxOut,colLength,approxTemp,colLength);
535 }
536 verbatim_copy(approxTemp,colLength,sigOut, colLength);
537 free(matrixInTemp);
538 free(approxTemp);
539 free(approxOut);
540
541 return;
542 }
543
544
545
iswt_input2(double * matrixApproxIn,double * matrixDetailIn,int rowLength,int colLength,double * sigOut,int sigOutLength,double * filterLow,double * filterHi,int filterLength)546 void iswt_input2 (double *matrixApproxIn, double *matrixDetailIn,
547 int rowLength, int colLength,
548 double *sigOut, int sigOutLength, double *filterLow,
549 double *filterHi, int filterLength)
550 {
551 int count, step;
552 double *approxTemp, *approxOut;
553 double *matrixApproxInTemp, *matrixDetailInTemp;
554
555 step = rowLength;
556 approxTemp = malloc(colLength*sizeof(double));
557 approxOut = malloc(colLength*sizeof(double));
558 matrixApproxInTemp = malloc(rowLength*colLength*sizeof(double));
559 matrixDetailInTemp = malloc(rowLength*colLength*sizeof(double));
560
561 matrix_tran(matrixApproxIn, colLength, rowLength,
562 matrixApproxInTemp,colLength, rowLength);
563 matrix_tran(matrixDetailIn, colLength, rowLength,
564 matrixDetailInTemp,colLength, rowLength);
565
566 verbatim_copy(matrixApproxInTemp+(step-1)*colLength,colLength,approxTemp,colLength);
567 for(count=0;count<step;count++)
568 {
569 iswt_conv_step(approxTemp,matrixDetailInTemp+(step-count-1)*colLength,colLength,
570 approxOut,colLength,filterLow,filterHi,filterLength,step-count);
571 verbatim_copy(approxOut,colLength,approxTemp,colLength);
572 }
573 verbatim_copy(approxTemp,colLength,sigOut, colLength);
574 free(approxTemp);
575 free(approxOut);
576 free(matrixApproxInTemp);
577 free(matrixDetailInTemp);
578
579 return;
580 }
581
582
iswt2(double * matrixInApprox,double * matrixInColDetail,double * matrixInRowDetail,double * matrixInDetail,int matrixInRow,int matrixInCol,double * matrixOut,int matrixOutRow,int matrixOutCol,double * filterLow,double * filterHi,int filterLength,int step)583 void iswt2(double *matrixInApprox, double *matrixInColDetail,
584 double *matrixInRowDetail, double *matrixInDetail,
585 int matrixInRow, int matrixInCol,
586 double *matrixOut, int matrixOutRow, int matrixOutCol,
587 double *filterLow, double *filterHi,
588 int filterLength, int step)
589 {
590 int row, col;//, count;
591 //int filterLen, twt;
592 double *matrixOutApproxPre, *matrixOutDetailPre;
593 double *matrixOutPre;
594 double *pMatrixOutApproxPre, *pMatrixOutDetailPre;
595
596
597 matrixOutApproxPre = malloc (matrixInRow *
598 matrixInCol * sizeof (double));
599 matrixOutDetailPre = malloc (matrixInRow *
600 matrixInCol * sizeof (double));
601
602
603 for (col=0;col<matrixInCol;col++)
604 {
605 iswt_conv_step (matrixInApprox+col*matrixInRow, matrixInColDetail+col*matrixInRow,
606 matrixInRow, matrixOutApproxPre+col*matrixInRow, matrixInRow,
607 filterLow, filterHi, filterLength, step);
608 iswt_conv_step (matrixInRowDetail+col*matrixInRow, matrixInDetail+col*matrixInRow,
609 matrixInRow, matrixOutDetailPre+col*matrixInRow, matrixInRow,
610 filterLow, filterHi, filterLength, step);
611 }
612
613
614 pMatrixOutApproxPre = malloc (matrixInRow *
615 matrixInCol * sizeof (double));
616 matrix_tran (matrixOutApproxPre, matrixInCol, matrixInRow,
617 pMatrixOutApproxPre, matrixInRow, matrixInCol);
618 free (matrixOutApproxPre);
619
620 pMatrixOutDetailPre = malloc (matrixInRow *
621 matrixInCol * sizeof (double));
622 matrix_tran (matrixOutDetailPre, matrixInCol, matrixInRow,
623 pMatrixOutDetailPre, matrixInRow, matrixInCol);
624
625 free (matrixOutDetailPre);
626
627 matrixOutPre = malloc(matrixInRow*matrixInCol*sizeof(double));
628 for (row = 0; row < matrixInRow; row++)
629 {
630 iswt_conv_step (pMatrixOutApproxPre+row*matrixInCol, pMatrixOutDetailPre+row*matrixInCol,
631 matrixInCol, matrixOutPre+row*matrixInCol, matrixInCol,
632 filterLow, filterHi, filterLength, step);
633 }
634 free (pMatrixOutApproxPre);
635 free (pMatrixOutDetailPre);
636
637 matrix_tran (matrixOutPre, matrixInRow, matrixInCol,
638 matrixOut, matrixInRow, matrixInCol);
639 free(matrixOutPre);
640
641 return;
642 }
643
644
iswt2_input4_step(double * matrixInApprox,double * matrixInColDetail,double * matrixInRowDetail,double * matrixInDetail,int matrixInRow,int matrixInCol,double * matrixOut,int matrixOutRow,int matrixOutCol,double * filterLow,double * filterHi,int filterLength,int step)645 void iswt2_input4_step(double *matrixInApprox, double *matrixInColDetail,
646 double *matrixInRowDetail, double *matrixInDetail,
647 int matrixInRow, int matrixInCol,
648 double *matrixOut, int matrixOutRow, int matrixOutCol,
649 double *filterLow, double *filterHi,
650 int filterLength, int step)
651 {
652 int count;
653 double *matrixApproxTemp, *matrixApproxOut;
654
655 matrixApproxTemp = malloc(matrixInRow*matrixInCol*sizeof(double));
656 matrixApproxOut = malloc(matrixInRow*matrixInCol*sizeof(double));
657
658 verbatim_copy(matrixInApprox+(step-1)*matrixInRow*matrixInCol,
659 matrixInRow*matrixInCol,
660 matrixApproxTemp, matrixInRow*matrixInCol);
661
662 for(count=0;count<step;count++)
663 {
664 iswt2(matrixApproxTemp, matrixInColDetail+(step-count-1)*matrixInRow*matrixInCol,
665 matrixInRowDetail+(step-count-1)*matrixInRow*matrixInCol,
666 matrixInDetail+(step-count-1)*matrixInRow*matrixInCol,
667 matrixInRow, matrixInCol, matrixApproxOut, matrixOutRow, matrixOutCol,
668 filterLow, filterHi, filterLength, step-count);
669 if (count!=(step-1))
670 {
671 verbatim_copy(matrixApproxOut, matrixInRow*matrixInCol,
672 matrixApproxTemp, matrixInRow*matrixInCol);
673 }
674 }
675 free(matrixApproxTemp);
676
677 verbatim_copy(matrixApproxOut, matrixInRow*matrixInCol,
678 matrixOut, matrixInRow*matrixInCol);
679 free(matrixApproxOut);
680
681 return;
682 }
683
iswt2_input1_step(double * matrixIn,int matrixInRow,int matrixInCol,double * matrixOut,int matrixOutRow,int matrixOutCol,double * filterLow,double * filterHi,int filterLength,int step)684 void iswt2_input1_step(double *matrixIn, int matrixInRow, int matrixInCol,
685 double *matrixOut, int matrixOutRow, int matrixOutCol,
686 double *filterLow, double *filterHi,
687 int filterLength, int step)
688 {
689 int count;
690 double *matrixApproxTemp, *matrixApproxOut;
691
692 matrixApproxTemp = malloc(matrixInRow*matrixInCol*sizeof(double));
693 matrixApproxOut = malloc(matrixInRow*matrixInCol*sizeof(double));
694
695 verbatim_copy(matrixIn+3*step*matrixInRow*matrixInCol,
696 matrixInRow*matrixInCol,
697 matrixApproxTemp, matrixInRow*matrixInCol);
698
699 for(count=0;count<step;count++)
700 {
701 iswt2(matrixApproxTemp, matrixIn+(step-count-1)*matrixInRow*matrixInCol,
702 matrixIn+(step-count-1)*matrixInRow*matrixInCol+step*matrixInRow*matrixInCol,
703 matrixIn+(step-count-1)*matrixInRow*matrixInCol+2*step*matrixInRow*matrixInCol,
704 matrixInRow, matrixInCol, matrixApproxOut, matrixOutRow, matrixOutCol,
705 filterLow, filterHi, filterLength, step-count);
706 if (count!=(step-1))
707 {
708 verbatim_copy(matrixApproxOut, matrixInRow*matrixInCol,
709 matrixApproxTemp, matrixInRow*matrixInCol);
710 }
711 }
712 free(matrixApproxTemp);
713
714 verbatim_copy(matrixApproxOut, matrixInRow*matrixInCol,
715 matrixOut, matrixInRow*matrixInCol);
716 free(matrixApproxOut);
717
718 return;
719 }
720