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

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

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