source: trunk/python/sbseparator.py@ 2711

Last change on this file since 2711 was 2711, checked in by Kana Sugimoto, 12 years ago

New Development: No

JIRA Issue: Yes (CAS-4141/CSV-1526)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: Added asap.sbseparator.set_lo1root() = SBSeparator::set_lo1root(),

and STSideBandSep::setLO1Root()

Test Programs:

Put in Release Notes: No

Module(s): sbseparator and STSideBandSep

Description:

Added a new python method asap.sbseparator.set_lo1root() to set MS
name to search for LO1 frequency.
Underlying c++ functions are SBSeparator::set_lo1root (interface) and STSideBandSep::setLO1Root().
Also developed a capability to get relevant table names from scantable header and
fill frequencies of image side band.
File size: 24.5 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 = SBSeparator()
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 lo1val = -1.
224 if isinstance(lo1, dict) and lo1["unit"] == "Hz":
225 lo1val = lo1["value"]
226 else:
227 lo1val = float(lo1)
228 if lo1val <= 0.:
229 asaplog.push("Got negative LO1 frequency. It will be ignored.")
230 asaplog.post("WARN")
231 else:
232 self.separator.set_lo1(lo1val)
233
234
235 def set_lo1root(self, name):
236 """
237 Set MS name which stores LO1 frequency of signal side band.
238 It is used to calculate frequency of image sideband.
239
240 name : MS name which contains 'ASDM_SPECTRALWINDOW' and
241 'ASDM_RECEIVER' tables.
242 """
243 self.separator.set_lo1root(name)
244
245
246 @asaplog_post_dec
247 def separate(self, outname="", overwrite=False):
248 """
249 Invoke sideband separation.
250
251 outname : a name of output scantable
252 overwrite : overwrite existing table
253 """
254 # List up valid scantables and IFNOs to convolve.
255 #self.separator._separate()
256 self._setup_shift()
257 #self._preprocess_tables()
258
259 nshift = len(self.tables)
260 signaltab = self._grid_outtable(self.tables[0])
261 if self.getboth:
262 imagetab = signaltab.copy()
263
264 rejrow = []
265 for irow in xrange(signaltab.nrow()):
266 currpol = signaltab.getpol(irow)
267 currbeam = signaltab.getbeam(irow)
268 currdir = signaltab.get_directionval(irow)
269 spec_array, tabidx = self._get_specarray(polid=currpol,\
270 beamid=currbeam,\
271 dir=currdir)
272 #if not spec_array:
273 if len(tabidx) == 0:
274 asaplog.post()
275 asaplog.push("skipping row %d" % irow)
276 rejrow.append(irow)
277 continue
278 signal = self._solve_signal(spec_array, tabidx)
279 signaltab.set_spectrum(signal, irow)
280
281 # Solve image side side band
282 if self.getboth:
283 image = self._solve_image(spec_array, tabidx)
284 imagetab.set_spectrum(image, irow)
285
286 # TODO: Need to remove rejrow form scantables here
287 signaltab.flag_row(rejrow)
288 if self.getboth:
289 imagetab.flag_row(rejrow)
290
291 if outname == "":
292 outname = "sbsepareted.asap"
293 signalname = outname + ".signalband"
294 if os.path.exists(signalname):
295 if not overwrite:
296 raise Exception, "Output file '%s' exists." % signalname
297 else:
298 shutil.rmtree(signalname)
299 signaltab.save(signalname)
300
301 if self.getboth:
302 # Warnings
303 asaplog.post()
304 asaplog.push("Saving IMAGE sideband.")
305 #asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")
306 #asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")
307 #asaplog.post("WARN")
308
309 imagename = outname + ".imageband"
310 if os.path.exists(imagename):
311 if not overwrite:
312 raise Exception, "Output file '%s' exists." % imagename
313 else:
314 shutil.rmtree(imagename)
315 # Update frequency information
316 self.separator.set_imgtable(imagetab)
317 self.separator.solve_imgfreq()
318 imagetab.save(imagename)
319
320
321 def _solve_signal(self, data, tabidx=None):
322 if not tabidx:
323 tabidx = range(len(data))
324
325 tempshift = []
326 dshift = []
327 if self.solveother:
328 for idx in tabidx:
329 tempshift.append(-self.imageShift[idx])
330 dshift.append(self.signalShift[idx] - self.imageShift[idx])
331 else:
332 for idx in tabidx:
333 tempshift.append(-self.signalShift[idx])
334 dshift.append(self.imageShift[idx] - self.signalShift[idx])
335
336 shiftdata = numpy.zeros(data.shape, numpy.float)
337 for i in range(len(data)):
338 shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
339 ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
340 result_image = self._combineResult(ifftdata)
341 if not self.solveother:
342 return result_image
343 result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
344 return result_signal
345
346
347 def _solve_image(self, data, tabidx=None):
348 if not tabidx:
349 tabidx = range(len(data))
350
351 tempshift = []
352 dshift = []
353 if self.solveother:
354 for idx in tabidx:
355 tempshift.append(-self.signalShift[idx])
356 dshift.append(self.imageShift[idx] - self.signalShift[idx])
357 else:
358 for idx in tabidx:
359 tempshift.append(-self.imageShift[idx])
360 dshift.append(self.signalShift[idx] - self.imageShift[idx])
361
362 shiftdata = numpy.zeros(data.shape, numpy.float)
363 for i in range(len(data)):
364 shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
365 ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
366 result_image = self._combineResult(ifftdata)
367 if not self.solveother:
368 return result_image
369 result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
370 return result_signal
371
372 @asaplog_post_dec
373 def _grid_outtable(self, table):
374 # Generate gridded table for output (Just to get rows)
375 gridder = asapgrid2(table)
376 gridder.setIF(self.baseif)
377
378 cellx = str(self.dirtol[0])+"rad"
379 celly = str(self.dirtol[1])+"rad"
380 dirarr = numpy.array(table.get_directionval()).transpose()
381 mapx = dirarr[0].max() - dirarr[0].min()
382 mapy = dirarr[1].max() - dirarr[1].min()
383 centy = 0.5 * (dirarr[1].max() + dirarr[1].min())
384 nx = max(1, numpy.ceil(mapx*numpy.cos(centy)/self.dirtol[0]))
385 ny = max(1, numpy.ceil(mapy/self.dirtol[0]))
386
387 asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
388 (cellx, celly))
389 gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
390 gridder.setFunc(func='box', convsupport=1)
391 gridder.setWeight(weightType='uniform')
392 gridder.grid()
393 return gridder.getResult()
394
395 @asaplog_post_dec
396 def _get_specarray(self, polid=None, beamid=None, dir=None):
397 ntable = len(self.tables)
398 spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
399 nspec = 0
400 asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
401 tabidx = []
402 for itab in range(ntable):
403 tab = self.tables[itab]
404 # Select rows by POLNO and BEAMNO
405 try:
406 tab.set_selection(pols=[polid], beams=[beamid])
407 if tab.nrow() > 0: tabidx.append(itab)
408 except: # no selection
409 asaplog.post()
410 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
411 asaplog.post("WARN")
412 continue
413
414 # Select rows by direction
415 spec = numpy.zeros(self.nchan, numpy.float)
416 selrow = []
417 for irow in range(tab.nrow()):
418 currdir = tab.get_directionval(irow)
419 if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
420 (abs(currdir[1] - dir[1]) > self.dirtol[1]):
421 continue
422 selrow.append(irow)
423 if len(selrow) == 0:
424 asaplog.post()
425 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
426 asaplog.post("WARN")
427 continue
428 else:
429 seltab = tab.copy()
430 seltab.set_selection(selector(rows=selrow))
431
432 if tab.nrow() > 1:
433 asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
434 tab = seltab.average_time(scanav=False, weight="tintsys")
435 else:
436 tab = seltab
437
438 spec_array[nspec] = tab._getspectrum()
439 nspec += 1
440
441 if nspec != ntable:
442 asaplog.post()
443 #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
444 asaplog.push("Could not find corresponding rows in some tables.")
445 asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
446 if nspec < 2:
447 asaplog.push("At least 2 spectra are necessary for convolution")
448 asaplog.post("ERROR")
449 return False, tabidx
450
451 return spec_array[0:nspec], tabidx
452
453
454 @asaplog_post_dec
455 def _setup_shift(self):
456 ### define self.tables, self.signalShift, and self.imageShift
457 if not self.intables:
458 asaplog.post()
459 raise RuntimeError, "Input data is not defined."
460 #if self.baseif < 0:
461 # asaplog.post()
462 # raise RuntimeError, "Reference IFNO is not defined."
463
464 byname = False
465 #if not self.intables:
466 if isinstance(self.intables[0], str):
467 # A list of file name is given
468 if not os.path.exists(self.intables[0]):
469 asaplog.post()
470 raise RuntimeError, "Could not find '%s'" % self.intables[0]
471
472 stab = scantable(self.intables[0],average=False)
473 ntab = len(self.intables)
474 byname = True
475 else:
476 stab = self.intables[0]
477 ntab = len(self.intables)
478
479 if len(stab.getbeamnos()) > 1:
480 asaplog.post()
481 asaplog.push("Mult-beam data is not supported by this module.")
482 asaplog.post("ERROR")
483 return
484
485 valid_ifs = stab.getifnos()
486 if self.baseif < 0:
487 self.baseif = valid_ifs[0]
488 asaplog.post()
489 asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
490
491 if not (self.baseif in valid_ifs):
492 asaplog.post()
493 errmsg = "IF%d does not exist in the first scantable" % \
494 self.baseif
495 raise RuntimeError, errmsg
496
497 asaplog.push("Start selecting tables and IFNOs to solve.")
498 asaplog.push("Checking frequency of the reference IF")
499 unit_org = stab.get_unit()
500 coord = stab._getcoordinfo()
501 frame_org = coord[1]
502 stab.set_unit("Hz")
503 if len(self.freqframe) > 0:
504 stab.set_freqframe(self.freqframe)
505 stab.set_selection(ifs=[self.baseif])
506 spx = stab._getabcissa()
507 stab.set_selection()
508 basech0 = spx[0]
509 baseinc = spx[1]-spx[0]
510 self.nchan = len(spx)
511 # frequency tolerance
512 if isinstance(self.freqtol, dict) and self.freqtol['unit'] == "Hz":
513 vftol = abs(self.freqtol['value'])
514 else:
515 vftol = abs(baseinc * float(self.freqtol))
516 self.freqtol = dict(value=vftol, unit="Hz")
517 # tolerance of frequency increment
518 inctol = abs(baseinc/float(self.nchan))
519 asaplog.push("Reference frequency setup (Table = 0, IFNO = %d): nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
520 asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
521 poltype0 = stab.poltype()
522
523 self.tables = []
524 self.signalShift = []
525 if self.dsbmode:
526 self.imageShift = []
527
528 for itab in range(ntab):
529 asaplog.push("Table %d:" % itab)
530 tab_selected = False
531 if itab > 0:
532 if byname:
533 stab = scantable(self.intables[itab],average=False)
534 else:
535 stab = self.intables[itab]
536 unit_org = stab.get_unit()
537 coord = stab._getcoordinfo()
538 frame_org = coord[1]
539 stab.set_unit("Hz")
540 if len(self.freqframe) > 0:
541 stab.set_freqframe(self.freqframe)
542
543 # Check POLTYPE should be equal to the first table.
544 if stab.poltype() != poltype0:
545 asaplog.post()
546 raise Exception, "POLTYPE should be equal to the first table."
547 # Multiple beam data may not handled properly
548 if len(stab.getbeamnos()) > 1:
549 asaplog.post()
550 asaplog.push("table contains multiple beams. It may not be handled properly.")
551 asaplog.push("WARN")
552
553 for ifno in stab.getifnos():
554 stab.set_selection(ifs=[ifno])
555 spx = stab._getabcissa()
556 if (len(spx) != self.nchan) or \
557 (abs(spx[0]-basech0) > vftol) or \
558 (abs(spx[1]-spx[0]-baseinc) > inctol):
559 continue
560 tab_selected = True
561 seltab = stab.copy()
562 seltab.set_unit(unit_org)
563 seltab.set_freqframe(frame_org)
564 self.tables.append(seltab)
565 self.signalShift.append((spx[0]-basech0)/baseinc)
566 if self.dsbmode:
567 self.imageShift.append(-self.signalShift[-1])
568 asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
569 stab.set_selection()
570 stab.set_unit(unit_org)
571 stab.set_freqframe(frame_org)
572 if not tab_selected:
573 asaplog.post()
574 asaplog.push("No data selected in Table %d" % itab)
575 asaplog.post("WARN")
576
577 asaplog.push("Total number of IFs selected = %d" % len(self.tables))
578 if len(self.tables) < 2:
579 asaplog.post()
580 raise RuntimeError, "At least 2 IFs are necessary for convolution!"
581
582 if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
583 asaplog.post()
584 errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
585 errmsg += "\nThe frequency tolerance (freqtol) you set may be too small."
586 raise RuntimeError, errmsg
587
588 self.signalShift = numpy.array(self.signalShift)
589 self.imageShift = numpy.array(self.imageShift)
590 self.nshift = len(self.tables)
591
592 @asaplog_post_dec
593 def _preprocess_tables(self):
594 ### temporary method to preprocess data
595 ### Do time averaging for now.
596 for itab in range(len(self.tables)):
597 self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
598
599# @asaplog_post_dec
600# def _setup_image_freq(self, table):
601# # User defined coordinate
602# # Get from associated MS
603# # Get from ASDM
604# lo1 = -1.
605# if self.lo1 > 0.:
606# asaplog.push("Using user defined LO1 frequency %e16.12 [Hz]" % self.lo1)
607# lo1 = self.lo1
608# else:
609# print "NOT IMPLEMENTED YET!!!"
610
611# def save(self, outfile, outform="ASAP", overwrite=False):
612# if not overwrite and os.path.exists(outfile):
613# raise RuntimeError, "Output file '%s' already exists" % outfile
614#
615# #self.separator._save(outfile, outform)
616
617# def done(self):
618# self.close()
619
620# def close(self):
621# pass
622# #del self.separator
623
624
625
626########################################################################
627 def _Deconvolution(self, data_array, shift, threshold=0.00000001):
628 FObs = []
629 Reject = 0
630 nshift, nchan = data_array.shape
631 nspec = nshift*(nshift-1)/2
632 ifftObs = numpy.zeros((nspec, nchan), numpy.float)
633 for i in range(nshift):
634 F = FFT.fft(data_array[i])
635 FObs.append(F)
636 z = 0
637 for i in range(nshift):
638 for j in range(i+1, nshift):
639 Fobs = (FObs[i]+FObs[j])/2.0
640 dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
641 #print 'dX,i,j=',dX,i,j
642 for k in range(1,self.nchan):
643 if math.fabs(math.sin(dX*k)) > threshold:
644 Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
645 else: Reject += 1
646 ifftObs[z] = FFT.ifft(Fobs)
647 z += 1
648 print 'Threshold=%s Reject=%d' % (threshold, Reject)
649 return ifftObs
650
651 def _combineResult(self, ifftObs):
652 nspec = len(ifftObs)
653 sum = ifftObs[0]
654 for i in range(1,nspec):
655 sum += ifftObs[i]
656 return(sum/float(nspec))
657
658 def _subtractOtherSide(self, data_array, shift, Data):
659 sum = numpy.zeros(len(Data), numpy.float)
660 numSP = len(data_array)
661 for i in range(numSP):
662 SPsub = data_array[i] - Data
663 sum += self._shiftSpectrum(SPsub, -shift[i])
664 return(sum/float(numSP))
665
666 def _shiftSpectrum(self, data, Shift):
667 Out = numpy.zeros(self.nchan, numpy.float)
668 w2 = Shift % 1
669 w1 = 1.0 - w2
670 for i in range(self.nchan):
671 c1 = int((Shift + i) % self.nchan)
672 c2 = (c1 + 1) % self.nchan
673 Out[c1] += data[i] * w1
674 Out[c2] += data[i] * w2
675 return Out.copy()
Note: See TracBrowser for help on using the repository browser.