Line data Source code
1 : //# FlagAgentTimeFreqCrop.cc: This file contains the implementation of the FlagAgentTimeFreqCrop class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #include <flagging/Flagging/FlagAgentTimeFreqCrop.h>
24 :
25 : // Polynomial-fitting classes
26 : #include <casacore/scimath/Functionals/Polynomial.h>
27 : #include <casacore/scimath/Fitting.h>
28 : #include <casacore/scimath/Fitting/LinearFit.h>
29 : #include <casacore/scimath/Fitting/GenericL2Fit.h>
30 :
31 : using namespace casacore;
32 : namespace casa { //# NAMESPACE CASA - BEGIN
33 :
34 : // TODO : Get rid of this macro ?
35 : #define MIN(a,b) ((a)<=(b) ? (a) : (b))
36 :
37 86 : FlagAgentTimeFreqCrop::FlagAgentTimeFreqCrop(FlagDataHandler *dh, Record config, Bool writePrivateFlagCube, Bool flag):
38 86 : FlagAgentBase(dh,config,ANTENNA_PAIRS,writePrivateFlagCube,flag)
39 : {
40 86 : setAgentParameters(config);
41 :
42 : // Request loading polarization map to FlagDataHandler
43 86 : flagDataHandler_p->setMapPolarizations(true);
44 86 : }
45 :
46 172 : FlagAgentTimeFreqCrop::~FlagAgentTimeFreqCrop()
47 : {
48 : // Compiler automagically calls FlagAgentBase::~FlagAgentBase()
49 172 : }
50 :
51 86 : void FlagAgentTimeFreqCrop::setAgentParameters(Record config)
52 : {
53 86 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
54 : int exists;
55 :
56 86 : exists = config.fieldNumber ("timecutoff");
57 86 : if (exists >= 0)
58 : {
59 78 : if( config.type(exists) != TpDouble && config.type(exists) != TpFloat && config.type(exists) != TpInt )
60 : {
61 0 : throw( AipsError ( "Parameter 'timecutoff' must be of type 'double'" ) );
62 : }
63 :
64 78 : time_cutoff_p = config.asDouble("timecutoff");
65 : }
66 : else
67 : {
68 8 : time_cutoff_p = 4.0;
69 : }
70 :
71 86 : *logger_p << logLevel_p << " timecutoff is " << time_cutoff_p << LogIO::POST;
72 :
73 86 : exists = config.fieldNumber ("freqcutoff");
74 86 : if (exists >= 0)
75 : {
76 78 : if( config.type(exists) != TpDouble && config.type(exists) != TpFloat && config.type(exists) != TpInt )
77 : {
78 0 : throw( AipsError ( "Parameter 'freqcutoff' must be of type 'double'" ) );
79 : }
80 :
81 78 : freq_cutoff_p = config.asDouble("freqcutoff");
82 : }
83 : else
84 : {
85 8 : freq_cutoff_p = 3.0;
86 : }
87 :
88 86 : *logger_p << logLevel_p << " freqcutoff is " << freq_cutoff_p << LogIO::POST;
89 :
90 86 : exists = config.fieldNumber ("maxnpieces");
91 86 : if (exists >= 0)
92 : {
93 78 : if( config.type(exists) != TpInt )
94 : {
95 0 : throw( AipsError ( "Parameter 'maxnpieces' must be of type 'int'" ) );
96 : }
97 :
98 78 : maxNPieces_p = config.asInt("maxnpieces");
99 :
100 78 : if ((maxNPieces_p<1) or (maxNPieces_p>9))
101 : {
102 0 : throw ( AipsError (" Unsupported maxnpieces : " + String::toString(maxNPieces_p) + ". Supported values: 1-9") );
103 : }
104 : }
105 : else
106 : {
107 8 : maxNPieces_p = 7;
108 : }
109 :
110 86 : *logger_p << logLevel_p << " maxnpieces is " << maxNPieces_p << LogIO::POST;
111 :
112 86 : exists = config.fieldNumber ("timefit");
113 86 : if (exists >= 0)
114 : {
115 78 : if( config.type(exists) != TpString )
116 : {
117 0 : throw( AipsError ( "Parameter 'timefit' must be of type 'string'" ) );
118 : }
119 :
120 78 : timeFitType_p = config.asString("timefit");
121 78 : if ((timeFitType_p.compare("line") != 0) and (timeFitType_p.compare("poly") != 0))
122 : {
123 0 : throw ( AipsError ("Unsupported timefit: " + timeFitType_p +". Supported modes: line,poly" ) );
124 : }
125 :
126 : }
127 : else
128 : {
129 8 : timeFitType_p = "line";
130 : }
131 :
132 86 : *logger_p << logLevel_p << " timefit is " << timeFitType_p << LogIO::POST;
133 :
134 86 : exists = config.fieldNumber ("freqfit");
135 86 : if (exists >= 0)
136 : {
137 78 : if( config.type(exists) != TpString )
138 : {
139 0 : throw( AipsError ( "Parameter 'timefit' must be of type 'string'" ) );
140 : }
141 :
142 78 : freqFitType_p = config.asString("freqfit");
143 :
144 78 : if ((freqFitType_p.compare("line") != 0) and (freqFitType_p.compare("poly") != 0))
145 : {
146 0 : throw ( AipsError ("Unsupported freqfit: " + freqFitType_p +". Supported modes: line,poly" ) );
147 : }
148 : }
149 : else
150 : {
151 8 : freqFitType_p = "poly";
152 : }
153 :
154 86 : *logger_p << logLevel_p << " freqfit is " << freqFitType_p << LogIO::POST;
155 :
156 86 : exists = config.fieldNumber ("flagdimension");
157 86 : if (exists >= 0)
158 : {
159 78 : if( config.type(exists) != TpString )
160 : {
161 0 : throw( AipsError ( "Parameter 'flagdimension' must be of type 'string'" ) );
162 : }
163 :
164 78 : flagDimension_p = config.asString("flagdimension");
165 :
166 156 : if ((flagDimension_p.compare("time") != 0) and (flagDimension_p.compare("freq") != 0)
167 156 : and (flagDimension_p.compare("timefreq") != 0) and (flagDimension_p.compare("freqtime") != 0))
168 : {
169 0 : throw ( AipsError ("Unsupported flagdimension: " + flagDimension_p +". Supported modes: time,freq,timefreq,freqtime" ) );
170 : }
171 : }
172 : else
173 : {
174 8 : flagDimension_p = "freqtime";
175 : }
176 :
177 86 : *logger_p << logLevel_p << " flagdimension is " << flagDimension_p << LogIO::POST;
178 :
179 86 : exists = config.fieldNumber ("halfwin");
180 86 : if (exists >= 0)
181 : {
182 78 : if( config.type(exists) != TpInt )
183 : {
184 0 : throw( AipsError ( "Parameter 'halfwin' must be of type 'int'" ) );
185 : }
186 :
187 78 : halfWin_p = config.asInt("halfwin");
188 :
189 78 : if ((halfWin_p < 1) or (halfWin_p > 3))
190 : {
191 0 : throw ( AipsError ("Unsupported halfwin: " + String::toString(halfWin_p) +". Supported values : 1,2,3" ) );
192 : }
193 : }
194 : else
195 : {
196 8 : halfWin_p = 1;
197 : }
198 :
199 86 : *logger_p << logLevel_p << " halfwin is " << halfWin_p << LogIO::POST;
200 :
201 86 : exists = config.fieldNumber ("usewindowstats");
202 86 : if (exists >= 0)
203 : {
204 78 : if( config.type(exists) != TpString )
205 : {
206 0 : throw( AipsError ( "Parameter 'usewindowstats' must be of type 'string'" ) );
207 : }
208 :
209 78 : winStats_p = config.asString("usewindowstats");
210 :
211 78 : if ((winStats_p.compare("none") != 0) and (winStats_p.compare("sum") != 0)
212 78 : and (winStats_p.compare("std") != 0) and (winStats_p.compare("both") != 0))
213 : {
214 0 : throw ( AipsError ("Unsupported usewindowstats: " + winStats_p +". Supported modes : none,sum,std,both") );
215 : }
216 : }
217 : else
218 : {
219 8 : winStats_p = "none";
220 : }
221 :
222 86 : *logger_p << logLevel_p << " usewindowstats is " << winStats_p << LogIO::POST;
223 :
224 :
225 86 : return;
226 : }
227 :
228 : bool
229 8841 : FlagAgentTimeFreqCrop::computeAntennaPairFlags(const vi::VisBuffer2 & /*visBuffer*/,
230 : VisMapper &visibilities,
231 : FlagMapper &flags,
232 : Int /*antenna1*/,
233 : Int /*antenna2*/,
234 : vector<uInt> & /*rows*/)
235 : {
236 : // Call 'fltBaseAndFlag' as specified by the user.
237 8841 : if(flagDimension_p == String("time"))
238 : {
239 0 : fitBaseAndFlag(timeFitType_p,String("time"),visibilities,flags);
240 : }
241 8841 : else if( flagDimension_p == String("freq") )
242 : {
243 0 : fitBaseAndFlag(freqFitType_p,String("freq"),visibilities,flags);
244 : }
245 8841 : else if( flagDimension_p == String("timefreq") )
246 : {
247 0 : fitBaseAndFlag(timeFitType_p,String("time"),visibilities,flags);
248 0 : fitBaseAndFlag(freqFitType_p,String("freq"),visibilities,flags);
249 : }
250 : else // freqtime (default)
251 : {
252 8841 : fitBaseAndFlag(freqFitType_p,String("freq"),visibilities,flags);
253 8841 : fitBaseAndFlag(timeFitType_p,String("time"),visibilities,flags);
254 : }
255 :
256 8841 : return false;
257 : }
258 :
259 : //----------------------------------------------------------------------------------------------------------
260 :
261 : // fitBaseAndFlag
262 : // Average the data along the axis not-specified by 'direction' to generate a vector along 'direction'.
263 : // Fit a piece-wise polynomial of type 'fittype' along the axis specified by 'direction'
264 : // Divide the unaveraged data by this fit, and iteratively flag outliers - use flags.applyFlag()
265 17682 : void FlagAgentTimeFreqCrop :: fitBaseAndFlag(String fittype, String direction, VisMapper &visibilities,FlagMapper &flags)
266 : {
267 : // Get shapes
268 17682 : IPosition flagCubeShape = visibilities.shape();
269 17682 : uInt nChannels = flagCubeShape(0);
270 17682 : uInt nTimes = flagCubeShape(1);
271 :
272 : // Work variables
273 17682 : Vector<Float> avgDat,avgFit;
274 17682 : Vector<Bool> avgFlag;
275 17682 : Vector<Int> mind(2);
276 17682 : Float tol=4.0,mn=1.0,sd=0,tpsd=0.0,tpsum=0.0, mval=0.0;
277 17682 : Int mcnt=0;
278 :
279 : // Decide which direction to fit the base polynomial
280 17682 : if(direction==String("freq"))
281 : {
282 8841 : mind[0]=nChannels; mind[1]=nTimes;
283 8841 : tol=freq_cutoff_p;
284 : }
285 8841 : else if(direction==String("time") )
286 : {
287 8841 : mind[0]=nTimes; mind[1]=nChannels;
288 8841 : tol=time_cutoff_p;
289 : }
290 : else
291 : {
292 0 : throw AipsError("Internal Error. Unrecognized axis direction for tfcrop : " + direction);
293 : }
294 :
295 : // ALLOC : Resize temp arrays to either nChannel or nTime
296 17682 : avgDat.resize(mind[0]); avgDat=0.0;
297 17682 : avgFit.resize(mind[0]); avgFit=0.0;
298 17682 : avgFlag.resize(mind[0]); avgFlag=false;
299 :
300 : // A way to tell if anything is non-zero in this piece.
301 17682 : Bool allzeros=true;
302 :
303 : // STEP 1 :
304 : // For each element in the dominant direction (axis0)
305 : // calculate the mean across the other direction (axis1).
306 : // The output average data and flags will later be sent for poly fitting...
307 768919 : for(int i0=0;i0<mind[0];i0++)
308 : {
309 : // Calc the mean across axes1 (with flags)
310 751237 : mval=0.0;mcnt=0;
311 24853765 : for(uInt i1=0;i1<(uInt) mind[1];i1++)
312 : {
313 24102528 : if(mind[0]==(Int) nChannels)// if i0 is channel, and i1 is time
314 : {
315 12051276 : if( ! ( flags.getModifiedFlags(i0,i1) ) ) //C// && usePreFlags_p ) )
316 : {
317 11931035 : mval += visibilities(i0,i1);
318 11931035 : mcnt++;
319 : }
320 : }
321 : else // if i1 is channel, and i0 is time
322 : {
323 12051252 : if( ! ( flags.getModifiedFlags(i1,i0) ) ) //C// && usePreFlags_p ) )
324 : {
325 11550091 : mval += visibilities(i1,i0);
326 11550091 : mcnt++;
327 : }
328 : }
329 : }//for i1
330 :
331 : // Fill in the mean value across axis1 for i0
332 751237 : avgDat[i0] = mcnt ? mval/mcnt : 0.0;
333 : // Fill in whether all visibilities are zero or not
334 751237 : avgFlag[i0] = mcnt ? false : true;
335 751237 : if(! avgFlag[i0] ) allzeros=false;
336 :
337 : }// for i0
338 :
339 :
340 : // STEP 2 :
341 : // If there are any non-zero unflagged values in the average,
342 : // fit polynomials and flag. Otherwise, return.
343 17682 : if(allzeros==false)
344 : {
345 : // Fit a piece-wise polynomial across axis0 (or do a line)
346 : // This fit is done to the average computed in the previous step.
347 16158 : if(fittype == String("poly")) // Piecewise Poly Fit
348 : {
349 8079 : fitPiecewisePoly(avgDat,avgFlag,avgFit,maxNPieces_p,4);
350 : }
351 : else // Line Fit
352 : {
353 8079 : fitPiecewisePoly(avgDat,avgFlag,avgFit,1,1);
354 : }
355 :
356 : // STEP 3 :
357 : // Now, iterate through the data again and flag.
358 : // This time, iterate in reverse order : for each i1, get all i0
359 : // This is because 'avgFit' is of size mind[0],
360 : // and window stats are to be computed along i0
361 717831 : for(uInt i1=0;i1<(uInt) mind[1];i1++)
362 : {
363 : // STEP 3A : Divide out the clean bandpass
364 701673 : avgDat=0,avgFlag=false;
365 24608569 : for(int i0=0;i0<mind[0];i0++)
366 : {
367 23906896 : if(mind[0]==(Int) nChannels)// if i0 is channel, and i1 is time
368 : {
369 11953460 : avgFlag[i0] = flags.getModifiedFlags(i0,i1); //C// && usePreFlags_p;
370 11953460 : if(avgFlag[i0]==false) avgDat[i0] = visibilities(i0,i1)/avgFit(i0);
371 : }
372 : else // if i1 is channel, i0 is time
373 : {
374 11953436 : avgFlag[i0] = flags.getModifiedFlags(i1,i0); //C// && usePreFlags_p;
375 11953436 : if(avgFlag[i0]==false) avgDat[i0] = visibilities(i1,i0)/avgFit(i0);
376 : }
377 : }//for i0
378 :
379 : // STEP 3B
380 : // Flag outliers based on absolute deviation from the model
381 : // Do this as a robust fit
382 701673 : Float temp=0;
383 1492785 : for(Int loop=0;loop<5;loop++)
384 : {
385 : // Calculate the standard-deviation of the normalized data w.r.to the mean
386 1427325 : sd = calcStd(avgDat,avgFlag,mn);
387 :
388 : // Flag if the data differs from mn=1 by N sd
389 52452330 : for(Int i0=0;i0<mind[0];i0++)
390 : {
391 51025005 : if(avgFlag[i0]==false && fabs(avgDat[i0]-mn) > tol*sd) avgFlag[i0]=true ;
392 : }
393 :
394 : // Stop iterating if the deviation of the normalized data from the mean is less than 10%
395 1427325 : if(fabs(temp-sd) < (Double)0.1)break;
396 : // else go on for 5 iterations
397 791112 : temp=sd;
398 : }//for loop
399 :
400 : // STEP 3C :
401 : // Additional flagging based on sliding-window statistics
402 : // Note : this step looks only at vis values, not existing flags.
403 701673 : if(halfWin_p>0 && (winStats_p != "none") )
404 : {
405 : // Start and end the sliding windows halfWin from the ends.
406 0 : for(Int i0=halfWin_p;i0<mind[0]-halfWin_p;i0++)
407 : {
408 0 : tpsum=0.0;tpsd=0.0;
409 :
410 : // Flag point i0 if average of N points around i0 crosses N sd
411 0 : if(winStats_p=="sum" || winStats_p=="both")
412 : {
413 0 : for(Int i=i0-halfWin_p; i<i0+halfWin_p+1; i++) tpsum += fabs(avgDat[i]-mn);
414 0 : if(tpsum/(2*halfWin_p+1.0) > tol*sd ) avgFlag[i0]=true;
415 : }
416 :
417 : // Flag point i0 if the N point std around i0 is larger then N sd
418 0 : if(winStats_p=="std" || winStats_p=="both")
419 : {
420 0 : for(Int i=i0-halfWin_p; i<i0+halfWin_p+1; i++) tpsd += (avgDat[i]-mn) * (avgDat[i]-mn) ;
421 0 : if(sqrt( tpsd / (2*halfWin_p+1.0) ) > tol*sd) avgFlag[i0]=true ;
422 : }
423 :
424 : }//for i0
425 : }// if winStats != none
426 :
427 :
428 : // STEP 3D :
429 : // Fill the flags into the FlagMapper
430 : // Note : To minimize copies, we can call flags.applyFlag() directly from STEPS 3B and 3C,
431 : // in addition to filling in avgFlag for STEP 3B.
432 : // However, the following ensures minimal calls to flags.applyFlag().
433 24608569 : for(Int i0=0;i0<mind[0];i0++)
434 : {
435 23906896 : if(mind[0]==(Int) nChannels) // if i0 is channel, and i1 is time
436 : {
437 11953460 : if(avgFlag[i0])
438 : {
439 404918 : flags.applyFlag(i0,i1);
440 404918 : visBufferFlags_p += 1;
441 : }
442 : }
443 : else //if i1 is channel, and i0 is time
444 : {
445 11953436 : if(avgFlag[i0])
446 : {
447 411415 : flags.applyFlag(i1,i0);
448 411415 : visBufferFlags_p += 1;
449 : }
450 : }
451 : }// for i0
452 :
453 : }//for i1
454 :
455 : }// if allzeros==false
456 :
457 :
458 35364 : return;
459 :
460 17682 : }// end fitBaseAndFlag
461 :
462 :
463 :
464 : /* Calculate the MEAN of 'vect' ignoring values flagged in 'flag' */
465 89351 : Float FlagAgentTimeFreqCrop :: calcMean(Vector<Float> &vect, Vector<Bool> &flag)
466 : {
467 89351 : Float mean=0;
468 89351 : Int cnt=0;
469 3416434 : for(uInt i=0;i<vect.nelements();i++)
470 3327083 : if(flag[i]==false)
471 : {
472 3022901 : mean += vect[i];
473 3022901 : cnt++;
474 : }
475 89351 : if(cnt==0) cnt=1;
476 89351 : return mean/cnt;
477 : }
478 :
479 :
480 : /* Calculate the variance of 'vect' w.r.to a given 'fit'
481 : * ignoring values flagged in 'flag'
482 : * Use median ( data - median(data) ) as the estimator of variance
483 : */
484 0 : Float FlagAgentTimeFreqCrop :: calcVar(Vector<Float> &vect, Vector<Bool> &flag, Vector<Float> &fit)
485 : {
486 : //// Float var=0;
487 0 : uInt n=0,cnt=0;
488 0 : n = vect.nelements() < fit.nelements() ? vect.nelements() : fit.nelements();
489 0 : for(uInt i=0;i<n;i++)
490 : {
491 0 : if(flag[i]==false)cnt++;
492 : }
493 :
494 0 : Vector<Float> validvals(cnt);
495 0 : cnt=0;
496 0 : for(uInt i=0;i<n;i++)
497 : {
498 0 : if(flag[i]==false)
499 : {
500 0 : validvals[cnt] = fit[i] < (Double)1e-6 ? (Double)0.0 : fabs( (vect[i] - fit[i])/fit[i] );
501 0 : cnt++;
502 : }
503 : }
504 :
505 0 : Float med=0.0;
506 :
507 0 : if(validvals.nelements())
508 : {
509 0 : med = median(validvals,false);
510 :
511 : //cout << "validvals : " << validvals << endl;
512 : //cout << "median : " << med << endl;
513 :
514 0 : for(uInt i=0;i<validvals.nelements();i++)
515 0 : validvals[i] = fabs( validvals[i]-med );
516 :
517 0 : med = median(validvals,false);
518 : //cout << "median(data-median(data)) : " << med << endl;
519 : }
520 :
521 0 : return med;
522 0 : }
523 :
524 :
525 :
526 : /* Calculate the STANDARD DEVN. of 'vect' w.r.to a given 'fit'
527 : * ignoring values flagged in 'flag' */
528 80790 : Float FlagAgentTimeFreqCrop :: calcStd(Vector<Float> &vect, Vector<Bool> &flag, Vector<Float> &fit)
529 : {
530 80790 : Float std=0;
531 80790 : uInt n=0,cnt=0;
532 80790 : n = vect.nelements() < fit.nelements() ? vect.nelements() : fit.nelements();
533 3589155 : for(uInt i=0;i<n;i++)
534 3508365 : if(flag[i]==false)
535 : {
536 3320992 : cnt++;
537 3320992 : std += (vect[i]-fit[i])*(vect[i]-fit[i]);
538 : }
539 80790 : if(cnt==0) cnt=1;
540 80790 : return sqrt(std/cnt);
541 : }
542 :
543 :
544 : /* Calculate the STANDARD DEVN. of 'vect' w.r.to a given mean
545 : * ignoring values flagged in 'flag' */
546 1516676 : Float FlagAgentTimeFreqCrop :: calcStd(Vector<Float> &vect, Vector<Bool> &flag, Float mean)
547 : {
548 1516676 : Float std=0;
549 1516676 : uInt cnt=0;
550 55868764 : for(uInt i=0;i<vect.nelements();i++)
551 54352088 : if(flag[i]==false)
552 : {
553 52489227 : cnt++;
554 52489227 : std += (vect[i]-mean)*(vect[i]-mean);
555 : }
556 1516676 : return sqrt(std/cnt);
557 : }
558 :
559 : /* Fit Piecewise polynomials to 'data' and get the 'fit' */
560 16158 : void FlagAgentTimeFreqCrop :: fitPiecewisePoly(Vector<Float> &data, Vector<Bool> &flag, Vector<Float> &fit, uInt maxnpieces, uInt maxdeg)
561 : {
562 16158 : Int deg=0;//,start=0;
563 16158 : Int left=0,right=0;
564 16158 : Float sd,TOL=3;
565 16158 : Vector<Float> tdata;
566 :
567 : // ALLOC : Another temp array to hold modified data
568 16158 : tdata.resize(data.nelements());
569 16158 : tdata = data;
570 :
571 16158 : AlwaysAssert(data.nelements()==flag.nelements(), AipsError);
572 :
573 : /* replace empty data values by adjacent values */
574 717831 : for(uInt i=0;i<tdata.nelements();i++)
575 : {
576 701673 : if(tdata[i]==0)
577 : {
578 20800 : if(i==0)// find first non-zero value and set to that.
579 : {
580 510 : Int ind=0;
581 20620 : for(uInt j=1;j<tdata.nelements();j++)
582 20332 : if(tdata[j]!=0){ind=j;break;}
583 510 : if(ind==0) tdata[i]=0;
584 222 : else tdata[i]=tdata[ind];
585 : }
586 : else// find next non-zero value and interpolate.
587 : {
588 20290 : Int indr=0;
589 720534 : for(uInt j=i+1;j<tdata.nelements();j++)
590 703439 : if(tdata[j]!=0){indr=j;break;}
591 20290 : Int indl=-1;
592 663188 : for(int j = i ; j >= 0 ; j--)
593 646196 : if(tdata[j]!=0){indl=j;break;}
594 :
595 20290 : if(indl==-1 && indr==0) tdata[i]=0;
596 20290 : if(indl>-1 && indr==0) tdata[i]=tdata[indl];
597 20290 : if(indl==-1 && indr>0) tdata[i]=tdata[indr];
598 20290 : if(indl>-1 && indr>0) tdata[i]=(tdata[indl]+tdata[indr])/2.0;
599 : }
600 : }
601 : }
602 :
603 :
604 : /* If there still are empty points (entire spectrum is flagged) flag it. */
605 717831 : for(uInt i=0;i<tdata.nelements();i++)
606 701673 : if(tdata[i]==0)
607 : {
608 17280 : flag[i]=true;
609 : }
610 :
611 16158 : fit = tdata;
612 :
613 16158 : Int psize=1;
614 16158 : Int leftover=1,leftover_front=0,npieces=1;
615 :
616 16158 : deg=1;
617 16158 : npieces=1;
618 :
619 96948 : for(uInt j=0;j<5;j++)
620 : {
621 80790 : npieces = MIN(2*j+1, maxnpieces);
622 80790 : if(j>1) {deg=2;}
623 80790 : if(j>2) {deg=3;}
624 80790 : deg = MIN(deg,(Int) maxdeg);
625 :
626 80790 : psize = (int)(tdata.nelements()/npieces);
627 : // cout << "Iter : " << j << " with Deg : " << deg << " and Piece-size : " << psize << endl;
628 :
629 80790 : leftover = (int)(tdata.nelements() % npieces);
630 :
631 80790 : leftover_front = (int)(leftover/2.0);
632 :
633 80790 : left=0; right=tdata.nelements()-1;
634 307002 : for(Int p=0;p<npieces;p++)
635 : {
636 226212 : if(npieces>1)
637 : {
638 177738 : left = leftover_front + p*psize;
639 177738 : right = left + psize;
640 :
641 177738 : if(p==0) {left = 0; right = leftover_front + psize;}
642 177738 : if(p==npieces-1) {right = tdata.nelements()-1;}
643 : }
644 226212 : if(deg==1)
645 72711 : lineFit(tdata,flag,fit,left,right);
646 : else
647 : //lineFit(tdata,flag,fit,left,right);
648 153501 : polyFit(tdata,flag,fit,left,right,deg);
649 : }
650 :
651 : /* Now, smooth the fit - make this nicer later */
652 :
653 80790 : int winstart=0,winend=0;
654 80790 : float winsum=0.0;
655 80790 : int offset=2;
656 3274440 : for(uInt i=offset;i<tdata.nelements()-offset;i++)
657 : {
658 3193865 : winstart = i-offset;
659 3193865 : winend = i+offset;
660 3193865 : if(winstart<0)winstart=0;
661 3193865 : if(winend>=(int) tdata.nelements())winend=tdata.nelements()-1;
662 3193865 : if(winend <= winstart) break;
663 3193650 : winsum=0.0;
664 19161900 : for(uInt wi=winstart;wi<=(uInt) winend;++wi)
665 15968250 : winsum += fit[wi];
666 3193650 : fit[i] = winsum/(winend-winstart+1);
667 : }
668 :
669 :
670 : /* Calculate the STD of the fit */
671 80790 : sd = calcStd(tdata,flag,fit);
672 80790 : if(j>=2) TOL=2;
673 32316 : else TOL=3;
674 :
675 :
676 : /* Detect outliers */
677 3589155 : for(uInt i=0;i<tdata.nelements();i++)
678 : {
679 3508365 : if(tdata[i]-fit[i] > TOL*sd)
680 144278 : flag[i]=true;
681 : }
682 :
683 : } // for j
684 :
685 16158 : } // end of fitPiecewisePoly
686 :
687 :
688 :
689 : /* Fit a polynomial to 'data' from lim1 to lim2, of given degree 'deg',
690 : * taking care of flags in 'flag', and returning the fitted values in 'fit' */
691 153501 : void FlagAgentTimeFreqCrop :: polyFit(Vector<Float> &data,Vector<Bool> &flag, Vector<Float> &fit, uInt lim1, uInt lim2,uInt deg)
692 : {
693 153501 : Vector<Double> x;
694 153501 : Vector<Double> y;
695 153501 : Vector<Double> sig;
696 153501 : Vector<Double> solution;
697 :
698 153501 : uInt cnt=0;
699 1737075 : for(uInt i=lim1;i<=lim2;i++)
700 1583574 : if(flag[i]==false) cnt++;
701 :
702 153501 : if(cnt <= deg)
703 : {
704 16640 : lineFit(data,flag,fit,lim1,lim2);
705 16640 : return;
706 : }
707 :
708 :
709 136861 : LinearFit<Double> fitter;
710 136861 : Polynomial<Double> combination(deg);
711 :
712 136861 : combination.setCoefficient(0,0.0);
713 136861 : if (deg >= 1) combination.setCoefficient(1, 0.0);
714 136861 : if (deg >= 2) combination.setCoefficient(2, 0.0);
715 136861 : if (deg >= 3) combination.setCoefficient(3, 0.0);
716 136861 : if (deg >= 4) combination.setCoefficient(4, 0.0);
717 :
718 : // ALLOC : Resize to the length of each 'piece' in the piecewise fit.
719 136861 : x.resize(lim2-lim1+1);
720 136861 : y.resize(lim2-lim1+1);
721 136861 : sig.resize(lim2-lim1+1);
722 136861 : solution.resize(deg+1);
723 :
724 1660117 : for(uInt i=lim1;i<=lim2;i++)
725 : {
726 1523256 : x[i-lim1] = i+1;
727 1523256 : y[i-lim1] = data[i];
728 1523256 : sig[i-lim1] = (flag[i]==true)?0:1;
729 : }
730 :
731 136861 : fitter.asWeight(true);
732 136861 : fitter.setFunction(combination);
733 136861 : solution = fitter.fit(x,y,sig);
734 :
735 1660117 : for(uInt i=lim1;i<=lim2;i++)
736 : {
737 1523256 : fit[i]=0;
738 7114096 : for(uInt j=0;j<deg+1;j++)
739 5590840 : fit[i] += solution[j]*pow((double)(x[i-lim1]),(double)j);
740 : }
741 :
742 203421 : }
743 :
744 :
745 :
746 : /* Fit a LINE to 'data' from lim1 to lim2, taking care of flags in
747 : * 'flag', and returning the fitted values in 'fit' */
748 89351 : void FlagAgentTimeFreqCrop :: lineFit(Vector<Float> &data, Vector<Bool> &flag, Vector<Float> &fit, uInt lim1, uInt lim2)
749 : {
750 89351 : float Sx = 0, Sy = 0, Sxx = 0, Sxy = 0, S = 0, a, b, sd, mn;
751 :
752 89351 : mn = calcMean(data, flag);
753 89351 : sd = calcStd (data, flag, mn);
754 :
755 2219882 : for (uInt i = lim1; i <= lim2; i++)
756 : {
757 2130531 : if (flag[i] == false) // if unflagged
758 : {
759 2007172 : S += 1 / (sd * sd);
760 2007172 : Sx += i / (sd * sd);
761 2007172 : Sy += data[i] / (sd * sd);
762 2007172 : Sxx += (i * i) / (sd * sd);
763 2007172 : Sxy += (i * data[i]) / (sd * sd);
764 : }
765 : }
766 89351 : a = (Sxx * Sy - Sx * Sxy) / (S * Sxx - Sx * Sx);
767 89351 : b = (S * Sxy - Sx * Sy) / (S * Sxx - Sx * Sx);
768 :
769 2219882 : for (uInt i = lim1; i <= lim2; i++)
770 2130531 : fit[i] = a + b * i;
771 :
772 89351 : }
773 :
774 :
775 :
776 : } //# NAMESPACE CASA - END
777 :
778 :
|
