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

Last change on this file since 1762 was 1762, checked in by aubourg, 24 years ago

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