Changeset 3592 in Sophya for trunk/AddOn/TAcq/mfits2spec.cc

Timestamp:

Apr 6, 2009, 11:57:39 PM (16 years ago)

Author:

ansari

Message:

Amelioration programme mfits2pec.cc - Reza 06/04/2009

File:

• trunk/AddOn/TAcq/mfits2spec.cc (modified) (4 diffs)

trunk/AddOn/TAcq/mfits2spec.cc

@@ -2 +2 @@
 #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++
@@ -26 +33 @@
 
 
-// include mini-fits lib 
+// include mini-fits lib , et structure paquet BAORadio
 #include "minifits.h"
-
-
-inline r_4 Zmod2(complex<r_4> z) 
-{ return (z.real()*z.real()+z.imag()*z.imag()); }
-
+#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 < 2) {
-    cout << " Usage: mfits2spec file1.fits [file2.fits ...] " << endl;
+  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;
   }
-  cout << " ---------- mfits2spec.cc Start ------------- " << endl;
 
   TArrayInitiator  _inia;
@@ -45 +61 @@
   int rc = 0;
   try {
-        ResourceUsage resu;
-        TVector<float> spectre;
-        sa_size_t nzm = 0;   // Nb de spectres moyennes
-    for(int ifile=1; ifile<narg; ifile++) {
-      string ffname = arg[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;      
-    }   
-    cout << "---- Ecriture spectre moyenne ds mspectre.ppf " << endl;      
-    spectre /= (r_4)(nzm);
-    POutPersist po("mspectre.ppf");
-    po << spectre;
+    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 ; 
   }
@@ -109 +92 @@
 
 }
+
+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;
+}