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

Last change on this file since 1462 was 1462, checked in by cmv, 24 years ago

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