source: Sophya/trunk/ArchTOIPipe/ProcWSophya/simtoipr.cc@ 1546

#include "array.h"
#include "simtoipr.h"
#include <math.h>
#include "toimanager.h"
#include "pexceptions.h"
#include "ctimer.h"
#include "fftpserver.h"

SimpleDeglitcher::SimpleDeglitcher(int wsz, double ns, int maxnpt, int minnpt)
{
  SetWSize(wsz);
  SetDetectionParam(ns, ns/2., maxnpt, minnpt);
  SetRange(-9.e39, 9.e39);
  RepBadSamples(true, true);

  totnscount = glnscount = glcount = out_range_nscount = 0;
  deglitchdone = false;

  cout << "SimpleDeglitcher::SimpleDeglitcher() WSize= " << wsz 
       << " nSig=" << ns << " maxNPt=" << maxnpt << endl;
}

SimpleDeglitcher::~SimpleDeglitcher()
{
}

void SimpleDeglitcher::SetDetectionParam(double ns, double ns2, int maxnpt, 
                                         int minnpt, int wszrec)
{
  nsig = (ns > 0.01) ? ns : 1.;
  nsig2 = (ns2 > 0.01) ? ns2 : nsig/2.;
  maxpoints = ((maxnpt > 0) && (maxnpt <= wsize)) ? maxnpt : 5;
  minpoints = ((minnpt > 0) && (minnpt <= maxpoints)) ? minnpt : maxpoints/2;
  wrecsize = ((wszrec > 0) && (wszrec <= wsize)) ? wszrec : 2*maxpoints;
}

inline char * _Bool2YesNo(bool fg) 
{
  if (fg) return "YES" ;
  else return "NO" ;
}

void SimpleDeglitcher::PrintStatus(ostream & os)
{
  os << "\n ------------------------------------------------------ \n" 
     << " SimpleDeglitcher::PrintStatus() - WindowSize=" << WSize() 
     << " MaxPoints=" << MaxPoints() << " MinPoints=" << MinPoints() << endl;
  os << " NbSigmas=" << NbSigmas() 
     << "  NbSigmas2=" << NbSigmas2() << "  WRecSize= " << WRecSize() 
     << "  Range_Min= " << range_min << " Range_Max= " << range_max << endl;
  os << "  RepOutOfRangeSamples: " << _Bool2YesNo(rec_out_range_samples) 
     << "  RepGlitchSamples: " << _Bool2YesNo(rec_gl_samples) 
     << "  UseWRec: " << _Bool2YesNo(rec_use_wrec) << endl;
  TOIProcessor::PrintStatus(os);
  if (deglitchdone) os << " Deglitching performed " << endl;
  else os << " NO deglitching done " << endl;
  double nst = (ProcessedSampleCount() > 0) ? ProcessedSampleCount() : 1.;
  os << " ProcessedSampleCount=" << ProcessedSampleCount() 
     << " OutOfRangeSampleCount=" << OutOfRangeSampleCount() << endl;
  os << " GlitchCount= " << GlitchCount() 
     << " GlitchSampleCount=" << GlitchSampleCount() 
     << "( " << (double)GlitchSampleCount()*100./nst << " % )" << endl;
  os << " ------------------------------------------------------ " << endl;
}

void SimpleDeglitcher::init() {
  cout << "SimpleDeglitcher::init" << endl;
  declareInput("in");
  declareOutput("out");
  declareOutput("mean");
  declareOutput("sigma");
  declareOutput("incopie");
  name = "SimpleDeglitcher";
  //  upExtra = 1;   A quoi ca sert ?
}


#define FG_OUTOFRANGE 1
#define FG_GLITCH 2
void SimpleDeglitcher::run() {

  //  TOIManager* mgr = TOIManager::getManager();
  int snb = getMinIn();
  int sne = getMaxIn();

  bool fgout = checkOutputTOIIndex(0);
  bool fgmean = checkOutputTOIIndex(1);
  bool fgsigma = checkOutputTOIIndex(2);
  bool fgincopie = checkOutputTOIIndex(3);
  
  if (!checkInputTOIIndex(0)) {
    cerr << " SimpleDeglitcher::run() - Input TOI (in) not connected! "
         << endl;
    throw ParmError("SimpleDeglitcher::run() Input TOI (in) not connected!");
  }
  if (!fgout && !fgmean && !fgsigma &&!fgincopie) {
    cerr << " SimpleDeglitcher::run() - No Output TOI connected! "
         << endl;
    throw ParmError("SimpleDeglitcher::run() No output TOI connected!");
  }
  
  if (!fgout) {
    cout << "Warning:  SimpleDeglitcher::run() - No TOI connected to out \n"
         << " No deglitching would be performed !" << endl;
  }

  cout << " SimpleDeglitcher::run() SNRange=" << snb << " - " << sne << endl; 
  try {
    Timer tm("SimpleDeglitcher::run()");
    Vector vin(wsize);
    Vector vas(wsize);

    int wrecsize = maxpoints*2;
    Vector vrec(wrecsize);

    TVector<uint_8> vfg(wsize);
    int wsz2 = wsize/2;
    // Le debut 
    int k;
    for(k=0; k<wsz2; k++) 
      getData(0, k+snb, vin(k), vfg(k));

    int nokdebut = 0;
    double s = 0.;    
    double mean = 0.;
    double s2 = 0.;
    double sigma = 0.;
    for(k=0; k<wsz2; k++) {
      if ( vfg(k) != 0) continue;
      if ( (vin(k) < range_min) || (vin(k) > range_max) ) continue;
      s += vin(k);
      s2 += vin(k)*vin(k);
      nokdebut++;
    }
    if (nokdebut > 0) {
      mean = s/nokdebut;
      if (nokdebut > 1)  sigma = sqrt(s2/nokdebut-mean*mean);
    }
    for(k=wsz2; k<wsize; k++) {
      vin(k) = mean;
      vfg(k) = 0;
    } 
    for(k=0; k<wsize; k++) {
      vas(k) = mean;
      if ( vfg(k) != 0) continue;
      if ( (vin(k) < range_min) || (vin(k) > range_max) ) continue;
      vas(k) = vin(k);
    }

    for(k=0; k<wrecsize; k++) {
      if ( (vin(k) < range_min) || (vin(k) > range_max) ) vrec(k)=mean;
      else vrec(k)=vin(k);    
    }

    bool fgokdebut = false;

    int kgl = -1;
    int ii,lastput;
    bool fgglitch = false;
    double valcur,valsub,valadd;
    double lastvalok = mean;
    uint_8 fgcur;
    bool fgokcur=false;

    int sx_refresh_count = 0;
    int sx_refresh_count_max = 16*wsize;

    // Boucle sur les sampleNum 
    int knext;
    int kfin = sne-snb;
    for(k=0;k<=kfin;k++) {
      totnscount++;
//       if (k%10000 == 0) cout << " DBG: K=" << k << endl;
      knext = k+wsz2;
      // Calcul mean-sigma
      if (knext<=kfin)  {
        valsub = vas(knext%wsize);
        getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));
        valadd = vin(knext%wsize);
        if ( vfg(knext%wsize) || 
             (valadd < range_min) || (valadd > range_max) ) {
          vas(knext%wsize) = valadd = mean;
          fgokcur = false;
        }
        else { 
          vas(knext%wsize) = valadd = vin(knext%wsize);
          fgokcur = true; 
        }
        if (!fgokdebut && fgokcur) {
          s += valadd;
          s2 += valadd*valadd;
          nokdebut++;
          mean = s/nokdebut;
          if (nokdebut > 1)  sigma = sqrt(s2/nokdebut-mean*mean);
          if (nokdebut >= wsize)  { 
            fgokdebut = true;
            cout << " SimpleDeglitcher::DebugInfo -  nokdebut=" << nokdebut 
                 << " k=" << k << " knext=" << knext 
                 << "\n ...DebugInfo mean=" << mean 
                 << " sigma=" << sigma << " s=" << s << " s2=" << s2 << endl;
          }
        } 
        else {
          if (sx_refresh_count >= sx_refresh_count_max) {
            //  On recalcule la somme
            s = vas.Sum();
            s2 = vas.SumX2();
            sx_refresh_count = 0;
          }
          else {
            s += (valadd-valsub);
            s2 += (valadd*valadd-valsub*valsub);
            sx_refresh_count++;
          }
          mean = s/wsize;
          sigma = sqrt(s2/wsize-mean*mean);
        }
      }


      //  On gere les sorties Mean et Sigma
      if (fgmean) 
        putData(1, k+snb, mean, 0);
      if (fgsigma) 
        putData(2, k+snb, sigma, 0);
      if (fgincopie) 
        putData(3, k+snb, vin(k%wsize), vfg(k%wsize));

      valcur = vin(k%wsize);
      if ( (valcur < range_min) || (valcur > range_max) ) {
        valcur = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
        if (rec_out_range_samples) vin(k%wsize) = valcur;
        vfg(k%wsize) |= FG_OUTOFRANGE;
        out_range_nscount++;
      }
      valcur = vin(k%wsize);
      fgcur = vfg(k%wsize);
        
      if (!fgout)   continue;  // Pas de sortie out (deglitche)

      
//             if (k<100) { 
//              cout << "DBG-A-Deglitch[" << k << "] mean=" 
//                   << mean << " sigma=" << sigma << " valcur=" 
//                   << valcur << " Kgl=" << kgl ;
//              if (fgglitch) cout << " In Glitch" ;
//              cout << endl;
//            }

      double curnsig = (fgglitch) ? nsig2 : nsig;
        
      if (valcur < mean+curnsig*sigma) {  // inferieur au seuil
        if (fgglitch) {  
          if ( (k-kgl <= maxpoints) && (k-kgl >= minpoints) ) {   
            // On vient de detecter un glitch
            glcount++;
            if (rec_gl_samples) { // On change la valeur des samples
              double recval = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
              for(ii=kgl; ii<k; ii++) {
                putData(0, ii+snb, recval, vfg(ii%wsize)|FG_GLITCH);
                glnscount++; 
              }
            }
            else {   // On ne fait que flagger les echantillons
              for(ii=kgl; ii<k; ii++) {
                putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize)|FG_GLITCH);
                glnscount++; 
              }
            }
            lastput = snb+k-1;
          }  // - Fin de detection de glitch
          else {  // Trop long ou trop court - ce n'est pas un glitch ...
            for(ii=kgl; ii<k; ii++) {
              putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
            }
            lastput = snb+k-1;  
          } 
        }
        putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
        lastput = snb+k;
        kgl = -1;  fgglitch = false;
        vrec(k%wrecsize) = lastvalok = valcur;
      }
      else {  // Superieur au seuil
        if (fgglitch) {
          if (k-kgl+1 > maxpoints) {  // serie de points > seuil
            for(ii=kgl; ii<=k; ii++)   // -> Donc pas glitch
              putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
            lastput = snb+k;
            fgglitch = false;  
            vrec(k%wrecsize) = lastvalok = valcur;          
          }
        }
        else {
          if (kgl < 0) {  // debut possible de glitch
            fgglitch = true;  kgl = k;
          }
          else { // On est toujours dans une serie > seuil
            putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
            lastput = snb+k;   lastvalok = valcur;
          }
          vrec(k%wrecsize) = lastvalok;     
        }
      }

//           if (k%5000 == 0) cout << " ---DBG2: K=" << k << " glcount=" 
//                              << glcount << " LastPut= " << lastput << endl;
    }  // Fin de Boucle sur les num-sample

    //DBG    cout << " La fin lastput=" << lastput << " SNE=" << sne; 
    //DBG    for(k=lastput-snb+1; k<sne-snb; k++) 
    //DBG      putData(0, k+snb, vin(k%wsize), 0);
    //    cout << " DBG3- OUT of try bloc ! " << endl;
    cout << " SimpleDeglitcher::run() - End of processing " 
    << " ProcessedSampleCount=" << ProcessedSampleCount() 
    <<  " GlitchCount= " << GlitchCount() << endl;
    if (fgout) deglitchdone = true;
  }  // Bloc try 
  catch (PException & exc) {
    cerr << "SimpleDeglitcher: Catched Exception " << (string)typeid(exc).name()
         << "\n .... Msg= " << exc.Msg() << endl;
  }                                                                             
}


// -------------------------------------------------------------------
//   Classe SimpleFilter : Filtre simple ds le domaine temporel
// -------------------------------------------------------------------

string SimpleFilter::FilterKind2String(FilterKind fk)
{
  switch (fk) {
    case UserFilter :
      return ("UserFilter");
      break;
    case MeanFilter :
      return ("MeanFilter");
      break;
    case SumFilter :
      return ("SumFilter");
      break;
    case GaussFilter :
      return ("GaussFilter");
      break;
    case DiffFilter :
      return ("DiffFilter");
      break;
    default :
      return ("ErrorFilterKind");
      break;
  }
  return("");
}

SimpleFilter::SimpleFilter(int wsz, FilterKind fk, double a, double s)
{
  if (wsz < 3) wsz = 3;
  if (wsz%2 == 0) wsz++;
  cout << "SimpleFilter::SimpleFilter() wsz= " << wsz 
       << " FilterKind=" << FilterKind2String(fk) << endl; 
  wsize = wsz;
  totnscount = 0;
  coef = new double[wsz];
  for(int k=0; k<wsz; k++) coef[k] = 0.;
  int kk;
  switch (fk) {
    case UserFilter :
      throw ParmError("SimpleFilter: Error in filter Kind (UserFilter)!");
      //      break;
    case MeanFilter :
      for(kk=0; kk<wsz; kk++) 
        coef[kk] = a/wsize;
      break;
    case SumFilter :
      for(kk=0; kk<wsz; kk++) 
        coef[kk] = a;
      break;
    case GaussFilter :
      for(kk=-(wsz/2); kk<=(wsz/2); kk++) 
        coef[kk+(wsz/2)] = a*exp(-(double)(kk*kk)/(2*s*s));
      break;
    case DiffFilter :
      for(kk=0; kk<wsz; kk++) 
        coef[kk] = -a/wsize;
      coef[wsz/2+1] += 1.;
      break;
    default :
      throw ParmError("SimpleFilter: Error in filter Kind (UnknownFilter)!");
      //      break;
  }
}

SimpleFilter::SimpleFilter(Vector const & vc)
{
  int wsz = vc.Size();
  if (wsz < 3) wsz = 3;
  if (wsz%2 == 0) wsz++;
  FilterKind fk = UserFilter;
  cout << "SimpleFilter::SimpleFilter(Vector & vc) vc.Size()= " 
       << vc.Size() << " WSize=" << wsz 
       << " FilterKind=" << FilterKind2String(fk) << endl; 
  wsize = wsz;
  totnscount = 0;
  coef = new double[wsz];
  int kk;
  for(kk=0; kk<vc.Size(); kk++) 
    coef[kk] = vc(kk);
  for(kk=vc.Size(); kk<wsz; kk++) 
    coef[kk] = 0.;
}

SimpleFilter::~SimpleFilter()
{
  delete[] coef;
}

void SimpleFilter::PrintStatus(ostream & os)
{
  os << "\n ------------------------------------------------------ \n" 
     << " SimpleFilter::PrintStatus() - WindowSize=" << WSize() 
     << " FilterKind= " << Type() << endl;
  TOIProcessor::PrintStatus(os);
  os << " Coeff= " ; 
  for(int k=0; k<wsize; k++) os << coef[k] << " " ;
  os << endl;
  os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
  os << " ------------------------------------------------------ " << endl;
}

void SimpleFilter::init() {
  cout << "SimpleFilter::init" << endl;
  declareInput("in");
  declareOutput("out");
  declareOutput("incopie");
  name = "SimpleFilter";
  //  upExtra = 1;
}

void SimpleFilter::run() {
  //  TOIManager* mgr = TOIManager::getManager();
  int snb = getMinIn();
  int sne = getMaxIn();

  bool fgout = checkOutputTOIIndex(0);
  bool fgincopie = checkOutputTOIIndex(1);

  if (!checkInputTOIIndex(0)) {
    cerr << " SimpleFilter::run() - Input TOI (in) not connected! "
         << endl;
    throw ParmError("SimpleFilter::run() Input TOI (in) not connected!");
  }
  if (!fgout) {
    cerr << " SimpleFilter::run() - No Output TOI connected! "
         << endl;
    throw ParmError("SimpleFilter::run() No output TOI connected!");
  }

  cout << " SimpleFilter::run() SNRange=" << snb << " - " << sne << endl; 
  

  try {
    Timer tm("SimpleFilter::run()");
    // Le debut
    int wsz2 = wsize/2;
    Vector vin(wsize);
    TVector<uint_8> vfg(wsize);
    int k;
    for(k=0; k<wsize; k++) 
      getData(0, k+snb, vin(k%wsize), vfg(k%wsize));
    
    double mean = vin.Sum()/wsize;
    for(k=wsz2+1; k<wsize; k++) { 
      vin(k) = mean; 
      vfg(k) = 0;
    }
    int knext;
    bool fgfin = false;
    // Boucle sur les sampleNum
    for(k=0;k<=sne-snb;k++) {
      double sortie = 0;
      for(int ii=-wsz2; ii<=wsz2; ii++) {
        sortie += vin((ii+k+wsize)%wsize)*coef[ii+wsz2];
      }
      putData(0,k+snb,sortie,vfg(k%wsize));
      if (fgincopie) 
        putData(1, k+snb, vin(k%wsize), vfg(k%wsize));
      knext = k+wsz2+1;
      if (knext<=(sne-snb)) 
        getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));

      else {
        if (!fgfin) { 
          mean = vin.Sum()/wsize;
          fgfin = true;
        }
        vin(knext%wsize) = mean;
        vfg(knext%wsize) = 0;
      }
      totnscount++;
    }  // Boucle sur les num-sample
    cout << " SimpleFilter::run() - End of processing " << endl;
  }  // Bloc try 

  catch (PException & exc) {
    cerr << "SimpleFilter: Catched Exception " << (string)typeid(exc).name()
         << "\n .... Msg= " << exc.Msg() << endl;
  }                                                                             
}

// ---------------------------------------------------------------
// -------------------- Classe SimpleAdder -----------------------
// ---------------------------------------------------------------

SimpleAdder::SimpleAdder(int nbinput)
  : gains(nbinput)
{
  if (nbinput < 1) 
    throw ParmError("SimpleAdder::SimpleAdder() NbInput < 1 !");
  nb_input = nbinput;
  for(int k=0; k<nb_input; k++) gains(k) = 1.;
  totnscount = 0;
}

SimpleAdder::~SimpleAdder()
{
}

void SimpleAdder::SetGain(int num, double g)
{
  if ((num < 0) || (num >= nb_input))
    throw RangeCheckError("SimpleAdder::SetGain() Out of range input number!");
  gains(num) = g;
  return;
}

double SimpleAdder::Gain(int num)
{
  if ((num < 0) || (num >= nb_input))
    throw RangeCheckError("SimpleAdder::Gain() Out of range input number!");
  return gains(num);
}

void SimpleAdder::PrintStatus(ostream & os)
{
  os << "\n ------------------------------------------------------ \n" 
     << " SimpleAdder::PrintStatus() - NbInput=" << NbInput() << endl; 
  TOIProcessor::PrintStatus(os);
  os << " Gains= " ; 
  for(int k=0; k<nb_input; k++) os << gains(k) << " " ;
  os << endl;
  os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
  os << " ------------------------------------------------------ " << endl;
}

void SimpleAdder::init() {
  cout << "SimpleAdder::init NbInput=" << nb_input << endl;
  char buff[32];
  for(int k=0; k<nb_input; k++) {
    sprintf(buff,"in%d", k);
    declareInput(buff);
  }

  declareOutput("out");
  name = "SimpleAdder";
  //  upExtra = 1;
}

void SimpleAdder::run() {
  //  TOIManager* mgr = TOIManager::getManager();
  int snb = getMinIn();
  int sne = getMaxIn();

  bool fgout = checkOutputTOIIndex(0);
  string msg_err;
  for(int ki=0;ki<nb_input;ki++) {
    if (!checkInputTOIIndex(ki)) {
      msg_err = "SimpleAdder::run() - Input TOI (" + getInName(ki) +
                " not connected!";
      cerr << msg_err << endl;
      throw ParmError(msg_err);
    }
  }
  if (!fgout) {
    cerr << " SimpleAdder::run() - No Output TOI connected! "
         << endl;
    throw ParmError("SimpleAdder::run() No output TOI connected!");
  }

  cout << " SimpleAdder::run() SNRange=" << snb << " - " << sne << endl; 
  

  try {
    Timer tm("SimpleAdder::run()");
    int k,i;
    double out = 0.;
    double valin = 0.;
    uint_8 fgin = 0;
    uint_8 fgout = 0;
    for(k=snb;k<=sne;k++) {
      out = 0;
      fgout = 0;
      for(i=0; i<nb_input; i++) {
        getData(i, k, valin, fgin);
        out += gains(i)*valin;
        fgout = fgout | fgin;
      }
    putData(0,k,out,fgout);
    totnscount++;
    }  // Boucle sur les num-sample
    cout << " SimpleAdder::run() - End of processing " << endl;
  }  // Bloc try 

  catch (PException & exc) {
    cerr << "SimpleAdder: Catched Exception " << (string)typeid(exc).name()
         << "\n .... Msg= " << exc.Msg() << endl;
  }                                                                             
}

// ----------------------------------------------------------------------
//   Classe SimpleFourierFilter : Filtre simple ds le domaine de Fourier
// ----------------------------------------------------------------------

SimpleFourierFilter::SimpleFourierFilter(Vector const & vc)
{
  ffcoef = vc;
  wsize = (ffcoef.Size()-1)*2;
  if (wsize < 16) 
    throw ParmError("SimpleFourierFilter::SimpleFourierFilter() WSize<16 !");
  KeepSpectra("spectra.ppf", 0);
  ComputeMeanSpectra(false);
  totnscount = 0;
  totnbblock = 0;
}

SimpleFourierFilter::~SimpleFourierFilter()
{
}


void SimpleFourierFilter::PrintStatus(ostream & os)
{
  os << "\n ------------------------------------------------------ \n" 
     << " SimpleFourierFilter::PrintStatus() - WindowSize=" 
     << WSize() << endl; 
  TOIProcessor::PrintStatus(os);
  os << " Coeff (Size= " << ffcoef.Size() << "): " << endl;  
  for(int i=0; i<16; i++) {
    os << ffcoef(i) << "  " ; 
    if (i == 7) os << endl;
  }
  os << "  .... " << endl;
  os << " ProcessedSampleCount=" << ProcessedSampleCount() 
     << " NbFFTBlocks= " << totnbblock << endl;
  os << " ------------------------------------------------------ " << endl;
}

void SimpleFourierFilter::init() {
  cout << "SimpleFourierFFilter::init" << endl;
  declareInput("in");
  declareOutput("out");
  declareOutput("incopie");
  name = "SimpleFourierFilter";
  //  upExtra = 1;
}


void SimpleFourierFilter::run() {
  //  TOIManager* mgr = TOIManager::getManager();
  int snb = getMinIn();
  int sne = getMaxIn();

  bool fgout = checkOutputTOIIndex(0);
  bool fgincopie = checkOutputTOIIndex(1);

  if (!checkInputTOIIndex(0)) {
    cerr << " SimpleFourierFilter::run() - Input TOI (in) not connected! "
         << endl;
    throw ParmError("SimpleFourierFilter::run() Input TOI (in) not connected!");
  }

  // ---- On peut utiliser cette classe pour calculer un spectre de Fourier ----
  //  if (!fgout) {
  //    cerr << " SimpleFourierFilter::run() - No Output TOI connected! "
  //     << endl;
  //    throw ParmError("SimpleFourierFilter::run() No output TOI connected!");
  //  }

  cout << " SimpleFourierFilter::run() SNRange=" << snb << " - " << sne << endl; 
  

  try {
    Timer tm("SimpleFourierFilter::run()");
    // Le debut
    Vector vin(wsize);
    Vector vout(wsize);
    TVector<uint_8> vfg(wsize);
    TVector< complex<r_8> > vfft, vfftmean;
    Vector meanpowerspectra;
    TVector< complex<r_8> > zcoef(ffcoef.Size());
    zcoef = ffcoef;


    FFTPackServer ffts;
    ffts.setNormalize(true);

    POutPersist pout(outppfname);

    int k,i,klast;
    int nks = 0;
    klast = snb-1;
    totnbblock = 0;
    // Boucle sur les sampleNum
    // 1er partie, on traite par paquets de wsize
    for(k=snb;k<=sne-wsize+1;k+=wsize) {
      for(i=0; i<wsize; i++) 
        getData(0, k+i, vin(i), vfg(i));
      ffts.FFTForward(vin, vfft);
      if (c_meanspectra) { // Compute mean-spectra
        if (totnbblock == 0)  { 
          vfftmean = vfft;
          meanpowerspectra.ReSize(vfft.Size());
          for(i=0; i<meanpowerspectra.Size(); i++) 
            meanpowerspectra(i) = sqrt(vfft(i).real()*vfft(i).real() + 
                                    vfft(i).imag()*vfft(i).imag() );
        }
        else { 
          vfftmean += vfft;
          for(i=0; i<meanpowerspectra.Size(); i++) 
            meanpowerspectra(i) += sqrt(vfft(i).real()*vfft(i).real() + 
                                        vfft(i).imag()*vfft(i).imag() );
        }
      }
      totnbblock++;
      if (nks < nb_keep) {
        TVector< complex<r_8> > vfftcopie;
        vfftcopie = vfft;
        string nomvfft = "spectra" + (string)MuTyV(nks);
        pout.PutObject(vfftcopie, nomvfft); 
        nks++;
      }
      if (fgout) { 
        vfft.MulElt(zcoef);
        ffts.FFTBackward(vfft, vout);
      }
      for(i=0; i<wsize; i++) {
        if (fgout) 
          putData(0,k+i,vout(i),vfg(i));
        if (fgincopie)  
          putData(1, k+i, vin(i), vfg(i));
      }
      klast+=wsize;
      totnscount+=wsize;
    }
    

    // 2eme partie, on traite la fin du bloc d'echantillons si necessaire
    double inval;
    uint_8 inflg;
    if (klast < sne) 
      for(k=klast+1; k<=sne; k++) {
        getData(0, k, inval, inflg);
        if (fgout) putData(0, k, inval, inflg);
        if (fgincopie) 
          putData(1, k, inval, inflg);
        totnscount++;
      }

    if (c_meanspectra) { 
      vfftmean /= complex<r_8>((r_8)totnbblock, 0.);    
      pout.PutObject(vfftmean, "meanspectra");
      meanpowerspectra /= (r_8)totnbblock;
      pout.PutObject(vfftmean, "meanpowerspectra");
    }
    pout.PutObject(ffcoef, "fourierfilter");

    cout << " SimpleFourierFilter::run() - End of processing " 
         << " NbFFTBlocks= " << totnbblock << endl;
  }  // Bloc try 

  catch (PException & exc) {
    cerr << "SimpleFourierFilter: Catched Exception " << (string)typeid(exc).name()
         << "\n .... Msg= " << exc.Msg() << endl;
  }                                                                             
}

// ---------------------------------------------------------------
// -------------------- Classe SimpleFanOut -----------------------
// ---------------------------------------------------------------

SimpleFanOut::SimpleFanOut(int nbinput, int mfanout)
{
  if (nbinput < 1) 
    throw ParmError("SimpleFanOut::SimpleFanOut() NbInput < 1 !");
  if (mfanout < 1) 
    throw ParmError("SimpleFanOut::SimpleFanOut() M_FanOut < 1 !");

  nb_input = nbinput;
  m_fanout = mfanout;
  totnscount = 0;
}

SimpleFanOut::~SimpleFanOut()
{
}


void SimpleFanOut::PrintStatus(ostream & os)
{
  os << "\n ------------------------------------------------------ \n" 
     << " SimpleFanOut::PrintStatus() - NbInput=" << NbInput() 
     << " M_FanOut=" << MFanOut() << endl; 
  TOIProcessor::PrintStatus(os);
  os << endl;
  os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
  os << " ------------------------------------------------------ " << endl;
}

void SimpleFanOut::init() {
  cout << "SimpleFanOut::init NbInput=" << nb_input << endl;
  char buff[64];
  for(int k=0; k<nb_input; k++) {
    sprintf(buff,"in%d", k);
    declareInput(buff);
    for(int j=0; j<m_fanout; j++) {
      sprintf(buff,"out%d_%d", k, j);
      declareOutput(buff);
    }
  }

  name = "SimpleFanOut";
  //  upExtra = 1;
}

void SimpleFanOut::run() {
  //  TOIManager* mgr = TOIManager::getManager();
  int snb = getMinIn();
  int sne = getMaxIn();

  TVector<int_4> in_index(nb_input);
  TMatrix<int_4> out_index(nb_input, m_fanout);
  in_index = -1;
  out_index = -1;
  int nbconin = 0;
  char buff[64];
  for(int ki=0;ki<nb_input;ki++) {
    sprintf(buff,"in%d", ki);
    int idx = getInputTOIIndex(buff);
    if (!checkInputTOIIndex(idx)) continue;
    nbconin++;
    in_index(ki) = idx;
    bool fgout = false;
    for(int jo=0; jo<m_fanout; jo++) {
      sprintf(buff,"out%d_%d", ki, jo);
      int odx = getOutputTOIIndex(buff);
      if (checkOutputTOIIndex(odx)) {
        out_index(ki, jo) = odx;
        fgout = true;
      }
    }
    if (!fgout) {
      string msg_err = 
        "SimpleFanOut::run() - No connected Output for Input TOI (" 
        + getInName(ki) + ") !";
      cerr << msg_err << endl;
      throw ParmError(msg_err);
    }
  }
  if (nbconin == 0) {
    cerr << " SimpleFanOut::run() - No Input TOI connected! "
         << endl;
    throw ParmError("SimpleFanOut::run() No Inout TOI connected!");
  }

  /*
  for(int ki=0;ki<nb_input;ki++) {
    cout << " SimpleFanOut::run() In(" << ki << ") Index=" << in_index(ki) 
         << " Name=" << getInName(in_index(ki)) << endl;
    for(int jo=0; jo<m_fanout; jo++) 
      cout << " .... Out(" << ki << "," << jo << ") Index=" << out_index(ki, jo)
           << " Name=" << getOutName(out_index(ki, jo)) << endl;
  }
  */
  cout << " SimpleFanOut::run() SNRange=" << snb << " - " << sne << endl; 
  

  try {
    Timer tm("SimpleFanOut::run()");
    double valin = 0.;
    uint_8 fgin = 0;
    for(int k=snb;k<=sne;k++) {
      for(int i=0;i<nb_input;i++) {
        if (in_index(i) < 0) continue;
        valin = 0;
        fgin = 0;
        getData(in_index(i), k, valin, fgin);
        for(int j=0; j<m_fanout; j++) {
          if (out_index(i, j) < 0)  continue;
          putData(out_index(i, j), k, valin, fgin);
        }
      }
      totnscount++;
    }  // Boucle sur les num-sample
    cout << " SimpleFanOut::run() - End of processing " 
         << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
  }  // Bloc try 

  catch (PException & exc) {
    cerr << "SimpleFanOut: Catched Exception " << (string)typeid(exc).name()
         << "\n .... Msg= " << exc.Msg() << endl;
  }                                                                             
}