source: trunk/python/sbseparator.py @ 2707

Last change on this file since 2707 was 2707, checked in by Kana Sugimoto, 11 years ago

New Development: Yes

JIRA Issue: Yes (CAS-4141)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

  • a new c++ class, STSideBandSep and its python interface.
  • a new method set_lo1 in the sbseparator module.

Test Programs: manual tests by CSV scientists.

Put in Release Notes: No

Module(s): STSidebandSep (a new class) and sbseparator

Description:

Added a new c++ class, STSideBandSep, and its python interface (python_STSideBandSep.cpp).
The class works on calculating frequency information of image side band and fill scantable.
Also added a new python method, set_lo1(double), in sbseparator module. This is to allow
user to set LO1 frequency to solve frequencies of image side band manually.
File size: 24.2 KB
Line 
1import os, shutil
2import numpy
3import numpy.fft as FFT
4import math
5
6from asap.scantable import scantable
7from asap.parameters import rcParams
8from asap.logging import asaplog, asaplog_post_dec
9from asap.selector import selector
10from asap.asapgrid import asapgrid2
11from asap._asap import SBSeparator
12
13class sbseparator:
14    """
15    The sbseparator class is defined to separate SIGNAL and IMAGE
16    sideband spectra observed by frequency-switching technique.
17    It also helps supressing emmission of IMAGE sideband.
18    *** WARNING *** THIS MODULE IS EXPERIMENTAL
19    Known issues:
20    - Frequencies of IMAGE sideband cannot be reconstructed from
21      information in scantable in sideband sparation. Frequency
22      setting of SIGNAL sideband is stored in output table for now.
23    - Flag information (stored in FLAGTRA) is ignored.
24
25    Example:
26        # Create sideband separator instance whith 3 input data
27        sbsep = sbseparator(['test1.asap', 'test2.asap', 'test3.asap'])
28        # Set reference IFNO and tolerance to select data
29        sbsep.set_frequency(5, 30, frame='TOPO')
30        # Set direction tolerance to select data in unit of radian
31        sbsep.set_dirtol(1.e-5)
32        # Set rejection limit of solution
33        sbsep.set_limit(0.2)
34        # Solve image sideband as well
35        sbsep.set_both(True)
36        # Invoke sideband separation
37        sbsep.separate('testout.asap', overwrite = True)
38    """
39    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        self.lo1 = 0.
52
53        self.tables = []
54        self.nshift = -1
55        self.nchan = -1
56
57        self.set_data(infiles)
58       
59        #self.separator = sidebandsep()
60
61    @asaplog_post_dec
62    def set_data(self, infiles):
63        """
64        Set data to be processed.
65
66        infiles  : a list of filenames or scantables
67        """
68        if not (type(infiles) in (list, tuple, numpy.ndarray)):
69            infiles = [infiles]
70        if isinstance(infiles[0], scantable):
71            # a list of scantable
72            for stab in infiles:
73                if not isinstance(stab, scantable):
74                    asaplog.post()
75                    raise TypeError, "Input data is not a list of scantables."
76
77            #self.separator._setdata(infiles)
78            self._reset_data()
79            self.intables = infiles
80        else:
81            # a list of filenames
82            for name in infiles:
83                if not os.path.exists(name):
84                    asaplog.post()
85                    raise ValueError, "Could not find input file '%s'" % name
86           
87            #self.separator._setdataname(infiles)
88            self._reset_data()
89            self.intables = infiles
90
91        asaplog.push("%d files are set to process" % len(self.intables))
92
93
94    def _reset_data(self):
95        del self.intables
96        self.intables = None
97        self.signalShift = []
98        #self.imageShift = []
99        self.tables = []
100        self.nshift = -1
101        self.nchan = -1
102
103    @asaplog_post_dec
104    def set_frequency(self, baseif, freqtol, frame=""):
105        """
106        Set IFNO and frequency tolerance to select data to process.
107
108        Parameters:
109          - reference IFNO to process in the first table in the list
110          - frequency tolerance from reference IF to select data
111          frame  : frequency frame to select IF
112        """
113        self._reset_if()
114        self.baseif = baseif
115        if isinstance(freqtol,dict) and freqtol["unit"] == "Hz":
116            if freqtol['value'] > 0.:
117                self.freqtol = freqtol
118            else:
119                asaplog.post()
120                asaplog.push("Frequency tolerance should be positive value.")
121                asaplog.post("ERROR")
122                return
123        else:
124            # torelance in channel unit
125            if freqtol > 0:
126                self.freqtol = float(freqtol)
127            else:
128                asaplog.post()
129                asaplog.push("Frequency tolerance should be positive value.")
130                asaplog.post("ERROR")
131                return
132        self.freqframe = frame
133
134    def _reset_if(self):
135        self.baseif = -1
136        self.freqtol = 0
137        self.freqframe = ""
138        self.signalShift = []
139        #self.imageShift = []
140        self.tables = []
141        self.nshift = 0
142        self.nchan = -1
143
144    @asaplog_post_dec
145    def set_dirtol(self, dirtol=[1.e-5,1.e-5]):
146        """
147        Set tolerance of direction to select data
148        """
149        # direction tolerance in rad
150        if not (type(dirtol) in [list, tuple, numpy.ndarray]):
151            dirtol = [dirtol, dirtol]
152        if len(dirtol) == 1:
153            dirtol = [dirtol[0], dirtol[0]]
154        if len(dirtol) > 1:
155            self.dirtol = dirtol[0:2]
156        else:
157            asaplog.post()
158            asaplog.push("Invalid direction tolerance. Should be a list of float in unit radian")
159            asaplog.post("ERROR")
160            return
161        asaplog.post("Set direction tolerance [%f, %f] (rad)" % \
162                     (self.dirtol[0], self.dirtol[1]))
163
164    @asaplog_post_dec
165    def set_shift(self, mode="DSB", imageshift=None):
166        """
167        Set shift mode and channel shift of image band.
168
169        mode       : shift mode ['DSB'|'SSB'(='2SB')]
170                     When mode='DSB', imageshift is assumed to be equal
171                     to the shift of signal sideband but in opposite direction.
172        imageshift : a list of number of channel shift in image band of
173                     each scantable. valid only mode='SSB'
174        """
175        if mode.upper().startswith("D"):
176            # DSB mode
177            self.dsbmode = True
178            self.imageShift = []
179        else:
180            if not imageshift:
181                raise ValueError, "Need to set shift value of image sideband"
182            self.dsbmode = False
183            self.imageShift = imageshift
184            asaplog.push("Image sideband shift is set manually: %s" % str(self.imageShift))
185
186    @asaplog_post_dec
187    def set_both(self, flag=False):
188        """
189        Resolve both image and signal sideband when True.
190        """
191        self.getboth = flag
192        if self.getboth:
193            asaplog.push("Both signal and image sidebands are solved and output as separate tables.")
194        else:
195            asaplog.push("Only signal sideband is solved and output as an table.")
196
197    @asaplog_post_dec
198    def set_limit(self, threshold=0.2):
199        """
200        Set rejection limit of solution.
201        """
202        #self.separator._setlimit(abs(threshold))
203        self.rejlimit = threshold
204        asaplog.push("The threshold of rejection is set to %f" % self.rejlimit)
205
206
207    @asaplog_post_dec
208    def set_solve_other(self, flag=False):
209        """
210        Calculate spectra by subtracting the solution of the other sideband
211        when True.
212        """
213        self.solveother = flag
214        if flag:
215            asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.")
216
217    def set_lo1(self,lo1):
218        """
219        Set LO1 frequency to calculate frequency of image sideband.
220
221        lo1 : LO1 frequency in float
222        """
223        if isinstance(lo1, dict) and lo1["unit"] == "Hz":
224            self.lo1 = lo1["value"]
225        else:
226            self.lo1 = float(lo1)
227        if self.lo1 <= 0.:
228            asaplog.push("Got negative LO1 frequency. It will be ignored.")
229            asaplog.post("WARN")
230
231
232    @asaplog_post_dec
233    def separate(self, outname="", overwrite=False):
234        """
235        Invoke sideband separation.
236
237        outname   : a name of output scantable
238        overwrite : overwrite existing table
239        """
240        # List up valid scantables and IFNOs to convolve.
241        #self.separator._separate()
242        self._setup_shift()
243        #self._preprocess_tables()
244
245        nshift = len(self.tables)
246        signaltab = self._grid_outtable(self.tables[0])
247        if self.getboth:
248            imagetab = signaltab.copy()
249
250        rejrow = []
251        for irow in xrange(signaltab.nrow()):
252            currpol = signaltab.getpol(irow)
253            currbeam = signaltab.getbeam(irow)
254            currdir = signaltab.get_directionval(irow)
255            spec_array, tabidx = self._get_specarray(polid=currpol,\
256                                                     beamid=currbeam,\
257                                                     dir=currdir)
258            #if not spec_array:
259            if len(tabidx) == 0:
260                asaplog.post()
261                asaplog.push("skipping row %d" % irow)
262                rejrow.append(irow)
263                continue
264            signal = self._solve_signal(spec_array, tabidx)
265            signaltab.set_spectrum(signal, irow)
266
267            # Solve image side side band
268            if self.getboth:
269                image = self._solve_image(spec_array, tabidx)
270                imagetab.set_spectrum(image, irow)
271
272        # TODO: Need to remove rejrow form scantables here
273        signaltab.flag_row(rejrow)
274        if self.getboth:
275            imagetab.flag_row(rejrow)
276       
277        if outname == "":
278            outname = "sbsepareted.asap"
279        signalname = outname + ".signalband"
280        if os.path.exists(signalname):
281            if not overwrite:
282                raise Exception, "Output file '%s' exists." % signalname
283            else:
284                shutil.rmtree(signalname)
285        signaltab.save(signalname)
286
287        if self.getboth:
288            # Warnings
289            asaplog.post()
290            asaplog.push("Saving IMAGE sideband.")
291            asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")
292            asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")
293            asaplog.post("WARN")
294
295            imagename = outname + ".imageband"
296            if os.path.exists(imagename):
297                if not overwrite:
298                    raise Exception, "Output file '%s' exists." % imagename
299                else:
300                    shutil.rmtree(imagename)
301            # Update frequency information
302            sbsep = SBSeparator()
303            sbsep.set_imgtable(imagetab)
304            if self.lo1 > 0.:
305                sbsep.set_lo1(self.lo1)
306            sbsep.solve_imgfreq()
307            imagetab.save(imagename)
308            del sbsep
309
310
311    def _solve_signal(self, data, tabidx=None):
312        if not tabidx:
313            tabidx = range(len(data))
314
315        tempshift = []
316        dshift = []
317        if self.solveother:
318            for idx in tabidx:
319                tempshift.append(-self.imageShift[idx])
320                dshift.append(self.signalShift[idx] - self.imageShift[idx])
321        else:
322            for idx in tabidx:
323                tempshift.append(-self.signalShift[idx])
324                dshift.append(self.imageShift[idx] - self.signalShift[idx])
325
326        shiftdata = numpy.zeros(data.shape, numpy.float)
327        for i in range(len(data)):
328            shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
329        ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
330        result_image = self._combineResult(ifftdata)
331        if not self.solveother:
332            return result_image
333        result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
334        return result_signal
335
336
337    def _solve_image(self, data, tabidx=None):
338        if not tabidx:
339            tabidx = range(len(data))
340
341        tempshift = []
342        dshift = []
343        if self.solveother:
344            for idx in tabidx:
345                tempshift.append(-self.signalShift[idx])
346                dshift.append(self.imageShift[idx] - self.signalShift[idx])
347        else:
348            for idx in tabidx:
349                tempshift.append(-self.imageShift[idx])
350                dshift.append(self.signalShift[idx] - self.imageShift[idx])
351
352        shiftdata = numpy.zeros(data.shape, numpy.float)
353        for i in range(len(data)):
354            shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
355        ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
356        result_image = self._combineResult(ifftdata)
357        if not self.solveother:
358            return result_image
359        result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
360        return result_signal
361
362    @asaplog_post_dec
363    def _grid_outtable(self, table):
364        # Generate gridded table for output (Just to get rows)
365        gridder = asapgrid2(table)
366        gridder.setIF(self.baseif)
367
368        cellx = str(self.dirtol[0])+"rad"
369        celly = str(self.dirtol[1])+"rad"
370        dirarr = numpy.array(table.get_directionval()).transpose()
371        mapx = dirarr[0].max() - dirarr[0].min()
372        mapy = dirarr[1].max() - dirarr[1].min()
373        centy = 0.5 * (dirarr[1].max() + dirarr[1].min())
374        nx = max(1, numpy.ceil(mapx*numpy.cos(centy)/self.dirtol[0]))
375        ny = max(1, numpy.ceil(mapy/self.dirtol[0]))
376
377        asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
378                     (cellx, celly))
379        gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
380        gridder.setFunc(func='box', convsupport=1)
381        gridder.setWeight(weightType='uniform')
382        gridder.grid()
383        return gridder.getResult()
384
385    @asaplog_post_dec
386    def _get_specarray(self, polid=None, beamid=None, dir=None):
387        ntable = len(self.tables)
388        spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
389        nspec = 0
390        asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
391        tabidx = []
392        for itab in range(ntable):
393            tab = self.tables[itab]
394            # Select rows by POLNO and BEAMNO
395            try:
396                tab.set_selection(pols=[polid], beams=[beamid])
397                if tab.nrow() > 0: tabidx.append(itab)
398            except: # no selection
399                asaplog.post()
400                asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
401                asaplog.post("WARN")
402                continue
403
404            # Select rows by direction
405            spec = numpy.zeros(self.nchan, numpy.float)
406            selrow = []
407            for irow in range(tab.nrow()):
408                currdir = tab.get_directionval(irow)
409                if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
410                   (abs(currdir[1] - dir[1]) > self.dirtol[1]):
411                    continue
412                selrow.append(irow)
413            if len(selrow) == 0:
414                asaplog.post()
415                asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
416                asaplog.post("WARN")
417                continue
418            else:
419                seltab = tab.copy()
420                seltab.set_selection(selector(rows=selrow))
421
422            if tab.nrow() > 1:
423                asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
424                tab = seltab.average_time(scanav=False, weight="tintsys")
425            else:
426                tab = seltab
427
428            spec_array[nspec] = tab._getspectrum()
429            nspec += 1
430
431        if nspec != ntable:
432            asaplog.post()
433            #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
434            asaplog.push("Could not find corresponding rows in some tables.")
435            asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
436            if nspec < 2:
437                asaplog.push("At least 2 spectra are necessary for convolution")
438                asaplog.post("ERROR")
439                return False, tabidx
440
441        return spec_array[0:nspec], tabidx
442
443
444    @asaplog_post_dec
445    def _setup_shift(self):
446        ### define self.tables, self.signalShift, and self.imageShift
447        if not self.intables:
448            asaplog.post()
449            raise RuntimeError, "Input data is not defined."
450        #if self.baseif < 0:
451        #    asaplog.post()
452        #    raise RuntimeError, "Reference IFNO is not defined."
453       
454        byname = False
455        #if not self.intables:
456        if isinstance(self.intables[0], str):
457            # A list of file name is given
458            if not os.path.exists(self.intables[0]):
459                asaplog.post()
460                raise RuntimeError, "Could not find '%s'" % self.intables[0]
461           
462            stab = scantable(self.intables[0],average=False)
463            ntab = len(self.intables)
464            byname = True
465        else:
466            stab = self.intables[0]
467            ntab = len(self.intables)
468
469        if len(stab.getbeamnos()) > 1:
470            asaplog.post()
471            asaplog.push("Mult-beam data is not supported by this module.")
472            asaplog.post("ERROR")
473            return
474
475        valid_ifs = stab.getifnos()
476        if self.baseif < 0:
477            self.baseif = valid_ifs[0]
478            asaplog.post()
479            asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
480
481        if not (self.baseif in valid_ifs):
482            asaplog.post()
483            errmsg = "IF%d does not exist in the first scantable" %  \
484                     self.baseif
485            raise RuntimeError, errmsg
486
487        asaplog.push("Start selecting tables and IFNOs to solve.")
488        asaplog.push("Checking frequency of the reference IF")
489        unit_org = stab.get_unit()
490        coord = stab._getcoordinfo()
491        frame_org = coord[1]
492        stab.set_unit("Hz")
493        if len(self.freqframe) > 0:
494            stab.set_freqframe(self.freqframe)
495        stab.set_selection(ifs=[self.baseif])
496        spx = stab._getabcissa()
497        stab.set_selection()
498        basech0 = spx[0]
499        baseinc = spx[1]-spx[0]
500        self.nchan = len(spx)
501        # frequency tolerance
502        if isinstance(self.freqtol, dict) and self.freqtol['unit'] == "Hz":
503            vftol = abs(self.freqtol['value'])
504        else:
505            vftol = abs(baseinc * float(self.freqtol))
506            self.freqtol = dict(value=vftol, unit="Hz")
507        # tolerance of frequency increment
508        inctol = abs(baseinc/float(self.nchan))
509        asaplog.push("Reference frequency setup (Table = 0, IFNO = %d):  nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
510        asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
511        poltype0 = stab.poltype()
512
513        self.tables = []
514        self.signalShift = []
515        if self.dsbmode:
516            self.imageShift = []
517
518        for itab in range(ntab):
519            asaplog.push("Table %d:" % itab)
520            tab_selected = False
521            if itab > 0:
522                if byname:
523                    stab = scantable(self.intables[itab],average=False)
524                else:
525                    stab = self.intables[itab]
526                unit_org = stab.get_unit()
527                coord = stab._getcoordinfo()
528                frame_org = coord[1]
529                stab.set_unit("Hz")
530                if len(self.freqframe) > 0:
531                    stab.set_freqframe(self.freqframe)
532
533            # Check POLTYPE should be equal to the first table.
534            if stab.poltype() != poltype0:
535                asaplog.post()
536                raise Exception, "POLTYPE should be equal to the first table."
537            # Multiple beam data may not handled properly
538            if len(stab.getbeamnos()) > 1:
539                asaplog.post()
540                asaplog.push("table contains multiple beams. It may not be handled properly.")
541                asaplog.push("WARN")
542
543            for ifno in stab.getifnos():
544                stab.set_selection(ifs=[ifno])
545                spx = stab._getabcissa()
546                if (len(spx) != self.nchan) or \
547                   (abs(spx[0]-basech0) > vftol) or \
548                   (abs(spx[1]-spx[0]-baseinc) > inctol):
549                    continue
550                tab_selected = True
551                seltab = stab.copy()
552                seltab.set_unit(unit_org)
553                seltab.set_freqframe(frame_org)
554                self.tables.append(seltab)
555                self.signalShift.append((spx[0]-basech0)/baseinc)
556                if self.dsbmode:
557                    self.imageShift.append(-self.signalShift[-1])
558                asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
559            stab.set_selection()
560            stab.set_unit(unit_org)
561            stab.set_freqframe(frame_org)
562            if not tab_selected:
563                asaplog.post()
564                asaplog.push("No data selected in Table %d" % itab)
565                asaplog.post("WARN")
566
567        asaplog.push("Total number of IFs selected = %d" % len(self.tables))
568        if len(self.tables) < 2:
569            asaplog.post()
570            raise RuntimeError, "At least 2 IFs are necessary for convolution!"
571
572        if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
573            asaplog.post()
574            errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
575            errmsg += "\nThe frequency tolerance (freqtol) you set may be too small."
576            raise RuntimeError, errmsg
577
578        self.signalShift = numpy.array(self.signalShift)
579        self.imageShift = numpy.array(self.imageShift)
580        self.nshift = len(self.tables)
581
582    @asaplog_post_dec
583    def _preprocess_tables(self):
584        ### temporary method to preprocess data
585        ### Do time averaging for now.
586        for itab in range(len(self.tables)):
587            self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
588
589#     @asaplog_post_dec
590#     def _setup_image_freq(self, table):
591#         # User defined coordinate
592#         # Get from associated MS
593#         # Get from ASDM
594#         lo1 = -1.
595#         if self.lo1 > 0.:
596#             asaplog.push("Using user defined LO1 frequency %e16.12 [Hz]" % self.lo1)
597#             lo1 = self.lo1
598#         else:
599#             print "NOT IMPLEMENTED YET!!!"
600
601#     def save(self, outfile, outform="ASAP", overwrite=False):
602#         if not overwrite and os.path.exists(outfile):
603#             raise RuntimeError, "Output file '%s' already exists" % outfile
604#
605#         #self.separator._save(outfile, outform)
606
607#     def done(self):
608#         self.close()
609
610#     def close(self):
611#         pass
612#         #del self.separator
613   
614
615
616########################################################################
617    def _Deconvolution(self, data_array, shift, threshold=0.00000001):
618        FObs = []
619        Reject = 0
620        nshift, nchan = data_array.shape
621        nspec = nshift*(nshift-1)/2
622        ifftObs  = numpy.zeros((nspec, nchan), numpy.float)
623        for i in range(nshift):
624           F = FFT.fft(data_array[i])
625           FObs.append(F)
626        z = 0
627        for i in range(nshift):
628            for j in range(i+1, nshift):
629                Fobs = (FObs[i]+FObs[j])/2.0
630                dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
631                #print 'dX,i,j=',dX,i,j
632                for k in range(1,self.nchan):
633                    if math.fabs(math.sin(dX*k)) > threshold:
634                        Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
635                    else: Reject += 1
636                ifftObs[z] = FFT.ifft(Fobs)
637                z += 1
638        print 'Threshold=%s Reject=%d' % (threshold, Reject)
639        return ifftObs
640
641    def _combineResult(self, ifftObs):
642        nspec = len(ifftObs)
643        sum = ifftObs[0]
644        for i in range(1,nspec):
645            sum += ifftObs[i]
646        return(sum/float(nspec))
647
648    def _subtractOtherSide(self, data_array, shift, Data):
649        sum = numpy.zeros(len(Data), numpy.float)
650        numSP = len(data_array)
651        for i in range(numSP):
652            SPsub = data_array[i] - Data
653            sum += self._shiftSpectrum(SPsub, -shift[i])
654        return(sum/float(numSP))
655
656    def _shiftSpectrum(self, data, Shift):
657        Out = numpy.zeros(self.nchan, numpy.float)
658        w2 = Shift % 1
659        w1 = 1.0 - w2
660        for i in range(self.nchan):
661            c1 = int((Shift + i) % self.nchan)
662            c2 = (c1 + 1) % self.nchan
663            Out[c1] += data[i] * w1
664            Out[c2] += data[i] * w2
665        return Out.copy()
Note: See TracBrowser for help on using the repository browser.