do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ;  \
      FREEUP(fout) ; FREEUP(cxar) ; FREEUP(tar)  ;  \
      FREEUP(fxar) ; FREEUP(fyar) ; FREEUP(dtr)  ;  \
    } while(0)

00096 {
00097    PLUGIN_interface * plint ;     /* will be the output of this routine */
00098 
00099    if( ncall > 1 ) return NULL ;  /* two interfaces */
00100 
00101 #ifdef ALLOW_TESTING
00102    if( ncall == 1 ) return TEST_init() ;
00103 #else
00104    if( ncall == 1 ) return NULL ;
00105 #endif
00106 
00107    /*---------------- set titles and call point ----------------*/
00108 
00109    plint = PLUTO_new_interface( "Power Spectrum" ,
00110                                 "Power Spectrum of a 3D+time Dataset" ,
00111                                 helpstring ,
00112                                 PLUGIN_CALL_VIA_MENU , POWER_main  ) ;
00113 
00114    PLUTO_add_hint( plint , "Power Spectrum of a 3D+time Dataset" ) ;
00115 
00116    PLUTO_set_sequence( plint , "A:newdset:statistics" ) ;
00117 
00118    /*--------- 1st line: Input dataset ---------*/
00119 
00120    PLUTO_add_option( plint ,
00121                      "Input" ,  /* label at left of input line */
00122                      "Input" ,  /* tag to return to plugin */
00123                      TRUE       /* is this mandatory? */
00124                    ) ;
00125 
00126    PLUTO_add_dataset(  plint ,
00127                        "---->>" ,         /* label next to button   */
00128                        ANAT_ALL_MASK ,    /* take any anat datasets */
00129                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
00130                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00131                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00132                     ) ;
00133 
00134    /*---------- 2nd line: Output dataset ----------*/
00135 
00136    PLUTO_add_option( plint ,
00137                      "Output" ,  /* label at left of input line */
00138                      "Output" ,  /* tag to return to plugin */
00139                      TRUE        /* is this mandatory? */
00140                    ) ;
00141 
00142    PLUTO_add_string(   plint ,
00143                        "Prefix" ,  /* label next to textfield */
00144                        0,NULL ,    /* no fixed strings to choose among */
00145                        19          /* 19 spaces for typing in value */
00146                    ) ;
00147 
00148    PLUTO_add_string(   plint ,
00149                        "Datum" ,          /* label next to chooser button */
00150                        NUM_TYPE_STRINGS , /* number of strings to choose among */
00151                        type_strings ,     /* list of strings to choose among */
00152                        0                  /* index of default string */
00153                    ) ;
00154 
00155    /*--------- Other lines: Parameters ---------*/
00156 
00157    PLUTO_add_option( plint , "Ignore" , "Ignore" , TRUE ) ;
00158 
00159    PLUTO_add_number( plint ,
00160                      "Count" ,   /* label next to chooser */
00161                      0 ,         /* smallest possible value */
00162                      999 ,       /* largest possible value */
00163                      0 ,         /* decimal shift (none in this case) */
00164                      4 ,         /* default value */
00165                      TRUE        /* allow user to edit value? */
00166                    ) ;
00167 
00168    PLUTO_add_option( plint , "Taper" , "Taper" , TRUE ) ;
00169 
00170    PLUTO_add_number( plint ,
00171                      "Percent" ,    /* label next to chooser */
00172                      0 ,            /* smallest possible value */
00173                      10 ,           /* largest possible value */
00174                      -1 ,           /* decimal shift (1 right == 0 to 100) */
00175                      0 ,            /* default value (with shift == 0) */
00176                      FALSE          /* allow user to edit value? */
00177                    ) ;
00178 
00179    PLUTO_add_option( plint , "FFT" , "FFT" , TRUE ) ;
00180 
00181    PLUTO_add_string( plint ,
00182                      "Length" ,         /* label next to chooser */
00183                      NUM_FFT_STRINGS ,  /* number of strings to choose among */
00184                      fft_strings ,      /* list of strings to choose among */
00185                      0                  /* index of default string */
00186                    ) ;
00187 
00188    /*--------- done with interface setup ---------*/
00189 
00190    return plint ;
00191 }

00213 {
00214    MCW_idcode * idc ;                          /* input dataset idcode */
00215    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00216    char * new_prefix , * str ;                 /* strings from user */
00217    int   new_datum , ignore , nfft , ninp ,    /* control parameters */
00218          old_datum , nuse , ntaper , ktbot ;
00219    float taper ;
00220 
00221    byte   ** bptr  = NULL ;  /* one of these will be the array of */
00222    short  ** sptr  = NULL ;  /* pointers to input dataset sub-bricks */
00223    float  ** fptr  = NULL ;  /* (depending on input datum type) */
00224 
00225    complex * cxar  = NULL ;  /* will be array of data to FFT */
00226    float   * fxar  = NULL ;  /* array loaded from input dataset */
00227    float   * fyar  = NULL ;  /* array loaded from input dataset */
00228    float  ** fout  = NULL ;  /* will be array of output floats */
00229 
00230    float   * tar   = NULL ;  /* will be array of taper coefficients */
00231    float   * dtr   = NULL ;  /* will be array of detrending coeff */
00232 
00233    float dfreq , pfact , phi , xr,xi , yr,yi ;
00234    float x0,x1 , y0,y1 , d0fac,d1fac ;
00235    int   nfreq , nvox , perc , new_units ;
00236    int   istr , ii,iip , ibot,itop , kk , icx ;       /* temp variables */
00237 
00238    /*--------------------------------------------------------------------*/
00239    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00240 
00241    /*--------- go to first input line ---------*/
00242 
00243    PLUTO_next_option(plint) ;
00244 
00245    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
00246    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00247    if( old_dset == NULL )
00248       return "*************************\n"
00249              "Cannot find Input Dataset\n"
00250              "*************************"  ;
00251 
00252    /*--------- go to second input line ---------*/
00253 
00254    PLUTO_next_option(plint) ;
00255 
00256    new_prefix = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
00257    if( ! PLUTO_prefix_ok(new_prefix) )      /* check if it is OK */
00258       return "************************\n"
00259              "Output Prefix is illegal\n"
00260              "************************"  ;
00261 
00262    str  = PLUTO_get_string(plint) ;              /* get string item (the datum type) */
00263    istr = PLUTO_string_index( str ,              /* find it in the list it came from */
00264                               NUM_TYPE_STRINGS ,
00265                               type_strings ) ;
00266    switch( istr ){
00267       default:
00268       case 0:
00269          new_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;  /* use old dataset type */
00270       break ;
00271 
00272       case 1: new_datum = MRI_byte  ; break ;  /* assign type of user's choice */
00273       case 2: new_datum = MRI_short ; break ;
00274       case 3: new_datum = MRI_float ; break ;
00275    }
00276 
00277    /*--------- go to next input lines ---------*/
00278 
00279    PLUTO_next_option(plint) ;                 /* skip to next line */
00280    ignore = PLUTO_get_number(plint) ;         /* get number item (ignore) */
00281 
00282    PLUTO_next_option(plint) ;                 /* skip to next line */
00283    taper  = PLUTO_get_number(plint) * 0.01 ;  /* get number item (taper %) */
00284 
00285    /* compute FFT length to use */
00286 
00287    PLUTO_next_option(plint) ;          /* skip to next line */
00288 
00289    str  = PLUTO_get_string(plint) ;    /* get string item for FFT count */
00290    ninp = DSET_NUM_TIMES(old_dset) ;   /* number of values in input */
00291    nuse = ninp - ignore ;              /* number of values to actually use */
00292 
00293    if( nuse < 4 )
00294       return "*****************************\n"
00295              "Not enough time points to FFT\n"
00296              "*****************************"  ;
00297 
00298    if( strcmp(str,fft_strings[0]) == 0 ){
00299 
00300       /*-- get next larger power-of-2 --*/
00301 #if 0
00302       for( nfft=32 ; nfft < nuse ; nfft *= 2 ) ; /* loop until nfft >= nuse */
00303 #else
00304       nfft = csfft_nextup_one35(nuse) ;
00305 #endif
00306 
00307    } else {
00308       nfft = strtol( str , NULL , 10 ) ;  /* just convert string to integer */
00309    }
00310 
00311    /* if the input FFT length is less than the data length,
00312       tell the user and truncate the amount of data to use */
00313 
00314    if( nfft < nuse ){
00315       char str[256] ;
00316 
00317       sprintf( str , "******************************\n"
00318                      "Warning:\n"
00319                      " Number of points in FFT =%4d\n"
00320                      " is less than available data\n"
00321                      " in time series = %d\n"
00322                      "******************************" ,
00323                nfft , nuse ) ;
00324 
00325       PLUTO_popup_transient( plint , str ) ;
00326 
00327       nuse = nfft ;  /* can't use more data than the FFT length */
00328    }
00329 
00330    /* compute the number of output points and the output grid spacing */
00331 
00332    nfreq = nfft / 2 ;                                 /* # frequencies */
00333    dfreq = 1.0 / (nfft * DSET_TIMESTEP(old_dset) ) ;  /* frequency grid */
00334 
00335    switch( DSET_TIMEUNITS(old_dset) ){
00336       case UNITS_MSEC_TYPE: dfreq *= 1000.0 ; new_units = UNITS_HZ_TYPE ; break;
00337       case UNITS_SEC_TYPE:                    new_units = UNITS_HZ_TYPE ; break;
00338       case UNITS_HZ_TYPE:                     new_units = UNITS_SEC_TYPE; break;
00339 
00340       default: new_units = DSET_TIMEUNITS(old_dset) ; break ; /* shouldn't happen */
00341    }
00342 
00343    /*------------------------------------------------------*/
00344    /*---------- At this point, the inputs are OK ----------*/
00345 
00346    PLUTO_popup_meter( plint ) ;  /* popup a progress meter */
00347 
00348    /*--------- set up pointers to each sub-brick in the input dataset ---------*/
00349 
00350    DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */
00351 
00352    old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
00353 
00354    switch( old_datum ){  /* pointer type depends on input datum type */
00355 
00356       default:
00357          return "******************************\n"
00358                 "Illegal datum in Input Dataset\n"
00359                 "******************************"  ;
00360 
00361       /** create array of pointers into old dataset sub-bricks **/
00362       /** Note that we skip the first 'ignore' sub-bricks here **/
00363 
00364       /*--------- input is bytes ----------*/
00365       /* voxel #i at time #k is bptr[k][i] */
00366       /* for i=0..nvox-1 and k=0..nuse-1.  */
00367 
00368       case MRI_byte:
00369          bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
00370          if( bptr == NULL ) return "Malloc\nFailure!\n [bptr]" ;
00371          for( kk=0 ; kk < nuse ; kk++ )
00372             bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
00373       break ;
00374 
00375       /*--------- input is shorts ---------*/
00376       /* voxel #i at time #k is sptr[k][i] */
00377       /* for i=0..nvox-1 and k=0..nuse-1.  */
00378 
00379       case MRI_short:
00380          sptr = (short **) malloc( sizeof(short *) * nuse ) ;
00381          if( sptr == NULL ) return "Malloc\nFailure!\n [sptr]" ;
00382          for( kk=0 ; kk < nuse ; kk++ )
00383             sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
00384       break ;
00385 
00386       /*--------- input is floats ---------*/
00387       /* voxel #i at time #k is fptr[k][i] */
00388       /* for i=0..nvox-1 and k=0..nuse-1.  */
00389 
00390       case MRI_float:
00391          fptr = (float **) malloc( sizeof(float *) * nuse ) ;
00392          if( fptr == NULL ) return "Malloc\nFailure!\n [fptr]" ;
00393          for( kk=0 ; kk < nuse ; kk++ )
00394             fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
00395       break ;
00396 
00397    } /* end of switch on input type */
00398 
00399    /*---- allocate space for 2 voxel timeseries and 1 FFT ----*/
00400 
00401    cxar = (complex *) malloc( sizeof(complex) * nfft ) ; /* FFT */
00402    fxar = (float *)   malloc( sizeof(float) * nuse ) ;   /* input */
00403    fyar = (float *)   malloc( sizeof(float) * nuse ) ;   /* input */
00404    if( cxar == NULL || fxar == NULL || fyar == NULL ){
00405       FREE_WORKSPACE ;
00406       return "Malloc\nFailure!\n [cxar]" ;
00407    }
00408 
00409    /*--------- make space for taper coefficient array ---------*/
00410 
00411    tar = (float *) malloc( sizeof(float) * MAX(nuse,nfreq) ) ;
00412    dtr = (float *) malloc( sizeof(float) * nuse ) ;
00413 
00414    if( tar == NULL || dtr == NULL ){
00415       FREE_WORKSPACE ;
00416       return "Malloc\nFailure!\n [tar]" ;
00417    }
00418 
00419    ntaper = (int)(0.5 * taper * nuse + 0.49) ; /* will taper data over */
00420    phi    = PI / MAX(ntaper,1) ;               /* kk=0..ntaper-1 on left */
00421    ktbot  = nuse - ntaper ;                    /* kk=ktbot..nuse-1 on right */
00422    pfact  = 0.0 ;                              /* sum of taper**2 */
00423 
00424    for( kk=0 ; kk < nuse ; kk++ ){                       /* Hamming-ize */
00425       if( kk < ntaper )
00426          tar[kk] = 0.54 - 0.46 * cos(kk*phi) ;           /* ramp up */
00427       else if( kk >= ktbot )
00428          tar[kk] = 0.54 + 0.46 * cos((kk-ktbot+1)*phi) ; /* ramp down */
00429       else
00430          tar[kk] = 1.0 ;                                 /* in the middle */
00431 
00432       pfact  += tar[kk] * tar[kk] ;
00433 
00434       dtr[kk] = kk - 0.5 * (nuse-1) ;  /* factors for linear detrending */
00435    }
00436 
00437    d0fac = 1.0 / nuse ;
00438    d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
00439 
00440    /*--- compute factor to go from |FFT|**2 to PSD;
00441          includes the scaling needed for loss of energy with tapering ---*/
00442 
00443    pfact = DSET_TIMESTEP(old_dset) / pfact ;
00444 
00445    /*--- include scaling factors for sub-bricks, if any ---*/
00446 
00447    for( kk=0 ; kk < nuse ; kk++ )
00448       if( DSET_BRICK_FACTOR(old_dset,kk+ignore) > 0.0 )
00449          tar[kk] *= DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
00450 
00451    /*---------------------- make a new dataset ----------------------*/
00452 
00453    new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
00454 
00455    { char * his = PLUTO_commandstring(plint) ;
00456      tross_Copy_History( old_dset , new_dset ) ;
00457      tross_Append_History( new_dset , his ) ; free(his) ;
00458    }
00459 
00460    /*-- edit some of its internal parameters --*/
00461 
00462    ii = EDIT_dset_items(
00463            new_dset ,
00464               ADN_prefix      , new_prefix ,           /* filename prefix */
00465               ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
00466               ADN_datum_all   , new_datum ,            /* atomic datum */
00467               ADN_nvals       , nfreq ,                /* # sub-bricks */
00468               ADN_ntt         , nfreq ,                /* # time points */
00469               ADN_ttorg       , dfreq ,                /* time origin */
00470               ADN_ttdel       , dfreq ,                /* time step */
00471               ADN_ttdur       , dfreq ,                /* time duration */
00472               ADN_nsl         , 0 ,                    /* z-axis time slicing */
00473               ADN_tunits      , new_units ,            /* time units */
00474            ADN_none ) ;
00475 
00476    if( ii != 0 ){
00477       THD_delete_3dim_dataset( new_dset , False ) ;
00478       FREE_WORKSPACE ;
00479       return "***********************************\n"
00480              "Error while creating output dataset\n"
00481              "***********************************"  ;
00482    }
00483 
00484    /*------ make floating point output sub-bricks
00485             (only at the end will scale to byte or shorts)
00486 
00487             Output #ii at freq #kk will go into fout[kk][ii],
00488             for kk=0..nfreq-1, and for ii=0..nvox-1.          ------*/
00489 
00490    nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
00491 
00492    fout = (float **) malloc( sizeof(float *) * nfreq ) ;  /* ptrs to sub-bricks */
00493 
00494    if( fout == NULL ){
00495       THD_delete_3dim_dataset( new_dset , False ) ;
00496       FREE_WORKSPACE ;
00497       return "Malloc\nFailure!\n [fout]" ;
00498    }
00499 
00500    for( kk=0 ; kk < nfreq ; kk++ ){
00501       fout[kk] = (float *) malloc( sizeof(float) * nvox ) ; /* sub-brick # kk */
00502       if( fout[kk] == NULL ) break ;
00503    }
00504 
00505    if( kk < nfreq ){
00506       for( ; kk >= 0 ; kk-- ) FREEUP(fout[kk]) ;   /* free all we did get */
00507       THD_delete_3dim_dataset( new_dset , False ) ;
00508       FREE_WORKSPACE ;
00509       return "Malloc\nFailure!\n [arrays]" ;
00510    }
00511 
00512    { char buf[128] ;
00513      ii = (nfreq * nvox * sizeof(float)) / (1024*1024) ;
00514      sprintf( buf , "  \n"
00515                     "*** 3D+time Power Spectral Density:\n"
00516                     "*** Using %d MBytes of workspace,\n "
00517                     "*** with FFT length = %d\n" , ii,nfft ) ;
00518      PLUTO_popup_transient( plint , buf ) ;
00519    }
00520 
00521    /*----------------------------------------------------*/
00522    /*----- Setup has ended.  Now do some real work. -----*/
00523 
00524    /***** loop over voxels *****/
00525 
00526    for( ii=0 ; ii < nvox ; ii += 2 ){  /* 2 time series at a time */
00527 
00528       iip = (ii+1) % nvox ;           /* voxel after ii */
00529 
00530       /*** load data from input dataset, depending on type ***/
00531 
00532       switch( old_datum ){
00533 
00534          /*** input = bytes ***/
00535 
00536          case MRI_byte:
00537             for( kk=0 ; kk < nuse ; kk++ ){
00538                fxar[kk] = bptr[kk][ii] ;
00539                fyar[kk] = bptr[kk][iip] ;
00540             }
00541          break ;
00542 
00543          /*** input = shorts ***/
00544 
00545          case MRI_short:
00546             for( kk=0 ; kk < nuse ; kk++ ){
00547                fxar[kk] = sptr[kk][ii] ;
00548                fyar[kk] = sptr[kk][iip] ;
00549             }
00550          break ;
00551 
00552          /*** input = floats ***/
00553 
00554          case MRI_float:
00555             for( kk=0 ; kk < nuse ; kk++ ){
00556                fxar[kk] = fptr[kk][ii] ;
00557                fyar[kk] = fptr[kk][iip] ;
00558             }
00559          break ;
00560 
00561       } /* end of switch over input type */
00562 
00563       /*** detrend:
00564              x0 = sum( fxar[kk] )
00565              x1 = sum( fxar[kk] * (kk-0.5*(N-1)) )
00566            x0 is used to remove the mean of fxar
00567            x1 is used to remove the linear trend of fxar ***/
00568 
00569       x0 = x1 = y0 = y1 = 0.0 ;
00570       for( kk=0 ; kk < nuse ; kk++ ){
00571          x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
00572          y0 += fyar[kk] ; y1 += fyar[kk] * dtr[kk] ;
00573       }
00574 
00575       x0 *= d0fac ; x1 *= d1fac ;  /* factors to remove mean and trend */
00576       y0 *= d0fac ; y1 *= d1fac ;
00577 
00578       for( kk=0 ; kk < nuse ; kk++ ){
00579          fxar[kk] -= (x0 + x1 * dtr[kk]) ;  /* remove mean and trend here! */
00580          fyar[kk] -= (y0 + y1 * dtr[kk]) ;
00581       }
00582 
00583       /*** taper, scale, and put into cxar array ***/
00584 
00585       for( kk=0 ; kk < nuse ; kk++ ){
00586          cxar[kk].r = fxar[kk] * tar[kk] ;
00587          cxar[kk].i = fyar[kk] * tar[kk] ;
00588       }
00589 
00590       /*** load zeros after where data was put ***/
00591 
00592       for( kk=nuse ; kk < nfft ; kk++ ) cxar[kk].r = cxar[kk].i = 0.0 ;
00593 
00594       /***** do the FFT (at long last) *****/
00595 
00596       csfft_cox( -1 , nfft , cxar ) ;
00597 
00598       /***** now compute output into corresponding voxels in fout *****/
00599 
00600       /*--- Let x = fxar (1st real time series)
00601                 y = fyar (2nd real time series)
00602                 z = cxar (complex time series) = x + i y
00603                 N = nfft (length of FFT)
00604 
00605             Then after FFT, since x and y are real, we have
00606               zhat[k]  = xhat[k] + i yhat[k]  > for k=1..N/2
00607             zhat[N-k]* = xhat[k] - i yhat[k]
00608 
00609             so we can untangle the FFTs of x and y by
00610               xhat[k] = 0.5 ( zhat[k] + zhat[N-k]* )
00611               yhat[k] = 0.5 ( zhat[k] - zhat[N-k]* ) / i
00612 
00613             This is the basis for doing 2 time series at once. ---*/
00614 
00615       for( kk=1 ; kk <= nfreq ; kk++ ){
00616          xr = 0.5 * ( cxar[kk].r + cxar[nfft-kk].r ) ; /* Re xhat[kk] */
00617          xi = 0.5 * ( cxar[kk].i - cxar[nfft-kk].i ) ; /* Im xhat[kk] */
00618          yr = 0.5 * ( cxar[kk].i + cxar[nfft-kk].i ) ; /* Re yhat[kk] */
00619          yi = 0.5 * ( cxar[kk].r - cxar[nfft-kk].r ) ; /*-Im yhat[kk] */
00620 
00621          fout[kk-1][ii]  = pfact * (xr*xr + xi*xi) ;
00622          fout[kk-1][iip] = pfact * (yr*yr + yi*yi) ;
00623       }
00624 
00625       perc = (100 * ii) / nvox ;        /* display percentage done */
00626       PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
00627 
00628    } /* end of outer loop over 2 voxels at a time */
00629 
00630    DSET_unload( old_dset ) ;  /* don't need this no more */
00631 
00632    /*------------------------------------------------------------*/
00633    /*------- The output is now in fout[kk][ii],
00634              for kk=0..nfreq-1 , ii=0..nvox-1.
00635              We must now put this into the output dataset -------*/
00636 
00637    switch( new_datum ){
00638 
00639       /*** output is floats is the simplest:
00640            we just have to attach the fout bricks to the dataset ***/
00641 
00642       case MRI_float:
00643          for( kk=0 ; kk < nfreq ; kk++ )
00644             EDIT_substitute_brick( new_dset , kk , MRI_float , fout[kk] ) ;
00645       break ;
00646 
00647       /*** output is shorts:
00648            we have to create a scaled sub-brick from fout ***/
00649 
00650       case MRI_short:{
00651          short * bout ;
00652          float fac ;
00653 
00654          for( kk=0 ; kk < nfreq ; kk++ ){  /* loop over sub-bricks */
00655 
00656             /*-- get output sub-brick --*/
00657 
00658             bout = (short *) malloc( sizeof(short) * nvox ) ;
00659             if( bout == NULL ){
00660                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
00661                return("\nFinal malloc error in plug_power!\n") ;
00662                /* EXIT(1) ;*/
00663             }
00664 
00665             /*-- find scaling and then scale --*/
00666 
00667             fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
00668             if( fac > 0.0 ){
00669                fac = 32767.0 / fac ;
00670                EDIT_coerce_scale_type( nvox,fac ,
00671                                        MRI_float,fout[kk] , MRI_short,bout ) ;
00672                fac = 1.0 / fac ;
00673             }
00674 
00675             free( fout[kk] ) ;  /* don't need this anymore */
00676 
00677             /*-- put output brick into dataset, and store scale factor --*/
00678 
00679             EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
00680             tar[kk] = fac ;
00681 
00682             perc = (100 * kk) / nfreq ;
00683             PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
00684          }
00685 
00686          /*-- save scale factor array into dataset --*/
00687 
00688          EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
00689 
00690       }
00691       break ;
00692 
00693       /*** output is bytes (byte = unsigned char)
00694            we have to create a scaled sub-brick from fout ***/
00695 
00696       case MRI_byte:{
00697          byte * bout ;
00698          float fac ;
00699 
00700          for( kk=0 ; kk < nfreq ; kk++ ){  /* loop over sub-bricks */
00701 
00702             /*-- get output sub-brick --*/
00703 
00704             bout = (byte *) malloc( sizeof(byte) * nvox ) ;
00705             if( bout == NULL ){
00706                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
00707                return("\nFinal malloc error in plug_power!\n\a") ;
00708                /* EXIT(1) ;*/
00709             }
00710 
00711             /*-- find scaling and then scale --*/
00712 
00713             fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
00714             if( fac > 0.0 ){
00715                fac = 255.0 / fac ;
00716                EDIT_coerce_scale_type( nvox,fac ,
00717                                        MRI_float,fout[kk] , MRI_byte,bout ) ;
00718                fac = 1.0 / fac ;
00719             }
00720 
00721             free( fout[kk] ) ;  /* don't need this anymore */
00722 
00723             /*-- put output brick into dataset, and store scale factor --*/
00724 
00725             EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
00726             tar[kk] = fac ;
00727 
00728             perc = (100 * kk) / nfreq ;
00729             PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
00730          }
00731 
00732          /*-- save scale factor array into dataset --*/
00733 
00734          EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
00735       }
00736       break ;
00737 
00738    } /* end of switch on output data type */
00739 
00740    /*-------------- Cleanup and go home ----------------*/
00741 
00742    PLUTO_set_meter( plint , 100 ) ;  /* set progress meter to 100% */
00743 
00744    PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00745 
00746    FREE_WORKSPACE ;
00747    return NULL ;  /* null string returned means all was OK */
00748 }

   
   { "shortest", "32" ,  "48" ,   "60" ,   "64" ,   "80" ,
                         "96" ,  "120" ,  "128" ,  "160" ,
                        "192" ,  "240" ,  "256" ,  "320" ,
                        "384" ,  "480" ,  "512" ,  "640" ,
                        "768" ,  "960" , "1024" , "1280" ,
                       "1536" , "1920" , "2048"            }

  " Purpose: Compute 'Power Spectrum' of a 3D+time dataset.\n"
  " Input items are:\n"
  "   Input = 3D+time dataset to analyze\n"
  "\n"
  "   Output: Prefix = Filename prefix for new dataset\n"
  "           Datum  = How to store results\n"
  "\n"
  "   Ignore Count   = How many points to ignore at start\n"
  "   Taper Percent  = Amount of data to taper (Hamming)\n"
  "   FFT Length     = Fourier transform size to use [N]\n"
  "                    (If N > size of data, data will be zero)\n"
  "                    (padded. 'shortest' means to use N just)\n"
  "                    (above the length of the time series.  )\n"
  "\n"
  " The output dataset will be stored in the 3D+time format, with\n"
  " the 'time' index actually being frequency.  The frequency grid\n"
  " spacing will be 1/(N*dt), where N=FFT length and dt = input\n"
  " dataset time spacing.\n"
  "\n"
  " The method used is the simplest known: squared periodogram.\n"
  " A single FFT is done (i.e., each point has DOF=2.)\n"