00165 {
00166    int ii , jj ;
00167    double val ;
00168 
00169    if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) Syntax() ;
00170 
00171    machdep(); if(argc > 1) AFNI_logger("waver",argc,argv) ;
00172 
00173    Process_Options( argc , argv ) ;
00174 
00175    /*---- compute duration of sampled waveform ----*/
00176 
00177    switch( waveform_type ){
00178 
00179       /* simple for the Cox function */
00180 
00181       default:
00182       case WAV_TYPE:
00183          WAV_duration = WAV_delay_time + WAV_rise_time + WAV_fall_time
00184                        + ( (WAV_undershoot != 0.0) ? WAV_restore_time : 0.0 ) ;
00185       break ;
00186 
00187       /* slightly complicated for the gamma variate */
00188 
00189       case GAM_TYPE:{
00190          double bal = 5.5/GAM_power + 1.0 ;
00191          double al  = bal ;
00192          double lal = log(al) ;
00193 
00194          while( al < bal + lal ){ al = bal + 1.1*lal ; lal = log(al) ; }
00195 
00196          WAV_duration = al * GAM_power * GAM_time ;
00197       }
00198       break ;
00199 
00200       /* 01 Aug 2001: yet more complicated for the EXPR type */
00201 
00202 #define FPASS 500
00203 #define SPASS 10000
00204 #define STHR  5
00205       case EXPR_TYPE:{
00206          double val , vtop=0.0 , vthr ;
00207          int itop=-1 , icount ;
00208          for( ii=0 ; ii < FPASS ; ii++ ){
00209             val = waveform_EXPR( WAV_dt * ii ) ; val = fabs(val) ;
00210             if( val > vtop ){ vtop = val ; itop = ii ; }
00211          }
00212          if( itop < 0 ){
00213             fprintf(stderr,"** -EXPR is 0 for 1st %d points!\n",FPASS);
00214             exit(1) ;
00215          }
00216          vthr = 0.01 * vtop ;
00217          for( icount=0,ii=itop+1 ; ii < SPASS && icount < STHR ; ii++ ){
00218             val = waveform_EXPR( WAV_dt * ii ) ; val = fabs(val) ;
00219             if( val <= vthr ) icount++ ;
00220             else              icount=0 ;
00221          }
00222          if( ii == SPASS && icount < STHR ){
00223             fprintf(stderr,"** -EXPR doesn't decay away in %d points!\n",SPASS);
00224             exit(1) ;
00225          }
00226          WAV_duration = WAV_dt * ii ;
00227 
00228          if( WAV_peak != 0.0 ) EXPR_fac = WAV_peak / vtop ;
00229          WAV_peak = 1.0 ;
00230       }
00231       break ;
00232    }
00233 
00234    /*---- compute sampled waveform ----*/
00235 
00236    WAV_npts = (int)(1 + ceil(WAV_duration/WAV_dt)) ;
00237    WAV_ts   = (double *) malloc( sizeof(double) * WAV_npts ) ;
00238 
00239    for( ii=0 ; ii < WAV_npts ; ii++ )
00240       WAV_ts[ii] = WAV_peak * waveform( WAV_dt * ii ) ;
00241 
00242    /*---- if no input timeseries, just output waveform ----*/
00243 
00244    if( IN_npts < 1 && IN_num_tstim < 1 ){
00245       int top = WAV_npts ;
00246       if( OUT_numout > 0 ) top = OUT_numout ;
00247       if( OUT_xy ){
00248          for( ii=0 ; ii < top ; ii++ )
00249             printf( "%g %g\n" , WAV_dt * ii , WAV_ts[ii] ) ;
00250       } else {
00251          for( ii=0 ; ii < top ; ii++ )
00252             printf( "%g\n" , WAV_ts[ii] ) ;
00253       }
00254       exit(0) ;
00255    }
00256 
00257    /*---- must convolve input with waveform ----*/
00258 
00259    if( IN_npts > 0 ){
00260       OUT_npts = IN_npts + WAV_npts ;
00261       if( OUT_numout > 0 ) OUT_npts = OUT_numout ;
00262       OUT_ts   = (double *) malloc( sizeof(double) * OUT_npts ) ;
00263       for( ii=0 ; ii < OUT_npts ; ii++ ) OUT_ts[ii] = 0.0 ;
00264 
00265       for( jj=0 ; jj < IN_npts ; jj++ ){
00266          val = IN_ts[jj] ;
00267          if( val == 0.0 || fabs(val) >= 33333.0 ) continue ;
00268          for( ii=0 ; ii < WAV_npts && ii+jj < OUT_npts ; ii++ )
00269             OUT_ts[ii+jj] += val * WAV_ts[ii] ;
00270       }
00271 
00272       for( jj=0 ; jj < IN_npts ; jj++ ){
00273          val = IN_ts[jj] ;
00274          if( fabs(val) >= 33333.0 ) OUT_ts[jj] = 99999.0 ;
00275       }
00276 
00277    } else if( IN_num_tstim > 0 ){  /* 16 May 2001 */
00278 #undef  TSTEP
00279 #define TSTEP 10              /* # expansion steps per WAV_dt */
00280       int ibot,itop , kk ;
00281       int nts ;
00282       double dts = WAV_dt/TSTEP , dur , aa ;
00283       double *tst , *ast ;
00284 
00285       /* setup the output array */
00286 
00287       OUT_npts = ceil(IN_top_tstim/WAV_dt) + WAV_npts ;
00288       if( OUT_numout > 0 ) OUT_npts = OUT_numout ;
00289       OUT_ts   = (double *) malloc( sizeof(double) * OUT_npts ) ;
00290       for( ii=0 ; ii < OUT_npts ; ii++ ) OUT_ts[ii] = 0.0 ;
00291 
00292       /* 12 May 2003: compute how many steps to expand to */
00293 
00294       nts = 0 ;
00295       for( jj=0 ; jj < IN_num_tstim ; jj++ ){
00296         dur = IN_tstim_b[jj] - IN_tstim_a[jj] ; dur = MAX(dur,0.0) ;
00297         ii  = ((int)ceil(dur/dts)) + 1 ;
00298         nts += ii ;
00299       }
00300 
00301       /* 12 May 2003: create expansion arrays */
00302 
00303       tst = (double *) malloc( sizeof(double) * nts ) ;
00304       ast = (double *) malloc( sizeof(double) * nts ) ;
00305       nts = 0 ;
00306       for( jj=0 ; jj < IN_num_tstim ; jj++ ){
00307         dur = IN_tstim_b[jj] - IN_tstim_a[jj] ; dur = MAX(dur,0.0) ;
00308         ii  = ((int)ceil(dur/dts)) + 1 ;
00309         if( ii == 1 ){              /* instantaneous impulse */
00310           tst[nts] = IN_tstim_a[jj] ; ast[nts] = IN_tstim_c[jj] ;
00311         } else {
00312           aa  = dur/(WAV_dt*ii) * IN_tstim_c[jj]; /* amplitude of each impulse */
00313           dur = dur / (ii-1) ;                    /* interval between impulses */
00314           for( kk=0 ; kk < ii ; kk++ ){
00315             tst[nts+kk] = IN_tstim_a[jj]+kk*dur ;
00316             ast[nts+kk] = aa ;
00317           }
00318         }
00319         nts += ii ;
00320       }
00321 
00322       /* Plop down copies of the waveform at each tst[] time */
00323 
00324       for( jj=0 ; jj < nts ; jj++ ){
00325         ibot = (int) (tst[jj]/WAV_dt) ;    /* may be 1 too early */
00326         itop = ibot + WAV_npts ;
00327         if( itop > OUT_npts ) itop = OUT_npts ;
00328         for( ii=ibot ; ii < itop ; ii++ ){
00329            val = WAV_peak * ast[jj] * waveform( WAV_dt * ii - tst[jj] ) ;
00330            OUT_ts[ii] += val ;
00331         }
00332       }
00333       free(ast); free(tst);
00334    }
00335 
00336    if( OUT_xy ){
00337       for( ii=0 ; ii < OUT_npts ; ii++ )
00338             printf( "%g %g\n" , WAV_dt * ii , OUT_ts[ii] ) ;
00339    } else {
00340       for( ii=0 ; ii < OUT_npts ; ii++ )
00341             printf( "%g\n" , OUT_ts[ii] ) ;
00342    }
00343 
00344    exit(0) ;
00345 }

00483 {
00484    int nopt = 1 ;
00485 
00486    while( nopt < argc ){
00487 
00488       if( strncmp(argv[nopt],"-GAM",4) == 0 ){
00489          waveform_type = GAM_TYPE ;
00490          nopt++ ; continue ;
00491       }
00492 
00493       if( strncmp(argv[nopt],"-WAV",4) == 0 ){
00494          waveform_type = WAV_TYPE ;
00495          nopt++ ; continue ;
00496       }
00497 
00498       if( strncmp(argv[nopt],"-EXPR",4) == 0 ){  /* 01 Aug 2001 */
00499          waveform_type = EXPR_TYPE ;
00500          if( EXPR_pcode != NULL ){
00501             fprintf(stderr,"** Can't have 2 -EXPR options!\n") ;
00502             exit(1) ;
00503          }
00504          nopt++ ;
00505          if( nopt >= argc ){
00506             fprintf(stderr,"** -EXPR needs an argument!\n") ; exit(1) ;
00507          }
00508          EXPR_pcode = PARSER_generate_code( argv[nopt] ) ;  /* compile */
00509          if( EXPR_pcode == NULL ){
00510             fprintf(stderr,"** Illegal -EXPR expression!\n") ;
00511             exit(1) ;
00512          }
00513          if( !PARSER_has_symbol("T",EXPR_pcode) ){
00514             fprintf(stderr,"** -EXPR expression doesn't use variable 't'!\n");
00515             exit(1) ;
00516          }
00517          nopt++ ; continue ;
00518       }
00519 
00520       if( strncmp(argv[nopt],"-gamb",5) == 0 ){
00521          if( nopt+1 >= argc ) ERROR ;
00522          GAM_power = strtod(argv[nopt+1],NULL) ;
00523          if( GAM_power <= 0.0 ) ERROR ;
00524          waveform_type = GAM_TYPE ;
00525          nopt++ ; nopt++ ; continue ;
00526       }
00527 
00528       if( strncmp(argv[nopt],"-gamc",5) == 0 ){
00529          if( nopt+1 >= argc ) ERROR ;
00530          GAM_time = strtod(argv[nopt+1],NULL) ;
00531          if( GAM_time <= 0.0 ) ERROR ;
00532          waveform_type = GAM_TYPE ;
00533          nopt++ ; nopt++ ; continue ;
00534       }
00535 
00536       if( strncmp(argv[nopt],"-gamd",5) == 0 ){
00537          if( nopt+1 >= argc ) ERROR ;
00538          GAM_delay_time = strtod(argv[nopt+1],NULL) ;
00539          /*if( GAM_time <= 0.0 ) ERROR ;*/
00540          waveform_type = GAM_TYPE ;
00541          nopt++ ; nopt++ ; continue ;
00542       }
00543 
00544       if( strncmp(argv[nopt],"-del",4) == 0 ){
00545          if( nopt+1 >= argc ) ERROR ;
00546          WAV_delay_time = strtod(argv[nopt+1],NULL) ;
00547          if( WAV_delay_time < 0.0 ) ERROR ;
00548          waveform_type = WAV_TYPE ;
00549          nopt++ ; nopt++ ; continue ;
00550       }
00551 
00552       if( strncmp(argv[nopt],"-ris",4) == 0 ){
00553          if( nopt+1 >= argc ) ERROR ;
00554          WAV_rise_time = strtod(argv[nopt+1],NULL) ;
00555          if( WAV_rise_time <= 0.0 ) ERROR ;
00556          waveform_type = WAV_TYPE ;
00557          nopt++ ; nopt++ ; continue ;
00558       }
00559 
00560       if( strncmp(argv[nopt],"-fal",4) == 0 ){
00561          if( nopt+1 >= argc ) ERROR ;
00562          WAV_fall_time = strtod(argv[nopt+1],NULL) ;
00563          if( WAV_fall_time <= 0.0 ) ERROR ;
00564          waveform_type = WAV_TYPE ;
00565          nopt++ ; nopt++ ; continue ;
00566       }
00567 
00568       if( strncmp(argv[nopt],"-und",4) == 0 ){
00569          if( nopt+1 >= argc ) ERROR ;
00570          WAV_undershoot = strtod(argv[nopt+1],NULL) ;
00571          waveform_type = WAV_TYPE ;
00572          nopt++ ; nopt++ ; continue ;
00573       }
00574 
00575       if( strncmp(argv[nopt],"-pea",4) == 0 ){
00576          if( nopt+1 >= argc ) ERROR ;
00577          WAV_peak = strtod(argv[nopt+1],NULL) ;
00578          nopt++ ; nopt++ ; continue ;
00579       }
00580 
00581       if( strncmp(argv[nopt],"-res",4) == 0 ){
00582          if( nopt+1 >= argc ) ERROR ;
00583          WAV_restore_time = strtod(argv[nopt+1],NULL) ;
00584          if( WAV_restore_time <= 0.0 ) ERROR ;
00585          waveform_type = WAV_TYPE ;
00586          nopt++ ; nopt++ ; continue ;
00587       }
00588 
00589       if( strncmp(argv[nopt],"-dt",3) == 0 || strncmp(argv[nopt],"-TR",3) == 0 ){
00590          if( nopt+1 >= argc ) ERROR ;
00591          WAV_dt = strtod(argv[nopt+1],NULL) ;
00592          if( WAV_dt <= 0.0 ) ERROR ;
00593          nopt++ ; nopt++ ; continue ;
00594       }
00595 
00596       if( strncmp(argv[nopt],"-xyo",4) == 0 ){
00597          OUT_xy = 1 ;
00598          nopt++ ; continue ;
00599       }
00600 
00601       if( strncmp(argv[nopt],"-inp",4) == 0 ){
00602          MRI_IMAGE * tsim ;
00603          float * tsar ;
00604          int ii ;
00605 
00606          if( IN_npts > 0 || IN_num_tstim > 0 ){
00607             fprintf(stderr,"Cannot input two timeseries!\n") ;
00608             exit(1) ;
00609          }
00610 
00611          if( nopt+1 >= argc ) ERROR ;
00612          tsim = mri_read_1D( argv[nopt+1] ) ;
00613          if( tsim == NULL ) ERROR ;
00614 
00615          IN_npts = tsim->nx ;
00616          IN_ts   = (double *) malloc( sizeof(double) * IN_npts ) ;
00617          tsar    = MRI_FLOAT_PTR(tsim) ;
00618          for( ii=0 ; ii < IN_npts ; ii++ ) IN_ts[ii] = tsar[ii] ;
00619          mri_free(tsim) ;
00620          nopt++ ; nopt++ ; continue ;
00621       }
00622 
00623       if( strcmp(argv[nopt],"-tstim") == 0 ){  /* 16 May 2001 */
00624         int iopt , nnn , zero_valc=0 ;
00625         float value , valb , valc ;
00626         char *cpt , *dpt ;
00627 
00628         if( IN_num_tstim > 0 || IN_npts > 0 ){
00629           fprintf(stderr,"Cannot input two timeseries!\n"); exit(1);
00630         }
00631         if( nopt+1 >= argc ) ERROR ;
00632 
00633         iopt         = nopt+1 ;
00634         IN_num_tstim = 0 ;
00635         IN_tstim_a   = (double *) malloc( sizeof(double) ) ;
00636         IN_tstim_b   = (double *) malloc( sizeof(double) ) ;
00637         IN_tstim_c   = (double *) malloc( sizeof(double) ) ;
00638 
00639         /* loop over argv until get a '-' or get to end of args */
00640 
00641         while( iopt < argc && argv[iopt][0] != '-' ){
00642 
00643           if( isspace(argv[iopt][0]) ){   /* skip if starts with blank */
00644             fprintf(stderr,               /* (usually from Microsoft!) */
00645                     "** Skipping -tstim value #%d that starts with whitespace!\n",
00646                     IN_num_tstim ) ;
00647             iopt++; continue;
00648           }
00649 
00650           /* formats:  start            ==> instantaneous impulse
00651                        start%duration   ==> extended duration
00652                        start:end        ==> extended duration
00653              can also have "/amplitude" afterwards to scale result */
00654 
00655           valb = 0.0 ; valc = 1.0 ;
00656           if( strchr(argv[iopt],'%') != NULL ){                       /* 12 May 2003 */
00657             nnn = sscanf( argv[iopt] , "%f%%%f" , &value , &valb ) ;
00658             if( nnn == 2 ) valb += value ;
00659           } else if( strchr(argv[iopt],':') != NULL ){
00660             nnn = sscanf( argv[iopt] , "%f:%f"  , &value , &valb ) ;
00661           } else {
00662             nnn = sscanf( argv[iopt] , "%f"     , &value ) ;
00663           }
00664           if( nnn < 1 || value < 0.0 ){
00665             fprintf(stderr,"** Weird value after -tstim: argv='%s'\n",argv[iopt]  ) ;
00666             fprintf(stderr,"**                  previous argv='%s'\n",argv[iopt-1]) ;
00667             fprintf(stderr,"** ==> Skipping this value!\n") ;
00668             iopt++; continue;
00669           }
00670           if( nnn == 1 || valb < value ) valb = value ;  /* 12 May 2003 */
00671 
00672           /* 13 May 2005: check for amplitude that follows a 'x' */
00673 
00674           cpt = strchr(argv[iopt],'x') ;
00675           if( cpt != NULL ){
00676             cpt++ ; valc = strtod( cpt , &dpt ) ;
00677             if( valc == 0.0 && dpt == cpt ) valc = 1.0 ;
00678             if( valc == 0.0 ) zero_valc++ ;
00679           }
00680 
00681           IN_tstim_a = (double *)realloc(IN_tstim_a,sizeof(double)*(IN_num_tstim+1));
00682           IN_tstim_b = (double *)realloc(IN_tstim_b,sizeof(double)*(IN_num_tstim+1));
00683           IN_tstim_c = (double *)realloc(IN_tstim_c,sizeof(double)*(IN_num_tstim+1));
00684           IN_tstim_a[IN_num_tstim] = value ;   /* start time */
00685           IN_tstim_b[IN_num_tstim] = valb  ;   /* end time */
00686           IN_tstim_c[IN_num_tstim] = valc  ;   /* amplitude */
00687           IN_num_tstim++ ;
00688           if( valb > IN_top_tstim ) IN_top_tstim = valb ;
00689           iopt++ ;
00690         }
00691         if( zero_valc == IN_num_tstim )
00692           fprintf(stderr,
00693                   "** WARNING: all '/' amplitudes in 'waver -tstim' are zero!\n") ;
00694 
00695         nopt = iopt ; continue ;  /* end of -tstim */
00696       }
00697 
00698       if( strncmp(argv[nopt],"-inl",4) == 0 ){
00699          int iopt , count , nnn ;
00700          float value ;
00701          char sep ;
00702 
00703          if( IN_npts > 0 || IN_num_tstim > 0 ){
00704             fprintf(stderr,"Cannot input two timeseries!\n") ;
00705             exit(1) ;
00706          }
00707 
00708          if( nopt+1 >= argc ) ERROR ;
00709          iopt    = nopt+1 ;
00710          IN_npts = 0 ;
00711          IN_ts   = (double *) malloc( sizeof(double) ) ;
00712          while( iopt < argc && argv[iopt][0] != '-' ){
00713 
00714             if( strstr(argv[iopt],"@") != NULL ||    /* if has one of the    */
00715                 strstr(argv[iopt],"x") != NULL ||    /* allowed separator    */
00716                 strstr(argv[iopt],"X") != NULL ||    /* characters, then     */
00717                 strstr(argv[iopt],"*") != NULL   ){  /* scan for count@value */
00718 
00719                nnn = sscanf( argv[iopt] , "%d%c%f" , &count , &sep , &value ) ;
00720                if( nnn != 3 || count < 1 ){
00721                   fprintf(stderr,"Illegal value after -inline: %s\n",argv[iopt]) ;
00722                   exit(1) ;
00723                }
00724 
00725             } else {                                 /* just scan for value */
00726                count = 1 ;
00727                nnn   = sscanf( argv[iopt] , "%f" , &value ) ;
00728                if( nnn != 1 ){
00729                   fprintf(stderr,"Illegal value after -inline: %s\n",argv[iopt]) ;
00730                   exit(1) ;
00731                }
00732             }
00733 
00734             IN_ts = (double *) realloc( IN_ts , sizeof(double) * (IN_npts+count) ) ;
00735             for( nnn=0 ; nnn < count ; nnn++ )
00736                IN_ts[nnn+IN_npts] = value ;
00737 
00738             IN_npts += count ; iopt++ ;
00739          }
00740          nopt = iopt ; continue ;
00741       }
00742 
00743       if( strcmp(argv[nopt],"-when") == 0 ){   /* 08 Apr 2002 */
00744          int iopt , bot,top , nn , nbt,*bt , count=0 , ii,kk ;
00745          float value ;
00746          char sep ;
00747 
00748          if( IN_npts > 0 || IN_num_tstim > 0 ){
00749             fprintf(stderr,"Cannot input two timeseries!\n") ;
00750             exit(1) ;
00751          }
00752 
00753          if( nopt+1 >= argc ) ERROR ;
00754          iopt = nopt+1 ;
00755          nbt  = 0 ;
00756          bt   = (int *) malloc(sizeof(int)) ;
00757          while( iopt < argc && argv[iopt][0] != '-' ){
00758 
00759             /* scan for value..value */
00760 
00761             bot = top = -1 ;
00762             nn = sscanf( argv[iopt],"%d..%d",&bot,&top) ;
00763             if( nn < 1 || bot < 0 ){
00764               fprintf(stderr,"Illegal value after -when: %s\n",argv[iopt]) ;
00765               exit(1) ;
00766             }
00767             if( nn == 1 ){
00768               top = bot ;
00769             } else if( top < bot ){
00770               fprintf(stderr,"Illegal value after -when: %s\n",argv[iopt]) ;
00771               exit(1) ;
00772             }
00773 
00774             /* save (bot,top) pairs in bt */
00775 
00776             bt = (int *) realloc( bt , sizeof(int)*(nbt+1)*2 ) ;
00777             bt[2*nbt  ] = bot ;
00778             bt[2*nbt+1] = top ; nbt++ ;
00779             if( count < top ) count = top ;
00780             iopt++ ;
00781          }
00782 
00783          if( nbt < 1 ){
00784             fprintf(stderr,"No ranges after -when?\n") ; exit(1) ;
00785          }
00786 
00787          IN_npts = count+1 ;
00788          IN_ts   = (double *) malloc( sizeof(double) * IN_npts ) ;
00789          for( ii=0 ; ii < IN_npts ; ii++ ) IN_ts[ii] = 0.0 ;
00790 
00791          for( kk=0 ; kk < nbt ; kk++ ){
00792             bot = bt[2*kk] ; top = bt[2*kk+1] ;
00793             for( ii=bot ; ii <= top ; ii++ ) IN_ts[ii] = 1.0 ;
00794          }
00795 
00796          free(bt) ; nopt = iopt ; continue ;
00797       }
00798 
00799       if( strcmp(argv[nopt],"-numout") == 0 ){   /* 08 Apr 2002 */
00800          int val = -1 ;
00801          if( nopt+1 >= argc ) ERROR ;
00802          sscanf(argv[nopt+1],"%d",&val) ;
00803          if( val <= 1 ){
00804            fprintf(stderr,"Illegal value after -numout: %s\n",argv[nopt]);
00805            exit(1);
00806          }
00807          OUT_numout = val ;
00808          nopt++ ; nopt++ ; continue ;
00809       }
00810 
00811       ERROR ;
00812    }
00813 
00814    if( WAV_peak == 0.0 && waveform_type != EXPR_TYPE ){
00815       fprintf(stderr,"** Illegal -peak 0 for non-EXPR waveform type!\n") ;
00816       exit(1) ;
00817    }
00818 
00819    return ;
00820 }

00348 {
00349    printf(
00350     "Usage: waver [options] > output_filename\n"
00351     "Creates an ideal waveform timeseries file.\n"
00352     "The output goes to stdout, and normally would be redirected to a file.\n"
00353     "\n"
00354     "Options: (# refers to a number; [xx] is the default value)\n"
00355     "  -WAV = Sets waveform to Cox special                    [default]\n"
00356     "           (cf. AFNI FAQ list for formulas)\n"
00357     "  -GAM = Sets waveform to form t^b * exp(-t/c)\n"
00358     "           (cf. Mark Cohen)\n"
00359     "\n"
00360     "  -EXPR \"expression\" = Sets waveform to the expression given,\n"
00361     "                         which should depend on the variable 't'.\n"
00362     "     e.g.: -EXPR \"step(t-2)*step(12-t)*(t-2)*(12-t)\"\n"
00363     "     N.B.: The peak value of the expression on the '-dt' grid will\n"
00364     "           be scaled to the value given by '-peak'; if this is not\n"
00365     "           desired, set '-peak 0', and the 'natural' peak value of\n"
00366     "           the expression will be used.\n"
00367     "\n"
00368     "These options set parameters for the -WAV waveform.\n"
00369     "  -delaytime #   = Sets delay time to # seconds                [2]\n"
00370     "  -risetime #    = Sets rise time to # seconds                 [4]\n"
00371     "  -falltime #    = Sets fall time to # seconds                 [6]\n"
00372     "  -undershoot #  = Sets undershoot to # times the peak         [0.2]\n"
00373     "                     (this should be a nonnegative factor)\n"
00374     "  -restoretime # = Sets time to restore from undershoot        [2]\n"
00375     "\n"
00376     "These options set parameters for the -GAM waveform:\n"
00377     "  -gamb #        = Sets the parameter 'b' to #                 [8.6]\n"
00378     "  -gamc #        = Sets the parameter 'c' to #                 [0.547]\n"
00379     "  -gamd #        = Sets the delay time to # seconds            [0.0]\n"
00380     "\n"
00381     "These options apply to all waveform types:\n"
00382     "  -peak #        = Sets peak value to #                        [100]\n"
00383     "  -dt #          = Sets time step of output AND input          [0.1]\n"
00384     "  -TR #          = '-TR' is equivalent to '-dt'\n"
00385     "\n"
00386     "The default is just to output the waveform defined by the parameters\n"
00387     "above.  If an input file is specified by one the options below, then\n"
00388     "the timeseries defined by that file will be convolved with the ideal\n"
00389     "waveform defined above -- that is, each nonzero point in the input\n"
00390     "timeseries will generate a copy of the waveform starting at that point\n"
00391     "in time, with the amplitude scaled by the input timeseries value.\n"
00392     "\n"
00393     "  -xyout         = Output data in 2 columns:\n"
00394     "                     1=time 2=waveform (useful for graphing)\n"
00395     "                     [default is 1 column=waveform]\n"
00396     "\n"
00397     "  -input infile  = Read timeseries from *.1D formatted 'infile';\n"
00398     "                     convolve with waveform to produce output\n"
00399     "              N.B.: you can use a sub-vector selector to choose\n"
00400     "                    a particular column of infile, as in\n"
00401     "                      -input 'fred.1D[3]'\n"
00402     "\n"
00403     "  -inline DATA   = Read timeseries from command line DATA;\n"
00404     "                     convolve with waveform to produce output\n"
00405     "                     DATA is in the form of numbers and\n"
00406     "                     count@value, as in\n"
00407     "                     -inline 20@0.0 5@1.0 30@0.0 1.0 20@0.0 2.0\n"
00408     "     which means a timeseries with 20 zeros, then 5 ones, then 30 zeros,\n"
00409     "     a single 1, 20 more zeros, and a final 2.\n"
00410     "     [The '@' character may actually be any of: '@', '*', 'x', 'X'.\n"
00411     "      Note that * must be typed as \\* to prevent the shell from\n"
00412     "      trying to interpret it as a filename wildcard.]\n"
00413     "\n"
00414     "  -tstim DATA    = Read discrete stimulation times from the command line\n"
00415     "                     and convolve the waveform with delta-functions at\n"
00416     "                     those times.  In this input format, the times do\n"
00417     "                     NOT have to be at intervals of '-dt'.  For example\n"
00418     "                       -dt 2.0 -tstim 5.6 9.3 13.7 16.4\n"
00419     "                     specifies a TR of 2 s and stimuli at 4 times\n"
00420     "                     (5.6 s, etc.) that do not correspond to integer\n"
00421     "                     multiples of TR.  DATA values cannot be negative.\n"
00422     "                   If the DATA is stored in a file, you can read it\n"
00423     "                     onto the command line using something like\n"
00424     "                       -tstim `cat filename`\n"
00425     "                     where using the backward-single-quote operator\n"
00426     "                     of the usual Unix shells.\n"
00427     "   ** 12 May 2003: The times after '-tstim' can now also be specified\n"
00428     "                     in the format 'a:b', indicating a continuous ON\n"
00429     "                     period from time 'a' to time 'b'.  For example,\n"
00430     "                       -dt 2.0 -tstim 13.2:15.7 20.3:25.3\n"
00431     "                     The amplitude of a response of duration equal to\n"
00432     "                     'dt' is equal the the amplitude of a single impulse\n"
00433     "                     response (which is the special case a=b).  N.B.: This\n"
00434     "                     means that something like '5:5.01' is very different\n"
00435     "                     from '5' (='5:5').  The former will have a small amplitude\n"
00436     "                     because of the small duration, but the latter will have\n"
00437     "                     a large amplitude because the case of an instantaneous\n"
00438     "                     input is special.  It is probably best NOT to mix the\n"
00439     "                     two types of input to '-tstim' for this reason.\n"
00440     "                     Compare the graphs from the 2 commands below:\n"
00441     "                       waver -dt 1.0 -tstim 5:5.1 | 1dplot -stdin\n"
00442     "                       waver -dt 1.0 -tstim 5     | 1dplot -stdin\n"
00443     "                     If you prefer, you can use the form 'a%%c' to indicate\n"
00444     "                     an ON interval from time=a to time=a+c.\n"
00445     "   ** 13 May 2005: You can now add an amplitude to each response individually.\n"
00446     "                     For example\n"
00447     "                       waver -dt 1.0 -peak 1.0 -tstim 3.2 17.9x2.0 23.1x-0.5\n"
00448     "                     puts the default response amplitude at time 3.2,\n"
00449     "                     2.0 times the default at time 17.9, and -0.5 times\n"
00450     "                     the default at time 23.1.\n"
00451     "\n"
00452     "  -when DATA     = Read time blocks when stimulus is 'on' (=1) from the\n"
00453     "                     command line and convolve the waveform with with\n"
00454     "                     a zero-one input.  For example:\n"
00455     "                       -when 20..40 60..80\n"
00456     "                     means that the stimulus function is 1.0 for time\n"
00457     "                     steps number 20 to 40, and 60 to 80 (inclusive),\n"
00458     "                     and zero otherwise.  (The first time step is\n"
00459     "                     numbered 0.)\n"
00460     "\n"
00461     "  -numout NN     = Output a timeseries with NN points; if this option\n"
00462     "                     is not given, then enough points are output to\n"
00463     "                     let the result tail back down to zero.\n"
00464     "\n"
00465     "At least one option is required, or the program will just print this message\n"
00466     "to stdout.  Only one of the 3 timeseries input options above can be used.\n"
00467     "\n"
00468     "If you have the 'xmgr' graphing program, then a useful way to preview the\n"
00469     "results of this program is through a command pipe like\n"
00470     "   waver -dt 0.25 -xyout -inline 16@1 40@0 16@1 40@0 | xmgr -source stdin\n"
00471     "Using the cruder AFNI package program 1dplot, you can do something like:\n"
00472     "   waver -GAM -tstim 0 7.7 | 1dplot -stdin\n"
00473     "\n"
00474     "If a square wave is desired, see the 'sqwave' program.\n"
00475    ) ;
00476    exit(0) ;
00477 }

00076 {
00077    switch( waveform_type ){
00078 
00079       default:
00080       case WAV_TYPE:  return waveform_WAV(t) ;
00081 
00082       case GAM_TYPE:  return waveform_GAM(t) ;
00083 
00084       case EXPR_TYPE: return waveform_EXPR(t);  /* 01 Aug 2001 */
00085    }
00086    return 0.0 ;  /* unreachable */
00087 }

00091 {
00092    static int first=1 ;
00093    static double atoz[26] ;
00094 
00095    if( t < 0.0 ) return 0.0 ;     /* 02 Aug 2001: oops */
00096 
00097    if( first ){
00098       int ii ;
00099       for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ;
00100       first =0 ;
00101    }
00102    atoz[TT] = t ;
00103    return (EXPR_fac * PARSER_evaluate_one(EXPR_pcode,atoz) ) ;
00104 }

00107 {
00108    if( GAM_ampl == 0.0 ){
00109       GAM_ampl = exp(GAM_power) / pow(GAM_power*GAM_time,GAM_power) ;
00110    }
00111 
00112    if( t-GAM_delay_time <= 0.0 ) return 0.0 ;
00113 
00114    return GAM_ampl * pow((t-GAM_delay_time),GAM_power) * exp(-(t-GAM_delay_time)/GAM_time) ;
00115 }

00118 {
00119    if( WAV_rise_start < 0.0 ){
00120       WAV_rise_start   = WAV_delay_time ;
00121       WAV_fall_start   = WAV_rise_start + WAV_rise_time ;
00122       WAV_fall_end     = WAV_fall_start + WAV_fall_time ;
00123       WAV_restore_end  = WAV_fall_end   + WAV_restore_time ;
00124    }
00125 
00126    if( t < WAV_rise_start )
00127       return 0.0 ;
00128 
00129    if( t < WAV_fall_start )
00130       return ztone( (t-WAV_rise_start)/WAV_rise_time ) ;
00131 
00132    if( t < WAV_fall_end )
00133       return (1.0+WAV_undershoot) * ztone( (WAV_fall_end-t)/WAV_fall_time )
00134              - WAV_undershoot ;
00135 
00136    if( t < WAV_restore_end )
00137       return -WAV_undershoot * ztone( (WAV_restore_end-t)/WAV_restore_time ) ;
00138 
00139    return 0.0 ;
00140 }

00035 {
00036    double y , z ;
00037 
00038    if( x <= 0.0 ) return 0.0 ;
00039    if( x >= 1.0 ) return 1.0 ;
00040 
00041    y = (0.5*PI) * ( 1.6 * x - 0.8 ) ;
00042    z = ZT_FAC * ( tanh(tan(y)) + ZT_ADD ) ;
00043    return z ;
00044 }