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

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

correction petit bug flag TOISrc - Reza 15/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.19 2002-05-15 22:44:51 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 uint_8 flg_src = 0;
282 if (k-kgl > maxpoints) flg_src = FlgToiSource; // Si trop long
283 for(ii=kgl; ii<k; ii++)
284 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize)|flg_src);
285
286 lastput = snb+k-1;
287 }
288 }
289 putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
290 lastput = snb+k;
291 kgl = -1; fgglitch = false;
292 vrec(k%wrecsize) = lastvalok = valcur;
293 }
294 else { // Superieur au seuil
295 if (fgglitch) {
296 if (k-kgl+1 > maxpoints) { // serie de points > seuil
297 for(ii=kgl; ii<=k; ii++) // -> Donc pas glitch
298 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
299 lastput = snb+k;
300 fgglitch = false;
301 vrec(k%wrecsize) = lastvalok = valcur;
302 }
303 }
304 else {
305 if (kgl < 0) { // debut possible de glitch
306 fgglitch = true; kgl = k;
307 }
308 else { // On est toujours dans une serie > seuil
309 putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
310 lastput = snb+k; lastvalok = valcur;
311 }
312 vrec(k%wrecsize) = lastvalok;
313 }
314 }
315
316// if (k%5000 == 0) cout << " ---DBG2: K=" << k << " glcount="
317// << glcount << " LastPut= " << lastput << endl;
318 } // Fin de Boucle sur les num-sample
319
320 //DBG cout << " La fin lastput=" << lastput << " SNE=" << sne;
321 //DBG for(k=lastput-snb+1; k<sne-snb; k++)
322 //DBG putData(0, k+snb, vin(k%wsize), 0);
323 // cout << " DBG3- OUT of try bloc ! " << endl;
324 cout << " SimpleDeglitcher::run() - End of processing "
325 << " ProcessedSampleCount=" << ProcessedSampleCount()
326 << " GlitchCount= " << GlitchCount() << endl;
327 if (fgout) deglitchdone = true;
328 } // Bloc try
329 catch (PException & exc) {
330 cerr << "SimpleDeglitcher: Catched Exception " << (string)typeid(exc).name()
331 << "\n .... Msg= " << exc.Msg() << endl;
332 }
333}
334
335
336// -------------------------------------------------------------------
337// Classe SimpleFilter : Filtre simple ds le domaine temporel
338// -------------------------------------------------------------------
339
340string SimpleFilter::FilterKind2String(FilterKind fk)
341{
342 switch (fk) {
343 case UserFilter :
344 return ("UserFilter");
345 break;
346 case MeanFilter :
347 return ("MeanFilter");
348 break;
349 case SumFilter :
350 return ("SumFilter");
351 break;
352 case GaussFilter :
353 return ("GaussFilter");
354 break;
355 case DiffFilter :
356 return ("DiffFilter");
357 break;
358 default :
359 return ("ErrorFilterKind");
360 break;
361 }
362 return("");
363}
364
365SimpleFilter::SimpleFilter(int wsz, FilterKind fk, double a, double s)
366{
367 if (wsz < 3) wsz = 3;
368 if (wsz%2 == 0) wsz++;
369 cout << "SimpleFilter::SimpleFilter() wsz= " << wsz
370 << " FilterKind=" << FilterKind2String(fk) << endl;
371 wsize = wsz;
372 totnscount = 0;
373 coef = new double[wsz];
374 for(int k=0; k<wsz; k++) coef[k] = 0.;
375 int kk;
376 switch (fk) {
377 case UserFilter :
378 throw ParmError("SimpleFilter: Error in filter Kind (UserFilter)!");
379 // break;
380 case MeanFilter :
381 for(kk=0; kk<wsz; kk++)
382 coef[kk] = a/wsize;
383 break;
384 case SumFilter :
385 for(kk=0; kk<wsz; kk++)
386 coef[kk] = a;
387 break;
388 case GaussFilter :
389 for(kk=-(wsz/2); kk<=(wsz/2); kk++)
390 coef[kk+(wsz/2)] = a*exp(-(double)(kk*kk)/(2*s*s));
391 break;
392 case DiffFilter :
393 for(kk=0; kk<wsz; kk++)
394 coef[kk] = -a/wsize;
395 coef[wsz/2+1] += 1.;
396 break;
397 default :
398 throw ParmError("SimpleFilter: Error in filter Kind (UnknownFilter)!");
399 // break;
400 }
401}
402
403SimpleFilter::SimpleFilter(Vector const & vc)
404{
405 int wsz = vc.Size();
406 if (wsz < 3) wsz = 3;
407 if (wsz%2 == 0) wsz++;
408 FilterKind fk = UserFilter;
409 cout << "SimpleFilter::SimpleFilter(Vector & vc) vc.Size()= "
410 << vc.Size() << " WSize=" << wsz
411 << " FilterKind=" << FilterKind2String(fk) << endl;
412 wsize = wsz;
413 totnscount = 0;
414 coef = new double[wsz];
415 int kk;
416 for(kk=0; kk<vc.Size(); kk++)
417 coef[kk] = vc(kk);
418 for(kk=vc.Size(); kk<wsz; kk++)
419 coef[kk] = 0.;
420}
421
422SimpleFilter::~SimpleFilter()
423{
424 delete[] coef;
425}
426
427void SimpleFilter::PrintStatus(::ostream & os)
428{
429 os << "\n ------------------------------------------------------ \n"
430 << " SimpleFilter::PrintStatus() - WindowSize=" << WSize()
431 << " FilterKind= " << Type() << endl;
432 TOIProcessor::PrintStatus(os);
433 os << " Coeff= " ;
434 for(int k=0; k<wsize; k++) os << coef[k] << " " ;
435 os << endl;
436 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
437 os << " ------------------------------------------------------ " << endl;
438}
439
440void SimpleFilter::init() {
441 cout << "SimpleFilter::init" << endl;
442 declareInput("in");
443 declareOutput("out");
444 declareOutput("incopie");
445 name = "SimpleFilter";
446 // upExtra = 1;
447}
448
449void SimpleFilter::run() {
450 // TOIManager* mgr = TOIManager::getManager();
451 int snb = getMinIn();
452 int sne = getMaxIn();
453
454 bool fgout = checkOutputTOIIndex(0);
455 bool fgincopie = checkOutputTOIIndex(1);
456
457 if (!checkInputTOIIndex(0)) {
458 cerr << " SimpleFilter::run() - Input TOI (in) not connected! "
459 << endl;
460 throw ParmError("SimpleFilter::run() Input TOI (in) not connected!");
461 }
462 if (!fgout) {
463 cerr << " SimpleFilter::run() - No Output TOI connected! "
464 << endl;
465 throw ParmError("SimpleFilter::run() No output TOI connected!");
466 }
467
468 cout << " SimpleFilter::run() SNRange=" << snb << " - " << sne << endl;
469
470
471 try {
472 Timer tm("SimpleFilter::run()");
473 // Le debut
474 int wsz2 = wsize/2;
475 Vector vin(wsize);
476 TVector<uint_8> vfg(wsize);
477 int k;
478 for(k=0; k<wsize; k++)
479 getData(0, k+snb, vin(k%wsize), vfg(k%wsize));
480
481 double mean = vin.Sum()/wsize;
482 for(k=wsz2+1; k<wsize; k++) {
483 vin(k) = mean;
484 vfg(k) = 0;
485 }
486 int knext;
487 bool fgfin = false;
488 // Boucle sur les sampleNum
489 for(k=0;k<=sne-snb;k++) {
490 double sortie = 0;
491 for(int ii=-wsz2; ii<=wsz2; ii++) {
492 sortie += vin((ii+k+wsize)%wsize)*coef[ii+wsz2];
493 }
494 putData(0,k+snb,sortie,vfg(k%wsize));
495 if (fgincopie)
496 putData(1, k+snb, vin(k%wsize), vfg(k%wsize));
497 knext = k+wsz2+1;
498 if (knext<=(sne-snb))
499 getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));
500
501 else {
502 if (!fgfin) {
503 mean = vin.Sum()/wsize;
504 fgfin = true;
505 }
506 vin(knext%wsize) = mean;
507 vfg(knext%wsize) = 0;
508 }
509 totnscount++;
510 } // Boucle sur les num-sample
511 cout << " SimpleFilter::run() - End of processing " << endl;
512 } // Bloc try
513
514 catch (PException & exc) {
515 cerr << "SimpleFilter: Catched Exception " << (string)typeid(exc).name()
516 << "\n .... Msg= " << exc.Msg() << endl;
517 }
518}
519
520// ---------------------------------------------------------------
521// -------------------- Classe SimpleAdder -----------------------
522// ---------------------------------------------------------------
523
524SimpleAdder::SimpleAdder(int nbinput)
525 : gains(nbinput)
526{
527 if (nbinput < 1)
528 throw ParmError("SimpleAdder::SimpleAdder() NbInput < 1 !");
529 nb_input = nbinput;
530 for(int k=0; k<nb_input; k++) gains(k) = 1.;
531 totnscount = 0;
532}
533
534SimpleAdder::~SimpleAdder()
535{
536}
537
538void SimpleAdder::SetGain(int num, double g)
539{
540 if ((num < 0) || (num >= nb_input))
541 throw RangeCheckError("SimpleAdder::SetGain() Out of range input number!");
542 gains(num) = g;
543 return;
544}
545
546double SimpleAdder::Gain(int num)
547{
548 if ((num < 0) || (num >= nb_input))
549 throw RangeCheckError("SimpleAdder::Gain() Out of range input number!");
550 return gains(num);
551}
552
553void SimpleAdder::PrintStatus(::ostream & os)
554{
555 os << "\n ------------------------------------------------------ \n"
556 << " SimpleAdder::PrintStatus() - NbInput=" << NbInput() << endl;
557 TOIProcessor::PrintStatus(os);
558 os << " Gains= " ;
559 for(int k=0; k<nb_input; k++) os << gains(k) << " " ;
560 os << endl;
561 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
562 os << " ------------------------------------------------------ " << endl;
563}
564
565void SimpleAdder::init() {
566 cout << "SimpleAdder::init NbInput=" << nb_input << endl;
567 char buff[32];
568 for(int k=0; k<nb_input; k++) {
569 sprintf(buff,"in%d", k);
570 declareInput(buff);
571 }
572
573 declareOutput("out");
574 name = "SimpleAdder";
575 // upExtra = 1;
576}
577
578void SimpleAdder::run() {
579 // TOIManager* mgr = TOIManager::getManager();
580 int snb = getMinIn();
581 int sne = getMaxIn();
582
583 bool fgout = checkOutputTOIIndex(0);
584 string msg_err;
585 for(int ki=0;ki<nb_input;ki++) {
586 if (!checkInputTOIIndex(ki)) {
587 msg_err = "SimpleAdder::run() - Input TOI (" + getInName(ki) +
588 " not connected!";
589 cerr << msg_err << endl;
590 throw ParmError(msg_err);
591 }
592 }
593 if (!fgout) {
594 cerr << " SimpleAdder::run() - No Output TOI connected! "
595 << endl;
596 throw ParmError("SimpleAdder::run() No output TOI connected!");
597 }
598
599 cout << " SimpleAdder::run() SNRange=" << snb << " - " << sne << endl;
600
601
602 try {
603 Timer tm("SimpleAdder::run()");
604 int k,i;
605 double out = 0.;
606 double valin = 0.;
607 uint_8 fgin = 0;
608 uint_8 fgout = 0;
609 for(k=snb;k<=sne;k++) {
610 out = 0;
611 fgout = 0;
612 for(i=0; i<nb_input; i++) {
613 getData(i, k, valin, fgin);
614 out += gains(i)*valin;
615 fgout = fgout | fgin;
616 }
617 putData(0,k,out,fgout);
618 totnscount++;
619 } // Boucle sur les num-sample
620 cout << " SimpleAdder::run() - End of processing " << endl;
621 } // Bloc try
622
623 catch (PException & exc) {
624 cerr << "SimpleAdder: Catched Exception " << (string)typeid(exc).name()
625 << "\n .... Msg= " << exc.Msg() << endl;
626 }
627}
628
629// ----------------------------------------------------------------------
630// Classe SimpleFourierFilter : Filtre simple ds le domaine de Fourier
631// ----------------------------------------------------------------------
632
633SimpleFourierFilter::SimpleFourierFilter(Vector const & vc)
634{
635 ffcoef = vc;
636 wsize = (ffcoef.Size()-1)*2;
637 if (wsize < 16)
638 throw ParmError("SimpleFourierFilter::SimpleFourierFilter() WSize<16 !");
639 KeepSpectra("spectra.ppf", 0);
640 ComputeMeanSpectra(false);
641 totnscount = 0;
642 totnbblock = 0;
643}
644
645SimpleFourierFilter::~SimpleFourierFilter()
646{
647}
648
649
650void SimpleFourierFilter::PrintStatus(::ostream & os)
651{
652 os << "\n ------------------------------------------------------ \n"
653 << " SimpleFourierFilter::PrintStatus() - WindowSize="
654 << WSize() << endl;
655 TOIProcessor::PrintStatus(os);
656 os << " Coeff (Size= " << ffcoef.Size() << "): " << endl;
657 for(int i=0; i<16; i++) {
658 os << ffcoef(i) << " " ;
659 if (i == 7) os << endl;
660 }
661 os << " .... " << endl;
662 os << " ProcessedSampleCount=" << ProcessedSampleCount()
663 << " NbFFTBlocks= " << totnbblock << endl;
664 os << " ------------------------------------------------------ " << endl;
665}
666
667void SimpleFourierFilter::init() {
668 cout << "SimpleFourierFFilter::init" << endl;
669 declareInput("in");
670 declareOutput("out");
671 declareOutput("incopie");
672 name = "SimpleFourierFilter";
673 // upExtra = 1;
674}
675
676
677void SimpleFourierFilter::run() {
678 // TOIManager* mgr = TOIManager::getManager();
679 int snb = getMinIn();
680 int sne = getMaxIn();
681
682 bool fgout = checkOutputTOIIndex(0);
683 bool fgincopie = checkOutputTOIIndex(1);
684
685 if (!checkInputTOIIndex(0)) {
686 cerr << " SimpleFourierFilter::run() - Input TOI (in) not connected! "
687 << endl;
688 throw ParmError("SimpleFourierFilter::run() Input TOI (in) not connected!");
689 }
690
691 // ---- On peut utiliser cette classe pour calculer un spectre de Fourier ----
692 // if (!fgout) {
693 // cerr << " SimpleFourierFilter::run() - No Output TOI connected! "
694 // << endl;
695 // throw ParmError("SimpleFourierFilter::run() No output TOI connected!");
696 // }
697
698 cout << " SimpleFourierFilter::run() SNRange=" << snb << " - " << sne << endl;
699
700
701 try {
702 Timer tm("SimpleFourierFilter::run()");
703 // Le debut
704 Vector vin(wsize);
705 Vector vout(wsize);
706 TVector<uint_8> vfg(wsize);
707 TVector< complex<r_8> > vfft, vfftmean;
708 Vector meanpowerspectra;
709 TVector< complex<r_8> > zcoef(ffcoef.Size());
710 zcoef = ffcoef;
711
712
713 FFTPackServer ffts;
714 ffts.setNormalize(true);
715
716 POutPersist pout(outppfname);
717
718 int k,i,klast;
719 int nks = 0;
720 klast = snb-1;
721 totnbblock = 0;
722 // Boucle sur les sampleNum
723 // 1er partie, on traite par paquets de wsize
724 for(k=snb;k<=sne-wsize+1;k+=wsize) {
725 for(i=0; i<wsize; i++)
726 getData(0, k+i, vin(i), vfg(i));
727 ffts.FFTForward(vin, vfft);
728 if (c_meanspectra) { // Compute mean-spectra
729 if (totnbblock == 0) {
730 vfftmean = vfft;
731 meanpowerspectra.ReSize(vfft.Size());
732 for(i=0; i<meanpowerspectra.Size(); i++)
733 meanpowerspectra(i) = sqrt(vfft(i).real()*vfft(i).real() +
734 vfft(i).imag()*vfft(i).imag() );
735 }
736 else {
737 vfftmean += vfft;
738 for(i=0; i<meanpowerspectra.Size(); i++)
739 meanpowerspectra(i) += sqrt(vfft(i).real()*vfft(i).real() +
740 vfft(i).imag()*vfft(i).imag() );
741 }
742 }
743 totnbblock++;
744 if (nks < nb_keep) {
745 TVector< complex<r_8> > vfftcopie;
746 vfftcopie = vfft;
747 string nomvfft = "spectra" + (string)MuTyV(nks);
748 pout.PutObject(vfftcopie, nomvfft);
749 nks++;
750 }
751 if (fgout) {
752 vfft.MulElt(zcoef);
753 ffts.FFTBackward(vfft, vout);
754 }
755 for(i=0; i<wsize; i++) {
756 if (fgout)
757 putData(0,k+i,vout(i),vfg(i));
758 if (fgincopie)
759 putData(1, k+i, vin(i), vfg(i));
760 }
761 klast+=wsize;
762 totnscount+=wsize;
763 }
764
765
766 // 2eme partie, on traite la fin du bloc d'echantillons si necessaire
767 double inval;
768 uint_8 inflg;
769 if (klast < sne)
770 for(k=klast+1; k<=sne; k++) {
771 getData(0, k, inval, inflg);
772 if (fgout) putData(0, k, inval, inflg);
773 if (fgincopie)
774 putData(1, k, inval, inflg);
775 totnscount++;
776 }
777
778 if (c_meanspectra) {
779 vfftmean /= complex<r_8>((r_8)totnbblock, 0.);
780 pout.PutObject(vfftmean, "meanspectra");
781 meanpowerspectra /= (r_8)totnbblock;
782 pout.PutObject(vfftmean, "meanpowerspectra");
783 }
784 pout.PutObject(ffcoef, "fourierfilter");
785
786 cout << " SimpleFourierFilter::run() - End of processing "
787 << " NbFFTBlocks= " << totnbblock << endl;
788 } // Bloc try
789
790 catch (PException & exc) {
791 cerr << "SimpleFourierFilter: Catched Exception " << (string)typeid(exc).name()
792 << "\n .... Msg= " << exc.Msg() << endl;
793 }
794}
795
796// ---------------------------------------------------------------
797// -------------------- Classe SimpleFanOut -----------------------
798// ---------------------------------------------------------------
799
800SimpleFanOut::SimpleFanOut(int nbinput, int mfanout)
801{
802 if (nbinput < 1)
803 throw ParmError("SimpleFanOut::SimpleFanOut() NbInput < 1 !");
804 if (mfanout < 1)
805 throw ParmError("SimpleFanOut::SimpleFanOut() M_FanOut < 1 !");
806
807 nb_input = nbinput;
808 m_fanout = mfanout;
809 totnscount = 0;
810}
811
812SimpleFanOut::~SimpleFanOut()
813{
814}
815
816
817void SimpleFanOut::PrintStatus(::ostream & os)
818{
819 os << "\n ------------------------------------------------------ \n"
820 << " SimpleFanOut::PrintStatus() - NbInput=" << NbInput()
821 << " M_FanOut=" << MFanOut() << endl;
822 TOIProcessor::PrintStatus(os);
823 os << endl;
824 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
825 os << " ------------------------------------------------------ " << endl;
826}
827
828void SimpleFanOut::init() {
829 cout << "SimpleFanOut::init NbInput=" << nb_input << endl;
830 char buff[64];
831 for(int k=0; k<nb_input; k++) {
832 sprintf(buff,"in%d", k);
833 declareInput(buff);
834 for(int j=0; j<m_fanout; j++) {
835 sprintf(buff,"out%d_%d", k, j);
836 declareOutput(buff);
837 }
838 }
839
840 name = "SimpleFanOut";
841 // upExtra = 1;
842}
843
844void SimpleFanOut::run() {
845 // TOIManager* mgr = TOIManager::getManager();
846 int snb = getMinIn();
847 int sne = getMaxIn();
848
849 TVector<int_4> in_index(nb_input);
850 TMatrix<int_4> out_index(nb_input, m_fanout);
851 in_index = -1;
852 out_index = -1;
853 int nbconin = 0;
854 char buff[64];
855 for(int ki=0;ki<nb_input;ki++) {
856 sprintf(buff,"in%d", ki);
857 int idx = getInputTOIIndex(buff);
858 if (!checkInputTOIIndex(idx)) continue;
859 nbconin++;
860 in_index(ki) = idx;
861 bool fgout = false;
862 for(int jo=0; jo<m_fanout; jo++) {
863 sprintf(buff,"out%d_%d", ki, jo);
864 int odx = getOutputTOIIndex(buff);
865 if (checkOutputTOIIndex(odx)) {
866 out_index(ki, jo) = odx;
867 fgout = true;
868 }
869 }
870 if (!fgout) {
871 string msg_err =
872 "SimpleFanOut::run() - No connected Output for Input TOI ("
873 + getInName(ki) + ") !";
874 cerr << msg_err << endl;
875 throw ParmError(msg_err);
876 }
877 }
878 if (nbconin == 0) {
879 cerr << " SimpleFanOut::run() - No Input TOI connected! "
880 << endl;
881 throw ParmError("SimpleFanOut::run() No Inout TOI connected!");
882 }
883
884 /*
885 for(int ki=0;ki<nb_input;ki++) {
886 cout << " SimpleFanOut::run() In(" << ki << ") Index=" << in_index(ki)
887 << " Name=" << getInName(in_index(ki)) << endl;
888 for(int jo=0; jo<m_fanout; jo++)
889 cout << " .... Out(" << ki << "," << jo << ") Index=" << out_index(ki, jo)
890 << " Name=" << getOutName(out_index(ki, jo)) << endl;
891 }
892 */
893 cout << " SimpleFanOut::run() SNRange=" << snb << " - " << sne << endl;
894
895
896 try {
897 Timer tm("SimpleFanOut::run()");
898 double valin = 0.;
899 uint_8 fgin = 0;
900 for(int k=snb;k<=sne;k++) {
901 for(int i=0;i<nb_input;i++) {
902 if (in_index(i) < 0) continue;
903 valin = 0;
904 fgin = 0;
905 getData(in_index(i), k, valin, fgin);
906 for(int j=0; j<m_fanout; j++) {
907 if (out_index(i, j) < 0) continue;
908 putData(out_index(i, j), k, valin, fgin);
909 }
910 }
911 totnscount++;
912 } // Boucle sur les num-sample
913 cout << " SimpleFanOut::run() - End of processing "
914 << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
915 } // Bloc try
916
917 catch (PException & exc) {
918 cerr << "SimpleFanOut: Catched Exception " << (string)typeid(exc).name()
919 << "\n .... Msg= " << exc.Msg() << endl;
920 }
921}
922
923
924
Note: See TracBrowser for help on using the repository browser.