1 #include "mrilib.h"
2
main(int argc,char * argv[])3 int main( int argc , char * argv[] )
4 {
5 MRI_IMAGE *inim , *ortim=NULL ;
6 float **vec , **ort=NULL ; int iv , nx,ny , nopt, nort=0 , do_norm=0 ;
7 int qdet=2 ; float dt=1.0f , fbot,ftop ;
8
9 /*-- help? --*/
10
11 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
12 printf("Usage: 1dBandpass [options] fbot ftop infile ~1~\n"
13 "\n"
14 " * infile is an AFNI *.1D file; each column is processed\n"
15 " * fbot = lowest frequency in the passband, in Hz\n"
16 " [can be 0 if you want to do a lowpass filter only,]\n"
17 " but the mean and Nyquist freq are always removed ]\n"
18 " * ftop = highest frequency in the passband (must be > fbot)\n"
19 " [if ftop > Nyquist freq, then we have a highpass filter only]\n"
20 " * You cannot construct a 'notch' filter with this program!\n"
21 " * Output vectors appear on stdout; redirect as desired\n"
22 " * Program will fail if fbot and ftop are too close for comfort\n"
23 " * The actual FFT length used will be printed, and may be larger\n"
24 " than the input time series length for the sake of efficiency.\n"
25 "\n"
26 "Options: ~1~\n"
27 " -dt dd = set time step to 'dd' sec [default = 1.0]\n"
28 " -ort f.1D = Also orthogonalize input to columns in f.1D\n"
29 " [only one '-ort' option is allowed]\n"
30 " -nodetrend = Skip the quadratic detrending of the input\n"
31 " -norm = Make output time series have L2 norm = 1\n"
32 "\n"
33 "Example: ~1~\n"
34 " 1deval -num 1000 -expr 'gran(0,1)' > r1000.1D\n"
35 " 1dBandpass 0.025 0.20 r1000.1D > f1000.1D\n"
36 " 1dfft f1000.1D - | 1dplot -del 0.000977 -stdin -plabel 'Filtered |FFT|'\n"
37 "\n"
38 "Goal: ~1~\n"
39 " * Mostly to test the functions in thd_bandpass.c -- RWCox -- May 2009\n"
40 ) ;
41 PRINT_COMPILE_DATE ; exit(0) ;
42 }
43
44 machdep() ;
45
46 nopt = 1 ;
47 while( nopt < argc && argv[nopt][0] == '-' ){
48
49 if( strcmp(argv[nopt],"-norm") == 0 ){
50 do_norm = 1 ; nopt++ ; continue ;
51 }
52
53 if( strcmp(argv[nopt],"-ort") == 0 ){
54 if( ++nopt >= argc ) ERROR_exit("need an argument after -ort!") ;
55 if( ortim != NULL ) ERROR_exit("can't have 2 -ort options!") ;
56 ortim = mri_read_1D( argv[nopt] ) ;
57 if( ortim == NULL ) ERROR_exit("can't read from -ort '%s'",argv[nopt]) ;
58 nopt++ ; continue ;
59 }
60
61 if( strncmp(argv[nopt],"-nodetrend",6) == 0 ){
62 qdet = 0 ; nopt++ ; continue ;
63 }
64
65 if( strcmp(argv[nopt],"-dt") == 0 ){
66 if( ++nopt >= argc ) ERROR_exit("need an argument after -dt!") ;
67 dt = (float)strtod(argv[nopt],NULL) ;
68 if( dt <= 0.0f ){
69 WARNING_message("value after -dt illegal: setting time step to 1.0") ;
70 dt = 1.0f ;
71 }
72 nopt++ ; continue ;
73 }
74
75 fprintf(stderr,"** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
76 }
77
78 /** check inputs for reasonablositiness **/
79
80 if( nopt+2 >= argc ) ERROR_exit("Need more arguments on command line!") ;
81
82 fbot = (float)strtod(argv[nopt++],NULL) ;
83 if( fbot < 0.0f ) ERROR_exit("fbot value can't be negative!") ;
84
85 ftop = (float)strtod(argv[nopt++],NULL) ;
86 if( ftop <= fbot ) ERROR_exit("ftop value must be greater than fbot value!") ;
87
88 inim = mri_read_1D(argv[nopt]) ;
89 if( inim == NULL ) ERROR_exit("Can't read 1D file '%s'",argv[nopt]) ;
90 nx = inim->nx ; ny = inim->ny ;
91 if( nx < 9 ) ERROR_exit("1D file '%s' has only %d lines",argv[nopt],nx) ;
92
93 if( !THD_bandpass_OK(nx,dt,fbot,ftop,1) ) ERROR_exit("Can't continue!") ;
94
95 /* make list of pointers to start of each input vector */
96
97 vec = (float **)malloc(sizeof(float *)*ny) ;
98 for( iv=0 ; iv < ny ; iv++ ) vec[iv] = MRI_FLOAT_PTR(inim) + iv*nx ;
99
100 /* similarly for the ort vectors */
101
102 if( ortim != NULL ){
103 if( ortim->nx != nx )
104 ERROR_exit("-ort file and input 1D file differ in column lengths!") ;
105 nort = ortim->ny ;
106 ort = (float **)malloc(sizeof(float *)*nort) ;
107 for( iv=0 ; iv < nort ; iv++ ) ort[iv] = MRI_FLOAT_PTR(ortim) + iv*nx ;
108 }
109
110 /* all the real work now */
111
112 (void)THD_bandpass_vectors( nx, ny, vec, dt, fbot,ftop, qdet, nort, ort ) ;
113
114 if( do_norm ){
115 for( iv=0 ; iv < ny ; iv++ ) THD_normalize( nx , vec[iv] ) ;
116 }
117
118 /* write to stdout */
119
120 mri_write_1D( "-" , inim ) ;
121 exit(0) ;
122 }
123