source: Sophya/trunk/AddOn/TAcq/mfits2spec.cc@ 3631

// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ------------------------------------------------------------------ 
   Programme de calcul de spectre moyenne a partir des fichiers fits  
   d'acquisition de BAORadio - donnees brutes ou fft 
   R. Ansari, C. Magneville
   V1 : Juillet 2008  , V2 : Avril 2009
   ------------------------------------------------------------------ */

// include standard c/c++
#include <math.h>
#include <stdio.h>

#include <iostream>
#include <string>

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

#include "FFTW/fftw3.h"

// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"


// include mini-fits lib , et structure paquet BAORadio
#include "minifits.h"
#include "brpaqu.h"


// Definition de type pour les tableaux de float 
#define DTF  r_4 

//---- Declaration des fonctions de calcul ----
// Fonction d'analyse 1ere version, pas d'entete ds le fichier, 1 voie
int ana_data_0(vector<string>& infiles, string& outfile);
// Fonction d'analyse 2eme version , 1 voie / paquet
int ana_data_1(vector<string>& infiles, string& oufile);
// Fonction d'analyse 2eme version , 2 voies / paquet
int ana_data_2(vector<string>& infiles, string& oufile);
// Fonction d'analyse 2eme version , 2 voies / paquet
int ana_data_fft_2(vector<string>& infiles, string& oufile);

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if (narg < 4) {
    cout << " ---Mean spectrum computation from BAORadio FITS files" << endl;
    cout << " Usage: mfits2spec ACT OutPPF file1 [file2 file3 ...] " << endl;
    cout << "  Or :  mfits2spec ACT OutPPF -infits DirName Imin Imax " << endl;
    cout << " ACT=-0,-1,-2,-fft2 \n" 
         << "   -0: Nancay-July2008 \n" 
         << "   -1,-2 :  Raw data  1 ou 2 channels / packet(or frame) \n" 
         << "   -fft2 :  FFT data 2 channels / packet " << endl;
    cout << " OutPPF : Output PPF file name " << endl;  
    cout << " file1,file2 ... : Input FITS files " << endl;  
    cout << " -infiles DirName Imin Imax : Input fits directory and sequence numbers \n" 
         << "    FileNames=DirName/signalII.fits Imin<=II<=Imax " << endl;  
    return 1;
  }

  TArrayInitiator  _inia;

  int rc = 0;
  try {
    string act = arg[1];
    string outppf = arg[2];
    vector<string> infiles;
    if (strcmp(arg[3],"-infits")==0) {
      if (narg<7)  { 
        cout << " mfits2spec/Error arguments - type mfits2spec for help" << endl;
        return 2;
      }
      char nbuff[1024];
      char* dirname = arg[4]; 
      int imin = atoi(arg[5]);
      int imax = atoi(arg[6]);
      for(int ii=imin; ii<=imax; ii++) {
        sprintf(nbuff,"%s/signal%d.fits",dirname,ii);
        infiles.push_back(nbuff);
      }
    }
    else // Liste de noms de fichiers fournis directement 
      for(int i=3; i<narg; i++) infiles.push_back(arg[i]);

    cout << " ---------- mfits2spec.cc Start - ACT= " << act << " ------------- " << endl;
    ResourceUsage resu;
    if (act == "-0") ana_data_0(infiles, outppf);
    else if (act == "-1") ana_data_1(infiles, outppf);
    else if (act == "-2") ana_data_2(infiles, outppf);
    else if (act == "-fft2") ana_data_fft_2(infiles, outppf);
    else cout << " mfits2spec.cc / Bad argument ACT=" << act << " -> exit" << endl;
    cout << resu ; 
  }
  catch (MiniFITSException& exc) {
    cerr << " mfits2spec.cc catched MiniFITSException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n mfits2spec.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " mfits2spec.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

  cout << ">>>> mfits2spec.cc ------- FIN ----------- RC=" << rc << endl;
  return rc;

}

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

/*--Nouvelle-Fonction--*/
int ana_data_0(vector<string>& infiles, string& outfile)
// Calcul spectre moyen 1 voie, donnees brutes, Donnees Juillet 2008 
{
  TVector<float> spectre;
  sa_size_t nzm = 0;   // Nb de spectres moyennes
  for(int ifile=0; ifile<infiles.size(); ifile++) {
    string ffname = infiles[ifile];
// -------------- Lecture de bytes      
    cout << ifile <<"-Ouverture/lecture fichier " << ffname << endl;
    MiniFITSFile mff(ffname, MF_Read);
    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1() 
         << " NAxis2=" << mff.NAxis2() << endl;
    if (mff.DataType() != MF_Byte) {
      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
    }
    size_t sx = mff.NAxis1(); 
    size_t sy = mff.NAxis2();

    Byte* data = new Byte[sx];
    TVector<DTF> vx(sx);
    TVector< complex<DTF> > cfour;
    FFTPackServer ffts;
      
    for(int j=0; j<sy; j++) {
      mff.ReadB(data, sx, j*sx);
      // On convertit en float 
      for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (DTF)(data[ix]);
      // On fait le FFT
      ffts.FFTForward(vx, cfour);
      if (!spectre.IsAllocated())  spectre.SetSize(cfour.Size());
        
      // On cumule les modules carres
      for(sa_size_t jf=1; jf<spectre.Size(); jf++) 
        spectre(jf) += Zmod2(cfour(jf)); 
      nzm++;

//      cout << " j=" << j << " data[0...3]=" << (int)data[0] << " , " 
//           << (int)data[1] << " , " <<  (int)data[2] << endl;
      }
    cout << "---- FIN lecture " << ffname << endl;      
    delete[] data;
  }   
  cout << "---- Ecriture spectre moyenne ds " << outfile << endl;      
  spectre /= (DTF)(nzm);
  POutPersist po(outfile);
  po << spectre;
  return 0;
}

/*--Nouvelle-Fonction--*/
int ana_data_1(vector<string>& infiles, string& outfile)
// Calcul spectre moyen 1 voie, donnees brutes  
{
  TVector<float> spectre;
  float freq0 = 0;  
  int paqsz = 0;
  int nfileok=0;
  sa_size_t nzm = 0;   // Nb de spectres moyennes
  Byte* data = NULL;
  FFTPackServer ffts;

  Timer tm("ana_data_1");
  size_t totnbytesrd = 0;

  for(int ifile=0; ifile<infiles.size(); ifile++) {
    string ffname = infiles[ifile];
// -------------- Lecture de bytes      
    cout << "ana_data_1[" << ifile <<"]Ouverture/lecture fichier " << ffname << endl;
    MiniFITSFile mff(ffname, MF_Read);
    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1() 
         << " NAxis2=" << mff.NAxis2() << endl;
    if (mff.DataType() != MF_Byte) {
      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
      continue;
    }
// Les fichier FITS contiennent l'entet (24 bytes), mais pas le trailer (16 bytes) ...
    if (paqsz == 0)  {  // premier passage, on fixe la taille de paquet et on alloue le buffer
      paqsz = mff.NAxis1()+16;
      data = new Byte[paqsz];
      for(int ib=0; ib<paqsz; ib++) data[ib]=0;
    }
    else {
      if (paqsz != mff.NAxis1()+16) {
      cout << " PB : paqsz=" << paqsz << " != mff.NAxis1()+16 --> skipping " << endl;
      continue;
      }
    }
    size_t sx = mff.NAxis1(); 
    size_t sy = mff.NAxis2();
    BRPaquet paq(NULL, data, paqsz);
    TVector<DTF> vx(paq.DataSize());
    TVector< complex<DTF> > cfour;
      
    for(int j=0; j<sy; j++) {
      mff.ReadB(data, sx, j*sx);
      // On convertit en float 
      for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (DTF)(paq.Data1()[ix])-127.5;
      // On fait le FFT
      ffts.FFTForward(vx, cfour);
      if (!spectre.IsAllocated())  spectre.SetSize(cfour.Size());
        
      // On cumule les modules carres  - en sautant la frequence zero 
      for(sa_size_t jf=1; jf<spectre.Size(); jf++) 
        spectre(jf) += Zmod2(cfour(jf)); 
      nzm++;  freq0 += cfour(0).real();
      }
    nfileok++;   
    size_t nbytesrd = sx*sy;
    totnbytesrd += nbytesrd;
    tm.Split();
    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl; 
    cout << "  TotalDiskRead= " << totnbytesrd/(1024*1024) << " MBytes Disk-Read rate= " 
         << (double)(nbytesrd)/(1024.*1024.)/tm.PartialCPUTime() << " MB/s" << endl;    
  }   
  if (data) delete[] data;
  if (nzm <= 0) {
    cout << " ana_data_1/ERROR : nzm=" << nzm << " spectres moyennes !" << endl;
    return nzm;
  } 
  cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;      
  spectre /= (DTF)(nzm);
  spectre.Info().Comment() = " SpectreMoyenne (Moyenne module^2) ";
  spectre.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
  spectre.Info()["Freq0"] = freq0;
  POutPersist po(outfile);
  po << PPFNameTag("specV1") << spectre;
  return nfileok;
}

/*--Nouvelle-Fonction--*/
int ana_data_2(vector<string>& infiles, string& outfile)
// Calcul spectres moyens 2 voies, donnees brutes  
{
  TVector<float> specV1, specV2;
  TVector< complex<float> > cxspecV12;
  float freq0v1 = 0;  
  float freq0v2 = 0;  
  int paqsz = 0;
  int nfileok=0;
  sa_size_t nzm = 0;   // Nb de spectres moyennes
  Byte* data = NULL;
  FFTPackServer ffts;

  Timer tm("ana_data_2");
  size_t totnbytesrd = 0;

  for(int ifile=0; ifile<infiles.size(); ifile++) {
    string ffname = infiles[ifile];
// -------------- Lecture de bytes      
    cout << "ana_data_2[" << ifile <<"]Ouverture/lecture fichier " << ffname << endl;
    MiniFITSFile mff(ffname, MF_Read);
    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1() 
         << " NAxis2=" << mff.NAxis2() << endl;
    if (mff.DataType() != MF_Byte) {
      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
      continue;
    }
// Les fichier FITS contiennent l'entet (24 bytes), mais pas le trailer (16 bytes) ...
    if (paqsz == 0)  {  // premier passage, on fixe la taille de paquet et on alloue le buffer
      paqsz = mff.NAxis1()+16;
      cout << " ana_data_2/ Allocating data , PaqSz=" << paqsz << endl;
      data = new Byte[paqsz];
      for(int ib=0; ib<paqsz; ib++) data[ib]=0;
    }
    else {
      if (paqsz != mff.NAxis1()+16) {
      cout << " PB : paqsz=" << paqsz << " != mff.NAxis1()+16 --> skipping " << endl;
      continue;
      }
    }
    size_t sx = mff.NAxis1(); 
    size_t sy = mff.NAxis2();
    BRPaquet paq(NULL, data, paqsz);
    TVector<DTF> vx1(paq.DataSize()/2);
    TVector<DTF> vx2(paq.DataSize()/2);
    TVector< complex<DTF> > cfour1,cfour2;
      
    for(int j=0; j<sy; j++) {
      mff.ReadB(data, sx, j*sx);
     // if (j%50 == 0) cout << " *DBG* mff.ReadB() , j=" << j << endl;
     // On convertit en float 
      for(sa_size_t ix=0; ix<vx1.Size(); ix++)  vx1(ix) = (DTF)(paq.Data1()[ix])-127.5;
      for(sa_size_t ix=0; ix<vx2.Size(); ix++)  vx2(ix) = (DTF)(paq.Data2()[ix])-127.5;
      // On fait le FFT
      ffts.FFTForward(vx1, cfour1);
      ffts.FFTForward(vx2, cfour2);
      if (!specV1.IsAllocated())  specV1.SetSize(cfour1.Size());
      if (!specV2.IsAllocated())  specV2.SetSize(cfour2.Size());
      if (!cxspecV12.IsAllocated())  cxspecV12.SetSize(cfour1.Size());
        
      // On cumule les modules carres  - en sautant la frequence zero 
      for(sa_size_t jf=1; jf<specV1.Size(); jf++) {
        specV1(jf) += Zmod2(cfour1(jf)); 
        specV2(jf) += Zmod2(cfour2(jf)); 
        cxspecV12(jf) += cfour1(jf)*cfour2(jf);
      }
      nzm++;  freq0v1 += cfour1(0).real();  freq0v2 += cfour2(0).real(); 
      }
    nfileok++;
    size_t nbytesrd = sx*sy;
    totnbytesrd += nbytesrd;
    tm.Split();
    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl; 
    cout << "  TotalDiskRead= " << totnbytesrd/(1024*1024) << " MBytes Disk-Read rate= " 
         << (double)(nbytesrd)/(1024.*1024.)/tm.PartialCPUTime() << " MB/s" << endl;    
  }  
  if (data) delete[] data;
  if (nzm <= 0) {
    cout << " ana_data_2/ERROR : nzm=" << nzm << " spectres moyennes !" << endl;
    return nzm;
  } 
  cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;      
  specV1 /= (DTF)(nzm);
  specV2 /= (DTF)(nzm);
  cxspecV12 /= complex<DTF>((DTF)(nzm),0.);
  specV1.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - Raw-Voie 1 (/2)";
  specV2.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - Raw-Voie 2 (/2)";
  specV1.Info()["NMOY"] = specV2.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
  specV1.Info()["Freq0"] = freq0v1;
  specV2.Info()["Freq0"] = freq0v2;
  POutPersist po(outfile);
  po << PPFNameTag("specV1") << specV1;
  po << PPFNameTag("specV2") << specV2;
  po << PPFNameTag("cxspecV12") << cxspecV12;
  return 0;
}


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

/*--Nouvelle-Fonction--*/
int ana_data_fft_2(vector<string>& infiles, string& outfile)
// Calcul spectres moyens 2 voies, donnees brutes  
{
  TVector<float> specV1, specV2;
  TVector< complex<float> > cxspecV12;
  float freq0v1 = 0;  
  float freq0v2 = 0;  
  int paqsz = 0;
  int nfileok=0;
  sa_size_t nzm = 0;   // Nb de spectres moyennes
  Byte* data = NULL;

  Timer tm("ana_data_fft_2");
  size_t totnbytesrd = 0;

  for(int ifile=0; ifile<infiles.size(); ifile++) {
    string ffname = infiles[ifile];
// -------------- Lecture de bytes      
    cout << "ana_data_fft_2[" << ifile <<"]Ouverture/lecture fichier " << ffname << endl;
    MiniFITSFile mff(ffname, MF_Read);
    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1() 
         << " NAxis2=" << mff.NAxis2() << endl;
    if (mff.DataType() != MF_Byte) {
      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
      continue;
    }
// Les fichier FITS contiennent l'entet (24 bytes), mais pas le trailer (16 bytes) ...
    if (paqsz == 0)  {  // premier passage, on fixe la taille de paquet et on alloue le buffer
      paqsz = mff.NAxis1()+16;
      cout << " ana_data_fft_2/ Allocating data , PaqSz=" << paqsz << endl;
      data = new Byte[paqsz];
      for(int ib=0; ib<paqsz; ib++) data[ib]=0;
    }
    else {
      if (paqsz != mff.NAxis1()+16) {
      cout << " PB : paqsz=" << paqsz << " != mff.NAxis1()+16 --> skipping " << endl;
      continue;
      }
    }
    size_t sx = mff.NAxis1(); 
    size_t sy = mff.NAxis2();
    BRPaquet paq(NULL, data, paqsz);
    TVector<DTF> vx1(paq.DataSize()/2);
    TVector<DTF> vx2(paq.DataSize()/2);
    TVector< complex<DTF> > cfour1,cfour2;

    if (!specV1.IsAllocated())  specV1.SetSize(paq.DataSize()/4+1);
    if (!specV2.IsAllocated())  specV2.SetSize(paq.DataSize()/4+1);
    if (!cxspecV12.IsAllocated())  cxspecV12.SetSize(paq.DataSize()/4+1);
      
    for(int j=0; j<sy; j++) {
      mff.ReadB(data, sx, j*sx);

      TwoByteComplex* zz1 = (TwoByteComplex*)paq.Data1();
      TwoByteComplex* zz2 = (TwoByteComplex*)paq.Data2();
      // Composante continue, partie reelle uniquement
      specV1(0) += zz1[0].realD()*zz1[0].realD();  
      specV2(0) += zz2[0].realD()*zz2[0].realD();  
      cxspecV12(0) += zz1[0].realD()*zz2[0].realD();

      for(sa_size_t jf=1; jf<specV1.Size()-1; jf++)  { 
        specV1(jf) += Zmod2TwoByte(zz1[jf]);
        specV2(jf) += Zmod2TwoByte(zz2[jf]);
        cxspecV12(jf) += (DTF)(zz1[jf].realD()*zz2[jf].realD()+zz1[jf].imagD()*zz2[jf].imagD());
      }
      // Freq. Nyquist a N/2  
      specV1(specV1.Size()-1) += zz1[0].imagD()*zz1[0].imagD();  // Freq. Nyquist a N/2 
      specV2(specV2.Size()-1) += zz2[0].imagD()*zz2[0].imagD();  // Freq. Nyquist a N/2 
      cxspecV12(cxspecV12.Size()-1) += zz1[0].imagD()*zz2[0].imagD();

      nzm++;  freq0v1 += zz1[0].realD();  freq0v2 += zz2[0].realD();
      }
    nfileok++;
    size_t nbytesrd = sx*sy;
    totnbytesrd += nbytesrd;
    tm.Split();
    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl; 
    cout << "  TotalDiskRead= " << totnbytesrd/(1024*1024) << " MBytes Disk-Read rate= " 
         << (double)(nbytesrd)/(1024.*1024.)/tm.PartialCPUTime() << " MB/s" << endl;    
  }  
  if (data) delete[] data;
  if (nzm <= 0) {
    cout << " ana_data_fft_2/ERROR : nzm=" << nzm << " spectres moyennes !" << endl;
    return nzm;
  } 
  cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;      
  specV1 /= (DTF)(nzm);
  specV2 /= (DTF)(nzm);
  cxspecV12 /= complex<DTF>((DTF)(nzm),0.);
  specV1.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - FFT-Voie 1 (/2)";
  specV2.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - FFT-Voie 2 (/2)";
  specV1.Info()["NMOY"] = specV2.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
  specV1.Info()["Freq0"] = freq0v1;
  specV2.Info()["Freq0"] = freq0v2;
  POutPersist po(outfile);
  po << PPFNameTag("specV1") << specV1;
  po << PPFNameTag("specV2") << specV2;
  po << PPFNameTag("cxspecV12") << cxspecV12;
  return 0;
}