00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #undef powerof2  /* needed for Mac OS X */
00008 
00009 /*---------------------------------------------------------------------------*/
00010 /*
00011   This file contains routines for performing wavelet analysis of
00012   time series data.
00013 
00014   File:    Wavelets.c
00015   Author:  B. Douglas Ward
00016   Date:    28 March 2000
00017 */
00018 
00019 /*---------------------------------------------------------------------------*/
00020 /*
00021   Include code for fast wavelet transforms.
00022 */
00023 
00024 #include "Haar.c"
00025 #include "Daubechies.c"
00026 
00027 
00028 /*---------------------------------------------------------------------------*/
00029 /*
00030    Print error message and stop.
00031 */
00032 
00033 void WA_error (char * message)
00034 {
00035   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00036   exit(1);
00037 }
00038 
00039 
00040 /*---------------------------------------------------------------------------*/
00041 /*
00042   Print time series data to screen.
00043 */
00044 
00045 void ts_print (int npts, float * data)
00046 {
00047   int i;
00048 
00049 
00050   for (i = 0;  i < npts;  i++)
00051     {
00052       printf ("%12.4f  ", data[i]);
00053       if (8*((i+1)/8) == i+1)  printf (" \n");
00054     }
00055       printf (" \n");
00056 
00057 }
00058 
00059 
00060 /*---------------------------------------------------------------------------*/
00061 /*
00062   Write time series data to specified file.
00063 */
00064 
00065 void ts_fprint (char * filename, int npts, float * data)
00066 {
00067   int i;
00068   FILE * outfile = NULL;
00069 
00070 
00071   outfile = fopen (filename, "w");
00072 
00073 
00074   for (i = 0;  i < npts;  i++)
00075     {
00076       fprintf (outfile, "%f", data[i]);
00077       fprintf (outfile, " \n");
00078     }
00079 
00080   fclose (outfile);
00081 }
00082 
00083 
00084 /*---------------------------------------------------------------------------*/
00085 /*
00086   Calculate integer power of 2.
00087 */
00088 
00089 int powerof2 (int n)
00090 {
00091   int i, j;
00092 
00093   j = 1;
00094 
00095   if (n > 0)
00096     for (i = 0;  i < n;  i++)
00097       j *= 2;
00098   else if (n < 0)
00099     j = 0;
00100 
00101   return (j);
00102 }
00103 
00104 
00105 /*---------------------------------------------------------------------------*/
00106 /*
00107   Calculate integer log base 2.
00108 */
00109 
00110 int my_log2 (int n)
00111 {
00112   int i;
00113    
00114   i = floor (log(n)/log(2.0) + 1.0e-10);
00115 
00116   return (i);
00117 }
00118 
00119 /*---------------------------------------------------------------------------*/
00120 /*
00121   Apply filter to wavelet coefficients.
00122 */
00123 
00124 void FWT_1d_filter 
00125 (
00126   float * filter,         /* array of filter coefficients */
00127   int N,                  /* log2(NPTS) */
00128   float * s               /* array of wavelet coefficients */ 
00129 )
00130 
00131 {
00132   int NPTS;        /* number of usable data points from input data */
00133   int ipts;        /* wavelet coefficient index */
00134 
00135 
00136   NPTS = powerof2 (N);
00137 
00138   for (ipts = 0;  ipts < NPTS;  ipts++)
00139     s[ipts] *= filter[ipts];
00140 
00141 }
00142 
00143 
00144 /*---------------------------------------------------------------------------*/
00145 /*
00146   Set up array indicating which wavelet coefficients to set to zero.
00147 */
00148 
00149 float * FWT_1d_stop_filter 
00150 (
00151   int num_stop_filters,   /* number of wavelet stop filters */
00152   int * stop_band,        /* wavelet filter stop band */ 
00153   int * stop_mintr,       /* min. time for stop band */
00154   int * stop_maxtr,       /* max. time for stop band */
00155   int NFirst,             /* first image from input 3d+time dataset to use */
00156   int NPTS                /* number of usable data points from input data */
00157 )
00158 
00159 {
00160   int N;                       /* log2(NPTS) */
00161   int ipts;                    /* wavelet coefficient index */
00162   int band;                    /* frequency band */
00163   int mintr;                   /* min. value for time window (in TR) */
00164   int maxtr;                   /* max. value for time window (in TR) */
00165   int ifilter;                 /* filter index */
00166   float * stop_filter = NULL;  /* select wavelet coefficients to stop */
00167 
00168 
00169   N = my_log2(NPTS);
00170   stop_filter = (float *) malloc (sizeof(float) * NPTS);   MTEST (stop_filter);
00171   
00172 
00173   for (ipts = 0;  ipts < NPTS;  ipts++)
00174     {
00175       if (ipts == 0)
00176         {
00177           band = -1;
00178           mintr = 0;
00179           maxtr = NPTS-1;
00180         }
00181       else
00182         {
00183           band = my_log2(ipts);
00184           mintr = (ipts - powerof2(band)) * powerof2(N-band);
00185           maxtr = mintr + powerof2(N-band) - 1;
00186         }
00187 
00188       mintr += NFirst;
00189       maxtr += NFirst;
00190 
00191       stop_filter[ipts] = 1.0;
00192       for (ifilter = 0;  ifilter < num_stop_filters;  ifilter++)
00193         {
00194           if (band == stop_band[ifilter])
00195             {
00196               if ((mintr >= stop_mintr[ifilter]) 
00197                   && (maxtr <= stop_maxtr[ifilter]))
00198                 stop_filter[ipts] = 0.0;
00199             }
00200         }
00201 
00202     }
00203 
00204   
00205   return (stop_filter);
00206 
00207 }
00208 
00209 
00210 /*---------------------------------------------------------------------------*/
00211 /*
00212   Set up array indicating which wavelet coefficients to include in the model.
00213 */
00214 
00215 float * FWT_1d_pass_filter
00216 (
00217   int num_pass_filters,   /* number of wavelet pass filters */
00218   int * pass_band,        /* wavelet filter pass band */ 
00219   int * pass_mintr,       /* min. time for pass band */
00220   int * pass_maxtr,       /* max. time for pass band */
00221   int NFirst,             /* first image from input 3d+time dataset to use */
00222   int NPTS                /* number of usable data points from input data */
00223 )
00224 
00225 {
00226   int N;                       /* log2(NPTS) */
00227   int ipts;                    /* wavelet coefficient index */
00228   int band;                    /* frequency band */
00229   int mintr;                   /* min. value for time window (in TR) */
00230   int maxtr;                   /* max. value for time window (in TR) */
00231   int ifilter;                 /* filter index */
00232   float * pass_filter = NULL;  /* select wavelet coefficients to pass */
00233 
00234 
00235   N = my_log2 (NPTS);
00236   pass_filter = (float *) malloc (sizeof(float) * NPTS);   MTEST (pass_filter);
00237   
00238 
00239   for (ipts = 0;  ipts < NPTS;  ipts++)
00240     {
00241       if (ipts == 0)
00242         {
00243           band = -1;
00244           mintr = 0;
00245           maxtr = NPTS-1;
00246         }
00247       else
00248         {
00249           band = my_log2(ipts);
00250           mintr = (ipts - powerof2(band)) * powerof2(N-band);
00251           maxtr = mintr + powerof2(N-band) - 1;
00252         }
00253 
00254       mintr += NFirst;
00255       maxtr += NFirst;
00256 
00257       pass_filter[ipts] = 0.0;
00258       for (ifilter = 0;  ifilter < num_pass_filters;  ifilter++)
00259         {
00260           if (band == pass_band[ifilter])
00261             {
00262               if ((mintr >= pass_mintr[ifilter]) 
00263                   && (maxtr <= pass_maxtr[ifilter]))
00264                 pass_filter[ipts] = 1.0;
00265             }
00266         }
00267 
00268     }
00269 
00270   
00271   return (pass_filter);
00272 
00273 }
00274 
00275 
00276 /*---------------------------------------------------------------------------*/
00277 /*
00278   Calculate the error sum of squares.
00279 */
00280 
00281 float calc_sse 
00282 (
00283   int NPTS,         /* number of usable data points from input data */
00284   float * trueData,     /* actual time series data */
00285   float * est       /* estimated time series data */
00286 )
00287 
00288 {
00289   int ipt;          /* time point index */
00290   float diff;       /* difference between actual and estimated data */
00291   float sse;        /* estimation error sum of squares */
00292 
00293   sse = 0.0;
00294   for (ipt = 0;  ipt < NPTS;  ipt++)
00295     {
00296       diff = trueData[ipt] - est[ipt];
00297       sse += diff * diff;
00298     }
00299   
00300   return (sse);
00301 }
00302 
00303 
00304 /*---------------------------------------------------------------------------*/
00305 /*
00306   Calculate the F-statistic for significance of the regression.
00307 */
00308 
00309 float calc_freg
00310 (
00311   int n,                      /* number of data points */
00312   int p,                      /* number of parameters in the full model */
00313   int q,                      /* number of parameters in the rdcd model */
00314   float ssef,                 /* error sum of squares from full model */
00315   float sser                  /* error sum of squares from reduced model */
00316 )
00317 
00318 {
00319   const float MAXF = 1000.0;         /* maximum value for F-statistic */
00320   const float EPSILON = 1.0e-2;      /* protection against divide by zero */
00321   float msef;                 /* mean square error for the full model */
00322   float msreg;                /* mean square due to the regression */
00323   float freg;                 /* F-statistic for the full regression model */
00324 
00325 
00326   /*----- Order of reduced model = order of full model ??? -----*/
00327   if (p <= q)   return (0.0);
00328 
00329 
00330   /*----- Calculate numerator and denominator of F-statistic -----*/
00331   msreg = (sser - ssef) / (p - q);    if (msreg < 0.0)  msreg = 0.0;
00332   msef   = ssef / (n - p);            if (msef  < 0.0)  msef  = 0.0;
00333 
00334 
00335   if (msef < EPSILON)
00336     freg = 0.0;
00337   else
00338     if (msreg > MAXF*msef)  freg = MAXF;
00339     else                    freg = msreg / msef;
00340 
00341 
00342   /*----- Limit range of values for F-statistic -----*/
00343   if (freg < 0.0)   freg = 0.0;
00344   if (freg > MAXF)  freg = MAXF;
00345 
00346 
00347   /*----- Return F-statistic for significance of the regression -----*/
00348   return (freg);
00349 
00350 }
00351 
00352 
00353 /*---------------------------------------------------------------------------*/
00354 /*
00355   Calculate the coefficient of multiple determination R^2.
00356 */
00357 
00358 float calc_rsqr 
00359 (
00360   float ssef,                 /* error sum of squares from full model */
00361   float sser                  /* error sum of squares from reduced model */
00362 )
00363 
00364 {
00365   const float EPSILON = 1.0e-2;     /* protection against divide by zero */
00366   float rsqr;                       /* coeff. of multiple determination R^2  */
00367 
00368 
00369   /*----- coefficient of multiple determination R^2 -----*/
00370   if (sser < EPSILON)
00371     rsqr = 0.0;
00372   else
00373     rsqr = (sser - ssef) / sser;
00374 
00375 
00376   /*----- Limit range of values for R^2 -----*/
00377   if (rsqr < 0.0)   rsqr = 0.0;
00378   if (rsqr > 1.0)   rsqr = 1.0;
00379 
00380 
00381   /*----- Return coefficient of multiple determination R^2 -----*/
00382   return (rsqr);
00383 }
00384 
00385 
00386 /*---------------------------------------------------------------------------*/
00387 /*
00388   Perform wavelet analysis on a single input time series.
00389 */
00390 
00391 void wavelet_analysis 
00392 (
00393   int wavelet_type,         /* indicates type of wavelet */   
00394   int f,                    /* number of parameters removed by the filter */
00395   float * stop_filter,      /* select wavelet coefs. to stop */
00396   int q,                    /* number of parameters in the baseline model */
00397   float * base_filter,      /* select wavelet coefs. for baseline */
00398   int p,                    /* number of parameters in the full model */
00399   float * full_filter,      /* select wavelet coefs. for full model */
00400   int NPTS,                 /* number of usable data points from input data */
00401   float * ts_array,         /* array of measured data for one voxel */
00402 
00403   float * coef,             /* full model wavelet coefficients */
00404   float * sse_base,         /* baseline model error sum of squares */
00405   float * sse_full,         /* full model error sum of squares */
00406   float * ffull,            /* full model F-statistic */
00407   float * rfull,            /* full model R^2 stats. */
00408 
00409   float * coefts,           /* filtered FWT coefficients */
00410   float * fitts,            /* filterd time series */
00411   float * sgnlts,           /* signal model fitted time series */
00412   float * errts             /* residual error time series */
00413 )
00414 
00415 {
00416   int N;                    /* log2(NPTS) */
00417   int it;                   /* time point index */
00418   int ip;                   /* full model parameter index */
00419   float * filtts = NULL;    /* stop filtered time series */
00420   float * basefwt = NULL;   /* baseline model FWT coefficients */
00421   float * basets = NULL;    /* baseline model time series */
00422   float * fullfwt = NULL;   /* full model FWT coefficients */
00423   float * fullts = NULL;    /* full model time series */
00424 
00425     
00426   /*----- Initialize local variables -----*/
00427   N = my_log2(NPTS);
00428 
00429 
00430   /*----- Allocate memory for arrays -----*/
00431   filtts     = (float *) malloc (sizeof(float) * NPTS);    MTEST (filtts);
00432   basefwt    = (float *) malloc (sizeof(float) * NPTS);    MTEST (basefwt);
00433   basets     = (float *) malloc (sizeof(float) * NPTS);    MTEST (basets);
00434   fullfwt    = (float *) malloc (sizeof(float) * NPTS);    MTEST (fullfwt);
00435   fullts     = (float *) malloc (sizeof(float) * NPTS);    MTEST (fullts);
00436 
00437 
00438   /*----- Forward Fast Wavelet Transform -----*/
00439   for (it = 0;  it < NPTS;  it++)
00440     coefts[it] = ts_array[it];
00441   switch (wavelet_type)
00442     {
00443     case WA_HAAR:
00444       Haar_forward_FWT_1d (N, coefts);
00445       break;
00446 
00447     case WA_DAUBECHIES:
00448       Daubechies_forward_FWT_1d (N, coefts);
00449       break;
00450     }
00451 
00452 
00453   /*----- Apply stop filter to wavelet transform coefficients -----*/
00454   FWT_1d_filter (stop_filter, N, coefts);
00455 
00456 
00457   /*----- Inverse Fast Wavelet Transform of filtered FWT -----*/
00458   for (it = 0;  it < NPTS;  it++)
00459     filtts[it] = coefts[it];
00460   switch (wavelet_type)
00461     {
00462     case WA_HAAR:
00463       Haar_inverse_FWT_1d (N, filtts);
00464       break;
00465 
00466     case WA_DAUBECHIES:
00467       Daubechies_inverse_FWT_1d (N, filtts);
00468       break;
00469     }
00470 
00471 
00472   /*----- Apply baseline pass filter to wavelet transform coefficients -----*/
00473   for (it = 0;  it < NPTS;  it++)
00474     basefwt[it] = coefts[it];
00475   FWT_1d_filter (base_filter, N, basefwt);
00476 
00477 
00478  /*----- Inverse Fast Wavelet Transform of baseline FWT -----*/
00479   for (it = 0;  it < NPTS;  it++)
00480     basets[it] = basefwt[it];
00481   switch (wavelet_type)
00482     {
00483     case WA_HAAR:
00484       Haar_inverse_FWT_1d (N, basets);
00485       break;
00486 
00487     case WA_DAUBECHIES:
00488       Daubechies_inverse_FWT_1d (N, basets);
00489       break;
00490     }
00491 
00492 
00493   /*----- Calculate error sum of squares (SSE) for baseline model -----*/
00494   *sse_base = calc_sse (NPTS, filtts, basets);
00495 
00496 
00497   /*----- Apply full model filter to wavelet transform coefficients -----*/
00498   for (it = 0;  it < NPTS;  it++)
00499     fullfwt[it] = coefts[it];
00500   FWT_1d_filter (full_filter, N, fullfwt);
00501 
00502 
00503   /*----- Save full model wavelet coefficients -----*/
00504   ip = 0;
00505   for (it = 0;  it < NPTS;  it++)
00506     if (full_filter[it] == 1.0)
00507       {
00508         coef[ip] = fullfwt[it];
00509         ip++;
00510         if (ip >= p)  break;
00511       }
00512       
00513 
00514  /*----- Inverse Fast Wavelet Transform of full model FWT -----*/
00515   for (it = 0;  it < NPTS;  it++)
00516     fullts[it] = fullfwt[it];
00517   switch (wavelet_type)
00518     {
00519     case WA_HAAR:
00520       Haar_inverse_FWT_1d (N, fullts);
00521       break;
00522 
00523     case WA_DAUBECHIES:
00524       Daubechies_inverse_FWT_1d (N, fullts);
00525       break;
00526     }
00527 
00528 
00529   /*----- Calculate error sum of squares (SSE) for signal model -----*/
00530   *sse_full = calc_sse (NPTS, filtts, fullts);
00531 
00532 
00533   /*----- Calculate coefficient of multiple determination R^2 -----*/
00534   *rfull = calc_rsqr (*sse_full, *sse_base);
00535 
00536 
00537   /*----- Calculate F-statistic for significance of the signal model -----*/
00538   *ffull = calc_freg (NPTS-f, p, q, *sse_full, *sse_base);
00539 
00540 
00541   /*----- Calculate residual errors -----*/
00542   for (it = 0;  it < NPTS;  it++)
00543     {
00544       if (p == 0)
00545         errts[it] = ts_array[it] - filtts[it];
00546       else
00547         errts[it] = filtts[it] - fullts[it];
00548     }
00549 
00550 
00551   /*----- Calculate "pure" signal time series -----*/
00552   for (it = 0;  it < NPTS;  it++)
00553     sgnlts[it] = fullts[it] - basets[it];
00554 
00555 
00556   /*----- Save the fitted time series -----*/
00557   for (it = 0;  it < NPTS;  it++)
00558     {
00559       if (p == 0)
00560         fitts[it] = filtts[it];
00561       else
00562         fitts[it] = fullts[it];
00563     }
00564 
00565 
00566   /*----- Release memory -----*/
00567   free (filtts);      filtts = NULL;
00568   free (basefwt);     basefwt = NULL;
00569   free (basets);      basets = NULL;
00570   free (fullfwt);     fullfwt = NULL;
00571   free (fullts);      fullts = NULL;
00572 
00573 
00574   return;
00575 
00576 }
00577 
00578 
00579 /*---------------------------------------------------------------------------*/
00580 /*
00581   Convert F-value to p-value.  
00582   This routine was copied from: mri_stats.c
00583 */
00584 
00585 
00586 double fstat_t2p( double ff , double dofnum , double dofden )
00587 {
00588    int which , status ;
00589    double p , q , f , dfn , dfd , bound ;
00590 
00591    which  = 1 ;
00592    p      = 0.0 ;
00593    q      = 0.0 ;
00594    f      = ff ;
00595    dfn    = dofnum ;
00596    dfd    = dofden ;
00597 
00598    cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
00599 
00600    if( status == 0 ) return q ;
00601    else              return 1.0 ;
00602 }
00603 
00604 /*---------------------------------------------------------------------------*/
00605 /*
00606   Report the results of wavelet analysis for a single time series.
00607 */
00608 
00609 static char lbuf[4096];   /* character string containing statistical summary */
00610 static char sbuf[256];
00611 
00612 void report_results 
00613 (
00614   int N,                /* number of usable data points from input data */
00615   int NFirst,           /* first image from input 3d+time dataset to use */
00616   int f,                /* number of parameters removed by the filter */
00617   int p,                /* number of parameters in the full model */
00618   int q,                /* number of parameters in the baseline model */
00619   float * base_filter,  /* select wavelet coefs. for baseline */
00620   float * sgnl_filter,  /* select wavelet coefs. for full model */
00621   float * coef,         /* full model wavelet coefficients */
00622   float sse_base,       /* baseline model error sum of squares */
00623   float sse_full,       /* full model error sum of squares */
00624   float ffull,          /* full model F-statistic */
00625   float rfull,          /* full model R^2 stat. */
00626   char ** label         /* statistical summary for ouput display */
00627 )
00628 
00629 {
00630   int it;               /* time index */
00631   int icoef;            /* coefficient index */
00632   double pvalue;        /* p-value */
00633   
00634 
00635   /** 22 Apr 1997: create label if desired by AFNI         **/
00636   /** [This is in static storage, since AFNI will copy it] **/
00637   
00638   if( label != NULL ){  /* assemble this 1 line at a time from sbuf */
00639     
00640     lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00641     
00642     /** for each reference, make a string into sbuf **/
00643     
00644 
00645   /*----- Display wavelet coefficients for full model -----*/
00646   icoef = 0;
00647   for (it = 0;  it < N;  it++)        
00648     {                      
00649       if (sgnl_filter[it])
00650         {                             
00651           {
00652             int band, mintr, maxtr;
00653 
00654             if (it == 0)
00655               {
00656                 band = -1;
00657                 mintr = 0;
00658                 maxtr = N-1;
00659               }
00660             else
00661               {
00662                 band = my_log2(it);
00663                 mintr = (it - powerof2(band)) * powerof2(my_log2(N)-band);
00664                 maxtr = mintr + powerof2(my_log2(N)-band) - 1;
00665               }
00666             
00667             mintr += NFirst;
00668             maxtr += NFirst;
00669             
00670             if (base_filter[it])
00671               sprintf (sbuf, "B(%2d)[%3d,%3d] = %f \n", 
00672                        band, mintr, maxtr, coef[icoef]);
00673             else
00674               sprintf (sbuf, "S(%2d)[%3d,%3d] = %f \n", 
00675                        band, mintr, maxtr, coef[icoef]);
00676           }
00677 
00678           strcat(lbuf,sbuf);
00679 
00680           icoef++;
00681         }
00682             
00683     }  /* End loop over Wavelet coefficients */
00684 
00685 
00686     /*----- Statistical results for baseline fit -----*/ 
00687     sprintf (sbuf, "\nBaseline: \n");
00688     strcat(lbuf,sbuf);
00689     sprintf (sbuf, "# params  = %4d \n", q);
00690     strcat (lbuf, sbuf);
00691     sprintf (sbuf, "SSE       = %10.3f \n", sse_base);
00692     strcat (lbuf, sbuf);
00693     sprintf (sbuf, "MSE       = %10.3f \n", sse_base/(N-f-q));
00694     strcat (lbuf, sbuf);
00695  
00696 
00697      /*----- Statistical results for full model -----*/
00698     sprintf (sbuf, "\nFull Model: \n");
00699     strcat (lbuf, sbuf);
00700     sprintf (sbuf, "# params  = %4d \n", p);
00701     strcat (lbuf, sbuf);
00702     sprintf (sbuf, "SSE       = %10.3f \n", sse_full);
00703     strcat (lbuf, sbuf);
00704     sprintf (sbuf, "MSE       = %10.3f \n", sse_full/(N-f-p));
00705     strcat (lbuf, sbuf);
00706     
00707 
00708      /*----- Summary statistics -----*/
00709     sprintf (sbuf, "\nSummary Statistics: \n");
00710     strcat (lbuf, sbuf);
00711 
00712     sprintf (sbuf, "R^2       = %10.3f \n", rfull);
00713     strcat (lbuf, sbuf);
00714     
00715     sprintf (sbuf, "F[%2d,%3d] = %10.3f \n", p-q, N-f-p, ffull);
00716     strcat (lbuf, sbuf);
00717 
00718     pvalue = fstat_t2p ( (double) ffull, (double) p-q, (double) N-f-p);
00719     sprintf (sbuf, "p-value   = %e  \n", pvalue);
00720     strcat (lbuf, sbuf);
00721     
00722     
00723     *label = lbuf ;  /* send address of lbuf back in what label points to */
00724   }
00725   
00726 }
00727 
00728 
00729 /*---------------------------------------------------------------------------*/
00730 
00731 
00732 
00733 
00734 
00735 
00736