1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #undef powerof2 /* needed for Mac OS X */
8
9 /*---------------------------------------------------------------------------*/
10 /*
11 This file contains routines for performing wavelet analysis of
12 time series data.
13
14 File: Wavelets.c
15 Author: B. Douglas Ward
16 Date: 28 March 2000
17 */
18
19 /*---------------------------------------------------------------------------*/
20 /*
21 Include code for fast wavelet transforms.
22 */
23
24 #include "Haar.c"
25 #include "Daubechies.c"
26
27
28 /*---------------------------------------------------------------------------*/
29 /*
30 Print error message and stop.
31 */
32
WA_error(char * message)33 void WA_error (char * message)
34 {
35 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
36 exit(1);
37 }
38
39
40 /*---------------------------------------------------------------------------*/
41 /*
42 Print time series data to screen.
43 */
44
ts_print(int npts,float * data)45 void ts_print (int npts, float * data)
46 {
47 int i;
48
49
50 for (i = 0; i < npts; i++)
51 {
52 printf ("%12.4f ", data[i]);
53 if (8*((i+1)/8) == i+1) printf (" \n");
54 }
55 printf (" \n");
56
57 }
58
59
60 /*---------------------------------------------------------------------------*/
61 /*
62 Write time series data to specified file.
63 */
64
ts_fprint(char * filename,int npts,float * data)65 void ts_fprint (char * filename, int npts, float * data)
66 {
67 int i;
68 FILE * outfile = NULL;
69
70
71 outfile = fopen (filename, "w");
72
73
74 for (i = 0; i < npts; i++)
75 {
76 fprintf (outfile, "%f", data[i]);
77 fprintf (outfile, " \n");
78 }
79
80 fclose (outfile);
81 }
82
83
84 /*---------------------------------------------------------------------------*/
85 /*
86 Calculate integer power of 2.
87 */
88
powerof2(int n)89 int powerof2 (int n)
90 {
91 int i, j;
92
93 j = 1;
94
95 if (n > 0)
96 for (i = 0; i < n; i++)
97 j *= 2;
98 else if (n < 0)
99 j = 0;
100
101 return (j);
102 }
103
104
105 /*---------------------------------------------------------------------------*/
106 /*
107 Calculate integer log base 2.
108 */
109
my_log2(int n)110 int my_log2 (int n)
111 {
112 int i;
113
114 i = floor (log(n)/log(2.0) + 1.0e-10);
115
116 return (i);
117 }
118
119 /*---------------------------------------------------------------------------*/
120 /*
121 Apply filter to wavelet coefficients.
122 */
123
FWT_1d_filter(float * filter,int N,float * s)124 void FWT_1d_filter
125 (
126 float * filter, /* array of filter coefficients */
127 int N, /* log2(NPTS) */
128 float * s /* array of wavelet coefficients */
129 )
130
131 {
132 int NPTS; /* number of usable data points from input data */
133 int ipts; /* wavelet coefficient index */
134
135
136 NPTS = powerof2 (N);
137
138 for (ipts = 0; ipts < NPTS; ipts++)
139 s[ipts] *= filter[ipts];
140
141 }
142
143
144 /*---------------------------------------------------------------------------*/
145 /*
146 Set up array indicating which wavelet coefficients to set to zero.
147 */
148
FWT_1d_stop_filter(int num_stop_filters,int * stop_band,int * stop_mintr,int * stop_maxtr,int NFirst,int NPTS)149 float * FWT_1d_stop_filter
150 (
151 int num_stop_filters, /* number of wavelet stop filters */
152 int * stop_band, /* wavelet filter stop band */
153 int * stop_mintr, /* min. time for stop band */
154 int * stop_maxtr, /* max. time for stop band */
155 int NFirst, /* first image from input 3d+time dataset to use */
156 int NPTS /* number of usable data points from input data */
157 )
158
159 {
160 int N; /* log2(NPTS) */
161 int ipts; /* wavelet coefficient index */
162 int band; /* frequency band */
163 int mintr; /* min. value for time window (in TR) */
164 int maxtr; /* max. value for time window (in TR) */
165 int ifilter; /* filter index */
166 float * stop_filter = NULL; /* select wavelet coefficients to stop */
167
168
169 N = my_log2(NPTS);
170 stop_filter = (float *) malloc (sizeof(float) * NPTS); MTEST (stop_filter);
171
172
173 for (ipts = 0; ipts < NPTS; ipts++)
174 {
175 if (ipts == 0)
176 {
177 band = -1;
178 mintr = 0;
179 maxtr = NPTS-1;
180 }
181 else
182 {
183 band = my_log2(ipts);
184 mintr = (ipts - powerof2(band)) * powerof2(N-band);
185 maxtr = mintr + powerof2(N-band) - 1;
186 }
187
188 mintr += NFirst;
189 maxtr += NFirst;
190
191 stop_filter[ipts] = 1.0;
192 for (ifilter = 0; ifilter < num_stop_filters; ifilter++)
193 {
194 if (band == stop_band[ifilter])
195 {
196 if ((mintr >= stop_mintr[ifilter])
197 && (maxtr <= stop_maxtr[ifilter]))
198 stop_filter[ipts] = 0.0;
199 }
200 }
201
202 }
203
204
205 return (stop_filter);
206
207 }
208
209
210 /*---------------------------------------------------------------------------*/
211 /*
212 Set up array indicating which wavelet coefficients to include in the model.
213 */
214
FWT_1d_pass_filter(int num_pass_filters,int * pass_band,int * pass_mintr,int * pass_maxtr,int NFirst,int NPTS)215 float * FWT_1d_pass_filter
216 (
217 int num_pass_filters, /* number of wavelet pass filters */
218 int * pass_band, /* wavelet filter pass band */
219 int * pass_mintr, /* min. time for pass band */
220 int * pass_maxtr, /* max. time for pass band */
221 int NFirst, /* first image from input 3d+time dataset to use */
222 int NPTS /* number of usable data points from input data */
223 )
224
225 {
226 int N; /* log2(NPTS) */
227 int ipts; /* wavelet coefficient index */
228 int band; /* frequency band */
229 int mintr; /* min. value for time window (in TR) */
230 int maxtr; /* max. value for time window (in TR) */
231 int ifilter; /* filter index */
232 float * pass_filter = NULL; /* select wavelet coefficients to pass */
233
234
235 N = my_log2 (NPTS);
236 pass_filter = (float *) malloc (sizeof(float) * NPTS); MTEST (pass_filter);
237
238
239 for (ipts = 0; ipts < NPTS; ipts++)
240 {
241 if (ipts == 0)
242 {
243 band = -1;
244 mintr = 0;
245 maxtr = NPTS-1;
246 }
247 else
248 {
249 band = my_log2(ipts);
250 mintr = (ipts - powerof2(band)) * powerof2(N-band);
251 maxtr = mintr + powerof2(N-band) - 1;
252 }
253
254 mintr += NFirst;
255 maxtr += NFirst;
256
257 pass_filter[ipts] = 0.0;
258 for (ifilter = 0; ifilter < num_pass_filters; ifilter++)
259 {
260 if (band == pass_band[ifilter])
261 {
262 if ((mintr >= pass_mintr[ifilter])
263 && (maxtr <= pass_maxtr[ifilter]))
264 pass_filter[ipts] = 1.0;
265 }
266 }
267
268 }
269
270
271 return (pass_filter);
272
273 }
274
275
276 /*---------------------------------------------------------------------------*/
277 /*
278 Calculate the error sum of squares.
279 */
280
calc_sse(int NPTS,float * trueData,float * est)281 float calc_sse
282 (
283 int NPTS, /* number of usable data points from input data */
284 float * trueData, /* actual time series data */
285 float * est /* estimated time series data */
286 )
287
288 {
289 int ipt; /* time point index */
290 float diff; /* difference between actual and estimated data */
291 float sse; /* estimation error sum of squares */
292
293 sse = 0.0;
294 for (ipt = 0; ipt < NPTS; ipt++)
295 {
296 diff = trueData[ipt] - est[ipt];
297 sse += diff * diff;
298 }
299
300 return (sse);
301 }
302
303
304 /*---------------------------------------------------------------------------*/
305 /*
306 Calculate the F-statistic for significance of the regression.
307 */
308
calc_freg(int n,int p,int q,float ssef,float sser)309 float calc_freg
310 (
311 int n, /* number of data points */
312 int p, /* number of parameters in the full model */
313 int q, /* number of parameters in the rdcd model */
314 float ssef, /* error sum of squares from full model */
315 float sser /* error sum of squares from reduced model */
316 )
317
318 {
319 const float MAXF = 1000.0; /* maximum value for F-statistic */
320 const float EPSILON = 1.0e-2; /* protection against divide by zero */
321 float msef; /* mean square error for the full model */
322 float msreg; /* mean square due to the regression */
323 float freg; /* F-statistic for the full regression model */
324
325
326 /*----- Order of reduced model = order of full model ??? -----*/
327 if (p <= q) return (0.0);
328
329
330 /*----- Calculate numerator and denominator of F-statistic -----*/
331 msreg = (sser - ssef) / (p - q); if (msreg < 0.0) msreg = 0.0;
332 msef = ssef / (n - p); if (msef < 0.0) msef = 0.0;
333
334
335 if (msef < EPSILON)
336 freg = 0.0;
337 else
338 if (msreg > MAXF*msef) freg = MAXF;
339 else freg = msreg / msef;
340
341
342 /*----- Limit range of values for F-statistic -----*/
343 if (freg < 0.0) freg = 0.0;
344 if (freg > MAXF) freg = MAXF;
345
346
347 /*----- Return F-statistic for significance of the regression -----*/
348 return (freg);
349
350 }
351
352
353 /*---------------------------------------------------------------------------*/
354 /*
355 Calculate the coefficient of multiple determination R^2.
356 */
357
calc_rsqr(float ssef,float sser)358 float calc_rsqr
359 (
360 float ssef, /* error sum of squares from full model */
361 float sser /* error sum of squares from reduced model */
362 )
363
364 {
365 const float EPSILON = 1.0e-2; /* protection against divide by zero */
366 float rsqr; /* coeff. of multiple determination R^2 */
367
368
369 /*----- coefficient of multiple determination R^2 -----*/
370 if (sser < EPSILON)
371 rsqr = 0.0;
372 else
373 rsqr = (sser - ssef) / sser;
374
375
376 /*----- Limit range of values for R^2 -----*/
377 if (rsqr < 0.0) rsqr = 0.0;
378 if (rsqr > 1.0) rsqr = 1.0;
379
380
381 /*----- Return coefficient of multiple determination R^2 -----*/
382 return (rsqr);
383 }
384
385
386 /*---------------------------------------------------------------------------*/
387 /*
388 Perform wavelet analysis on a single input time series.
389 */
390
wavelet_analysis(int wavelet_type,int f,float * stop_filter,int q,float * base_filter,int p,float * full_filter,int NPTS,float * ts_array,float * coef,float * sse_base,float * sse_full,float * ffull,float * rfull,float * coefts,float * fitts,float * sgnlts,float * errts)391 void wavelet_analysis
392 (
393 int wavelet_type, /* indicates type of wavelet */
394 int f, /* number of parameters removed by the filter */
395 float * stop_filter, /* select wavelet coefs. to stop */
396 int q, /* number of parameters in the baseline model */
397 float * base_filter, /* select wavelet coefs. for baseline */
398 int p, /* number of parameters in the full model */
399 float * full_filter, /* select wavelet coefs. for full model */
400 int NPTS, /* number of usable data points from input data */
401 float * ts_array, /* array of measured data for one voxel */
402
403 float * coef, /* full model wavelet coefficients */
404 float * sse_base, /* baseline model error sum of squares */
405 float * sse_full, /* full model error sum of squares */
406 float * ffull, /* full model F-statistic */
407 float * rfull, /* full model R^2 stats. */
408
409 float * coefts, /* filtered FWT coefficients */
410 float * fitts, /* filterd time series */
411 float * sgnlts, /* signal model fitted time series */
412 float * errts /* residual error time series */
413 )
414
415 {
416 int N; /* log2(NPTS) */
417 int it; /* time point index */
418 int ip; /* full model parameter index */
419 float * filtts = NULL; /* stop filtered time series */
420 float * basefwt = NULL; /* baseline model FWT coefficients */
421 float * basets = NULL; /* baseline model time series */
422 float * fullfwt = NULL; /* full model FWT coefficients */
423 float * fullts = NULL; /* full model time series */
424
425
426 /*----- Initialize local variables -----*/
427 N = my_log2(NPTS);
428
429
430 /*----- Allocate memory for arrays -----*/
431 filtts = (float *) malloc (sizeof(float) * NPTS); MTEST (filtts);
432 basefwt = (float *) malloc (sizeof(float) * NPTS); MTEST (basefwt);
433 basets = (float *) malloc (sizeof(float) * NPTS); MTEST (basets);
434 fullfwt = (float *) malloc (sizeof(float) * NPTS); MTEST (fullfwt);
435 fullts = (float *) malloc (sizeof(float) * NPTS); MTEST (fullts);
436
437
438 /*----- Forward Fast Wavelet Transform -----*/
439 for (it = 0; it < NPTS; it++)
440 coefts[it] = ts_array[it];
441 switch (wavelet_type)
442 {
443 case WA_HAAR:
444 Haar_forward_FWT_1d (N, coefts);
445 break;
446
447 case WA_DAUBECHIES:
448 Daubechies_forward_FWT_1d (N, coefts);
449 break;
450 }
451
452
453 /*----- Apply stop filter to wavelet transform coefficients -----*/
454 FWT_1d_filter (stop_filter, N, coefts);
455
456
457 /*----- Inverse Fast Wavelet Transform of filtered FWT -----*/
458 for (it = 0; it < NPTS; it++)
459 filtts[it] = coefts[it];
460 switch (wavelet_type)
461 {
462 case WA_HAAR:
463 Haar_inverse_FWT_1d (N, filtts);
464 break;
465
466 case WA_DAUBECHIES:
467 Daubechies_inverse_FWT_1d (N, filtts);
468 break;
469 }
470
471
472 /*----- Apply baseline pass filter to wavelet transform coefficients -----*/
473 for (it = 0; it < NPTS; it++)
474 basefwt[it] = coefts[it];
475 FWT_1d_filter (base_filter, N, basefwt);
476
477
478 /*----- Inverse Fast Wavelet Transform of baseline FWT -----*/
479 for (it = 0; it < NPTS; it++)
480 basets[it] = basefwt[it];
481 switch (wavelet_type)
482 {
483 case WA_HAAR:
484 Haar_inverse_FWT_1d (N, basets);
485 break;
486
487 case WA_DAUBECHIES:
488 Daubechies_inverse_FWT_1d (N, basets);
489 break;
490 }
491
492
493 /*----- Calculate error sum of squares (SSE) for baseline model -----*/
494 *sse_base = calc_sse (NPTS, filtts, basets);
495
496
497 /*----- Apply full model filter to wavelet transform coefficients -----*/
498 for (it = 0; it < NPTS; it++)
499 fullfwt[it] = coefts[it];
500 FWT_1d_filter (full_filter, N, fullfwt);
501
502
503 /*----- Save full model wavelet coefficients -----*/
504 ip = 0;
505 for (it = 0; it < NPTS; it++)
506 if (full_filter[it] == 1.0)
507 {
508 coef[ip] = fullfwt[it];
509 ip++;
510 if (ip >= p) break;
511 }
512
513
514 /*----- Inverse Fast Wavelet Transform of full model FWT -----*/
515 for (it = 0; it < NPTS; it++)
516 fullts[it] = fullfwt[it];
517 switch (wavelet_type)
518 {
519 case WA_HAAR:
520 Haar_inverse_FWT_1d (N, fullts);
521 break;
522
523 case WA_DAUBECHIES:
524 Daubechies_inverse_FWT_1d (N, fullts);
525 break;
526 }
527
528
529 /*----- Calculate error sum of squares (SSE) for signal model -----*/
530 *sse_full = calc_sse (NPTS, filtts, fullts);
531
532
533 /*----- Calculate coefficient of multiple determination R^2 -----*/
534 *rfull = calc_rsqr (*sse_full, *sse_base);
535
536
537 /*----- Calculate F-statistic for significance of the signal model -----*/
538 *ffull = calc_freg (NPTS-f, p, q, *sse_full, *sse_base);
539
540
541 /*----- Calculate residual errors -----*/
542 for (it = 0; it < NPTS; it++)
543 {
544 if (p == 0)
545 errts[it] = ts_array[it] - filtts[it];
546 else
547 errts[it] = filtts[it] - fullts[it];
548 }
549
550
551 /*----- Calculate "pure" signal time series -----*/
552 for (it = 0; it < NPTS; it++)
553 sgnlts[it] = fullts[it] - basets[it];
554
555
556 /*----- Save the fitted time series -----*/
557 for (it = 0; it < NPTS; it++)
558 {
559 if (p == 0)
560 fitts[it] = filtts[it];
561 else
562 fitts[it] = fullts[it];
563 }
564
565
566 /*----- Release memory -----*/
567 free (filtts); filtts = NULL;
568 free (basefwt); basefwt = NULL;
569 free (basets); basets = NULL;
570 free (fullfwt); fullfwt = NULL;
571 free (fullts); fullts = NULL;
572
573
574 return;
575
576 }
577
578
579 /*---------------------------------------------------------------------------*/
580 /*
581 Convert F-value to p-value.
582 This routine was copied from: mri_stats.c - linked with libmri.a?
583 */
584
585 #if 0
586 double fstat_t2p( double ff , double dofnum , double dofden )
587 {
588 int which , status ;
589 double p , q , f , dfn , dfd , bound ;
590
591 which = 1 ;
592 p = 0.0 ;
593 q = 0.0 ;
594 f = ff ;
595 dfn = dofnum ;
596 dfd = dofden ;
597
598 cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
599
600 if( status == 0 ) return q ;
601 else return 1.0 ;
602 }
603 #endif
604
605 /*---------------------------------------------------------------------------*/
606 /*
607 Report the results of wavelet analysis for a single time series.
608 */
609
610 static char lbuf[4096]; /* character string containing statistical summary */
611 static char sbuf[256];
612
report_results(int N,int NFirst,int f,int p,int q,float * base_filter,float * sgnl_filter,float * coef,float sse_base,float sse_full,float ffull,float rfull,char ** label)613 void report_results
614 (
615 int N, /* number of usable data points from input data */
616 int NFirst, /* first image from input 3d+time dataset to use */
617 int f, /* number of parameters removed by the filter */
618 int p, /* number of parameters in the full model */
619 int q, /* number of parameters in the baseline model */
620 float * base_filter, /* select wavelet coefs. for baseline */
621 float * sgnl_filter, /* select wavelet coefs. for full model */
622 float * coef, /* full model wavelet coefficients */
623 float sse_base, /* baseline model error sum of squares */
624 float sse_full, /* full model error sum of squares */
625 float ffull, /* full model F-statistic */
626 float rfull, /* full model R^2 stat. */
627 char ** label /* statistical summary for ouput display */
628 )
629
630 {
631 int it; /* time index */
632 int icoef; /* coefficient index */
633 double pvalue; /* p-value */
634
635
636 /** 22 Apr 1997: create label if desired by AFNI **/
637 /** [This is in static storage, since AFNI will copy it] **/
638
639 if( label != NULL ){ /* assemble this 1 line at a time from sbuf */
640
641 lbuf[0] = '\0' ; /* make this a 0 length string to start */
642
643 /** for each reference, make a string into sbuf **/
644
645
646 /*----- Display wavelet coefficients for full model -----*/
647 icoef = 0;
648 for (it = 0; it < N; it++)
649 {
650 if (sgnl_filter[it])
651 {
652 {
653 int band, mintr, maxtr;
654
655 if (it == 0)
656 {
657 band = -1;
658 mintr = 0;
659 maxtr = N-1;
660 }
661 else
662 {
663 band = my_log2(it);
664 mintr = (it - powerof2(band)) * powerof2(my_log2(N)-band);
665 maxtr = mintr + powerof2(my_log2(N)-band) - 1;
666 }
667
668 mintr += NFirst;
669 maxtr += NFirst;
670
671 if (base_filter[it])
672 sprintf (sbuf, "B(%2d)[%3d,%3d] = %f \n",
673 band, mintr, maxtr, coef[icoef]);
674 else
675 sprintf (sbuf, "S(%2d)[%3d,%3d] = %f \n",
676 band, mintr, maxtr, coef[icoef]);
677 }
678
679 strcat(lbuf,sbuf);
680
681 icoef++;
682 }
683
684 } /* End loop over Wavelet coefficients */
685
686
687 /*----- Statistical results for baseline fit -----*/
688 sprintf (sbuf, "\nBaseline: \n");
689 strcat(lbuf,sbuf);
690 sprintf (sbuf, "# params = %4d \n", q);
691 strcat (lbuf, sbuf);
692 sprintf (sbuf, "SSE = %10.3f \n", sse_base);
693 strcat (lbuf, sbuf);
694 sprintf (sbuf, "MSE = %10.3f \n", sse_base/(N-f-q));
695 strcat (lbuf, sbuf);
696
697
698 /*----- Statistical results for full model -----*/
699 sprintf (sbuf, "\nFull Model: \n");
700 strcat (lbuf, sbuf);
701 sprintf (sbuf, "# params = %4d \n", p);
702 strcat (lbuf, sbuf);
703 sprintf (sbuf, "SSE = %10.3f \n", sse_full);
704 strcat (lbuf, sbuf);
705 sprintf (sbuf, "MSE = %10.3f \n", sse_full/(N-f-p));
706 strcat (lbuf, sbuf);
707
708
709 /*----- Summary statistics -----*/
710 sprintf (sbuf, "\nSummary Statistics: \n");
711 strcat (lbuf, sbuf);
712
713 sprintf (sbuf, "R^2 = %10.3f \n", rfull);
714 strcat (lbuf, sbuf);
715
716 sprintf (sbuf, "F[%2d,%3d] = %10.3f \n", p-q, N-f-p, ffull);
717 strcat (lbuf, sbuf);
718
719 pvalue = fstat_t2p ( (double) ffull, (double) p-q, (double) N-f-p);
720 sprintf (sbuf, "p-value = %e \n", pvalue);
721 strcat (lbuf, sbuf);
722
723
724 *label = lbuf ; /* send address of lbuf back in what label points to */
725 }
726
727 }
728
729
730 /*---------------------------------------------------------------------------*/
731
732
733
734
735
736
737
738