source: Sophya/trunk/ArchTOIPipe/ProcWSophya/simtoipr.cc@ 1454

Last change on this file since 1454 was 1454, checked in by ansari, 24 years ago

Modifs effectuees a St Pierre de Charetreuse - Reza 9/4/2001
File size: 14.0 KB
Line 
1#include "simtoipr.h"
2#include <math.h>
3#include "toimanager.h"
4#include "pexceptions.h"
5#include "ctimer.h"
6
7SimpleDeglitcher::SimpleDeglitcher(int wsz, double ns, int maxnpt)
8{
9 if (wsz < 5) wsz = 5;
10 if (wsz%2 == 0) wsz++;
11 if (maxnpt > wsz) maxnpt = wsz;
12
13 cout << "SimpleDeglitcher::SimpleDeglitcher() WSize= " << wsz
14 << " nSig=" << ns << " maxNPt=" << maxnpt << endl;
15
16 wsize = wsz;
17 nsig = ns;
18 maxpoints = maxnpt;
19 totnscount = glnscount = glcount = out_range_nscount = 0;
20 deglitchdone = false;
21 SetRange(-9.e39, 9.e39);
22}
23
24SimpleDeglitcher::~SimpleDeglitcher()
25{
26}
27
28
29void SimpleDeglitcher::PrintStatus(ostream & os)
30{
31 os << "\n ------------------------------------------------------ \n"
32 << " SimpleDeglitcher::PrintStatus() - WindowSize=" << WSize()
33 << " NbSigmas=" << NbSigmas() << " MaxPoints=" << MaxPoints() << endl;
34 os << " Range_Min= " << range_min << " Range_Max= " << range_max << endl;
35 TOIProcessor::PrintStatus(os);
36 if (deglitchdone) os << " Deglitching performed " << endl;
37 else os << " NO deglitching done " << endl;
38 double nst = (ProcessedSampleCount() > 0) ? ProcessedSampleCount() : 1.;
39 os << " ProcessedSampleCount=" << ProcessedSampleCount()
40 << " OutOfRangeSampleCount=" << OutOfRangeSampleCount() << endl;
41 os << " GlitchCount= " << GlitchCount()
42 << " GlitchSampleCount=" << GlitchSampleCount()
43 << "( " << (double)GlitchSampleCount()*100./nst << " % )" << endl;
44 os << " ------------------------------------------------------ " << endl;
45}
46
47void SimpleDeglitcher::init() {
48 cout << "SimpleDeglitcher::init" << endl;
49 declareInput("in");
50 declareOutput("out");
51 declareOutput("mean");
52 declareOutput("sigma");
53 declareOutput("incopie");
54 name = "SimpleDeglitcher";
55 // upExtra = 1; A quoi ca sert ?
56}
57
58void SimpleDeglitcher::run() {
59
60 // TOIManager* mgr = TOIManager::getManager();
61 int snb = getMinIn();
62 int sne = getMaxIn();
63
64 bool fgout = checkOutputTOIIndex(0);
65 bool fgmean = checkOutputTOIIndex(1);
66 bool fgsigma = checkOutputTOIIndex(2);
67 bool fgincopie = checkOutputTOIIndex(3);
68
69 if (!checkInputTOIIndex(0)) {
70 cerr << " SimpleDeglitcher::run() - Input TOI (in) not connected! "
71 << endl;
72 throw ParmError("SimpleDeglitcher::run() Input TOI (in) not connected!");
73 }
74 if (!fgout && !fgmean && !fgsigma &&!fgincopie) {
75 cerr << " SimpleDeglitcher::run() - No Output TOI connected! "
76 << endl;
77 throw ParmError("SimpleDeglitcher::run() No output TOI connected!");
78 }
79
80 if (!fgout) {
81 cout << "Warning: SimpleDeglitcher::run() - No TOI connected to out \n"
82 << " No deglitching would be performed !" << endl;
83 }
84
85 cout << " SimpleDeglitcher::run() SNRange=" << snb << " - " << sne << endl;
86 try {
87 Timer tm("SimpleDeglitcher::run()");
88 // Le debut
89 Vector vin(wsize);
90 TVector<int_8> vfg(wsize);
91 int k;
92 for(k=0; k<wsize; k++) {
93 vin(k) = getData(0, k+snb);
94 vfg(k) = getFlag(0, k+snb);
95 }
96 double s = vin.Sum();
97 double mean = s/wsize;
98 for(k=0; k<wsize; k++) {
99 if ( (vin(k) < range_min) || (vin(k) > range_max) ) {
100 vin(k) = mean;
101 vfg(k) |= 2;
102 out_range_nscount++;
103 }
104 }
105 s = vin.Sum();
106 mean = s/wsize;
107
108 double s2 = vin.SumX2();
109 double sigma = sqrt(s2/wsize-mean*mean);
110 int kgl = -1;
111 int ii,lastput;
112 bool fgglitch = false;
113 double valcur,valsub,valadd;
114 uint_8 fgcur;
115
116 // Boucle sur les sampleNum
117
118 for(k=0;k<=sne-snb;k++) {
119 totnscount++;
120// if (k%10000 == 0) cout << " DBG: K=" << k << endl;
121 // Calcul mean-sigma
122 if (k>=wsize) {
123 valsub = vin(k%wsize);
124 vin(k%wsize) = valadd = getData(0, k+snb);
125 vfg(k%wsize) = getFlag(0, k+snb);
126 if ( (valadd < range_min) || (valadd > range_max) )
127 valadd = mean;
128 s += (valadd-valsub);
129 s2 += (valadd*valadd-valsub*valsub);
130 mean = s/wsize;
131 sigma = sqrt(s2/wsize-mean*mean);
132 }
133
134 // On gere les sorties Mean et Sigma
135 if (fgmean)
136 putData(1, k+snb, mean, 0);
137 if (fgsigma)
138 putData(2, k+snb, sigma, 0);
139 if (fgincopie)
140 putData(3, k+snb, vin(k%wsize), vfg(k%wsize));
141
142 valcur = vin(k%wsize);
143 if ( (valcur < range_min) || (valcur > range_max) ) {
144 vin(k%wsize) = valcur = mean;
145 vfg(k%wsize) |= 2;
146 out_range_nscount++;
147 }
148
149 if (!fgout) continue; // Pas de sortie out (deglitche)
150
151
152// if (k<100) {
153// cout << "DBG-A-Deglitch[" << k << "] mean="
154// << mean << " sigma=" << sigma << " valcur="
155// << valcur << " Kgl=" << kgl ;
156// if (fgglitch) cout << " In Glitch" ;
157// cout << endl;
158// }
159
160 if (valcur < mean+nsig*sigma) { // inferieur au seuil
161 if (fgglitch) {
162 if (k-kgl < maxpoints) { // On vient de detecter un glitch
163 glcount++;
164 for(ii=kgl; ii<k; ii++) {
165 putData(0, ii+snb, mean, vfg(ii%wsize)|5);
166 glnscount++;
167 }
168 lastput = snb+k-1;
169 }
170 else {
171 for(ii=kgl; ii<k; ii++) {
172 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
173 }
174 lastput = snb+k-1;
175 }
176 }
177 putData(0, k+snb, valcur, 0);
178 lastput = snb+k;
179 kgl = -1; fgglitch = false;
180 }
181 else { // Superieur au seuil
182 if (fgglitch) {
183 if (k-kgl+1 >= maxpoints) { // serie de points > seuil
184 for(ii=kgl; ii<=k; ii++) // -> Donc pas glitch
185 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
186 lastput = snb+k;
187 fgglitch = false;
188 }
189 }
190 else {
191 if (kgl < 0) { // debut possible de glitch
192 fgglitch = true; kgl = k;
193 }
194 else { // On est toujours dans une serie > seuil
195 putData(0, k+snb, valcur, fgcur);
196 lastput = snb+k;
197 }
198 }
199 }
200
201// if (k%5000 == 0) cout << " ---DBG2: K=" << k << " glcount="
202// << glcount << " LastPut= " << lastput << endl;
203 } // Fin de Boucle sur les num-sample
204
205 //DBG cout << " La fin lastput=" << lastput << " SNE=" << sne;
206 //DBG for(k=lastput-snb+1; k<sne-snb; k++)
207 //DBG putData(0, k+snb, vin(k%wsize), 0);
208 // cout << " DBG3- OUT of try bloc ! " << endl;
209 cout << " SimpleDeglitcher::run() - End of processing "
210 << " ProcessedSampleCount=" << ProcessedSampleCount()
211 << " GlitchCount= " << GlitchCount() << endl;
212 if (fgout) deglitchdone = true;
213 } // Bloc try
214 catch (PException & exc) {
215 cerr << "SimpleDeglitcher: Catched Exception " << (string)typeid(exc).name()
216 << "\n .... Msg= " << exc.Msg() << endl;
217 }
218}
219
220
221// -------------------------------------------------------------------
222// Classe SimpleFilter : Filtre simple ds le domaine temporel
223// -------------------------------------------------------------------
224
225string SimpleFilter::FilterKind2String(FilterKind fk)
226{
227 switch (fk) {
228 case UserFilter :
229 return ("UserFilter");
230 break;
231 case MeanFilter :
232 return ("MeanFilter");
233 break;
234 case SumFilter :
235 return ("SumFilter");
236 break;
237 case GaussFilter :
238 return ("GaussFilter");
239 break;
240 case DiffFilter :
241 return ("DiffFilter");
242 break;
243 default :
244 return ("ErrorFilterKind");
245 break;
246 }
247 return("");
248}
249
250SimpleFilter::SimpleFilter(int wsz, FilterKind fk, double a, double s)
251{
252 if (wsz < 3) wsz = 3;
253 if (wsz%2 == 0) wsz++;
254 cout << "SimpleFilter::SimpleFilter() wsz= " << wsz
255 << " FilterKind=" << FilterKind2String(fk) << endl;
256 wsize = wsz;
257 totnscount = 0;
258 coef = new double[wsz];
259 switch (fk) {
260 case UserFilter :
261 throw ParmError("SimpleFilter: Error in filter Kind (UserFilter)!");
262 break;
263 case MeanFilter :
264 for(int kk=0; kk<wsz; kk++)
265 coef[kk] = a/wsize;
266 break;
267 case SumFilter :
268 for(int kk=0; kk<wsz; kk++)
269 coef[kk] = a;
270 break;
271 case GaussFilter :
272 for(int kk=-(wsz/2); kk<=(wsz/2); kk++)
273 coef[kk] = a*exp((double)(kk*kk)/(s*s));
274 break;
275 case DiffFilter :
276 for(int kk=0; kk<wsz; kk++)
277 coef[kk] = -a/wsize;
278 coef[wsz/2+1] += 1.;
279 break;
280 default :
281 throw ParmError("SimpleFilter: Error in filter Kind (UnknownFilter)!");
282 break;
283 }
284}
285
286SimpleFilter::SimpleFilter(Vector const & vc)
287{
288 int wsz = vc.Size();
289 if (wsz < 3) wsz = 3;
290 if (wsz%2 == 0) wsz++;
291 FilterKind fk = UserFilter;
292 cout << "SimpleFilter::SimpleFilter(Vector & vc) vc.Size()= "
293 << vc.Size() << " WSize=" << wsz
294 << " FilterKind=" << FilterKind2String(fk) << endl;
295 wsize = wsz;
296 totnscount = 0;
297 coef = new double[wsz];
298 for(int kk=0; kk<vc.Size(); kk++)
299 coef[kk] = vc(kk);
300 for(int kk=vc.Size(); kk<wsz; kk++)
301 coef[kk] = 0.;
302}
303
304SimpleFilter::~SimpleFilter()
305{
306 delete[] coef;
307}
308
309void SimpleFilter::PrintStatus(ostream & os)
310{
311 os << "\n ------------------------------------------------------ \n"
312 << " SimpleFilter::PrintStatus() - WindowSize=" << WSize()
313 << " FilterKind= " << Type() << endl;
314 TOIProcessor::PrintStatus(os);
315 os << " Coeff= " ;
316 for(int k=0; k<wsize; k++) os << coef[k] << " " ;
317 os << endl;
318 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
319 os << " ------------------------------------------------------ " << endl;
320}
321
322void SimpleFilter::init() {
323 cout << "SimpleFilter::init" << endl;
324 declareInput("in");
325 declareOutput("out");
326 declareOutput("incopie");
327 name = "SimpleFilter";
328 // upExtra = 1;
329}
330
331void SimpleFilter::run() {
332 // TOIManager* mgr = TOIManager::getManager();
333 int snb = getMinIn();
334 int sne = getMaxIn();
335
336 bool fgout = checkOutputTOIIndex(0);
337 bool fgincopie = checkOutputTOIIndex(1);
338
339 if (!checkInputTOIIndex(0)) {
340 cerr << " SimpleFilter::run() - Input TOI (in) not connected! "
341 << endl;
342 throw ParmError("SimpleFilter::run() Input TOI (in) not connected!");
343 }
344 if (!fgout) {
345 cerr << " SimpleFilter::run() - No Output TOI connected! "
346 << endl;
347 throw ParmError("SimpleFilter::run() No output TOI connected!");
348 }
349
350 cout << " SimpleFilter::run() SNRange=" << snb << " - " << sne << endl;
351
352
353 try {
354 Timer tm("SimpleFilter::run()");
355 // Le debut
356 int wsz2 = wsize/2;
357 Vector vin(wsize);
358 TVector<int_8> vfg(wsize);
359 int k;
360 for(k=0; k<wsize; k++) {
361 vin(k%wsize) = getData(0, k+snb);
362 vfg(k%wsize) = getFlag(0, k+snb);
363 }
364 double mean = vin.Sum()/wsize;
365 for(k=wsz2+1; k<wsize; k++) {
366 vin(k) = mean;
367 vfg(k) = 0;
368 }
369 int knext;
370 bool fgfin = false;
371 // Boucle sur les sampleNum
372 for(k=0;k<=sne-snb;k++) {
373 double sortie = 0;
374 for(int ii=0; ii<wsize; ii++) {
375 sortie += vin(ii)*coef[(ii+wsz2)%wsize];
376 }
377 putData(0,k+snb,sortie,vfg(k%wsize));
378 if (fgincopie)
379 putData(1, k+snb, vin(k%wsize), vfg(k%wsize));
380 knext = k+wsz2+1;
381 if (knext<=(sne-snb)) {
382 vin(knext%wsize) = getData(0, knext+snb);
383 vfg(knext%wsize) = getFlag(0, knext+snb);
384 }
385 else {
386 if (!fgfin) {
387 mean = vin.Sum()/wsize;
388 fgfin = true;
389 }
390 vin(knext%wsize) = mean;
391 vfg(knext%wsize) = 0;
392 }
393 totnscount++;
394 } // Boucle sur les num-sample
395 cout << " SimpleFilter::run() - End of processing " << endl;
396 } // Bloc try
397
398 catch (PException & exc) {
399 cerr << "SimpleFilter: Catched Exception " << (string)typeid(exc).name()
400 << "\n .... Msg= " << exc.Msg() << endl;
401 }
402}
403
404// ---------------------------------------------------------------
405// -------------------- Classe SimpleAdder -----------------------
406// ---------------------------------------------------------------
407
408SimpleAdder::SimpleAdder(int nbinput)
409 : gains(nbinput)
410{
411 if (nbinput < 1)
412 throw ParmError("SimpleAdder::SimpleAdder() NbInput < 1 !");
413 nb_input = nbinput;
414 for(int k=0; k<nb_input; k++) gains(k) = 1.;
415 totnscount = 0;
416}
417
418SimpleAdder::~SimpleAdder()
419{
420}
421
422void SimpleAdder::SetGain(int num, double g)
423{
424 if ((num < 0) || (num >= nb_input))
425 throw RangeCheckError("SimpleAdder::SetGain() Out of range input number!");
426 gains(num) = g;
427 return;
428}
429
430double SimpleAdder::Gain(int num)
431{
432 if ((num < 0) || (num >= nb_input))
433 throw RangeCheckError("SimpleAdder::Gain() Out of range input number!");
434 return gains(num);
435}
436
437void SimpleAdder::PrintStatus(ostream & os)
438{
439 os << "\n ------------------------------------------------------ \n"
440 << " SimpleAdder::PrintStatus() - NbInput=" << NbInput() << endl;
441 TOIProcessor::PrintStatus(os);
442 os << " Gains= " ;
443 for(int k=0; k<nb_input; k++) os << gains(k) << " " ;
444 os << endl;
445 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
446 os << " ------------------------------------------------------ " << endl;
447}
448
449void SimpleAdder::init() {
450 cout << "SimpleAdder::init NbInput=" << nb_input << endl;
451 char buff[32];
452 for(int k=0; k<nb_input; k++) {
453 sprintf(buff,"in%d", k);
454 declareInput(buff);
455 }
456
457 declareOutput("out");
458 name = "SimpleAdder";
459 // upExtra = 1;
460}
461
462void SimpleAdder::run() {
463 // TOIManager* mgr = TOIManager::getManager();
464 int snb = getMinIn();
465 int sne = getMaxIn();
466
467 bool fgout = checkOutputTOIIndex(0);
468 string msg_err;
469 for(int ki=0;ki<nb_input;ki++) {
470 if (!checkInputTOIIndex(ki)) {
471 msg_err = "SimpleAdder::run() - Input TOI (" + getInName(ki) +
472 " not connected!";
473 cerr << msg_err << endl;
474 throw ParmError(msg_err);
475 }
476 }
477 if (!fgout) {
478 cerr << " SimpleAdder::run() - No Output TOI connected! "
479 << endl;
480 throw ParmError("SimpleAdder::run() No output TOI connected!");
481 }
482
483 cout << " SimpleAdder::run() SNRange=" << snb << " - " << sne << endl;
484
485
486 try {
487 Timer tm("SimpleAdder::run()");
488 int k,i;
489 double out = 0.;
490 unsigned long fgout = 0;
491 for(k=snb;k<=sne;k++) {
492 out = 0;
493 fgout = 0;
494 for(i=0; i<nb_input; i++) {
495 out += gains(i)*getData(i, k);
496 fgout = fgout | (unsigned int)getFlag(i,k);
497 }
498 putData(0,k,out,fgout);
499 totnscount++;
500 } // Boucle sur les num-sample
501 cout << " SimpleAdder::run() - End of processing " << endl;
502 } // Bloc try
503
504 catch (PException & exc) {
505 cerr << "SimpleAdder: Catched Exception " << (string)typeid(exc).name()
506 << "\n .... Msg= " << exc.Msg() << endl;
507 }
508}
509
510
511
Note: See TracBrowser for help on using the repository browser.