Changeset 2794 for trunk

Timestamp:

03/15/13 18:23:04 (12 years ago)

Author:

Kana Sugimoto

Message:

New Development: Yes

JIRA Issue: Yes (CAS-4141/TRAC-290)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sbseparator module

Description: Pushed python codes down to cpp.

Location:

trunk

Files:

• python/sbseparator.py (modified) (2 diffs)
• src/STSideBandSep.cpp (modified) (14 diffs)
• src/STSideBandSep.h (modified) (8 diffs)
• src/python_STSideBandSep.cpp (modified) (1 diff)

trunk/python/sbseparator.py

@@ -38 +38 @@
     """
     def __init__(self, infiles):
-        self.intables = None
-        self.signalShift = []
-        self.imageShift = []
-        self.dsbmode = True
-        self.getboth = False
-        self.rejlimit = 0.2
-        self.baseif = -1
-        self.freqtol = 10.
-        self.freqframe = ""
-        self.solveother = False
-        self.dirtol = [1.e-5, 1.e-5] # direction tolerance in rad (2 arcsec)
-        #self.lo1 = 0.
-
-        self.tables = []
-        self.nshift = -1
-        self.nchan = -1
-
-        self.set_data(infiles)
-        
         self._separator = SBSeparator(infiles)
 
-    @asaplog_post_dec
-    def set_data(self, infiles):
-        """
-        Set data to be processed.
-
-        infiles  : a list of filenames or scantables
-        """
-        if not (type(infiles) in (list, tuple, numpy.ndarray)):
-            infiles = [infiles]
-        if isinstance(infiles[0], scantable):
-            # a list of scantable
-            for stab in infiles:
-                if not isinstance(stab, scantable):
-                    asaplog.post()
-                    raise TypeError, "Input data is not a list of scantables."
-
-            #self._separator._setdata(infiles)
-            self._reset_data()
-            self.intables = infiles
-        else:
-            # a list of filenames
-            for name in infiles:
-                if not os.path.exists(name):
-                    asaplog.post()
-                    raise ValueError, "Could not find input file '%s'" % name
-            
-            #self._separator._setdataname(infiles)
-            self._reset_data()
-            self.intables = infiles
-
-        asaplog.push("%d files are set to process" % len(self.intables))
-
-
-    def _reset_data(self):
-        del self.intables
-        self.intables = None
-        self.signalShift = []
-        #self.imageShift = []
-        self.tables = []
-        self.nshift = -1
-        self.nchan = -1
-
-    @asaplog_post_dec
+################## temp functions #################
+#     def rfft(self, invec):
+#         print "Python method cppRfft"
+#         self._separator._cpprfft(invec)
+#         return 
+################## temp functions #################
+
     def set_frequency(self, baseif, freqtol, frame=""):
         """
@@ -108 +53 @@
         Parameters:
           - reference IFNO to process in the first table in the list
-          - frequency tolerance from reference IF to select data
+          - frequency tolerance from reference IF to select data (string)
           frame  : frequency frame to select IF
         """
-        self._reset_if()
-        self.baseif = baseif
-        if isinstance(freqtol,dict) and freqtol["unit"] == "Hz":
-            if freqtol['value'] > 0.:
-                self.freqtol = freqtol
+        if type(freqtol) in (float, int):
+            freqtol = str(freqtol)
+        elif isinstance(freqtol, dict):
+            try:
+                freqtol = str(freqtol['value']) + freqtol['unit']
+            except:
+                raise ValueError, "Invalid frequency tolerance."
+        self._separator.set_freq(baseif, freqtol, frame)
+
+
+    def set_dirtol(self, dirtol=["2arcsec", "2arcsec"]):
+        """
+        Set tolerance of direction to select data
+        """
+        if isinstance(dirtol, str):
+            dirtol = [dirtol]
+
+        self._separator.set_dirtol(dirtol)
+    
+            
+    def set_shift(self, imageshift=[]):
+        """
+        Set shift mode and channel shift of image band.
+
+        imageshift : a list of number of channels shifted in image
+                     side band of each scantable.
+                     If the shifts are not set, they are assumed to be
+                     equal to those of signal side band, but in opposite
+                     direction as usual by LO1 offsetting of DSB receivers.
+        """
+        if not imageshift:
+            imageshift = []
+        self._separator.set_shift(imageshift)
+
+
+    @asaplog_post_dec
+    def set_both(self, flag=False):
+        """
+        Resolve both image and signal sideband when True.
+        """
+        self._separator.solve_both(flag)
+        if flag:
+            asaplog.push("Both signal and image sidebands will be solved and stored in separate tables.")
+        else:
+            asaplog.push("Only signal sideband will be solved and stored in an table.")
+
+    @asaplog_post_dec
+    def set_limit(self, threshold=0.2):
+        """
+        Set rejection limit of solution.
+        """
+        self._separator.set_limit(threshold)
+
+
+    @asaplog_post_dec
+    def set_solve_other(self, flag=False):
+        """
+        Calculate spectra by subtracting the solution of the other sideband
+        when True.
+        """
+        self._separator.subtract_other(flag)
+        if flag:
+            asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.")
+
+
+    def set_lo1(self,lo1, frame="TOPO", reftime=-1, refdir=""):
+        """
+        Set LO1 frequency to calculate frequency of image sideband.
+
+        lo1     : LO1 frequency
+        frame   : the frequency frame of LO1
+        reftime : the reference time used in frequency frame conversion.
+        refdir  : the reference direction used in frequency frame conversion.
+        """
+        self._separator.set_lo1(lo1val, frame, reftime, refdir)
+
+
+    def set_lo1root(self, name):
+        """
+        Set MS name which stores LO1 frequency of signal side band.
+        It is used to calculate frequency of image sideband.
+
+        name : MS name which contains 'ASDM_SPECTRALWINDOW' and
+               'ASDM_RECEIVER' tables.
+        """
+        self._separator.set_lo1root(name)
+
+    @asaplog_post_dec
+    def separate(self, outname="", overwrite=False):
+        """
+        Invoke sideband separation.
+
+        outname   : a name of output scantable
+        overwrite : overwrite existing table
+        """
+        out_default = "sbseparated.asap"
+        if len(outname) == 0:
+            outname = out_default
+            asaplog.post()
+            asaplog.push("The output file name is not specified.")
+            asaplog.push("Using default name '%s'" % outname)
+            asaplog.post("WARN")
+
+        if os.path.exists(outname):
+            if overwrite:
+                asaplog.push("removing the old file '%s'" % outname)
+                shutil.rmtree(outname)
             else:
                 asaplog.post()
-                asaplog.push("Frequency tolerance should be positive value.")
+                asaplog.push("Output file '%s' already exists." % outname)
                 asaplog.post("ERROR")
-                return
-        else:
-            # torelance in channel unit
-            if freqtol > 0:
-                self.freqtol = float(freqtol)
-            else:
-                asaplog.post()
-                asaplog.push("Frequency tolerance should be positive value.")
-                asaplog.post("ERROR")
-                return
-        self.freqframe = frame
-
-    def _reset_if(self):
-        self.baseif = -1
-        self.freqtol = 0
-        self.freqframe = ""
-        self.signalShift = []
-        #self.imageShift = []
-        self.tables = []
-        self.nshift = 0
-        self.nchan = -1
-
-    @asaplog_post_dec
-    def set_dirtol(self, dirtol=[1.e-5,1.e-5]):
-        """
-        Set tolerance of direction to select data
-        """
-        # direction tolerance in rad
-        if not (type(dirtol) in [list, tuple, numpy.ndarray]):
-            dirtol = [dirtol, dirtol]
-        if len(dirtol) == 1:
-            dirtol = [dirtol[0], dirtol[0]]
-        if len(dirtol) > 1:
-            self.dirtol = dirtol[0:2]
-        else:
-            asaplog.post()
-            asaplog.push("Invalid direction tolerance. Should be a list of float in unit radian")
-            asaplog.post("ERROR")
-            return
-        asaplog.post("Set direction tolerance [%f, %f] (rad)" % \
-                     (self.dirtol[0], self.dirtol[1]))
-
-    @asaplog_post_dec
-    def set_shift(self, mode="DSB", imageshift=None):
-        """
-        Set shift mode and channel shift of image band.
-
-        mode       : shift mode ['DSB'|'SSB'(='2SB')]
-                     When mode='DSB', imageshift is assumed to be equal
-                     to the shift of signal sideband but in opposite direction.
-        imageshift : a list of number of channel shift in image band of
-                     each scantable. valid only mode='SSB'
-        """
-        if mode.upper().startswith("D"):
-            # DSB mode
-            self.dsbmode = True
-            self.imageShift = []
-        else:
-            if not imageshift:
-                raise ValueError, "Need to set shift value of image sideband"
-            self.dsbmode = False
-            self.imageShift = imageshift
-            asaplog.push("Image sideband shift is set manually: %s" % str(self.imageShift))
-
-    @asaplog_post_dec
-    def set_both(self, flag=False):
-        """
-        Resolve both image and signal sideband when True.
-        """
-        self.getboth = flag
-        if self.getboth:
-            asaplog.push("Both signal and image sidebands are solved and output as separate tables.")
-        else:
-            asaplog.push("Only signal sideband is solved and output as an table.")
-
-    @asaplog_post_dec
-    def set_limit(self, threshold=0.2):
-        """
-        Set rejection limit of solution.
-        """
-        #self._separator._setlimit(abs(threshold))
-        self.rejlimit = threshold
-        asaplog.push("The threshold of rejection is set to %f" % self.rejlimit)
-
-
-    @asaplog_post_dec
-    def set_solve_other(self, flag=False):
-        """
-        Calculate spectra by subtracting the solution of the other sideband
-        when True.
-        """
-        self.solveother = flag
-        if flag:
-            asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.")
-
-    def set_lo1(self,lo1):
-        """
-        Set LO1 frequency to calculate frequency of image sideband.
-
-        lo1 : LO1 frequency in float
-        """
-        lo1val = -1.
-        if isinstance(lo1, dict) and lo1["unit"] == "Hz":
-            lo1val = lo1["value"]
-        else:
-            lo1val = float(lo1)
-        if lo1val <= 0.:
-            asaplog.push("Got negative LO1 frequency. It will be ignored.")
-            asaplog.post("WARN")
-        else:
-            self._separator.set_lo1(lo1val)
-
-
-    def set_lo1root(self, name):
-        """
-        Set MS name which stores LO1 frequency of signal side band.
-        It is used to calculate frequency of image sideband.
-
-        name : MS name which contains 'ASDM_SPECTRALWINDOW' and
-               'ASDM_RECEIVER' tables.
-        """
-        self._separator.set_lo1root(name)
-
-
-    @asaplog_post_dec
-    def separate(self, outname="", overwrite=False):
-        """
-        Invoke sideband separation.
-
-        outname   : a name of output scantable
-        overwrite : overwrite existing table
-        """
-        # List up valid scantables and IFNOs to convolve.
-        #self._separator._separate()
-        self._setup_shift()
-        #self._preprocess_tables()
-
-        ### TEMPORAL ###
-        self._separator._get_asistb_from_scantb(self.tables[0])
-        ################
-
-        nshift = len(self.tables)
-        signaltab = self._grid_outtable(self.tables[0])
-        if self.getboth:
-            imagetab = signaltab.copy()
-
-        rejrow = []
-        for irow in xrange(signaltab.nrow()):
-            currpol = signaltab.getpol(irow)
-            currbeam = signaltab.getbeam(irow)
-            currdir = signaltab.get_directionval(irow)
-            spec_array, tabidx = self._get_specarray(polid=currpol,\
-                                                     beamid=currbeam,\
-                                                     dir=currdir)
-            #if not spec_array:
-            if len(tabidx) == 0:
-                asaplog.post()
-                asaplog.push("skipping row %d" % irow)
-                rejrow.append(irow)
-                continue
-            signal = self._solve_signal(spec_array, tabidx)
-            signaltab.set_spectrum(signal, irow)
-
-            # Solve image side side band
-            if self.getboth:
-                image = self._solve_image(spec_array, tabidx)
-                imagetab.set_spectrum(image, irow)
-
-        # TODO: Need to remove rejrow form scantables here
-        signaltab.flag_row(rejrow)
-        if self.getboth:
-            imagetab.flag_row(rejrow)
+                return False
+
+        self._separator.separate(outname)
+
+######################################################################
+#     @asaplog_post_dec
+#     def separate0(self, outname="", overwrite=False):
+#         """
+#         Invoke sideband separation.
+
+#         outname   : a name of output scantable
+#         overwrite : overwrite existing table
+#         """
+#         # List up valid scantables and IFNOs to convolve.
+#         #self._separator.separate()
+#         self._setup_shift()
+#         #self._preprocess_tables()
+
+#         ### TEMPORAL ###
+#         self._separator._get_asistb_from_scantb(self.tables[0])
+#         ################
+
+#         nshift = len(self.tables)
+#         signaltab = self._grid_outtable(self.tables[0])
+#         if self.getboth:
+#             imagetab = signaltab.copy()
+
+#         rejrow = []
+#         for irow in xrange(signaltab.nrow()):
+#             currpol = signaltab.getpol(irow)
+#             currbeam = signaltab.getbeam(irow)
+#             currdir = signaltab.get_directionval(irow)
+#             spec_array, tabidx = self._get_specarray(polid=currpol,\
+#                                                      beamid=currbeam,\
+#                                                      dir=currdir)
+#             #if not spec_array:
+#             if len(tabidx) == 0:
+#                 asaplog.post()
+#                 asaplog.push("skipping row %d" % irow)
+#                 rejrow.append(irow)
+#                 continue
+#             signal = self._solve_signal(spec_array, tabidx)
+#             signaltab.set_spectrum(signal, irow)
+
+#             # Solve image side side band
+#             if self.getboth:
+#                 image = self._solve_image(spec_array, tabidx)
+#                 imagetab.set_spectrum(image, irow)
+
+#         # TODO: Need to remove rejrow form scantables here
+#         signaltab.flag_row(rejrow)
+#         if self.getboth:
+#             imagetab.flag_row(rejrow)
         
-        if outname == "":
-            outname = "sbsepareted.asap"
-        signalname = outname + ".signalband"
-        if os.path.exists(signalname):
-            if not overwrite:
-                raise Exception, "Output file '%s' exists." % signalname
-            else:
-                shutil.rmtree(signalname)
-        signaltab.save(signalname)
-
-        if self.getboth:
-            # Warnings
-            asaplog.post()
-            asaplog.push("Saving IMAGE sideband.")
-            #asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")
-            #asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")
-            #asaplog.post("WARN")
-
-            imagename = outname + ".imageband"
-            if os.path.exists(imagename):
-                if not overwrite:
-                    raise Exception, "Output file '%s' exists." % imagename
-                else:
-                    shutil.rmtree(imagename)
-            # Update frequency information
-            self._separator.set_imgtable(imagetab)
-            self._separator.solve_imgfreq()
-            imagetab.save(imagename)
-
-
-    def _solve_signal(self, data, tabidx=None):
-        if not tabidx:
-            tabidx = range(len(data))
-
-        tempshift = []
-        dshift = []
-        if self.solveother:
-            for idx in tabidx:
-                tempshift.append(-self.imageShift[idx])
-                dshift.append(self.signalShift[idx] - self.imageShift[idx])
-        else:
-            for idx in tabidx:
-                tempshift.append(-self.signalShift[idx])
-                dshift.append(self.imageShift[idx] - self.signalShift[idx])
-
-        shiftdata = numpy.zeros(data.shape, numpy.float)
-        for i in range(len(data)):
-            shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
-        ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
-        result_image = self._combineResult(ifftdata)
-        if not self.solveother:
-            return result_image
-        result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
-        return result_signal
-
-
-    def _solve_image(self, data, tabidx=None):
-        if not tabidx:
-            tabidx = range(len(data))
-
-        tempshift = []
-        dshift = []
-        if self.solveother:
-            for idx in tabidx:
-                tempshift.append(-self.signalShift[idx])
-                dshift.append(self.imageShift[idx] - self.signalShift[idx])
-        else:
-            for idx in tabidx:
-                tempshift.append(-self.imageShift[idx])
-                dshift.append(self.signalShift[idx] - self.imageShift[idx])
-
-        shiftdata = numpy.zeros(data.shape, numpy.float)
-        for i in range(len(data)):
-            shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
-        ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
-        result_image = self._combineResult(ifftdata)
-        if not self.solveother:
-            return result_image
-        result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
-        return result_signal
-
-    @asaplog_post_dec
-    def _grid_outtable(self, table):
-        # Generate gridded table for output (Just to get rows)
-        gridder = asapgrid2(table)
-        gridder.setIF(self.baseif)
-
-        cellx = str(self.dirtol[0])+"rad"
-        celly = str(self.dirtol[1])+"rad"
-        dirarr = numpy.array(table.get_directionval()).transpose()
-        mapx = dirarr[0].max() - dirarr[0].min()
-        mapy = dirarr[1].max() - dirarr[1].min()
-        centy = 0.5 * (dirarr[1].max() + dirarr[1].min())
-        nx = max(1, numpy.ceil(mapx*numpy.cos(centy)/self.dirtol[0]))
-        ny = max(1, numpy.ceil(mapy/self.dirtol[0]))
-
-        asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
-                     (cellx, celly))
-        gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
-        gridder.setFunc(func='box', convsupport=1)
-        gridder.setWeight(weightType='uniform')
-        gridder.grid()
-        return gridder.getResult()
-
-    @asaplog_post_dec
-    def _get_specarray(self, polid=None, beamid=None, dir=None):
-        ntable = len(self.tables)
-        spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
-        nspec = 0
-        asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
-        tabidx = []
-        for itab in range(ntable):
-            tab = self.tables[itab]
-            # Select rows by POLNO and BEAMNO
-            try:
-                tab.set_selection(pols=[polid], beams=[beamid])
-                if tab.nrow() > 0: tabidx.append(itab)
-            except: # no selection
-                asaplog.post()
-                asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
-                asaplog.post("WARN")
-                continue
-
-            # Select rows by direction
-            spec = numpy.zeros(self.nchan, numpy.float)
-            selrow = []
-            for irow in range(tab.nrow()):
-                currdir = tab.get_directionval(irow)
-                if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
-                   (abs(currdir[1] - dir[1]) > self.dirtol[1]):
-                    continue
-                selrow.append(irow)
-            if len(selrow) == 0:
-                asaplog.post()
-                asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
-                asaplog.post("WARN")
-                continue
-            else:
-                seltab = tab.copy()
-                seltab.set_selection(selector(rows=selrow))
-
-            if tab.nrow() > 1:
-                asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
-                tab = seltab.average_time(scanav=False, weight="tintsys")
-            else:
-                tab = seltab
-
-            spec_array[nspec] = tab._getspectrum()
-            nspec += 1
-
-        if nspec != ntable:
-            asaplog.post()
-            #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
-            asaplog.push("Could not find corresponding rows in some tables.")
-            asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
-            if nspec < 2:
-                asaplog.push("At least 2 spectra are necessary for convolution")
-                asaplog.post("ERROR")
-                return False, tabidx
-
-        return spec_array[0:nspec], tabidx
-
-
-    @asaplog_post_dec
-    def _setup_shift(self):
-        ### define self.tables, self.signalShift, and self.imageShift
-        if not self.intables:
-            asaplog.post()
-            raise RuntimeError, "Input data is not defined."
-        #if self.baseif < 0:
-        #    asaplog.post()
-        #    raise RuntimeError, "Reference IFNO is not defined."
+#         if outname == "":
+#             outname = "sbsepareted.asap"
+#         signalname = outname + ".signalband"
+#         if os.path.exists(signalname):
+#             if not overwrite:
+#                 raise Exception, "Output file '%s' exists." % signalname
+#             else:
+#                 shutil.rmtree(signalname)
+#         signaltab.save(signalname)
+
+#         if self.getboth:
+#             # Warnings
+#             asaplog.post()
+#             asaplog.push("Saving IMAGE sideband.")
+#             #asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")
+#             #asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")
+#             #asaplog.post("WARN")
+
+#             imagename = outname + ".imageband"
+#             if os.path.exists(imagename):
+#                 if not overwrite:
+#                     raise Exception, "Output file '%s' exists." % imagename
+#                 else:
+#                     shutil.rmtree(imagename)
+#             # Update frequency information
+#             self._separator.set_imgtable(imagetab)
+#             self._separator.solve_imgfreq()
+#             imagetab.save(imagename)
+
+
+#     def _solve_signal(self, data, tabidx=None):
+#         if not tabidx:
+#             tabidx = range(len(data))
+
+#         tempshift = []
+#         dshift = []
+#         if self.solveother:
+#             for idx in tabidx:
+#                 tempshift.append(-self.imageShift[idx])
+#                 dshift.append(self.signalShift[idx] - self.imageShift[idx])
+#         else:
+#             for idx in tabidx:
+#                 tempshift.append(-self.signalShift[idx])
+#                 dshift.append(self.imageShift[idx] - self.signalShift[idx])
+
+#         shiftdata = numpy.zeros(data.shape, numpy.float)
+#         for i in range(len(data)):
+#             shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
+#         ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
+#         result_image = self._combineResult(ifftdata)
+#         if not self.solveother:
+#             return result_image
+#         result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
+#         return result_signal
+
+
+#     def _solve_image(self, data, tabidx=None):
+#         if not tabidx:
+#             tabidx = range(len(data))
+
+#         tempshift = []
+#         dshift = []
+#         if self.solveother:
+#             for idx in tabidx:
+#                 tempshift.append(-self.signalShift[idx])
+#                 dshift.append(self.imageShift[idx] - self.signalShift[idx])
+#         else:
+#             for idx in tabidx:
+#                 tempshift.append(-self.imageShift[idx])
+#                 dshift.append(self.signalShift[idx] - self.imageShift[idx])
+
+#         shiftdata = numpy.zeros(data.shape, numpy.float)
+#         for i in range(len(data)):
+#             shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
+#         ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
+#         result_image = self._combineResult(ifftdata)
+#         if not self.solveother:
+#             return result_image
+#         result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
+#         return result_signal
+
+#     @asaplog_post_dec
+#     def _grid_outtable(self, table):
+#         # Generate gridded table for output (Just to get rows)
+#         gridder = asapgrid2(table)
+#         gridder.setIF(self.baseif)
+
+#         cellx = str(self.dirtol[0])+"rad"
+#         celly = str(self.dirtol[1])+"rad"
+#         dirarr = numpy.array(table.get_directionval()).transpose()
+#         mapx = dirarr[0].max() - dirarr[0].min()
+#         mapy = dirarr[1].max() - dirarr[1].min()
+#         centy = 0.5 * (dirarr[1].max() + dirarr[1].min())
+#         nx = max(1, numpy.ceil(mapx*numpy.cos(centy)/self.dirtol[0]))
+#         ny = max(1, numpy.ceil(mapy/self.dirtol[0]))
+
+#         asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
+#                      (cellx, celly))
+#         gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
+#         gridder.setFunc(func='box', convsupport=1)
+#         gridder.setWeight(weightType='uniform')
+#         gridder.grid()
+#         return gridder.getResult()
+
+#     @asaplog_post_dec
+#     def _get_specarray(self, polid=None, beamid=None, dir=None):
+#         ntable = len(self.tables)
+#         spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
+#         nspec = 0
+#         asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
+#         tabidx = []
+#         for itab in range(ntable):
+#             tab = self.tables[itab]
+#             # Select rows by POLNO and BEAMNO
+#             try:
+#                 tab.set_selection(pols=[polid], beams=[beamid])
+#                 if tab.nrow() > 0: tabidx.append(itab)
+#             except: # no selection
+#                 asaplog.post()
+#                 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
+#                 asaplog.post("WARN")
+#                 continue
+
+#             # Select rows by direction
+#             spec = numpy.zeros(self.nchan, numpy.float)
+#             selrow = []
+#             for irow in range(tab.nrow()):
+#                 currdir = tab.get_directionval(irow)
+#                 if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
+#                    (abs(currdir[1] - dir[1]) > self.dirtol[1]):
+#                     continue
+#                 selrow.append(irow)
+#             if len(selrow) == 0:
+#                 asaplog.post()
+#                 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
+#                 asaplog.post("WARN")
+#                 continue
+#             else:
+#                 seltab = tab.copy()
+#                 seltab.set_selection(selector(rows=selrow))
+
+#             if seltab.nrow() > 1:
+#                 asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
+#                 tab = seltab.average_time(scanav=False, weight="tintsys")
+#             else:
+#                 tab = seltab
+
+#             spec_array[nspec] = tab._getspectrum()
+#             nspec += 1
+
+#         if nspec != ntable:
+#             asaplog.post()
+#             #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
+#             asaplog.push("Could not find corresponding rows in some tables.")
+#             asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
+#             if nspec < 2:
+#                 asaplog.push("At least 2 spectra are necessary for convolution")
+#                 asaplog.post("ERROR")
+#                 return False, tabidx
+
+#         return spec_array[0:nspec], tabidx
+
+
+#     @asaplog_post_dec
+#     def _setup_shift(self):
+#         ### define self.tables, self.signalShift, and self.imageShift
+#         if not self.intables:
+#             asaplog.post()
+#             raise RuntimeError, "Input data is not defined."
+#         #if self.baseif < 0:
+#         #    asaplog.post()
+#         #    raise RuntimeError, "Reference IFNO is not defined."
         
-        byname = False
-        #if not self.intables:
-        if isinstance(self.intables[0], str):
-            # A list of file name is given
-            if not os.path.exists(self.intables[0]):
-                asaplog.post()
-                raise RuntimeError, "Could not find '%s'" % self.intables[0]
+#         byname = False
+#         #if not self.intables:
+#         if isinstance(self.intables[0], str):
+#             # A list of file name is given
+#             if not os.path.exists(self.intables[0]):
+#                 asaplog.post()
+#                 raise RuntimeError, "Could not find '%s'" % self.intables[0]
             
-            stab = scantable(self.intables[0],average=False)
-            ntab = len(self.intables)
-            byname = True
-        else:
-            stab = self.intables[0]
-            ntab = len(self.intables)
-
-        if len(stab.getbeamnos()) > 1:
-            asaplog.post()
-            asaplog.push("Mult-beam data is not supported by this module.")
-            asaplog.post("ERROR")
-            return
-
-        valid_ifs = stab.getifnos()
-        if self.baseif < 0:
-            self.baseif = valid_ifs[0]
-            asaplog.post()
-            asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
-
-        if not (self.baseif in valid_ifs):
-            asaplog.post()
-            errmsg = "IF%d does not exist in the first scantable" %  \
-                     self.baseif
-            raise RuntimeError, errmsg
-
-        asaplog.push("Start selecting tables and IFNOs to solve.")
-        asaplog.push("Checking frequency of the reference IF")
-        unit_org = stab.get_unit()
-        coord = stab._getcoordinfo()
-        frame_org = coord[1]
-        stab.set_unit("Hz")
-        if len(self.freqframe) > 0:
-            stab.set_freqframe(self.freqframe)
-        stab.set_selection(ifs=[self.baseif])
-        spx = stab._getabcissa()
-        stab.set_selection()
-        basech0 = spx[0]
-        baseinc = spx[1]-spx[0]
-        self.nchan = len(spx)
-        # frequency tolerance
-        if isinstance(self.freqtol, dict) and self.freqtol['unit'] == "Hz":
-            vftol = abs(self.freqtol['value'])
-        else:
-            vftol = abs(baseinc * float(self.freqtol))
-            self.freqtol = dict(value=vftol, unit="Hz")
-        # tolerance of frequency increment
-        inctol = abs(baseinc/float(self.nchan))
-        asaplog.push("Reference frequency setup (Table = 0, IFNO = %d):  nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
-        asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
-        poltype0 = stab.poltype()
-
-        self.tables = []
-        self.signalShift = []
-        if self.dsbmode:
-            self.imageShift = []
-
-        for itab in range(ntab):
-            asaplog.push("Table %d:" % itab)
-            tab_selected = False
-            if itab > 0:
-                if byname:
-                    stab = scantable(self.intables[itab],average=False)
-                else:
-                    stab = self.intables[itab]
-                unit_org = stab.get_unit()
-                coord = stab._getcoordinfo()
-                frame_org = coord[1]
-                stab.set_unit("Hz")
-                if len(self.freqframe) > 0:
-                    stab.set_freqframe(self.freqframe)
-
-            # Check POLTYPE should be equal to the first table.
-            if stab.poltype() != poltype0:
-                asaplog.post()
-                raise Exception, "POLTYPE should be equal to the first table."
-            # Multiple beam data may not handled properly
-            if len(stab.getbeamnos()) > 1:
-                asaplog.post()
-                asaplog.push("table contains multiple beams. It may not be handled properly.")
-                asaplog.push("WARN")
-
-            for ifno in stab.getifnos():
-                stab.set_selection(ifs=[ifno])
-                spx = stab._getabcissa()
-                if (len(spx) != self.nchan) or \
-                   (abs(spx[0]-basech0) > vftol) or \
-                   (abs(spx[1]-spx[0]-baseinc) > inctol):
-                    continue
-                tab_selected = True
-                seltab = stab.copy()
-                seltab.set_unit(unit_org)
-                seltab.set_freqframe(frame_org)
-                self.tables.append(seltab)
-                self.signalShift.append((spx[0]-basech0)/baseinc)
-                if self.dsbmode:
-                    self.imageShift.append(-self.signalShift[-1])
-                asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
-            stab.set_selection()
-            stab.set_unit(unit_org)
-            stab.set_freqframe(frame_org)
-            if not tab_selected:
-                asaplog.post()
-                asaplog.push("No data selected in Table %d" % itab)
-                asaplog.post("WARN")
-
-        asaplog.push("Total number of IFs selected = %d" % len(self.tables))
-        if len(self.tables) < 2:
-            asaplog.post()
-            raise RuntimeError, "At least 2 IFs are necessary for convolution!"
-
-        if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
-            asaplog.post()
-            errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
-            errmsg += "\nThe frequency tolerance (freqtol) you set may be too small."
-            raise RuntimeError, errmsg
-
-        self.signalShift = numpy.array(self.signalShift)
-        self.imageShift = numpy.array(self.imageShift)
-        self.nshift = len(self.tables)
-
-    @asaplog_post_dec
-    def _preprocess_tables(self):
-        ### temporary method to preprocess data
-        ### Do time averaging for now.
-        for itab in range(len(self.tables)):
-            self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
+#             stab = scantable(self.intables[0],average=False)
+#             ntab = len(self.intables)
+#             byname = True
+#         else:
+#             stab = self.intables[0]
+#             ntab = len(self.intables)
+
+#         if len(stab.getbeamnos()) > 1:
+#             asaplog.post()
+#             asaplog.push("Mult-beam data is not supported by this module.")
+#             asaplog.post("ERROR")
+#             return
+
+#         valid_ifs = stab.getifnos()
+#         if self.baseif < 0:
+#             self.baseif = valid_ifs[0]
+#             asaplog.post()
+#             asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
+
+#         if not (self.baseif in valid_ifs):
+#             asaplog.post()
+#             errmsg = "IF%d does not exist in the first scantable" %  \
+#                      self.baseif
+#             raise RuntimeError, errmsg
+
+#         asaplog.push("Start selecting tables and IFNOs to solve.")
+#         asaplog.push("Checking frequency of the reference IF")
+#         unit_org = stab.get_unit()
+#         coord = stab._getcoordinfo()
+#         frame_org = coord[1]
+#         stab.set_unit("Hz")
+#         if len(self.freqframe) > 0:
+#             stab.set_freqframe(self.freqframe)
+#         stab.set_selection(ifs=[self.baseif])
+#         spx = stab._getabcissa()
+#         stab.set_selection()
+#         basech0 = spx[0]
+#         baseinc = spx[1]-spx[0]
+#         self.nchan = len(spx)
+#         # frequency tolerance
+#         if isinstance(self.freqtol, dict) and self.freqtol['unit'] == "Hz":
+#             vftol = abs(self.freqtol['value'])
+#         else:
+#             vftol = abs(baseinc * float(self.freqtol))
+#             self.freqtol = dict(value=vftol, unit="Hz")
+#         # tolerance of frequency increment
+#         inctol = abs(baseinc/float(self.nchan))
+#         asaplog.push("Reference frequency setup (Table = 0, IFNO = %d):  nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
+#         asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
+#         poltype0 = stab.poltype()
+
+#         self.tables = []
+#         self.signalShift = []
+#         if self.dsbmode:
+#             self.imageShift = []
+
+#         for itab in range(ntab):
+#             asaplog.push("Table %d:" % itab)
+#             tab_selected = False
+#             if itab > 0:
+#                 if byname:
+#                     stab = scantable(self.intables[itab],average=False)
+#                 else:
+#                     stab = self.intables[itab]
+#                 unit_org = stab.get_unit()
+#                 coord = stab._getcoordinfo()
+#                 frame_org = coord[1]
+#                 stab.set_unit("Hz")
+#                 if len(self.freqframe) > 0:
+#                     stab.set_freqframe(self.freqframe)
+
+#             # Check POLTYPE should be equal to the first table.
+#             if stab.poltype() != poltype0:
+#                 asaplog.post()
+#                 raise Exception, "POLTYPE should be equal to the first table."
+#             # Multiple beam data may not handled properly
+#             if len(stab.getbeamnos()) > 1:
+#                 asaplog.post()
+#                 asaplog.push("table contains multiple beams. It may not be handled properly.")
+#                 asaplog.push("WARN")
+
+#             for ifno in stab.getifnos():
+#                 stab.set_selection(ifs=[ifno])
+#                 spx = stab._getabcissa()
+#                 if (len(spx) != self.nchan) or \
+#                    (abs(spx[0]-basech0) > vftol) or \
+#                    (abs(spx[1]-spx[0]-baseinc) > inctol):
+#                     continue
+#                 tab_selected = True
+#                 seltab = stab.copy()
+#                 seltab.set_unit(unit_org)
+#                 seltab.set_freqframe(frame_org)
+#                 self.tables.append(seltab)
+#                 self.signalShift.append((spx[0]-basech0)/baseinc)
+#                 if self.dsbmode:
+#                     self.imageShift.append(-self.signalShift[-1])
+#                 asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
+#             stab.set_selection()
+#             stab.set_unit(unit_org)
+#             stab.set_freqframe(frame_org)
+#             if not tab_selected:
+#                 asaplog.post()
+#                 asaplog.push("No data selected in Table %d" % itab)
+#                 asaplog.post("WARN")
+
+#         asaplog.push("Total number of IFs selected = %d" % len(self.tables))
+#         if len(self.tables) < 2:
+#             asaplog.post()
+#             raise RuntimeError, "At least 2 IFs are necessary for convolution!"
+
+#         if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
+#             asaplog.post()
+#             errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
+#             errmsg += "\nThe frequency tolerance (freqtol) you set may be too small."
+#             raise RuntimeError, errmsg
+
+#         self.signalShift = numpy.array(self.signalShift)
+#         self.imageShift = numpy.array(self.imageShift)
+#         self.nshift = len(self.tables)
 
 #     @asaplog_post_dec
-#     def _setup_image_freq(self, table):
-#         # User defined coordinate
-#         # Get from associated MS
-#         # Get from ASDM
-#         lo1 = -1.
-#         if self.lo1 > 0.:
-#             asaplog.push("Using user defined LO1 frequency %e16.12 [Hz]" % self.lo1)
-#             lo1 = self.lo1
-#         else:
-#             print "NOT IMPLEMENTED YET!!!"
-
-#     def save(self, outfile, outform="ASAP", overwrite=False):
-#         if not overwrite and os.path.exists(outfile):
-#             raise RuntimeError, "Output file '%s' already exists" % outfile
-# 
-#         #self._separator._save(outfile, outform)
-
-#     def done(self):
-#         self.close()
-
-#     def close(self):
-#         pass
-#         #del self._separator
-    
-
-
-########################################################################
-    def _Deconvolution(self, data_array, shift, threshold=0.00000001):
-        FObs = []
-        Reject = 0
-        nshift, nchan = data_array.shape
-        nspec = nshift*(nshift-1)/2
-        ifftObs  = numpy.zeros((nspec, nchan), numpy.float)
-        for i in range(nshift):
-           F = FFT.fft(data_array[i])
-           FObs.append(F)
-        z = 0
-        for i in range(nshift):
-            for j in range(i+1, nshift):
-                Fobs = (FObs[i]+FObs[j])/2.0
-                dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
-                #print 'dX,i,j=',dX,i,j
-                for k in range(1,self.nchan):
-                    if math.fabs(math.sin(dX*k)) > threshold:
-                        Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
-                    else: Reject += 1
-                ifftObs[z] = FFT.ifft(Fobs)
-                z += 1
-        print 'Threshold=%s Reject=%d' % (threshold, Reject)
-        return ifftObs
-
-    def _combineResult(self, ifftObs):
-        nspec = len(ifftObs)
-        sum = ifftObs[0]
-        for i in range(1,nspec):
-            sum += ifftObs[i]
-        return(sum/float(nspec))
-
-    def _subtractOtherSide(self, data_array, shift, Data):
-        sum = numpy.zeros(len(Data), numpy.float)
-        numSP = len(data_array)
-        for i in range(numSP):
-            SPsub = data_array[i] - Data
-            sum += self._shiftSpectrum(SPsub, -shift[i])
-        return(sum/float(numSP))
-
-    def _shiftSpectrum(self, data, Shift):
-        Out = numpy.zeros(self.nchan, numpy.float)
-        w2 = Shift % 1
-        w1 = 1.0 - w2
-        for i in range(self.nchan):
-            c1 = int((Shift + i) % self.nchan)
-            c2 = (c1 + 1) % self.nchan
-            Out[c1] += data[i] * w1
-            Out[c2] += data[i] * w2
-        return Out.copy()
+#     def _preprocess_tables(self):
+#         ### temporary method to preprocess data
+#         ### Do time averaging for now.
+#         for itab in range(len(self.tables)):
+#             self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
+
+
+# ########################################################################
+#     def _Deconvolution(self, data_array, shift, threshold=0.00000001):
+#         FObs = []
+#         Reject = 0
+#         nshift, nchan = data_array.shape
+#         nspec = nshift*(nshift-1)/2
+#         ifftObs  = numpy.zeros((nspec, nchan), numpy.float)
+#         for i in range(nshift):
+#            F = FFT.fft(data_array[i])
+#            FObs.append(F)
+#         z = 0
+#         for i in range(nshift):
+#             for j in range(i+1, nshift):
+#                 Fobs = (FObs[i]+FObs[j])/2.0
+#                 dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
+#                 #print 'dX,i,j=',dX,i,j
+#                 for k in range(1,self.nchan):
+#                     if math.fabs(math.sin(dX*k)) > threshold:
+#                         Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
+#                     else: Reject += 1
+#                 ifftObs[z] = FFT.ifft(Fobs)
+#                 z += 1
+#         print 'Threshold=%s Reject=%d' % (threshold, Reject)
+#         return ifftObs
+
+#     def _combineResult(self, ifftObs):
+#         nspec = len(ifftObs)
+#         sum = ifftObs[0]
+#         for i in range(1,nspec):
+#             sum += ifftObs[i]
+#         return(sum/float(nspec))
+
+#     def _subtractOtherSide(self, data_array, shift, Data):
+#         sum = numpy.zeros(len(Data), numpy.float)
+#         numSP = len(data_array)
+#         for i in range(numSP):
+#             SPsub = data_array[i] - Data
+#             sum += self._shiftSpectrum(SPsub, -shift[i])
+#         return(sum/float(numSP))
+
+#     def _shiftSpectrum(self, data, Shift):
+#         Out = numpy.zeros(self.nchan, numpy.float)
+#         w2 = Shift % 1
+#         w1 = 1.0 - w2
+#         for i in range(self.nchan):
+#             c1 = int((Shift + i) % self.nchan)
+#             c2 = (c1 + 1) % self.nchan
+#             Out[c1] += data[i] * w1
+#             Out[c2] += data[i] * w2
+#         return Out.copy()

trunk/src/STSideBandSep.cpp

@@ -26 +26 @@
 
 // asap
+#include "STGrid.h"
+#include "STMath.h"
+#include "MathUtils.h"
 #include "STSideBandSep.h"
 
@@ -32 +35 @@
 using namespace asap ;
 
-#ifndef KS_DEBUG
-#define KS_DEBUG
-#endif
+// #ifndef KS_DEBUG
+// #define KS_DEBUG
+// #endif
 
 namespace asap {
@@ -78 +81 @@
 
   init();
+  tp_ = intabList_[0]->table().tableType();
 
   os << ntable_ << " tables are set." << LogIO::POST;
@@ -92 +96 @@
 {
   // frequency setup
-  sigIfno_= 0;
+  sigIfno_= -1;
   ftol_ = -1;
   solFrame_ = MFrequency::N_Types;
@@ -109 +113 @@
   // Default LO frame is TOPO
   loFrame_ = MFrequency::TOPO;
+  // scantable storage
+  tp_ = Table::Memory;
 };
 
@@ -130 +136 @@
 
   // IFNO
-  sigIfno_ = ifno;
+  sigIfno_ = (int) ifno;
 
   // Frequency tolerance
@@ -219 +225 @@
       if (i != imgShift_.size()-1) os << " , ";
     }
-    os << " ) [channel]" << LogIO::POST;
+    os << " ) [channels]" << LogIO::POST;
   }
 };
@@ -227 +233 @@
   LogIO os(LogOrigin("STSideBandSep","setThreshold()", WHERE));
   if (limit < 0)
-    throw( AipsError("Rejection limit should be positive number.") );
+    throw( AipsError("Rejection limit should be a positive number.") );
 
   rejlimit_ = limit;
-  os << "Rejection limit is set to " << rejlimit_;
-};
-
-// Temporal function to set scantable of image sideband
-void STSideBandSep::setImageTable(const ScantableWrapper &s)
-{
-#ifdef KS_DEBUG
-  cout << "STSideBandSep::setImageTable" << endl;
-  cout << "got image table nrow = " << s.nrow() << endl;
-#endif
-  imgTab_p = s.getCP();
-  AlwaysAssert(!imgTab_p.null(),AipsError);
-
-};
-
-//
-void STSideBandSep::setLO1(const double lo1, const string frame,
+  os << "Rejection limit is set to " << rejlimit_ << LogIO::POST;
+};
+
+void STSideBandSep::separate(string outname)
+{
+#ifdef KS_DEBUG
+  cout << "STSideBandSep::separate" << endl;
+#endif
+  LogIO os(LogOrigin("STSideBandSep","separate()", WHERE));
+  if (outname.empty())
+    outname = "sbseparated.asap";
+
+  // Set up a goup of IFNOs in the list of scantables within
+  // the frequency tolerance and make them a list.
+  nshift_ = setupShift();
+  if (nshift_ < 2)
+    throw( AipsError("At least 2 IFs are necessary for convolution.") );
+  // Grid scantable and generate output tables
+  ScantableWrapper gridst = gridTable();
+  sigTab_p = gridst.getCP();
+  if (doboth_)
+    imgTab_p = gridst.copy().getCP();
+  vector<unsigned int> remRowIds;
+  remRowIds.resize(0);
+  Matrix<float> specMat(nchan_, nshift_);
+  vector<float> sigSpec(nchan_), imgSpec(nchan_);
+  vector<uInt> tabIdvec;
+
+  //Generate FFTServer
+  fftsf.resize(IPosition(1, nchan_), FFTEnums::REALTOCOMPLEX);
+  fftsi.resize(IPosition(1, nchan_), FFTEnums::COMPLEXTOREAL);
+
+  /// Loop over sigTab_p and separate sideband
+  for (int irow = 0; irow < sigTab_p->nrow(); irow++){
+    tabIdvec.resize(0);
+    const int polId = sigTab_p->getPol(irow);
+    const int beamId = sigTab_p->getBeam(irow);
+    const vector<double> dir = sigTab_p->getDirectionVector(irow);
+    // Get a set of spectra to solve
+    if (!getSpectraToSolve(polId, beamId, dir[0], dir[1],
+                           specMat, tabIdvec)){
+      remRowIds.push_back(irow);
+#ifdef KS_DEBUG
+      cout << "no matching row found. skipping row = " << irow << endl;
+#endif
+      continue;
+    }
+    // Solve signal sideband
+    sigSpec = solve(specMat, tabIdvec, true);
+    sigTab_p->setSpectrum(sigSpec, irow);
+
+    // Solve image sideband
+    if (doboth_) {
+      imgSpec = solve(specMat, tabIdvec, false);
+      imgTab_p->setSpectrum(imgSpec, irow);
+    }
+  } // end of row loop
+
+  // Remove or flag rows without relevant data from gridded tables
+  if (remRowIds.size() > 0) {
+    const size_t nrem = remRowIds.size();
+    if ( sigTab_p->table().canRemoveRow() ) {
+      sigTab_p->table().removeRow(remRowIds);
+      os << "Removing " << nrem << " rows from the signal band table"
+         << LogIO::POST;
+    } else {
+      sigTab_p->flagRow(remRowIds, false);
+      os << "Cannot remove rows from the signal band table. Flagging "
+         << nrem << " rows" << LogIO::POST;
+    }
+
+    if (doboth_) {
+      if ( imgTab_p->table().canRemoveRow() ) {
+        imgTab_p->table().removeRow(remRowIds);
+        os << "Removing " << nrem << " rows from the image band table"
+           << LogIO::POST;
+      } else {
+        imgTab_p->flagRow(remRowIds, false);
+        os << "Cannot remove rows from the image band table. Flagging "
+           << nrem << " rows" << LogIO::POST;
+      }
+    }
+  }
+
+  // Finally, save tables on disk
+  if (outname.size() ==0)
+    outname = "sbseparated.asap";
+  const string sigName = outname + ".signalband";
+  os << "Saving SIGNAL sideband table: " << sigName << LogIO::POST;
+  sigTab_p->makePersistent(sigName);
+  if (doboth_) {
+    solveImageFrequency();
+    const string imgName = outname + ".imageband";
+    os << "Saving IMAGE sideband table: " << sigName << LogIO::POST;
+    imgTab_p->makePersistent(imgName);
+  }
+
+};
+
+unsigned int STSideBandSep::setupShift()
+{
+  LogIO os(LogOrigin("STSideBandSep","setupShift()", WHERE));
+  if (infileList_.size() == 0 && intabList_.size() == 0)
+    throw( AipsError("No scantable has been set. Set a list of scantables first.") );
+
+  const bool byname = (intabList_.size() == 0);
+  // Make sure sigIfno_ exists in the first table.
+  CountedPtr<Scantable> stab;
+  vector<string> coordsav;
+  vector<string> coordtmp(3);
+  os << "Checking IFNO in the first table." << LogIO::POST;
+  if (byname) {
+    if (!checkFile(infileList_[0], "d"))
+      os << LogIO::SEVERE << "Could not find scantable '" << infileList_[0]
+         << "'" << LogIO::EXCEPTION;
+    stab = CountedPtr<Scantable>(new Scantable(infileList_[0]));
+  } else {
+    stab = intabList_[0];
+  }
+  if (sigIfno_ < 0) {
+    sigIfno_ = (int) stab->getIF(0);
+    os << "IFNO to process has not been set. Using the first IF = "
+       << sigIfno_ << LogIO::POST;
+  }
+
+  unsigned int basench;
+  double basech0, baseinc, ftolval, inctolval;
+  coordsav = stab->getCoordInfo();
+  const string stfframe = coordsav[1];
+  coordtmp[0] = "Hz";
+  coordtmp[1] = ( (solFrame_ == MFrequency::N_Types) ?
+                  stfframe :
+                  MFrequency::showType(solFrame_) );
+  coordtmp[2] = coordsav[2];
+  stab->setCoordInfo(coordtmp);
+  if (!getFreqInfo(stab, (unsigned int) sigIfno_, basech0, baseinc, basench)) {
+    os << LogIO::SEVERE << "No data with IFNO=" << sigIfno_
+       << " found in the first table." << LogIO::EXCEPTION;
+  }
+  else {
+    os << "Found IFNO = " << sigIfno_
+       << " in the first table." << LogIO::POST;
+  }
+  if (ftol_.getUnit().empty()) {
+    // tolerance in unit of channels
+    ftolval = ftol_.getValue() * baseinc;
+  }
+  else {
+    ftolval = ftol_.getValue("Hz");
+  }
+  inctolval = abs(baseinc/(double) basench);
+  const string poltype0 = stab->getPolType();
+
+  // Initialize shift values
+  initshift();
+
+  const bool setImg = ( doboth_ && (imgShift_.size() == 0) );
+  // Select IFs
+  for (unsigned int itab = 0; itab < ntable_; itab++ ){
+    os << "Table " << itab << LogIO::POST;
+    if (itab > 0) {
+      if (byname) {
+        if (!checkFile(infileList_[itab], "d"))
+          os << LogIO::SEVERE << "Could not find scantable '"
+             << infileList_[itab] << "'" << LogIO::EXCEPTION;
+        stab = CountedPtr<Scantable>(new Scantable(infileList_[itab]));
+      } else {
+        stab = intabList_[itab];
+      }
+      //POLTYPE should be the same.
+      if (stab->getPolType() != poltype0 ) {
+        os << LogIO::WARN << "POLTYPE differs from the first table."
+           << " Skipping the table" << LogIO::POST;
+        continue;
+      }
+      // Multi beam data may not handled properly
+      if (stab->nbeam() > 1)
+        os <<  LogIO::WARN << "Table contains multiple beams. "
+           << "It may not be handled properly."  << LogIO::POST;
+
+      coordsav = stab->getCoordInfo();
+      coordtmp[2] = coordsav[2];
+      stab->setCoordInfo(coordtmp);
+    }
+    bool selected = false;
+    vector<uint> ifnos = stab->getIFNos();
+    vector<uint>::iterator iter;
+    const STSelector& basesel = stab->getSelection();
+    for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
+      unsigned int nch;
+      double freq0, incr;
+      if ( getFreqInfo(stab, *iter, freq0, incr, nch) ){
+        if ( (nch == basench) && (abs(freq0-basech0) < ftolval)
+             && (abs(incr-baseinc) < inctolval) ){
+          //Found
+          STSelector sel(basesel);
+          sel.setIFs(vector<int>(1,(int) *iter));
+          stab->setSelection(sel);
+          CountedPtr<Scantable> seltab = ( new Scantable(*stab, false) );
+          stab->setSelection(basesel);
+          seltab->setCoordInfo(coordsav);
+          const double chShift = (freq0 - basech0) / baseinc;
+          tableList_.push_back(seltab);
+          sigShift_.push_back(chShift);
+          if (setImg)
+            imgShift_.push_back(-chShift);
+
+          selected = true;
+          os << "- IF" << *iter << " selected: sideband shift = "
+             << chShift << " channels" << LogIO::POST;
+        }
+      }
+    } // ifno loop
+    stab->setCoordInfo(coordsav);
+    if (!selected)
+      os << LogIO::WARN << "No data selected in Table "
+         << itab << LogIO::POST;
+  } // table loop
+  nchan_ = basench;
+
+  os << "Total number of IFs selected = " << tableList_.size()
+     << LogIO::POST;
+  if ( setImg && (sigShift_.size() != imgShift_.size()) ){
+      os << LogIO::SEVERE
+         << "User defined channel shift of image sideband has "
+         << imgShift_.size()
+         << " elements, while selected IFNOs are " << sigShift_.size()
+         << "\nThe frequency tolerance (freqtol) may be too small."
+         << LogIO::EXCEPTION;
+  }
+
+  return tableList_.size();
+};
+
+bool STSideBandSep::getFreqInfo(const CountedPtr<Scantable> &stab,
+                                const unsigned int &ifno,
+                                double &freq0, double &incr,
+                                unsigned int &nchan)
+{
+    vector<uint> ifnos = stab->getIFNos();
+    bool found = false;
+    vector<uint>::iterator iter;
+    for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
+      if (*iter == ifno) {
+        found = true;
+        break;
+      }
+    }
+    if (!found)
+      return false;
+
+    const STSelector& basesel = stab->getSelection();
+    STSelector sel(basesel);
+    sel.setIFs(vector<int>(1,(int) ifno));
+    stab->setSelection(sel);
+    vector<double> freqs;
+    freqs = stab->getAbcissa(0);
+    freq0 = freqs[0];
+    incr = freqs[1] - freqs[0];
+    nchan = freqs.size();
+    stab->setSelection(basesel);
+    return true;
+};
+
+ScantableWrapper STSideBandSep::gridTable()
+{
+  LogIO os(LogOrigin("STSideBandSep","gridTable()", WHERE));
+  if (tableList_.size() == 0)
+    throw( AipsError("Internal error. No scantable has been set to grid.") );
+  Double xmin, xmax, ymin, ymax;
+  mapExtent(tableList_, xmin, xmax, ymin, ymax);
+  const Double centx = 0.5 * (xmin + xmax);
+  const Double centy = 0.5 * (ymin + ymax);
+  const int nx = max(1, (int) ceil( (xmax - xmin) * cos(centy) /xtol_ ) );
+  const int ny = max(1, (int) ceil( (ymax - ymin) / ytol_ ) );
+
+  string scellx, scelly;
+  {
+    ostringstream oss;
+    oss << xtol_ << "rad" ;
+    scellx = oss.str();
+  }
+  {
+    ostringstream oss;
+    oss << ytol_ << "rad" ;
+    scelly = oss.str();
+  }
+
+  ScantableWrapper stab0;
+  if (intabList_.size() > 0)
+    stab0 = ScantableWrapper(intabList_[0]);
+  else
+    stab0 = ScantableWrapper(infileList_[0]);
+
+  string scenter;
+  {
+    ostringstream oss;
+    oss << stab0.getCP()->getDirectionRefString() << " "
+        << centx << "rad" << " " << centy << "rad";
+    scenter = oss.str();
+  } 
+  
+  STGrid2 gridder = STGrid2(stab0);
+  gridder.setIF(sigIfno_);
+  gridder.defineImage(nx, ny, scellx, scelly, scenter);
+  gridder.setFunc("box", 1);
+  gridder.setWeight("uniform");
+#ifdef KS_DEBUG
+  cout << "Grid parameter summary: " << endl;
+  cout << "- IF = " << sigIfno_ << endl;
+  cout << "- center = " << scenter << ")\n"
+       << "- npix = (" << nx << ", " << ny << ")\n"
+       << "- cell = (" << scellx << ", " << scelly << endl;
+#endif
+  gridder.grid();
+  const int itp = (tp_ == Table::Memory ? 0 : 1);
+  ScantableWrapper gtab = gridder.getResultAsScantable(itp);
+  return gtab;
+};
+
+void STSideBandSep::mapExtent(vector< CountedPtr<Scantable> > &tablist,
+                              Double &xmin, Double &xmax, 
+                              Double &ymin, Double &ymax)
+{
+  ROArrayColumn<Double> dirCol_;
+  dirCol_.attach( tablist[0]->table(), "DIRECTION" );
+  Matrix<Double> direction = dirCol_.getColumn();
+  Vector<Double> ra( direction.row(0) );
+  mathutil::rotateRA(ra);
+  minMax( xmin, xmax, ra );
+  minMax( ymin, ymax, direction.row(1) );
+  Double amin, amax, bmin, bmax;
+  const uInt ntab = tablist.size();
+  for ( uInt i = 1 ; i < ntab ; i++ ) {
+    dirCol_.attach( tablist[i]->table(), "DIRECTION" );
+    direction.assign( dirCol_.getColumn() );
+    ra.assign( direction.row(0) );
+    mathutil::rotateRA(ra);
+    minMax( amin, amax, ra );
+    minMax( bmin, bmax, direction.row(1) );
+    xmin = min(xmin, amin);
+    xmax = max(xmax, amax);
+    ymin = min(ymin, bmin);
+    ymax = max(ymax, bmax);
+  }
+};
+
+bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
+                                      const double dirX, const double dirY,
+                                      Matrix<float> &specMat, vector<uInt> &tabIdvec)
+{
+  LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE));
+
+  tabIdvec.resize(0);
+  specMat.resize(nchan_, nshift_);
+  Vector<float> spec;
+  uInt nspec = 0;
+  STMath stm(false); // insitu has no effect for average.
+  for (uInt itab = 0 ; itab < nshift_ ; itab++) {
+    CountedPtr<Scantable> currtab_p = tableList_[itab];
+    // Selection by POLNO and BEAMNO
+    const STSelector& basesel = currtab_p->getSelection();
+    STSelector sel(basesel);
+    sel.setPolarizations(vector<int>(1, polId));
+    sel.setBeams(vector<int>(1, beamId));
+    try {
+      currtab_p->setSelection(sel);
+    } catch (...) {
+#ifdef KS_DEBUG
+      cout << "Table " << itab << " - No spectrum found. skipping the table."
+           << endl;
+#endif
+      continue;
+    }
+    // Selection by direction;
+    vector<int> selrow(0);
+    vector<double> currDir(2, 0.);
+    const int nselrow = currtab_p->nrow();
+    for (int irow = 0 ; irow < nselrow ; irow++) {
+      currDir = currtab_p->getDirectionVector(irow);
+      if ( (abs(currDir[0]-dirX) > xtol_) ||
+           (abs(currDir[1]-dirY) > ytol_) )
+        continue;
+      // within direction tolerance
+      selrow.push_back(irow);
+    } // end of row loop
+
+    if (selrow.size() < 1) {
+      currtab_p->setSelection(basesel);
+
+#ifdef KS_DEBUG
+      cout << "Table " << itab << " - No spectrum found. skipping the table."
+           << endl;
+#endif
+
+      continue;
+    }
+
+    // At least a spectrum is selected in this table
+    CountedPtr<Scantable> seltab_p = ( new Scantable(*currtab_p, false) );
+    currtab_p->setSelection(basesel);
+    STSelector sel2(seltab_p->getSelection());
+    sel2.setRows(selrow);
+    seltab_p->setSelection(sel2);
+    CountedPtr<Scantable> avetab_p;
+    vector<bool> mask;
+    if (seltab_p->nrow() > 1) {
+      avetab_p = stm.average(vector< CountedPtr<Scantable> >(1, seltab_p),
+                             vector<bool>(), "TINTSYS", "NONE");
+#ifdef KS_DEBUG
+      cout << "Table " << itab << " - more than a spectrum is selected. averaging rows..."
+           << endl;
+#endif
+      if (avetab_p->nrow() > 1)
+        throw( AipsError("Averaged table has more than a row. Somethigs went wrong.") );
+    } else {
+      avetab_p = seltab_p;
+    }
+    spec.reference(specMat.column(nspec));
+    spec = avetab_p->getSpectrum(0);
+    tabIdvec.push_back((uInt) itab);
+    nspec++;
+  } // end of table loop
+  if (nspec != nshift_){
+    //shiftSpecmat.resize(nchan_, nspec, true);
+#ifdef KS_DEBUG
+      cout << "Could not find corresponding rows in some tables."
+           << endl;
+      cout << "Number of spectra selected = " << nspec << endl;
+#endif
+  }
+  if (nspec < 2) {
+#ifdef KS_DEBUG
+      cout << "At least 2 spectra are necessary for convolution"
+           << endl;
+#endif
+      return false;
+  }
+  return true;
+};
+
+vector<float> STSideBandSep::solve(const Matrix<float> &specmat,
+                                   const vector<uInt> &tabIdvec,
+                                   const bool signal)
+{
+  LogIO os(LogOrigin("STSideBandSep","solve()", WHERE));
+  if (tabIdvec.size() == 0)
+    throw(AipsError("Internal error. Table index is not defined."));
+  if (specmat.ncolumn() != tabIdvec.size())
+    throw(AipsError("Internal error. The row number of input matrix is not conformant."));
+  if (specmat.nrow() != nchan_)
+    throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
+  
+
+#ifdef KS_DEBUG
+  cout << "Solving " << (signal ? "SIGNAL" : "IMAGE") << " sideband."
+     << endl;
+#endif
+
+  const size_t nspec = tabIdvec.size();
+  vector<double> *thisShift, *otherShift;
+  if (signal == otherside_) {
+    // (solve signal && solveother = T) OR (solve image && solveother = F)
+    thisShift = &imgShift_;
+    otherShift = &sigShift_;
+#ifdef KS_DEBUG
+    cout << "Image sideband will be deconvolved." << endl;
+#endif
+  } else {
+    // (solve signal && solveother = F) OR (solve image && solveother = T)
+    thisShift =  &sigShift_;
+    otherShift = &imgShift_; 
+#ifdef KS_DEBUG
+    cout << "Signal sideband will be deconvolved." << endl;
+#endif
+ }
+
+  vector<double> spshift(nspec);
+  Matrix<float> shiftSpecmat(nchan_, nspec, 0.);
+  double tempshift;
+  Vector<float> shiftspvec;
+  for (uInt i = 0 ; i < nspec; i++) {
+    spshift[i] = otherShift->at(i) - thisShift->at(i);
+    tempshift = - thisShift->at(i);
+    shiftspvec.reference(shiftSpecmat.column(i));
+    shiftSpectrum(specmat.column(i), tempshift, shiftspvec);
+  }
+
+  Matrix<float> convmat(nchan_, nspec*(nspec-1)/2, 0.);
+  vector<float> thisvec(nchan_, 0.);
+
+  float minval, maxval;
+  minMax(minval, maxval, shiftSpecmat);
+#ifdef KS_DEBUG
+  cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
+#endif
+
+#ifdef KS_DEBUG
+  cout << "starting deconvolution" << endl;
+#endif
+  deconvolve(shiftSpecmat, spshift, rejlimit_, convmat);
+#ifdef KS_DEBUG
+  cout << "finished deconvolution" << endl;
+#endif
+
+  minMax(minval, maxval, convmat);
+#ifdef KS_DEBUG
+  cout << "Max/Min of output Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
+#endif
+
+  aggregateMat(convmat, thisvec);
+
+  if (!otherside_) return thisvec;
+
+  // subtract from the other side band.
+  vector<float> othervec(nchan_, 0.);
+  subtractFromOther(shiftSpecmat, thisvec, spshift, othervec);
+  return othervec;
+};
+
+
+void STSideBandSep::shiftSpectrum(const Vector<float> &invec,
+                                  double shift,
+                                  Vector<float> &outvec)
+{
+  LogIO os(LogOrigin("STSideBandSep","shiftSpectrum()", WHERE));
+  if (invec.size() != nchan_)
+    throw(AipsError("Internal error. The length of input vector differs from nchan_"));
+  if (outvec.size() != nchan_)
+    throw(AipsError("Internal error. The length of output vector differs from nchan_"));
+
+#ifdef KS_DEBUG
+  cout << "Start shifting spectrum for " << shift << "channels" << endl;
+#endif
+
+  // tweak shift to be in 0 ~ nchan_-1
+  if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
+  if (shift < 0.) shift += nchan_;
+  double rweight = fmod(shift, 1.);
+  if (rweight < 0.) rweight += 1.;
+  double lweight = 1. - rweight;
+  uInt lchan, rchan;
+
+  outvec = 0;
+  for (uInt i = 0 ; i < nchan_ ; i++) {
+    lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
+    if (lchan < 0.) lchan += nchan_;
+    rchan = ( (lchan + 1) % nchan_ );
+    outvec(lchan) += invec(i) * lweight;
+    outvec(rchan) += invec(i) * rweight;
+  }
+};
+
+void STSideBandSep::deconvolve(Matrix<float> &specmat,
+                               const vector<double> shiftvec,
+                               const double threshold,
+                               Matrix<float> &outmat)
+{
+  LogIO os(LogOrigin("STSideBandSep","deconvolve()", WHERE));
+  if (specmat.nrow() != nchan_)
+    throw(AipsError("Internal error. The length of input matrix differs from nchan_"));
+  if (specmat.ncolumn() != shiftvec.size())
+    throw(AipsError("Internal error. The number of input shifts and spectrum  differs."));
+
+  float minval, maxval;
+#ifdef KS_DEBUG
+  minMax(minval, maxval, specmat);
+  cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
+#endif
+
+  uInt ninsp = shiftvec.size();
+  outmat.resize(nchan_, ninsp*(ninsp-1)/2, 0.);
+  Matrix<Complex> fftspmat(nchan_/2+1, ninsp, 0.);
+  Vector<float> rvecref(nchan_, 0.);
+  Vector<Complex> cvecref(nchan_/2+1, 0.);
+  uInt icol = 0;
+  unsigned int nreject = 0;
+
+#ifdef KS_DEBUG
+  cout << "Starting initial FFT. The number of input spectra = " << ninsp << endl;
+  cout << "out matrix has ncolumn = " << outmat.ncolumn() << endl;
+#endif
+
+  for (uInt isp = 0 ; isp < ninsp ; isp++) {
+    rvecref.reference( specmat.column(isp) );
+    cvecref.reference( fftspmat.column(isp) );
+
+#ifdef KS_DEBUG
+    minMax(minval, maxval, rvecref);
+    cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
+#endif
+
+    fftsf.fft0(cvecref, rvecref, true);
+
+#ifdef KS_DEBUG
+    double maxr=cvecref[0].real(), minr=cvecref[0].real(),
+      maxi=cvecref[0].imag(), mini=cvecref[0].imag();
+    for (uInt i = 1; i<cvecref.size();i++){
+      maxr = max(maxr, cvecref[i].real());
+      maxi = max(maxi, cvecref[i].imag());
+      minr = min(minr, cvecref[i].real());
+      mini = min(mini, cvecref[i].imag());
+    }
+    cout << "Max/Min of inv FFTed Matrix (size=" << cvecref.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
+#endif
+  }
+
+  //Liberate from reference
+  rvecref.unique();
+
+  Vector<Complex> cspec(nchan_/2+1, 0.);
+  const double PI = 6.0 * asin(0.5);
+  const double nchani = 1. / (float) nchan_;
+  const Complex trans(0., 1.);
+#ifdef KS_DEBUG
+  cout << "starting actual deconvolution" << endl;
+#endif
+  for (uInt j = 0 ; j < ninsp ; j++) {
+    for (uInt k = j+1 ; k < ninsp ; k++) {
+      const double dx = (shiftvec[k] - shiftvec[j]) * 2. * PI * nchani;
+
+#ifdef KS_DEBUG
+      cout << "icol = " << icol << endl;
+#endif
+
+      for (uInt ichan = 0 ; ichan < cspec.size() ; ichan++){
+        cspec[ichan] = ( fftspmat(ichan, j) + fftspmat(ichan, k) )*0.5;
+        double phase = dx*ichan;
+        if ( fabs( sin(phase) ) > threshold){
+          cspec[ichan] += ( fftspmat(ichan, j) - fftspmat(ichan, k) ) * 0.5
+            * trans * sin(phase) / ( 1. - cos(phase) );
+        } else {
+          nreject++;
+        }
+      } // end of channel loop
+
+#ifdef KS_DEBUG
+      cout << "done calculation of cspec" << endl;
+#endif
+
+      Vector<Float> rspec;
+      rspec.reference( outmat.column(icol) );
+
+#ifdef KS_DEBUG
+      cout << "Starting inverse FFT. icol = " << icol << endl;
+      //cout << "- size of complex vector = " << cspec.size() << endl;
+      double maxr=cspec[0].real(), minr=cspec[0].real(),
+        maxi=cspec[0].imag(), mini=cspec[0].imag();
+      for (uInt i = 1; i<cspec.size();i++){
+        maxr = max(maxr, cspec[i].real());
+        maxi = max(maxi, cspec[i].imag());
+        minr = min(minr, cspec[i].real());
+        mini = min(mini, cspec[i].imag());
+      }
+      cout << "Max/Min of conv vector (size=" << cspec.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
+#endif
+
+      fftsi.fft0(rspec, cspec, false);
+
+#ifdef KS_DEBUG
+      //cout << "- size of inversed real vector = " << rspec.size() << endl;
+      minMax(minval, maxval, rspec);
+      cout << "Max/Min of inv FFTed Vector (size=" << rspec.size() << ") = (max: " << maxval << ", min: " << minval << ")" << endl;
+      //cout << "Done inverse FFT. icol = " << icol << endl;
+#endif
+
+      icol++;
+    }
+  }
+
+#ifdef KS_DEBUG
+  minMax(minval, maxval, outmat);
+  cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
+#endif
+
+  os << "Threshold = " << threshold << ", Rejected channels = " << nreject << endl;
+};
+
+////////////////////////////////////////////////////////////////////
+// void STSideBandSep::cpprfft(std::vector<float> invec)
+// {
+//   cout << "c++ method cpprfft" << endl;
+//   const unsigned int len = invec.size();
+//   Vector<Complex> carr(len/2+1, 0.);
+//   Vector<float> inarr = Vector<float>(invec);
+//   Vector <float> outarr(len, 0.);
+//   FFTServer<Float, Complex> fftsf, fftsi;
+//   fftsf.resize(IPosition(1, len), FFTEnums::REALTOCOMPLEX);
+//   fftsi.resize(IPosition(1, invec.size()), FFTEnums::COMPLEXTOREAL);
+//   cout << "Input float array (length = " << len << "):" << endl;
+//   for (uInt i = 0 ; i < len ; i++){
+//     cout << (i == 0 ? "( " : " ") << inarr[i] << (i == len-1 ? ")" : ",");
+//   }
+//   cout << endl;
+//   cout << "R->C transform" << endl;
+//   fftsf.fft0(carr, inarr, true);
+//   cout << "FFTed complex array (length = " << carr.size() << "):" << endl;
+//   for (uInt i = 0 ; i < carr.size() ; i++){
+//     cout << (i == 0 ? "( " : " ") << carr[i] << ( (i == carr.size()-1) ? ")" : "," );
+//   }
+//   cout << endl;
+//   cout << "C->R transform" << endl;
+//   fftsi.fft0(outarr, carr, false);
+//   cout << "invFFTed float array (length = " << outarr.size() << "):" << endl;
+//   for (uInt i = 0 ; i < outarr.size() ; i++){
+//     cout << (i == 0 ? "( " : " ") << outarr[i] << ( (i == outarr.size()-1) ? ")" : "," );
+//   }
+//   cout << endl;
+// };
+////////////////////////////////////////////////////////////////////
+
+
+void STSideBandSep::aggregateMat(Matrix<float> &inmat,
+                                 vector<float> &outvec)
+{
+  LogIO os(LogOrigin("STSideBandSep","aggregateMat()", WHERE));
+  if (inmat.nrow() != nchan_)
+    throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
+//   if (outvec.size() != nchan_)
+//     throw(AipsError("Internal error. The size of output vector should be equal to nchan_"));
+
+  os << "Averaging " << inmat.ncolumn() << " spectra in the input matrix."
+     << LogIO::POST;
+
+  const uInt nspec = inmat.ncolumn();
+  const double scale = 1./( (double) nspec );
+  // initialize values with 0
+  outvec.assign(nchan_, 0);
+  for (uInt isp = 0 ; isp < nspec ; isp++) {
+    for (uInt ich = 0 ; ich < nchan_ ; ich++) {
+      outvec[ich] += inmat(ich, isp);
+    }
+  }
+
+  vector<float>::iterator iter;
+  for (iter = outvec.begin(); iter != outvec.end(); iter++){
+    *iter *= scale;
+  }
+};
+
+void STSideBandSep::subtractFromOther(const Matrix<float> &shiftmat,
+                                      const vector<float> &invec,
+                                      const vector<double> &shift,
+                                      vector<float> &outvec)
+{
+  LogIO os(LogOrigin("STSideBandSep","subtractFromOther()", WHERE));
+  if (shiftmat.nrow() != nchan_)
+    throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
+  if (invec.size() != nchan_)
+    throw(AipsError("Internal error. The length of input vector should be nchan_"));
+  if (shift.size() != shiftmat.ncolumn())
+    throw(AipsError("Internal error. The column numbers of input matrix != the number of elements in shift"));
+
+  const uInt nspec = shiftmat.ncolumn();
+  Vector<float> subsp(nchan_, 0.), shiftsub;
+  Matrix<float> submat(nchan_, nspec, 0.);
+  vector<float>::iterator iter;
+  for (uInt isp = 0 ; isp < nspec ; isp++) {
+    for (uInt ich = 0; ich < nchan_ ; ich++) {
+      subsp(ich) = shiftmat(ich, isp) - invec[ich];
+    }
+    shiftsub.reference(submat.column(isp));
+    shiftSpectrum(subsp, shift[isp], shiftsub);
+  }
+
+  aggregateMat(submat, outvec);
+};
+
+
+void STSideBandSep::setLO1(const string lo1, const string frame,
                            const double reftime, const string refdir)
 {
-  lo1Freq_ = lo1;
+  Quantum<Double> qfreq;
+  readQuantity(qfreq, String(lo1));
+  lo1Freq_ = qfreq.getValue("Hz");
   MFrequency::getType(loFrame_, frame);
   loTime_ = reftime;
@@ -297 +1058 @@
 };
 
-///// TEMPORAL FUNCTION!!! /////
-void STSideBandSep::setScanTb0(const ScantableWrapper &s){
-  st0_ = s.getCP();
-};
-////////////////////////////////
-
-void STSideBandSep::solveImageFreqency()
-{
-#ifdef KS_DEBUG
-  cout << "STSideBandSep::solveImageFrequency" << endl;
-#endif
+
+void STSideBandSep::solveImageFrequency()
+{
   LogIO os(LogOrigin("STSideBandSep","solveImageFreqency()", WHERE));
   os << "Start calculating frequencies of image side band" << LogIO::POST;
@@ -383 +1136 @@
     sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue();
 
-
   // Check for the availability of LO1
   if (lo1Freq_ > 0.) {
@@ -402 +1154 @@
     // Try getting ASDM name from scantable header
     os << "Try getting information from scantable header" << LogIO::POST;
-    if (!getLo1FromScanTab(st0_, sigrefval, refpix, increment, nChan)) {
+    if (!getLo1FromScanTab(tableList_[0], sigrefval, refpix, increment, nChan)) {
       //throw AipsError("Failed to get LO1 frequency from asis table");
       os << LogIO::WARN << "Failed to get LO1 frequency using information in scantable." << LogIO::POST;
@@ -410 +1162 @@
     }
   }
+
   // LO1 should now be ready.
   if (lo1Freq_ < 0.)
     throw(AipsError("Got negative LO1 Frequency"));
 
-  // Print summary (LO1)
-  Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue();
-  os << "[LO1 settings]" << LogIO::POST;
-  os << "- Frequency: " << lo1Freq_ << " [Hz] ("
-     << MFrequency::showType(loFrame_) << ")" << LogIO::POST;
-  os << "- Reference time: " << me.get(Unit(String("d"))).getValue()
-     << " [day]" << LogIO::POST;
-  os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1]
-     << "] (" << md.getRefString() << ") " << LogIO::POST;
-
-  //Print summary (signal)
-  os << "[Signal side band]" << LogIO::POST;
-  os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")"
-     << LogIO::POST;
-  os << "- Reference value: " << refval << " [Hz] ("
-     << MFrequency::showType(tabframe) << ") = "
-     << sigrefval << " [Hz] (" <<  MFrequency::showType(loFrame_)
-     << ")" << LogIO::POST;
-  os << "- Reference pixel: " << refpix  << LogIO::POST;
-  os << "- Increment: " << increment << " [Hz]" << LogIO::POST;
+  // Print summary
+  {
+    // LO1
+    Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue();
+    os << "[LO1 settings]" << LogIO::POST;
+    os << "- Frequency: " << lo1Freq_ << " [Hz] ("
+       << MFrequency::showType(loFrame_) << ")" << LogIO::POST;
+    os << "- Reference time: " << me.get(Unit(String("d"))).getValue()
+       << " [day]" << LogIO::POST;
+    os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1]
+       << "] (" << md.getRefString() << ") " << LogIO::POST;
+
+    // signal sideband
+    os << "[Signal side band]" << LogIO::POST;
+    os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")"
+       << LogIO::POST;
+    os << "- Reference value: " << refval << " [Hz] ("
+       << MFrequency::showType(tabframe) << ") = "
+       << sigrefval << " [Hz] (" <<  MFrequency::showType(loFrame_)
+       << ")" << LogIO::POST;
+    os << "- Reference pixel: " << refpix  << LogIO::POST;
+    os << "- Increment: " << increment << " [Hz]" << LogIO::POST;
+  }
 
   // Calculate image band incr and refval in loFrame_
@@ -447 +1203 @@
   freqIdVec = fIDnew;
 
-  // Print summary (Image side band)
-  os << "[Image side band]" << LogIO::POST;
-  os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0)
-     << ")" << LogIO::POST;
-  os << "- Reference value: " << imgrefval << " [Hz] ("
-     << MFrequency::showType(tabframe) << ") = "
-     << imgrefval_tab << " [Hz] (" <<  MFrequency::showType(loFrame_)
-     << ")" << LogIO::POST;
-  os << "- Reference pixel: " << refpix  << LogIO::POST;
-  os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST;
+  // Print summary (Image sideband)
+  {
+    os << "[Image side band]" << LogIO::POST;
+    os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0)
+       << ")" << LogIO::POST;
+    os << "- Reference value: " << imgrefval << " [Hz] ("
+       << MFrequency::showType(tabframe) << ") = "
+       << imgrefval_tab << " [Hz] (" <<  MFrequency::showType(loFrame_)
+       << ")" << LogIO::POST;
+    os << "- Reference pixel: " << refpix  << LogIO::POST;
+    os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST;
+  }
 };
 
@@ -619 +1377 @@
 };
 
-// String STSideBandSep::()
-// {
-
-// };
-
 } //namespace asap

trunk/src/STSideBandSep.h

@@ -22 +22 @@
 #include <coordinates/Coordinates/DirectionCoordinate.h>
 #include <coordinates/Coordinates/SpectralCoordinate.h>
+#include <scimath/Mathematics/FFTServer.h>
 // asap
 #include "ScantableWrapper.h"
@@ -40 +41 @@
   explicit STSideBandSep(const vector<ScantableWrapper> &tables);
   virtual ~STSideBandSep();
+
+  ///////////// temp functions //////////////////////
+  //void cpprfft(std::vector<float> invec);
+  //////////////////////////////////////////////////
+
+  /**
+   * Separate side bands
+   **/
+  void separate(string outname);
 
   /**
@@ -83 +93 @@
    * Set additional information to fill frequencies of image sideband
    **/
-  void setLO1(const double lo1, const string frame="TOPO",
+  void setLO1(const string lo1, const string frame="TOPO",
               const double reftime=-1, string refdir="");
   void setLO1Root(const string name);
-
-  /**
-   * Actual calculation of frequencies of image sideband
-   **/
-  void solveImageFreqency();
 
 private:
@@ -99 +104 @@
   /** Return if the path exists (optionally, check file type) **/
   Bool checkFile(const string name, string type="");
+
+  /** **/
+  unsigned int setupShift();
+  bool getFreqInfo(const CountedPtr<Scantable> &stab, const unsigned int &ifno,
+                   double &freq0, double &incr, unsigned int &nchan);
+
+  /** Grid scantable **/
+  ScantableWrapper gridTable();
+  void mapExtent(vector< CountedPtr<Scantable> > &tablist,
+                 Double &xmin, Double &xmax,
+                 Double &ymin, Double &ymax);
+
+  /**
+   * Actual calculation of frequencies of image sideband
+   **/
+  void solveImageFrequency();
 
   /** 
@@ -112 +133 @@
                          const double refval, const double refpix,
                          const double increment, const int nChan);
+  bool getSpectraToSolve(const int polId, const int beamId,
+                         const double dirX, const double dirY,
+                         Matrix<float> &specmat, vector<uInt> &tabIdvec);
+
+  vector<float> solve(const Matrix<float> &specmat,
+                      const vector<uInt> &tabIdvec,
+                      const bool signal = true);
+
+  void shiftSpectrum(const Vector<float> &invec, double shift,
+                     Vector<float> &outvec);
+
+  void deconvolve(Matrix<float> &specmat, const vector<double> shiftvec,
+                  const double threshold, Matrix<float> &outmat);
+
+  void aggregateMat(Matrix<float> &inmat, vector<float> &outvec);
+
+  void subtractFromOther(const Matrix<float> &shiftmat,
+                         const vector<float> &invec,
+                         const vector<double> &shift,
+                         vector<float> &outvec);
+
+
 
   /** Member variables **/
@@ -119 +162 @@
   unsigned int ntable_;
   // frequency and direction setup to select data.
-  unsigned int sigIfno_;
+  int sigIfno_;
   Quantum<Double> ftol_;
   MFrequency::Types solFrame_;
@@ -130 +173 @@
   double rejlimit_;
   // LO1
-  double lo1Freq_;
+  double lo1Freq_; // in Hz
   MFrequency::Types loFrame_;
   double loTime_;
@@ -136 +179 @@
   string asdmName_, asisName_;
 
+  //CountedPtr<Scantable> imgTab_p, sigTab_p;
   CountedPtr<Scantable> imgTab_p, sigTab_p;
+  Table::TableType tp_;
+  FFTServer<Float, Complex> fftsf, fftsi;
   // TEMPORAL member
   CountedPtr<Scantable> st0_;

trunk/src/python_STSideBandSep.cpp

@@ -34 +34 @@
            boost::python::arg("refdir")="") )
     .def( "set_lo1root", &STSideBandSep::setLO1Root )
-    // temporal methods
-    .def( "set_imgtable", &STSideBandSep::setImageTable )
-    .def( "solve_imgfreq", &STSideBandSep::solveImageFreqency )
-    .def( "_get_asistb_from_scantb", &STSideBandSep::setScanTb0 )
+    .def( "separate", &STSideBandSep::separate )
+    //// temporal methods
+    //.def( "_cpprfft", &STSideBandSep::cpprfft )
+    //.def( "solve_imgfreq", &STSideBandSep::solveImageFreqency )
+    //.def( "_get_asistb_from_scantb", &STSideBandSep::setScanTb0 )
   ;
 };