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

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

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