Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/sbseparator.py

    r2685 r2859  
    99from asap.selector import selector
    1010from asap.asapgrid import asapgrid2
    11 #from asap._asap import sidebandsep
     11from asap._asap import SBSeparator
    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.
    2423
    2524    Example:
    … …  
    3837    """
    3938    def __init__(self, infiles):
    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))
     39        self._separator = SBSeparator(infiles)
    9140
    9241
    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
    10342    def set_frequency(self, baseif, freqtol, frame=""):
    10443        """
    … …  
    10746        Parameters:
    10847          - reference IFNO to process in the first table in the list
    109           - frequency tolerance from reference IF to select data
     48          - frequency tolerance from reference IF to select data (string)
    11049          frame  : frequency frame to select IF
    11150        """
    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
     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)
    13259
    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
    14260
    143     @asaplog_post_dec
    144     def set_dirtol(self, dirtol=[1.e-5,1.e-5]):
     61    def set_dirtol(self, dirtol=["2arcsec", "2arcsec"]):
    14562        """
    14663        Set tolerance of direction to select data
    14764        """
    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]))
     65        if isinstance(dirtol, str):
     66            dirtol = [dirtol]
    16267
    163     @asaplog_post_dec
    164     def set_shift(self, mode="DSB", imageshift=None):
     68        self._separator.set_dirtol(dirtol)
     69   
     70           
     71    def set_shift(self, imageshift=[]):
    16572        """
    16673        Set shift mode and channel shift of image band.
    16774
    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'
     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.
    17380        """
    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 = []
     81        if not imageshift:
     82            imageshift = []
     83        self._separator.set_shift(imageshift)
     84
    18485
    18586    @asaplog_post_dec
    … …  
    18889        Resolve both image and signal sideband when True.
    18990        """
    190         self.getboth = flag
    191         if self.getboth:
    192             asaplog.push("Both signal and image sidebands are solved and output as separate tables.")
     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.")
    19394        else:
    194             asaplog.push("Only signal sideband is solved and output as an table.")
     95            asaplog.push("Only signal sideband will be solved and stored in an table.")
    19596
    19697    @asaplog_post_dec
    … …  
    199100        Set rejection limit of solution.
    200101        """
    201         #self.separator._setlimit(abs(threshold))
    202         self.rejlimit = threshold
    203         asaplog.push("The threshold of rejection is set to %f" % self.rejlimit)
     102        self._separator.set_limit(threshold)
    204103
    205104
    … …  
    210109        when True.
    211110        """
    212         self.solveother = flag
     111        self._separator.subtract_other(flag)
    213112        if flag:
    214113            asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.")
    215114
     115
     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)
    216137
    217138    @asaplog_post_dec
    … …  
    223144        overwrite : overwrite existing table
    224145        """
    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
     146        out_default = "sbseparated.asap"
     147        if len(outname) == 0:
     148            outname = out_default
    271149            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.")
     150            asaplog.push("The output file name is not specified.")
     151            asaplog.push("Using default name '%s'" % outname)
    275152            asaplog.post("WARN")
    276153
    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)
     154        if os.path.exists(outname):
     155            if overwrite:
     156                asaplog.push("removing the old file '%s'" % outname)
     157                shutil.rmtree(outname)
     158            else:
     159                asaplog.post()
     160                asaplog.push("Output file '%s' already exists." % outname)
     161                asaplog.post("ERROR")
     162                return False
    284163
     164        self._separator.separate(outname)
    285165
    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
    392             else:
    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")
    412                 asaplog.post("ERROR")
    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.