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

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

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