IGBdft.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.
// ----------------------------------------------------------------------------

#include <fftw3.h>
#include <math.h>
#include <string>
#include <iostream>
#include <algorithm>
#include "IGBheader.h"
#include "dft_cmdline.h"
#include "Filter.h"
#ifdef _OPENMP
#include <omp.h>
#endif

using namespace std;

#define creal(A)  (A[0])
#define cimag(A)  (A[1])
#define cmag(A)   sqrt(A[0]*A[0]+A[1]*A[1])

void
do_phase( double **in, int n, fftw_plan fftp, fftw_plan bftp,
                                   fftw_complex **out, int t )
{
#pragma omp parallel for
  for( int i=0; i<n; i++ ) {
#ifdef _OPENMP
    int thr = omp_get_thread_num();
#else
    int thr = 0;
#endif
    memset( out[thr], 0, sizeof(fftw_complex)*t );
    fftw_execute_dft_r2c( fftp, in[i], out[thr] );
    creal(out[thr][0]) = cimag(out[thr][0]) = 0.; // remove dc
    // this is an inverse DFT which ignores the negative components
    fftw_execute_dft( bftp, out[thr], out[thr] );
    for( int j=0; j<t; j++ )
      in[i][j] = atan2( cimag(out[thr][j]), creal(out[thr][j]) );
  }
}

void
bp_filt( double **in, int n, fftw_plan fftp, fftw_plan bftp,
               fftw_complex **out, int t, float dt, float low, float high )
{
  int mini = low*dt*(t/2)/500.;
  int maxi = high ? high*dt*(t/2)/500. : t;

#pragma omp parallel for
  for( int i=0; i<n; i++ ) {
#ifdef _OPENMP
    int thr = omp_get_thread_num();
#else
    int thr = 0;
#endif
    fftw_execute_dft_r2c( fftp, in[i], out[thr] );
    for( int j=0; j<=mini; j++ ) 
      creal(out[thr][j]) = cimag(out[thr][j]) = 0.;
    for( int j=maxi; j<=t/2; j++ ) 
      creal(out[thr][j]) = cimag(out[thr][j]) = 0.;
    fftw_execute_dft_c2r( bftp, out[thr], in[i] );
    // doing the transform and its inverse scales by t
    for( int j=0; j<t; j++ )
      in[i][j] /= t;
  }
}

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 notchi = notch*dt*(t/2)/500.;

#pragma omp parallel for
  for( int i=0; i<n; i++ ) {
#ifdef _OPENMP
    int thr = omp_get_thread_num();
#else
    int thr = 0;
#endif
    fftw_execute_dft_r2c( fftp, in[i], out[thr] );
    for( int j=notchi-nwidth; j<=notchi+nwidth; j++ ) 
      creal(out[thr][j]) = cimag(out[thr][j]) = 0.;
    fftw_execute_dft_c2r( bftp, out[thr], in[i] );
    // doing the transform and its inverse scales by t
    for( int j=0; j<t; j++ )
      in[i][j] /= t;
  }
}


void
do_PSD( double **in, int n, fftw_plan ftp, 
        fftw_complex **out, int t, float dt )
{
#pragma omp parallel for
  for( int i=0; i<n; i++ ) {
#ifdef _OPENMP
    int thr = omp_get_thread_num();
#else
    int thr = 0;
#endif
    fftw_execute_dft_r2c( ftp, in[i], out[thr] );
    for( int j=0; j<=t/2; j++ ) {
      in[i][j] = creal(out[thr][j])*creal(out[thr][j])+
                     cimag(out[thr][j])*cimag(out[thr][j]);
    }
  }
}


void
do_domF( double **in, int n, fftw_plan ftp, 
        fftw_complex **out, float t, float dt )
{
#pragma omp parallel for
  for( int i=0; i<n; i++ ) {
#ifdef _OPENMP
    int thr = omp_get_thread_num();
#else
    int thr = 0;
#endif

    fftw_execute_dft_r2c( ftp, in[i], out[thr] );

    int maxf = 1;
    float maxmag = creal(out[thr][1])*creal(out[thr][1])+
                   cimag(out[thr][1])*cimag(out[thr][1]);
    for( int j=2; j<t/2; j++ ) {
      float tmag = creal(out[thr][j])*creal(out[thr][j])+
                   cimag(out[thr][j])*cimag(out[thr][j]);
      if( tmag > maxmag ) {
        maxmag = tmag;
        maxf = j;
      }
    }
    in[i][0] = 500./dt/(t/2.)*maxf;
  }
}


template<class T>
void
write_traces( T **data, int nr, int n, IGBheader *ohead )
{
  FILE *fp = static_cast<FILE *>(ohead->fileptr());

  for( long t=0; t<ohead->t(); t++ ) {
    fseek( fp, 1024+(t*ohead->slice_sz()+n)*ohead->data_size(), SEEK_SET );
    for( int i=0; i<nr; i++ ){
      float d = data[i][t];
      fwrite( &d, sizeof(float), 1, fp );
    }
  }
}


int main(int argc, char *argv[] )
{
#ifdef _OPENMP
  if(!getenv("OMP_NUM_THREADS") ){
    cerr << "unst" << endl;
    omp_set_num_threads(1);
  }
#endif

  gengetopt_args_info opts;

  if (cmdline_parser(argc, argv, &opts) != 0)
    exit(1);

  if( !opts.inputs_num ) {
    cmdline_parser_print_help();
    exit(0);
  }

  gzFile in = gzopen((!opts.inputs_num||!strcmp(opts.inputs[0],"-"))?"/dev/stdin":opts.inputs[0],"r");
  if( !in ) {
    cerr << "Cannot open input file: " << opts.inputs[0] << endl;
    exit(1);
  }
  IGBheader orig(in,true);

  FILE* dftout;
  char tname[] = "/tmp/igbdftXXXXXX";
  string out_fname;
  if( !strcmp( opts.output_arg, "." ) ||
      !strcmp( opts.output_arg, "-") ) {
    mkstemp( tname );
    out_fname = tname;
  } else
    out_fname = opts.output_arg;

  dftout = fopen( out_fname.c_str(), "w" );
  if( !dftout ) {
    cerr << "Cannot open output file: " << out_fname << endl;
    exit(1);
  }

  IGBheader dfthead = orig;
  dfthead.fileptr( dftout );
  dfthead.type(IGB_FLOAT);

  int      npar         = min(opts.numpar_arg,(int)orig.slice_sz());
  double **data         = new double *[npar];
  int      aligned_size =((orig.t()+7)/8)*8;   // align for SIMD
  double  *bbuf         = fftw_alloc_real(npar*aligned_size);
  for( int i=0; i<npar; i++ ) 
    data[i] = bbuf+aligned_size*i;

  fftw_plan     ftp, iftp;
  IIR_Filt     *filt;
  double       *outf;

#ifdef _OPENMP
  int numthr=omp_get_max_threads();
#else
  int numthr=1;
#endif
  fftw_complex **out = new fftw_complex*[numthr];

  switch( opts.op_arg ) {
      case op_arg_DFT_bp :  
      case op_arg_notch :  
        for( int i=0; i<numthr; i++ )
          out[i] = fftw_alloc_complex(orig.t()/2+1);
        if( !(ftp = fftw_plan_dft_r2c_1d( orig.t(), data[0], out[0], 0)) ) {
          cerr << "Cannot create fftw3 plan!" << endl;
          exit(1);
        }
        iftp = fftw_plan_dft_c2r_1d( orig.t(), out[0], data[0], 0);
        break;
      case op_arg_PSD :  
        dfthead.unites_t("Hz");
        dfthead.t(orig.t()/2+1);
        dfthead.org_t(0);
        dfthead.inc_t(1000./orig.t()/orig.inc_t()); 
        for( int i=0; i<numthr; i++ )
          out[i] = fftw_alloc_complex(orig.t()/2+1);
        if( !(ftp = fftw_plan_dft_r2c_1d( orig.t(), data[0], out[0], 0)) ) {
          cerr << "Cannot create fftw3 plan!" << endl;
          exit(1);
        }
        break;
      case op_arg_DomF :
        dfthead.unites("Hz");
        dfthead.t(1);
        for( int i=0; i<numthr; i++ )
          out[i] = fftw_alloc_complex(orig.t()/2+1);
        if( !(ftp = fftw_plan_dft_r2c_1d( orig.t(), data[0], out[0], 0)) ) {
          cerr << "Cannot create fftw3 plan!" << endl;
          exit(1);
        }
        break;
      case op_arg_phase :
        dfthead.unites("radians");
        for( int i=0; i<numthr; i++ )
          out[i]  = fftw_alloc_complex(orig.t());
        ftp  = fftw_plan_dft_r2c_1d( orig.t(), data[0], out[0], 0 );
        iftp = fftw_plan_dft_1d( orig.t(), out[0], out[0], FFTW_BACKWARD, 0);
        break;
      case op_arg_Butter_bp :
        filt = new Butterworth(  opts.order_arg, opts.low_arg*orig.inc_t()/500.,
                             opts.upper_arg*orig.inc_t()/500. );
        outf = new double[orig.t()];
        break;
  }
  dfthead.write();

  double *slice    = new double[orig.slice_sz()];
  float  *dftslice = new float[orig.slice_sz()];

  for( int n=0; n<orig.slice_sz(); n+=npar ) {
    int nr = orig.slice_sz()-n>=npar?npar:orig.slice_sz()-n;
    cerr << "Reading nodes ["<< n << "-" << (n+nr-1) << "]"<< endl; 
    gzseek( in, 1024, SEEK_SET );
    for( int i=0; i<orig.t(); i++ ) {
      orig.read_data( slice );
      for( int j=0; j<nr; j++ )
        data[j][i] = slice[n+j];
    }
    switch( opts.op_arg ) {
        case op_arg_PSD :
          do_PSD( data, nr, ftp, out, orig.t(), orig.inc_t() );
          break;
        case op_arg_DomF :
          do_domF( data, nr, ftp, out, orig.t(), orig.inc_t() );
          break;
        case op_arg_phase :
          do_phase( data, nr, ftp, iftp, out, orig.t() );
          break;
        case op_arg_DFT_bp:
          bp_filt( data, nr, ftp, iftp, out, orig.t(), orig.inc_t(), 
                   opts.low_arg, opts.upper_arg );
          break;
        case op_arg_notch:
          bp_notch( data, nr, ftp, iftp, out, orig.t(), orig.inc_t(), 
                   opts.low_arg, opts.notch_width_arg );
          break;
        case op_arg_Butter_bp :
          filter_all( data, nr, filt, outf, orig.t() );
          break;
    }
    write_traces( data, nr, n, &dfthead ); 
  }
  gzclose( in );

  fclose( dftout );

  // overwrite original file
  if( !strcmp( opts.output_arg, "." ) )
    rename( out_fname.c_str(),  opts.inputs[0] );
  else if( !strcmp( opts.output_arg, "-" ) ) {
    dftout = fopen( out_fname.c_str(), "r" );
    int c;
    while( (c=getc(dftout)) != EOF )
      putchar(c);
    fclose( dftout );
    //remove( out_fname.c_str() );
  }

  return 0;
}