/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#include "mrilib.h"
#define UMIN 4
#define TOCX 1
#define FROMCX 2
int main( int argc , char * argv[] )
{
MRI_IMAGE * inim , * outim ;
int ii , jj , nx,nfft=0,ny , nopt,nby2 , ignore=0,nxi , use=0 ;
complex * cxar ;
float * iar , * oar , * far ;
int nodetrend=0 , cxop=0 ; /* 29 Nov 1999 */
int hilbert=0 ; /* 09 Dec 1999 */
/*-- help? --*/
if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
printf("Usage: 1dfft [options] infile outfile\n"
"where infile is an AFNI *.1D file (ASCII list of numbers arranged\n"
"in columns); outfile will be a similar file, with the absolute\n"
"value of the FFT of the input columns. The length of the file\n"
"will be 1+(FFT length)/2.\n"
"\n"
"Options:\n"
" -ignore sss = Skip the first 'sss' lines in the input file.\n"
" [default = no skipping]\n"
" -use uuu = Use only 'uuu' lines of the input file.\n"
" [default = use them all, Frank]\n"
" -nfft nnn = Set FFT length to 'nnn'.\n"
" [default = length of data (# of lines used)]\n"
" -tocx = Save Re and Im parts of transform in 2 columns.\n"
" -fromcx = Convert 2 column complex input into 1 column\n"
" real output.\n"
" -hilbert = When -fromcx is used, the inverse FFT will\n"
" do the Hilbert transform instead.\n"
" -nodetrend = Skip the detrending of the input.\n"
"\n"
"Nota Bene:\n"
" * Each input time series has any quadratic trend of the\n"
" form 'a+b*t+c*t*t' removed before the FFT, where 't'\n"
" is the line number.\n"
" * The FFT length will be a power-of-2 times at most one\n"
" factor of 3 and one factor of 5. The smallest such\n"
" length >= to the specified FFT length will be used.\n"
" * If the FFT length is longer than the file length, the\n"
" data is zero-padded to make up the difference.\n"
" * Do NOT call the output of this program the Power Spectrum!\n"
" That is something else entirely.\n"
" * If 'outfile' is '-', the output appears on stdout.\n"
) ;
exit(0) ;
}
machdep() ;
nopt = 1 ;
while( nopt < argc && argv[nopt][0] == '-' ){
if( strncmp(argv[nopt],"-nodetrend",6) == 0 ){ /* 29 Nov 1999 */
nodetrend++ ;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-hilbert") == 0 ){ /* 09 Dec 1999 */
hilbert = 1 ;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-tocx") == 0 ){ /* 29 Nov 1999 */
cxop = TOCX ;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-fromcx") == 0 ){ /* 29 Nov 1999 */
cxop = FROMCX ;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-nfft") == 0 ){
if( ++nopt >= argc ){
fprintf(stderr,"** Need argument after -nfft!\n");exit(1);
}
nfft = strtol( argv[nopt] , NULL , 10 ) ;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-ignore") == 0 ){
if( ++nopt >= argc ){
fprintf(stderr,"** Need argument after -ignore!\n");exit(1);
}
ignore = strtol( argv[nopt] , NULL , 10 ) ;
if( ignore < 0 ){
fprintf(stderr,"** Illegal value after -ignore!\n");exit(1);
}
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-use") == 0 ){
if( ++nopt >= argc ){
fprintf(stderr,"** Need argument after -use!\n");exit(1);
}
use = strtol( argv[nopt] , NULL , 10 ) ;
if( use < UMIN ){
fprintf(stderr,"** Illegal value after -use!\n");exit(1);
}
nopt++ ; continue ;
}
fprintf(stderr,"** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
}
if( nopt >= argc ){
fprintf(stderr,"** Need input filenames!\n");exit(1);
}
if( argc > nopt+1 && !THD_filename_ok(argv[nopt+1]) ){
fprintf(stderr,"** Illegal output filename!\n"); exit(1);
}
if( argc > nopt+1 && strcmp(argv[nopt+1],"-") != 0 && THD_is_file(argv[nopt+1]) ){
fprintf(stderr,"** Output file already exists!\n"); exit(1);
}
if( hilbert && cxop != FROMCX ){
fprintf(stderr,"** -hilbert is illegal without -fromcx!\n"); exit(1);
}
/* read input file */
inim = mri_read_1D( argv[nopt] ) ;
if( inim == NULL ){
fprintf(stderr,"** Can't read input file!\n"); exit(1);
}
nx = inim->nx ; ny = inim->ny ;
nxi = nx - ignore ;
if( use > 0 ){
if( use > nxi ){
fprintf(stderr,
"** Not enough data in file for ignore=%d and use=%d\n",
ignore , use ) ;
exit(1) ;
}
nxi = use ;
} else if( nxi < UMIN ){
fprintf(stderr,"** Input file has too few rows!\n"); exit(1);
}
/*- 29 Nov 1999: complex operations are finicky -*/
switch( cxop ){
case TOCX:
if( ny != 1 ){
fprintf(stderr,
"** -tocx option only works on 1 column files!\n") ;
exit(1) ;
}
break ;
case FROMCX:
if( ny != 2 ){
fprintf(stderr,
"** -fromcx option only works on 2 column files!\n") ;
exit(1) ;
}
if( nx != nxi ){
fprintf(stderr,
"** -fromcx option doesn't allow -ignore or -use!\n") ;
exit(1) ;
}
jj = csfft_nextup_one35( 2*(nx-1) ) ;
if( jj != 2*(nx-1) ){
fprintf(stderr,
"** -fromcx only works with certain file lengths!\n") ;
exit(1) ;
}
break ;
}
/* 29 Nov 1999: two possible paths */
if( cxop != FROMCX ){ /* real input */
if( nfft < nxi ) nfft = nxi ;
nfft = csfft_nextup_one35(nfft) ;
fprintf(stderr,"++ FFT length = %d\n",nfft) ;
nby2 = nfft/2 ;
cxar = (complex *) malloc( sizeof(complex) * nfft ) ;
if( cxop == TOCX )
outim = mri_new( nby2+1 , 2 , MRI_float ) ; /* complex output */
else
outim = mri_new( nby2+1 , ny , MRI_float ) ; /* abs() output */
iar = MRI_FLOAT_PTR(inim) ;
oar = MRI_FLOAT_PTR(outim) ;
for( jj=0 ; jj < ny ; jj++ ){
far = iar + jj*nx ;
if( !nodetrend ){
float f0,f1,f2 ;
THD_quadratic_detrend( nxi , far+ignore , &f0,&f1,&f2 ) ;
fprintf(stderr,"++ quadratic trend: %g + %g * i + %g * i*i\n"
" mid = %g end = %g\n" ,
f0,f1,f2 , f0+0.5*nxi*f1+0.25*nxi*nxi*f2 ,
f0+(nxi-1)*f1+(nxi-1)*(nxi-1)*f2 ) ;
}
for( ii=0 ; ii < nxi ; ii++ ){
cxar[ii].r = far[ii+ignore]; cxar[ii].i = 0.0;
}
for( ii=nxi ; ii < nfft ; ii++ ){ cxar[ii].r = cxar[ii].i = 0.0; }
csfft_cox( -1 , nfft , cxar ) ;
far = oar + jj*(nby2+1) ;
if( cxop == TOCX )
for( ii=0 ; ii <= nby2 ; ii++ ){ far[ii] = cxar[ii].r ;
far[ii+(nby2+1)] = cxar[ii].i ; }
else
for( ii=0 ; ii <= nby2 ; ii++ ){ far[ii] = CABS(cxar[ii]) ; }
}
} else { /* complex input */
nfft = 2*(nx-1) ;
nby2 = nfft/2 ;
fprintf(stderr,"++ FFT length = %d\n",nfft) ;
cxar = (complex *) malloc( sizeof(complex) * nfft ) ;
outim = mri_new( nfft , 1 , MRI_float ) ;
iar = MRI_FLOAT_PTR(inim) ;
oar = MRI_FLOAT_PTR(outim) ;
cxar[0].r = iar[0] ; cxar[0].i = 0.0 ;
cxar[nby2].r = iar[nx-1] ; cxar[nby2].i = 0.0 ;
for( ii=1 ; ii < nby2 ; ii++ ){
cxar[ii].r = iar[ii] ; cxar[ii].i = iar[ii+nx] ;
cxar[nfft-ii].r = iar[ii] ; cxar[nfft-ii].i = -iar[ii+nx] ;
}
if( hilbert ){
float val ;
cxar[0].r = cxar[0].i = cxar[nby2].r = cxar[nby2].i = 0.0 ;
for( ii=1 ; ii < nby2 ; ii++ ){
val = cxar[ii].r ; cxar[ii].r = -cxar[ii].i ;
cxar[ii].i = val ;
val = cxar[nfft-ii].r ; cxar[nfft-ii].r = cxar[nfft-ii].i ;
cxar[nfft-ii].i = -val ;
}
}
csfft_cox( 1 , nfft , cxar ) ;
for( ii=0 ; ii < nfft ; ii++ ) oar[ii] = cxar[ii].r / nfft ;
}
mri_write_1D( (argc > nopt+1) ? argv[nopt+1] : "-" , outim ) ;
exit(0) ;
}
syntax highlighted by Code2HTML, v. 0.9.1