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

Last change on this file since 1999 was 1999, checked in by ansari, 23 years ago

Ajout flags, correction classe espion RzProcSampleCounter - Reza 13/5/2002
File size: 25.9 KB
Line 
1// ArchTOIPipe (C) CEA/DAPNIA/SPP IN2P3/LAL
2// Eric Aubourg
3// Christophe Magneville
4// Reza Ansari
5// $Id: simtoipr.cc,v 1.18 2002-05-13 15:50:50 ansari Exp $
6
7#include "config.h"
8
9#include "array.h"
10#include "simtoipr.h"
11#include <math.h>
12#include "toimanager.h"
13#include "pexceptions.h"
14#include "ctimer.h"
15#include "fftpserver.h"
16
17#include "flagtoidef.h"
18
19SimpleDeglitcher::SimpleDeglitcher(int wsz, double ns, int maxnpt, int minnpt)
20{
21 SetWSize(wsz);
22 SetDetectionParam(ns, ns/2., maxnpt, minnpt);
23 SetRange(-9.e39, 9.e39);
24 RepBadSamples(true, true);
25
26 totnscount = glnscount = glcount = out_range_nscount = 0;
27 deglitchdone = false;
28
29 cout << "SimpleDeglitcher::SimpleDeglitcher() WSize= " << wsz
30 << " nSig=" << ns << " maxNPt=" << maxnpt << endl;
31}
32
33SimpleDeglitcher::~SimpleDeglitcher()
34{
35}
36
37void SimpleDeglitcher::SetDetectionParam(double ns, double ns2, int maxnpt,
38 int minnpt, int wszrec)
39{
40 nsig = (ns > 0.01) ? ns : 1.;
41 nsig2 = (ns2 > 0.01) ? ns2 : nsig/2.;
42 maxpoints = ((maxnpt > 0) && (maxnpt <= wsize)) ? maxnpt : 5;
43 minpoints = ((minnpt > 0) && (minnpt <= maxpoints)) ? minnpt : maxpoints/2;
44 wrecsize = ((wszrec > 0) && (wszrec <= wsize)) ? wszrec : 2*maxpoints;
45}
46
47inline char * _Bool2YesNo(bool fg)
48{
49 if (fg) return "YES" ;
50 else return "NO" ;
51}
52
53void SimpleDeglitcher::PrintStatus(::ostream & os)
54{
55 os << "\n ------------------------------------------------------ \n"
56 << " SimpleDeglitcher::PrintStatus() - WindowSize=" << WSize()
57 << " MaxPoints=" << MaxPoints() << " MinPoints=" << MinPoints() << endl;
58 os << " NbSigmas=" << NbSigmas()
59 << " NbSigmas2=" << NbSigmas2() << " WRecSize= " << WRecSize()
60 << " Range_Min= " << range_min << " Range_Max= " << range_max << endl;
61 os << " RepOutOfRangeSamples: " << _Bool2YesNo(rec_out_range_samples)
62 << " RepGlitchSamples: " << _Bool2YesNo(rec_gl_samples)
63 << " UseWRec: " << _Bool2YesNo(rec_use_wrec) << endl;
64 TOIProcessor::PrintStatus(os);
65 if (deglitchdone) os << " Deglitching performed " << endl;
66 else os << " NO deglitching done " << endl;
67 double nst = (ProcessedSampleCount() > 0) ? ProcessedSampleCount() : 1.;
68 os << " ProcessedSampleCount=" << ProcessedSampleCount()
69 << " OutOfRangeSampleCount=" << OutOfRangeSampleCount() << endl;
70 os << " GlitchCount= " << GlitchCount()
71 << " GlitchSampleCount=" << GlitchSampleCount()
72 << "( " << (double)GlitchSampleCount()*100./nst << " % )" << endl;
73 os << " ------------------------------------------------------ " << endl;
74}
75
76void SimpleDeglitcher::init() {
77 cout << "SimpleDeglitcher::init" << endl;
78 declareInput("in");
79 declareOutput("out");
80 declareOutput("mean");
81 declareOutput("sigma");
82 declareOutput("incopie");
83 name = "SimpleDeglitcher";
84 // upExtra = 1; A quoi ca sert ?
85}
86
87
88void SimpleDeglitcher::run() {
89
90 // TOIManager* mgr = TOIManager::getManager();
91 int snb = getMinIn();
92 int sne = getMaxIn();
93
94 bool fgout = checkOutputTOIIndex(0);
95 bool fgmean = checkOutputTOIIndex(1);
96 bool fgsigma = checkOutputTOIIndex(2);
97 bool fgincopie = checkOutputTOIIndex(3);
98
99 if (!checkInputTOIIndex(0)) {
100 cerr << " SimpleDeglitcher::run() - Input TOI (in) not connected! "
101 << endl;
102 throw ParmError("SimpleDeglitcher::run() Input TOI (in) not connected!");
103 }
104 if (!fgout && !fgmean && !fgsigma &&!fgincopie) {
105 cerr << " SimpleDeglitcher::run() - No Output TOI connected! "
106 << endl;
107 throw ParmError("SimpleDeglitcher::run() No output TOI connected!");
108 }
109
110 if (!fgout) {
111 cout << "Warning: SimpleDeglitcher::run() - No TOI connected to out \n"
112 << " No deglitching would be performed !" << endl;
113 }
114
115 cout << " SimpleDeglitcher::run() SNRange=" << snb << " - " << sne << endl;
116 try {
117 Timer tm("SimpleDeglitcher::run()");
118 Vector vin(wsize);
119 Vector vas(wsize);
120
121 int wrecsize = maxpoints*2;
122 Vector vrec(wrecsize);
123
124 TVector<uint_8> vfg(wsize);
125 int wsz2 = wsize/2;
126 // Le debut
127 int k;
128 for(k=0; k<wsz2; k++)
129 getData(0, k+snb, vin(k), vfg(k));
130
131 int nokdebut = 0;
132 double s = 0.;
133 double mean = 0.;
134 double s2 = 0.;
135 double sigma = 0.;
136 for(k=0; k<wsz2; k++) {
137 if ( vfg(k) != 0) continue;
138 if ( (vin(k) < range_min) || (vin(k) > range_max) ) continue;
139 s += vin(k);
140 s2 += vin(k)*vin(k);
141 nokdebut++;
142 }
143 if (nokdebut > 0) {
144 mean = s/nokdebut;
145 if (nokdebut > 1) sigma = sqrt(s2/nokdebut-mean*mean);
146 }
147 for(k=wsz2; k<wsize; k++) {
148 vin(k) = mean;
149 vfg(k) = 0;
150 }
151 for(k=0; k<wsize; k++) {
152 vas(k) = mean;
153 if ( vfg(k) != 0) continue;
154 if ( (vin(k) < range_min) || (vin(k) > range_max) ) continue;
155 vas(k) = vin(k);
156 }
157
158 for(k=0; k<wrecsize; k++) {
159 if ( (vin(k) < range_min) || (vin(k) > range_max) ) vrec(k)=mean;
160 else vrec(k)=vin(k);
161 }
162
163 bool fgokdebut = false;
164
165 int kgl = -1;
166 int ii,lastput;
167 bool fgglitch = false;
168 double valcur,valsub,valadd;
169 double lastvalok = mean;
170 uint_8 fgcur;
171 bool fgokcur=false;
172
173 int sx_refresh_count = 0;
174 int sx_refresh_count_max = 16*wsize;
175
176 // Boucle sur les sampleNum
177 int knext;
178 int kfin = sne-snb;
179 for(k=0;k<=kfin;k++) {
180 totnscount++;
181// if (k%10000 == 0) cout << " DBG: K=" << k << endl;
182 knext = k+wsz2;
183 // Calcul mean-sigma
184 if (knext<=kfin) {
185 valsub = vas(knext%wsize);
186 getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));
187 valadd = vin(knext%wsize);
188 if ( vfg(knext%wsize) ||
189 (valadd < range_min) || (valadd > range_max) ) {
190 vas(knext%wsize) = valadd = mean;
191 fgokcur = false;
192 }
193 else {
194 vas(knext%wsize) = valadd = vin(knext%wsize);
195 fgokcur = true;
196 }
197 if (!fgokdebut && fgokcur) {
198 s += valadd;
199 s2 += valadd*valadd;
200 nokdebut++;
201 mean = s/nokdebut;
202 if (nokdebut > 1) sigma = sqrt(s2/nokdebut-mean*mean);
203 if (nokdebut >= wsize) {
204 fgokdebut = true;
205 cout << " SimpleDeglitcher::DebugInfo - nokdebut=" << nokdebut
206 << " k=" << k << " knext=" << knext
207 << "\n ...DebugInfo mean=" << mean
208 << " sigma=" << sigma << " s=" << s << " s2=" << s2 << endl;
209 }
210 }
211 else {
212 if (sx_refresh_count >= sx_refresh_count_max) {
213 // On recalcule la somme
214 s = vas.Sum();
215 s2 = vas.SumX2();
216 sx_refresh_count = 0;
217 }
218 else {
219 s += (valadd-valsub);
220 s2 += (valadd*valadd-valsub*valsub);
221 sx_refresh_count++;
222 }
223 mean = s/wsize;
224 sigma = sqrt(s2/wsize-mean*mean);
225 }
226 }
227
228
229 // On gere les sorties Mean et Sigma
230 if (fgmean)
231 putData(1, k+snb, mean, 0);
232 if (fgsigma)
233 putData(2, k+snb, sigma, 0);
234 if (fgincopie)
235 putData(3, k+snb, vin(k%wsize), vfg(k%wsize));
236
237 valcur = vin(k%wsize);
238 if ( (valcur < range_min) || (valcur > range_max) ) {
239 valcur = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
240 if (rec_out_range_samples) vin(k%wsize) = valcur;
241 vfg(k%wsize) |= FlgToiOut;
242 out_range_nscount++;
243 }
244 valcur = vin(k%wsize);
245 fgcur = vfg(k%wsize);
246
247 if (!fgout) continue; // Pas de sortie out (deglitche)
248
249
250// if (k<100) {
251// cout << "DBG-A-Deglitch[" << k << "] mean="
252// << mean << " sigma=" << sigma << " valcur="
253// << valcur << " Kgl=" << kgl ;
254// if (fgglitch) cout << " In Glitch" ;
255// cout << endl;
256// }
257
258 double curnsig = (fgglitch) ? nsig2 : nsig;
259
260 if (valcur < mean+curnsig*sigma) { // inferieur au seuil
261 if (fgglitch) {
262 if ( (k-kgl <= maxpoints) && (k-kgl >= minpoints) ) {
263 // On vient de detecter un glitch
264 glcount++;
265 if (rec_gl_samples) { // On change la valeur des samples
266 double recval = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
267 for(ii=kgl; ii<k; ii++) {
268 putData(0, ii+snb, recval, vfg(ii%wsize)|FlgToiSpike);
269 glnscount++;
270 }
271 }
272 else { // On ne fait que flagger les echantillons
273 for(ii=kgl; ii<k; ii++) {
274 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize)|FlgToiSpike);
275 glnscount++;
276 }
277 }
278 lastput = snb+k-1;
279 } // - Fin de detection de glitch
280 else { // Trop long ou trop court - ce n'est pas un glitch ...
281 for(ii=kgl; ii<k; ii++) {
282 if (k-kgl > maxpoints) // Si trop long
283 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize)|FlgToiSource);
284 else // Trop court
285 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
286 }
287 lastput = snb+k-1;
288 }
289 }
290 putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
291 lastput = snb+k;
292 kgl = -1; fgglitch = false;
293 vrec(k%wrecsize) = lastvalok = valcur;
294 }
295 else { // Superieur au seuil
296 if (fgglitch) {
297 if (k-kgl+1 > maxpoints) { // serie de points > seuil
298 for(ii=kgl; ii<=k; ii++) // -> Donc pas glitch
299 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
300 lastput = snb+k;
301 fgglitch = false;
302 vrec(k%wrecsize) = lastvalok = valcur;
303 }
304 }
305 else {
306 if (kgl < 0) { // debut possible de glitch
307 fgglitch = true; kgl = k;
308 }
309 else { // On est toujours dans une serie > seuil
310 putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
311 lastput = snb+k; lastvalok = valcur;
312 }
313 vrec(k%wrecsize) = lastvalok;
314 }
315 }
316
317// if (k%5000 == 0) cout << " ---DBG2: K=" << k << " glcount="
318// << glcount << " LastPut= " << lastput << endl;
319 } // Fin de Boucle sur les num-sample
320
321 //DBG cout << " La fin lastput=" << lastput << " SNE=" << sne;
322 //DBG for(k=lastput-snb+1; k<sne-snb; k++)
323 //DBG putData(0, k+snb, vin(k%wsize), 0);
324 // cout << " DBG3- OUT of try bloc ! " << endl;
325 cout << " SimpleDeglitcher::run() - End of processing "
326 << " ProcessedSampleCount=" << ProcessedSampleCount()
327 << " GlitchCount= " << GlitchCount() << endl;
328 if (fgout) deglitchdone = true;
329 } // Bloc try
330 catch (PException & exc) {
331 cerr << "SimpleDeglitcher: Catched Exception " << (string)typeid(exc).name()
332 << "\n .... Msg= " << exc.Msg() << endl;
333 }
334}
335
336
337// -------------------------------------------------------------------
338// Classe SimpleFilter : Filtre simple ds le domaine temporel
339// -------------------------------------------------------------------
340
341string SimpleFilter::FilterKind2String(FilterKind fk)
342{
343 switch (fk) {
344 case UserFilter :
345 return ("UserFilter");
346 break;
347 case MeanFilter :
348 return ("MeanFilter");
349 break;
350 case SumFilter :
351 return ("SumFilter");
352 break;
353 case GaussFilter :
354 return ("GaussFilter");
355 break;
356 case DiffFilter :
357 return ("DiffFilter");
358 break;
359 default :
360 return ("ErrorFilterKind");
361 break;
362 }
363 return("");
364}
365
366SimpleFilter::SimpleFilter(int wsz, FilterKind fk, double a, double s)
367{
368 if (wsz < 3) wsz = 3;
369 if (wsz%2 == 0) wsz++;
370 cout << "SimpleFilter::SimpleFilter() wsz= " << wsz
371 << " FilterKind=" << FilterKind2String(fk) << endl;
372 wsize = wsz;
373 totnscount = 0;
374 coef = new double[wsz];
375 for(int k=0; k<wsz; k++) coef[k] = 0.;
376 int kk;
377 switch (fk) {
378 case UserFilter :
379 throw ParmError("SimpleFilter: Error in filter Kind (UserFilter)!");
380 // break;
381 case MeanFilter :
382 for(kk=0; kk<wsz; kk++)
383 coef[kk] = a/wsize;
384 break;
385 case SumFilter :
386 for(kk=0; kk<wsz; kk++)
387 coef[kk] = a;
388 break;
389 case GaussFilter :
390 for(kk=-(wsz/2); kk<=(wsz/2); kk++)
391 coef[kk+(wsz/2)] = a*exp(-(double)(kk*kk)/(2*s*s));
392 break;
393 case DiffFilter :
394 for(kk=0; kk<wsz; kk++)
395 coef[kk] = -a/wsize;
396 coef[wsz/2+1] += 1.;
397 break;
398 default :
399 throw ParmError("SimpleFilter: Error in filter Kind (UnknownFilter)!");
400 // break;
401 }
402}
403
404SimpleFilter::SimpleFilter(Vector const & vc)
405{
406 int wsz = vc.Size();
407 if (wsz < 3) wsz = 3;
408 if (wsz%2 == 0) wsz++;
409 FilterKind fk = UserFilter;
410 cout << "SimpleFilter::SimpleFilter(Vector & vc) vc.Size()= "
411 << vc.Size() << " WSize=" << wsz
412 << " FilterKind=" << FilterKind2String(fk) << endl;
413 wsize = wsz;
414 totnscount = 0;
415 coef = new double[wsz];
416 int kk;
417 for(kk=0; kk<vc.Size(); kk++)
418 coef[kk] = vc(kk);
419 for(kk=vc.Size(); kk<wsz; kk++)
420 coef[kk] = 0.;
421}
422
423SimpleFilter::~SimpleFilter()
424{
425 delete[] coef;
426}
427
428void SimpleFilter::PrintStatus(::ostream & os)
429{
430 os << "\n ------------------------------------------------------ \n"
431 << " SimpleFilter::PrintStatus() - WindowSize=" << WSize()
432 << " FilterKind= " << Type() << endl;
433 TOIProcessor::PrintStatus(os);
434 os << " Coeff= " ;
435 for(int k=0; k<wsize; k++) os << coef[k] << " " ;
436 os << endl;
437 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
438 os << " ------------------------------------------------------ " << endl;
439}
440
441void SimpleFilter::init() {
442 cout << "SimpleFilter::init" << endl;
443 declareInput("in");
444 declareOutput("out");
445 declareOutput("incopie");
446 name = "SimpleFilter";
447 // upExtra = 1;
448}
449
450void SimpleFilter::run() {
451 // TOIManager* mgr = TOIManager::getManager();
452 int snb = getMinIn();
453 int sne = getMaxIn();
454
455 bool fgout = checkOutputTOIIndex(0);
456 bool fgincopie = checkOutputTOIIndex(1);
457
458 if (!checkInputTOIIndex(0)) {
459 cerr << " SimpleFilter::run() - Input TOI (in) not connected! "
460 << endl;
461 throw ParmError("SimpleFilter::run() Input TOI (in) not connected!");
462 }
463 if (!fgout) {
464 cerr << " SimpleFilter::run() - No Output TOI connected! "
465 << endl;
466 throw ParmError("SimpleFilter::run() No output TOI connected!");
467 }
468
469 cout << " SimpleFilter::run() SNRange=" << snb << " - " << sne << endl;
470
471
472 try {
473 Timer tm("SimpleFilter::run()");
474 // Le debut
475 int wsz2 = wsize/2;
476 Vector vin(wsize);
477 TVector<uint_8> vfg(wsize);
478 int k;
479 for(k=0; k<wsize; k++)
480 getData(0, k+snb, vin(k%wsize), vfg(k%wsize));
481
482 double mean = vin.Sum()/wsize;
483 for(k=wsz2+1; k<wsize; k++) {
484 vin(k) = mean;
485 vfg(k) = 0;
486 }
487 int knext;
488 bool fgfin = false;
489 // Boucle sur les sampleNum
490 for(k=0;k<=sne-snb;k++) {
491 double sortie = 0;
492 for(int ii=-wsz2; ii<=wsz2; ii++) {
493 sortie += vin((ii+k+wsize)%wsize)*coef[ii+wsz2];
494 }
495 putData(0,k+snb,sortie,vfg(k%wsize));
496 if (fgincopie)
497 putData(1, k+snb, vin(k%wsize), vfg(k%wsize));
498 knext = k+wsz2+1;
499 if (knext<=(sne-snb))
500 getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));
501
502 else {
503 if (!fgfin) {
504 mean = vin.Sum()/wsize;
505 fgfin = true;
506 }
507 vin(knext%wsize) = mean;
508 vfg(knext%wsize) = 0;
509 }
510 totnscount++;
511 } // Boucle sur les num-sample
512 cout << " SimpleFilter::run() - End of processing " << endl;
513 } // Bloc try
514
515 catch (PException & exc) {
516 cerr << "SimpleFilter: Catched Exception " << (string)typeid(exc).name()
517 << "\n .... Msg= " << exc.Msg() << endl;
518 }
519}
520
521// ---------------------------------------------------------------
522// -------------------- Classe SimpleAdder -----------------------
523// ---------------------------------------------------------------
524
525SimpleAdder::SimpleAdder(int nbinput)
526 : gains(nbinput)
527{
528 if (nbinput < 1)
529 throw ParmError("SimpleAdder::SimpleAdder() NbInput < 1 !");
530 nb_input = nbinput;
531 for(int k=0; k<nb_input; k++) gains(k) = 1.;
532 totnscount = 0;
533}
534
535SimpleAdder::~SimpleAdder()
536{
537}
538
539void SimpleAdder::SetGain(int num, double g)
540{
541 if ((num < 0) || (num >= nb_input))
542 throw RangeCheckError("SimpleAdder::SetGain() Out of range input number!");
543 gains(num) = g;
544 return;
545}
546
547double SimpleAdder::Gain(int num)
548{
549 if ((num < 0) || (num >= nb_input))
550 throw RangeCheckError("SimpleAdder::Gain() Out of range input number!");
551 return gains(num);
552}
553
554void SimpleAdder::PrintStatus(::ostream & os)
555{
556 os << "\n ------------------------------------------------------ \n"
557 << " SimpleAdder::PrintStatus() - NbInput=" << NbInput() << endl;
558 TOIProcessor::PrintStatus(os);
559 os << " Gains= " ;
560 for(int k=0; k<nb_input; k++) os << gains(k) << " " ;
561 os << endl;
562 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
563 os << " ------------------------------------------------------ " << endl;
564}
565
566void SimpleAdder::init() {
567 cout << "SimpleAdder::init NbInput=" << nb_input << endl;
568 char buff[32];
569 for(int k=0; k<nb_input; k++) {
570 sprintf(buff,"in%d", k);
571 declareInput(buff);
572 }
573
574 declareOutput("out");
575 name = "SimpleAdder";
576 // upExtra = 1;
577}
578
579void SimpleAdder::run() {
580 // TOIManager* mgr = TOIManager::getManager();
581 int snb = getMinIn();
582 int sne = getMaxIn();
583
584 bool fgout = checkOutputTOIIndex(0);
585 string msg_err;
586 for(int ki=0;ki<nb_input;ki++) {
587 if (!checkInputTOIIndex(ki)) {
588 msg_err = "SimpleAdder::run() - Input TOI (" + getInName(ki) +
589 " not connected!";
590 cerr << msg_err << endl;
591 throw ParmError(msg_err);
592 }
593 }
594 if (!fgout) {
595 cerr << " SimpleAdder::run() - No Output TOI connected! "
596 << endl;
597 throw ParmError("SimpleAdder::run() No output TOI connected!");
598 }
599
600 cout << " SimpleAdder::run() SNRange=" << snb << " - " << sne << endl;
601
602
603 try {
604 Timer tm("SimpleAdder::run()");
605 int k,i;
606 double out = 0.;
607 double valin = 0.;
608 uint_8 fgin = 0;
609 uint_8 fgout = 0;
610 for(k=snb;k<=sne;k++) {
611 out = 0;
612 fgout = 0;
613 for(i=0; i<nb_input; i++) {
614 getData(i, k, valin, fgin);
615 out += gains(i)*valin;
616 fgout = fgout | fgin;
617 }
618 putData(0,k,out,fgout);
619 totnscount++;
620 } // Boucle sur les num-sample
621 cout << " SimpleAdder::run() - End of processing " << endl;
622 } // Bloc try
623
624 catch (PException & exc) {
625 cerr << "SimpleAdder: Catched Exception " << (string)typeid(exc).name()
626 << "\n .... Msg= " << exc.Msg() << endl;
627 }
628}
629
630// ----------------------------------------------------------------------
631// Classe SimpleFourierFilter : Filtre simple ds le domaine de Fourier
632// ----------------------------------------------------------------------
633
634SimpleFourierFilter::SimpleFourierFilter(Vector const & vc)
635{
636 ffcoef = vc;
637 wsize = (ffcoef.Size()-1)*2;
638 if (wsize < 16)
639 throw ParmError("SimpleFourierFilter::SimpleFourierFilter() WSize<16 !");
640 KeepSpectra("spectra.ppf", 0);
641 ComputeMeanSpectra(false);
642 totnscount = 0;
643 totnbblock = 0;
644}
645
646SimpleFourierFilter::~SimpleFourierFilter()
647{
648}
649
650
651void SimpleFourierFilter::PrintStatus(::ostream & os)
652{
653 os << "\n ------------------------------------------------------ \n"
654 << " SimpleFourierFilter::PrintStatus() - WindowSize="
655 << WSize() << endl;
656 TOIProcessor::PrintStatus(os);
657 os << " Coeff (Size= " << ffcoef.Size() << "): " << endl;
658 for(int i=0; i<16; i++) {
659 os << ffcoef(i) << " " ;
660 if (i == 7) os << endl;
661 }
662 os << " .... " << endl;
663 os << " ProcessedSampleCount=" << ProcessedSampleCount()
664 << " NbFFTBlocks= " << totnbblock << endl;
665 os << " ------------------------------------------------------ " << endl;
666}
667
668void SimpleFourierFilter::init() {
669 cout << "SimpleFourierFFilter::init" << endl;
670 declareInput("in");
671 declareOutput("out");
672 declareOutput("incopie");
673 name = "SimpleFourierFilter";
674 // upExtra = 1;
675}
676
677
678void SimpleFourierFilter::run() {
679 // TOIManager* mgr = TOIManager::getManager();
680 int snb = getMinIn();
681 int sne = getMaxIn();
682
683 bool fgout = checkOutputTOIIndex(0);
684 bool fgincopie = checkOutputTOIIndex(1);
685
686 if (!checkInputTOIIndex(0)) {
687 cerr << " SimpleFourierFilter::run() - Input TOI (in) not connected! "
688 << endl;
689 throw ParmError("SimpleFourierFilter::run() Input TOI (in) not connected!");
690 }
691
692 // ---- On peut utiliser cette classe pour calculer un spectre de Fourier ----
693 // if (!fgout) {
694 // cerr << " SimpleFourierFilter::run() - No Output TOI connected! "
695 // << endl;
696 // throw ParmError("SimpleFourierFilter::run() No output TOI connected!");
697 // }
698
699 cout << " SimpleFourierFilter::run() SNRange=" << snb << " - " << sne << endl;
700
701
702 try {
703 Timer tm("SimpleFourierFilter::run()");
704 // Le debut
705 Vector vin(wsize);
706 Vector vout(wsize);
707 TVector<uint_8> vfg(wsize);
708 TVector< complex<r_8> > vfft, vfftmean;
709 Vector meanpowerspectra;
710 TVector< complex<r_8> > zcoef(ffcoef.Size());
711 zcoef = ffcoef;
712
713
714 FFTPackServer ffts;
715 ffts.setNormalize(true);
716
717 POutPersist pout(outppfname);
718
719 int k,i,klast;
720 int nks = 0;
721 klast = snb-1;
722 totnbblock = 0;
723 // Boucle sur les sampleNum
724 // 1er partie, on traite par paquets de wsize
725 for(k=snb;k<=sne-wsize+1;k+=wsize) {
726 for(i=0; i<wsize; i++)
727 getData(0, k+i, vin(i), vfg(i));
728 ffts.FFTForward(vin, vfft);
729 if (c_meanspectra) { // Compute mean-spectra
730 if (totnbblock == 0) {
731 vfftmean = vfft;
732 meanpowerspectra.ReSize(vfft.Size());
733 for(i=0; i<meanpowerspectra.Size(); i++)
734 meanpowerspectra(i) = sqrt(vfft(i).real()*vfft(i).real() +
735 vfft(i).imag()*vfft(i).imag() );
736 }
737 else {
738 vfftmean += vfft;
739 for(i=0; i<meanpowerspectra.Size(); i++)
740 meanpowerspectra(i) += sqrt(vfft(i).real()*vfft(i).real() +
741 vfft(i).imag()*vfft(i).imag() );
742 }
743 }
744 totnbblock++;
745 if (nks < nb_keep) {
746 TVector< complex<r_8> > vfftcopie;
747 vfftcopie = vfft;
748 string nomvfft = "spectra" + (string)MuTyV(nks);
749 pout.PutObject(vfftcopie, nomvfft);
750 nks++;
751 }
752 if (fgout) {
753 vfft.MulElt(zcoef);
754 ffts.FFTBackward(vfft, vout);
755 }
756 for(i=0; i<wsize; i++) {
757 if (fgout)
758 putData(0,k+i,vout(i),vfg(i));
759 if (fgincopie)
760 putData(1, k+i, vin(i), vfg(i));
761 }
762 klast+=wsize;
763 totnscount+=wsize;
764 }
765
766
767 // 2eme partie, on traite la fin du bloc d'echantillons si necessaire
768 double inval;
769 uint_8 inflg;
770 if (klast < sne)
771 for(k=klast+1; k<=sne; k++) {
772 getData(0, k, inval, inflg);
773 if (fgout) putData(0, k, inval, inflg);
774 if (fgincopie)
775 putData(1, k, inval, inflg);
776 totnscount++;
777 }
778
779 if (c_meanspectra) {
780 vfftmean /= complex<r_8>((r_8)totnbblock, 0.);
781 pout.PutObject(vfftmean, "meanspectra");
782 meanpowerspectra /= (r_8)totnbblock;
783 pout.PutObject(vfftmean, "meanpowerspectra");
784 }
785 pout.PutObject(ffcoef, "fourierfilter");
786
787 cout << " SimpleFourierFilter::run() - End of processing "
788 << " NbFFTBlocks= " << totnbblock << endl;
789 } // Bloc try
790
791 catch (PException & exc) {
792 cerr << "SimpleFourierFilter: Catched Exception " << (string)typeid(exc).name()
793 << "\n .... Msg= " << exc.Msg() << endl;
794 }
795}
796
797// ---------------------------------------------------------------
798// -------------------- Classe SimpleFanOut -----------------------
799// ---------------------------------------------------------------
800
801SimpleFanOut::SimpleFanOut(int nbinput, int mfanout)
802{
803 if (nbinput < 1)
804 throw ParmError("SimpleFanOut::SimpleFanOut() NbInput < 1 !");
805 if (mfanout < 1)
806 throw ParmError("SimpleFanOut::SimpleFanOut() M_FanOut < 1 !");
807
808 nb_input = nbinput;
809 m_fanout = mfanout;
810 totnscount = 0;
811}
812
813SimpleFanOut::~SimpleFanOut()
814{
815}
816
817
818void SimpleFanOut::PrintStatus(::ostream & os)
819{
820 os << "\n ------------------------------------------------------ \n"
821 << " SimpleFanOut::PrintStatus() - NbInput=" << NbInput()
822 << " M_FanOut=" << MFanOut() << endl;
823 TOIProcessor::PrintStatus(os);
824 os << endl;
825 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
826 os << " ------------------------------------------------------ " << endl;
827}
828
829void SimpleFanOut::init() {
830 cout << "SimpleFanOut::init NbInput=" << nb_input << endl;
831 char buff[64];
832 for(int k=0; k<nb_input; k++) {
833 sprintf(buff,"in%d", k);
834 declareInput(buff);
835 for(int j=0; j<m_fanout; j++) {
836 sprintf(buff,"out%d_%d", k, j);
837 declareOutput(buff);
838 }
839 }
840
841 name = "SimpleFanOut";
842 // upExtra = 1;
843}
844
845void SimpleFanOut::run() {
846 // TOIManager* mgr = TOIManager::getManager();
847 int snb = getMinIn();
848 int sne = getMaxIn();
849
850 TVector<int_4> in_index(nb_input);
851 TMatrix<int_4> out_index(nb_input, m_fanout);
852 in_index = -1;
853 out_index = -1;
854 int nbconin = 0;
855 char buff[64];
856 for(int ki=0;ki<nb_input;ki++) {
857 sprintf(buff,"in%d", ki);
858 int idx = getInputTOIIndex(buff);
859 if (!checkInputTOIIndex(idx)) continue;
860 nbconin++;
861 in_index(ki) = idx;
862 bool fgout = false;
863 for(int jo=0; jo<m_fanout; jo++) {
864 sprintf(buff,"out%d_%d", ki, jo);
865 int odx = getOutputTOIIndex(buff);
866 if (checkOutputTOIIndex(odx)) {
867 out_index(ki, jo) = odx;
868 fgout = true;
869 }
870 }
871 if (!fgout) {
872 string msg_err =
873 "SimpleFanOut::run() - No connected Output for Input TOI ("
874 + getInName(ki) + ") !";
875 cerr << msg_err << endl;
876 throw ParmError(msg_err);
877 }
878 }
879 if (nbconin == 0) {
880 cerr << " SimpleFanOut::run() - No Input TOI connected! "
881 << endl;
882 throw ParmError("SimpleFanOut::run() No Inout TOI connected!");
883 }
884
885 /*
886 for(int ki=0;ki<nb_input;ki++) {
887 cout << " SimpleFanOut::run() In(" << ki << ") Index=" << in_index(ki)
888 << " Name=" << getInName(in_index(ki)) << endl;
889 for(int jo=0; jo<m_fanout; jo++)
890 cout << " .... Out(" << ki << "," << jo << ") Index=" << out_index(ki, jo)
891 << " Name=" << getOutName(out_index(ki, jo)) << endl;
892 }
893 */
894 cout << " SimpleFanOut::run() SNRange=" << snb << " - " << sne << endl;
895
896
897 try {
898 Timer tm("SimpleFanOut::run()");
899 double valin = 0.;
900 uint_8 fgin = 0;
901 for(int k=snb;k<=sne;k++) {
902 for(int i=0;i<nb_input;i++) {
903 if (in_index(i) < 0) continue;
904 valin = 0;
905 fgin = 0;
906 getData(in_index(i), k, valin, fgin);
907 for(int j=0; j<m_fanout; j++) {
908 if (out_index(i, j) < 0) continue;
909 putData(out_index(i, j), k, valin, fgin);
910 }
911 }
912 totnscount++;
913 } // Boucle sur les num-sample
914 cout << " SimpleFanOut::run() - End of processing "
915 << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
916 } // Bloc try
917
918 catch (PException & exc) {
919 cerr << "SimpleFanOut: Catched Exception " << (string)typeid(exc).name()
920 << "\n .... Msg= " << exc.Msg() << endl;
921 }
922}
923
924
925
Note: See TracBrowser for help on using the repository browser.