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

// 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 (non 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"

//---- 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);

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if (narg < 4) {
    cout << " ---Calcul spectres moyennes a partir de fits BAORadio " << endl;
    cout << " Usage: mfits2spec ACT OutPPF file1 [file2 file3 ...] " << endl;
    cout << " ACT=-0,-1,-2 ==> 0: Nancay-Juil2008, - 1,2 : 1/2 voies / paquet " << endl;
    cout << " OutPPF : Output PPF file name, file1,file2 ... Input FITS files " << endl;  
    return 1;
  }

  TArrayInitiator  _inia;

  int rc = 0;
  try {
    string act = arg[1];
    string outppf = arg[2];
    vector<string> infiles;
    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 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 r_4 Zmod2(complex<r_4> z) 
{ return (z.real()*z.real()+z.imag()*z.imag()); }

/*--Nouvelle-Fonction--*/
int ana_data_0(vector<string>& infiles, string& outfile)
{
  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<r_4> vx(sx);
    TVector< complex<r_4> > 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) = (r_4)(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 /= (r_4)(nzm);
  POutPersist po(outfile);
  po << spectre;
  return 0;
}

/*--Nouvelle-Fonction--*/
int ana_data_1(vector<string>& infiles, string& outfile)
{
  TVector<float> spectre;
  float freq0 = 0;  
  int paqsz = 0;
  int nfileok;
  sa_size_t nzm = 0;   // Nb de spectres moyennes
  Byte* data = NULL;
  FFTPackServer ffts;
  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<r_4> vx(paq.DataSize());
    TVector< complex<r_4> > 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) = (r_4)(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++;
    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << 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 /= (r_4)(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)
{
  TVector<float> specV1, specV2;
  TVector< complex<float> > cxspecV12;
  float freq0v1 = 0;  
  float freq0v2 = 0;  
  int paqsz = 0;
  int nfileok;
  sa_size_t nzm = 0;   // Nb de spectres moyennes
  Byte* data = NULL;
  FFTPackServer ffts;
  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<r_4> vx1(paq.DataSize()/2);
    TVector<r_4> vx2(paq.DataSize()/2);
    TVector< complex<r_4> > 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) = (r_4)(paq.Data1()[ix])-127.5;
      for(sa_size_t ix=0; ix<vx2.Size(); ix++)  vx2(ix) = (r_4)(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++;
    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << 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 /= (r_4)(nzm);
  specV2 /= (r_4)(nzm);
  cxspecV12 /= complex<r_4>((r_4)(nzm),0.);
  specV1.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - Voie 1 (/2)";
  specV2.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - 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;
}