Ajout programme visi2tmfreq.cc , calcul de cartes temps(Y)-frequence(X) pour...

Ajout programme visi2tmfreq.cc , calcul de cartes temps(Y)-frequence(X) pour les 8 auto-correlations, protection pour g(nu)<GAINTHR=0.01 ds visi2ntac.cc visi2dtacx.cc, Reza 19/06/2015

README

@@ -16,3 +16,10 @@ List of files:
   with total power as a function of time for the 8 PAON4 auto-correlations 
     Updated April 2015: reads the gain g(nu) file produced by rdvisip4.cc , corrects the 
     auto-correlation signals before integrating in the specified bands 
+- visi2dtacx.cc : a main program to read visibilitiy files and produce a datatable 
+  with auto-correlation and 6 cross correlations (1H-2H 1H-3H 1H-4H 2H-3H 2H-4H 3H-4H) 
+  as a function of time for three frequency bins - Can also make a time-frequency map of 
+  the 6 cross-correlations. The visibilities are corrected for the electronic gain, using the 
+  gain curve g(nu) computed by rdvisip4.cc 
+- visi2tmfreq.cc : a main program to read visibilitiy files and produce time-frequency maps 
+  for the 8 auto-correlations, corrected for the gain, using the gain curve g(nu) computed by rdvisip4.cc 

makefile

@@ -3,12 +3,13 @@
 include $(SOPHYABASE)/include/sophyamake.inc
 
 #  Define our target list 
-all : Objs/rdvisip4  Objs/visi2ntac  Objs/visi2dtacx
+all : Objs/rdvisip4  Objs/visi2ntac  Objs/visi2dtacx  Objs/visi2tmfreq 
 
 clean :
 	rm -f Objs/* 
 
 ######
+## programme de remplissage de DataTable avec cross-correlation fct du temps 
 Objs/visi2dtacx : Objs/visi2dtacx.o Objs/visip4reader.o 
 	$(CXXLINK) -o Objs/visi2dtacx Objs/visi2dtacx.o Objs/visip4reader.o $(SOPHYAEXTSLBLIST) 
 
@@ -16,6 +17,7 @@ Objs/visi2dtacx.o : visi2dtacx.cc visip4reader.h
 	$(CXXCOMPILE)  -o Objs/visi2dtacx.o visi2dtacx.cc
 
 ######
+## programme de remplissage de DataTable (NTuple) avec auto-correlation fct du temps 
 Objs/visi2ntac : Objs/visi2ntac.o Objs/visip4reader.o 
 	$(CXXLINK) -o Objs/visi2ntac Objs/visi2ntac.o Objs/visip4reader.o $(SOPHYAEXTSLBLIST) 
 
@@ -23,6 +25,15 @@ Objs/visi2ntac.o : visi2ntac.cc visip4reader.h
 	$(CXXCOMPILE)  -o Objs/visi2ntac.o visi2ntac.cc
 
 ######
+## programme de calcul de matrice Visibilites V_ij(nu) moyennee et calcul de gain g(nu)  
+Objs/visi2tmfreq : Objs/visi2tmfreq.o Objs/visip4reader.o 
+	$(CXXLINK) -o Objs/visi2tmfreq Objs/visi2tmfreq.o Objs/visip4reader.o $(SOPHYAEXTSLBLIST) 
+
+Objs/visi2tmfreq.o : visi2tmfreq.cc visip4reader.h
+	$(CXXCOMPILE)  -o Objs/visi2tmfreq.o visi2tmfreq.cc
+
+######
+## programme de calcul de matrice Visibilites V_ij(nu) moyennee et calcul de gain g(nu)  
 Objs/rdvisip4 : Objs/rdvisip4.o Objs/visip4reader.o 
 	$(CXXLINK) -o Objs/rdvisip4 Objs/rdvisip4.o Objs/visip4reader.o $(SOPHYAEXTSLBLIST) 
 

visi2dtacx.cc

@@ -154,6 +154,7 @@ int main(int narg, char* arg[])
     cout << " visi2dtacx/Info: reading input gain file"<<gainfile<<" ..."; 
     Matrix gains;
     Vector gn;
+    double GAINTHR = 0.005;  // limite inferieur du gain avant faire 1/G
     {
       PInPersist pin(gainfile);
       pin>>PPFNameTag("gains")>>gains;
@@ -161,7 +162,7 @@ int main(int narg, char* arg[])
       for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
       // on transforme gains en facteur multiplicatif
       for(sa_size_t r=0;r<gains.NRows();r++) {
-	for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=1./gains(r,c);
+	for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=((gains(r,c)>GAINTHR)?1./gains(r,c):0.);
       }
       cout<<endl;
     }

visi2ntac.cc

@@ -100,6 +100,7 @@ int main(int narg, char* arg[])
   cout << " visi2ntac/Info: reading input gain file"<<gainfile<<endl;
   Matrix gains;
   Vector gn;
+  double GAINTHR = 0.005;  // limite inferieur du gain avant faire 1/G
   {
     PInPersist pin(gainfile);
     pin>>PPFNameTag("gains")>>gains;
@@ -107,7 +108,7 @@ int main(int narg, char* arg[])
     for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
     // on transforme gains en facteur multiplicatif
     for(sa_size_t r=0;r<gains.NRows();r++) {
-      for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=1./gains(r,c);
+      for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=((gains(r,c)>GAINTHR)?1./gains(r,c):0.);
     }
     cout<<endl;
   }

visi2tmfreq.cc

+// Utilisation de SOPHYA pour faciliter les tests ...
+#include "sopnamsp.h"
+#include "machdefs.h"
+
+/* ---------------------------------------------------------- 
+   Projet BAORadio/PAON4 - (C) LAL/IRFU  2008-2015 
+
+   visi2tmfreq.cc: programme de lecture des fichiers matrices de 
+   visibilites de PAON4 produits mfacq, creation d'une tableau 
+   temps-frequence pour chaque auto-correlation 
+   R. Ansari, T. Etourneau   -  LAL
+   V1 : Juin 2015 
+   ---------------------------------------------------------- */
+
+// include standard c/c++
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <iostream>
+#include <string>
+
+#include "pexceptions.h"
+#include "tvector.h"
+#include "fioarr.h"
+// #include "tarrinit.h"
+#include "ntuple.h" 
+#include "datatable.h" 
+#include "histinit.h" 
+#include "matharr.h" 
+#include "timestamp.h"
+#include <utilarr.h>
+
+// include sophya mesure ressource CPU/memoire ...
+#include "resusage.h"
+#include "ctimer.h"
+#include "timing.h"
+
+// include lecteur de fichiers visibilites 
+#include "visip4reader.h"
+
+ 
+//--------------------------- Fonctions de ce fichier   ------------------- 
+int Usage(bool fgshort=true);
+// int DecodeArgs(int narg, char* arg[]);
+
+/* --Fonction-- */
+int Usage(bool fgea)
+{
+  cout << " --- visi2tmfreq.cc : Read PPF files produced by mfacq to create time-frequency maps\n" << endl;
+  if (fgea) cout << " rdvisip4: Missing or wrong argument ! " << endl;
+  cout << " Usage: rdvisip4 InPathBAO5 InPathBAO6 Imin,Imax  GainPPFFile OutPPFFile [DeltaIAvg=60] [NBinFreq=4] [PrintLev=0] \n"
+       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
+       << "  Imin,Imax : read files for Imin<=iii<=Imax  \n"
+       << "  GainPPFFile: Input gains PPF file name \n"
+       << "  OutPPFFile: Output PPF file name which will contain the 8 auto-correlation time-frequency arrays \n"
+       << "  DeltaIAvg: compute average power every DeltaI input vismtx files, def=60 \n"
+       << "  NBinFreq: The binning in frequency in the computed time-frequency maps (def=4 -> 250 kHz) \n"
+       << endl;
+  /*
+  if (fgshort) {
+    cout << " rdvisip4 -h for detailed instructions" << endl;
+    return 1;
+  }
+  */
+  return 1;
+}
+
+//----------------------------------------------------
+//----------------------------------------------------
+int main(int narg, char* arg[])
+{
+  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
+  if (narg<6) return Usage(true);
+  string path5 = arg[1];
+  string path6 = arg[2];
+  int Imin,Imax,Istep=1; 
+  sscanf(arg[3],"%d,%d",&Imin,&Imax);
+  string gainfile=arg[4];
+  string outfile=arg[5];
+  Istep=1;
+
+  // taille de fenetre de moyenne en temps
+  int deltaIavg=60;
+  if ((narg>6)&&(strcmp(arg[6],"-")!=0))  deltaIavg=atoi(arg[6]);
+  // taille de fenetre de moyenne en frequence (en bin de frequence, cad 61 kHz) 
+  sa_size_t TFMfbin=4;
+  if ((narg>7)&&(strcmp(arg[7],"-")!=0))  TFMfbin=atoi(arg[7]);
+  int prtlev=0;
+  if (narg>8) prtlev=atoi(arg[8]);
+
+  cout << " visi2tmfreq/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
+       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<'\n'
+       << " GainFile: "<<gainfile<<" OutFile: "<<outfile
+       << " DeltaIAvg="<<deltaIavg<<" NBinFreq="<<TFMfbin<<" PrtLev="<<prtlev<<endl;
+  HiStatsInitiator _inia;
+  //   TArrayInitiator  _inia;
+
+  int rc = 0;
+  try {
+    ResourceUsage resu;
+    //  Numero des lignes des auto-correlations dans la matrice des visibilites 
+    sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
+ 
+    cout << " visi2tmfreq/Info: reading input gain file"<<gainfile<<" ..."; 
+    Matrix gains;
+    Vector gn;
+    double GAINTHR = 0.005;  // limite inferieur du gain avant faire 1/G
+    {
+      PInPersist pin(gainfile);
+      pin>>PPFNameTag("gains")>>gains;
+      pin>>PPFNameTag("gn")>>gn;
+      for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
+      // on transforme gains en facteur multiplicatif
+      for(sa_size_t r=0;r<gains.NRows();r++) {
+	for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=((gains(r,c)>GAINTHR)?1./gains(r,c):0.);
+      }
+      cout<<endl;
+    }
+    cout << " DONE"<<endl;
+
+    // Creation de l'instance de la classe lecteur des fichiers visimtxXXX.ppf 
+    vector<string> paths; 
+    paths.push_back(path5); 
+    paths.push_back(path6);
+    VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
+    vreader.setPrintLevel(prtlev);
+
+    bool fgok=true;
+    TMatrix< complex<r_4> > vismtx;
+    TMatrix< complex<r_4> > acsum;
+
+    double acdt[32];   // les 4*8=32 valeurs d'autocorrelation pour remplissage dans la table 
+    complex<double> cxdt[18];   // les 3*6=18 valeurs de cross-correlations pour remplissage dans la table 
+
+    TimeStamp dateobs, cfdate;
+    TimeStamp dateorg(2015,1,1,12,0,0.);  // Date origine 1 jan 2015
+    double mttag;
+    int cnt=0, cntnt=0, pcntnt=0;
+    int I=0; 
+    //----- 8 TimeFrequency maps of auto-correlation 
+    vector< TArray< r_4 > > vtfm;
+    
+    //---- for the time-freqency map filling    
+    sa_size_t TFMtmidx=0;
+    sa_size_t tfmSX, tfmSY;
+
+    while (fgok) {  // boucle de lecture des fichiers d'entree 
+      fgok=vreader.ReadNext(vismtx, cfdate, mttag);
+      if (fgok)  {
+	if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
+	  acsum.SetSize(8, vismtx.NCols());
+	  // Allocating time-frequency maps 
+	  tfmSX=vismtx.NCols()/TFMfbin;   // frequence selon l'axe X 
+	  tfmSY=(Imax-Imin)/Istep/deltaIavg;  // temps selon l'axe Y 
+	    cout << " visi2tmfreq/Info: allocating Time-Freqncy maps : Freq->NX= "<<tfmSX<<" x Time->NY="<<tfmSY<<endl;
+	    for(int k=0; k<6; k++) vtfm.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 
+	  }
+	}
+	if (I==0) {   // start filling a new DataTable row 
+	  dateobs=cfdate;
+	  if (prtlev>0) 
+	    cout << "visi2tmfreq/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
+	  acsum = complex<r_4>(0.,0.);
+	}
+	//   sum (integration) along the time axis 
+	for(int k=0; k<8; k++)    acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
+
+      //   sum (integration) along the time axis 
+      for(int k=0; k<8; k++)    acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
+      I++;    // incrementing DeltaTime counter 
+      if (I==deltaIavg) {
+	Vector vac(acsum.NCols());
+	// facteur de normalisation pour que notre carte temps-frequence corresponde a la moyenne et pas la somme 
+	double facnorm=1./(double)(TFMfbin*deltaIavg);  
+	for(int k=0; k<8; k++)  {  // loop over 8 auto-corr
+	  // correct for gain(nu)
+	  for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=acsum(k,c).real()*gains(k,c);
+	  //filling time-frequency maps  
+	  TArray< r_4 > & tfmap = vtfm[k];
+	  for(sa_size_t ix=0; ix<tfmSX; ix++) 
+	    tfmap(ix, TFMtmidx) = (vac( Range(ix*TFMfbin, (ix+1)*TFMfbin-1) ).Sum())*facnorm;
+	}  //---- end of loop over 8 auto-corr
+	TFMtmidx++;
+      }
+    }    //--- fin de boucle de lecture des fichiers d'entree 
+
+    cout<<" visi2tmfreq/Info: count="<<cnt<<" visimtx read "<<endl;
+    cout<<" visi2tmfreq/Info: Saving time-frequency maps to PPF file "<<outfile<<endl;
+    POutPersist potfm(outfile);
+    const char* tfm_names[8]={"TFM_AC0", "TFM_AC1", "TFM_AC2", "TFM_AC3", "TFM_AC4", "TFM_AC5", "TFM_AC6", "TFM_AC7"};
+    for(int k=0; k<8; k++)  {   // loop over the 8 auto-cor 
+      TArray< r_4 > & tfmap = vtfm[k];
+      potfm << PPFNameTag(tfm_names[k]) << tfmap;
+    }
+    //    resu.Update();
+    cout << resu;   // Update est fait lors du print
+  }
+  catch (PException& exc) {
+    cerr << " visi2tmfreq.cc catched PException " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n visi2tmfreq.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+    rc = 78;
+  }
+  catch (...) {
+    cerr << " visi2tmfreq.cc catched unknown (...) exception  " << endl; 
+    rc = 79; 
+  } 
+
+  cout << ">>>> visi2tmfreq.cc ------- END ----------- RC=" << rc << endl;
+  return rc;
+
+}
+