source: trunk/python/sbseparator.py@ 2770

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

New Development: No

JIRA Issue: Yes (CAS-4141)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: added temporal methods to set scantable to STSideBandSep for a workaround.

Test Programs:

Put in Release Notes: No

Module(s): STSideBandSep and sbseparator

Description:

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