00373 {
00374   float * ts_array;        /* array of measured data for one voxel */
00375   float * coef;            /* regression parameters */
00376   float sse_base;          /* baseline model error sum of squares */
00377   float sse_full;          /* full model error sum of squares */
00378   float ffull;             /* full model F-statistic */
00379   float rfull;             /* full model R^2 stat. */
00380 
00381   int N;                   /* number of data points */
00382   int f;                   /* number of parameters removed by the filter */
00383   int p;                   /* number of parameters in the full model */
00384   int q;                   /* number of parameters in the baseline model */
00385   int n;                   /* data point index */
00386   int i;                   /* data point index */
00387   int ok = 1;
00388 
00389 
00390   /*----- Check initialization flag -----*/
00391   if (! plug_initialize)  return (0);
00392 
00393 
00394   /*----- Initialize local variables -----*/
00395   *NFirst = plug_NFirst;
00396 
00397   *NLast = plug_NLast;
00398   if (*NLast > nt-1)  *NLast = nt-1;
00399 
00400   N = *NLast - *NFirst + 1;
00401   N = powerof2(my_log2(N));
00402   *NLast = N + *NFirst - 1;
00403 
00404 
00405   plug_stop_filter = 
00406     FWT_1d_stop_filter (plug_num_stop_filters, plug_stop_band,
00407                         plug_stop_mintr, plug_stop_maxtr, *NFirst, N);
00408 
00409   plug_base_filter = 
00410     FWT_1d_pass_filter (plug_num_base_filters, plug_base_band,
00411                         plug_base_mintr, plug_base_maxtr, *NFirst, N);
00412 
00413   plug_sgnl_filter = 
00414     FWT_1d_pass_filter (plug_num_sgnl_filters, plug_sgnl_band,
00415                         plug_sgnl_mintr, plug_sgnl_maxtr, *NFirst, N);
00416 
00417   f = 0;
00418   for (i = 0;  i < N;  i++)
00419     if (plug_stop_filter[i] == 0.0)
00420       {
00421         f++;
00422         plug_base_filter[i] = 0.0;
00423         plug_sgnl_filter[i] = 0.0;
00424       }
00425 
00426   q = 0;
00427   for (i = 0;  i < N;  i++)
00428     if (plug_base_filter[i] == 1.0)
00429       {
00430         q++;
00431         plug_sgnl_filter[i] = 1.0;
00432       }
00433 
00434   p = 0;
00435   for (i = 0;  i < N;  i++)
00436     if (plug_sgnl_filter[i] == 1.0)
00437       {
00438         p++;
00439       }    
00440 
00441 
00442   /*----- Allocate memory for fitted time series and residuals -----*/
00443    coef     = (float *) malloc (sizeof(float) * p);    MTEST (coef);
00444   *coefts   = (float *) malloc (sizeof(float) * N);    MTEST (coefts);
00445   *fitts    = (float *) malloc (sizeof(float) * N);    MTEST (fitts);
00446   *sgnlts   = (float *) malloc (sizeof(float) * N);    MTEST (sgnlts);
00447   *errts    = (float *) malloc (sizeof(float) * N);    MTEST (errts);
00448 
00449 
00450   /*----- Extract data for this voxel -----*/
00451   ts_array = vec + *NFirst;
00452       
00453       
00454   /*----- Perform the wavelet analysis for this voxel-----*/
00455   wavelet_analysis (plug_wavelet_type, 
00456                     f, plug_stop_filter,
00457                     q, plug_base_filter, 
00458                     p, plug_sgnl_filter,
00459                     N, ts_array, coef,
00460                     &sse_base, &sse_full, &ffull, &rfull,
00461                     *coefts, *fitts, *sgnlts, *errts);
00462 
00463       
00464   /*----- Report results for this voxel -----*/
00465   printf ("\nResults for Voxel: \n");
00466   report_results (N, *NFirst, f, p, q, 
00467                   plug_base_filter, plug_sgnl_filter, 
00468                   coef, sse_base, sse_full, ffull, rfull, 
00469                   label);  
00470   printf ("%s \n", *label);
00471 
00472   plug_prev_nt = nt;
00473    
00474 
00475   /*----- Release memory -----*/
00476   free (plug_stop_filter);   plug_stop_filter = NULL;
00477   free (plug_base_filter);   plug_base_filter = NULL;
00478   free (plug_sgnl_filter);   plug_sgnl_filter = NULL;
00479   free (coef);               coef     = NULL;
00480 
00481 
00482   return (ok);
00483 }

00135 {
00136   PLUGIN_interface * plint ;     /* will be the output of this routine */
00137   int is;                        /* filter index */
00138 
00139 
00140   if( ncall > 0 ) return NULL ;  /* generate interface for ncall 0 */
00141 
00142 
00143   /***** do interface #0 *****/
00144 
00145   /*---------------- set titles and call point ----------------*/
00146 
00147   plint = PLUTO_new_interface ("Wavelets" ,
00148                                "Wavelet Analysis of Time Series Data" ,
00149                                helpstring, PLUGIN_CALL_VIA_MENU, WA_main);
00150 
00151   global_plint = plint ;  /* make global copy */
00152   
00153   PLUTO_add_hint (plint, 
00154                   "Control Wavelet Analysis Functions");
00155 
00156   PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00157 
00158   
00159   /*----- Initialize Global Variables -----*/
00160   for (is =0;  is < MAX_FILTERS;  is++)
00161     {
00162       plug_stop_band[is]  = 0;
00163       plug_stop_mintr[is] = 0.0;
00164       plug_stop_maxtr[is] = 0.0;
00165       plug_base_band[is]  = 0;
00166       plug_base_mintr[is] = 0.0;
00167       plug_base_maxtr[is] = 0.0;
00168       plug_sgnl_band[is]  = 0;
00169       plug_sgnl_mintr[is] = 0.0;
00170       plug_sgnl_maxtr[is] = 0.0;
00171     }
00172    
00173 
00174   /*----- Parameters -----*/
00175   PLUTO_add_option (plint, "Control",   "Control", TRUE);
00176   PLUTO_add_string (plint, "Wavelet",   MAX_WAVELET_TYPE,  WAVELET_TYPE_name, 
00177                     plug_wavelet_type);
00178   PLUTO_add_number (plint, "NFirst",    0, 32767, 0, plug_NFirst, TRUE);
00179   PLUTO_add_number (plint, "NLast",     0, 32767, 0, plug_NLast,  TRUE);
00180 
00181 
00182   /*----- Stop Filters -----*/
00183   for (is = 0;  is < MAX_FILTERS;  is++)
00184     {
00185       PLUTO_add_option (plint, "Filter", "Filter", FALSE);
00186       PLUTO_add_string (plint, "Type",   NFILTER,  filter_strings, 
00187                         plug_filter_type);
00188       PLUTO_add_number (plint, "Band",  -1, MAX_BAND, 0, 0, TRUE);
00189       PLUTO_add_number (plint, "Min TR", 0, 10000, 0, 0, TRUE);
00190       PLUTO_add_number (plint, "Max TR", 0, 10000, 0, 0, TRUE);
00191     }
00192 
00193 
00194   /*--------- done with interface setup ---------*/
00195   PLUTO_register_1D_funcstr ("WA_FWT",  WA_fwt);
00196   PLUTO_register_1D_funcstr ("WA_Fit",  WA_fit);
00197   PLUTO_register_1D_funcstr ("WA_Sgnl", WA_sgnl);
00198   PLUTO_register_1D_funcstr ("WA_Err",  WA_err);
00199 
00200 
00201   return plint ;
00202 }

00649 {
00650   int NFirst;              /* first image from input 3d+time dataset to use */
00651   int NLast;               /* last image from input 3d+time dataset to use */
00652   int n;                   /* time index */
00653   int ok;                  /* Boolean for successful calculation */
00654   float * coefts  = NULL;     /* forward wavelet transform coefficients */
00655   float * fitts   = NULL;     /* wavelet filtered time series */
00656   float * sgnlts  = NULL;     /* signal model fitted time series */
00657   float * errts   = NULL;     /* residual error time series */
00658 
00659 
00660   /*----- Calculate the wavelet filtered time series data -----*/
00661   ok = calculate_results (nt, vec, &NFirst, &NLast, label, 
00662                           &coefts, &fitts, &sgnlts, &errts);
00663   if (!ok)
00664     {
00665       for (n = 0;  n < nt;  n++)  vec[n] = 0.0;
00666       return;
00667     }
00668 
00669 
00670   /*----- Store the filtered time series residuals back into the array -----*/
00671   for (n = NFirst;  n <= NLast;  n++)
00672     {
00673       vec[n] = errts[n-NFirst];
00674     }
00675 
00676   for (n = 0;  n < NFirst;  n++)
00677     vec[n] = 0.0;
00678 
00679   for (n = NLast+1;  n < nt;  n++)
00680     vec[n] = 0.0;
00681     
00682 
00683   /*----- Deallocate memory -----*/
00684   free (coefts);   coefts = NULL;
00685   free (fitts);    fitts  = NULL;
00686   free (sgnlts);   sgnlts = NULL;
00687   free (errts);    errts  = NULL;
00688 
00689 
00690   return;
00691 }

00545 {
00546   int NFirst;              /* first image from input 3d+time dataset to use */
00547   int NLast;               /* last image from input 3d+time dataset to use */
00548   int n;                   /* time index */
00549   int ok;                  /* Boolean for successful calculation */
00550   float * coefts  = NULL;     /* forward wavelet transform coefficients */
00551   float * fitts   = NULL;     /* wavelet filtered time series */
00552   float * sgnlts  = NULL;     /* signal model fitted time series */
00553   float * errts   = NULL;     /* residual error time series */
00554 
00555 
00556   /*----- Calculate the wavelet filtered time series data -----*/
00557   ok = calculate_results (nt, vec, &NFirst, &NLast, label, 
00558                           &coefts, &fitts, &sgnlts, &errts);
00559   if (!ok)
00560     {
00561       for (n = 0;  n < nt;  n++)  vec[n] = 0.0;
00562       return;
00563     }
00564 
00565 
00566   /*----- Store the filtered time series data back into the array -----*/
00567   for (n = NFirst;  n <= NLast;  n++)
00568     {
00569       vec[n] = fitts[n-NFirst];
00570     }
00571 
00572   for (n = 0;  n < NFirst;  n++)
00573     vec[n] = vec[NFirst];
00574 
00575   for (n = NLast+1;  n < nt;  n++)
00576     vec[n] = vec[NLast];
00577     
00578 
00579   /*----- Deallocate memory -----*/
00580   free (coefts);   coefts = NULL;
00581   free (fitts);    fitts  = NULL;
00582   free (sgnlts);   sgnlts = NULL;
00583   free (errts);    errts  = NULL;
00584 
00585 
00586   return;
00587 }

00493 {
00494   int NFirst;              /* first image from input 3d+time dataset to use */
00495   int NLast;               /* last image from input 3d+time dataset to use */
00496   int n;                   /* time index */
00497   int ok;                  /* Boolean for successful calculation */
00498   float * coefts  = NULL;     /* forward wavelet transform coefficients */
00499   float * fitts   = NULL;     /* wavelet filtered time series */
00500   float * sgnlts  = NULL;     /* signal model fitted time series */
00501   float * errts   = NULL;     /* residual error time series */
00502 
00503 
00504   /*----- Calculate the wavelet filtered time series data -----*/
00505   ok = calculate_results (nt, vec, &NFirst, &NLast, label, 
00506                           &coefts, &fitts, &sgnlts, &errts);
00507   if (!ok)
00508     {
00509       for (n = 0;  n < nt;  n++)  vec[n] = 0.0;
00510       return;
00511     }
00512 
00513 
00514   /*----- Store the filtered time series data back into the array -----*/
00515   for (n = NFirst;  n <= NLast;  n++)
00516     {
00517       vec[n] = coefts[n-NFirst];
00518     }
00519 
00520   for (n = 0;  n < NFirst;  n++)
00521     vec[n] = 0.0;
00522 
00523   for (n = NLast+1;  n < nt;  n++)
00524     vec[n] = 0.0;
00525     
00526 
00527   /*----- Deallocate memory -----*/
00528   free (coefts);   coefts = NULL;
00529   free (fitts);    fitts  = NULL;
00530   free (sgnlts);   sgnlts = NULL;
00531   free (errts);    errts  = NULL;
00532 
00533 
00534   return;
00535 }

00213 {
00214   char * str;                           /* input string */
00215   int is;                               /* filter index */
00216   
00217 
00218   /*----- reset flag for successful initialization -----*/
00219   plug_initialize = 0;
00220     
00221 
00222   /*--------- go to Control input line ---------*/
00223   PLUTO_next_option (plint);
00224   str    = PLUTO_get_string (plint);
00225   plug_wavelet_type = PLUTO_string_index (str, MAX_WAVELET_TYPE, 
00226                                           WAVELET_TYPE_name);
00227   plug_NFirst = PLUTO_get_number (plint);
00228   plug_NLast  = PLUTO_get_number (plint);
00229 
00230 
00231   /*------ read input line(s) -----*/
00232   plug_num_stop_filters = 0;
00233   plug_num_base_filters = 0;
00234   plug_num_sgnl_filters = 0;
00235 
00236   do
00237     {
00238       str = PLUTO_get_optiontag(plint); 
00239       if (str == NULL)  break;
00240       if (strcmp (str, "Filter") != 0)
00241         return "************************\n"
00242                "Illegal optiontag found!\n"
00243                "************************";
00244      
00245 
00246       /*----- Read Filter Specification Line -----*/
00247       str    = PLUTO_get_string (plint);
00248       plug_filter_type = PLUTO_string_index (str, NFILTER, filter_strings);
00249       
00250       switch (plug_filter_type)
00251         {
00252         case 0:
00253           {
00254             plug_stop_band[plug_num_stop_filters] = PLUTO_get_number(plint);   
00255             plug_stop_mintr[plug_num_stop_filters] 
00256               = PLUTO_get_number(plint);
00257             plug_stop_maxtr[plug_num_stop_filters] 
00258               = PLUTO_get_number(plint);
00259           
00260             if (plug_stop_mintr[plug_num_stop_filters] > 
00261                 plug_stop_maxtr[plug_num_stop_filters])
00262               return "*************************\n"
00263                 "Require Min TR <= Max TR \n"
00264                 "*************************"  ;
00265             
00266             plug_num_stop_filters++;
00267             break;
00268           }
00269 
00270         case 1:
00271           {
00272             plug_base_band[plug_num_base_filters] = PLUTO_get_number(plint);   
00273             plug_base_mintr[plug_num_base_filters] 
00274               = PLUTO_get_number(plint);
00275             plug_base_maxtr[plug_num_base_filters] 
00276               = PLUTO_get_number(plint);
00277           
00278             if (plug_base_mintr[plug_num_base_filters] > 
00279                 plug_base_maxtr[plug_num_base_filters])
00280               return "*************************\n"
00281                 "Require Min TR <= Max TR \n"
00282                 "*************************"  ;
00283             
00284             plug_num_base_filters++;
00285             break;
00286           }
00287 
00288         case 2:
00289           {
00290             plug_sgnl_band[plug_num_sgnl_filters]=PLUTO_get_number(plint); 
00291             plug_sgnl_mintr[plug_num_sgnl_filters]
00292               = PLUTO_get_number(plint);
00293             plug_sgnl_maxtr[plug_num_sgnl_filters]
00294               = PLUTO_get_number(plint);
00295           
00296             if (plug_sgnl_mintr[plug_num_sgnl_filters] > 
00297                 plug_sgnl_maxtr[plug_num_sgnl_filters])
00298               return "*************************\n"
00299                 "Require Min TR <= Max TR \n"
00300                 "*************************"  ;
00301             
00302             plug_num_sgnl_filters++;
00303             break;
00304           }
00305 
00306         }
00307     }
00308   while (1);
00309 
00310 
00311   /*----- Identify software -----*/
00312   printf ("\n\n");
00313   printf ("Program: %s \n", PROGRAM_NAME);
00314   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
00315   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00316   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00317   printf ("\n");
00318 
00319   /*----- show current input options -----*/
00320   printf ("\nControls: \n");
00321   printf ("Wavelet Type = %10s \n", WAVELET_TYPE_name[plug_wavelet_type]);
00322   printf ("NFirst       = %10d \n", plug_NFirst);
00323   printf ("NLast        = %10d \n", plug_NLast);
00324 
00325   for (is = 0;  is < plug_num_stop_filters;  is++)
00326     {
00327       printf ("\nStop Filter:       Band = %4d   ", plug_stop_band[is]);
00328       printf ("Min. TR = %4d   Max. TR = %4d \n", 
00329               plug_stop_mintr[is], plug_stop_maxtr[is]);
00330     }
00331  
00332   for (is = 0;  is < plug_num_base_filters;  is++)
00333     {
00334       printf ("\nBaseline Filter:   Band = %4d   ", plug_base_band[is]);
00335       printf ("Min. TR = %4d   Max. TR = %4d \n", 
00336               plug_base_mintr[is], plug_base_maxtr[is]);
00337     }
00338  
00339   for (is = 0;  is < plug_num_sgnl_filters;  is++)
00340     {
00341       printf ("\nSignal Filter:     Band = %4d   ", plug_sgnl_band[is]);
00342       printf ("Min. TR = %4d   Max. TR = %4d \n", 
00343               plug_sgnl_mintr[is], plug_sgnl_maxtr[is]);
00344     }
00345  
00346 
00347   /*--- nothing left to do until data arrives ---*/
00348   plug_initialize = 1 ;  /* successful initialization */
00349   plug_prev_nt = 0;      /* previous time series length */
00350   
00351   return NULL ;
00352 }

00597 {
00598   int NFirst;              /* first image from input 3d+time dataset to use */
00599   int NLast;               /* last image from input 3d+time dataset to use */
00600   int n;                   /* time index */
00601   int ok;                  /* Boolean for successful calculation */
00602   float * coefts  = NULL;     /* forward wavelet transform coefficients */
00603   float * fitts   = NULL;     /* wavelet filtered time series */
00604   float * sgnlts  = NULL;     /* signal model fitted time series */
00605   float * errts   = NULL;     /* residual error time series */
00606 
00607 
00608   /*----- Calculate the wavelet filtered time series data -----*/
00609   ok = calculate_results (nt, vec, &NFirst, &NLast, label, 
00610                           &coefts, &fitts, &sgnlts, &errts);
00611   if (!ok)
00612     {
00613       for (n = 0;  n < nt;  n++)  vec[n] = 0.0;
00614       return;
00615     }
00616 
00617 
00618   /*----- Store the filtered time series residuals back into the array -----*/
00619   for (n = NFirst;  n <= NLast;  n++)
00620     {
00621       vec[n] = sgnlts[n-NFirst];
00622     }
00623 
00624   for (n = 0;  n < NFirst;  n++)
00625     vec[n] = 0.0;
00626 
00627   for (n = NLast+1;  n < nt;  n++)
00628     vec[n] = 0.0;
00629     
00630 
00631   /*----- Deallocate memory -----*/
00632   free (coefts);   coefts = NULL;
00633   free (fitts);    fitts  = NULL;
00634   free (sgnlts);   sgnlts = NULL;
00635   free (errts);    errts  = NULL;
00636 
00637 
00638   return;
00639 }

   " Purpose:    Wavelet Analysis of FMRI time series data.                 \n"
   "                                                                        \n"
   " Control:    Wavelet    = Type of wavelet to be used in analysis        \n"
   "             NFirst     = Number of first time series point to use.     \n"
   "             NLast      = Number of last time series point to use.      \n"
   "                                                                        \n"
   " Filter:     Type       = Type of filtering to apply in this            \n"
   "                          time-frequency window:                        \n"
   "               Stop     = set wavelet coefficients to zero              \n"
   "               Baseline = assign wavelet coefficients to baseline model \n"
   "               Signal   = assign wavelet coefficients to signal model   \n"
   "                                                                        \n"
   "             Band       = Frequency band at which to apply filtering.   \n"
   "             Min TR     = Minimum value for time window (in TR)         \n"
   "             Max TR     = Maximum value for time window (in TR)         \n"
   "                                                                        \n"
   "For more information, see 'Wavelet Analysis of FMRI Time Series Data'   \n"
   "which is contained in file 3dWavelets.ps of the AFNI distribution.      \n"