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 /*---------------------------------------------------------------------------*/
00008 /*
00009   Plugin to perform wavelet analysis of time series data.
00010 
00011   File:    plug_wavelet.c
00012   Author:  B. Douglas Ward
00013   Date:    28 March 2000
00014 
00015   Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
00016   Date:    02 December 2002
00017 */
00018 
00019 
00020 /*---------------------------------------------------------------------------*/
00021 /*
00022   Software identification.
00023 */
00024 
00025 #define PROGRAM_NAME "plug_wavelets"                 /* name of this program */
00026 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00027 #define PROGRAM_INITIAL "28 March 2000"      /* initial program release date */
00028 #define PROGRAM_LATEST  "02 December  2002"   /* latest program revision date*/
00029 
00030 /*---------------------------------------------------------------------------*/
00031 /*
00032   Include header files and source code.
00033 */
00034 
00035 #include "afni.h"
00036 #include "Wavelets.h"
00037 #include "Wavelets.c"
00038 
00039 /*---------------------------------------------------------------------------*/
00040 /*
00041   Global constants. 
00042 */
00043 
00044 #define MAX_NAME_LENGTH THD_MAX_NAME    /* max. string length for file names */
00045 #define MAX_FILTERS 20                  /* max. number of blocking filters */
00046 #define MAX_BAND 20                     /* max. frequency band */
00047 
00048 /*---------------------------------------------------------------------------*/
00049 
00050 #ifndef ALLOW_PLUGINS
00051 #  error "Plugins not properly set up -- see machdep.h"
00052 #endif
00053 
00054 
00055 /*------------- string to 'help' the user -------------*/
00056 
00057 static char helpstring[] =
00058    " Purpose:    Wavelet Analysis of FMRI time series data.                 \n"
00059    "                                                                        \n"
00060    " Control:    Wavelet    = Type of wavelet to be used in analysis        \n"
00061    "             NFirst     = Number of first time series point to use.     \n"
00062    "             NLast      = Number of last time series point to use.      \n"
00063    "                                                                        \n"
00064    " Filter:     Type       = Type of filtering to apply in this            \n"
00065    "                          time-frequency window:                        \n"
00066    "               Stop     = set wavelet coefficients to zero              \n"
00067    "               Baseline = assign wavelet coefficients to baseline model \n"
00068    "               Signal   = assign wavelet coefficients to signal model   \n"
00069    "                                                                        \n"
00070    "             Band       = Frequency band at which to apply filtering.   \n"
00071    "             Min TR     = Minimum value for time window (in TR)         \n"
00072    "             Max TR     = Maximum value for time window (in TR)         \n"
00073    "                                                                        \n"
00074    "For more information, see 'Wavelet Analysis of FMRI Time Series Data'   \n"
00075    "which is contained in file 3dWavelets.ps of the AFNI distribution.      \n"
00076 ;
00077 
00078 
00079 /*------- Strings for filter type ------*/
00080 #define NFILTER 3
00081 static char * filter_strings[NFILTER] = {"Stop",  "Baseline", "Signal"};
00082 
00083 
00084 /*--------------- prototypes for internal routines ---------------*/
00085 
00086 char * WA_main( PLUGIN_interface * ) ;  /* the entry point */
00087 
00088 void WA_fwt  (int nt, double to, double dt, float * vec, char ** label);
00089 void WA_fit  (int nt, double to, double dt, float * vec, char ** label);
00090 void WA_sgnl (int nt, double to, double dt, float * vec, char ** label);
00091 void WA_err  (int nt, double to, double dt, float * vec, char ** label);
00092 
00093 
00094 
00095 /*---------------- global data -------------------*/
00096 
00097 static PLUGIN_interface * global_plint = NULL;
00098 
00099 static int plug_wavelet_type=0;        /* type of wavelet to use in analysis */
00100 static int plug_NFirst=0;   /* first image from input 3d+time dataset to use */
00101 static int plug_NLast=32767; /* last image from input 3d+time dataset to use */
00102 static int plug_initialize=0;             /* flag for perform initialization */
00103 static int plug_prev_nt=0;                    /* previous time series length */
00104 static int plug_filter_type=0;          /* type of filter to use in analysis */
00105 
00106 static int plug_num_stop_filters = 0;    /* number of stop filters */
00107 static int plug_stop_band[MAX_FILTERS];  /* freq band for stop filter */
00108 static int plug_stop_mintr[MAX_FILTERS]; /* min time for stop filter */
00109 static int plug_stop_maxtr[MAX_FILTERS]; /* max time for stop filter */
00110 static float * plug_stop_filter = NULL;  /* select wavelet coefs. to stop */
00111 
00112 static int plug_num_base_filters = 0;    /* number of base filters */
00113 static int plug_base_band[MAX_FILTERS];  /* freq band for base filter */
00114 static int plug_base_mintr[MAX_FILTERS]; /* min time for base filter */
00115 static int plug_base_maxtr[MAX_FILTERS]; /* max time for base filter */
00116 static float * plug_base_filter = NULL;  /* select wavelet coefs. 
00117                                             for baseline */
00118 
00119 static int plug_num_sgnl_filters = 0;    /* number of signal filters */
00120 static int plug_sgnl_band[MAX_FILTERS];  /* freq band for signal filter */
00121 static int plug_sgnl_mintr[MAX_FILTERS]; /* min time for signal filter */
00122 static int plug_sgnl_maxtr[MAX_FILTERS]; /* max time for signal filter */
00123 static float * plug_sgnl_filter = NULL;  /* select wavelet coefs. for signal */
00124 
00125 
00126 /*---------------------------------------------------------------------------*/
00127 /*
00128    Set up the interface to the user.
00129 */
00130 
00131 
00132 DEFINE_PLUGIN_PROTOTYPE
00133 
00134 PLUGIN_interface * PLUGIN_init( int ncall )
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 }
00203 
00204 
00205 /*---------------------------------------------------------------------------*/
00206 /*
00207   Main routine for this plugin (will be called from AFNI).
00208   If the return string is not NULL, some error transpired, and
00209   AFNI will popup the return string in a message box.
00210 */
00211 
00212 char * WA_main( PLUGIN_interface * plint )
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 }
00353 
00354 
00355 /*---------------------------------------------------------------------------*/
00356 /*
00357   Perform the wavelet analysis.
00358 */
00359 
00360 int calculate_results 
00361 (
00362   int nt,               /* number of time points */
00363   float * vec,          /* input measured time series data */
00364   int * NFirst,         /* first image from input 3d+time dataset to use */
00365   int * NLast,          /* last image from input 3d+time dataset to use */
00366   char ** label,        /* string containing statistical summary of results */
00367   float ** coefts,      /* forward wavelet transform coefficients */
00368   float ** fitts,       /* full model fit to the time series data */
00369   float ** sgnlts,      /* signal model fit to the time series data */
00370   float ** errts        /* residual error time series */
00371 )
00372   
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 }
00484 
00485 
00486 /*---------------------------------------------------------------------------*/
00487 /*
00488   Calculate the forward wavelet transform coefficients.
00489 */
00490 
00491 void WA_fwt (int nt, double to, double dt, float * vec, char ** label)
00492 
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 }
00536 
00537 
00538 /*---------------------------------------------------------------------------*/
00539 /*
00540   Calculate the full model wavelet fit to the time series data.
00541 */
00542 
00543 void WA_fit (int nt, double to, double dt, float * vec, char ** label)
00544 
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 }
00588 
00589 
00590 /*---------------------------------------------------------------------------*/
00591 /*
00592   Calculate the signal model wavelet fit to the time series data.
00593 */
00594 
00595 void WA_sgnl (int nt, double to, double dt, float * vec, char ** label)
00596 
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 }
00640 
00641 
00642 /*---------------------------------------------------------------------------*/
00643 /*
00644   Calculate the residual error from the full model wavelet fit.
00645 */
00646 
00647 void WA_err (int nt, double to, double dt, float * vec, char ** label)
00648 
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 }
00692 
00693 
00694 /*---------------------------------------------------------------------------*/
00695 
00696 
00697 
00698 
00699 
00700 
00701 
00702 
00703 
00704 
00705 
00706 
00707 
00708