Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/python/sbseparator.py

    r2859 r2685  
    99from asap.selector import selector
    1010from asap.asapgrid import asapgrid2
    11 from asap._asap import SBSeparator
     11#from asap._asap import sidebandsep
    1212
    1313class sbseparator:
     
    2121      information in scantable in sideband sparation. Frequency
    2222      setting of SIGNAL sideband is stored in output table for now.
     23    - Flag information (stored in FLAGTRA) is ignored.
    2324
    2425    Example:
     
    3738    """
    3839    def __init__(self, infiles):
    39         self._separator = SBSeparator(infiles)
    40 
    41 
     40        self.intables = None
     41        self.signalShift = []
     42        self.imageShift = []
     43        self.dsbmode = True
     44        self.getboth = False
     45        self.rejlimit = 0.2
     46        self.baseif = -1
     47        self.freqtol = 10.
     48        self.freqframe = ""
     49        self.solveother = False
     50        self.dirtol = [1.e-5, 1.e-5] # direction tolerance in rad (2 arcsec)
     51
     52        self.tables = []
     53        self.nshift = -1
     54        self.nchan = -1
     55
     56        self.set_data(infiles)
     57       
     58        #self.separator = sidebandsep()
     59
     60    @asaplog_post_dec
     61    def set_data(self, infiles):
     62        """
     63        Set data to be processed.
     64
     65        infiles  : a list of filenames or scantables
     66        """
     67        if not (type(infiles) in (list, tuple, numpy.ndarray)):
     68            infiles = [infiles]
     69        if isinstance(infiles[0], scantable):
     70            # a list of scantable
     71            for stab in infiles:
     72                if not isinstance(stab, scantable):
     73                    asaplog.post()
     74                    raise TypeError, "Input data is not a list of scantables."
     75
     76            #self.separator._setdata(infiles)
     77            self._reset_data()
     78            self.intables = infiles
     79        else:
     80            # a list of filenames
     81            for name in infiles:
     82                if not os.path.exists(name):
     83                    asaplog.post()
     84                    raise ValueError, "Could not find input file '%s'" % name
     85           
     86            #self.separator._setdataname(infiles)
     87            self._reset_data()
     88            self.intables = infiles
     89
     90        asaplog.push("%d files are set to process" % len(self.intables))
     91
     92
     93    def _reset_data(self):
     94        del self.intables
     95        self.intables = None
     96        self.signalShift = []
     97        #self.imageShift = []
     98        self.tables = []
     99        self.nshift = -1
     100        self.nchan = -1
     101
     102    @asaplog_post_dec
    42103    def set_frequency(self, baseif, freqtol, frame=""):
    43104        """
     
    46107        Parameters:
    47108          - reference IFNO to process in the first table in the list
    48           - frequency tolerance from reference IF to select data (string)
     109          - frequency tolerance from reference IF to select data
    49110          frame  : frequency frame to select IF
    50111        """
    51         if type(freqtol) in (float, int):
    52             freqtol = str(freqtol)
    53         elif isinstance(freqtol, dict):
    54             try:
    55                 freqtol = str(freqtol['value']) + freqtol['unit']
    56             except:
    57                 raise ValueError, "Invalid frequency tolerance."
    58         self._separator.set_freq(baseif, freqtol, frame)
    59 
    60 
    61     def set_dirtol(self, dirtol=["2arcsec", "2arcsec"]):
     112        self._reset_if()
     113        self.baseif = baseif
     114        if isinstance(freqtol,dict) and freqtol["unit"] == "Hz":
     115            if freqtol['value'] > 0.:
     116                self.freqtol = freqtol
     117            else:
     118                asaplog.post()
     119                asaplog.push("Frequency tolerance should be positive value.")
     120                asaplog.post("ERROR")
     121                return
     122        else:
     123            # torelance in channel unit
     124            if freqtol > 0:
     125                self.freqtol = float(freqtol)
     126            else:
     127                asaplog.post()
     128                asaplog.push("Frequency tolerance should be positive value.")
     129                asaplog.post("ERROR")
     130                return
     131        self.freqframe = frame
     132
     133    def _reset_if(self):
     134        self.baseif = -1
     135        self.freqtol = 0
     136        self.freqframe = ""
     137        self.signalShift = []
     138        #self.imageShift = []
     139        self.tables = []
     140        self.nshift = 0
     141        self.nchan = -1
     142
     143    @asaplog_post_dec
     144    def set_dirtol(self, dirtol=[1.e-5,1.e-5]):
    62145        """
    63146        Set tolerance of direction to select data
    64147        """
    65         if isinstance(dirtol, str):
    66             dirtol = [dirtol]
    67 
    68         self._separator.set_dirtol(dirtol)
    69    
    70            
    71     def set_shift(self, imageshift=[]):
     148        # direction tolerance in rad
     149        if not (type(dirtol) in [list, tuple, numpy.ndarray]):
     150            dirtol = [dirtol, dirtol]
     151        if len(dirtol) == 1:
     152            dirtol = [dirtol[0], dirtol[0]]
     153        if len(dirtol) > 1:
     154            self.dirtol = dirtol[0:2]
     155        else:
     156            asaplog.post()
     157            asaplog.push("Invalid direction tolerance. Should be a list of float in unit radian")
     158            asaplog.post("ERROR")
     159            return
     160        asaplog.post("Set direction tolerance [%f, %f] (rad)" % \
     161                     (self.dirtol[0], self.dirtol[1]))
     162
     163    @asaplog_post_dec
     164    def set_shift(self, mode="DSB", imageshift=None):
    72165        """
    73166        Set shift mode and channel shift of image band.
    74167
    75         imageshift : a list of number of channels shifted in image
    76                      side band of each scantable.
    77                      If the shifts are not set, they are assumed to be
    78                      equal to those of signal side band, but in opposite
    79                      direction as usual by LO1 offsetting of DSB receivers.
    80         """
    81         if not imageshift:
    82             imageshift = []
    83         self._separator.set_shift(imageshift)
    84 
     168        mode       : shift mode ['DSB'|'SSB']
     169                     When mode='DSB', imageshift is assumed to be equal
     170                     to the shift of signal sideband but in opposite direction.
     171        imageshift : a list of number of channel shift in image band of
     172                     each scantable. valid only mode='SSB'
     173        """
     174        if mode.upper().startswith("S"):
     175            if not imageshift:
     176                raise ValueError, "Need to set shift value of image sideband"
     177            self.dsbmode = False
     178            self.imageShift = imageshift
     179            asaplog.push("Image sideband shift is set manually: %s" % str(self.imageShift))
     180        else:
     181            # DSB mode
     182            self.dsbmode = True
     183            self.imageShift = []
    85184
    86185    @asaplog_post_dec
     
    89188        Resolve both image and signal sideband when True.
    90189        """
    91         self._separator.solve_both(flag)
    92         if flag:
    93             asaplog.push("Both signal and image sidebands will be solved and stored in separate tables.")
    94         else:
    95             asaplog.push("Only signal sideband will be solved and stored in an table.")
     190        self.getboth = flag
     191        if self.getboth:
     192            asaplog.push("Both signal and image sidebands are solved and output as separate tables.")
     193        else:
     194            asaplog.push("Only signal sideband is solved and output as an table.")
    96195
    97196    @asaplog_post_dec
     
    100199        Set rejection limit of solution.
    101200        """
    102         self._separator.set_limit(threshold)
     201        #self.separator._setlimit(abs(threshold))
     202        self.rejlimit = threshold
     203        asaplog.push("The threshold of rejection is set to %f" % self.rejlimit)
    103204
    104205
     
    109210        when True.
    110211        """
    111         self._separator.subtract_other(flag)
     212        self.solveother = flag
    112213        if flag:
    113214            asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.")
    114215
    115216
    116     def set_lo1(self, lo1, frame="TOPO", reftime=-1, refdir=""):
    117         """
    118         Set LO1 frequency to calculate frequency of image sideband.
    119 
    120         lo1     : LO1 frequency
    121         frame   : the frequency frame of LO1
    122         reftime : the reference time used in frequency frame conversion.
    123         refdir  : the reference direction used in frequency frame conversion.
    124         """
    125         self._separator.set_lo1(lo1, frame, reftime, refdir)
    126 
    127 
    128     def set_lo1root(self, name):
    129         """
    130         Set MS name which stores LO1 frequency of signal side band.
    131         It is used to calculate frequency of image sideband.
    132 
    133         name : MS name which contains 'ASDM_SPECTRALWINDOW' and
    134                'ASDM_RECEIVER' tables.
    135         """
    136         self._separator.set_lo1root(name)
    137 
    138217    @asaplog_post_dec
    139218    def separate(self, outname="", overwrite=False):
     
    144223        overwrite : overwrite existing table
    145224        """
    146         out_default = "sbseparated.asap"
    147         if len(outname) == 0:
    148             outname = out_default
    149             asaplog.post()
    150             asaplog.push("The output file name is not specified.")
    151             asaplog.push("Using default name '%s'" % outname)
     225        # List up valid scantables and IFNOs to convolve.
     226        #self.separator._separate()
     227        self._setup_shift()
     228        #self._preprocess_tables()
     229
     230        nshift = len(self.tables)
     231        signaltab = self._grid_outtable(self.tables[0].copy())
     232        if self.getboth:
     233            imagetab = signaltab.copy()
     234
     235        rejrow = []
     236        for irow in xrange(signaltab.nrow()):
     237            currpol = signaltab.getpol(irow)
     238            currbeam = signaltab.getbeam(irow)
     239            currdir = signaltab.get_directionval(irow)
     240            spec_array, tabidx = self._get_specarray(polid=currpol,\
     241                                                     beamid=currbeam,\
     242                                                     dir=currdir)
     243            #if not spec_array:
     244            if len(tabidx) == 0:
     245                asaplog.post()
     246                asaplog.push("skipping row %d" % irow)
     247                rejrow.append(irow)
     248                continue
     249            signal = self._solve_signal(spec_array, tabidx)
     250            signaltab.set_spectrum(signal, irow)
     251            if self.getboth:
     252                image = self._solve_image(spec_array, tabidx)
     253                imagetab.set_spectrum(image, irow)
     254       
     255        # TODO: Need to remove rejrow form scantables here
     256        signaltab.flag_row(rejrow)
     257        if self.getboth:
     258            imagetab.flag_row(rejrow)
     259       
     260        if outname == "":
     261            outname = "sbsepareted.asap"
     262        signalname = outname + ".signalband"
     263        if os.path.exists(signalname):
     264            if not overwrite:
     265                raise Exception, "Output file '%s' exists." % signalname
     266            else:
     267                shutil.rmtree(signalname)
     268        signaltab.save(signalname)
     269        if self.getboth:
     270            # Warnings
     271            asaplog.post()
     272            asaplog.push("Saving IMAGE sideband.")
     273            asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")
     274            asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")
    152275            asaplog.post("WARN")
    153276
    154         if os.path.exists(outname):
    155             if overwrite:
    156                 asaplog.push("removing the old file '%s'" % outname)
    157                 shutil.rmtree(outname)
     277            imagename = outname + ".imageband"
     278            if os.path.exists(imagename):
     279                if not overwrite:
     280                    raise Exception, "Output file '%s' exists." % imagename
     281                else:
     282                    shutil.rmtree(imagename)
     283            imagetab.save(imagename)
     284
     285
     286    def _solve_signal(self, data, tabidx=None):
     287        if not tabidx:
     288            tabidx = range(len(data))
     289
     290        tempshift = []
     291        dshift = []
     292        if self.solveother:
     293            for idx in tabidx:
     294                tempshift.append(-self.imageShift[idx])
     295                dshift.append(self.signalShift[idx] - self.imageShift[idx])
     296        else:
     297            for idx in tabidx:
     298                tempshift.append(-self.signalShift[idx])
     299                dshift.append(self.imageShift[idx] - self.signalShift[idx])
     300
     301        shiftdata = numpy.zeros(data.shape, numpy.float)
     302        for i in range(len(data)):
     303            shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
     304        ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
     305        result_image = self._combineResult(ifftdata)
     306        if not self.solveother:
     307            return result_image
     308        result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
     309        return result_signal
     310
     311
     312    def _solve_image(self, data, tabidx=None):
     313        if not tabidx:
     314            tabidx = range(len(data))
     315
     316        tempshift = []
     317        dshift = []
     318        if self.solveother:
     319            for idx in tabidx:
     320                tempshift.append(-self.signalShift[idx])
     321                dshift.append(self.imageShift[idx] - self.signalShift[idx])
     322        else:
     323            for idx in tabidx:
     324                tempshift.append(-self.imageShift[idx])
     325                dshift.append(self.signalShift[idx] - self.imageShift[idx])
     326
     327        shiftdata = numpy.zeros(data.shape, numpy.float)
     328        for i in range(len(data)):
     329            shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
     330        ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
     331        result_image = self._combineResult(ifftdata)
     332        if not self.solveother:
     333            return result_image
     334        result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
     335        return result_signal
     336
     337    @asaplog_post_dec
     338    def _grid_outtable(self, table):
     339        # Generate gridded table for output (Just to get rows)
     340        gridder = asapgrid2(table)
     341        gridder.setIF(self.baseif)
     342       
     343        cellx = str(self.dirtol[0])+"rad"
     344        celly = str(self.dirtol[1])+"rad"
     345        dirarr = numpy.array(table.get_directionval()).transpose()
     346        mapx = dirarr[0].max() - dirarr[0].min()
     347        mapy = dirarr[1].max() - dirarr[1].min()
     348        nx = max(1, numpy.ceil(mapx/self.dirtol[0]))
     349        ny = max(1, numpy.ceil(mapy/self.dirtol[0]))
     350
     351        asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
     352                     (cellx, celly))
     353        gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
     354        gridder.setFunc(func='box', convsupport=1)
     355        gridder.setWeight(weightType='uniform')
     356        gridder.grid()
     357        return gridder.getResult()
     358
     359    @asaplog_post_dec
     360    def _get_specarray(self, polid=None, beamid=None, dir=None):
     361        ntable = len(self.tables)
     362        spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
     363        nspec = 0
     364        asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
     365        tabidx = []
     366        for itab in range(ntable):
     367            tab = self.tables[itab]
     368            # Select rows by POLNO and BEAMNO
     369            try:
     370                tab.set_selection(pols=[polid], beams=[beamid])
     371                if tab.nrow() > 0: tabidx.append(itab)
     372            except: # no selection
     373                asaplog.post()
     374                asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
     375                asaplog.post("WARN")
     376                continue
     377
     378            # Select rows by direction
     379            spec = numpy.zeros(self.nchan, numpy.float)
     380            selrow = []
     381            for irow in range(tab.nrow()):
     382                currdir = tab.get_directionval(irow)
     383                if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
     384                   (abs(currdir[1] - dir[1]) > self.dirtol[1]):
     385                    continue
     386                selrow.append(irow)
     387            if len(selrow) == 0:
     388                asaplog.post()
     389                asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
     390                asaplog.post("WARN")
     391                continue
    158392            else:
    159                 asaplog.post()
    160                 asaplog.push("Output file '%s' already exists." % outname)
     393                seltab = tab.copy()
     394                seltab.set_selection(selector(rows=selrow))
     395           
     396            if tab.nrow() > 1:
     397                asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
     398                tab = seltab.average_time(scanav=False, weight="tintsys")
     399            else:
     400                tab = seltab
     401
     402            spec_array[nspec] = tab._getspectrum()
     403            nspec += 1
     404
     405        if nspec != ntable:
     406            asaplog.post()
     407            #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
     408            asaplog.push("Could not find corresponding rows in some tables.")
     409            asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
     410            if nspec < 2:
     411                asaplog.push("At least 2 spectra are necessary for convolution")
    161412                asaplog.post("ERROR")
    162                 return False
    163 
    164         self._separator.separate(outname)
    165 
     413                return False, tabidx
     414
     415        return spec_array[0:nspec], tabidx
     416           
     417
     418    @asaplog_post_dec
     419    def _setup_shift(self):
     420        ### define self.tables, self.signalShift, and self.imageShift
     421        if not self.intables:
     422            asaplog.post()
     423            raise RuntimeError, "Input data is not defined."
     424        #if self.baseif < 0:
     425        #    asaplog.post()
     426        #    raise RuntimeError, "Reference IFNO is not defined."
     427       
     428        byname = False
     429        #if not self.intables:
     430        if isinstance(self.intables[0], str):
     431            # A list of file name is given
     432            if not os.path.exists(self.intables[0]):
     433                asaplog.post()
     434                raise RuntimeError, "Could not find '%s'" % self.intables[0]
     435           
     436            stab = scantable(self.intables[0],average=False)
     437            ntab = len(self.intables)
     438            byname = True
     439        else:
     440            stab = self.intables[0]
     441            ntab = len(self.intables)
     442
     443        if len(stab.getbeamnos()) > 1:
     444            asaplog.post()
     445            asaplog.push("Mult-beam data is not supported by this module.")
     446            asaplog.post("ERROR")
     447            return
     448
     449        valid_ifs = stab.getifnos()
     450        if self.baseif < 0:
     451            self.baseif = valid_ifs[0]
     452            asaplog.post()
     453            asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
     454       
     455        if not (self.baseif in valid_ifs):
     456            asaplog.post()
     457            errmsg = "IF%d does not exist in the first scantable" %  \
     458                     self.baseif
     459            raise RuntimeError, errmsg
     460
     461        asaplog.push("Start selecting tables and IFNOs to solve.")
     462        asaplog.push("Cheching frequency of the reference IF")
     463        unit_org = stab.get_unit()
     464        coord = stab._getcoordinfo()
     465        frame_org = coord[1]
     466        stab.set_unit("Hz")
     467        if len(self.freqframe) > 0:
     468            stab.set_freqframe(self.freqframe)
     469        stab.set_selection(ifs=[self.baseif])
     470        spx = stab._getabcissa()
     471        stab.set_selection()
     472        basech0 = spx[0]
     473        baseinc = spx[1]-spx[0]
     474        self.nchan = len(spx)
     475        if isinstance(self.freqtol, float):
     476            vftol = abs(baseinc * self.freqtol)
     477            self.freqtol = dict(value=vftol, unit="Hz")
     478        else:
     479            vftol = abs(self.freqtol['value'])
     480        inctol = abs(baseinc/float(self.nchan))
     481        asaplog.push("Reference frequency setup (Table = 0, IFNO = %d):  nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
     482        asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
     483        poltype0 = stab.poltype()
     484       
     485        self.tables = []
     486        self.signalShift = []
     487        if self.dsbmode:
     488            self.imageShift = []
     489
     490        for itab in range(ntab):
     491            asaplog.push("Table %d:" % itab)
     492            tab_selected = False
     493            if itab > 0:
     494                if byname:
     495                    stab = scantable(self.intables[itab],average=False)
     496                else:
     497                    stab = self.intables[itab]
     498                unit_org = stab.get_unit()
     499                coord = stab._getcoordinfo()
     500                frame_org = coord[1]
     501                stab.set_unit("Hz")
     502                if len(self.freqframe) > 0:
     503                    stab.set_freqframe(self.freqframe)
     504
     505            # Check POLTYPE should be equal to the first table.
     506            if stab.poltype() != poltype0:
     507                asaplog.post()
     508                raise Exception, "POLTYPE should be equal to the first table."
     509            # Multiple beam data may not handled properly
     510            if len(stab.getbeamnos()) > 1:
     511                asaplog.post()
     512                asaplog.push("table contains multiple beams. It may not be handled properly.")
     513                asaplog.push("WARN")
     514
     515            for ifno in stab.getifnos():
     516                stab.set_selection(ifs=[ifno])
     517                spx = stab._getabcissa()
     518                if (len(spx) != self.nchan) or \
     519                   (abs(spx[0]-basech0) > vftol) or \
     520                   (abs(spx[1]-spx[0]-baseinc) > inctol):
     521                    continue
     522                tab_selected = True
     523                seltab = stab.copy()
     524                seltab.set_unit(unit_org)
     525                seltab.set_freqframe(frame_org)
     526                self.tables.append(seltab)
     527                self.signalShift.append((spx[0]-basech0)/baseinc)
     528                if self.dsbmode:
     529                    self.imageShift.append(-self.signalShift[-1])
     530                asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
     531            stab.set_selection()
     532            stab.set_unit(unit_org)
     533            stab.set_freqframe(frame_org)
     534            if not tab_selected:
     535                asaplog.post()
     536                asaplog.push("No data selected in Table %d" % itab)
     537                asaplog.post("WARN")
     538
     539        asaplog.push("Total number of IFs selected = %d" % len(self.tables))
     540        if len(self.tables) < 2:
     541            asaplog.post()
     542            raise RuntimeError, "At least 2 IFs are necessary for convolution!"
     543
     544        if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
     545            asaplog.post()
     546            errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
     547            raise RuntimeError, errmsg
     548
     549        self.signalShift = numpy.array(self.signalShift)
     550        self.imageShift = numpy.array(self.imageShift)
     551        self.nshift = len(self.tables)
     552
     553    @asaplog_post_dec
     554    def _preprocess_tables(self):
     555        ### temporary method to preprocess data
     556        ### Do time averaging for now.
     557        for itab in range(len(self.tables)):
     558            self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
     559       
     560
     561#     def save(self, outfile, outform="ASAP", overwrite=False):
     562#         if not overwrite and os.path.exists(outfile):
     563#             raise RuntimeError, "Output file '%s' already exists" % outfile
     564#
     565#         #self.separator._save(outfile, outform)
     566
     567#     def done(self):
     568#         self.close()
     569
     570#     def close(self):
     571#         pass
     572#         #del self.separator
     573   
     574
     575
     576########################################################################
     577    def _Deconvolution(self, data_array, shift, threshold=0.00000001):
     578        FObs = []
     579        Reject = 0
     580        nshift, nchan = data_array.shape
     581        nspec = nshift*(nshift-1)/2
     582        ifftObs  = numpy.zeros((nspec, nchan), numpy.float)
     583        for i in range(nshift):
     584           F = FFT.fft(data_array[i])
     585           FObs.append(F)
     586        z = 0
     587        for i in range(nshift):
     588            for j in range(i+1, nshift):
     589                Fobs = (FObs[i]+FObs[j])/2.0
     590                dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
     591                #print 'dX,i,j=',dX,i,j
     592                for k in range(1,self.nchan):
     593                    if math.fabs(math.sin(dX*k)) > threshold:
     594                        Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
     595                    else: Reject += 1
     596                ifftObs[z] = FFT.ifft(Fobs)
     597                z += 1
     598        print 'Threshold=%s Reject=%d' % (threshold, Reject)
     599        return ifftObs
     600
     601    def _combineResult(self, ifftObs):
     602        nspec = len(ifftObs)
     603        sum = ifftObs[0]
     604        for i in range(1,nspec):
     605            sum += ifftObs[i]
     606        return(sum/float(nspec))
     607
     608    def _subtractOtherSide(self, data_array, shift, Data):
     609        sum = numpy.zeros(len(Data), numpy.float)
     610        numSP = len(data_array)
     611        for i in range(numSP):
     612            SPsub = data_array[i] - Data
     613            sum += self._shiftSpectrum(SPsub, -shift[i])
     614        return(sum/float(numSP))
     615
     616    def _shiftSpectrum(self, data, Shift):
     617        Out = numpy.zeros(self.nchan, numpy.float)
     618        w2 = Shift % 1
     619        w1 = 1.0 - w2
     620        for i in range(self.nchan):
     621            c1 = int((Shift + i) % self.nchan)
     622            c2 = (c1 + 1) % self.nchan
     623            Out[c1] += data[i] * w1
     624            Out[c2] += data[i] * w2
     625        return Out.copy()
Note: See TracChangeset for help on using the changeset viewer.