26 #include "dft_cmdline.h" 34 #define creal(A) (A[0]) 35 #define cimag(A) (A[1]) 36 #define cmag(A) sqrt(A[0]*A[0]+A[1]*A[1]) 50 do_phase(
double **in,
int n, fftw_plan fftp, fftw_plan bftp,
51 fftw_complex **out,
int t )
53 #pragma omp parallel for 54 for(
int i=0; i<n; i++ ) {
56 int thr = omp_get_thread_num();
60 memset( out[thr], 0,
sizeof(fftw_complex)*t );
61 fftw_execute_dft_r2c( fftp, in[i], out[thr] );
64 fftw_execute_dft( bftp, out[thr], out[thr] );
65 for(
int j=0; j<t; j++ )
66 in[i][j] = atan2(
cimag(out[thr][j]),
creal(out[thr][j]) );
85 bp_filt(
double **in,
int n, fftw_plan fftp, fftw_plan bftp,
86 fftw_complex **out,
int t,
float dt,
float low,
float high )
88 int mini = low*dt*(t/2)/500.;
89 int maxi = high ? high*dt*(t/2)/500. : t;
91 #pragma omp parallel for 92 for(
int i=0; i<n; i++ ) {
94 int thr = omp_get_thread_num();
98 fftw_execute_dft_r2c( fftp, in[i], out[thr] );
99 for(
int j=0; j<=mini; j++ )
101 for(
int j=maxi; j<=t/2; j++ )
103 fftw_execute_dft_c2r( bftp, out[thr], in[i] );
105 for(
int j=0; j<t; j++ )
125 bp_notch(
double **in,
int n, fftw_plan fftp, fftw_plan bftp,
126 fftw_complex **out,
int t,
float dt,
float notch,
int nwidth )
128 int notchi = notch*dt*(t/2)/500.;
130 #pragma omp parallel for 131 for(
int i=0; i<n; i++ ) {
133 int thr = omp_get_thread_num();
137 fftw_execute_dft_r2c( fftp, in[i], out[thr] );
138 for(
int j=notchi-nwidth; j<=notchi+nwidth; j++ )
140 fftw_execute_dft_c2r( bftp, out[thr], in[i] );
142 for(
int j=0; j<t; j++ )
158 do_PSD(
double **in,
int n, fftw_plan ftp,
159 fftw_complex **out,
int t,
float dt )
161 #pragma omp parallel for 162 for(
int i=0; i<n; i++ ) {
164 int thr = omp_get_thread_num();
168 fftw_execute_dft_r2c( ftp, in[i], out[thr] );
169 for(
int j=0; j<=t/2; j++ ) {
170 in[i][j] =
creal(out[thr][j])*
creal(out[thr][j])+
188 fftw_complex **out,
float t,
float dt )
190 #pragma omp parallel for 191 for(
int i=0; i<n; i++ ) {
193 int thr = omp_get_thread_num();
198 fftw_execute_dft_r2c( ftp, in[i], out[thr] );
201 float maxmag =
creal(out[thr][1])*
creal(out[thr][1])+
203 for(
int j=2; j<t/2; j++ ) {
204 float tmag =
creal(out[thr][j])*
creal(out[thr][j])+
206 if( tmag > maxmag ) {
211 in[i][0] = 500./dt/(t/2.)*maxf;
228 FILE *fp =
static_cast<FILE *
>(ohead->
fileptr());
230 for(
long t=0; t<ohead->
t(); t++ ) {
232 for(
int i=0; i<nr; i++ ){
233 float d = data[i][t];
234 fwrite( &d,
sizeof(
float), 1, fp );
240 int main(
int argc,
char *argv[] )
243 if(!getenv(
"OMP_NUM_THREADS") ){
244 cerr <<
"unst" << endl;
245 omp_set_num_threads(1);
249 gengetopt_args_info opts;
251 if (cmdline_parser(argc, argv, &opts) != 0)
254 if( !opts.inputs_num ) {
255 cmdline_parser_print_help();
259 gzFile in = gzopen((!opts.inputs_num||!strcmp(opts.inputs[0],
"-"))?
"/dev/stdin":opts.inputs[0],
"r");
261 cerr <<
"Cannot open input file: " << opts.inputs[0] << endl;
267 char tname[] =
"/tmp/igbdftXXXXXX";
269 if( !strcmp( opts.output_arg,
"." ) ||
270 !strcmp( opts.output_arg,
"-") ) {
274 out_fname = opts.output_arg;
276 dftout = fopen( out_fname.c_str(),
"w" );
278 cerr <<
"Cannot open output file: " << out_fname << endl;
286 int npar =
min(opts.numpar_arg,(
int)orig.
slice_sz());
287 double **data =
new double *[npar];
288 int aligned_size =((orig.
t()+7)/8)*8;
289 double *bbuf = fftw_alloc_real(npar*aligned_size);
290 for(
int i=0; i<npar; i++ )
291 data[i] = bbuf+aligned_size*i;
298 int numthr=omp_get_max_threads();
302 fftw_complex **out =
new fftw_complex*[numthr];
304 switch( opts.op_arg ) {
307 for(
int i=0; i<numthr; i++ )
308 out[i] = fftw_alloc_complex(orig.
t()/2+1);
309 if( !(ftp = fftw_plan_dft_r2c_1d( orig.
t(), data[0], out[0], 0)) ) {
310 cerr <<
"Cannot create fftw3 plan!" << endl;
313 iftp = fftw_plan_dft_c2r_1d( orig.
t(), out[0], data[0], 0);
317 dfthead.
t(orig.
t()/2+1);
320 for(
int i=0; i<numthr; i++ )
321 out[i] = fftw_alloc_complex(orig.
t()/2+1);
322 if( !(ftp = fftw_plan_dft_r2c_1d( orig.
t(), data[0], out[0], 0)) ) {
323 cerr <<
"Cannot create fftw3 plan!" << endl;
330 for(
int i=0; i<numthr; i++ )
331 out[i] = fftw_alloc_complex(orig.
t()/2+1);
332 if( !(ftp = fftw_plan_dft_r2c_1d( orig.
t(), data[0], out[0], 0)) ) {
333 cerr <<
"Cannot create fftw3 plan!" << endl;
338 dfthead.
unites(
"radians");
339 for(
int i=0; i<numthr; i++ )
340 out[i] = fftw_alloc_complex(orig.
t());
341 ftp = fftw_plan_dft_r2c_1d( orig.
t(), data[0], out[0], 0 );
342 iftp = fftw_plan_dft_1d( orig.
t(), out[0], out[0], FFTW_BACKWARD, 0);
344 case op_arg_Butter_bp :
346 opts.upper_arg*orig.
inc_t()/500. );
347 outf =
new double[orig.
t()];
352 double *slice =
new double[orig.
slice_sz()];
353 float *dftslice =
new float[orig.
slice_sz()];
355 for(
int n=0; n<orig.
slice_sz(); n+=npar ) {
357 cerr <<
"Reading nodes ["<< n <<
"-" << (n+nr-1) <<
"]"<< endl;
358 gzseek( in, 1024, SEEK_SET );
359 for(
int i=0; i<orig.
t(); i++ ) {
361 for(
int j=0; j<nr; j++ )
362 data[j][i] = slice[n+j];
364 switch( opts.op_arg ) {
372 do_phase( data, nr, ftp, iftp, out, orig.
t() );
376 opts.low_arg, opts.upper_arg );
380 opts.low_arg, opts.notch_width_arg );
382 case op_arg_Butter_bp :
393 if( !strcmp( opts.output_arg,
"." ) )
394 rename( out_fname.c_str(), opts.inputs[0] );
395 else if( !strcmp( opts.output_arg,
"-" ) ) {
396 dftout = fopen( out_fname.c_str(),
"r" );
398 while( (c=getc(dftout)) != EOF )
constexpr T min(T a, T b)
void bp_notch(double **in, int n, fftw_plan fftp, fftw_plan bftp, fftw_complex **out, int t, float dt, float notch, int nwidth)
int main(int argc, char *argv[])
void filter_all(double **in, int nr, IIR_Filt *filt, double *out, int t)
void bp_filt(double **in, int n, fftw_plan fftp, fftw_plan bftp, fftw_complex **out, int t, float dt, float low, float high)
void do_domF(double **in, int n, fftw_plan ftp, fftw_complex **out, float t, float dt)
void do_PSD(double **in, int n, fftw_plan ftp, fftw_complex **out, int t, float dt)
void write_traces(T **data, int nr, int n, IGBheader *ohead)
void do_phase(double **in, int n, fftw_plan fftp, fftw_plan bftp, fftw_complex **out, int t)