Doxygen Source Code Documentation
1dfft.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | UMIN 4 |
#define | TOCX 1 |
#define | FROMCX 2 |
Functions | |
int | main (int argc, char *argv[]) |
Define Documentation
|
Definition at line 12 of file 1dfft.c. Referenced by main(). |
|
Definition at line 11 of file 1dfft.c. Referenced by main(). |
|
Definition at line 9 of file 1dfft.c. Referenced by main(). |
Function Documentation
|
convert three sub-briks to a raw dataset with consecutive triplets Definition at line 14 of file 1dfft.c. References argc, CABS, csfft_cox(), csfft_nextup_one35(), far, FROMCX, complex::i, machdep(), malloc, MRI_FLOAT_PTR, mri_new(), mri_read_1D(), mri_write_1D(), MRI_IMAGE::nx, MRI_IMAGE::ny, complex::r, THD_filename_ok(), THD_is_file(), THD_quadratic_detrend(), TOCX, and UMIN.
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 } |