Changeset 1479 in Sophya

Timestamp:

Apr 27, 2001, 2:52:15 AM (24 years ago)

Author:

ansari

Message:

debug deglitcher - Ajout Filtre Fourier (SimpleFourierFilter) - Reza 27/4/2001

Location:

trunk/ArchTOIPipe

Files:

• ProcWSophya/simtoipr.cc (modified) (9 diffs)
• ProcWSophya/simtoipr.h (modified) (6 diffs)
• TestPipes/simtst.cc (modified) (15 diffs)

trunk/ArchTOIPipe/ProcWSophya/simtoipr.cc

@@ -5 +5 @@
 #include "pexceptions.h"
 #include "ctimer.h"
-
-SimpleDeglitcher::SimpleDeglitcher(int wsz, double ns, int maxnpt)
+#include "fftpserver.h"
+
+SimpleDeglitcher::SimpleDeglitcher(int wsz, double ns, int maxnpt, int minnpt)
 {
   SetWSize(wsz);
-  SetDetectionParam(ns, ns/2., maxnpt);
+  SetDetectionParam(ns, ns/2., maxnpt, minnpt);
   SetRange(-9.e39, 9.e39);
   RepBadSamples(true, true);
@@ -24 +25 @@
 }
 
-void SimpleDeglitcher::SetDetectionParam(double ns, double ns2, int maxnpt, int wszrec)
+void SimpleDeglitcher::SetDetectionParam(double ns, double ns2, int maxnpt, 
+                                         int minnpt, int wszrec)
 {
   nsig = (ns > 0.01) ? ns : 1.;
   nsig2 = (ns2 > 0.01) ? ns2 : nsig/2.;
   maxpoints = ((maxnpt > 0) && (maxnpt <= wsize)) ? maxnpt : 5;
+  minpoints = ((minnpt > 0) && (minnpt <= maxpoints)) ? minnpt : maxpoints/2;
   wrecsize = ((wszrec > 0) && (wszrec <= wsize)) ? wszrec : 2*maxpoints;
 }
@@ -42 +45 @@
   os << "\n ------------------------------------------------------ \n" 
      << " SimpleDeglitcher::PrintStatus() - WindowSize=" << WSize() 
-     << "  NbSigmas=" << NbSigmas() << " MaxPoints=" << MaxPoints() << endl;
-  os << "  NbSigmas2=" << NbSigmas2() << "  WRecSize= " << WRecSize() 
+     << " MaxPoints=" << MaxPoints() << " MinPoints=" << MinPoints() << endl;
+  os << " NbSigmas=" << NbSigmas() 
+     << "  NbSigmas2=" << NbSigmas2() << "  WRecSize= " << WRecSize() 
      << "  Range_Min= " << range_min << " Range_Max= " << range_max << endl;
   os << "  RepOutOfRangeSamples: " << _Bool2YesNo(rec_out_range_samples) 
@@ -117 +121 @@
       getData(0, k+snb, vin(k), vfg(k));
 
-    int nokdebut = 0.;
+    int nokdebut = 0;
     double s = 0.;    
     double mean = 0.;
@@ -236 +240 @@
       if (valcur < mean+curnsig*sigma) {  // inferieur au seuil
         if (fgglitch) {  
-          if (k-kgl < maxpoints) {   // On vient de detecter un glitch
+          if ( (k-kgl <= maxpoints) && (k-kgl >= minpoints) ) {   
+            // On vient de detecter un glitch
             glcount++;
             if (rec_gl_samples) { // On change la valeur des samples
@@ -253 +258 @@
             lastput = snb+k-1;
           }  // - Fin de detection de glitch
-          else {  // Trop long - ce n'est pas un glitch ...
+          else {  // Trop long ou trop court - ce n'est pas un glitch ...
             for(ii=kgl; ii<k; ii++) {
               putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
@@ -267 +272 @@
       else {  // Superieur au seuil
         if (fgglitch) {
-          if (k-kgl+1 >= maxpoints) {  // serie de points > seuil
+          if (k-kgl+1 > maxpoints) {  // serie de points > seuil
             for(ii=kgl; ii<=k; ii++)   // -> Donc pas glitch
               putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
@@ -600 +605 @@
 }
 
+// ----------------------------------------------------------------------
+//   Classe SimpleFourierFilter : Filtre simple ds le domaine de Fourier
+// ----------------------------------------------------------------------
+
+SimpleFourierFilter::SimpleFourierFilter(Vector const & vc)
+{
+  ffcoef = vc;
+  wsize = (ffcoef.Size()-1)*2;
+  if (wsize < 16) 
+    throw ParmError("SimpleFourierFilter::SimpleFourierFilter() WSize<16 !");
+  KeepSpectra(0);
+}
+
+SimpleFourierFilter::~SimpleFourierFilter()
+{
+}
+
+
+void SimpleFourierFilter::PrintStatus(ostream & os)
+{
+  os << "\n ------------------------------------------------------ \n" 
+     << " SimpleFourierFilter::PrintStatus() - WindowSize=" 
+     << WSize() << endl; 
+  TOIProcessor::PrintStatus(os);
+  os << " Coeff= " << ffcoef << endl;
+  os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
+  os << " ------------------------------------------------------ " << endl;
+}
+
+void SimpleFourierFilter::init() {
+  cout << "SimpleFourierFFilter::init" << endl;
+  declareInput("in");
+  declareOutput("out");
+  declareOutput("incopie");
+  name = "SimpleFourierFilter";
+  //  upExtra = 1;
+}
+
+
+void SimpleFourierFilter::run() {
+  //  TOIManager* mgr = TOIManager::getManager();
+  int snb = getMinIn();
+  int sne = getMaxIn();
+
+  bool fgout = checkOutputTOIIndex(0);
+  bool fgincopie = checkOutputTOIIndex(1);
+
+  if (!checkInputTOIIndex(0)) {
+    cerr << " SimpleFourierFilter::run() - Input TOI (in) not connected! "
+         << endl;
+    throw ParmError("SimpleFourierFilter::run() Input TOI (in) not connected!");
+  }
+  if (!fgout) {
+    cerr << " SimpleFourierFilter::run() - No Output TOI connected! "
+         << endl;
+    throw ParmError("SimpleFourierFilter::run() No output TOI connected!");
+  }
+
+  cout << " SimpleFourierFilter::run() SNRange=" << snb << " - " << sne << endl; 
+  
+
+  try {
+    Timer tm("SimpleFourierFilter::run()");
+    // Le debut
+    Vector vin(wsize);
+    Vector vout(wsize);
+    TVector<int_8> vfg(wsize);
+    TVector< complex<r_8> > vfft;
+    TVector< complex<r_8> > zcoef(ffcoef.Size());
+    zcoef = ffcoef;
+
+
+    FFTPackServer ffts;
+    ffts.setNormalize(true);
+    // Boucle sur les sampleNum
+
+    int k,i,klast;
+    int nks = 0;
+
+    klast = snb-1;
+    // 1er partie, on traite par paquets de wsize
+    for(k=snb;k<=sne-wsize+1;k+=wsize) {
+      for(i=0; i<wsize; i++) 
+        getData(0, k+i, vin(i), vfg(i));
+      ffts.FFTForward(vin, vfft);
+      if (nks < nb_keep) {
+        string dbgfile = "ffiltdbg" + (string)MuTyV(nks) + ".ppf";
+        POutPersist po(dbgfile);
+        po << vfft; 
+        nks++;
+      }
+      vfft.MulElt(zcoef);
+      ffts.FFTBackward(vfft, vout);
+      for(i=0; i<wsize; i++) 
+        putData(0,k+i,vout(i),vfg(i));
+      if (fgincopie) 
+        for(i=0; i<wsize; i++) 
+          putData(1, k+i, vin(i), vfg(i));
+      klast+=wsize;
+      totnscount+=wsize;
+    }
+    
+
+    // 2eme partie, on traite la fin du bloc d'echantillons si necessaire
+    double inval;
+    int_8 inflg;
+    if (klast < sne) 
+      for(k=klast+1; k<=sne; k++) {
+        getData(0, k, inval, inflg);
+        putData(0, k, inval, inflg);
+        if (fgincopie) 
+          putData(1, k, inval, inflg);
+        totnscount++;
+      }
+    cout << " SimpleFourierFilter::run() - End of processing " << endl;
+  }  // Bloc try 
+
+  catch (PException & exc) {
+    cerr << "SimpleFourierFilter: Catched Exception " << (string)typeid(exc).name()
+         << "\n .... Msg= " << exc.Msg() << endl;
+  }                                                                             
+}
 
 // ---------------------------------------------------------------
@@ -659 +786 @@
   out_index = -1;
   int nbconin = 0;
-  bool fggin = false;
   char buff[64];
   for(int ki=0;ki<nb_input;ki++) {

trunk/ArchTOIPipe/ProcWSophya/simtoipr.h

@@ -14 +14 @@
 class SimpleDeglitcher : public TOIProcessor {
 public:
-                SimpleDeglitcher(int wsz=64, double ns=3, int maxnpt=5);
+                SimpleDeglitcher(int wsz=64, double ns=3, 
+                                 int maxnpt=5, int minnpt=2);
   virtual       ~SimpleDeglitcher();
 
@@ -26 +27 @@
 
 
-  void          SetDetectionParam(double ns, double ns2, int maxnpt, int wszrec=0);
+  void          SetDetectionParam(double ns, double ns2, int maxnpt, 
+                                  int minnpt, int wszrec=0);
     
   inline void   RepBadSamples(bool gl_samples, bool out_range_samples, bool use_wrec=true)
@@ -40 +42 @@
   inline double NbSigmas2() const { return nsig2; }
   inline int    MaxPoints() const { return maxpoints; }
+  inline int    MinPoints() const { return minpoints; }
   
   inline int_8  ProcessedSampleCount() const { return totnscount; }
@@ -61 +64 @@
   double nsig2;     // Seuil en nb de sigmas, pour les points suivants le 1er
   int maxpoints;    // Nb maxi de points > ns sigmas
+  int minpoints;    // Nb mini de points > ns sigmas pour avoir un glitch
   double range_min, range_max;  // Range acceptable pour in
 
@@ -136 +140 @@
 
 
+//  Un filtre simple, dans le domaine de Fourier
+//  InverseFFT ( FFT(Vecteur(in)) * FilterCoefficient )
+
+class SimpleFourierFilter : public TOIProcessor {
+public:
+                SimpleFourierFilter(Vector const & vc);  
+                ~SimpleFourierFilter();
+
+  inline int    WSize() const { return wsize; }
+  inline int_8  ProcessedSampleCount() const { return totnscount; }
+  inline Vector FilterCoefficients() const
+                { Vector rcv; rcv = ffcoef; return(rcv); } 
+
+  virtual void  PrintStatus(ostream & os) ; // const plus tard
+
+  virtual void  init();  
+  virtual void  run();
+
+  inline  void  KeepSpectra(int nb=0) 
+                { nb_keep = nb; }
+protected:
+  int_8 totnscount;   // Nombre total d'echantillon processe
+  int wsize;         // Taille de fenetre de travail
+  Vector ffcoef;     // Coefficients du filtre
+  int nb_keep;
+};
+
+
 //  Classe SimpleFanOut
-//  Transforme recopie chaque entree sur M lignes de sortie
+//  Recopie chaque entree sur M lignes de sortie
 
 class SimpleFanOut : public TOIProcessor {
@@ -159 +191 @@
 };
 
+
 #endif

trunk/ArchTOIPipe/TestPipes/simtst.cc

@@ -33 +33 @@
   else {
     cout << "\n Usage : simtst [-dbg] [-start snb] [-end sne] \n" 
-         << "         [-wtoi sz] [-w2 sz] inFitsName outFitsName ppfFileName \n" 
+         << "         [-wtoi sz] [-wdegli sz] [-wfft sz] [-keepfft n] \n"
+         << "         inFitsName outFitsName  \n" 
          << "   -dbg : sets TOISeqBuffered debug level to 1 \n"
          << "   -start snb : sets the start sample num \n"
@@ -42 +43 @@
          << "   -wtoi sz : sets TOISeqBuff buffer size (def= 8192)\n"
          << "   -wdegli sz : sets deglitcher window size (def= 512) \n"
+         << "   -wfft sz : Activate Fourier filter and sets its width \n"
+         << "   -keepfft n : Keeps n spectra (if Fourier filter activated) \n" 
          << endl;
     exit(0);
@@ -54 +57 @@
 
   bool fgdbg = false;
-  int w1 = 8192;
-  int w2 = 512;
+  int wtoi = 8192;
+  int wdegli = 512;
+  int wfft = 4096;
+  int keepfft = 0;
   int nmax = 10;
-  int istart = 104121000+w1*5;
+  int istart = 104121000+wtoi*5;
   int iend = 0;
   double range_min = -16000;
@@ -63 +68 @@
   string infile;
   string outfile;
-  string ppffile;
   string intoi = "boloMuV_27";
+  bool fg_f_filt = false;
 
   if (narg < 4) Usage(true);
@@ -80 +85 @@
     else if (strcmp(arg[ia],"-wtoi") == 0) {
       if (ia == narg-1) Usage(true);  
-      w1 = atoi(arg[ia+1]); ia++;
+      wtoi = atoi(arg[ia+1]); ia++;
     }    
     else if (strcmp(arg[ia],"-wdegli") == 0) {
       if (ia == narg-1) Usage(true);  
-      w2 = atoi(arg[ia+1]); ia++;
+      wdegli = atoi(arg[ia+1]); ia++;
+    }    
+    else if (strcmp(arg[ia],"-wfft") == 0) {
+      if (ia == narg-1) Usage(true);  
+      wfft = atoi(arg[ia+1]); ia++;
+      fg_f_filt = true;
+    }    
+    else if (strcmp(arg[ia],"-keepfft") == 0) {
+      if (ia == narg-1) Usage(true);  
+      keepfft = atoi(arg[ia+1]); ia++;
     }    
     else if (strcmp(arg[ia],"-range") == 0) {
@@ -100 +114 @@
   }
 
-  if (iend < istart) iend = istart+w1*(nmax+5);
-  if ((narg-ko) < 3)  Usage(true); 
+  if (iend < istart) iend = istart+wtoi*(nmax+5);
+  if ((narg-ko) < 2)  Usage(true); 
   infile = arg[ko];
   outfile = arg[ko+1];
-  ppffile = arg[ko+2];
 
   cout << " Initializing SOPHYA ... " << endl;
@@ -110 +123 @@
   InitTim();
 
-  cout << ">>>> tstrztoi: Infile= " << infile << " outFile=" << outfile 
-       << " ppfFile= " << ppffile << endl;
-  cout << ">>>> Window Size W1= " << w1 << " W2= " << w2
+  cout << ">>>> tstrztoi: Infile= " << infile << " outFile=" 
+       << outfile << endl; 
+  cout << ">>>> Window Size WTOI= " << wtoi << " WDEGLI= " << wdegli
        << "  iStart= " << istart << " iEnd= " << iend << endl;
   cout << ">>>> InTOIName= " << intoi << endl;
+
   try {
     TOIManager* mgr = TOIManager::getManager();
@@ -131 +145 @@
 
 
+    int w1 = wtoi;
+    int w2 = (fg_f_filt && (wfft > w1)) ? wfft : w1; 
+
     //    char * colname[5] = {"MJD", "UTC","boloMuV_11","magnFlux","pivot"};
     //    char * toiname[5] = {"MJD", "UTC","bolo11","magneto","pivot"};
@@ -139 +156 @@
     TOISeqBuffered * toimean = new TOISeqBuffered("mean", w1);
     if (fgdbg) toimean->setDebugLevel(1);
-    TOISeqBuffered * toisig = new TOISeqBuffered("sigma", w1);
+    TOISeqBuffered * toisig = new TOISeqBuffered("sigma", w2);
     if (fgdbg) toisig->setDebugLevel(1);
     TOISeqBuffered * toiincopie = new TOISeqBuffered("incopie", w1);
@@ -160 +177 @@
 
     cout << " Creating SimpleDeglitcher() " << endl;
-    SimpleDeglitcher degl(w2);
+    SimpleDeglitcher degl(wdegli);
+    degl.SetDetectionParam(3.,1.5, 4, 2, 10);
     cout << " Setting Range for deglitcher: " << range_min 
          << " - " << range_max << endl;
+
     degl.SetRange(range_min, range_max);
     degl.addInput("in", toiin);
@@ -174 +193 @@
     TOISeqBuffered * toidegli0 = new TOISeqBuffered("degli0", w1);
     if (fgdbg) toidegli0->setDebugLevel(1);
-    TOISeqBuffered * toidegli1 = new TOISeqBuffered("degli1", w1);
+    TOISeqBuffered * toidegli1 = new TOISeqBuffered("degli1", w2);
 
     TOISeqBuffered * toimean0 = new TOISeqBuffered("mean0", w1);
     if (fgdbg) toimean0->setDebugLevel(1);
-    TOISeqBuffered * toimean1 = new TOISeqBuffered("mean1", w1);
+    TOISeqBuffered * toimean1 = new TOISeqBuffered("mean1", w2);
     if (fgdbg) toimean1->setDebugLevel(1);
     
@@ -206 +225 @@
     TOISeqBuffered * toiout = new TOISeqBuffered("out", w1);
     if (fgdbg) toiout->setDebugLevel(1);
-    TOISeqBuffered * toideglioffcopie = new TOISeqBuffered("deglioffcopie", w1);
+    TOISeqBuffered * toideglioffcopie = new TOISeqBuffered("deglioffcopie", w2);
     if (fgdbg) toideglioffcopie->setDebugLevel(1);
     filt.addOutput("out", toiout);
     filt.addOutput("incopie", toideglioffcopie);
-    
+
+    Vector fcoef(wfft/2);
+    fcoef = RegularSequence(1., -2./wfft);
+    cout << " Creating Fourier Filter ... " << endl;
+    SimpleFourierFilter sfft(fcoef);
+    sfft.KeepSpectra(keepfft);
+    TOISeqBuffered * toifft = NULL;
+    TOISeqBuffered * toifiltcopie = NULL;
+
+    if (fg_f_filt) {
+      cout << " Connecting Fourier Filter ... " << endl;
+      toifft = new TOISeqBuffered("fftout", w1);
+      toifiltcopie = new TOISeqBuffered("filtcopie", w1);
+      sfft.addInput("in",toiout);
+      sfft.addOutput("out", toifft);
+      sfft.addOutput("incopie", toifiltcopie);
+    }
 
     cout << " Connecting to output (FitsWriter) " << endl;
@@ -216 +251 @@
     w.setOutFlags(true);
     w.addInput("in", toiincopie);
-    w.addInput("out", toiout);
+    if (fg_f_filt) {
+      w.addInput("filtout", toifiltcopie);
+      w.addInput("fftout", toifft);
+    } 
+    else   w.addInput("out", toiout);
     w.addInput("degli", toidegli1);
     w.addInput("deglioff", toideglioffcopie);
@@ -239 +278 @@
     adder.start();
     filt.start();
+    if (fg_f_filt) sfft.start();
     w.start();
 
@@ -268 +308 @@
     cout << filt;
     cout << adder;
+    if (fg_f_filt) cout << sfft ; 
 
   }