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 #define UMIN 4
00010 
00011 #define TOCX   1
00012 #define FROMCX 2
00013 
00014 int main( int argc , char * argv[] )
00015 {
00016    MRI_IMAGE * inim , * outim ;
00017    int ii , jj , nx,nfft=0,ny , nopt,nby2 , ignore=0,nxi , use=0 ;
00018    complex * cxar ;
00019    float * iar , * oar , * far ;
00020    int nodetrend=0 , cxop=0 ;     /* 29 Nov 1999 */
00021    int hilbert=0 ;                /* 09 Dec 1999 */
00022 
00023    /*-- help? --*/
00024 
00025    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00026      printf("Usage: 1dfft [options] infile outfile\n"
00027             "where infile is an AFNI *.1D file (ASCII list of numbers arranged\n"
00028             "in columns); outfile will be a similar file, with the absolute\n"
00029             "value of the FFT of the input columns.  The length of the file\n"
00030             "will be 1+(FFT length)/2.\n"
00031             "\n"
00032             "Options:\n"
00033             "  -ignore sss = Skip the first 'sss' lines in the input file.\n"
00034             "                [default = no skipping]\n"
00035             "  -use uuu    = Use only 'uuu' lines of the input file.\n"
00036             "                [default = use them all, Frank]\n"
00037             "  -nfft nnn   = Set FFT length to 'nnn'.\n"
00038             "                [default = length of data (# of lines used)]\n"
00039             "  -tocx       = Save Re and Im parts of transform in 2 columns.\n"
00040             "  -fromcx     = Convert 2 column complex input into 1 column\n"
00041             "                  real output.\n"
00042             "  -hilbert    = When -fromcx is used, the inverse FFT will\n"
00043             "                  do the Hilbert transform instead.\n"
00044             "  -nodetrend  = Skip the detrending of the input.\n"
00045             "\n"
00046             "Nota Bene:\n"
00047             " * Each input time series has any quadratic trend of the\n"
00048             "     form 'a+b*t+c*t*t' removed before the FFT, where 't'\n"
00049             "     is the line number.\n"
00050             " * The FFT length will be a power-of-2 times at most one\n"
00051             "     factor of 3 and one factor of 5.  The smallest such\n"
00052             "     length >= to the specified FFT length will be used.\n"
00053             " * If the FFT length is longer than the file length, the\n"
00054             "     data is zero-padded to make up the difference.\n"
00055             " * Do NOT call the output of this program the Power Spectrum!\n"
00056             "     That is something else entirely.\n"
00057             " * If 'outfile' is '-', the output appears on stdout.\n"
00058            ) ;
00059       exit(0) ;
00060    }
00061 
00062    machdep() ;
00063 
00064    nopt = 1 ;
00065    while( nopt < argc && argv[nopt][0] == '-' ){
00066 
00067       if( strncmp(argv[nopt],"-nodetrend",6) == 0 ){ /* 29 Nov 1999 */
00068          nodetrend++ ;
00069          nopt++ ; continue ;
00070       }
00071 
00072       if( strcmp(argv[nopt],"-hilbert") == 0 ){      /* 09 Dec 1999 */
00073          hilbert = 1 ;
00074          nopt++ ; continue ;
00075       }
00076 
00077       if( strcmp(argv[nopt],"-tocx") == 0 ){         /* 29 Nov 1999 */
00078          cxop = TOCX ;
00079          nopt++ ; continue ;
00080       }
00081 
00082       if( strcmp(argv[nopt],"-fromcx") == 0 ){       /* 29 Nov 1999 */
00083          cxop = FROMCX ;
00084          nopt++ ; continue ;
00085       }
00086 
00087       if( strcmp(argv[nopt],"-nfft") == 0 ){
00088          if( ++nopt >= argc ){
00089             fprintf(stderr,"** Need argument after -nfft!\n");exit(1);
00090          }
00091          nfft = strtol( argv[nopt] , NULL , 10 ) ;
00092          nopt++ ; continue ;
00093       }
00094 
00095       if( strcmp(argv[nopt],"-ignore") == 0 ){
00096          if( ++nopt >= argc ){
00097             fprintf(stderr,"** Need argument after -ignore!\n");exit(1);
00098          }
00099          ignore = strtol( argv[nopt] , NULL , 10 ) ;
00100          if( ignore < 0 ){
00101             fprintf(stderr,"** Illegal value after -ignore!\n");exit(1);
00102          }
00103          nopt++ ; continue ;
00104       }
00105 
00106       if( strcmp(argv[nopt],"-use") == 0 ){
00107          if( ++nopt >= argc ){
00108             fprintf(stderr,"** Need argument after -use!\n");exit(1);
00109          }
00110          use = strtol( argv[nopt] , NULL , 10 ) ;
00111          if( use < UMIN ){
00112             fprintf(stderr,"** Illegal value after -use!\n");exit(1);
00113          }
00114          nopt++ ; continue ;
00115       }
00116 
00117       fprintf(stderr,"** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00118    }
00119 
00120    if( nopt >= argc ){
00121       fprintf(stderr,"** Need input filenames!\n");exit(1);
00122    }
00123 
00124    if( argc > nopt+1 && !THD_filename_ok(argv[nopt+1]) ){
00125      fprintf(stderr,"** Illegal output filename!\n"); exit(1);
00126    }
00127    if( argc > nopt+1 && strcmp(argv[nopt+1],"-") != 0 && THD_is_file(argv[nopt+1]) ){
00128      fprintf(stderr,"** Output file already exists!\n"); exit(1);
00129    }
00130 
00131    if( hilbert && cxop != FROMCX ){
00132       fprintf(stderr,"** -hilbert is illegal without -fromcx!\n"); exit(1);
00133    }
00134 
00135    /* read input file */
00136 
00137    inim = mri_read_1D( argv[nopt] ) ;
00138    if( inim == NULL ){
00139       fprintf(stderr,"** Can't read input file!\n"); exit(1);
00140    }
00141 
00142    nx = inim->nx ; ny = inim->ny ;
00143    nxi = nx - ignore ;
00144    if( use > 0 ){
00145       if( use > nxi ){
00146          fprintf(stderr,
00147                  "** Not enough data in file for ignore=%d and use=%d\n",
00148                  ignore , use ) ;
00149          exit(1) ;
00150       }
00151       nxi = use ;
00152    } else if( nxi < UMIN ){
00153       fprintf(stderr,"** Input file has too few rows!\n"); exit(1);
00154    }
00155 
00156    /*- 29 Nov 1999: complex operations are finicky -*/
00157 
00158    switch( cxop ){
00159       case TOCX:
00160          if( ny != 1 ){
00161             fprintf(stderr,
00162                     "** -tocx option only works on 1 column files!\n") ;
00163             exit(1) ;
00164          }
00165       break ;
00166 
00167       case FROMCX:
00168          if( ny != 2 ){
00169             fprintf(stderr,
00170                     "** -fromcx option only works on 2 column files!\n") ;
00171             exit(1) ;
00172          }
00173          if( nx != nxi ){
00174             fprintf(stderr,
00175                    "** -fromcx option doesn't allow -ignore or -use!\n") ;
00176             exit(1) ;
00177          }
00178          jj = csfft_nextup_one35( 2*(nx-1) ) ;
00179          if( jj != 2*(nx-1) ){
00180             fprintf(stderr,
00181                     "** -fromcx only works with certain file lengths!\n") ;
00182             exit(1) ;
00183          }
00184       break ;
00185    }
00186 
00187    /* 29 Nov 1999: two possible paths */
00188 
00189    if( cxop != FROMCX ){                     /* real input */
00190       if( nfft < nxi ) nfft = nxi ;
00191       nfft = csfft_nextup_one35(nfft) ;
00192       fprintf(stderr,"++ FFT length = %d\n",nfft) ;
00193       nby2 = nfft/2 ;
00194 
00195       cxar = (complex *) malloc( sizeof(complex) * nfft ) ;
00196       if( cxop == TOCX )
00197          outim = mri_new( nby2+1 , 2 , MRI_float ) ;   /* complex output */
00198       else
00199          outim = mri_new( nby2+1 , ny , MRI_float ) ;  /* abs() output */
00200 
00201       iar = MRI_FLOAT_PTR(inim) ;
00202       oar = MRI_FLOAT_PTR(outim) ;
00203 
00204       for( jj=0 ; jj < ny ; jj++ ){
00205          far = iar + jj*nx ;
00206          if( !nodetrend ){
00207             float f0,f1,f2 ;
00208             THD_quadratic_detrend( nxi , far+ignore , &f0,&f1,&f2 ) ;
00209             fprintf(stderr,"++ quadratic trend: %g + %g * i + %g * i*i\n"
00210                            "                    mid = %g  end = %g\n" ,
00211                     f0,f1,f2 , f0+0.5*nxi*f1+0.25*nxi*nxi*f2 ,
00212                                f0+(nxi-1)*f1+(nxi-1)*(nxi-1)*f2 ) ;
00213          }
00214          for( ii=0 ; ii < nxi ; ii++ ){
00215             cxar[ii].r = far[ii+ignore]; cxar[ii].i = 0.0;
00216          }
00217          for( ii=nxi ; ii < nfft ; ii++ ){ cxar[ii].r = cxar[ii].i = 0.0; }
00218          csfft_cox( -1 , nfft , cxar ) ;
00219          far = oar + jj*(nby2+1) ;
00220          if( cxop == TOCX )
00221             for( ii=0 ; ii <= nby2 ; ii++ ){ far[ii]          = cxar[ii].r ;
00222                                              far[ii+(nby2+1)] = cxar[ii].i ; }
00223          else
00224             for( ii=0 ; ii <= nby2 ; ii++ ){ far[ii] = CABS(cxar[ii]) ; }
00225       }
00226 
00227    } else {   /* complex input */
00228       nfft = 2*(nx-1) ;
00229       nby2 = nfft/2 ;
00230       fprintf(stderr,"++ FFT length = %d\n",nfft) ;
00231       cxar = (complex *) malloc( sizeof(complex) * nfft ) ;
00232 
00233       outim = mri_new( nfft , 1 , MRI_float ) ;
00234 
00235       iar = MRI_FLOAT_PTR(inim) ;
00236       oar = MRI_FLOAT_PTR(outim) ;
00237 
00238       cxar[0].r    = iar[0]    ; cxar[0].i    = 0.0 ;
00239       cxar[nby2].r = iar[nx-1] ; cxar[nby2].i = 0.0 ;
00240       for( ii=1 ; ii < nby2 ; ii++ ){
00241          cxar[ii].r      = iar[ii] ; cxar[ii].i      =  iar[ii+nx] ;
00242          cxar[nfft-ii].r = iar[ii] ; cxar[nfft-ii].i = -iar[ii+nx] ;
00243       }
00244 
00245       if( hilbert ){
00246          float val ;
00247          cxar[0].r = cxar[0].i = cxar[nby2].r = cxar[nby2].i = 0.0 ;
00248          for( ii=1 ; ii < nby2 ; ii++ ){
00249             val = cxar[ii].r      ; cxar[ii].r      = -cxar[ii].i ;
00250                                     cxar[ii].i      =  val ;
00251             val = cxar[nfft-ii].r ; cxar[nfft-ii].r =  cxar[nfft-ii].i ;
00252                                     cxar[nfft-ii].i = -val ;
00253          }
00254       }
00255       csfft_cox( 1 , nfft , cxar ) ;
00256       for( ii=0 ; ii < nfft ; ii++ ) oar[ii] = cxar[ii].r / nfft ;
00257    }
00258 
00259    mri_write_1D( (argc > nopt+1) ? argv[nopt+1] : "-" , outim ) ;
00260    exit(0) ;
00261 }