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