Changeset 2186

Timestamp:

06/07/11 23:49:53 (14 years ago)

Author:

WataruKawasaki

Message:

New Development: Yes

JIRA Issue: Yes CAS-3149

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: scantable.*sinusoid_baseline() params

Test Programs:

Put in Release Notes: Yes

Module(s):

Description: (1) Implemented an automated sinusoidal fitting functionality

(2) FFT available with scantable.fft()
(3) fixed a bug of parsing 'edge' param used by linefinder.
(4) a function to show progress status for row-based iterations.

Location:

trunk

Files:

• python/scantable.py (modified) (22 diffs)
• src/MathUtils.cpp (modified) (1 diff)
• src/MathUtils.h (modified) (4 diffs)
• src/STMath.cpp (modified) (4 diffs)
• src/STMath.h (modified) (2 diffs)
• src/Scantable.cpp (modified) (21 diffs)
• src/Scantable.h (modified) (6 diffs)
• src/ScantableWrapper.h (modified) (2 diffs)
• src/python_Scantable.cpp (modified) (1 diff)

trunk/python/scantable.py

@@ -76 +76 @@
     else:
         return False
+
+def normalise_edge_param(edge):
+    """\
+    Convert a given edge value to a one-dimensional array that can be
+    given to baseline-fitting/subtraction functions.
+    The length of the output value will be an even because values for
+    the both sides of spectra are to be contained for each IF. When
+    the length is 2, the values will be applied to all IFs. If the length
+    is larger than 2, it will be 2*ifnos(). 
+    Accepted format of edge include:
+            * an integer - will be used for both sides of spectra of all IFs.
+                  e.g. 10 is converted to [10,10]
+            * an empty list/tuple [] - converted to [0, 0] and used for all IFs. 
+            * a list/tuple containing an integer - same as the above case.
+                  e.g. [10] is converted to [10,10]
+            * a list/tuple containing two integers - will be used for all IFs. 
+                  e.g. [5,10] is output as it is. no need to convert.
+            * a list/tuple of lists/tuples containing TWO integers -
+                  each element of edge will be used for each IF.
+                  e.g. [[5,10],[15,20]] - [5,10] for IF[0] and [15,20] for IF[1].
+                  
+                  If an element contains the same integer values, the input 'edge'
+                  parameter can be given in a simpler shape in the following cases:
+                      ** when len(edge)!=2
+                          any elements containing the same values can be replaced
+                          to single integers.
+                          e.g. [[15,15]] can be simplified to [15] (or [15,15] or 15 also). 
+                          e.g. [[1,1],[2,2],[3,3]] can be simplified to [1,2,3]. 
+                      ** when len(edge)=2
+                          care is needed for this case: ONLY ONE of the
+                          elements can be a single integer,
+                          e.g. [[5,5],[10,10]] can be simplified to [5,[10,10]]
+                               or [[5,5],10], but can NOT be simplified to [5,10].
+                               when [5,10] given, it is interpreted as
+                               [[5,10],[5,10],[5,10],...] instead, as shown before. 
+    """
+    from asap import _is_sequence_or_number as _is_valid
+    if isinstance(edge, list) or isinstance(edge, tuple):
+        for edgepar in edge:
+            if not _is_valid(edgepar, int):
+                raise ValueError, "Each element of the 'edge' tuple has \
+                                   to be a pair of integers or an integer."
+            if isinstance(edgepar, list) or isinstance(edgepar, tuple):
+                if len(edgepar) != 2:
+                    raise ValueError, "Each element of the 'edge' tuple has \
+                                       to be a pair of integers or an integer."
+    else:
+        if not _is_valid(edge, int):
+            raise ValueError, "Parameter 'edge' has to be an integer or a \
+                               pair of integers specified as a tuple. \
+                               Nested tuples are allowed \
+                               to make individual selection for different IFs."
+        
+
+    if isinstance(edge, int):
+        edge = [ edge, edge ]                 # e.g. 3   => [3,3]
+    elif isinstance(edge, list) or isinstance(edge, tuple):
+        if len(edge) == 0:
+            edge = [0, 0]                     # e.g. []  => [0,0]
+        elif len(edge) == 1:
+            if isinstance(edge[0], int):
+                edge = [ edge[0], edge[0] ]   # e.g. [1] => [1,1]
+
+    commonedge = True
+    if len(edge) > 2: commonedge = False
+    else:
+        for edgepar in edge:
+            if isinstance(edgepar, list) or isinstance(edgepar, tuple):
+                commonedge = False
+                break
+
+    if commonedge:
+        if len(edge) > 1:
+            norm_edge = edge
+        else:
+            norm_edge = edge + edge
+    else:
+        norm_edge = []
+        for edgepar in edge:
+            if isinstance(edgepar, int):
+                norm_edge += [edgepar, edgepar]
+            else:
+                norm_edge += edgepar
+
+    return norm_edge
+
+def raise_fitting_failure_exception(e):
+    msg = "The fit failed, possibly because it didn't converge."
+    if rcParams["verbose"]:
+        asaplog.push(str(e))
+        asaplog.push(str(msg))
+    else:
+        raise RuntimeError(str(e)+'\n'+msg)
     
+
 class scantable(Scantable):
     """\
@@ -1120 +1214 @@
 
     @asaplog_post_dec
-    def fft(self, getrealimag=False, rowno=[]):
+    def fft(self, rowno=[], mask=[], getrealimag=False):
         """\
         Apply FFT to the spectra.
@@ -1126 +1220 @@
 
         Parameters:
+
+            rowno:          The row number(s) to be processed. int, list 
+                            and tuple are accepted. By default ([]), FFT 
+                            is applied to the whole data.
+
+            mask:           Auxiliary channel mask(s). Given as a boolean
+                            list, it is applied to all specified rows.
+                            A list of boolean lists can also be used to 
+                            apply different masks. In the latter case, the
+                            length of 'mask' must be the same as that of
+                            'rowno'. The default is []. 
         
             getrealimag:    If True, returns the real and imaginary part 
@@ -1133 +1238 @@
                             phase (argument, in unit of radian).
 
-            rowno:          The row number(s) to be processed. int, list 
-                            and tuple are accepted. By default ([]), FFT 
-                            is applied to the whole data.
-
         Returns:
 
-            A dictionary containing two keys, i.e., 'real' and 'imag' for
-            getrealimag = True, or 'ampl' and 'phase' for getrealimag = False,
+            A list of dictionaries containing the results for each spectrum.
+            Each dictionary contains two values, the real and the imaginary 
+            parts when getrealimag = True, or the amplitude(absolute value)
+            and the phase(argument) when getrealimag = False. The key for
+            these values are 'real' and 'imag', or 'ampl' and 'phase',
             respectively.
-            The value for each key is a list of lists containing the FFT
-            results from each spectrum.
-            
-        """
-        if getrealimag:
-            res = {"real":[], "imag":[]}  # return real and imaginary values
-        else:
-            res = {"ampl":[], "phase":[]} # return amplitude and phase(argument)
-        
+        """
         if isinstance(rowno, int):
             rowno = [rowno]
         elif not (isinstance(rowno, list) or isinstance(rowno, tuple)):
-            raise TypeError("The row number(s) should be int, list or tuple.")
+            raise TypeError("The row number(s) must be int, list or tuple.")
+
+        if len(rowno) == 0: rowno = [i for i in xrange(self.nrow())]
+
+        if not (isinstance(mask, list) or isinstance(mask, tuple)):
+            raise TypeError("The mask must be a boolean list or a list of boolean list.")
+        if len(mask) == 0: mask = [True for i in xrange(self.nchan())]
+        if isinstance(mask[0], bool): mask = [mask]
+        elif not (isinstance(mask[0], list) or isinstance(mask[0], tuple)):
+            raise TypeError("The mask must be a boolean list or a list of boolean list.")
+
+        usecommonmask = (len(mask) == 1)
+        if not usecommonmask:
+            if len(mask) != len(rowno):
+                raise ValueError("When specifying masks for each spectrum, the numbers of them must be identical.")
+        for amask in mask:
+            if len(amask) != self.nchan():
+                raise ValueError("The spectra and the mask have different length.")
         
-        nrow = len(rowno)
-        if nrow == 0: nrow = self.nrow()  # by default, calculate for all rows. 
-        
-        fspec = self._math._fft(self, rowno, getrealimag)
-        nspec = len(fspec)/nrow
-        
-        i = 0
-        while (i < nrow):
-            v1 = []
-            v2 = []
+        res = []
+
+        imask = 0
+        for whichrow in rowno:
+            fspec = self._fft(whichrow, mask[imask], getrealimag)
+            nspec = len(fspec)
             
-            j = 0
-            while (j < nspec):
-                k = i*nspec + j
-                v1.append(fspec[k])
-                v2.append(fspec[k+1])
-                j += 2
+            i=0
+            v1=[]
+            v2=[]
+            reselem = {"real":[],"imag":[]} if getrealimag else {"ampl":[],"phase":[]}
+            
+            while (i < nspec):
+                v1.append(fspec[i])
+                v2.append(fspec[i+1])
+                i+=2
             
             if getrealimag:
-                res["real"].append(v1)
-                res["imag"].append(v2)
+                reselem["real"]  += v1
+                reselem["imag"]  += v2
             else:
-                res["ampl"].append(v1)
-                res["phase"].append(v2)
+                reselem["ampl"]  += v1
+                reselem["phase"] += v2
             
-            i += 1
-
+            res.append(reselem)
+            
+            if not usecommonmask: imask += 1
+        
         return res
 
@@ -2131 +2245 @@
         else: return s
 
-
-    @asaplog_post_dec
-    def sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
-                          clipthresh=None, clipniter=None, plot=None, getresidual=None, outlog=None, blfile=None):
+    @asaplog_post_dec
+    def _parse_wn(self, wn):
+        if isinstance(wn, list) or isinstance(wn, tuple):
+            return wn
+        elif isinstance(wn, int):
+            return [ wn ]
+        elif isinstance(wn, str):
+            if '-' in wn:                            # case 'a-b' : return [a,a+1,...,b-1,b]
+                val = wn.split('-')
+                val = [int(val[0]), int(val[1])]
+                val.sort()
+                res = [i for i in xrange(val[0], val[1]+1)]
+            elif wn[:2] == '<=' or wn[:2] == '=<':   # cases '<=a','=<a' : return [0,1,...,a-1,a]
+                val = int(wn[2:])+1
+                res = [i for i in xrange(val)]
+            elif wn[-2:] == '>=' or wn[-2:] == '=>': # cases 'a>=','a=>' : return [0,1,...,a-1,a]
+                val = int(wn[:-2])+1
+                res = [i for i in xrange(val)]
+            elif wn[0] == '<':                       # case '<a' :         return [0,1,...,a-2,a-1]
+                val = int(wn[1:])
+                res = [i for i in xrange(val)]
+            elif wn[-1] == '>':                      # case 'a>' :         return [0,1,...,a-2,a-1]
+                val = int(wn[:-1])
+                res = [i for i in xrange(val)]
+            elif wn[:2] == '>=' or wn[:2] == '=>':   # cases '>=a','=>a' : return [a,a+1,...,a_nyq]
+                val = int(wn[2:])
+                res = [i for i in xrange(val, self.nchan()/2+1)]
+            elif wn[-2:] == '<=' or wn[-2:] == '=<': # cases 'a<=','a=<' : return [a,a+1,...,a_nyq]
+                val = int(wn[:-2])
+                res = [i for i in xrange(val, self.nchan()/2+1)]
+            elif wn[0] == '>':                       # case '>a' :         return [a+1,a+2,...,a_nyq]
+                val = int(wn[1:])+1
+                res = [i for i in xrange(val, self.nchan()/2+1)]
+            elif wn[-1] == '<':                      # case 'a<' :         return [a+1,a+2,...,a_nyq]
+                val = int(wn[:-1])+1
+                res = [i for i in xrange(val, self.nchan()/2+1)]
+
+            return res
+        else:
+            msg = 'wrong value given for addwn/rejwn'
+            raise RuntimeError(msg)
+
+
+    @asaplog_post_dec
+    def sinusoid_baseline(self, insitu=None, mask=None, applyfft=None, fftmethod=None, fftthresh=None,
+                          addwn=None, rejwn=None, clipthresh=None, clipniter=None, plot=None,
+                          getresidual=None, outlog=None, blfile=None):
         """\
         Return a scan which has been baselined (all rows) with sinusoidal functions.
         Parameters:
-            insitu:        If False a new scantable is returned.
+            insitu:        if False a new scantable is returned.
                            Otherwise, the scaling is done in-situ
                            The default is taken from .asaprc (False)
-            mask:          An optional mask
-            nwave:         the maximum wave number of sinusoids within
-                           maxwavelength * (spectral range).
-                           The default is 3 (i.e., sinusoids with wave 
-                           number of 0(=constant), 1, 2, and 3 are 
-                           used for fitting). Also it is possible to
-                           explicitly specify all the wave numbers to
-                           be used, by giving a list including them 
-                           (e.g. [0,1,2,15,16]).
-            maxwavelength: the longest sinusoidal wavelength. The
-                           default is 1.0 (unit: spectral range). 
+            mask:          an optional mask
+            applyfft:      if True use some method, such as FFT, to find
+                           strongest sinusoidal components in the wavenumber
+                           domain to be used for baseline fitting.
+                           default is True.
+            fftmethod:     method to find the strong sinusoidal components.
+                           now only 'fft' is available and it is the default.
+            fftthresh:     the threshold to select wave numbers to be used for
+                           fitting from the distribution of amplitudes in the
+                           wavenumber domain. 
+                           both float and string values accepted. 
+                           given a float value, the unit is set to sigma.
+                           for string values, allowed formats include:
+                               'xsigma' or 'x' (= x-sigma level. e.g., '3sigma'), or
+                               'topx' (= the x strongest ones, e.g. 'top5'). 
+                           default is 3.0 (unit: sigma). 
+            addwn:         the additional wave numbers to be used for fitting.
+                           list or integer value is accepted to specify every
+                           wave numbers. also string value can be used in case
+                           you need to specify wave numbers in a certain range, 
+                           e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
+                                 '<a'  (= 0,1,...,a-2,a-1),
+                                 '>=a' (= a, a+1, ... up to the maximum wave
+                                        number corresponding to the Nyquist
+                                        frequency for the case of FFT).
+                           default is []. 
+            rejwn:         the wave numbers NOT to be used for fitting. 
+                           can be set just as addwn but has higher priority: 
+                           wave numbers which are specified both in addwn
+                           and rejwn will NOT be used. default is []. 
             clipthresh:    Clipping threshold. (default is 3.0, unit: sigma)
             clipniter:     maximum number of iteration of 'clipthresh'-sigma clipping (default is 0)
@@ -2164 +2340 @@
             blfile:        Name of a text file in which the best-fit
                            parameter values to be written
-                           (default is "": no file/logger output)
+                           (default is '': no file/logger output)
 
         Example:
@@ -2170 +2346 @@
             # wave numbers in spectral window up to 10, 
             # also with 3-sigma clipping, iteration up to 4 times
-            bscan = scan.sinusoid_baseline(nwave=10,clipthresh=3.0,clipniter=4)
+            bscan = scan.sinusoid_baseline(addwn='<=10',clipthresh=3.0,clipniter=4)
 
         Note:
@@ -2177 +2353 @@
         """
         
-        varlist = vars()
+        try:
+            varlist = vars()
         
-        if insitu is None: insitu = rcParams["insitu"]
-        if insitu:
-            workscan = self
-        else:
-            workscan = self.copy()
-
-        nchan = workscan.nchan()
-        
-        if mask          is None: mask          = [True for i in xrange(nchan)]
-        if nwave         is None: nwave         = 3
-        if maxwavelength is None: maxwavelength = 1.0
-        if clipthresh    is None: clipthresh    = 3.0
-        if clipniter     is None: clipniter     = 0
-        if plot          is None: plot          = False
-        if getresidual   is None: getresidual   = True
-        if outlog        is None: outlog        = False
-        if blfile        is None: blfile        = ""
-
-        if isinstance(nwave, int):
-            in_nwave = nwave
-            nwave = []
-            for i in xrange(in_nwave+1): nwave.append(i)
-        
-        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
-        
-        try:
+            if insitu is None: insitu = rcParams['insitu']
+            if insitu:
+                workscan = self
+            else:
+                workscan = self.copy()
+            
+            if mask          is None: mask          = [True for i in xrange(workscan.nchan())]
+            if applyfft      is None: applyfft      = True
+            if fftmethod     is None: fftmethod     = 'fft'
+            if fftthresh     is None: fftthresh     = 3.0
+            if addwn         is None: addwn         = []
+            if rejwn         is None: rejwn         = []
+            if clipthresh    is None: clipthresh    = 3.0
+            if clipniter     is None: clipniter     = 0
+            if plot          is None: plot          = False
+            if getresidual   is None: getresidual   = True
+            if outlog        is None: outlog        = False
+            if blfile        is None: blfile        = ''
+
             #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method. 
-            workscan._sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile)
-            
-            workscan._add_history("sinusoid_baseline", varlist)
+            workscan._sinusoid_baseline(mask, applyfft, fftmethod.lower(), str(fftthresh).lower(), workscan._parse_wn(addwn), workscan._parse_wn(rejwn), clipthresh, clipniter, getresidual, outlog, blfile)
+            workscan._add_history('sinusoid_baseline', varlist)
             
             if insitu:
@@ -2216 +2385 @@
             
         except RuntimeError, e:
-            msg = "The fit failed, possibly because it didn't converge."
-            if rcParams["verbose"]:
-                asaplog.push(str(e))
-                asaplog.push(str(msg))
-                return
-            else:
-                raise RuntimeError(str(e)+'\n'+msg)
-
-
-    def auto_sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
-                               clipthresh=None, clipniter=None, edge=None, threshold=None,
+            raise_fitting_failure_exception(e)
+
+
+    @asaplog_post_dec
+    def auto_sinusoid_baseline(self, insitu=None, mask=None, applyfft=None, fftmethod=None, fftthresh=None,
+                               addwn=None, rejwn=None, clipthresh=None, clipniter=None, edge=None, threshold=None,
                                chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None):
         """\
@@ -2238 +2402 @@
                            The default is taken from .asaprc (False)
             mask:          an optional mask retreived from scantable
-            nwave:         the maximum wave number of sinusoids within
-                           maxwavelength * (spectral range).
-                           The default is 3 (i.e., sinusoids with wave 
-                           number of 0(=constant), 1, 2, and 3 are 
-                           used for fitting). Also it is possible to
-                           explicitly specify all the wave numbers to
-                           be used, by giving a list including them 
-                           (e.g. [0,1,2,15,16]).
-            maxwavelength: the longest sinusoidal wavelength. The
-                           default is 1.0 (unit: spectral range). 
+            applyfft:      if True use some method, such as FFT, to find
+                           strongest sinusoidal components in the wavenumber
+                           domain to be used for baseline fitting.
+                           default is True.
+            fftmethod:     method to find the strong sinusoidal components.
+                           now only 'fft' is available and it is the default.
+            fftthresh:     the threshold to select wave numbers to be used for
+                           fitting from the distribution of amplitudes in the
+                           wavenumber domain. 
+                           both float and string values accepted. 
+                           given a float value, the unit is set to sigma.
+                           for string values, allowed formats include:
+                               'xsigma' or 'x' (= x-sigma level. e.g., '3sigma'), or
+                               'topx' (= the x strongest ones, e.g. 'top5'). 
+                           default is 3.0 (unit: sigma). 
+            addwn:         the additional wave numbers to be used for fitting.
+                           list or integer value is accepted to specify every
+                           wave numbers. also string value can be used in case
+                           you need to specify wave numbers in a certain range, 
+                           e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
+                                 '<a'  (= 0,1,...,a-2,a-1),
+                                 '>=a' (= a, a+1, ... up to the maximum wave
+                                        number corresponding to the Nyquist
+                                        frequency for the case of FFT).
+                           default is []. 
+            rejwn:         the wave numbers NOT to be used for fitting. 
+                           can be set just as addwn but has higher priority: 
+                           wave numbers which are specified both in addwn
+                           and rejwn will NOT be used. default is []. 
             clipthresh:    Clipping threshold. (default is 3.0, unit: sigma)
             clipniter:     maximum number of iteration of 'clipthresh'-sigma clipping (default is 0)
@@ -2283 +2466 @@
 
         Example:
-            bscan = scan.auto_sinusoid_baseline(nwave=10, insitu=False)
+            bscan = scan.auto_sinusoid_baseline(addwn='<=10', insitu=False)
         
         Note:
@@ -2290 +2473 @@
         """
 
-        varlist = vars()
-
-        if insitu is None: insitu = rcParams['insitu']
-        if insitu:
-            workscan = self
-        else:
-            workscan = self.copy()
-
-        nchan = workscan.nchan()
-        
-        if mask           is None: mask           = [True for i in xrange(nchan)]
-        if nwave          is None: nwave          = 3
-        if maxwavelength  is None: maxwavelength  = 1.0
-        if clipthresh     is None: clipthresh     = 3.0
-        if clipniter      is None: clipniter      = 0
-        if edge           is None: edge           = (0,0)
-        if threshold      is None: threshold      = 3
-        if chan_avg_limit is None: chan_avg_limit = 1
-        if plot           is None: plot           = False
-        if getresidual    is None: getresidual    = True
-        if outlog         is None: outlog         = False
-        if blfile         is None: blfile         = ""
-
-        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
-        
-        from asap.asaplinefind import linefinder
-        from asap import _is_sequence_or_number as _is_valid
-
-        if isinstance(nwave, int):
-            in_nwave = nwave
-            nwave = []
-            for i in xrange(in_nwave+1): nwave.append(i)
-
-        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
-        individualedge = False;
-        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
-
-        if individualedge:
-            for edgepar in edge:
-                if not _is_valid(edgepar, int):
-                    raise ValueError, "Each element of the 'edge' tuple has \
-                                       to be a pair of integers or an integer."
-        else:
-            if not _is_valid(edge, int):
-                raise ValueError, "Parameter 'edge' has to be an integer or a \
-                                   pair of integers specified as a tuple. \
-                                   Nested tuples are allowed \
-                                   to make individual selection for different IFs."
-
-            if len(edge) > 1:
-                curedge = edge
+        try:
+            varlist = vars()
+
+            if insitu is None: insitu = rcParams['insitu']
+            if insitu:
+                workscan = self
             else:
-                curedge = edge + edge
-
-        try:
+                workscan = self.copy()
+            
+            if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
+            if applyfft       is None: applyfft       = True
+            if fftmethod      is None: fftmethod      = 'fft'
+            if fftthresh      is None: fftthresh      = 3.0
+            if addwn          is None: addwn          = []
+            if rejwn          is None: rejwn          = []
+            if clipthresh     is None: clipthresh     = 3.0
+            if clipniter      is None: clipniter      = 0
+            if edge           is None: edge           = (0,0)
+            if threshold      is None: threshold      = 3
+            if chan_avg_limit is None: chan_avg_limit = 1
+            if plot           is None: plot           = False
+            if getresidual    is None: getresidual    = True
+            if outlog         is None: outlog         = False
+            if blfile         is None: blfile         = ''
+
             #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method. 
-            if individualedge:
-                curedge = []
-                for i in xrange(len(edge)):
-                    curedge += edge[i]
-                
-            workscan._auto_sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
-
+            workscan._auto_sinusoid_baseline(mask, applyfft, fftmethod.lower(), str(fftthresh).lower(), workscan._parse_wn(addwn), workscan._parse_wn(rejwn), clipthresh, clipniter, normalise_edge_param(edge), threshold, chan_avg_limit, getresidual, outlog, blfile)
             workscan._add_history("auto_sinusoid_baseline", varlist)
             
@@ -2361 +2508 @@
             
         except RuntimeError, e:
-            msg = "The fit failed, possibly because it didn't converge."
-            if rcParams["verbose"]:
-                asaplog.push(str(e))
-                asaplog.push(str(msg))
-                return
-            else:
-                raise RuntimeError(str(e)+'\n'+msg)
-
+            raise_fitting_failure_exception(e)
 
     @asaplog_post_dec
@@ -2405 +2545 @@
         """
         
-        varlist = vars()
-        
-        if insitu is None: insitu = rcParams["insitu"]
-        if insitu:
-            workscan = self
-        else:
-            workscan = self.copy()
-
-        nchan = workscan.nchan()
-        
-        if mask        is None: mask        = [True for i in xrange(nchan)]
-        if npiece      is None: npiece      = 2
-        if clipthresh  is None: clipthresh  = 3.0
-        if clipniter   is None: clipniter   = 0
-        if plot        is None: plot        = False
-        if getresidual is None: getresidual = True
-        if outlog      is None: outlog      = False
-        if blfile      is None: blfile      = ""
-
-        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
-        
         try:
+            varlist = vars()
+            
+            if insitu is None: insitu = rcParams['insitu']
+            if insitu:
+                workscan = self
+            else:
+                workscan = self.copy()
+
+            if mask        is None: mask        = [True for i in xrange(workscan.nchan())]
+            if npiece      is None: npiece      = 2
+            if clipthresh  is None: clipthresh  = 3.0
+            if clipniter   is None: clipniter   = 0
+            if plot        is None: plot        = False
+            if getresidual is None: getresidual = True
+            if outlog      is None: outlog      = False
+            if blfile      is None: blfile      = ''
+
             #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 
             workscan._cspline_baseline(mask, npiece, clipthresh, clipniter, getresidual, outlog, blfile)
-            
             workscan._add_history("cspline_baseline", varlist)
             
@@ -2438 +2573 @@
             
         except RuntimeError, e:
-            msg = "The fit failed, possibly because it didn't converge."
-            if rcParams["verbose"]:
-                asaplog.push(str(e))
-                asaplog.push(str(msg))
-                return
-            else:
-                raise RuntimeError(str(e)+'\n'+msg)
-
-
+            raise_fitting_failure_exception(e)
+
+    @asaplog_post_dec
     def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None,
                               clipniter=None, edge=None, threshold=None,
@@ -2505 +2634 @@
         """
 
-        varlist = vars()
-
-        if insitu is None: insitu = rcParams['insitu']
-        if insitu:
-            workscan = self
-        else:
-            workscan = self.copy()
-
-        nchan = workscan.nchan()
-        
-        if mask           is None: mask           = [True for i in xrange(nchan)]
-        if npiece         is None: npiece         = 2
-        if clipthresh     is None: clipthresh     = 3.0
-        if clipniter      is None: clipniter      = 0
-        if edge           is None: edge           = (0, 0)
-        if threshold      is None: threshold      = 3
-        if chan_avg_limit is None: chan_avg_limit = 1
-        if plot           is None: plot           = False
-        if getresidual    is None: getresidual    = True
-        if outlog         is None: outlog         = False
-        if blfile         is None: blfile         = ""
-
-        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
-        
-        from asap.asaplinefind import linefinder
-        from asap import _is_sequence_or_number as _is_valid
-
-        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
-        individualedge = False;
-        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
-
-        if individualedge:
-            for edgepar in edge:
-                if not _is_valid(edgepar, int):
-                    raise ValueError, "Each element of the 'edge' tuple has \
-                                       to be a pair of integers or an integer."
-        else:
-            if not _is_valid(edge, int):
-                raise ValueError, "Parameter 'edge' has to be an integer or a \
-                                   pair of integers specified as a tuple. \
-                                   Nested tuples are allowed \
-                                   to make individual selection for different IFs."
-
-            if len(edge) > 1:
-                curedge = edge
+        try:
+            varlist = vars()
+
+            if insitu is None: insitu = rcParams['insitu']
+            if insitu:
+                workscan = self
             else:
-                curedge = edge + edge
-
-        try:
-            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 
-            if individualedge:
-                curedge = []
-                for i in xrange(len(edge)):
-                    curedge += edge[i]
-                
-            workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
-
+                workscan = self.copy()
+            
+            if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
+            if npiece         is None: npiece         = 2
+            if clipthresh     is None: clipthresh     = 3.0
+            if clipniter      is None: clipniter      = 0
+            if edge           is None: edge           = (0, 0)
+            if threshold      is None: threshold      = 3
+            if chan_avg_limit is None: chan_avg_limit = 1
+            if plot           is None: plot           = False
+            if getresidual    is None: getresidual    = True
+            if outlog         is None: outlog         = False
+            if blfile         is None: blfile         = ''
+
+            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
+            workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, normalise_edge_param(edge), threshold, chan_avg_limit, getresidual, outlog, blfile)
             workscan._add_history("auto_cspline_baseline", varlist)
             
@@ -2570 +2665 @@
             
         except RuntimeError, e:
-            msg = "The fit failed, possibly because it didn't converge."
-            if rcParams["verbose"]:
-                asaplog.push(str(e))
-                asaplog.push(str(msg))
-                return
-            else:
-                raise RuntimeError(str(e)+'\n'+msg)
-
+            raise_fitting_failure_exception(e)
 
     @asaplog_post_dec
@@ -2606 +2694 @@
         """
         
-        varlist = vars()
+        try:
+            varlist = vars()
         
-        if insitu is None: insitu = rcParams["insitu"]
-        if insitu:
-            workscan = self
-        else:
-            workscan = self.copy()
-
-        nchan = workscan.nchan()
-        
-        if mask        is None: mask        = [True for i in xrange(nchan)]
-        if order       is None: order       = 0
-        if plot        is None: plot        = False
-        if getresidual is None: getresidual = True
-        if outlog      is None: outlog      = False
-        if blfile      is None: blfile      = ""
-
-        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
-        
-        try:
-            rows = xrange(workscan.nrow())
-            
-            #if len(rows) > 0: workscan._init_blinfo()
+            if insitu is None: insitu = rcParams["insitu"]
+            if insitu:
+                workscan = self
+            else:
+                workscan = self.copy()
+
+            if mask        is None: mask        = [True for i in xrange(workscan.nchan())]
+            if order       is None: order       = 0
+            if plot        is None: plot        = False
+            if getresidual is None: getresidual = True
+            if outlog      is None: outlog      = False
+            if blfile      is None: blfile      = ""
 
             if plot:
+                outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
                 if outblfile: blf = open(blfile, "a")
                 
                 f = fitter()
                 f.set_function(lpoly=order)
+                
+                rows = xrange(workscan.nrow())
+                #if len(rows) > 0: workscan._init_blinfo()
+                
                 for r in rows:
                     f.x = workscan._getabcissa(r)
@@ -2673 +2758 @@
             
         except RuntimeError, e:
-            msg = "The fit failed, possibly because it didn't converge."
-            if rcParams["verbose"]:
-                asaplog.push(str(e))
-                asaplog.push(str(msg))
-                return
-            else:
-                raise RuntimeError(str(e)+'\n'+msg)
-
-
+            raise_fitting_failure_exception(e)
+
+    @asaplog_post_dec
     def auto_poly_baseline(self, insitu=None, mask=None, order=None, edge=None, threshold=None,
                            chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None):
@@ -2731 +2810 @@
         """
 
-        varlist = vars()
-
-        if insitu is None: insitu = rcParams['insitu']
-        if insitu:
-            workscan = self
-        else:
-            workscan = self.copy()
-
-        nchan = workscan.nchan()
-        
-        if mask           is None: mask           = [True for i in xrange(nchan)]
-        if order          is None: order          = 0
-        if edge           is None: edge           = (0, 0)
-        if threshold      is None: threshold      = 3
-        if chan_avg_limit is None: chan_avg_limit = 1
-        if plot           is None: plot           = False
-        if getresidual    is None: getresidual    = True
-        if outlog         is None: outlog         = False
-        if blfile         is None: blfile         = ""
-
-        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
-        
-        from asap.asaplinefind import linefinder
-        from asap import _is_sequence_or_number as _is_valid
-
-        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
-        individualedge = False;
-        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
-
-        if individualedge:
-            for edgepar in edge:
-                if not _is_valid(edgepar, int):
-                    raise ValueError, "Each element of the 'edge' tuple has \
-                                       to be a pair of integers or an integer."
-        else:
-            if not _is_valid(edge, int):
-                raise ValueError, "Parameter 'edge' has to be an integer or a \
-                                   pair of integers specified as a tuple. \
-                                   Nested tuples are allowed \
-                                   to make individual selection for different IFs."
-
-            if len(edge) > 1:
-                curedge = edge
+        try:
+            varlist = vars()
+
+            if insitu is None: insitu = rcParams['insitu']
+            if insitu:
+                workscan = self
             else:
-                curedge = edge + edge
-
-        try:
-            rows = xrange(workscan.nrow())
-            
-            #if len(rows) > 0: workscan._init_blinfo()
+                workscan = self.copy()
+
+            if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
+            if order          is None: order          = 0
+            if edge           is None: edge           = (0, 0)
+            if threshold      is None: threshold      = 3
+            if chan_avg_limit is None: chan_avg_limit = 1
+            if plot           is None: plot           = False
+            if getresidual    is None: getresidual    = True
+            if outlog         is None: outlog         = False
+            if blfile         is None: blfile         = ''
+
+            edge = normalise_edge_param(edge)
 
             if plot:
+                outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
                 if outblfile: blf = open(blfile, "a")
                 
+                from asap.asaplinefind import linefinder
                 fl = linefinder()
                 fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
                 fl.set_scan(workscan)
+                
                 f = fitter()
                 f.set_function(lpoly=order)
 
+                rows = xrange(workscan.nrow())
+                #if len(rows) > 0: workscan._init_blinfo()
+                
                 for r in rows:
-                    if individualedge:
-                        if len(edge) <= workscan.getif(r):
-                            raise RuntimeError, "Number of edge elements appear to " \
-                                  "be less than the number of IFs"
-                        else:
-                            curedge = edge[workscan.getif(r)]
-
-                    fl.find_lines(r, mask_and(mask, workscan._getmask(r)), curedge)  # (CAS-1434)
+                    idx = 2*workscan.getif(r)
+                    fl.find_lines(r, mask_and(mask, workscan._getmask(r)), edge[idx:idx+2])  # (CAS-1434)
 
                     f.x = workscan._getabcissa(r)
@@ -2827 +2876 @@
 
                 if outblfile: blf.close()
-                
             else:
-                if individualedge:
-                    curedge = []
-                    for i in xrange(len(edge)):
-                        curedge += edge[i]
-                
-                workscan._auto_poly_baseline(mask, order, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
+                workscan._auto_poly_baseline(mask, order, edge, threshold, chan_avg_limit, getresidual, outlog, blfile)
 
             workscan._add_history("auto_poly_baseline", varlist)
@@ -2844 +2887 @@
             
         except RuntimeError, e:
-            msg = "The fit failed, possibly because it didn't converge."
-            if rcParams["verbose"]:
-                asaplog.push(str(e))
-                asaplog.push(str(msg))
-                return
-            else:
-                raise RuntimeError(str(e)+'\n'+msg)
-
-
-    ### OBSOLETE ##################################################################
-    @asaplog_post_dec
-    def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
-        """
-        Return a scan which has been baselined (all rows) by a polynomial.
-        
-        Parameters:
-
-            mask:       an optional mask
-
-            order:      the order of the polynomial (default is 0)
-
-            plot:       plot the fit and the residual. In this each
-                        indivual fit has to be approved, by typing 'y'
-                        or 'n'
-
-            uselin:     use linear polynomial fit
-
-            insitu:     if False a new scantable is returned.
-                        Otherwise, the scaling is done in-situ
-                        The default is taken from .asaprc (False)
-
-            rows:       row numbers of spectra to be processed.
-                        (default is None: for all rows)
-        
-        Example:
-            # return a scan baselined by a third order polynomial,
-            # not using a mask
-            bscan = scan.poly_baseline(order=3)
-
-        """
-        if insitu is None: insitu = rcParams['insitu']
-        if not insitu:
-            workscan = self.copy()
-        else:
-            workscan = self
-        varlist = vars()
-        if mask is None:
-            mask = [True for i in xrange(self.nchan())]
-
-        try:
-            f = fitter()
-            if uselin:
-                f.set_function(lpoly=order)
-            else:
-                f.set_function(poly=order)
-
-            if rows == None:
-                rows = xrange(workscan.nrow())
-            elif isinstance(rows, int):
-                rows = [ rows ]
-            
-            if len(rows) > 0:
-                self.blpars = []
-                self.masklists = []
-                self.actualmask = []
-            
-            for r in rows:
-                f.x = workscan._getabcissa(r)
-                f.y = workscan._getspectrum(r)
-                f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
-                f.data = None
-                f.fit()
-                if plot:
-                    f.plot(residual=True)
-                    x = raw_input("Accept fit ( [y]/n ): ")
-                    if x.upper() == 'N':
-                        self.blpars.append(None)
-                        self.masklists.append(None)
-                        self.actualmask.append(None)
-                        continue
-                workscan._setspectrum(f.fitter.getresidual(), r)
-                self.blpars.append(f.get_parameters())
-                self.masklists.append(workscan.get_masklist(f.mask, row=r, silent=True))
-                self.actualmask.append(f.mask)
-
-            if plot:
-                f._p.unmap()
-                f._p = None
-            workscan._add_history("poly_baseline", varlist)
-            if insitu:
-                self._assign(workscan)
-            else:
-                return workscan
-        except RuntimeError:
-            msg = "The fit failed, possibly because it didn't converge."
-            raise RuntimeError(msg)
+            raise_fitting_failure_exception(e)
 
     def _init_blinfo(self):

trunk/src/MathUtils.cpp

@@ -248 +248 @@
   }
 }
+
+void mathutil::doZeroOrderInterpolation(casa::Vector<casa::Float>& data, 
+                                        std::vector<bool>& mask) {
+  int fstart = -1;
+  int fend = -1;
+  for (uInt i = 0; i < mask.size(); ++i) {
+    if (!mask[i]) {
+      fstart = i;
+      while (!mask[i] && i < mask.size()) {
+        fend = i;
+        i++;
+      }
+    }
+
+    // execute interpolation as the following criteria:
+    // (1) for a masked region inside the spectrum, replace the spectral 
+    //     values with the mean of those at the two channels just outside 
+    //     the both edges of the masked region. 
+    // (2) for a masked region at the spectral edge, replace the values 
+    //     with the one at the nearest non-masked channel. 
+    //     (ZOH, but bilateral)
+    Float interp = 0.0;
+    if (fstart-1 > 0) {
+      interp = data[fstart-1];
+      if (fend+1 < Int(data.nelements())) {
+        interp = (interp + data[fend+1]) / 2.0;
+      }
+    } else {
+      interp = data[fend+1];
+    }
+    if (fstart > -1 && fend > -1) {
+      for (int j = fstart; j <= fend; ++j) {
+        data[j] = interp;
+      }
+    }
+
+    fstart = -1;
+    fend = -1;
+  }
+}

trunk/src/MathUtils.h

@@ -68 +68 @@
  * @param hwidth half-width of the smoothing window
  */
- void runningMedian(casa::Vector<casa::Float>& out,
+void runningMedian(casa::Vector<casa::Float>& out,
                    casa::Vector<casa::Bool>& outflag,
                    const casa::Vector<casa::Float>& in,
@@ -74 +74 @@
                    float hwidth);
 
- void polyfit(casa::Vector<casa::Float>& out,
-                   casa::Vector<casa::Bool>& outmask,
-                   const casa::Vector<casa::Float>& in,
-                   const casa::Vector<casa::Bool>& mask,
-                   float hwidth, int order);
+void polyfit(casa::Vector<casa::Float>& out,
+             casa::Vector<casa::Bool>& outmask,
+             const casa::Vector<casa::Float>& in,
+             const casa::Vector<casa::Bool>& mask,
+             float hwidth, int order);
 
 // Generate specified statistic
@@ -85 +85 @@
 
 // Return a position of min or max value
- casa::IPosition minMaxPos(const casa::String& which,
+casa::IPosition minMaxPos(const casa::String& which,
                  const casa::MaskedArray<casa::Float>& data);
 
@@ -106 +106 @@
 casa::Vector<casa::String> toVectorString(const std::vector<std::string>& in);
 
+void doZeroOrderInterpolation(casa::Vector<casa::Float>& data, 
+                              std::vector<bool>& mask);
+
 }
 

trunk/src/STMath.cpp

@@ -49 +49 @@
 #include <atnf/PKSIO/SrcType.h>
 
-#include "MathUtils.h"
 #include "RowAccumulator.h"
 #include "STAttr.h"
@@ -2827 +2826 @@
   std::vector<float> res;
   Table tab = in->table();
-  ArrayColumn<Float> specCol(tab, "SPECTRA");
-  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
-  FFTServer<Float,Complex> ffts;
+  std::vector<bool> mask;
 
   if (whichrow.size() < 1) {          // for all rows (by default)
     int nrow = int(tab.nrow());
     for (int i = 0; i < nrow; ++i) {
-      Vector<Float> spec = specCol(i);
-      Vector<uChar> flag = flagCol(i);
-      doFFT(res, ffts, getRealImag, spec, flag);
+      res = in->execFFT(i, mask, getRealImag);
     }
   } else {                           // for specified rows
     for (uInt i = 0; i < whichrow.size(); ++i) {
-      Vector<Float> spec = specCol(whichrow[i]);
-      Vector<uChar> flag = flagCol(whichrow[i]);
-      doFFT(res, ffts, getRealImag, spec, flag);
+      res = in->execFFT(i, mask, getRealImag);
     }
   }
@@ -2849 +2842 @@
 }
 
-void asap::STMath::doFFT( std::vector<float>& out, 
-                          FFTServer<Float,Complex>& ffts, 
-                          bool getRealImag,
-                          Vector<Float>& spec, 
-                          Vector<uChar>& flag )
-{
-  doZeroOrderInterpolation(spec, flag);
-
-  Vector<Complex> fftout;
-  ffts.fft0(fftout, spec);
-
-  float norm = float(2.0/double(spec.size()));
-  if (getRealImag) {
-    for (uInt j = 0; j < fftout.size(); ++j) {
-      out.push_back(real(fftout[j])*norm);
-      out.push_back(imag(fftout[j])*norm);
-    }
-  } else {
-    for (uInt j = 0; j < fftout.size(); ++j) {
-      out.push_back(abs(fftout[j])*norm);
-      out.push_back(arg(fftout[j]));
-    }
-  }
-
-}
-
-void asap::STMath::doZeroOrderInterpolation( Vector<Float>& spec, 
-                                             Vector<uChar>& flag )
-{
-  int fstart = -1;
-  int fend = -1;
-  for (uInt i = 0; i < flag.nelements(); ++i) {
-    if (flag[i] > 0) {
-      fstart = i;
-      while (flag[i] > 0 && i < flag.nelements()) {
-        fend = i;
-        i++;
-      }
-    }
-
-    // execute interpolation as the following criteria:
-    // (1) for a masked region inside the spectrum, replace the spectral 
-    //     values with the mean of those at the two channels just outside 
-    //     the both edges of the masked region. 
-    // (2) for a masked region at the spectral edge, replace the values 
-    //     with the one at the nearest non-masked channel. 
-    //     (ZOH, but bilateral)
-    Float interp = 0.0;
-    if (fstart-1 > 0) {
-      interp = spec[fstart-1];
-      if (fend+1 < Int(spec.nelements())) {
-        interp = (interp + spec[fend+1]) / 2.0;
-      }
-    } else {
-      interp = spec[fend+1];
-    }
-    if (fstart > -1 && fend > -1) {
-      for (int j = fstart; j <= fend; ++j) {
-        spec[j] = interp;
-      }
-    }
-
-    fstart = -1;
-    fend = -1;
-  }
-}
 
 CountedPtr<Scantable>
@@ -2939 +2866 @@
       Vector<Float> spec = specCol(i);
       Vector<uChar> flag = flagCol(i);
-      doZeroOrderInterpolation(spec, flag);
+      std::vector<bool> mask;
+      for (uInt j = 0; j < flag.nelements(); ++j) {
+        mask.push_back(!(flag[j]>0));
+      }
+      mathutil::doZeroOrderInterpolation(spec, mask);
 
       Vector<Complex> lags;

trunk/src/STMath.h

@@ -21 +21 @@
 #include <casa/Utilities/CountedPtr.h>
 
-#include <scimath/Mathematics/FFTServer.h>
 #include <scimath/Mathematics/InterpolateArray1D.h>
 
@@ -433 +432 @@
                                         int index ) ;
   double getMJD( string strtime ) ;
-  void doFFT( std::vector<float>& out, 
-              casa::FFTServer< casa::Float, casa::Complex >& ffts, 
-              bool getRealImag, 
-              casa::Vector<casa::Float>& spec, 
-              casa::Vector<casa::uChar>& flag ) ;
-  void doZeroOrderInterpolation( casa::Vector<casa::Float>& spec, 
-                                 casa::Vector<casa::uChar>& flag) ;
 
   bool insitu_;

trunk/src/Scantable.cpp

@@ -11 +11 @@
 //
 #include <map>
-#include <fstream>
+
+#include <atnf/PKSIO/SrcType.h>
 
 #include <casa/aips.h>
+#include <casa/iomanip.h>
 #include <casa/iostream.h>
-#include <casa/iomanip.h>
+#include <casa/OS/File.h>
 #include <casa/OS/Path.h>
-#include <casa/OS/File.h>
 #include <casa/Arrays/Array.h>
+#include <casa/Arrays/ArrayAccessor.h>
+#include <casa/Arrays/ArrayLogical.h>
 #include <casa/Arrays/ArrayMath.h>
 #include <casa/Arrays/MaskArrMath.h>
-#include <casa/Arrays/ArrayLogical.h>
-#include <casa/Arrays/ArrayAccessor.h>
+#include <casa/Arrays/Slice.h>
 #include <casa/Arrays/Vector.h>
 #include <casa/Arrays/VectorSTLIterator.h>
-#include <casa/Arrays/Slice.h>
 #include <casa/BasicMath/Math.h>
 #include <casa/BasicSL/Constants.h>
+#include <casa/Containers/RecordField.h>
+#include <casa/Logging/LogIO.h>
 #include <casa/Quanta/MVAngle.h>
-#include <casa/Containers/RecordField.h>
+#include <casa/Quanta/MVTime.h>
 #include <casa/Utilities/GenSort.h>
-#include <casa/Logging/LogIO.h>
-
-#include <tables/Tables/TableParse.h>
-#include <tables/Tables/TableDesc.h>
-#include <tables/Tables/TableCopy.h>
-#include <tables/Tables/SetupNewTab.h>
-#include <tables/Tables/ScaColDesc.h>
-#include <tables/Tables/ArrColDesc.h>
-#include <tables/Tables/TableRow.h>
-#include <tables/Tables/TableVector.h>
-#include <tables/Tables/TableIter.h>
-
-#include <tables/Tables/ExprNode.h>
-#include <tables/Tables/TableRecord.h>
-#include <casa/Quanta/MVTime.h>
-#include <casa/Quanta/MVAngle.h>
-#include <measures/Measures/MeasRef.h>
-#include <measures/Measures/MeasTable.h>
+
+#include <coordinates/Coordinates/CoordinateUtil.h>
+
 // needed to avoid error in .tcc
 #include <measures/Measures/MCDirection.h>
 //
 #include <measures/Measures/MDirection.h>
+#include <measures/Measures/MEpoch.h>
 #include <measures/Measures/MFrequency.h>
-#include <measures/Measures/MEpoch.h>
+#include <measures/Measures/MeasRef.h>
+#include <measures/Measures/MeasTable.h>
+#include <measures/TableMeasures/ScalarMeasColumn.h>
+#include <measures/TableMeasures/TableMeasDesc.h>
 #include <measures/TableMeasures/TableMeasRefDesc.h>
 #include <measures/TableMeasures/TableMeasValueDesc.h>
-#include <measures/TableMeasures/TableMeasDesc.h>
-#include <measures/TableMeasures/ScalarMeasColumn.h>
-#include <coordinates/Coordinates/CoordinateUtil.h>
-
-#include <atnf/PKSIO/SrcType.h>
-#include "Scantable.h"
-#include "STPolLinear.h"
-#include "STPolCircular.h"
-#include "STPolStokes.h"
+
+#include <tables/Tables/ArrColDesc.h>
+#include <tables/Tables/ExprNode.h>
+#include <tables/Tables/ScaColDesc.h>
+#include <tables/Tables/SetupNewTab.h>
+#include <tables/Tables/TableCopy.h>
+#include <tables/Tables/TableDesc.h>
+#include <tables/Tables/TableIter.h>
+#include <tables/Tables/TableParse.h>
+#include <tables/Tables/TableRecord.h>
+#include <tables/Tables/TableRow.h>
+#include <tables/Tables/TableVector.h>
+
+#include "MathUtils.h"
 #include "STAttr.h"
 #include "STLineFinder.h"
-#include "MathUtils.h"
+#include "STPolCircular.h"
+#include "STPolLinear.h"
+#include "STPolStokes.h"
+#include "Scantable.h"
 
 using namespace casa;
@@ -1847 +1847 @@
     setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
+    showProgressOnTerminal(whichrow, nRow);
   }
 
@@ -1909 +1910 @@
 
     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
+    showProgressOnTerminal(whichrow, nRow);
   }
 
@@ -1950 +1952 @@
 
     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
+    showProgressOnTerminal(whichrow, nRow);
   }
 
@@ -2018 +2021 @@
 
     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
+    showProgressOnTerminal(whichrow, nRow);
   }
 
@@ -2243 +2247 @@
 }
 
-void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, bool getResidual, bool outLogger, const std::string& blfile)
+  void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
+{
+  nWaves.clear();
+
+  if (applyFFT) {
+    string fftThAttr;
+    float fftThSigma;
+    int fftThTop;
+    parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
+    doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
+  }
+
+  addAuxWaveNumbers(addNWaves, rejectNWaves, nWaves);
+}
+
+void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
+{
+  uInt idxSigma = fftThresh.find("sigma");
+  uInt idxTop   = fftThresh.find("top");
+
+  if (idxSigma == fftThresh.size() - 5) {
+    std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
+    is >> fftThSigma;
+    fftThAttr = "sigma";
+  } else if (idxTop == 0) {
+    std::istringstream is(fftThresh.substr(3));
+    is >> fftThTop;
+    fftThAttr = "top";
+  } else {
+    bool isNumber = true;
+    for (uInt i = 0; i < fftThresh.size()-1; ++i) {
+      char ch = (fftThresh.substr(i, 1).c_str())[0];
+      if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
+        isNumber = false;
+        break;
+      }
+    }
+    if (isNumber) {
+      std::istringstream is(fftThresh);
+      is >> fftThSigma;
+      fftThAttr = "sigma";
+    } else {
+      throw(AipsError("fftthresh has a wrong value"));
+    }
+  }
+}
+
+void Scantable::doSelectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftMethod, const float fftThSigma, const int fftThTop, const std::string& fftThAttr, std::vector<int>& nWaves)
+{
+  std::vector<float> fspec;
+  if (fftMethod == "fft") {
+    fspec = execFFT(whichrow, chanMask, false, true);
+  //} else if (fftMethod == "lsp") {
+  //  fspec = lombScarglePeriodogram(whichrow);
+  }
+
+  if (fftThAttr == "sigma") {
+    float mean  = 0.0;
+    float mean2 = 0.0;
+    for (uInt i = 0; i < fspec.size(); ++i) {
+      mean  += fspec[i];
+      mean2 += fspec[i]*fspec[i];
+    }
+    mean  /= float(fspec.size());
+    mean2 /= float(fspec.size());
+    float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
+
+    for (uInt i = 0; i < fspec.size(); ++i) {
+      if (fspec[i] >= thres) {
+        nWaves.push_back(i);
+      }
+    }
+
+  } else if (fftThAttr == "top") {
+    for (int i = 0; i < fftThTop; ++i) {
+      float max = 0.0;
+      int maxIdx = 0;
+      for (uInt j = 0; j < fspec.size(); ++j) {
+        if (fspec[j] > max) {
+          max = fspec[j];
+          maxIdx = j;
+        }
+      }
+      nWaves.push_back(maxIdx);
+      fspec[maxIdx] = 0.0;
+    }
+
+  }
+
+  if (nWaves.size() > 1) {
+    sort(nWaves.begin(), nWaves.end());
+  }
+}
+
+void Scantable::addAuxWaveNumbers(const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
+{
+  for (uInt i = 0; i < addNWaves.size(); ++i) {
+    bool found = false;
+    for (uInt j = 0; j < nWaves.size(); ++j) {
+      if (nWaves[j] == addNWaves[i]) {
+        found = true;
+        break;
+      }
+    }
+    if (!found) nWaves.push_back(addNWaves[i]);
+  }
+
+  for (uInt i = 0; i < rejectNWaves.size(); ++i) {
+    for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
+      if (*j == rejectNWaves[i]) {
+        j = nWaves.erase(j);
+      } else {
+        ++j;
+      }
+    }
+  }
+
+  if (nWaves.size() > 1) {
+    sort(nWaves.begin(), nWaves.end());
+    unique(nWaves.begin(), nWaves.end());
+  }
+}
+
+void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, bool outLogger, const std::string& blfile)
 {
   ofstream ofs;
@@ -2267 +2394 @@
   int nRow = nrow();
   std::vector<bool> chanMask;
+  std::vector<int> nWaves;
 
   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     chanMask = getCompositeChanMask(whichrow, mask);
+    selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
+
+    //FOR DEBUGGING------------
+    if (whichrow < 0) {// == nRow -1) {
+      cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
+      if (applyFFT) {
+          cout << "[ ";
+          for (uInt j = 0; j < nWaves.size(); ++j) {
+            cout << nWaves[j] << ", ";
+          }
+          cout << " ]    " << endl;
+      }
+      cout << flush;
+    }
+    //-------------------------
+
     //fitBaseline(chanMask, whichrow, fitter);
     //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     std::vector<float> params;
-    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
     setSpectrum(res, whichrow);
     //
 
     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
+    showProgressOnTerminal(whichrow, nRow);
   }
 
@@ -2283 +2428 @@
 }
 
-void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, bool outLogger, const std::string& blfile)
+void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, bool outLogger, const std::string& blfile)
 {
   ofstream ofs;
@@ -2307 +2452 @@
   int nRow = nrow();
   std::vector<bool> chanMask;
+  std::vector<int> nWaves;
+
   int minEdgeSize = getIFNos().size()*2;
   STLineFinder lineFinder = STLineFinder();
@@ -2336 +2483 @@
     //-------------------------------------------------------
 
+    selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
 
     //fitBaseline(chanMask, whichrow, fitter);
     //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     std::vector<float> params;
-    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
     setSpectrum(res, whichrow);
     //
 
     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
+    showProgressOnTerminal(whichrow, nRow);
   }
 
@@ -2350 +2499 @@
 }
 
-std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, float maxWaveLength, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
+std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
 {
   if (data.size() != mask.size()) {
@@ -2359 +2508 @@
   }
   if (waveNumbers.size() == 0) {
-    throw(AipsError("missing wave number info"));
+    throw(AipsError("no wave numbers given"));
   }
   std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
@@ -2390 +2539 @@
 
   const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
-  double baseXFactor = 2.0*PI/(double)maxWaveLength/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
+  double baseXFactor = 2.0*PI/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
 
   // xArray : contains elemental values for computing the least-square matrix. 
@@ -2571 +2720 @@
 std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
 {
-  std::vector<bool> chanMask = getMask(whichrow);
-  uInt chanMaskSize = chanMask.size();
-  if (chanMaskSize != inMask.size()) {
-    throw(AipsError("different mask sizes"));
-  }
-  for (uInt i = 0; i < chanMaskSize; ++i) {
-    chanMask[i] = chanMask[i] && inMask[i];
-  }
-
-  return chanMask;
+  std::vector<bool> mask = getMask(whichrow);
+  uInt maskSize = mask.size();
+  if (maskSize != inMask.size()) {
+    throw(AipsError("mask sizes are not the same."));
+  }
+  for (uInt i = 0; i < maskSize; ++i) {
+    mask[i] = mask[i] && inMask[i];
+  }
+
+  return mask;
 }
 
@@ -2610 +2759 @@
 
 /* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
-void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter) {
+void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter)
+{
   if (outLogger || outTextFile) {
     std::vector<float> params = fitter.getParameters();
@@ -2628 +2778 @@
 
 /* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
-void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params) {
+void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params)
+{
   if (outLogger || outTextFile) {
     float rms = getRms(chanMask, whichrow);
@@ -2646 +2797 @@
 
 /* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
-void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params) {
+void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params)
+{
   if (outLogger || outTextFile) {
     float rms = getRms(chanMask, whichrow);
@@ -2663 +2815 @@
 }
 
-float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
+void Scantable::showProgressOnTerminal(const int nProcessed, const int nTotal, const int nTotalThreshold)
+{
+  if (nTotal >= nTotalThreshold) {
+    int nInterval = int(floor(double(nTotal)/100.0));
+    if (nInterval == 0) nInterval++;
+
+    if (nProcessed == 0) {
+      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
+      printf("[  0%%]");
+      printf("\x1b[39m\x1b[0m");             //default attributes
+      fflush(NULL);
+    } else if (nProcessed % nInterval == 0) {
+      printf("\r\x1b[1C");                   //go to the 2nd column
+      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
+      printf("%3d", (int)(100.0*(double(nProcessed+1))/(double(nTotal))) );
+      printf("\x1b[39m\x1b[0m");             //default attributes
+      printf("\x1b[2C");                     //go to the end of line
+      fflush(NULL);
+    }
+    if (nProcessed == nTotal - 1) {
+      printf("\r\x1b[K");                    //clear
+      fflush(NULL);
+    }
+  }
+}
+
+std::vector<float> Scantable::execFFT(const int whichrow, const std::vector<bool>& inMask, bool getRealImag, bool getAmplitudeOnly)
+{
+  std::vector<bool>  mask = getMask(whichrow);
+
+  if (inMask.size() > 0) {
+    uInt maskSize = mask.size();
+    if (maskSize != inMask.size()) {
+      throw(AipsError("mask sizes are not the same."));
+    }
+    for (uInt i = 0; i < maskSize; ++i) {
+      mask[i] = mask[i] && inMask[i];
+    }
+  }
+
+  Vector<Float> spec = getSpectrum(whichrow);
+  mathutil::doZeroOrderInterpolation(spec, mask);
+
+  FFTServer<Float,Complex> ffts;
+  Vector<Complex> fftres;
+  ffts.fft0(fftres, spec);
+
+  std::vector<float> res;
+  float norm = float(2.0/double(spec.size()));
+
+  if (getRealImag) {
+    for (uInt i = 0; i < fftres.size(); ++i) {
+      res.push_back(real(fftres[i])*norm);
+      res.push_back(imag(fftres[i])*norm);
+    }
+  } else {
+    for (uInt i = 0; i < fftres.size(); ++i) {
+      res.push_back(abs(fftres[i])*norm);
+      if (!getAmplitudeOnly) res.push_back(arg(fftres[i]));
+    }
+  }
+
+  return res;
+}
+
+
+float Scantable::getRms(const std::vector<bool>& mask, int whichrow)
+{
   Vector<Float> spec;
   specCol_.get(whichrow, spec);
@@ -2685 +2904 @@
 
 
-std::string Scantable::formatBaselineParamsHeader(int whichrow, 
-                                                  const std::string& masklist, 
-                                                  bool verbose) const
+std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
 {
   ostringstream oss;
@@ -2722 +2939 @@
 }
 
-  std::string Scantable::formatBaselineParams(const std::vector<float>& params, 
-                                              const std::vector<bool>& fixed, 
-                                              float rms, 
-                                              const std::string& masklist, 
-                                              int whichrow, 
-                                              bool verbose, 
-                                              int start, int count,
-                                              bool resetparamid) const
+std::string Scantable::formatBaselineParams(const std::vector<float>& params, 
+                                            const std::vector<bool>& fixed, 
+                                            float rms, 
+                                            const std::string& masklist, 
+                                            int whichrow, 
+                                            bool verbose, 
+                                            int start, int count,
+                                            bool resetparamid) const
 {
   int nParam = (int)(params.size());
@@ -2904 +3121 @@
 }
 
-/*
-STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
-{
-  Fitter fitter = Fitter();
-  fitter.setExpression("poly", order);
-
-  std::vector<bool> fmask = getMask(rowno);
-  if (fmask.size() != mask.size()) {
-    throw(AipsError("different mask sizes"));
-  }
-  for (int i = 0; i < fmask.size(); ++i) {
-    fmask[i] = fmask[i] && mask[i];
-  }
-
-  fitBaseline(fmask, rowno, fitter);
-  setSpectrum(fitter.getResidual(), rowno);
-  return fitter.getFitEntry();
-}
-*/
 
 }

trunk/src/Scantable.h

@@ -14 +14 @@
 
 // STL
+#include <fstream>
+#include <iostream>
+#include <sstream>
 #include <string>
 #include <vector>
-#include <iostream>
-#include <fstream>
 // AIPS++
 #include <casa/aips.h>
+#include <casa/Arrays/MaskedArray.h>
+#include <casa/Arrays/Vector.h>
+#include <casa/BasicSL/String.h>
 #include <casa/Containers/Record.h>
-#include <casa/Arrays/MaskedArray.h>
-#include <casa/BasicSL/String.h>
+#include <casa/Exceptions/Error.h>
+#include <casa/Quanta/Quantum.h>
 #include <casa/Utilities/CountedPtr.h>
 
-#include <tables/Tables/Table.h>
+#include <coordinates/Coordinates/SpectralCoordinate.h>
+
+#include <measures/TableMeasures/ScalarMeasColumn.h>
+
+#include <scimath/Mathematics/FFTServer.h>
+
 #include <tables/Tables/ArrayColumn.h>
 #include <tables/Tables/ScalarColumn.h>
-
-#include <measures/TableMeasures/ScalarMeasColumn.h>
-
-#include <coordinates/Coordinates/SpectralCoordinate.h>
-
-#include <casa/Arrays/Vector.h>
-#include <casa/Quanta/Quantum.h>
-
-#include <casa/Exceptions/Error.h>
+#include <tables/Tables/Table.h>
 
 #include "Logger.h"
-#include "STHeader.h"
-#include "STFrequencies.h"
-#include "STWeather.h"
-#include "STFocus.h"
-#include "STTcal.h"
-#include "STMolecules.h"
-#include "STSelector.h"
-#include "STHistory.h"
-#include "STPol.h"
+#include "MathUtils.h"
 #include "STFit.h"
 #include "STFitEntry.h"
 #include "STFitter.h"
+#include "STFocus.h"
+#include "STFrequencies.h"
+#include "STHeader.h"
+#include "STHistory.h"
+#include "STMolecules.h"
+#include "STPol.h"
+#include "STSelector.h"
+#include "STTcal.h"
+#include "STWeather.h"
 
 namespace asap {
@@ -528 +530 @@
                                const std::string& blfile="");
   void sinusoidBaseline(const std::vector<bool>& mask,
-                        const std::vector<int>& nWaves,
-                        float maxWaveLength,
+                        const bool applyFFT, 
+                        const std::string& fftMethod, 
+                        const std::string& fftThresh, 
+                        const std::vector<int>& addNWaves, 
+                        const std::vector<int>& rejectNWaves, 
                         float thresClip, 
                         int nIterClip,
@@ -536 +541 @@
                         const std::string& blfile="");
   void autoSinusoidBaseline(const std::vector<bool>& mask,
-                            const std::vector<int>& nWaves,
-                            float maxWaveLength,
+                            const bool applyFFT, 
+                            const std::string& fftMethod, 
+                            const std::string& fftThresh, 
+                            const std::vector<int>& addNWaves, 
+                            const std::vector<int>& rejectNWaves, 
                             float thresClip, 
                             int nIterClip,
@@ -546 +554 @@
                             bool outLogger=false, 
                             const std::string& blfile="");
+  std::vector<float> execFFT(const int whichrow, 
+                             const std::vector<bool>& inMask,
+                             bool getRealImag=false,
+                             bool getAmplitudeOnly=false);
   float getRms(const std::vector<bool>& mask, int whichrow);
   std::string formatBaselineParams(const std::vector<float>& params, 
@@ -689 +701 @@
                                        const std::vector<bool>& mask,
                                        const std::vector<int>& waveNumbers,
-                                       float maxWaveLength, 
                                        std::vector<float>& params,
                                        float thresClip=3.0, 
                                        int nIterClip=1,
                                        bool getResidual=true);
+  void selectWaveNumbers(const int whichrow, 
+                         const std::vector<bool>& chanMask, 
+                         const bool applyFFT, 
+                         const std::string& fftMethod, 
+                         const std::string& fftThresh, 
+                         const std::vector<int>& addNWaves, 
+                         const std::vector<int>& rejectNWaves, 
+                         std::vector<int>& nWaves);
+  void parseThresholdExpression(const std::string& fftThresh,
+                                std::string& fftThAttr,
+                                float& fftThSigma,
+                                int& fftThTop);
+  void doSelectWaveNumbers(const int whichrow,
+                           const std::vector<bool>& chanMask, 
+                           const std::string& fftMethod, 
+                           const float fftThSigma, 
+                           const int fftThTop, 
+                           const std::string& fftThAttr, 
+                           std::vector<int>& nWaves);
+  void addAuxWaveNumbers(const std::vector<int>& addNWaves, 
+                         const std::vector<int>& rejectNWaves, 
+                         std::vector<int>& nWaves);
   bool hasSameNchanOverIFs();
   std::string getMaskRangeList(const std::vector<bool>& mask, 
@@ -708 +741 @@
   void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params);
   void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<float>& params);
+  void showProgressOnTerminal(const int nProcessed, const int nTotal, const int nTotalThreshold=1000);
 
   void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

trunk/src/ScantableWrapper.h

@@ -276 +276 @@
   { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
 
-  void sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, bool getresidual=true, bool outlog=false, const std::string& blfile="")
-  { table_->sinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile); }
-
-  void autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, bool outlog=false, const std::string& blfile="")
-  { table_->autoSinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
+  void sinusoidBaseline(const std::vector<bool>& mask, const bool applyfft, const std::string& fftmethod, const std::string& fftthresh, const std::vector<int>& addwn, const std::vector<int>& rejwn, float clipthresh, int clipniter, bool getresidual=true, bool outlog=false, const std::string& blfile="")
+  { table_->sinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, getresidual, outlog, blfile); }
+
+  void autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyfft, const std::string& fftmethod, const std::string& fftthresh, const std::vector<int>& addwn, const std::vector<int>& rejwn, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, bool outlog=false, const std::string& blfile="")
+  { table_->autoSinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
 
   float getRms(const std::vector<bool>& mask, int whichrow)
@@ -293 +293 @@
   bool getFlagtraFast(int whichrow=0) const
   { return table_->getFlagtraFast(casa::uInt(whichrow)); }
+
+  std::vector<float> execFFT(int whichrow, const std::vector<bool>& mask, bool getRealImag=false, bool getAmplitudeOnly=false)
+  { return table_->execFFT(whichrow, mask, getRealImag, getAmplitudeOnly); }
 
 

trunk/src/python_Scantable.cpp

@@ -154 +154 @@
     .def("_getflagtrafast", &ScantableWrapper::getFlagtraFast,
          (boost::python::arg("whichrow")=0) )
+    .def("_fft", &ScantableWrapper::execFFT)
     //.def("_sspline_baseline", &ScantableWrapper::smoothingSplineBaseline)
     //.def("_test_cin", &ScantableWrapper::testCin)