Changeset 3905 in Sophya for trunk/AddOn/TAcq/brproc.cc

Timestamp:

Oct 15, 2010, 11:02:36 PM (15 years ago)

Author:

ansari

Message:

implementation de la fenetre en temps pour les spectres ds BRMeanSpecCalculator, Reza 15/10/2010

File:

• trunk/AddOn/TAcq/brproc.cc (modified) (11 diffs)

trunk/AddOn/TAcq/brproc.cc

@@ -32 +32 @@
 
 //---------------------------------------------------------------------
-// Classe de traitement simple - calcul de spectres moyennes / voie
+// Classe de traitement de spectres - 
+// Calcul de spectres moyennes,variance / voie + nettoyage
 //---------------------------------------------------------------------
 /* --Methode-- */
@@ -39 +40 @@
   : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean),
     fgdatafft_(fgdatafft), fgsinglechannel_(fgsinglechan), 
-    clnflg_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()),
+    nbpaq4mean_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()), 
     nbadpaq_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()) 
 {
@@ -59 +60 @@
 
   numfile_=0;
-  nbpaq4mean_=0;
   totnbpaq_=0;
 
+  size_t nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres());
+  
+  for(size_t i=0; i<nchan; i++) {
+    nbpaq4mean_[i]=nbadpaq_[i]=0;
+  }
+
+  // Definition des tailles de fenetres de spectres, etc ...
+  SetSpectraWindowSize();
+  SetMaxNbSepcWinFiles();
+  nbtot_specwin_=0;
   SetVarianceLimits();
 
@@ -69 +79 @@
   ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create);
   dtp_ = new SwFitsDataTable(*ofsdtp_,1024,true);
-  int nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres());
   char cnom[32];
   for(int i=0; i<nchan; i++) {
-    sprintf(cnom,"var%d",i);
+    sprintf(cnom,"variance%d",i);
     dtp_->AddFloatColumn(cnom);    
   }
+  /*
   for(int i=0; i<nchan; i++) {
-    sprintf(cnom,"varnorm%d",i);
+    sprintf(cnom,"sigma%d",i);
     dtp_->AddFloatColumn(cnom);    
   }
+  */
   xnt_=new double[nchan*2];
+
 }
 
@@ -85 +97 @@
 BRMeanSpecCalculator::~BRMeanSpecCalculator()
 {
-  if (nbpaq4mean_>1)  SaveSpectra();
+  uint_8 npqm=0;
+  for(size_t i=0; i<nbpaq4mean_.size(); i++)  npqm+=nbpaq4mean_[i];
+  if (npqm>nmean_*nbpaq4mean_.size()/10)  SaveMeanSpectra();
   cout << " ---------------- BRMeanSpecCalculator()_Finalizing -------------------- " << endl;
   for(size_t i=0; i<nbadpaq_.size(); i++) {
@@ -100 +114 @@
 
 /* --Methode-- */
-void BRMeanSpecCalculator::ReadGainFitsFile(string filename)
+void BRMeanSpecCalculator::SetSpectraWindowSize(uint_4 winsz, uint_4 wszext)
+{
+  if (winsz < 3) {
+    winsz=1;  wszext=0;
+  }  
+  if (wszext>=winsz/2) wszext=winsz/2;
+  sa_size_t sz[5]={0,0,0,0,0};
+  sz[0]=mspecmtx_.NCols();
+  sz[1]=mspecmtx_.NRows();
+  sz[2]=winsz+2*wszext;
+  spec_window_.SetSize(3,sz);
+  spwin_ext_sz_=wszext;
+  sz[0]=mspecmtx_.NRows();
+  sz[1]=winsz+2*wszext;  
+  clnflg_.SetSize(2,sz);
+  cout << "BRMeanSpecCalculator::SetSpectraWindowSize()/Info: SpectraWindowSize()=" << GetSpectraWindowSize()
+       << " ExtensionSize=" << GetSpecWinExtensionSize() << " Overlap=" << GetSpecWinOverlapSize() 
+       << " ArraySize=" << spec_window_.SizeZ() << endl;
+
+  paqnum_w_start=spwin_ext_sz_;  // premiere initialisation du numero de paquet 
+  return;
+}
+
+/* --Methode-- */
+void BRMeanSpecCalculator::ReadGainFitsFile(string filename, bool fgapp)
 {
   cout << " BRMeanSpecCalculator::ReadGainFitsFile() - reading file " << filename;
   FitsInOutFile fis(filename, FitsInOutFile::Fits_RO);
   fis >> sgain_;
-  cout << " MeanGain=" << sgain_.Sum()/sgain_.Size() << endl;
+  fg_apply_gains_=fgapp;
+  cout << " MeanGain=" << sgain_.Sum()/sgain_.Size() << " ApplyGains=" 
+       << ((fg_apply_gains_)?"true":"false") << endl;
 }
 
@@ -111 +151 @@
 { return (z.real()*z.real()+z.imag()*z.imag()); }
 
+
+
+
 /* --Methode-- */
 int BRMeanSpecCalculator::Process()
 {
-  if (nbpaq4mean_==nmean_)  SaveSpectra();
-  FlagBadPackets();
-    
+  // Cette methode remplit le tableau spec_window_ avec les spectres (renormalise avec 
+  // les gains si demande) et appelle la methode du traitement de la fenetre temporelle 
+  // des spectres le cas echeant ProcSpecWin()
+
+  int_8 nbpaqdec = (int_8)totnbpaq_-(int_8)GetSpecWinOverlapSize();
+  if ((nbpaqdec>0)&&(nbpaqdec%GetSpectraWindowSize()==0)) {
+    paqnum_w_end=totnbpaq_-GetSpecWinExtensionSize();
+    ProcSpecWin(paqnum_w_start, paqnum_w_end);
+    paqnum_w_start=totnbpaq_-GetSpecWinExtensionSize();
+  }
+
   if (fgdatafft_) {  // Donnees firmware FFT 
-    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) {
-      if (clnflg_[i])  { nbadpaq_[i]++;  continue; }   // si le paquet a ete flagge mauvais ( clnflg_[i] <> 0 ) 
+    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
       TwoByteComplex* zp=NULL;
       if (fgsinglechannel_) {
@@ -128 +178 @@
         if (i%2==1)  zp=vpaq_[i/2].Data2C();
       }
-      TVector< r_4 > spec = mspecmtx_.Row(i); 
-      TVector< r_4 > sspec = sigspecmtx_.Row(i); 
-      for(sa_size_t f=1; f<spec.Size(); f++) {
-        r_4 m2zf=zp[f].module2F();;
-        spec(f) += m2zf;
-        sspec(f) += m2zf*m2zf;
-      }
+      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
+      for(sa_size_t f=0; f<spec_window_.SizeX(); f++) 
+        spec_window_(f,i,kz) = zp[f].module2F();
     }
   }
   else {  // Donnees RAW qui ont du etre processe par BRFFTCalculator
-    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) {
-      if (clnflg_[i])  { nbadpaq_[i]++;  continue; }   // si le paquet a ete flagge mauvais ( clnflg_[i] <> 0 ) 
+    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) {
       complex<ODT>* zp=NULL;
       if (fgsinglechannel_) {
@@ -148 +193 @@
         if (i%2==1)  zp= reinterpret_cast< complex<ODT>* >(vprocpaq_[i/2]+memgr_.ProcPaqSize()/2) ;
       }
+      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
+      for(sa_size_t f=0; f<spec_window_.SizeX(); f++) 
+        spec_window_(f,i,kz) = Zmod2(zp[f]); 
+    }
+  }
+  if (fg_apply_gains_) {   // Application des gains, si demande 
+    sa_size_t kz=PaqNumToArrayIndex(totnbpaq_);
+    for(sa_size_t i=0; i<spec_window_.SizeY(); i++) 
+      (spec_window_(Range::all(), Range(i), Range(kz))).Div(sgain_.Row(i));
+  }
+
+  totnbpaq_++;
+  return 0;
+}
+
+
+/* --Methode-- */
+void BRMeanSpecCalculator::ProcSpecWin(uint_8 numpaqstart, uint_8 numpaqend)
+{
+  //DBG  cout << "BRMeanSpecCalculator::ProcSpecWin()/Debug: numpaqstart=" << numpaqstart 
+  //DBG     << " numpaqend=" << numpaqend << endl;
+
+  // On appelle la routine de nettoyage qui doit flagger les mauvais paquets 
+  FlagBadPackets(numpaqstart, numpaqend);
+
+  // Boucle sur les numeros de paquets de la fenetre en temps
+  for (uint_8 jp=numpaqstart; jp<numpaqend; jp++) {
+    // On sauvegarde les spectres moyennes si necessaire 
+    if ((nbpaq4mean_[0]>0)&&(nbpaq4mean_[0]%nmean_ == 0))  SaveMeanSpectra();
+    //  On peut aussi acceder aux spectres et flags pour (jpmin -),(jpmax+) GetSpecWinExtensionSize()
+    sa_size_t kz=PaqNumToArrayIndex(jp); 
+    // Boucle sur les numeros de voie (canaux) 
+    for(sa_size_t i=0; i<spec_window_.SizeY(); i++)  {
+      if ( clnflg_(i,kz) != 0)  continue;
       TVector< r_4 > spec = mspecmtx_.Row(i); 
       TVector< r_4 > sspec = sigspecmtx_.Row(i); 
-      for(sa_size_t f=1; f<spec.Size(); f++)  {
-        r_4 m2zf=Zmod2(zp[f]);
-        spec(f) += m2zf;
-        sspec(f) += (m2zf*m2zf);        
-      }     
+      // Calcul de spectres moyennes et variance 
+      for(sa_size_t f=1; f<spec.Size(); f++)  {   // boucle sur les frequences 
+        spec(f) += spec_window_(f,i,kz);
+        sspec(f) += spec_window_(f,i,kz)*spec_window_(f,i,kz);  
+      }
+      nbpaq4mean_[i]++;  // compteur de paquets OK pour la moyenne 
     }
   }
-  nbpaq4mean_++;   totnbpaq_++;
-  return 0;
-}
-
-/* --Methode-- */
-void BRMeanSpecCalculator::FlagBadPackets()
-{
-  if (fgdatafft_) {  // Donnees firmware FFT 
-    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) {
-      TwoByteComplex* zp=NULL;
-      if (fgsinglechannel_) {
-        zp=vpaq_[i].Data1C();
-      }
-      else {
-        zp=vpaq_[i/2].Data1C();
-        if (i%2==1)  zp=vpaq_[i/2].Data2C();
-      }
+  if (nbtot_specwin_<nmaxfiles_specw_)  SaveSpectraWindow();
+  nbtot_specwin_++;
+  return;
+}
+
+/* --Methode-- */
+void BRMeanSpecCalculator::FlagBadPackets(uint_8 numpaqstart, uint_8 numpaqend)
+{
+  // Boucle sur les numeros de paquets de la fenetre en temps
+  for (uint_8 jp=numpaqstart; jp<numpaqend; jp++) {
+    //  On peut aussi acceder aux spectres et flags pour (jpmin -),(jpmax+) GetSpecWinExtensionSize()
+    sa_size_t kz=PaqNumToArrayIndex(jp); 
+    // Boucle sur les numeros de voie (canaux) 
+    for(sa_size_t i=0; i<spec_window_.SizeY(); i++)  {
+      double mean, sigma;
+      sa_size_t kz=PaqNumToArrayIndex(totnbpaq_); 
       double variance=0.;
-      double varnorm=0.;
-      for(sa_size_t f=1; f<mspecmtx_.NCols(); f++) {
-        double modsq=(double)zp[f].module2F();
-        variance+=modsq;
-        varnorm+=modsq/sgain_(i,f);
-      }
+      variance=spec_window_(Range(1,Range::lastIndex()), Range(i), Range(kz)).Sum();
       xnt_[i]=variance;
-      xnt_[i+mspecmtx_.NRows()]=varnorm;
-      clnflg_[i]=0;
-      if (varnorm<varmin_) clnflg_[i]=1;
-      else if (varnorm>varmax_) clnflg_[i]=2;
+      clnflg_(i,kz)=0;
+      if (variance<varmin_) { clnflg_(i,kz)=1;  nbadpaq_[i]++; }
+      else if (variance>varmax_) { clnflg_(i,kz)=2;  nbadpaq_[i]++; }
     }
-  }
-  else {  // Donnees RAW qui ont du etre processe par BRFFTCalculator
-    for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) {
-      complex<ODT>* zp=NULL;
-      if (fgsinglechannel_) {
-        zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i]);
-      }
-      else {
-        zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i/2]);
-        if (i%2==1)  zp= reinterpret_cast< complex<ODT>* >(vprocpaq_[i/2]+memgr_.ProcPaqSize()/2) ;
-      }
-      double variance=0.;
-      double varnorm=0.;
-      for(sa_size_t f=1; f<mspecmtx_.NCols(); f++) {
-        double modsq=(double)Zmod2(zp[f]);
-        variance+=modsq;
-        varnorm+=modsq/sgain_(i,f);
-      }
-      xnt_[i]=variance;
-      xnt_[i+mspecmtx_.NRows()]=varnorm;
-      clnflg_[i]=0;
-      if (varnorm<varmin_) clnflg_[i]=1;
-      else if (varnorm>varmax_) clnflg_[i]=2;
+  dtp_->AddRow(xnt_);
+  }
+  return;
+}
+
+/* --Methode-- */
+void BRMeanSpecCalculator::SaveMeanSpectra()
+{
+  for(sa_size_t ir=0; ir<mspecmtx_.NRows(); ir++) {
+    char buff[32];
+    sprintf(buff,"NPAQSUM_%d",(int)ir);
+    mspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_[0];
+    mspecmtx_.Info()[buff] = nbpaq4mean_[ir];
+    sigspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_[0];
+    sigspecmtx_.Info()[buff] = nbpaq4mean_[ir];
+    if (nbpaq4mean_[ir] > 0) {
+      mspecmtx_.Row(ir)  /= (r_4)nbpaq4mean_[ir];
+      sigspecmtx_.Row(ir) /= (r_4)nbpaq4mean_[ir];
+      sigspecmtx_.Row(ir) -= (mspecmtx_.Row(ir) && mspecmtx_.Row(ir));  // Mean(X^2) - [ Mean(X) ]^2
     }
   }
-  dtp_->AddRow(xnt_);
-  return;
-}
-
-/* --Methode-- */
-void BRMeanSpecCalculator::SaveSpectra()
-{
-  mspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_;
-  mspecmtx_ /= (double)nbpaq4mean_;
-  sigspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_;
-  sigspecmtx_ /= (double)nbpaq4mean_;
-  sigspecmtx_ -= (mspecmtx_ && mspecmtx_);  // Mean(X^2) - [ Mean(X) ]^2
   char nfile[64];
   string flnm;
-  /*
-  {
-  sprintf(nfile,"mspecmtx%d.ppf",numfile_);
-  flnm=outpath_+nfile;
-  POutPersist po(flnm);
-  po << mspecmtx_;
-  }
-  {
-  sprintf(nfile,"sigspecmtx%d.ppf",numfile_);
-  flnm=outpath_+nfile;
-  POutPersist po(flnm);
-  po << sigspecmtx_;
-  }
-  */
   {
   sprintf(nfile,"mspecmtx%d.fits",numfile_);
@@ -253 +296 @@
   }
 
-  cout << numfile_ << "-BRMeanSpecCalculator::SaveSpectra() NPaqProc="
+  cout << numfile_ << "-BRMeanSpecCalculator::SaveMeanSpectra() NPaqProc="
        << totnbpaq_ << "  -> Mean/Sig spectra Matrix in " << flnm << " /sigspec...ppf" << endl;
-  numfile_++;  nbpaq4mean_=0;
+  numfile_++;  
+
+  for(size_t i=0; i<nbpaq4mean_.size(); i++) nbpaq4mean_[i]=0;
   mspecmtx_ = (r_4)(0.);
   sigspecmtx_ = (r_4)(0.);
   return;
+}
+
+/* --Methode-- */
+void BRMeanSpecCalculator::SaveSpectraWindow()
+{
+  char nfile[64];
+  string flnm;
+  sprintf(nfile,"specwin%d.fits",nbtot_specwin_);
+  flnm="!"+outpath_+nfile;
+  FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create);
+  fos << spec_window_;
+  cout << " SaveSpectraWindow() " << nbtot_specwin_ << "- file " << nfile << " created " << endl;
 }
 
@@ -332 +389 @@
   for(uint_4 k=0; k<inp.Size(); k++)   inp(k)=(ODT)indata[k];
   fftwf_execute(myplan_); 
-  for(uint_4 k=0; k<outfc.Size(); k++)   ofc[k]=outfc(k);
+  for(uint_4 k=0; k<outfc.Size(); k++)   ofc[k]=outfc(k)/(ODT)sz_;   // on renormalise les coeff FFT ( / sz )
   return;
 }