source: Sophya/trunk/AddOn/TAcq/racquproc.cc@ 3629

#include "racquproc.h"

#include <stdlib.h>
#include <unistd.h>
#include <fstream>
#include <signal.h>

#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"
#include "timestamp.h"
#include "fftpserver.h"
#include "fftwserver.h"

#include "FFTW/fftw3.h"


#include "pciewrap.h"
#include "brpaqu.h"
#include "minifits.h"
struct sigaction act;

//-------------------------------------------------------
// Classe thread de traitement  avec 1 voie par frame
//-------------------------------------------------------

void DataProcSignal(int s)
{
  cout <<"............................................... receive CtrlC" << endl;

}

DataProc::DataProc(RAcqMemZoneMgr& mem, string& path, uint_4 nmean, uint_4 step, uint_4 nmax)
  :  memgr(mem) 
{
  nmax_ = nmax; 
  nmean_ = nmean;
  step_ = step;
  stop_ = false;        
  path_ = path; 
  act.sa_handler=DataProcSignal;
}

inline r_4 Zmod2(complex<r_4> z) 
{ return (z.real()*z.real()+z.imag()*z.imag()); }

void DataProc::Stop()
{
 stop_=true;
 //  cout <<" DataProc::Stop ... > STOP " << endl;
}

void DataProc::run()
{
  
  // sigaddset(&act.sa_mask,SIGINT);  // pour proteger le transfert DMA
  // sigaction(SIGINT,&act,NULL);       
  setRC(1);     
  try {
    TimeStamp ts; 
    cout << " DataProc::run() - Starting " << ts << " NMaxMemZones=" << nmax_ 
         << " NMean=" << nmean_ << " Step=" << step_ << endl;   
    char fname[512];
    sprintf(fname,"%s/proc.log",path_.c_str());
    ofstream filog(fname);
    filog << " DataProc::run() - starting log file " << ts << endl;                    
    filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl; 

// Initialisation pour clcul FFT 
    TVector< complex<r_4> > cfour;  // composant TF
    uint_4 paqsz = memgr.PaqSize();
    BRPaquet pq(NULL, NULL, paqsz); 
    TVector<r_4> vx(pq.DataSize());
    vx = (r_4)(0.);
    FFTPackServer ffts;
    ffts.FFTForward(vx, cfour);
    TVector<r_4> spectre;
    spectre.ReSize(cfour.Size());
    
    fftwf_plan plan = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(), 
                                       (fftwf_complex *)cfour.Data(), FFTW_ESTIMATE); 

    uint_4 ifile = 0;                           
    uint_4 nzm = 0;
    for (uint_4 kmz=0; kmz<nmax_; kmz++) {
      if (stop_) break;
      int mid = memgr.FindMemZoneId(MemZA_Proc);
      Byte* buff = memgr.GetMemZone(mid);
      if (buff == NULL) {
        cout << " DataProc::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
        setRC(2);       
        return; 
      }
      BRPaquet paq0(NULL, buff, paqsz); 
      for(uint_4 i=0; i<memgr.NbPaquets(); i+=step_) {
        BRPaquet paq(NULL, buff+i*paqsz, paqsz); 
        Byte min = 255;
        Byte max = 0;
            
        for(sa_size_t j=0; j<vx.Size(); j++) 
          vx(j) = (r_4)(*(paq.Data1()+j))-127.5;
        fftwf_execute(plan); 
        //        ffts_.FFTForward(vx, cfour_);
        for(sa_size_t j=0; j<spectre.Size(); j++) 
          spectre(j) += Zmod2(cfour(j)); 
        nzm++;
      }
      if ((nzm >= nmean_) ||(kmz==(nmax_-1))) {
        spectre /= (r_4)(nzm);
        spectre.Info()["NPaqMoy"] = nzm;
        sprintf(fname,"%s/spectre%d.ppf",path_.c_str(),(int)ifile);
        POutPersist po(fname);
        po << spectre; 
        spectre = (r_4)(0.);
        nzm = 0;  ifile++; 
        ts.SetNow();
        filog << ts << " :  proc file  " << fname << endl;                     
        cout << " DataProc::run() " << ts << " : created file  " << fname << endl;
      }
      
      memgr.FreeMemZone(mid, MemZS_Proc);
    }
  }
  catch (PException& exc) {
    cout << " DataProc::run()/catched PException " << exc.Msg() << endl;
    setRC(3);   
    return; 
  }
  catch(...) {
    cout << " DataProc::run()/catched unknown ... exception " << endl;
    setRC(4);   
    return; 
  }
  setRC(0);
  return;
}   

//---------------------------------------------------------------
// Classe thread de traitement donnees ADC avec 2 voies par frame
//---------------------------------------------------------------

DataProc2C::DataProc2C(RAcqMemZoneMgr& mem, string& path, uint_4 nmean, uint_4 step, uint_4 nmax)
  :  memgr(mem) 
{
  nmax_ = nmax; 
  nmean_ = nmean;
  step_ = step;
  stop_ = false;        
  path_ = path; 
  act.sa_handler=DataProcSignal;
}
void DataProc2C::Stop()
{
 stop_=true;
 // cout <<" DataProc2C::Stop ... > STOP " << endl;
}

void DataProc2C::run()
{
  //   sigaction(SIGINT,&act,NULL);     
  setRC(1);     
  try {
    TimeStamp ts; 
    cout << " DataProc2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_ 
         << " NMean=" << nmean_ << " Step=" << step_ << endl;   
    char fname[512];
    sprintf(fname,"%s/proc.log",path_.c_str());
    ofstream filog(fname);
    filog << " DataProc2C::run() - starting log file " << ts << endl;                  
    filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl; 

// Initialisation pour clcul FFT 
    TVector< complex<r_4> > cfour;  // composant TF
    uint_4 paqsz = memgr.PaqSize();
    BRPaquet pq(NULL, NULL, paqsz); 
    TVector<r_4> vx(pq.DataSize()/2);
    vx = (r_4)(0.);
    FFTPackServer ffts;
    ffts.FFTForward(vx, cfour);
    TVector<r_4> spectreV1, spectreV2;
    spectreV1.ReSize(cfour.Size());
    spectreV2.ReSize(cfour.Size());
    
    fftwf_plan plan = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(), 
                                       (fftwf_complex *)cfour.Data(), FFTW_ESTIMATE); 

    uint_4 ifile = 0;                           
    uint_4 nzm = 0;
    for (uint_4 kmz=0; kmz<nmax_; kmz++) {
      if (stop_) break;
      int mid = memgr.FindMemZoneId(MemZA_Proc);
      Byte* buff = memgr.GetMemZone(mid);
      if (buff == NULL) {
        cout << " DataProc2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
        setRC(2);       
        return; 
      }
      BRPaquet paq0(NULL, buff, paqsz); 
      for(uint_4 i=0; i<memgr.NbPaquets(); i+=step_) {
        BRPaquet paq(NULL, buff+i*paqsz, paqsz); 
        Byte min = 255;
        Byte max = 0;
            
        for(sa_size_t j=0; j<vx.Size(); j++) 
          vx(j) = (r_4)(*(paq.Data1()+j))-127.5;
        fftwf_execute(plan); 
        //   ffts_.FFTForward(vx, cfour_);
        for(sa_size_t j=0; j<spectreV1.Size(); j++) 
          spectreV1(j) += Zmod2(cfour(j));
 
        for(sa_size_t j=0; j<vx.Size(); j++) 
          vx(j) = (r_4)(*(paq.Data2()+j))-127.5;
        fftwf_execute(plan); 
        //    ffts_.FFTForward(vx, cfour_);
        for(sa_size_t j=0; j<spectreV2.Size(); j++) 
          spectreV2(j) += Zmod2(cfour(j)); 

        nzm++;
      }
      if ((nzm >= nmean_) ||(kmz==(nmax_-1))) {
        spectreV1 /= (r_4)(nzm);
        spectreV2 /= (r_4)(nzm);
        spectreV1.Info()["NPaqMoy"] = nzm;
        spectreV2.Info()["NPaqMoy"] = nzm;
        {
        sprintf(fname,"%s/spectre2C_%d.ppf",path_.c_str(),(int)ifile);
        POutPersist po(fname);
        po << PPFNameTag("specV1") << spectreV1; 
        po << PPFNameTag("specV2") << spectreV2; 
        }
        spectreV1 = (r_4)(0.);
        spectreV2 = (r_4)(0.);
        nzm = 0;  ifile++; 
        ts.SetNow();
        filog << ts << " :  proc file  " << fname << endl;                     
        cout << " DataProc2C::run() " << ts << " : created file  " << fname << endl;
      }
      
      memgr.FreeMemZone(mid, MemZS_Proc);
    }
  }
  catch (PException& exc) {
    cout << " DataProc2C::run()/catched PException " << exc.Msg() << endl;
    setRC(3);   
    return; 
  }
  catch(...) {
    cout << " DataProc2C::run()/catched unknown ... exception " << endl;
    setRC(4);   
    return; 
  }
  setRC(0);
  return;
}   




//---------------------------------------------------------------
// Classe thread de traitement donnees FFT avec 2 voies par frame
//---------------------------------------------------------------

inline r_4 Zmod2TwoByte(TwoByteComplex z) 
{ return (z.realD()*z.realD()+z.imagD()*z.imagD()); }

DataProcFFT2C::DataProcFFT2C(RAcqMemZoneMgr& mem, string& path, uint_4 nmean, uint_4 step, uint_4 nmax)
  :  memgr(mem) 
{
  nmax_ = nmax; 
  nmean_ = nmean;
  step_ = step;
  stop_ = false;        
  path_ = path; 
 act.sa_handler=DataProcSignal;
}
void DataProcFFT2C::Stop()
{
 stop_=true;
 // cout <<" DataProcFFT2C::Stop ... > STOP " << endl;
}

void DataProcFFT2C::run()
{
  // sigaction(SIGINT,&act,NULL);       
  setRC(1);     
  try {
    TimeStamp ts; 
    cout << " DataProcFFT2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_ 
         << " NMean=" << nmean_ << " Step=" << step_ << endl;   
    char fname[512];
    sprintf(fname,"%s/proc.log",path_.c_str());
    ofstream filog(fname);
    filog << " DataProcFFT2C::run() - starting log file " << ts << endl;                       
    filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl; 

// Initialisation pour clcul FFT 
    TVector< complex<r_4> > cfour;  // composant TF
    uint_4 paqsz = memgr.PaqSize();
    BRPaquet pq(NULL, NULL, paqsz); 
 
    TVector<r_4> spectreV1(pq.DataSize()/4+1), spectreV2(pq.DataSize()/4+1);    

    uint_4 ifile = 0;                           
    for (uint_4 kmz=0; kmz<nmax_; kmz++) {
      if (stop_ ) break;
      int mid = memgr.FindMemZoneId(MemZA_Proc);
      Byte* buff = memgr.GetMemZone(mid);
      if (buff == NULL) {
        cout << " DataProcFFT2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
        setRC(2);       
        return; 
      }
      BRPaquet paq0(NULL, buff, paqsz); 
      uint_4 nzm = 0;
      for(uint_4 i=0; i<memgr.NbPaquets(); i+=step_) {
        BRPaquet paq(NULL, buff+i*paqsz, paqsz); 
        Byte min = 255;
        Byte max = 0;
                
        TwoByteComplex* zz;
        // Traitement Voie 1
        zz = (TwoByteComplex*)paq.Data1();
        spectreV1(0) += zz[0].realD()*zz[0].realD();  // Composante continue, partie reelle uniquement
        for(sa_size_t j=1; j<spectreV1.Size()-1; j++)  spectreV1(j) += Zmod2TwoByte(zz[j]);

        spectreV1(spectreV1.Size()-1) += zz[0].imagD()*zz[0].imagD();  // Freq. Nyquist a N/2 

        // Traitement Voie 2
        zz = (TwoByteComplex*)paq.Data2();
        spectreV2(0) += zz[0].realD()*zz[0].realD();  // Composante continue, partie reelle uniquement
        for(sa_size_t j=1; j<spectreV2.Size()-1; j++)  spectreV2(j) += Zmod2TwoByte(zz[j]);

        spectreV2(spectreV2.Size()-1) += zz[0].imagD()*zz[0].imagD();  // Freq. Nyquist a N/2 

        nzm++;
      }
      if ((nzm >= nmean_) ||(kmz==(nmax_-1))) {
        spectreV1 /= (r_4)(nzm);
        spectreV2 /= (r_4)(nzm);
        spectreV1.Info()["NPaqMoy"] = nzm;
        spectreV2.Info()["NPaqMoy"] = nzm;
        {
        sprintf(fname,"%s/spectreFFT2C_%d.ppf",path_.c_str(),(int)ifile);
        POutPersist po(fname);
        po << PPFNameTag("specV1") << spectreV1; 
        po << PPFNameTag("specV2") << spectreV2; 
        }
        spectreV1 = (r_4)(0.);
        spectreV2 = (r_4)(0.);
        nzm = 0;  ifile++; 
        ts.SetNow();
        filog << ts << " :  proc file  " << fname << endl;                     
        cout << " DataProcFFT2C::run() " << ts << " : created file  " << fname << endl;
      }
      
      memgr.FreeMemZone(mid, MemZS_Proc);
    }
  }
  catch (PException& exc) {
    cout << " DataProcFFT2C::run()/catched PException " << exc.Msg() << endl;
    setRC(3);   
    return; 
  }
  catch(...) {
    cout << " DataProcFFT2C::run()/catched unknown ... exception " << endl;
    setRC(4);   
    return; 
  }
  setRC(0);
  return;
}