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 #include "mrilib.h"
00008 
00009 /*------- Adapted from plug_stats.c ---------*/
00010 
00011 #define METH_MEAN   0
00012 #define METH_SLOPE  1
00013 #define METH_SIGMA  2
00014 #define METH_CVAR   3
00015 
00016 #define METH_MEDIAN 4   /* 14 Feb 2001 */
00017 #define METH_MAD    5
00018 
00019 #define METH_MAX    6
00020 #define METH_MIN    7
00021 
00022 #define METH_DW     8   /* KRH 3 Dec 2002 */
00023 
00024 #define METH_SIGMA_NOD     9   /* KRH 27 Dec 2002 */
00025 #define METH_CVAR_NOD     10   /* KRH 27 Dec 2002 */
00026 
00027 #define METH_AUTOCORR     11  /* KRH 16 Jun 2004 */
00028 #define METH_AUTOREGP     12  /* KRH 16 Jun 2004 */
00029 
00030 #define METH_ABSMAX     13  /* KRH 15 Feb 2005 */
00031 
00032 #define METH_ARGMAX     14  /* KRH 4 Aug 2005 */
00033 #define METH_ARGMIN     15  /* KRH 4 Aug 2005 */
00034 #define METH_ARGABSMAX     16  /* KRH 4 Aug 2005 */
00035 
00036 #define MAX_NUM_OF_METHS 20
00037 static int meth[MAX_NUM_OF_METHS] = {METH_MEAN};
00038 static int nmeths = 0;
00039 static char prefix[THD_MAX_PREFIX] = "stat" ;
00040 static int datum                   = MRI_float ;
00041 static char meth_names[][20] = {"Mean","Slope","Std Dev","Coef of Var","Median",
00042                             "Med Abs Dev", "Max", "Min", "Durbin-Watson", "Std Dev (NOD)",
00043                             "Coef Var(NOD)","AutoCorr","AutoReg","Absolute Max",
00044                             "ArgMax","ArgMin","ArgAbsMax"};
00045 static void STATS_tsfunc( double tzero , double tdelta ,
00046                          int npts , float ts[] , double ts_mean ,
00047                          double ts_slope , void * ud , int nbriks, float * val ) ;
00048 
00049 static void autocorr( int npts, float ints[], int numVals, float outcoeff[] ) ;
00050 
00051 int main( int argc , char * argv[] )
00052 {
00053    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00054    int nopt, nbriks, ii ;
00055    int addBriks = 0;
00056    int numMultBriks,methIndex,brikIndex;
00057 
00058    /*----- Read command line -----*/
00059    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00060       printf("Usage: 3dTstat [options] dataset\n"
00061              "Computes one or more voxel-wise statistics for a 3D+time dataset\n"
00062              "and stores them in a bucket dataset.\n"
00063              "\n"
00064              "Options:\n"
00065              " -mean   = compute mean of input voxels [DEFAULT]\n"
00066              " -slope  = compute mean slope of input voxels vs. time\n"
00067              " -stdev  = compute standard deviation of input voxels\n"
00068              "             [N.B.: this is computed after    ]\n"
00069              "             [      the slope has been removed]\n"
00070              " -cvar   = compute coefficient of variation of input\n"
00071              "             voxels = stdev/fabs(mean)\n"
00072              "   **N.B.: You can add NOD to the end of the above 2\n"
00073              "           options to turn off detrending, as in\n"
00074              "             -stdevNOD or -cvarNOD\n"
00075              "\n"
00076              " -MAD    = compute MAD (median absolute deviation) of\n"
00077              "             input voxels = median(|voxel-median(voxel)|)\n"
00078              "             [N.B.: the trend is NOT removed for this]\n"
00079              " -DW    = compute Durbin-Watson Statistic of\n"
00080              "             input voxels\n"
00081              "             [N.B.: the trend is removed for this]\n"
00082              " -median = compute median of input voxels  [undetrended]\n"
00083              " -min    = compute minimum of input voxels [undetrended]\n"
00084              " -max    = compute maximum of input voxels [undetrended]\n"
00085              " -absmax    = compute absolute maximum of input voxels [undetrended]\n"
00086              " -argmin    = index of minimum of input voxels [undetrended]\n"
00087              " -argmax    = index of maximum of input voxels [undetrended]\n"
00088              " -argabsmax    = index of absolute maximum of input voxels [undetrended]\n"
00089              "\n"
00090              " -prefix p = use string 'p' for the prefix of the\n"
00091              "               output dataset [DEFAULT = 'stat']\n"
00092              " -datum d  = use data type 'd' for the type of storage\n"
00093              "               of the output, where 'd' is one of\n"
00094              "               'byte', 'short', or 'float' [DEFAULT=float]\n"
00095              " -autocorr n = compute autocorrelation function and return\n"
00096              "               first n coefficients\n"
00097              " -autoreg n = compute autoregression coefficients and return\n"
00098              "               first n coefficients\n"
00099              "    [N.B.: -autocorr 0 and/or -autoreg 0 will return coefficients\n"
00100              "           equal to the length of the input data\n"
00101              "\n"
00102              "The output is a bucket dataset.  The input dataset\n"
00103              "may use a sub-brick selection list, as in program 3dcalc.\n"
00104            ) ;
00105       printf("\n" MASTER_SHORTHELP_STRING ) ;
00106       exit(0) ;
00107    }
00108 
00109    mainENTRY("3dTstat main"); machdep(); AFNI_logger("3dTstat",argc,argv);
00110    PRINT_VERSION("3dTstat");
00111 
00112    nopt = 1 ;
00113    nbriks = 0 ;
00114    nmeths = 0 ;
00115    while( nopt < argc && argv[nopt][0] == '-' ){
00116 
00117       /*-- methods --*/
00118 
00119       if( strcmp(argv[nopt],"-median") == 0 ){
00120          meth[nmeths++] = METH_MEDIAN ;
00121          nbriks++ ;
00122          nopt++ ; continue ;
00123       }
00124 
00125       if( strcmp(argv[nopt],"-DW") == 0 ){
00126          meth[nmeths++] = METH_DW ;
00127          nbriks++ ;
00128          nopt++ ; continue ;
00129       }
00130 
00131       if( strcmp(argv[nopt],"-MAD") == 0 ){
00132          meth[nmeths++] = METH_MAD ;
00133          nbriks++ ;
00134          nopt++ ; continue ;
00135       }
00136 
00137       if( strcmp(argv[nopt],"-mean") == 0 ){
00138          meth[nmeths++] = METH_MEAN ;
00139          nbriks++ ;
00140          nopt++ ; continue ;
00141       }
00142 
00143       if( strcmp(argv[nopt],"-slope") == 0 ){
00144          meth[nmeths++] = METH_SLOPE ;
00145          nbriks++ ;
00146          nopt++ ; continue ;
00147       }
00148 
00149       if( strcmp(argv[nopt],"-stdev") == 0 ||
00150           strcmp(argv[nopt],"-sigma") == 0   ){
00151 
00152          meth[nmeths++] = METH_SIGMA ;
00153          nbriks++ ;
00154          nopt++ ; continue ;
00155       }
00156 
00157       if( strcmp(argv[nopt],"-cvar") == 0 ){
00158          meth[nmeths++] = METH_CVAR ;
00159          nbriks++ ;
00160          nopt++ ; continue ;
00161       }
00162 
00163       if( strcmp(argv[nopt],"-stdevNOD") == 0 ||
00164           strcmp(argv[nopt],"-sigmaNOD") == 0   ){  /* 07 Dec 2001 */
00165 
00166          meth[nmeths++] = METH_SIGMA_NOD ;
00167          nbriks++ ;
00168          nopt++ ; continue ;
00169       }
00170 
00171       if( strcmp(argv[nopt],"-cvarNOD") == 0 ){     /* 07 Dec 2001 */
00172          meth[nmeths++] = METH_CVAR_NOD ;
00173          nbriks++ ;
00174          nopt++ ; continue ;
00175       }
00176 
00177       if( strcmp(argv[nopt],"-min") == 0 ){
00178          meth[nmeths++] = METH_MIN ;
00179          nbriks++ ;
00180          nopt++ ; continue ;
00181       }
00182 
00183       if( strcmp(argv[nopt],"-max") == 0 ){
00184          meth[nmeths++] = METH_MAX ;
00185          nbriks++ ;
00186          nopt++ ; continue ;
00187       }
00188 
00189       if( strcmp(argv[nopt],"-absmax") == 0 ){
00190          meth[nmeths++] = METH_ABSMAX ;
00191          nbriks++ ;
00192          nopt++ ; continue ;
00193       }
00194 
00195       if( strcmp(argv[nopt],"-argmin") == 0 ){
00196          meth[nmeths++] = METH_ARGMIN ;
00197          nbriks++ ;
00198          nopt++ ; continue ;
00199       }
00200 
00201       if( strcmp(argv[nopt],"-argmax") == 0 ){
00202          meth[nmeths++] = METH_ARGMAX ;
00203          nbriks++ ;
00204          nopt++ ; continue ;
00205       }
00206 
00207       if( strcmp(argv[nopt],"-argabsmax") == 0 ){
00208          meth[nmeths++] = METH_ARGABSMAX ;
00209          nbriks++ ;
00210          nopt++ ; continue ;
00211       }
00212 
00213       if( strcmp(argv[nopt],"-autocorr") == 0 ){
00214          meth[nmeths++] = METH_AUTOCORR ;
00215          if( ++nopt >= argc ){
00216             fprintf(stderr,"*** -autocorr needs an argument!\n"); exit(1);
00217          }
00218          meth[nmeths++] = atoi(argv[nopt++]);
00219          if (meth[nmeths - 1] == 0) {
00220            addBriks++;
00221          } else {
00222            nbriks+=meth[nmeths - 1] ;
00223          }
00224          continue ;
00225       }
00226 
00227       if( strcmp(argv[nopt],"-autoreg") == 0 ){
00228          meth[nmeths++] = METH_AUTOREGP ;
00229          if( ++nopt >= argc ){
00230             fprintf(stderr,"*** -autoreg needs an argument!\n"); exit(1);
00231          }
00232          meth[nmeths++] = atoi(argv[nopt++]);
00233          if (meth[nmeths - 1] == 0) {
00234            addBriks++;
00235          } else {
00236            nbriks+=meth[nmeths - 1] ;
00237          }
00238          continue ;
00239       }
00240 
00241       /*-- prefix --*/
00242 
00243       if( strcmp(argv[nopt],"-prefix") == 0 ){
00244          if( ++nopt >= argc ){
00245             fprintf(stderr,"*** -prefix needs an argument!\n"); exit(1);
00246          }
00247          MCW_strncpy(prefix,argv[nopt],THD_MAX_PREFIX) ;
00248          if( !THD_filename_ok(prefix) ){
00249             fprintf(stderr,"*** %s is not a valid prefix!\n",prefix); exit(1);
00250          }
00251          nopt++ ; continue ;
00252       }
00253 
00254       /*-- datum --*/
00255 
00256       if( strcmp(argv[nopt],"-datum") == 0 ){
00257          if( ++nopt >= argc ){
00258             fprintf(stderr,"*** -datum needs an argument!\n"); exit(1);
00259          }
00260          if( strcmp(argv[nopt],"short") == 0 ){
00261             datum = MRI_short ;
00262          } else if( strcmp(argv[nopt],"float") == 0 ){
00263             datum = MRI_float ;
00264          } else if( strcmp(argv[nopt],"byte") == 0 ){
00265             datum = MRI_byte ;
00266          } else {
00267             fprintf(stderr,"-datum of type '%s' is not supported!\n",
00268                     argv[nopt] ) ;
00269             exit(1) ;
00270          }
00271          nopt++ ; continue ;
00272       }
00273 
00274       /*-- Quien sabe'? --*/
00275 
00276       fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00277    }
00278 
00279    /*--- If no options selected, default to single stat MEAN--KRH---*/
00280 
00281    if (nmeths == 0) nmeths = 1;
00282    if (nbriks == 0 && addBriks == 0) nbriks = 1;
00283 
00284    /*----- read input dataset -----*/
00285 
00286    if( nopt >= argc ){
00287       fprintf(stderr,"*** No input dataset!?\n"); exit(1);
00288    }
00289 
00290    old_dset = THD_open_dataset( argv[nopt] ) ;
00291    if( !ISVALID_DSET(old_dset) ){
00292       fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]); exit(1);
00293    }
00294 
00295    if( DSET_NVALS(old_dset) < 2 ){
00296       fprintf(stderr,"*** Can't use dataset with < 2 values per voxel!\n") ;
00297       exit(1) ;
00298    }
00299 
00300    if( DSET_NUM_TIMES(old_dset) < 2 ){
00301       fprintf(stderr,"--- Input dataset is not 3D+time!\n"
00302                      "--- Adding an artificial time axis with dt=1.0\n" ) ;
00303       EDIT_dset_items( old_dset ,
00304                           ADN_ntt    , DSET_NVALS(old_dset) ,
00305                           ADN_ttorg  , 0.0 ,
00306                           ADN_ttdel  , 1.0 ,
00307                           ADN_tunits , UNITS_SEC_TYPE ,
00308                        NULL ) ;
00309    }
00310 
00311    /* If one or more of the -autocorr/-autoreg options was called with */
00312    /* an argument of 0, then I'll now add extra BRIKs for the N-1 data */
00313    /* output points for each.                                          */
00314    nbriks += ((DSET_NVALS(old_dset)-1) * addBriks);
00315 
00316    /*------------- ready to compute new dataset -----------*/
00317 
00318    new_dset = MAKER_4D_to_typed_fbuc(
00319                  old_dset ,             /* input dataset */
00320                  prefix ,               /* output prefix */
00321                  datum ,                /* output datum  */
00322                  0 ,                    /* ignore count  */
00323                  0 ,              /* can't detrend in maker function  KRH 12/02*/
00324                  nbriks ,               /* number of briks */
00325                  STATS_tsfunc ,         /* timeseries processor */
00326                  NULL                   /* data for tsfunc */
00327               ) ;
00328 
00329    if( new_dset != NULL ){
00330       tross_Copy_History( old_dset , new_dset ) ;
00331       tross_Make_History( "3dTstat" , argc,argv , new_dset ) ;
00332       for (methIndex = 0,brikIndex = 0; methIndex < nmeths; methIndex++, brikIndex++) {
00333         if ((meth[methIndex] == METH_AUTOCORR) || (meth[methIndex] == METH_AUTOREGP)) {
00334           numMultBriks = meth[methIndex+1];
00335           if (numMultBriks == 0) numMultBriks = DSET_NVALS(old_dset);
00336           for (ii = 1; ii <= numMultBriks; ii++) {
00337             char tmpstr[25];
00338             if (meth[methIndex] == METH_AUTOREGP) {
00339               sprintf(tmpstr,"%s[%d](%d)",meth_names[meth[methIndex]],numMultBriks,ii);
00340             } else {
00341               sprintf(tmpstr,"%s(%d)",meth_names[meth[methIndex]],ii);
00342             }
00343             EDIT_BRICK_LABEL(new_dset, (brikIndex + ii - 1), tmpstr) ;
00344           }
00345           methIndex++;
00346           brikIndex += numMultBriks - 1;
00347         } else {
00348           EDIT_BRICK_LABEL(new_dset, brikIndex, meth_names[meth[methIndex]]) ;
00349         }
00350       }
00351       DSET_write( new_dset ) ;
00352       fprintf(stderr,"--- Output dataset %s\n",DSET_BRIKNAME(new_dset)) ;
00353    } else {
00354       fprintf(stderr,"*** Unable to compute output dataset!\n") ;
00355       exit(1) ;
00356    }
00357 
00358    exit(0) ;
00359 }
00360 
00361 /**********************************************************************
00362    Function that does the real work
00363 ***********************************************************************/
00364 
00365 static void STATS_tsfunc( double tzero, double tdelta ,
00366                           int npts, float ts[],
00367                           double ts_mean, double ts_slope,
00368                           void * ud, int nbriks, float * val          )
00369 {
00370    static int nvox , ncall ;
00371    int meth_index, ii , out_index;
00372    float* ts_det;
00373 
00374    /** is this a "notification"? **/
00375 
00376    if( val == NULL ){
00377 
00378       if( npts > 0 ){  /* the "start notification" */
00379 
00380          nvox  = npts ;                       /* keep track of   */
00381          ncall = 0 ;                          /* number of calls */
00382 
00383       } else {  /* the "end notification" */
00384 
00385          /* nothing to do here */
00386 
00387       }
00388       return ;
00389    }
00390 
00391    /* KRH make copy and detrend it right here.  Use ts_mean and
00392     * ts_slope to detrend */
00393 
00394    ts_det = (float*)calloc(npts, sizeof(float));
00395    memcpy( ts_det, ts, npts * sizeof(float));
00396    for( ii = 0; ii < npts; ii++) ts_det[ii] -= 
00397            (ts_mean - (ts_slope * (npts - 1) * tdelta/2) + ts_slope * tdelta * ii) ;
00398    
00399    /** OK, actually do some work **/
00400 
00401    /* This main loop now uses meth_index AND out_index as loop variables, mainly */
00402    /* to avoid rewriting code that worked.                                       */
00403 
00404    /* meth_index is an index into the static method array, which contains the    */
00405    /* list of methods to be executed (and will also contain an integer           */
00406    /* parameter specifying the number of return values if -autocorr n and/or     */
00407    /* -autoreg p have been requested).                                           */
00408 
00409    /* out_index is an index into the output array.                               */
00410    for (meth_index = 0, out_index = 0 ; meth_index < nmeths; meth_index++, out_index++) {
00411    switch( meth[meth_index] ){
00412 
00413       default:
00414       case METH_MEAN:  val[out_index] = ts_mean  ; break ;
00415 
00416       case METH_SLOPE: val[out_index] = ts_slope ; break ;
00417 
00418       case METH_CVAR_NOD:
00419       case METH_SIGMA_NOD:
00420       case METH_CVAR:
00421       case METH_SIGMA:{
00422          register int ii ;
00423          register double sum ;
00424 
00425          sum = 0.0 ;
00426          if((meth[meth_index] == METH_CVAR) || (meth[meth_index] == METH_SIGMA )){
00427            for( ii=0 ; ii < npts ; ii++ ) sum += ts_det[ii] * ts_det[ii] ;
00428          } else {
00429            for( ii=0 ; ii < npts ; ii++ ) sum += (ts[ii]-ts_mean)
00430                                                 *(ts[ii]-ts_mean) ;
00431          }
00432 
00433          sum = sqrt( sum/(npts-1) ) ;
00434 
00435          if((meth[meth_index] == METH_SIGMA) || (meth[meth_index] == METH_SIGMA_NOD))  val[out_index] = sum ;
00436          else if( ts_mean != 0.0 ) val[out_index] = sum / fabs(ts_mean) ;
00437          else                      val[out_index] = 0.0 ;
00438       }
00439       break ;
00440 
00441       /* 14 Feb 2000: these 2 new methods disturb the array ts[] */
00442       /* 18 Dec 2002: these 2 methods no longer disturb the array ts[] */
00443 
00444       case METH_MEDIAN:{
00445          float* ts_copy;
00446          ts_copy = (float*)calloc(npts, sizeof(float));
00447          memcpy( ts_copy, ts, npts * sizeof(float));
00448          val[out_index] = qmed_float( npts , ts_copy ) ;
00449          free(ts_copy);
00450       }
00451       break ;
00452 
00453       case METH_MAD:{
00454          float* ts_copy;
00455          register int ii ;
00456          register float vm ;
00457          ts_copy = (float*)calloc(npts, sizeof(float));
00458          memcpy( ts_copy, ts, npts * sizeof(float));
00459          vm = qmed_float( npts , ts_copy ) ;
00460          for( ii=0 ; ii < npts ; ii++ ) ts_copy[ii] = fabs(ts_copy[ii]-vm) ;
00461          val[out_index] = qmed_float( npts , ts_copy ) ;
00462          free(ts_copy);
00463       }
00464       break ;
00465 
00466       case METH_ARGMIN:
00467       case METH_MIN:{
00468          register int ii,outdex=0 ;
00469          register float vm=ts[0] ;
00470          for( ii=1 ; ii < npts ; ii++ ) if( ts[ii] < vm ) {
00471            vm = ts[ii] ;
00472            outdex = ii ;
00473          }
00474          if (meth[meth_index] == METH_MIN) {
00475            val[out_index] = vm ;
00476          } else {
00477            val[out_index] = outdex ;
00478          }
00479       }
00480       break ;
00481 
00482       case METH_DW:{
00483          register int ii ;
00484          register float den=ts[0]*ts[0] ;
00485          register float num=0 ;
00486          for( ii=1 ; ii < npts ; ii++ ) {
00487            num = num + (ts_det[ii] - ts_det[ii-1])
00488                        *(ts_det[ii] - ts_det[ii-1]);
00489            den = den + ts_det[ii] * ts_det[ii];
00490          }
00491          if (den == 0) {
00492            val[out_index] = 0 ;
00493          } else {
00494            val[out_index] = num/den ;
00495          }
00496       }
00497       break ;
00498 
00499       case METH_ABSMAX:
00500       case METH_ARGMAX:
00501       case METH_ARGABSMAX:
00502       case METH_MAX:{
00503          register int ii, outdex=0 ;
00504          register float vm=ts[0] ;
00505          if ((meth[meth_index] == METH_ABSMAX) || (meth[meth_index] == METH_ARGABSMAX)) {
00506            vm = fabs(vm) ;
00507            for( ii=1 ; ii < npts ; ii++ ) {
00508              if( fabs(ts[ii]) > vm ) {
00509                vm = fabs(ts[ii]) ;
00510                outdex = ii ;
00511              }
00512            }
00513          } else {
00514            for( ii=1 ; ii < npts ; ii++ ) {
00515              if( ts[ii] > vm ) {
00516                vm = ts[ii] ;
00517                outdex = ii ;
00518              }
00519            }
00520          }
00521 
00522          if ((meth[meth_index] == METH_ABSMAX) || (meth[meth_index] == METH_MAX)) {
00523            val[out_index] = vm ;
00524          } else {
00525            val[out_index] = outdex ;
00526          }
00527       }
00528       break ;
00529 
00530       case METH_AUTOCORR:{
00531         int numVals;
00532         float* ts_corr;
00533         /* for these new methods, the extra, needed integer */
00534         /* parameter is stored in the static array "meth",  */
00535         /* in the element right after the indicator for     */
00536         /* this method.  This parameter indicates the number*/
00537         /* of values to return, or 0 for the same length as */
00538         /* the input array.                                 */
00539         numVals = meth[meth_index + 1];
00540         if (numVals == 0) numVals = npts - 1;
00541         ts_corr = (float*)calloc(numVals,sizeof(float));
00542         /* Call the autocorrelation function, in this file. */
00543         autocorr(npts,ts,numVals,ts_corr);
00544         /* Copy the values into the output array val, which */
00545         /* will be returned to the fbuc MAKER caller to     */
00546         /* populate the appropriate BRIKs with the data.    */
00547         for( ii = 0; ii < numVals; ii++) {
00548           val[out_index + ii] = ts_corr[ii];
00549         }
00550         /* Although meth_index will be incremented by the   */
00551         /* for loop, it needs to be incremented an extra    */
00552         /* time to account for the extra parameter we       */
00553         /* pulled from the meth array.                      */
00554         meth_index++;
00555         /* Similarly, though the out_index will be incremented */
00556         /* by the for loop, we have added potentially several  */
00557         /* values to the output array, and we need to account  */
00558         /* for that here.                                      */
00559         out_index+=(numVals - 1);
00560         free(ts_corr);
00561       }
00562       break ;
00563 
00564       case METH_AUTOREGP:{
00565         int numVals,kk,mm;
00566         float *ts_corr, *y, *z;
00567         float alpha, beta, tmp;
00568 
00569         /* For these new methods, the extra, needed integer */
00570         /* parameter is stored in the static array "meth",  */
00571         /* in the element right after the indicator for     */
00572         /* this method.  This parameter indicates the number*/
00573         /* of values to return, or 0 for the same length as */
00574         /* the input array.                                 */
00575         numVals = meth[meth_index + 1];
00576         if (numVals == 0) numVals = npts - 1;
00577 
00578         /* Allocate space for the autocorrelation coefficients, */
00579         /* result, and a temp array for calculations.           */
00580         /* Correlation coeff array is larger, because we must   */
00581         /* disregard the r0, which is identically 1.            */
00582         /* Might fix this in autocorr function                  */
00583         ts_corr = (float*)malloc((numVals) * sizeof(float));
00584         y = (float*)malloc(numVals * sizeof(float));
00585         z = (float*)malloc(numVals * sizeof(float));
00586 
00587         /* Call autocorr function to generate the autocorrelation  */
00588         /* coefficients.                                           */
00589         autocorr(npts,ts,numVals,ts_corr);
00590 
00591         /* Durbin algorithm for solving Yule-Walker equations.        */
00592         /* The Yule-Walker equations specify the autoregressive       */
00593         /* coefficients in terms of the autocorrelation coefficients. */
00594         /* R*Phi = r, where r is vector of correlation coefficients,  */
00595         /* R is Toeplitz matrix formed from r, and Phi is a vector of */
00596         /* the autoregression coefficients.                           */
00597 
00598         /* In this implementation, 'y' is 'Phi' above and 'ts_corr' is 'r'    */
00599         
00600         y[0] = -ts_corr[0];
00601         alpha = -ts_corr[0];
00602         beta = 1;
00603         for (kk = 0 ; kk < (numVals - 1) ; kk++ ) {
00604           beta = (1 - alpha * alpha ) * beta ;
00605           tmp = 0;
00606           for ( mm = 0 ; mm <= kk ; mm++) {
00607             tmp = tmp + ts_corr[kk - mm] * y[mm];
00608           }
00609           alpha = - ( ts_corr[kk+1] + tmp ) / beta ;
00610           for ( mm = 0 ; mm <= kk ; mm++ ) {
00611             z[mm] = y[mm] + alpha * y[kk - mm];
00612           }
00613           for ( mm = 0 ; mm <= kk ; mm++ ) {
00614             y[mm] = z[mm];
00615           }
00616           y[kk+1] = alpha ;
00617         }
00618 
00619         /* Copy the values into the output array val, which */
00620         /* will be returned to the fbuc MAKER caller to     */
00621         /* populate the appropriate BRIKs with the data.    */
00622         for( ii = 0; ii < numVals; ii++) {
00623           val[out_index + ii] = y[ii];
00624           if (!finite(y[ii])) {
00625             fprintf(stderr,"BAD NUMBER y[%d] = %f, Call# %d\n",ii,y[ii],ncall);
00626           }
00627         }
00628         /* Although meth_index will be incremented by the   */
00629         /* for loop, it needs to be incremented an extra    */
00630         /* time to account for the extra parameter we       */
00631         /* pulled from the meth array.                      */
00632         meth_index++;
00633         /* Similarly, though the out_index will be incremented */
00634         /* by the for loop, we have added potentially several  */
00635         /* values to the output array, and we need to account  */
00636         /* for that here.                                      */
00637         out_index+=(numVals - 1);
00638         free(ts_corr);
00639         free(y);
00640         free(z);
00641       }
00642       break ;
00643    }
00644    }
00645 
00646    free(ts_det);
00647    ncall++ ; return ;
00648 }
00649 
00650 
00651 static void autocorr( int npts, float in_ts[], int numVals, float outcoeff[] )
00652 {
00653   /* autocorr is the inv fourier xform of the fourier xform abs squared */
00654 
00655   int ii,nfft;
00656   double scaler;
00657   complex * cxar = NULL;
00658 
00659   /* Calculate size for FFT, including padding for eliminating overlap  */
00660   /* from circular convolution */
00661   nfft = csfft_nextup_one35(npts * 2 - 1);
00662 /*  fprintf(stderr,"++ FFT length = %d\n",nfft) ; */
00663 
00664   cxar = (complex *) calloc( sizeof(complex) , nfft ) ;
00665 
00666   /* Populate complex array with input (real-only) time series */
00667   for( ii=0 ; ii < npts ; ii++ ){
00668     cxar[ii].r = in_ts[ii]; cxar[ii].i = 0.0;
00669   }
00670   /* Zero-pad input outside range of original time series */
00671   for( ii=npts ; ii < nfft ; ii++ ){ cxar[ii].r = cxar[ii].i = 0.0; }
00672 
00673   /* Calculate nfft-point FFT of input time series */
00674   /* First argument is "mode."  -1 value indicates */
00675   /* negative exponent in fourier sum, which means */
00676   /* this is an FFT and not an inverse FFT.        */
00677   /* Output will be complex and symmetrical.       */
00678   csfft_cox( -1 , nfft , cxar ) ;
00679 
00680   /* Use macro to calculate absolute square of FFT */
00681   for( ii=0 ; ii < nfft ; ii++ ){ 
00682     cxar[ii].r = CSQR(cxar[ii]) ; cxar[ii].i = 0 ; 
00683   }
00684 
00685   /* Take inverse FFT of result.  First function called */
00686   /* sets a static variable that the output should be   */
00687   /* scaled by 1/N.  This plus the mode argument of 1   */
00688   /* indicate that this is an inverse FFT.              */
00689   csfft_scale_inverse( 1 ) ;
00690   csfft_cox( 1 , nfft, cxar ) ;
00691 
00692   /* Copy the requested number of coefficients to the   */
00693   /* output array outcoeff.  These will be copied into  */
00694   /* the output BRIKs or used for further calculations  */
00695   /* in the calling function, STATS_tsfunc() above.     */
00696   /* The output coefficients are scaled by 1/(M-p) to   */
00697   /* provide an unbiased estimate, and also scaled by   */
00698   /* the final value of the zeroth coefficient.         */
00699   scaler = cxar[0].r/npts;
00700   for (ii = 0 ; ii < numVals ; ii++ ) {
00701     outcoeff[ii] = cxar[ii+1].r/((npts - (ii+1)) * scaler);
00702   }
00703   free(cxar);
00704   cxar = NULL;
00705 }