Index: /trunk/AddOn/TAcq/mfits2spec.cc
===================================================================
--- /trunk/AddOn/TAcq/mfits2spec.cc	(revision 3592)
+++ /trunk/AddOn/TAcq/mfits2spec.cc	(revision 3593)
@@ -5,5 +5,5 @@
 /* ------------------------------------------------------------------ 
    Programme de calcul de spectre moyenne a partir des fichiers fits  
-   d'acquisition de BAORadio - donnees brutes (non FFT)
+   d'acquisition de BAORadio - donnees brutes ou fft 
    R. Ansari, C. Magneville
    V1 : Juillet 2008  , V2 : Avril 2009
@@ -37,4 +37,8 @@
 #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
@@ -44,4 +48,6 @@
 // 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);
 
 //----------------------------------------------------
@@ -50,8 +56,15 @@
 {
   if (narg < 4) {
-    cout << " ---Calcul spectres moyennes a partir de fits BAORadio " << endl;
+    cout << " ---Mean spectrum computation from BAORadio FITS files" << 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;  
+    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;
   }
@@ -64,5 +77,21 @@
     string outppf = arg[2];
     vector<string> infiles;
-    for(int i=3; i<narg; i++) infiles.push_back(arg[i]);
+    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;
@@ -70,4 +99,5 @@
     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 ; 
@@ -93,9 +123,10 @@
 }
 
-inline r_4 Zmod2(complex<r_4> z) 
+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;
@@ -115,6 +146,6 @@
 
     Byte* data = new Byte[sx];
-    TVector<r_4> vx(sx);
-    TVector< complex<r_4> > cfour;
+    TVector<DTF> vx(sx);
+    TVector< complex<DTF> > cfour;
     FFTPackServer ffts;
       
@@ -122,5 +153,5 @@
       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]);
+      for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (DTF)(data[ix]);
       // On fait le FFT
       ffts.FFTForward(vx, cfour);
@@ -139,5 +170,5 @@
   }   
   cout << "---- Ecriture spectre moyenne ds " << outfile << endl;      
-  spectre /= (r_4)(nzm);
+  spectre /= (DTF)(nzm);
   POutPersist po(outfile);
   po << spectre;
@@ -147,12 +178,17 @@
 /*--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;
+  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];
@@ -181,11 +217,11 @@
     size_t sy = mff.NAxis2();
     BRPaquet paq(NULL, data, paqsz);
-    TVector<r_4> vx(paq.DataSize());
-    TVector< complex<r_4> > cfour;
+    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) = (r_4)(paq.Data1()[ix])-127.5;
+      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);
@@ -197,6 +233,11 @@
       nzm++;  freq0 += cfour(0).real();
       }
-    nfileok++;
-    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;      
+    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;
@@ -206,5 +247,5 @@
   } 
   cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;      
-  spectre /= (r_4)(nzm);
+  spectre /= (DTF)(nzm);
   spectre.Info().Comment() = " SpectreMoyenne (Moyenne module^2) ";
   spectre.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
@@ -217,4 +258,5 @@
 /*--Nouvelle-Fonction--*/
 int ana_data_2(vector<string>& infiles, string& outfile)
+// Calcul spectres moyens 2 voies, donnees brutes  
 {
   TVector<float> specV1, specV2;
@@ -223,8 +265,12 @@
   float freq0v2 = 0;  
   int paqsz = 0;
-  int nfileok;
+  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];
@@ -254,7 +300,7 @@
     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;
+    TVector<DTF> vx1(paq.DataSize()/2);
+    TVector<DTF> vx2(paq.DataSize()/2);
+    TVector< complex<DTF> > cfour1,cfour2;
       
     for(int j=0; j<sy; j++) {
@@ -262,6 +308,6 @@
      // 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;
+      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);
@@ -280,5 +326,10 @@
       }
     nfileok++;
-    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;      
+    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;
@@ -288,9 +339,9 @@
   } 
   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 /= (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;
@@ -302,2 +353,108 @@
   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;
+}
