Changeset 1478 in Sophya for trunk/ArchTOIPipe


Ignore:
Timestamp:
Apr 26, 2001, 7:22:02 PM (24 years ago)
Author:
ansari
Message:

amelioration deglitcher - Reza 26/4/2001

Location:
trunk/ArchTOIPipe
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/ArchTOIPipe/ProcWSophya/simtoipr.cc

    r1476 r1478  
    88SimpleDeglitcher::SimpleDeglitcher(int wsz, double ns, int maxnpt)
    99{
    10   if (wsz < 5) wsz = 5;
    11   if (wsz%2 == 0) wsz++;
    12   if (maxnpt > wsz) maxnpt = wsz;
     10  SetWSize(wsz);
     11  SetDetectionParam(ns, ns/2., maxnpt);
     12  SetRange(-9.e39, 9.e39);
     13  RepBadSamples(true, true);
     14
     15  totnscount = glnscount = glcount = out_range_nscount = 0;
     16  deglitchdone = false;
    1317
    1418  cout << "SimpleDeglitcher::SimpleDeglitcher() WSize= " << wsz
    1519       << " nSig=" << ns << " maxNPt=" << maxnpt << endl;
    16  
    17   wsize = wsz;
    18   nsig = ns;
    19   maxpoints = maxnpt;
    20   totnscount = glnscount = glcount = out_range_nscount = 0;
    21   deglitchdone = false;
    22   SetRange(-9.e39, 9.e39);
    2320}
    2421
    … …  
    2724}
    2825
     26void SimpleDeglitcher::SetDetectionParam(double ns, double ns2, int maxnpt, int wszrec)
     27{
     28  nsig = (ns > 0.01) ? ns : 1.;
     29  nsig2 = (ns2 > 0.01) ? ns2 : nsig/2.;
     30  maxpoints = ((maxnpt > 0) && (maxnpt <= wsize)) ? maxnpt : 5;
     31  wrecsize = ((wszrec > 0) && (wszrec <= wsize)) ? wszrec : 2*maxpoints;
     32}
     33
     34inline char * _Bool2YesNo(bool fg)
     35{
     36  if (fg) return "YES" ;
     37  else return "NO" ;
     38}
    2939
    3040void SimpleDeglitcher::PrintStatus(ostream & os)
    … …  
    3343     << " SimpleDeglitcher::PrintStatus() - WindowSize=" << WSize()
    3444     << "  NbSigmas=" << NbSigmas() << " MaxPoints=" << MaxPoints() << endl;
    35   os << "  Range_Min= " << range_min << " Range_Max= " << range_max << endl;
     45  os << "  NbSigmas2=" << NbSigmas2() << "  WRecSize= " << WRecSize()
     46     << "  Range_Min= " << range_min << " Range_Max= " << range_max << endl;
     47  os << "  RepOutOfRangeSamples: " << _Bool2YesNo(rec_out_range_samples)
     48     << "  RepGlitchSamples: " << _Bool2YesNo(rec_gl_samples)
     49     << "  UseWRec: " << _Bool2YesNo(rec_use_wrec) << endl;
    3650  TOIProcessor::PrintStatus(os);
    3751  if (deglitchdone) os << " Deglitching performed " << endl;
    … …  
    5771}
    5872
     73
     74#define FG_OUTOFRANGE 1
     75#define FG_GLITCH 2
    5976void SimpleDeglitcher::run() {
    6077
    … …  
    88105    Timer tm("SimpleDeglitcher::run()");
    89106    Vector vin(wsize);
    90 
    91     int wrec = maxpoints*2;
    92     Vector vrec(wrec);
     107    Vector vas(wsize);
     108
     109    int wrecsize = maxpoints*2;
     110    Vector vrec(wrecsize);
    93111
    94112    TVector<int_8> vfg(wsize);
    … …  
    119137      vfg(k) = 0;
    120138    }
    121 
    122     for(k=0; k<wrec; k++) {
     139    for(k=0; k<wsize; k++) {
     140      vas(k) = mean;
     141      if ( vfg(k) != 0) continue;
     142      if ( (vin(k) < range_min) || (vin(k) > range_max) ) continue;
     143      vas(k) = vin(k);
     144    }
     145
     146    for(k=0; k<wrecsize; k++) {
    123147      if ( (vin(k) < range_min) || (vin(k) > range_max) ) vrec(k)=mean;
    124148      else vrec(k)=vin(k);   
    … …  
    144168      // Calcul mean-sigma
    145169      if (knext<=kfin)  {
    146         valsub = vin(knext%wsize);
     170        valsub = vas(knext%wsize);
    147171        getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));
    148172        valadd = vin(knext%wsize);
    149         if ( (valadd < range_min) || (valadd > range_max) ) {
    150           valadd = mean;
     173        if ( vfg(knext%wsize) ||
     174             (valadd < range_min) || (valadd > range_max) ) {
     175          vas(knext%wsize) = valadd = mean;
    151176          fgokcur = false;
    152177        }
    153         else fgokcur = true;
    154         if ( (valsub < range_min) || (valsub > range_max) )
    155           valsub = mean;
     178        else {
     179          vas(knext%wsize) = valadd = vin(knext%wsize);
     180          fgokcur = true;
     181        }
    156182        if (!fgokdebut && fgokcur) {
    157183          s += valadd;
    … …  
    185211        putData(3, k+snb, vin(k%wsize), vfg(k%wsize));
    186212
     213      valcur = vin(k%wsize);
    187214      if ( (valcur < range_min) || (valcur > range_max) ) {
    188         vin(k%wsize) = valcur = mean;
    189         vfg(k%wsize) |= 2;
     215        valcur = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
     216        if (rec_out_range_samples) vin(k%wsize) = valcur;
     217        vfg(k%wsize) |= FG_OUTOFRANGE;
    190218        out_range_nscount++;
    191219      }
    … …  
    204232//            }
    205233
    206       double curnsig = nsig;
    207       if (fgglitch)  curnsig /= 2.;
     234      double curnsig = (fgglitch) ? nsig2 : nsig;
    208235       
    209236      if (valcur < mean+curnsig*sigma) {  // inferieur au seuil
    … …  
    211238          if (k-kgl < maxpoints) {   // On vient de detecter un glitch
    212239            glcount++;
    213             double recval = vrec.Sum()/wrec;
    214             for(ii=kgl; ii<k; ii++) {
    215               putData(0, ii+snb, recval, vfg(ii%wsize)|5);
    216               glnscount++;
     240            if (rec_gl_samples) { // On change la valeur des samples
     241              double recval = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
     242              for(ii=kgl; ii<k; ii++) {
     243                putData(0, ii+snb, recval, vfg(ii%wsize)|FG_GLITCH);
     244                glnscount++;
     245              }
     246            }
     247            else {   // On ne fait que flagger les echantillons
     248              for(ii=kgl; ii<k; ii++) {
     249                putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize)|FG_GLITCH);
     250                glnscount++;
     251              }
    217252            }
    218253            lastput = snb+k-1;
    219           }
    220           else {
     254          }  // - Fin de detection de glitch
     255          else {  // Trop long - ce n'est pas un glitch ...
    221256            for(ii=kgl; ii<k; ii++) {
    222257              putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
    … …  
    225260          }
    226261        }
    227         putData(0, k+snb, valcur, 0);
     262        putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
    228263        lastput = snb+k;
    229264        kgl = -1;  fgglitch = false;
    230         vrec(k%wrec) = lastvalok = valcur;
     265        vrec(k%wrecsize) = lastvalok = valcur;
    231266      }
    232267      else {  // Superieur au seuil
    … …  
    237272            lastput = snb+k;
    238273            fgglitch = false; 
    239             vrec(k%wrec) = lastvalok = valcur;     
     274            vrec(k%wrecsize) = lastvalok = valcur;         
    240275          }
    241276        }
    … …  
    245280          }
    246281          else { // On est toujours dans une serie > seuil
    247             putData(0, k+snb, valcur, fgcur);
    248             lastput = snb+k;
     282            putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
     283            lastput = snb+k;   lastvalok = valcur;
    249284          }
    250           vrec(k%wrec) = lastvalok;         
     285          vrec(k%wrecsize) = lastvalok;     
    251286        }
    252287      }

  • trunk/ArchTOIPipe/ProcWSophya/simtoipr.h

    r1467 r1478  
    2222    { min = range_min; max = range_max; }
    2323
     24  inline void   SetWSize(int wsz)
     25    { wsize = (wsz < 5) ? 5 : wsz; }
     26
     27
     28  void          SetDetectionParam(double ns, double ns2, int maxnpt, int wszrec=0);
     29   
     30  inline void   RepBadSamples(bool gl_samples, bool out_range_samples, bool use_wrec=true)
     31    {  rec_gl_samples = gl_samples;  rec_out_range_samples = out_range_samples;
     32    rec_use_wrec =  use_wrec; }
     33
    2434  virtual void  init(); 
    2535  virtual void  run();
    2636
    2737  inline int    WSize() const { return wsize; }
     38  inline int    WRecSize() const { return wrecsize; }
    2839  inline double NbSigmas() const { return nsig; }
     40  inline double NbSigmas2() const { return nsig2; }
    2941  inline int    MaxPoints() const { return maxpoints; }
    3042 
    … …  
    4456
    4557  int wsize;        // Taille de fenetre de travail
     58  int wrecsize;     // Taille de fenetre de calcul pour reconstruite les mauvaises valeur
     59                    // pour valeur de glitch
    4660  double nsig;      // Seuil en nb de sigmas
     61  double nsig2;     // Seuil en nb de sigmas, pour les points suivants le 1er
    4762  int maxpoints;    // Nb maxi de points > ns sigmas
    4863  double range_min, range_max;  // Range acceptable pour in
     64
     65  bool rec_gl_samples;        // if true, replace glitch sample values
     66  bool rec_out_range_samples; // if true, replace out of range sample values
     67  bool rec_use_wrec;          // if true, use Mean[Window(wrecsize)] to replace bad samples
     68                              // else use sliding mean value for
     69
    4970};
    5071

  • trunk/ArchTOIPipe/TestPipes/simtst.cc

    r1468 r1478  
     1/*   Test de processeurs ds simtoipr.cc   - Reza Avril 2001
     2
     3----------------   Exemple d'appel  ---------------------
     4csh> simtst -start 104385384 -end 104399964 -range -500,500
     5            -intoi boloMuV_26 -wtoi 8192 -wdegli 1024 inputbolo.fits degli.fits xx.ppf
     6*/
     7
     8
     9
    110#include "machdefs.h"
    211#include <math.h>
Note: See TracChangeset for help on using the changeset viewer.