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 }

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 }

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 }

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 }

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 }

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 }

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 }

00111 {
00112   int i;
00113    
00114   i = floor (log(n)/log(2.0) + 1.0e-10);
00115 
00116   return (i);
00117 }

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 }

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 }

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 }

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 }

00034 {
00035   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00036   exit(1);
00037 }

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 }