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

Last change on this file since 1483 was 1483, checked in by ansari, 24 years ago

Petites amelioration SimpleFourierFilter - Reza 30/4/2001
File size: 24.2 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<int_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 // Boucle sur les sampleNum
165
166 int knext;
167 int kfin = sne-snb;
168 for(k=0;k<=kfin;k++) {
169 totnscount++;
170// if (k%10000 == 0) cout << " DBG: K=" << k << endl;
171 knext = k+wsz2;
172 // Calcul mean-sigma
173 if (knext<=kfin) {
174 valsub = vas(knext%wsize);
175 getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));
176 valadd = vin(knext%wsize);
177 if ( vfg(knext%wsize) ||
178 (valadd < range_min) || (valadd > range_max) ) {
179 vas(knext%wsize) = valadd = mean;
180 fgokcur = false;
181 }
182 else {
183 vas(knext%wsize) = valadd = vin(knext%wsize);
184 fgokcur = true;
185 }
186 if (!fgokdebut && fgokcur) {
187 s += valadd;
188 s2 += valadd*valadd;
189 nokdebut++;
190 mean = s/nokdebut;
191 if (nokdebut > 1) sigma = sqrt(s2/nokdebut-mean*mean);
192 if (nokdebut >= wsize) {
193 fgokdebut = true;
194 cout << " SimpleDeglitcher::DebugInfo - nokdebut=" << nokdebut
195 << " k=" << k << " knext=" << knext
196 << "\n ...DebugInfo mean=" << mean
197 << " sigma=" << sigma << " s=" << s << " s2=" << s2 << endl;
198 }
199 }
200 else {
201 s += (valadd-valsub);
202 s2 += (valadd*valadd-valsub*valsub);
203 mean = s/wsize;
204 sigma = sqrt(s2/wsize-mean*mean);
205 }
206 }
207
208
209 // On gere les sorties Mean et Sigma
210 if (fgmean)
211 putData(1, k+snb, mean, 0);
212 if (fgsigma)
213 putData(2, k+snb, sigma, 0);
214 if (fgincopie)
215 putData(3, k+snb, vin(k%wsize), vfg(k%wsize));
216
217 valcur = vin(k%wsize);
218 if ( (valcur < range_min) || (valcur > range_max) ) {
219 valcur = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
220 if (rec_out_range_samples) vin(k%wsize) = valcur;
221 vfg(k%wsize) |= FG_OUTOFRANGE;
222 out_range_nscount++;
223 }
224 valcur = vin(k%wsize);
225 fgcur = vfg(k%wsize);
226
227 if (!fgout) continue; // Pas de sortie out (deglitche)
228
229
230// if (k<100) {
231// cout << "DBG-A-Deglitch[" << k << "] mean="
232// << mean << " sigma=" << sigma << " valcur="
233// << valcur << " Kgl=" << kgl ;
234// if (fgglitch) cout << " In Glitch" ;
235// cout << endl;
236// }
237
238 double curnsig = (fgglitch) ? nsig2 : nsig;
239
240 if (valcur < mean+curnsig*sigma) { // inferieur au seuil
241 if (fgglitch) {
242 if ( (k-kgl <= maxpoints) && (k-kgl >= minpoints) ) {
243 // On vient de detecter un glitch
244 glcount++;
245 if (rec_gl_samples) { // On change la valeur des samples
246 double recval = (rec_use_wrec) ? vrec.Sum()/wrecsize : mean;
247 for(ii=kgl; ii<k; ii++) {
248 putData(0, ii+snb, recval, vfg(ii%wsize)|FG_GLITCH);
249 glnscount++;
250 }
251 }
252 else { // On ne fait que flagger les echantillons
253 for(ii=kgl; ii<k; ii++) {
254 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize)|FG_GLITCH);
255 glnscount++;
256 }
257 }
258 lastput = snb+k-1;
259 } // - Fin de detection de glitch
260 else { // Trop long ou trop court - ce n'est pas un glitch ...
261 for(ii=kgl; ii<k; ii++) {
262 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
263 }
264 lastput = snb+k-1;
265 }
266 }
267 putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
268 lastput = snb+k;
269 kgl = -1; fgglitch = false;
270 vrec(k%wrecsize) = lastvalok = valcur;
271 }
272 else { // Superieur au seuil
273 if (fgglitch) {
274 if (k-kgl+1 > maxpoints) { // serie de points > seuil
275 for(ii=kgl; ii<=k; ii++) // -> Donc pas glitch
276 putData(0, ii+snb, vin(ii%wsize), vfg(ii%wsize));
277 lastput = snb+k;
278 fgglitch = false;
279 vrec(k%wrecsize) = lastvalok = valcur;
280 }
281 }
282 else {
283 if (kgl < 0) { // debut possible de glitch
284 fgglitch = true; kgl = k;
285 }
286 else { // On est toujours dans une serie > seuil
287 putData(0, k+snb, vin(k%wsize), vfg(k%wsize));
288 lastput = snb+k; lastvalok = valcur;
289 }
290 vrec(k%wrecsize) = lastvalok;
291 }
292 }
293
294// if (k%5000 == 0) cout << " ---DBG2: K=" << k << " glcount="
295// << glcount << " LastPut= " << lastput << endl;
296 } // Fin de Boucle sur les num-sample
297
298 //DBG cout << " La fin lastput=" << lastput << " SNE=" << sne;
299 //DBG for(k=lastput-snb+1; k<sne-snb; k++)
300 //DBG putData(0, k+snb, vin(k%wsize), 0);
301 // cout << " DBG3- OUT of try bloc ! " << endl;
302 cout << " SimpleDeglitcher::run() - End of processing "
303 << " ProcessedSampleCount=" << ProcessedSampleCount()
304 << " GlitchCount= " << GlitchCount() << endl;
305 if (fgout) deglitchdone = true;
306 } // Bloc try
307 catch (PException & exc) {
308 cerr << "SimpleDeglitcher: Catched Exception " << (string)typeid(exc).name()
309 << "\n .... Msg= " << exc.Msg() << endl;
310 }
311}
312
313
314// -------------------------------------------------------------------
315// Classe SimpleFilter : Filtre simple ds le domaine temporel
316// -------------------------------------------------------------------
317
318string SimpleFilter::FilterKind2String(FilterKind fk)
319{
320 switch (fk) {
321 case UserFilter :
322 return ("UserFilter");
323 break;
324 case MeanFilter :
325 return ("MeanFilter");
326 break;
327 case SumFilter :
328 return ("SumFilter");
329 break;
330 case GaussFilter :
331 return ("GaussFilter");
332 break;
333 case DiffFilter :
334 return ("DiffFilter");
335 break;
336 default :
337 return ("ErrorFilterKind");
338 break;
339 }
340 return("");
341}
342
343SimpleFilter::SimpleFilter(int wsz, FilterKind fk, double a, double s)
344{
345 if (wsz < 3) wsz = 3;
346 if (wsz%2 == 0) wsz++;
347 cout << "SimpleFilter::SimpleFilter() wsz= " << wsz
348 << " FilterKind=" << FilterKind2String(fk) << endl;
349 wsize = wsz;
350 totnscount = 0;
351 coef = new double[wsz];
352 for(int k=0; k<wsz; k++) coef[k] = 0.;
353 int kk;
354 switch (fk) {
355 case UserFilter :
356 throw ParmError("SimpleFilter: Error in filter Kind (UserFilter)!");
357 // break;
358 case MeanFilter :
359 for(kk=0; kk<wsz; kk++)
360 coef[kk] = a/wsize;
361 break;
362 case SumFilter :
363 for(kk=0; kk<wsz; kk++)
364 coef[kk] = a;
365 break;
366 case GaussFilter :
367 for(kk=-(wsz/2); kk<=(wsz/2); kk++)
368 coef[kk+(wsz/2)] = a*exp(-(double)(kk*kk)/(s*s));
369 break;
370 case DiffFilter :
371 for(kk=0; kk<wsz; kk++)
372 coef[kk] = -a/wsize;
373 coef[wsz/2+1] += 1.;
374 break;
375 default :
376 throw ParmError("SimpleFilter: Error in filter Kind (UnknownFilter)!");
377 // break;
378 }
379}
380
381SimpleFilter::SimpleFilter(Vector const & vc)
382{
383 int wsz = vc.Size();
384 if (wsz < 3) wsz = 3;
385 if (wsz%2 == 0) wsz++;
386 FilterKind fk = UserFilter;
387 cout << "SimpleFilter::SimpleFilter(Vector & vc) vc.Size()= "
388 << vc.Size() << " WSize=" << wsz
389 << " FilterKind=" << FilterKind2String(fk) << endl;
390 wsize = wsz;
391 totnscount = 0;
392 coef = new double[wsz];
393 int kk;
394 for(kk=0; kk<vc.Size(); kk++)
395 coef[kk] = vc(kk);
396 for(kk=vc.Size(); kk<wsz; kk++)
397 coef[kk] = 0.;
398}
399
400SimpleFilter::~SimpleFilter()
401{
402 delete[] coef;
403}
404
405void SimpleFilter::PrintStatus(ostream & os)
406{
407 os << "\n ------------------------------------------------------ \n"
408 << " SimpleFilter::PrintStatus() - WindowSize=" << WSize()
409 << " FilterKind= " << Type() << endl;
410 TOIProcessor::PrintStatus(os);
411 os << " Coeff= " ;
412 for(int k=0; k<wsize; k++) os << coef[k] << " " ;
413 os << endl;
414 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
415 os << " ------------------------------------------------------ " << endl;
416}
417
418void SimpleFilter::init() {
419 cout << "SimpleFilter::init" << endl;
420 declareInput("in");
421 declareOutput("out");
422 declareOutput("incopie");
423 name = "SimpleFilter";
424 // upExtra = 1;
425}
426
427void SimpleFilter::run() {
428 // TOIManager* mgr = TOIManager::getManager();
429 int snb = getMinIn();
430 int sne = getMaxIn();
431
432 bool fgout = checkOutputTOIIndex(0);
433 bool fgincopie = checkOutputTOIIndex(1);
434
435 if (!checkInputTOIIndex(0)) {
436 cerr << " SimpleFilter::run() - Input TOI (in) not connected! "
437 << endl;
438 throw ParmError("SimpleFilter::run() Input TOI (in) not connected!");
439 }
440 if (!fgout) {
441 cerr << " SimpleFilter::run() - No Output TOI connected! "
442 << endl;
443 throw ParmError("SimpleFilter::run() No output TOI connected!");
444 }
445
446 cout << " SimpleFilter::run() SNRange=" << snb << " - " << sne << endl;
447
448
449 try {
450 Timer tm("SimpleFilter::run()");
451 // Le debut
452 int wsz2 = wsize/2;
453 Vector vin(wsize);
454 TVector<int_8> vfg(wsize);
455 int k;
456 for(k=0; k<wsize; k++)
457 getData(0, k+snb, vin(k%wsize), vfg(k%wsize));
458
459 double mean = vin.Sum()/wsize;
460 for(k=wsz2+1; k<wsize; k++) {
461 vin(k) = mean;
462 vfg(k) = 0;
463 }
464 int knext;
465 bool fgfin = false;
466 // Boucle sur les sampleNum
467 for(k=0;k<=sne-snb;k++) {
468 double sortie = 0;
469 for(int ii=-wsz2; ii<=wsz2; ii++) {
470 sortie += vin((ii+k+wsize)%wsize)*coef[ii+wsz2];
471 }
472 putData(0,k+snb,sortie,vfg(k%wsize));
473 if (fgincopie)
474 putData(1, k+snb, vin(k%wsize), vfg(k%wsize));
475 knext = k+wsz2+1;
476 if (knext<=(sne-snb))
477 getData(0, knext+snb, vin(knext%wsize), vfg(knext%wsize));
478
479 else {
480 if (!fgfin) {
481 mean = vin.Sum()/wsize;
482 fgfin = true;
483 }
484 vin(knext%wsize) = mean;
485 vfg(knext%wsize) = 0;
486 }
487 totnscount++;
488 } // Boucle sur les num-sample
489 cout << " SimpleFilter::run() - End of processing " << endl;
490 } // Bloc try
491
492 catch (PException & exc) {
493 cerr << "SimpleFilter: Catched Exception " << (string)typeid(exc).name()
494 << "\n .... Msg= " << exc.Msg() << endl;
495 }
496}
497
498// ---------------------------------------------------------------
499// -------------------- Classe SimpleAdder -----------------------
500// ---------------------------------------------------------------
501
502SimpleAdder::SimpleAdder(int nbinput)
503 : gains(nbinput)
504{
505 if (nbinput < 1)
506 throw ParmError("SimpleAdder::SimpleAdder() NbInput < 1 !");
507 nb_input = nbinput;
508 for(int k=0; k<nb_input; k++) gains(k) = 1.;
509 totnscount = 0;
510}
511
512SimpleAdder::~SimpleAdder()
513{
514}
515
516void SimpleAdder::SetGain(int num, double g)
517{
518 if ((num < 0) || (num >= nb_input))
519 throw RangeCheckError("SimpleAdder::SetGain() Out of range input number!");
520 gains(num) = g;
521 return;
522}
523
524double SimpleAdder::Gain(int num)
525{
526 if ((num < 0) || (num >= nb_input))
527 throw RangeCheckError("SimpleAdder::Gain() Out of range input number!");
528 return gains(num);
529}
530
531void SimpleAdder::PrintStatus(ostream & os)
532{
533 os << "\n ------------------------------------------------------ \n"
534 << " SimpleAdder::PrintStatus() - NbInput=" << NbInput() << endl;
535 TOIProcessor::PrintStatus(os);
536 os << " Gains= " ;
537 for(int k=0; k<nb_input; k++) os << gains(k) << " " ;
538 os << endl;
539 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
540 os << " ------------------------------------------------------ " << endl;
541}
542
543void SimpleAdder::init() {
544 cout << "SimpleAdder::init NbInput=" << nb_input << endl;
545 char buff[32];
546 for(int k=0; k<nb_input; k++) {
547 sprintf(buff,"in%d", k);
548 declareInput(buff);
549 }
550
551 declareOutput("out");
552 name = "SimpleAdder";
553 // upExtra = 1;
554}
555
556void SimpleAdder::run() {
557 // TOIManager* mgr = TOIManager::getManager();
558 int snb = getMinIn();
559 int sne = getMaxIn();
560
561 bool fgout = checkOutputTOIIndex(0);
562 string msg_err;
563 for(int ki=0;ki<nb_input;ki++) {
564 if (!checkInputTOIIndex(ki)) {
565 msg_err = "SimpleAdder::run() - Input TOI (" + getInName(ki) +
566 " not connected!";
567 cerr << msg_err << endl;
568 throw ParmError(msg_err);
569 }
570 }
571 if (!fgout) {
572 cerr << " SimpleAdder::run() - No Output TOI connected! "
573 << endl;
574 throw ParmError("SimpleAdder::run() No output TOI connected!");
575 }
576
577 cout << " SimpleAdder::run() SNRange=" << snb << " - " << sne << endl;
578
579
580 try {
581 Timer tm("SimpleAdder::run()");
582 int k,i;
583 double out = 0.;
584 double valin = 0.;
585 int_8 fgin = 0;
586 int_8 fgout = 0;
587 for(k=snb;k<=sne;k++) {
588 out = 0;
589 fgout = 0;
590 for(i=0; i<nb_input; i++) {
591 getData(i, k, valin, fgin);
592 out += gains(i)*valin;
593 fgout = fgout | fgin;
594 }
595 putData(0,k,out,fgout);
596 totnscount++;
597 } // Boucle sur les num-sample
598 cout << " SimpleAdder::run() - End of processing " << endl;
599 } // Bloc try
600
601 catch (PException & exc) {
602 cerr << "SimpleAdder: Catched Exception " << (string)typeid(exc).name()
603 << "\n .... Msg= " << exc.Msg() << endl;
604 }
605}
606
607// ----------------------------------------------------------------------
608// Classe SimpleFourierFilter : Filtre simple ds le domaine de Fourier
609// ----------------------------------------------------------------------
610
611SimpleFourierFilter::SimpleFourierFilter(Vector const & vc)
612{
613 ffcoef = vc;
614 wsize = (ffcoef.Size()-1)*2;
615 if (wsize < 16)
616 throw ParmError("SimpleFourierFilter::SimpleFourierFilter() WSize<16 !");
617 KeepSpectra("spectra.ppf", 0);
618}
619
620SimpleFourierFilter::~SimpleFourierFilter()
621{
622}
623
624
625void SimpleFourierFilter::PrintStatus(ostream & os)
626{
627 os << "\n ------------------------------------------------------ \n"
628 << " SimpleFourierFilter::PrintStatus() - WindowSize="
629 << WSize() << endl;
630 TOIProcessor::PrintStatus(os);
631 os << " Coeff= " << ffcoef << endl;
632 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
633 os << " ------------------------------------------------------ " << endl;
634}
635
636void SimpleFourierFilter::init() {
637 cout << "SimpleFourierFFilter::init" << endl;
638 declareInput("in");
639 declareOutput("out");
640 declareOutput("incopie");
641 name = "SimpleFourierFilter";
642 // upExtra = 1;
643}
644
645
646void SimpleFourierFilter::run() {
647 // TOIManager* mgr = TOIManager::getManager();
648 int snb = getMinIn();
649 int sne = getMaxIn();
650
651 bool fgout = checkOutputTOIIndex(0);
652 bool fgincopie = checkOutputTOIIndex(1);
653
654 if (!checkInputTOIIndex(0)) {
655 cerr << " SimpleFourierFilter::run() - Input TOI (in) not connected! "
656 << endl;
657 throw ParmError("SimpleFourierFilter::run() Input TOI (in) not connected!");
658 }
659 if (!fgout) {
660 cerr << " SimpleFourierFilter::run() - No Output TOI connected! "
661 << endl;
662 throw ParmError("SimpleFourierFilter::run() No output TOI connected!");
663 }
664
665 cout << " SimpleFourierFilter::run() SNRange=" << snb << " - " << sne << endl;
666
667
668 try {
669 Timer tm("SimpleFourierFilter::run()");
670 // Le debut
671 Vector vin(wsize);
672 Vector vout(wsize);
673 TVector<int_8> vfg(wsize);
674 TVector< complex<r_8> > vfft, vfftmean;
675 TVector< complex<r_8> > zcoef(ffcoef.Size());
676 zcoef = ffcoef;
677
678
679 FFTPackServer ffts;
680 ffts.setNormalize(true);
681
682 POutPersist pout(outppfname);
683
684 int k,i,klast;
685 int nks = 0;
686 int nblk = 0;
687 klast = snb-1;
688 // Boucle sur les sampleNum
689 // 1er partie, on traite par paquets de wsize
690 for(k=snb;k<=sne-wsize+1;k+=wsize) {
691 for(i=0; i<wsize; i++)
692 getData(0, k+i, vin(i), vfg(i));
693 ffts.FFTForward(vin, vfft);
694 if (nblk == 0) vfftmean = vfft;
695 else vfftmean += vfft;
696 nblk++;
697 if (nks < nb_keep) {
698 TVector< complex<r_8> > vfftcopie;
699 vfftcopie = vfft;
700 string nomvfft = "spectra" + MuTyV(nks);
701 pout.PutObject(vfftcopie, nomvfft);
702 nks++;
703 }
704 vfft.MulElt(zcoef);
705 ffts.FFTBackward(vfft, vout);
706 for(i=0; i<wsize; i++)
707 putData(0,k+i,vout(i),vfg(i));
708 if (fgincopie)
709 for(i=0; i<wsize; i++)
710 putData(1, k+i, vin(i), vfg(i));
711 klast+=wsize;
712 totnscount+=wsize;
713 }
714
715
716 // 2eme partie, on traite la fin du bloc d'echantillons si necessaire
717 double inval;
718 int_8 inflg;
719 if (klast < sne)
720 for(k=klast+1; k<=sne; k++) {
721 getData(0, k, inval, inflg);
722 putData(0, k, inval, inflg);
723 if (fgincopie)
724 putData(1, k, inval, inflg);
725 totnscount++;
726 }
727
728
729 vfftmean /= complex<r_8>((r_8)nblk, 0.);
730 pout.PutObject(vfftmean, "meanspectra");
731
732 cout << " SimpleFourierFilter::run() - End of processing " << endl;
733 } // Bloc try
734
735 catch (PException & exc) {
736 cerr << "SimpleFourierFilter: Catched Exception " << (string)typeid(exc).name()
737 << "\n .... Msg= " << exc.Msg() << endl;
738 }
739}
740
741// ---------------------------------------------------------------
742// -------------------- Classe SimpleFanOut -----------------------
743// ---------------------------------------------------------------
744
745SimpleFanOut::SimpleFanOut(int nbinput, int mfanout)
746{
747 if (nbinput < 1)
748 throw ParmError("SimpleFanOut::SimpleFanOut() NbInput < 1 !");
749 if (mfanout < 1)
750 throw ParmError("SimpleFanOut::SimpleFanOut() M_FanOut < 1 !");
751
752 nb_input = nbinput;
753 m_fanout = mfanout;
754 totnscount = 0;
755}
756
757SimpleFanOut::~SimpleFanOut()
758{
759}
760
761
762void SimpleFanOut::PrintStatus(ostream & os)
763{
764 os << "\n ------------------------------------------------------ \n"
765 << " SimpleFanOut::PrintStatus() - NbInput=" << NbInput()
766 << " M_FanOut=" << MFanOut() << endl;
767 TOIProcessor::PrintStatus(os);
768 os << endl;
769 os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
770 os << " ------------------------------------------------------ " << endl;
771}
772
773void SimpleFanOut::init() {
774 cout << "SimpleFanOut::init NbInput=" << nb_input << endl;
775 char buff[64];
776 for(int k=0; k<nb_input; k++) {
777 sprintf(buff,"in%d", k);
778 declareInput(buff);
779 for(int j=0; j<m_fanout; j++) {
780 sprintf(buff,"out%d_%d", k, j);
781 declareOutput(buff);
782 }
783 }
784
785 name = "SimpleFanOut";
786 // upExtra = 1;
787}
788
789void SimpleFanOut::run() {
790 // TOIManager* mgr = TOIManager::getManager();
791 int snb = getMinIn();
792 int sne = getMaxIn();
793
794 TVector<int_4> in_index(nb_input);
795 TMatrix<int_4> out_index(nb_input, m_fanout);
796 in_index = -1;
797 out_index = -1;
798 int nbconin = 0;
799 char buff[64];
800 for(int ki=0;ki<nb_input;ki++) {
801 sprintf(buff,"in%d", ki);
802 int idx = getInputTOIIndex(buff);
803 if (!checkInputTOIIndex(idx)) continue;
804 nbconin++;
805 in_index(ki) = idx;
806 bool fgout = false;
807 for(int jo=0; jo<m_fanout; jo++) {
808 sprintf(buff,"out%d_%d", ki, jo);
809 int odx = getOutputTOIIndex(buff);
810 if (checkOutputTOIIndex(odx)) {
811 out_index(ki, jo) = odx;
812 fgout = true;
813 }
814 }
815 if (!fgout) {
816 string msg_err =
817 "SimpleFanOut::run() - No connected Output for Input TOI ("
818 + getInName(ki) + ") !";
819 cerr << msg_err << endl;
820 throw ParmError(msg_err);
821 }
822 }
823 if (nbconin == 0) {
824 cerr << " SimpleFanOut::run() - No Input TOI connected! "
825 << endl;
826 throw ParmError("SimpleFanOut::run() No Inout TOI connected!");
827 }
828
829 /*
830 for(int ki=0;ki<nb_input;ki++) {
831 cout << " SimpleFanOut::run() In(" << ki << ") Index=" << in_index(ki)
832 << " Name=" << getInName(in_index(ki)) << endl;
833 for(int jo=0; jo<m_fanout; jo++)
834 cout << " .... Out(" << ki << "," << jo << ") Index=" << out_index(ki, jo)
835 << " Name=" << getOutName(out_index(ki, jo)) << endl;
836 }
837 */
838 cout << " SimpleFanOut::run() SNRange=" << snb << " - " << sne << endl;
839
840
841 try {
842 Timer tm("SimpleFanOut::run()");
843 double valin = 0.;
844 int_8 fgin = 0;
845 for(int k=snb;k<=sne;k++) {
846 for(int i=0;i<nb_input;i++) {
847 if (in_index(i) < 0) continue;
848 valin = 0;
849 fgin = 0;
850 getData(in_index(i), k, valin, fgin);
851 for(int j=0; j<m_fanout; j++) {
852 if (out_index(i, j) < 0) continue;
853 putData(out_index(i, j), k, valin, fgin);
854 }
855 }
856 totnscount++;
857 } // Boucle sur les num-sample
858 cout << " SimpleFanOut::run() - End of processing "
859 << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
860 } // Bloc try
861
862 catch (PException & exc) {
863 cerr << "SimpleFanOut: Catched Exception " << (string)typeid(exc).name()
864 << "\n .... Msg= " << exc.Msg() << endl;
865 }
866}
867
868
869
Note: See TracBrowser for help on using the repository browser.