openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
IGBdft.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is free software: you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation, either version 3 of the License, or
9 // (at your option) any later version.
10 //
11 // This program is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with this program. If not, see <https://www.gnu.org/licenses/>.
18 // ----------------------------------------------------------------------------
19 
20 #include <fftw3.h>
21 #include <math.h>
22 #include <string>
23 #include <iostream>
24 #include <algorithm>
25 #include "IGBheader.h"
26 #include "dft_cmdline.h"
27 #include "Filter.h"
28 #ifdef _OPENMP
29 #include <omp.h>
30 #endif
31 
32 using namespace std;
33 
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])
37 
49 void
50 do_phase( double **in, int n, fftw_plan fftp, fftw_plan bftp,
51  fftw_complex **out, int t )
52 {
53 #pragma omp parallel for
54  for( int i=0; i<n; i++ ) {
55 #ifdef _OPENMP
56  int thr = omp_get_thread_num();
57 #else
58  int thr = 0;
59 #endif
60  memset( out[thr], 0, sizeof(fftw_complex)*t );
61  fftw_execute_dft_r2c( fftp, in[i], out[thr] );
62  creal(out[thr][0]) = cimag(out[thr][0]) = 0.; // remove dc
63  // this is an inverse DFT which ignores the negative components
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]) );
67  }
68 }
69 
84 void
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 )
87 {
88  int mini = low*dt*(t/2)/500.;
89  int maxi = high ? high*dt*(t/2)/500. : t;
90 
91 #pragma omp parallel for
92  for( int i=0; i<n; i++ ) {
93 #ifdef _OPENMP
94  int thr = omp_get_thread_num();
95 #else
96  int thr = 0;
97 #endif
98  fftw_execute_dft_r2c( fftp, in[i], out[thr] );
99  for( int j=0; j<=mini; j++ )
100  creal(out[thr][j]) = cimag(out[thr][j]) = 0.;
101  for( int j=maxi; j<=t/2; j++ )
102  creal(out[thr][j]) = cimag(out[thr][j]) = 0.;
103  fftw_execute_dft_c2r( bftp, out[thr], in[i] );
104  // doing the transform and its inverse scales by t
105  for( int j=0; j<t; j++ )
106  in[i][j] /= t;
107  }
108 }
109 
124 void
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 )
127 {
128  int notchi = notch*dt*(t/2)/500.;
129 
130 #pragma omp parallel for
131  for( int i=0; i<n; i++ ) {
132 #ifdef _OPENMP
133  int thr = omp_get_thread_num();
134 #else
135  int thr = 0;
136 #endif
137  fftw_execute_dft_r2c( fftp, in[i], out[thr] );
138  for( int j=notchi-nwidth; j<=notchi+nwidth; j++ )
139  creal(out[thr][j]) = cimag(out[thr][j]) = 0.;
140  fftw_execute_dft_c2r( bftp, out[thr], in[i] );
141  // doing the transform and its inverse scales by t
142  for( int j=0; j<t; j++ )
143  in[i][j] /= t;
144  }
145 }
146 
147 
157 void
158 do_PSD( double **in, int n, fftw_plan ftp,
159  fftw_complex **out, int t, float dt )
160 {
161 #pragma omp parallel for
162  for( int i=0; i<n; i++ ) {
163 #ifdef _OPENMP
164  int thr = omp_get_thread_num();
165 #else
166  int thr = 0;
167 #endif
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])+
171  cimag(out[thr][j])*cimag(out[thr][j]);
172  }
173  }
174 }
175 
176 
186 void
187 do_domF( double **in, int n, fftw_plan ftp,
188  fftw_complex **out, float t, float dt )
189 {
190 #pragma omp parallel for
191  for( int i=0; i<n; i++ ) {
192 #ifdef _OPENMP
193  int thr = omp_get_thread_num();
194 #else
195  int thr = 0;
196 #endif
197 
198  fftw_execute_dft_r2c( ftp, in[i], out[thr] );
199 
200  int maxf = 1;
201  float maxmag = creal(out[thr][1])*creal(out[thr][1])+
202  cimag(out[thr][1])*cimag(out[thr][1]);
203  for( int j=2; j<t/2; j++ ) {
204  float tmag = creal(out[thr][j])*creal(out[thr][j])+
205  cimag(out[thr][j])*cimag(out[thr][j]);
206  if( tmag > maxmag ) {
207  maxmag = tmag;
208  maxf = j;
209  }
210  }
211  in[i][0] = 500./dt/(t/2.)*maxf;
212  }
213 }
214 
215 
224 template<class T>
225 void
226 write_traces( T **data, int nr, int n, IGBheader *ohead )
227 {
228  FILE *fp = static_cast<FILE *>(ohead->fileptr());
229 
230  for( long t=0; t<ohead->t(); t++ ) {
231  fseek( fp, 1024+(t*ohead->slice_sz()+n)*ohead->data_size(), SEEK_SET );
232  for( int i=0; i<nr; i++ ){
233  float d = data[i][t];
234  fwrite( &d, sizeof(float), 1, fp );
235  }
236  }
237 }
238 
239 
240 int main(int argc, char *argv[] )
241 {
242 #ifdef _OPENMP
243  if(!getenv("OMP_NUM_THREADS") ){
244  cerr << "unst" << endl;
245  omp_set_num_threads(1);
246  }
247 #endif
248 
249  gengetopt_args_info opts;
250 
251  if (cmdline_parser(argc, argv, &opts) != 0)
252  exit(1);
253 
254  if( !opts.inputs_num ) {
255  cmdline_parser_print_help();
256  exit(0);
257  }
258 
259  gzFile in = gzopen((!opts.inputs_num||!strcmp(opts.inputs[0],"-"))?"/dev/stdin":opts.inputs[0],"r");
260  if( !in ) {
261  cerr << "Cannot open input file: " << opts.inputs[0] << endl;
262  exit(1);
263  }
264  IGBheader orig(in,true);
265 
266  FILE* dftout;
267  char tname[] = "/tmp/igbdftXXXXXX";
268  string out_fname;
269  if( !strcmp( opts.output_arg, "." ) ||
270  !strcmp( opts.output_arg, "-") ) {
271  mkstemp( tname );
272  out_fname = tname;
273  } else
274  out_fname = opts.output_arg;
275 
276  dftout = fopen( out_fname.c_str(), "w" );
277  if( !dftout ) {
278  cerr << "Cannot open output file: " << out_fname << endl;
279  exit(1);
280  }
281 
282  IGBheader dfthead = orig;
283  dfthead.fileptr( dftout );
284  dfthead.type(IGB_FLOAT);
285 
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; // align for SIMD
289  double *bbuf = fftw_alloc_real(npar*aligned_size);
290  for( int i=0; i<npar; i++ )
291  data[i] = bbuf+aligned_size*i;
292 
293  fftw_plan ftp, iftp;
294  IIR_Filt *filt;
295  double *outf;
296 
297 #ifdef _OPENMP
298  int numthr=omp_get_max_threads();
299 #else
300  int numthr=1;
301 #endif
302  fftw_complex **out = new fftw_complex*[numthr];
303 
304  switch( opts.op_arg ) {
305  case op_arg_DFT_bp :
306  case op_arg_notch :
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;
311  exit(1);
312  }
313  iftp = fftw_plan_dft_c2r_1d( orig.t(), out[0], data[0], 0);
314  break;
315  case op_arg_PSD :
316  dfthead.unites_t("Hz");
317  dfthead.t(orig.t()/2+1);
318  dfthead.org_t(0);
319  dfthead.inc_t(1000./orig.t()/orig.inc_t());
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;
324  exit(1);
325  }
326  break;
327  case op_arg_DomF :
328  dfthead.unites("Hz");
329  dfthead.t(1);
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;
334  exit(1);
335  }
336  break;
337  case op_arg_phase :
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);
343  break;
344  case op_arg_Butter_bp :
345  filt = new Butterworth( opts.order_arg, opts.low_arg*orig.inc_t()/500.,
346  opts.upper_arg*orig.inc_t()/500. );
347  outf = new double[orig.t()];
348  break;
349  }
350  dfthead.write();
351 
352  double *slice = new double[orig.slice_sz()];
353  float *dftslice = new float[orig.slice_sz()];
354 
355  for( int n=0; n<orig.slice_sz(); n+=npar ) {
356  int nr = orig.slice_sz()-n>=npar?npar:orig.slice_sz()-n;
357  cerr << "Reading nodes ["<< n << "-" << (n+nr-1) << "]"<< endl;
358  gzseek( in, 1024, SEEK_SET );
359  for( int i=0; i<orig.t(); i++ ) {
360  orig.read_data( slice );
361  for( int j=0; j<nr; j++ )
362  data[j][i] = slice[n+j];
363  }
364  switch( opts.op_arg ) {
365  case op_arg_PSD :
366  do_PSD( data, nr, ftp, out, orig.t(), orig.inc_t() );
367  break;
368  case op_arg_DomF :
369  do_domF( data, nr, ftp, out, orig.t(), orig.inc_t() );
370  break;
371  case op_arg_phase :
372  do_phase( data, nr, ftp, iftp, out, orig.t() );
373  break;
374  case op_arg_DFT_bp:
375  bp_filt( data, nr, ftp, iftp, out, orig.t(), orig.inc_t(),
376  opts.low_arg, opts.upper_arg );
377  break;
378  case op_arg_notch:
379  bp_notch( data, nr, ftp, iftp, out, orig.t(), orig.inc_t(),
380  opts.low_arg, opts.notch_width_arg );
381  break;
382  case op_arg_Butter_bp :
383  filter_all( data, nr, filt, outf, orig.t() );
384  break;
385  }
386  write_traces( data, nr, n, &dfthead );
387  }
388  gzclose( in );
389 
390  fclose( dftout );
391 
392  // overwrite original file
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" );
397  int c;
398  while( (c=getc(dftout)) != EOF )
399  putchar(c);
400  fclose( dftout );
401  //remove( out_fname.c_str() );
402  }
403 
404  return 0;
405 }
constexpr T min(T a, T b)
Definition: ion_type.h:33
int type(void)
Definition: IGBheader.h:280
int read_data(T *dp, size_t numt=1, char *buf=NULL)
Definition: IGBheader.h:431
void bp_notch(double **in, int n, fftw_plan fftp, fftw_plan bftp, fftw_complex **out, int t, float dt, float notch, int nwidth)
Definition: IGBdft.cc:125
#define IGB_FLOAT
Definition: IGBheader.h:60
size_t slice_sz()
Definition: IGBheader.h:263
int main(int argc, char *argv[])
Definition: IGBdft.cc:240
size_t t(void)
Definition: IGBheader.h:277
size_t data_size(void)
Definition: IGBheader.h:267
int write()
Definition: IGBheader.cc:343
void unites(const char *a)
Definition: IGBheader.h:365
void filter_all(double **in, int nr, IIR_Filt *filt, double *out, int t)
Definition: Filter.cc:76
void bp_filt(double **in, int n, fftw_plan fftp, fftw_plan bftp, fftw_complex **out, int t, float dt, float low, float high)
Definition: IGBdft.cc:85
float org_t(void)
Definition: IGBheader.h:317
void unites_t(const char *a)
Definition: IGBheader.h:362
void do_domF(double **in, int n, fftw_plan ftp, fftw_complex **out, float t, float dt)
Definition: IGBdft.cc:187
void inc_t(float a)
Definition: IGBheader.h:329
#define creal(A)
Definition: IGBdft.cc:34
void do_PSD(double **in, int n, fftw_plan ftp, fftw_complex **out, int t, float dt)
Definition: IGBdft.cc:158
#define cimag(A)
Definition: IGBdft.cc:35
void write_traces(T **data, int nr, int n, IGBheader *ohead)
Definition: IGBdft.cc:226
void fileptr(gzFile f)
Definition: IGBheader.cc:329
void do_phase(double **in, int n, fftw_plan fftp, fftw_plan bftp, fftw_complex **out, int t)
Definition: IGBdft.cc:50