1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 #define UMIN 4
10
11 #define TOCX 1
12 #define FROMCX 2
13
main(int argc,char * argv[])14 int main( int argc , char * argv[] )
15 {
16 MRI_IMAGE *inim , *outim ;
17 int ii , jj , nx,nfft=0,ny , nopt,nby2 , ignore=0,nxi , use=0 ;
18 complex *cxar ;
19 float *iar , *oar , *far ;
20 int nodetrend=0 , cxop=0 ; /* 29 Nov 1999 */
21 int hilbert=0 ; /* 09 Dec 1999 */
22
23 /*-- help? --*/
24
25 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
26 printf("\n"
27 "Usage: 1dfft [options] infile outfile\n"
28 "where infile is an AFNI *.1D file (ASCII list of numbers arranged\n"
29 "in columns); outfile will be a similar file, with the absolute\n"
30 "value of the FFT of the input columns. The length of the file\n"
31 "will be 1+(FFT length)/2.\n"
32 "\n"
33 "Options:\n"
34 " -ignore sss = Skip the first 'sss' lines in the input file.\n"
35 " [default = no skipping]\n"
36 " -use uuu = Use only 'uuu' lines of the input file.\n"
37 " [default = use them all, Frank]\n"
38 " -nfft nnn = Set FFT length to 'nnn'.\n"
39 " [default = length of data (# of lines used)]\n"
40 " -tocx = Save Re and Im parts of transform in 2 columns.\n"
41 " -fromcx = Convert 2 column complex input into 1 column\n"
42 " real output.\n"
43 " [-fromcx will not work if the original]\n"
44 " [data FFT length was an odd number! :(]\n"
45 " -hilbert = When -fromcx is used, the inverse FFT will\n"
46 " do the Hilbert transform instead.\n"
47 " -nodetrend = Skip the detrending of the input.\n"
48 "\n"
49 "Nota Bene:\n"
50 " * Each input time series has any quadratic trend of the\n"
51 " form 'a+b*t+c*t*t' removed before the FFT, where 't'\n"
52 " is the line number.\n"
53 #if 0
54 " * The FFT length will be a power-of-2 times at most one\n"
55 " factor of 3 and one factor of 5. The smallest such\n"
56 " length >= to the specified FFT length will be used.\n"
57 #else
58 " * The FFT length can be any positive even integer, but\n"
59 " the Fast Fourier Transform algorithm will be slower if\n"
60 " any prime factors of the FFT length are large (say > 997)\n"
61 " Unless you are applying this program to VERY long files,\n"
62 " this slowdown will probably not be appreciable.\n"
63 #endif
64 " * If the FFT length is longer than the file length, the\n"
65 " data is zero-padded to make up the difference.\n"
66 " * Do NOT call the output of this program the Power Spectrum!\n"
67 " That is something else entirely.\n"
68 " * If 'outfile' is '-' (or missing), the output appears on stdout.\n"
69 ) ;
70 PRINT_COMPILE_DATE ; exit(0) ;
71 }
72
73 machdep() ;
74
75 if( AFNI_yesenv("AFNI_USE_FFTN") ) csfft_force_fftn(1) ;
76
77 nopt = 1 ;
78 while( nopt < argc && argv[nopt][0] == '-' ){
79
80 if( strncmp(argv[nopt],"-nodetrend",6) == 0 ){ /* 29 Nov 1999 */
81 nodetrend++ ;
82 nopt++ ; continue ;
83 }
84
85 if( strcmp(argv[nopt],"-hilbert") == 0 ){ /* 09 Dec 1999 */
86 hilbert = 1 ;
87 nopt++ ; continue ;
88 }
89
90 if( strcmp(argv[nopt],"-tocx") == 0 ){ /* 29 Nov 1999 */
91 cxop = TOCX ;
92 nopt++ ; continue ;
93 }
94
95 if( strcmp(argv[nopt],"-fromcx") == 0 ){ /* 29 Nov 1999 */
96 cxop = FROMCX ;
97 nopt++ ; continue ;
98 }
99
100 if( strcmp(argv[nopt],"-nfft") == 0 ){
101 if( ++nopt >= argc ){
102 fprintf(stderr,"** Need argument after -nfft!\n");exit(1);
103 }
104 nfft = strtol( argv[nopt] , NULL , 10 ) ;
105 nopt++ ; continue ;
106 }
107
108 if( strcmp(argv[nopt],"-ignore") == 0 ){
109 if( ++nopt >= argc ){
110 fprintf(stderr,"** Need argument after -ignore!\n");exit(1);
111 }
112 ignore = strtol( argv[nopt] , NULL , 10 ) ;
113 if( ignore < 0 ){
114 fprintf(stderr,"** Illegal value after -ignore!\n");exit(1);
115 }
116 nopt++ ; continue ;
117 }
118
119 if( strcmp(argv[nopt],"-use") == 0 ){
120 if( ++nopt >= argc ){
121 fprintf(stderr,"** Need argument after -use!\n");exit(1);
122 }
123 use = strtol( argv[nopt] , NULL , 10 ) ;
124 if( use < UMIN ){
125 fprintf(stderr,"** Illegal value after -use!\n");exit(1);
126 }
127 nopt++ ; continue ;
128 }
129
130 fprintf(stderr,"** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
131 }
132
133 if( nopt >= argc ){
134 fprintf(stderr,"** Need input filenames!\n");exit(1);
135 }
136
137 if( argc > nopt+1 && !THD_filename_ok(argv[nopt+1]) ){
138 fprintf(stderr,"** Illegal output filename!\n"); exit(1);
139 }
140 if( argc > nopt+1 && strcmp(argv[nopt+1],"-") != 0 && THD_is_file(argv[nopt+1]) ){
141 fprintf(stderr,"** Output file already exists!\n"); exit(1);
142 }
143
144 if( hilbert && cxop != FROMCX ){
145 fprintf(stderr,"** -hilbert is illegal without -fromcx!\n"); exit(1);
146 }
147
148 /* read input file */
149
150 inim = mri_read_1D( argv[nopt] ) ;
151 if( inim == NULL ){
152 fprintf(stderr,"** Can't read input file!\n"); exit(1);
153 }
154
155 nx = inim->nx ; ny = inim->ny ;
156 nxi = nx - ignore ;
157 if( use > 0 ){
158 if( use > nxi ){
159 fprintf(stderr,
160 "** Not enough data in file for ignore=%d and use=%d\n",
161 ignore , use ) ;
162 exit(1) ;
163 }
164 nxi = use ;
165 } else if( nxi < UMIN ){
166 fprintf(stderr,"** Input file has too few rows!\n"); exit(1);
167 }
168
169 /*- 29 Nov 1999: complex operations are finicky -*/
170
171 switch( cxop ){
172 case TOCX:
173 if( ny != 1 ){
174 fprintf(stderr,
175 "** -tocx option only works on 1 column files!\n") ;
176 exit(1) ;
177 }
178 break ;
179
180 case FROMCX:
181 if( ny != 2 ){
182 fprintf(stderr,
183 "** -fromcx option only works on 2 column files!\n") ;
184 exit(1) ;
185 }
186 if( nx != nxi ){
187 fprintf(stderr,
188 "** -fromcx option doesn't allow -ignore or -use!\n") ;
189 exit(1) ;
190 }
191 jj = csfft_nextup_even( 2*(nx-1) ) ;
192 if( jj != 2*(nx-1) ){
193 fprintf(stderr,
194 "** -fromcx only works with certain file lengths!\n") ;
195 exit(1) ;
196 }
197 break ;
198 }
199
200 /* 29 Nov 1999: two possible paths */
201
202 if( cxop != FROMCX ){ /* real input */
203 if( nfft < nxi ) nfft = nxi ;
204 /* nfft = csfft_nextup_even(nfft) ; */
205 fprintf(stderr,"++ 1dfft length = %d\n",nfft) ;
206 nby2 = nfft/2 ;
207
208 cxar = (complex *) malloc( sizeof(complex) * nfft ) ;
209 if( cxop == TOCX )
210 outim = mri_new( nby2+1 , 2 , MRI_float ) ; /* complex output */
211 else
212 outim = mri_new( nby2+1 , ny , MRI_float ) ; /* abs() output */
213
214 iar = MRI_FLOAT_PTR(inim) ;
215 oar = MRI_FLOAT_PTR(outim) ;
216
217 for( jj=0 ; jj < ny ; jj++ ){
218 far = iar + jj*nx ;
219 if( !nodetrend ){
220 float f0,f1,f2 ;
221 THD_quadratic_detrend( nxi , far+ignore , &f0,&f1,&f2 ) ;
222 #if 0
223 fprintf(stderr,"++ quadratic trend: %g + %g * i + %g * i*i\n"
224 " mid = %g end = %g\n" ,
225 f0,f1,f2 , f0+0.5*nxi*f1+0.25*nxi*nxi*f2 ,
226 f0+(nxi-1)*f1+(nxi-1)*(nxi-1)*f2 ) ;
227 #endif
228 }
229 for( ii=0 ; ii < nxi ; ii++ ){
230 cxar[ii].r = far[ii+ignore]; cxar[ii].i = 0.0;
231 }
232 for( ii=nxi ; ii < nfft ; ii++ ){ cxar[ii].r = cxar[ii].i = 0.0; }
233 csfft_cox( -1 , nfft , cxar ) ;
234 far = oar + jj*(nby2+1) ;
235 if( cxop == TOCX )
236 for( ii=0 ; ii <= nby2 ; ii++ ){ far[ii] = cxar[ii].r ;
237 far[ii+(nby2+1)] = cxar[ii].i ; }
238 else
239 for( ii=0 ; ii <= nby2 ; ii++ ){ far[ii] = CABS(cxar[ii]) ; }
240 }
241
242 } else { /* complex input */
243 nfft = 2*(nx-1) ;
244 nby2 = nfft/2 ;
245 fprintf(stderr,"++ FFT length = %d\n",nfft) ;
246 cxar = (complex *) malloc( sizeof(complex) * nfft ) ;
247
248 outim = mri_new( nfft , 1 , MRI_float ) ;
249
250 iar = MRI_FLOAT_PTR(inim) ;
251 oar = MRI_FLOAT_PTR(outim) ;
252
253 cxar[0].r = iar[0] ; cxar[0].i = 0.0 ;
254 cxar[nby2].r = iar[nx-1] ; cxar[nby2].i = 0.0 ;
255 for( ii=1 ; ii < nby2 ; ii++ ){
256 cxar[ii].r = iar[ii] ; cxar[ii].i = iar[ii+nx] ;
257 cxar[nfft-ii].r = iar[ii] ; cxar[nfft-ii].i = -iar[ii+nx] ;
258 }
259
260 if( hilbert ){
261 float val ;
262 cxar[0].r = cxar[0].i = cxar[nby2].r = cxar[nby2].i = 0.0 ;
263 for( ii=1 ; ii < nby2 ; ii++ ){
264 val = cxar[ii].r ; cxar[ii].r = -cxar[ii].i ;
265 cxar[ii].i = val ;
266 val = cxar[nfft-ii].r ; cxar[nfft-ii].r = cxar[nfft-ii].i ;
267 cxar[nfft-ii].i = -val ;
268 }
269 }
270 csfft_cox( 1 , nfft , cxar ) ;
271 for( ii=0 ; ii < nfft ; ii++ ) oar[ii] = cxar[ii].r / nfft ;
272 }
273
274 mri_write_1D( (argc > nopt+1) ? argv[nopt+1] : "-" , outim ) ;
275 exit(0) ;
276 }
277