##// END OF EJS Templates
Fix pairsList in SpectraProc
jespinoza -
r1341:34c23d10b2ab
parent child
Show More
@@ -1,874 +1,876
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """Spectra processing Unit and operations
5 """Spectra processing Unit and operations
6
6
7 Here you will find the processing unit `SpectraProc` and several operations
7 Here you will find the processing unit `SpectraProc` and several operations
8 to work with Spectra data type
8 to work with Spectra data type
9 """
9 """
10
10
11 import time
11 import time
12 import itertools
12 import itertools
13
13
14 import numpy
14 import numpy
15
15
16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 from schainpy.model.data.jrodata import Spectra
17 from schainpy.model.data.jrodata import Spectra
18 from schainpy.model.data.jrodata import hildebrand_sekhon
18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 from schainpy.utils import log
19 from schainpy.utils import log
20
20
21
21
22 class SpectraProc(ProcessingUnit):
22 class SpectraProc(ProcessingUnit):
23
23
24 def __init__(self):
24 def __init__(self):
25
25
26 ProcessingUnit.__init__(self)
26 ProcessingUnit.__init__(self)
27
27
28 self.buffer = None
28 self.buffer = None
29 self.firstdatatime = None
29 self.firstdatatime = None
30 self.profIndex = 0
30 self.profIndex = 0
31 self.dataOut = Spectra()
31 self.dataOut = Spectra()
32 self.id_min = None
32 self.id_min = None
33 self.id_max = None
33 self.id_max = None
34 self.setupReq = False #Agregar a todas las unidades de proc
34 self.setupReq = False #Agregar a todas las unidades de proc
35
35
36 def __updateSpecFromVoltage(self):
36 def __updateSpecFromVoltage(self):
37
37
38 self.dataOut.timeZone = self.dataIn.timeZone
38 self.dataOut.timeZone = self.dataIn.timeZone
39 self.dataOut.dstFlag = self.dataIn.dstFlag
39 self.dataOut.dstFlag = self.dataIn.dstFlag
40 self.dataOut.errorCount = self.dataIn.errorCount
40 self.dataOut.errorCount = self.dataIn.errorCount
41 self.dataOut.useLocalTime = self.dataIn.useLocalTime
41 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 try:
42 try:
43 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
43 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
44 except:
44 except:
45 pass
45 pass
46 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
46 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
47
48 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
47 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
49 self.dataOut.channelList = self.dataIn.channelList
48 self.dataOut.channelList = self.dataIn.channelList
50 self.dataOut.heightList = self.dataIn.heightList
49 self.dataOut.heightList = self.dataIn.heightList
51 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
50 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
52 self.dataOut.nProfiles = self.dataOut.nFFTPoints
51 self.dataOut.nProfiles = self.dataOut.nFFTPoints
53 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
52 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
54 self.dataOut.utctime = self.firstdatatime
53 self.dataOut.utctime = self.firstdatatime
55 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
54 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
56 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
55 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
57 self.dataOut.flagShiftFFT = False
56 self.dataOut.flagShiftFFT = False
58 self.dataOut.nCohInt = self.dataIn.nCohInt
57 self.dataOut.nCohInt = self.dataIn.nCohInt
59 self.dataOut.nIncohInt = 1
58 self.dataOut.nIncohInt = 1
60 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
61 self.dataOut.frequency = self.dataIn.frequency
60 self.dataOut.frequency = self.dataIn.frequency
62 self.dataOut.realtime = self.dataIn.realtime
61 self.dataOut.realtime = self.dataIn.realtime
63 self.dataOut.azimuth = self.dataIn.azimuth
62 self.dataOut.azimuth = self.dataIn.azimuth
64 self.dataOut.zenith = self.dataIn.zenith
63 self.dataOut.zenith = self.dataIn.zenith
65 self.dataOut.beam.codeList = self.dataIn.beam.codeList
64 self.dataOut.beam.codeList = self.dataIn.beam.codeList
66 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
65 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
67 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
66 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
68
67
69 def __getFft(self):
68 def __getFft(self):
70 """
69 """
71 Convierte valores de Voltaje a Spectra
70 Convierte valores de Voltaje a Spectra
72
71
73 Affected:
72 Affected:
74 self.dataOut.data_spc
73 self.dataOut.data_spc
75 self.dataOut.data_cspc
74 self.dataOut.data_cspc
76 self.dataOut.data_dc
75 self.dataOut.data_dc
77 self.dataOut.heightList
76 self.dataOut.heightList
78 self.profIndex
77 self.profIndex
79 self.buffer
78 self.buffer
80 self.dataOut.flagNoData
79 self.dataOut.flagNoData
81 """
80 """
82 fft_volt = numpy.fft.fft(
81 fft_volt = numpy.fft.fft(
83 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
82 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
84 fft_volt = fft_volt.astype(numpy.dtype('complex'))
83 fft_volt = fft_volt.astype(numpy.dtype('complex'))
85 dc = fft_volt[:, 0, :]
84 dc = fft_volt[:, 0, :]
86
85
87 # calculo de self-spectra
86 # calculo de self-spectra
88 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
87 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
89 spc = fft_volt * numpy.conjugate(fft_volt)
88 spc = fft_volt * numpy.conjugate(fft_volt)
90 spc = spc.real
89 spc = spc.real
91
90
92 blocksize = 0
91 blocksize = 0
93 blocksize += dc.size
92 blocksize += dc.size
94 blocksize += spc.size
93 blocksize += spc.size
95
94
96 cspc = None
95 cspc = None
97 pairIndex = 0
96 pairIndex = 0
98 if self.dataOut.pairsList != None:
97 if self.dataOut.pairsList != None:
99 # calculo de cross-spectra
98 # calculo de cross-spectra
100 cspc = numpy.zeros(
99 cspc = numpy.zeros(
101 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
100 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
102 for pair in self.dataOut.pairsList:
101 for pair in self.dataOut.pairsList:
103 if pair[0] not in self.dataOut.channelList:
102 if pair[0] not in self.dataOut.channelList:
104 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
103 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
105 str(pair), str(self.dataOut.channelList)))
104 str(pair), str(self.dataOut.channelList)))
106 if pair[1] not in self.dataOut.channelList:
105 if pair[1] not in self.dataOut.channelList:
107 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
106 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
108 str(pair), str(self.dataOut.channelList)))
107 str(pair), str(self.dataOut.channelList)))
109
108
110 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
109 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
111 numpy.conjugate(fft_volt[pair[1], :, :])
110 numpy.conjugate(fft_volt[pair[1], :, :])
112 pairIndex += 1
111 pairIndex += 1
113 blocksize += cspc.size
112 blocksize += cspc.size
114
113
115 self.dataOut.data_spc = spc
114 self.dataOut.data_spc = spc
116 self.dataOut.data_cspc = cspc
115 self.dataOut.data_cspc = cspc
117 self.dataOut.data_dc = dc
116 self.dataOut.data_dc = dc
118 self.dataOut.blockSize = blocksize
117 self.dataOut.blockSize = blocksize
119 self.dataOut.flagShiftFFT = False
118 self.dataOut.flagShiftFFT = False
120
119
121 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
120 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
122
121
123 if self.dataIn.type == "Spectra":
122 if self.dataIn.type == "Spectra":
124 self.dataOut.copy(self.dataIn)
123 self.dataOut.copy(self.dataIn)
125 if shift_fft:
124 if shift_fft:
126 #desplaza a la derecha en el eje 2 determinadas posiciones
125 #desplaza a la derecha en el eje 2 determinadas posiciones
127 shift = int(self.dataOut.nFFTPoints/2)
126 shift = int(self.dataOut.nFFTPoints/2)
128 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
127 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
129
128
130 if self.dataOut.data_cspc is not None:
129 if self.dataOut.data_cspc is not None:
131 #desplaza a la derecha en el eje 2 determinadas posiciones
130 #desplaza a la derecha en el eje 2 determinadas posiciones
132 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
131 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
132 if pairsList:
133 self.__selectPairs(pairsList)
133
134
134 elif self.dataIn.type == "Voltage":
135 elif self.dataIn.type == "Voltage":
135
136
136 self.dataOut.flagNoData = True
137 self.dataOut.flagNoData = True
137
138
138 if nFFTPoints == None:
139 if nFFTPoints == None:
139 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
140 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
140
141
141 if nProfiles == None:
142 if nProfiles == None:
142 nProfiles = nFFTPoints
143 nProfiles = nFFTPoints
143
144
144 if ippFactor == None:
145 if ippFactor == None:
145 self.dataOut.ippFactor = 1
146 self.dataOut.ippFactor = 1
146
147
147 self.dataOut.nFFTPoints = nFFTPoints
148 self.dataOut.nFFTPoints = nFFTPoints
148
149
149 if self.buffer is None:
150 if self.buffer is None:
150 self.buffer = numpy.zeros((self.dataIn.nChannels,
151 self.buffer = numpy.zeros((self.dataIn.nChannels,
151 nProfiles,
152 nProfiles,
152 self.dataIn.nHeights),
153 self.dataIn.nHeights),
153 dtype='complex')
154 dtype='complex')
154
155
155 if self.dataIn.flagDataAsBlock:
156 if self.dataIn.flagDataAsBlock:
156 nVoltProfiles = self.dataIn.data.shape[1]
157 nVoltProfiles = self.dataIn.data.shape[1]
157
158
158 if nVoltProfiles == nProfiles:
159 if nVoltProfiles == nProfiles:
159 self.buffer = self.dataIn.data.copy()
160 self.buffer = self.dataIn.data.copy()
160 self.profIndex = nVoltProfiles
161 self.profIndex = nVoltProfiles
161
162
162 elif nVoltProfiles < nProfiles:
163 elif nVoltProfiles < nProfiles:
163
164
164 if self.profIndex == 0:
165 if self.profIndex == 0:
165 self.id_min = 0
166 self.id_min = 0
166 self.id_max = nVoltProfiles
167 self.id_max = nVoltProfiles
167
168
168 self.buffer[:, self.id_min:self.id_max,
169 self.buffer[:, self.id_min:self.id_max,
169 :] = self.dataIn.data
170 :] = self.dataIn.data
170 self.profIndex += nVoltProfiles
171 self.profIndex += nVoltProfiles
171 self.id_min += nVoltProfiles
172 self.id_min += nVoltProfiles
172 self.id_max += nVoltProfiles
173 self.id_max += nVoltProfiles
173 else:
174 else:
174 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
175 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
175 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
176 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
176 self.dataOut.flagNoData = True
177 self.dataOut.flagNoData = True
177 else:
178 else:
178 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
179 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
179 self.profIndex += 1
180 self.profIndex += 1
180
181
181 if self.firstdatatime == None:
182 if self.firstdatatime == None:
182 self.firstdatatime = self.dataIn.utctime
183 self.firstdatatime = self.dataIn.utctime
183
184
184 if self.profIndex == nProfiles:
185 if self.profIndex == nProfiles:
185 self.__updateSpecFromVoltage()
186 self.__updateSpecFromVoltage()
186 if pairsList == None:
187 if pairsList == None:
187 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
188 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
189 else:
190 self.dataOut.pairsList = pairsList
188 self.__getFft()
191 self.__getFft()
189
190 self.dataOut.flagNoData = False
192 self.dataOut.flagNoData = False
191 self.firstdatatime = None
193 self.firstdatatime = None
192 self.profIndex = 0
194 self.profIndex = 0
193 else:
195 else:
194 raise ValueError("The type of input object '%s' is not valid".format(
196 raise ValueError("The type of input object '%s' is not valid".format(
195 self.dataIn.type))
197 self.dataIn.type))
196
198
197 def __selectPairs(self, pairsList):
199 def __selectPairs(self, pairsList):
198
200
199 if not pairsList:
201 if not pairsList:
200 return
202 return
201
203
202 pairs = []
204 pairs = []
203 pairsIndex = []
205 pairsIndex = []
204
206
205 for pair in pairsList:
207 for pair in pairsList:
206 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
208 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
207 continue
209 continue
208 pairs.append(pair)
210 pairs.append(pair)
209 pairsIndex.append(pairs.index(pair))
211 pairsIndex.append(pairs.index(pair))
210
212
211 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
213 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
212 self.dataOut.pairsList = pairs
214 self.dataOut.pairsList = pairs
213
215
214 return
216 return
215
217
216 def selectFFTs(self, minFFT, maxFFT ):
218 def selectFFTs(self, minFFT, maxFFT ):
217 """
219 """
218 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
220 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
219 minFFT<= FFT <= maxFFT
221 minFFT<= FFT <= maxFFT
220 """
222 """
221
223
222 if (minFFT > maxFFT):
224 if (minFFT > maxFFT):
223 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
225 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
224
226
225 if (minFFT < self.dataOut.getFreqRange()[0]):
227 if (minFFT < self.dataOut.getFreqRange()[0]):
226 minFFT = self.dataOut.getFreqRange()[0]
228 minFFT = self.dataOut.getFreqRange()[0]
227
229
228 if (maxFFT > self.dataOut.getFreqRange()[-1]):
230 if (maxFFT > self.dataOut.getFreqRange()[-1]):
229 maxFFT = self.dataOut.getFreqRange()[-1]
231 maxFFT = self.dataOut.getFreqRange()[-1]
230
232
231 minIndex = 0
233 minIndex = 0
232 maxIndex = 0
234 maxIndex = 0
233 FFTs = self.dataOut.getFreqRange()
235 FFTs = self.dataOut.getFreqRange()
234
236
235 inda = numpy.where(FFTs >= minFFT)
237 inda = numpy.where(FFTs >= minFFT)
236 indb = numpy.where(FFTs <= maxFFT)
238 indb = numpy.where(FFTs <= maxFFT)
237
239
238 try:
240 try:
239 minIndex = inda[0][0]
241 minIndex = inda[0][0]
240 except:
242 except:
241 minIndex = 0
243 minIndex = 0
242
244
243 try:
245 try:
244 maxIndex = indb[0][-1]
246 maxIndex = indb[0][-1]
245 except:
247 except:
246 maxIndex = len(FFTs)
248 maxIndex = len(FFTs)
247
249
248 self.selectFFTsByIndex(minIndex, maxIndex)
250 self.selectFFTsByIndex(minIndex, maxIndex)
249
251
250 return 1
252 return 1
251
253
252 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
254 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
253 newheis = numpy.where(
255 newheis = numpy.where(
254 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
256 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
255
257
256 if hei_ref != None:
258 if hei_ref != None:
257 newheis = numpy.where(self.dataOut.heightList > hei_ref)
259 newheis = numpy.where(self.dataOut.heightList > hei_ref)
258
260
259 minIndex = min(newheis[0])
261 minIndex = min(newheis[0])
260 maxIndex = max(newheis[0])
262 maxIndex = max(newheis[0])
261 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
263 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
262 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
264 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
263
265
264 # determina indices
266 # determina indices
265 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
267 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
266 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
268 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
267 avg_dB = 10 * \
269 avg_dB = 10 * \
268 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
270 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
269 beacon_dB = numpy.sort(avg_dB)[-nheis:]
271 beacon_dB = numpy.sort(avg_dB)[-nheis:]
270 beacon_heiIndexList = []
272 beacon_heiIndexList = []
271 for val in avg_dB.tolist():
273 for val in avg_dB.tolist():
272 if val >= beacon_dB[0]:
274 if val >= beacon_dB[0]:
273 beacon_heiIndexList.append(avg_dB.tolist().index(val))
275 beacon_heiIndexList.append(avg_dB.tolist().index(val))
274
276
275 #data_spc = data_spc[:,:,beacon_heiIndexList]
277 #data_spc = data_spc[:,:,beacon_heiIndexList]
276 data_cspc = None
278 data_cspc = None
277 if self.dataOut.data_cspc is not None:
279 if self.dataOut.data_cspc is not None:
278 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
280 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
279 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
281 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
280
282
281 data_dc = None
283 data_dc = None
282 if self.dataOut.data_dc is not None:
284 if self.dataOut.data_dc is not None:
283 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
285 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
284 #data_dc = data_dc[:,beacon_heiIndexList]
286 #data_dc = data_dc[:,beacon_heiIndexList]
285
287
286 self.dataOut.data_spc = data_spc
288 self.dataOut.data_spc = data_spc
287 self.dataOut.data_cspc = data_cspc
289 self.dataOut.data_cspc = data_cspc
288 self.dataOut.data_dc = data_dc
290 self.dataOut.data_dc = data_dc
289 self.dataOut.heightList = heightList
291 self.dataOut.heightList = heightList
290 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
292 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
291
293
292 return 1
294 return 1
293
295
294 def selectFFTsByIndex(self, minIndex, maxIndex):
296 def selectFFTsByIndex(self, minIndex, maxIndex):
295 """
297 """
296
298
297 """
299 """
298
300
299 if (minIndex < 0) or (minIndex > maxIndex):
301 if (minIndex < 0) or (minIndex > maxIndex):
300 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
302 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
301
303
302 if (maxIndex >= self.dataOut.nProfiles):
304 if (maxIndex >= self.dataOut.nProfiles):
303 maxIndex = self.dataOut.nProfiles-1
305 maxIndex = self.dataOut.nProfiles-1
304
306
305 #Spectra
307 #Spectra
306 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
308 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
307
309
308 data_cspc = None
310 data_cspc = None
309 if self.dataOut.data_cspc is not None:
311 if self.dataOut.data_cspc is not None:
310 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
312 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
311
313
312 data_dc = None
314 data_dc = None
313 if self.dataOut.data_dc is not None:
315 if self.dataOut.data_dc is not None:
314 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
316 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
315
317
316 self.dataOut.data_spc = data_spc
318 self.dataOut.data_spc = data_spc
317 self.dataOut.data_cspc = data_cspc
319 self.dataOut.data_cspc = data_cspc
318 self.dataOut.data_dc = data_dc
320 self.dataOut.data_dc = data_dc
319
321
320 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
322 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
321 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
323 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
322 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
324 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
323
325
324 return 1
326 return 1
325
327
326 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
328 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
327 # validacion de rango
329 # validacion de rango
328 if minHei == None:
330 if minHei == None:
329 minHei = self.dataOut.heightList[0]
331 minHei = self.dataOut.heightList[0]
330
332
331 if maxHei == None:
333 if maxHei == None:
332 maxHei = self.dataOut.heightList[-1]
334 maxHei = self.dataOut.heightList[-1]
333
335
334 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
336 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
335 print('minHei: %.2f is out of the heights range' % (minHei))
337 print('minHei: %.2f is out of the heights range' % (minHei))
336 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
338 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
337 minHei = self.dataOut.heightList[0]
339 minHei = self.dataOut.heightList[0]
338
340
339 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
341 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
340 print('maxHei: %.2f is out of the heights range' % (maxHei))
342 print('maxHei: %.2f is out of the heights range' % (maxHei))
341 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
343 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
342 maxHei = self.dataOut.heightList[-1]
344 maxHei = self.dataOut.heightList[-1]
343
345
344 # validacion de velocidades
346 # validacion de velocidades
345 velrange = self.dataOut.getVelRange(1)
347 velrange = self.dataOut.getVelRange(1)
346
348
347 if minVel == None:
349 if minVel == None:
348 minVel = velrange[0]
350 minVel = velrange[0]
349
351
350 if maxVel == None:
352 if maxVel == None:
351 maxVel = velrange[-1]
353 maxVel = velrange[-1]
352
354
353 if (minVel < velrange[0]) or (minVel > maxVel):
355 if (minVel < velrange[0]) or (minVel > maxVel):
354 print('minVel: %.2f is out of the velocity range' % (minVel))
356 print('minVel: %.2f is out of the velocity range' % (minVel))
355 print('minVel is setting to %.2f' % (velrange[0]))
357 print('minVel is setting to %.2f' % (velrange[0]))
356 minVel = velrange[0]
358 minVel = velrange[0]
357
359
358 if (maxVel > velrange[-1]) or (maxVel < minVel):
360 if (maxVel > velrange[-1]) or (maxVel < minVel):
359 print('maxVel: %.2f is out of the velocity range' % (maxVel))
361 print('maxVel: %.2f is out of the velocity range' % (maxVel))
360 print('maxVel is setting to %.2f' % (velrange[-1]))
362 print('maxVel is setting to %.2f' % (velrange[-1]))
361 maxVel = velrange[-1]
363 maxVel = velrange[-1]
362
364
363 # seleccion de indices para rango
365 # seleccion de indices para rango
364 minIndex = 0
366 minIndex = 0
365 maxIndex = 0
367 maxIndex = 0
366 heights = self.dataOut.heightList
368 heights = self.dataOut.heightList
367
369
368 inda = numpy.where(heights >= minHei)
370 inda = numpy.where(heights >= minHei)
369 indb = numpy.where(heights <= maxHei)
371 indb = numpy.where(heights <= maxHei)
370
372
371 try:
373 try:
372 minIndex = inda[0][0]
374 minIndex = inda[0][0]
373 except:
375 except:
374 minIndex = 0
376 minIndex = 0
375
377
376 try:
378 try:
377 maxIndex = indb[0][-1]
379 maxIndex = indb[0][-1]
378 except:
380 except:
379 maxIndex = len(heights)
381 maxIndex = len(heights)
380
382
381 if (minIndex < 0) or (minIndex > maxIndex):
383 if (minIndex < 0) or (minIndex > maxIndex):
382 raise ValueError("some value in (%d,%d) is not valid" % (
384 raise ValueError("some value in (%d,%d) is not valid" % (
383 minIndex, maxIndex))
385 minIndex, maxIndex))
384
386
385 if (maxIndex >= self.dataOut.nHeights):
387 if (maxIndex >= self.dataOut.nHeights):
386 maxIndex = self.dataOut.nHeights - 1
388 maxIndex = self.dataOut.nHeights - 1
387
389
388 # seleccion de indices para velocidades
390 # seleccion de indices para velocidades
389 indminvel = numpy.where(velrange >= minVel)
391 indminvel = numpy.where(velrange >= minVel)
390 indmaxvel = numpy.where(velrange <= maxVel)
392 indmaxvel = numpy.where(velrange <= maxVel)
391 try:
393 try:
392 minIndexVel = indminvel[0][0]
394 minIndexVel = indminvel[0][0]
393 except:
395 except:
394 minIndexVel = 0
396 minIndexVel = 0
395
397
396 try:
398 try:
397 maxIndexVel = indmaxvel[0][-1]
399 maxIndexVel = indmaxvel[0][-1]
398 except:
400 except:
399 maxIndexVel = len(velrange)
401 maxIndexVel = len(velrange)
400
402
401 # seleccion del espectro
403 # seleccion del espectro
402 data_spc = self.dataOut.data_spc[:,
404 data_spc = self.dataOut.data_spc[:,
403 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
405 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
404 # estimacion de ruido
406 # estimacion de ruido
405 noise = numpy.zeros(self.dataOut.nChannels)
407 noise = numpy.zeros(self.dataOut.nChannels)
406
408
407 for channel in range(self.dataOut.nChannels):
409 for channel in range(self.dataOut.nChannels):
408 daux = data_spc[channel, :, :]
410 daux = data_spc[channel, :, :]
409 sortdata = numpy.sort(daux, axis=None)
411 sortdata = numpy.sort(daux, axis=None)
410 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
412 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
411
413
412 self.dataOut.noise_estimation = noise.copy()
414 self.dataOut.noise_estimation = noise.copy()
413
415
414 return 1
416 return 1
415
417
416 class removeDC(Operation):
418 class removeDC(Operation):
417
419
418 def run(self, dataOut, mode=2):
420 def run(self, dataOut, mode=2):
419 self.dataOut = dataOut
421 self.dataOut = dataOut
420 jspectra = self.dataOut.data_spc
422 jspectra = self.dataOut.data_spc
421 jcspectra = self.dataOut.data_cspc
423 jcspectra = self.dataOut.data_cspc
422
424
423 num_chan = jspectra.shape[0]
425 num_chan = jspectra.shape[0]
424 num_hei = jspectra.shape[2]
426 num_hei = jspectra.shape[2]
425
427
426 if jcspectra is not None:
428 if jcspectra is not None:
427 jcspectraExist = True
429 jcspectraExist = True
428 num_pairs = jcspectra.shape[0]
430 num_pairs = jcspectra.shape[0]
429 else:
431 else:
430 jcspectraExist = False
432 jcspectraExist = False
431
433
432 freq_dc = int(jspectra.shape[1] / 2)
434 freq_dc = int(jspectra.shape[1] / 2)
433 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
435 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
434 ind_vel = ind_vel.astype(int)
436 ind_vel = ind_vel.astype(int)
435
437
436 if ind_vel[0] < 0:
438 if ind_vel[0] < 0:
437 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
439 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
438
440
439 if mode == 1:
441 if mode == 1:
440 jspectra[:, freq_dc, :] = (
442 jspectra[:, freq_dc, :] = (
441 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
443 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
442
444
443 if jcspectraExist:
445 if jcspectraExist:
444 jcspectra[:, freq_dc, :] = (
446 jcspectra[:, freq_dc, :] = (
445 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
447 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
446
448
447 if mode == 2:
449 if mode == 2:
448
450
449 vel = numpy.array([-2, -1, 1, 2])
451 vel = numpy.array([-2, -1, 1, 2])
450 xx = numpy.zeros([4, 4])
452 xx = numpy.zeros([4, 4])
451
453
452 for fil in range(4):
454 for fil in range(4):
453 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
455 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
454
456
455 xx_inv = numpy.linalg.inv(xx)
457 xx_inv = numpy.linalg.inv(xx)
456 xx_aux = xx_inv[0, :]
458 xx_aux = xx_inv[0, :]
457
459
458 for ich in range(num_chan):
460 for ich in range(num_chan):
459 yy = jspectra[ich, ind_vel, :]
461 yy = jspectra[ich, ind_vel, :]
460 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
462 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
461
463
462 junkid = jspectra[ich, freq_dc, :] <= 0
464 junkid = jspectra[ich, freq_dc, :] <= 0
463 cjunkid = sum(junkid)
465 cjunkid = sum(junkid)
464
466
465 if cjunkid.any():
467 if cjunkid.any():
466 jspectra[ich, freq_dc, junkid.nonzero()] = (
468 jspectra[ich, freq_dc, junkid.nonzero()] = (
467 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
469 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
468
470
469 if jcspectraExist:
471 if jcspectraExist:
470 for ip in range(num_pairs):
472 for ip in range(num_pairs):
471 yy = jcspectra[ip, ind_vel, :]
473 yy = jcspectra[ip, ind_vel, :]
472 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
474 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
473
475
474 self.dataOut.data_spc = jspectra
476 self.dataOut.data_spc = jspectra
475 self.dataOut.data_cspc = jcspectra
477 self.dataOut.data_cspc = jcspectra
476
478
477 return self.dataOut
479 return self.dataOut
478
480
479 class removeInterference(Operation):
481 class removeInterference(Operation):
480
482
481 def removeInterference2(self):
483 def removeInterference2(self):
482
484
483 cspc = self.dataOut.data_cspc
485 cspc = self.dataOut.data_cspc
484 spc = self.dataOut.data_spc
486 spc = self.dataOut.data_spc
485 Heights = numpy.arange(cspc.shape[2])
487 Heights = numpy.arange(cspc.shape[2])
486 realCspc = numpy.abs(cspc)
488 realCspc = numpy.abs(cspc)
487
489
488 for i in range(cspc.shape[0]):
490 for i in range(cspc.shape[0]):
489 LinePower= numpy.sum(realCspc[i], axis=0)
491 LinePower= numpy.sum(realCspc[i], axis=0)
490 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
492 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
491 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
493 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
492 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
494 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
493 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
495 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
494 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
496 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
495
497
496
498
497 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
499 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
498 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
500 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
499 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
501 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
500 cspc[i,InterferenceRange,:] = numpy.NaN
502 cspc[i,InterferenceRange,:] = numpy.NaN
501
503
502 self.dataOut.data_cspc = cspc
504 self.dataOut.data_cspc = cspc
503
505
504 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
506 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
505
507
506 jspectra = self.dataOut.data_spc
508 jspectra = self.dataOut.data_spc
507 jcspectra = self.dataOut.data_cspc
509 jcspectra = self.dataOut.data_cspc
508 jnoise = self.dataOut.getNoise()
510 jnoise = self.dataOut.getNoise()
509 num_incoh = self.dataOut.nIncohInt
511 num_incoh = self.dataOut.nIncohInt
510
512
511 num_channel = jspectra.shape[0]
513 num_channel = jspectra.shape[0]
512 num_prof = jspectra.shape[1]
514 num_prof = jspectra.shape[1]
513 num_hei = jspectra.shape[2]
515 num_hei = jspectra.shape[2]
514
516
515 # hei_interf
517 # hei_interf
516 if hei_interf is None:
518 if hei_interf is None:
517 count_hei = int(num_hei / 2)
519 count_hei = int(num_hei / 2)
518 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
520 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
519 hei_interf = numpy.asarray(hei_interf)[0]
521 hei_interf = numpy.asarray(hei_interf)[0]
520 # nhei_interf
522 # nhei_interf
521 if (nhei_interf == None):
523 if (nhei_interf == None):
522 nhei_interf = 5
524 nhei_interf = 5
523 if (nhei_interf < 1):
525 if (nhei_interf < 1):
524 nhei_interf = 1
526 nhei_interf = 1
525 if (nhei_interf > count_hei):
527 if (nhei_interf > count_hei):
526 nhei_interf = count_hei
528 nhei_interf = count_hei
527 if (offhei_interf == None):
529 if (offhei_interf == None):
528 offhei_interf = 0
530 offhei_interf = 0
529
531
530 ind_hei = list(range(num_hei))
532 ind_hei = list(range(num_hei))
531 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
533 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
532 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
534 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
533 mask_prof = numpy.asarray(list(range(num_prof)))
535 mask_prof = numpy.asarray(list(range(num_prof)))
534 num_mask_prof = mask_prof.size
536 num_mask_prof = mask_prof.size
535 comp_mask_prof = [0, num_prof / 2]
537 comp_mask_prof = [0, num_prof / 2]
536
538
537 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
539 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
538 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
540 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
539 jnoise = numpy.nan
541 jnoise = numpy.nan
540 noise_exist = jnoise[0] < numpy.Inf
542 noise_exist = jnoise[0] < numpy.Inf
541
543
542 # Subrutina de Remocion de la Interferencia
544 # Subrutina de Remocion de la Interferencia
543 for ich in range(num_channel):
545 for ich in range(num_channel):
544 # Se ordena los espectros segun su potencia (menor a mayor)
546 # Se ordena los espectros segun su potencia (menor a mayor)
545 power = jspectra[ich, mask_prof, :]
547 power = jspectra[ich, mask_prof, :]
546 power = power[:, hei_interf]
548 power = power[:, hei_interf]
547 power = power.sum(axis=0)
549 power = power.sum(axis=0)
548 psort = power.ravel().argsort()
550 psort = power.ravel().argsort()
549
551
550 # Se estima la interferencia promedio en los Espectros de Potencia empleando
552 # Se estima la interferencia promedio en los Espectros de Potencia empleando
551 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
553 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
552 offhei_interf, nhei_interf + offhei_interf))]]]
554 offhei_interf, nhei_interf + offhei_interf))]]]
553
555
554 if noise_exist:
556 if noise_exist:
555 # tmp_noise = jnoise[ich] / num_prof
557 # tmp_noise = jnoise[ich] / num_prof
556 tmp_noise = jnoise[ich]
558 tmp_noise = jnoise[ich]
557 junkspc_interf = junkspc_interf - tmp_noise
559 junkspc_interf = junkspc_interf - tmp_noise
558 #junkspc_interf[:,comp_mask_prof] = 0
560 #junkspc_interf[:,comp_mask_prof] = 0
559
561
560 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
562 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
561 jspc_interf = jspc_interf.transpose()
563 jspc_interf = jspc_interf.transpose()
562 # Calculando el espectro de interferencia promedio
564 # Calculando el espectro de interferencia promedio
563 noiseid = numpy.where(
565 noiseid = numpy.where(
564 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
566 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
565 noiseid = noiseid[0]
567 noiseid = noiseid[0]
566 cnoiseid = noiseid.size
568 cnoiseid = noiseid.size
567 interfid = numpy.where(
569 interfid = numpy.where(
568 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
570 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
569 interfid = interfid[0]
571 interfid = interfid[0]
570 cinterfid = interfid.size
572 cinterfid = interfid.size
571
573
572 if (cnoiseid > 0):
574 if (cnoiseid > 0):
573 jspc_interf[noiseid] = 0
575 jspc_interf[noiseid] = 0
574
576
575 # Expandiendo los perfiles a limpiar
577 # Expandiendo los perfiles a limpiar
576 if (cinterfid > 0):
578 if (cinterfid > 0):
577 new_interfid = (
579 new_interfid = (
578 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
580 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
579 new_interfid = numpy.asarray(new_interfid)
581 new_interfid = numpy.asarray(new_interfid)
580 new_interfid = {x for x in new_interfid}
582 new_interfid = {x for x in new_interfid}
581 new_interfid = numpy.array(list(new_interfid))
583 new_interfid = numpy.array(list(new_interfid))
582 new_cinterfid = new_interfid.size
584 new_cinterfid = new_interfid.size
583 else:
585 else:
584 new_cinterfid = 0
586 new_cinterfid = 0
585
587
586 for ip in range(new_cinterfid):
588 for ip in range(new_cinterfid):
587 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
589 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
588 jspc_interf[new_interfid[ip]
590 jspc_interf[new_interfid[ip]
589 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
591 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
590
592
591 jspectra[ich, :, ind_hei] = jspectra[ich, :,
593 jspectra[ich, :, ind_hei] = jspectra[ich, :,
592 ind_hei] - jspc_interf # Corregir indices
594 ind_hei] - jspc_interf # Corregir indices
593
595
594 # Removiendo la interferencia del punto de mayor interferencia
596 # Removiendo la interferencia del punto de mayor interferencia
595 ListAux = jspc_interf[mask_prof].tolist()
597 ListAux = jspc_interf[mask_prof].tolist()
596 maxid = ListAux.index(max(ListAux))
598 maxid = ListAux.index(max(ListAux))
597
599
598 if cinterfid > 0:
600 if cinterfid > 0:
599 for ip in range(cinterfid * (interf == 2) - 1):
601 for ip in range(cinterfid * (interf == 2) - 1):
600 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
602 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
601 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
603 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
602 cind = len(ind)
604 cind = len(ind)
603
605
604 if (cind > 0):
606 if (cind > 0):
605 jspectra[ich, interfid[ip], ind] = tmp_noise * \
607 jspectra[ich, interfid[ip], ind] = tmp_noise * \
606 (1 + (numpy.random.uniform(cind) - 0.5) /
608 (1 + (numpy.random.uniform(cind) - 0.5) /
607 numpy.sqrt(num_incoh))
609 numpy.sqrt(num_incoh))
608
610
609 ind = numpy.array([-2, -1, 1, 2])
611 ind = numpy.array([-2, -1, 1, 2])
610 xx = numpy.zeros([4, 4])
612 xx = numpy.zeros([4, 4])
611
613
612 for id1 in range(4):
614 for id1 in range(4):
613 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
615 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
614
616
615 xx_inv = numpy.linalg.inv(xx)
617 xx_inv = numpy.linalg.inv(xx)
616 xx = xx_inv[:, 0]
618 xx = xx_inv[:, 0]
617 ind = (ind + maxid + num_mask_prof) % num_mask_prof
619 ind = (ind + maxid + num_mask_prof) % num_mask_prof
618 yy = jspectra[ich, mask_prof[ind], :]
620 yy = jspectra[ich, mask_prof[ind], :]
619 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
621 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
620 yy.transpose(), xx)
622 yy.transpose(), xx)
621
623
622 indAux = (jspectra[ich, :, :] < tmp_noise *
624 indAux = (jspectra[ich, :, :] < tmp_noise *
623 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
625 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
624 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
626 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
625 (1 - 1 / numpy.sqrt(num_incoh))
627 (1 - 1 / numpy.sqrt(num_incoh))
626
628
627 # Remocion de Interferencia en el Cross Spectra
629 # Remocion de Interferencia en el Cross Spectra
628 if jcspectra is None:
630 if jcspectra is None:
629 return jspectra, jcspectra
631 return jspectra, jcspectra
630 num_pairs = int(jcspectra.size / (num_prof * num_hei))
632 num_pairs = int(jcspectra.size / (num_prof * num_hei))
631 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
633 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
632
634
633 for ip in range(num_pairs):
635 for ip in range(num_pairs):
634
636
635 #-------------------------------------------
637 #-------------------------------------------
636
638
637 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
639 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
638 cspower = cspower[:, hei_interf]
640 cspower = cspower[:, hei_interf]
639 cspower = cspower.sum(axis=0)
641 cspower = cspower.sum(axis=0)
640
642
641 cspsort = cspower.ravel().argsort()
643 cspsort = cspower.ravel().argsort()
642 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
644 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
643 offhei_interf, nhei_interf + offhei_interf))]]]
645 offhei_interf, nhei_interf + offhei_interf))]]]
644 junkcspc_interf = junkcspc_interf.transpose()
646 junkcspc_interf = junkcspc_interf.transpose()
645 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
647 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
646
648
647 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
649 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
648
650
649 median_real = int(numpy.median(numpy.real(
651 median_real = int(numpy.median(numpy.real(
650 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
652 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
651 median_imag = int(numpy.median(numpy.imag(
653 median_imag = int(numpy.median(numpy.imag(
652 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
654 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
653 comp_mask_prof = [int(e) for e in comp_mask_prof]
655 comp_mask_prof = [int(e) for e in comp_mask_prof]
654 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
656 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
655 median_real, median_imag)
657 median_real, median_imag)
656
658
657 for iprof in range(num_prof):
659 for iprof in range(num_prof):
658 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
660 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
659 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
661 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
660
662
661 # Removiendo la Interferencia
663 # Removiendo la Interferencia
662 jcspectra[ip, :, ind_hei] = jcspectra[ip,
664 jcspectra[ip, :, ind_hei] = jcspectra[ip,
663 :, ind_hei] - jcspc_interf
665 :, ind_hei] - jcspc_interf
664
666
665 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
667 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
666 maxid = ListAux.index(max(ListAux))
668 maxid = ListAux.index(max(ListAux))
667
669
668 ind = numpy.array([-2, -1, 1, 2])
670 ind = numpy.array([-2, -1, 1, 2])
669 xx = numpy.zeros([4, 4])
671 xx = numpy.zeros([4, 4])
670
672
671 for id1 in range(4):
673 for id1 in range(4):
672 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
674 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
673
675
674 xx_inv = numpy.linalg.inv(xx)
676 xx_inv = numpy.linalg.inv(xx)
675 xx = xx_inv[:, 0]
677 xx = xx_inv[:, 0]
676
678
677 ind = (ind + maxid + num_mask_prof) % num_mask_prof
679 ind = (ind + maxid + num_mask_prof) % num_mask_prof
678 yy = jcspectra[ip, mask_prof[ind], :]
680 yy = jcspectra[ip, mask_prof[ind], :]
679 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
681 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
680
682
681 # Guardar Resultados
683 # Guardar Resultados
682 self.dataOut.data_spc = jspectra
684 self.dataOut.data_spc = jspectra
683 self.dataOut.data_cspc = jcspectra
685 self.dataOut.data_cspc = jcspectra
684
686
685 return 1
687 return 1
686
688
687 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
689 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
688
690
689 self.dataOut = dataOut
691 self.dataOut = dataOut
690
692
691 if mode == 1:
693 if mode == 1:
692 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
694 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
693 elif mode == 2:
695 elif mode == 2:
694 self.removeInterference2()
696 self.removeInterference2()
695
697
696 return self.dataOut
698 return self.dataOut
697
699
698
700
699 class IncohInt(Operation):
701 class IncohInt(Operation):
700
702
701 __profIndex = 0
703 __profIndex = 0
702 __withOverapping = False
704 __withOverapping = False
703
705
704 __byTime = False
706 __byTime = False
705 __initime = None
707 __initime = None
706 __lastdatatime = None
708 __lastdatatime = None
707 __integrationtime = None
709 __integrationtime = None
708
710
709 __buffer_spc = None
711 __buffer_spc = None
710 __buffer_cspc = None
712 __buffer_cspc = None
711 __buffer_dc = None
713 __buffer_dc = None
712
714
713 __dataReady = False
715 __dataReady = False
714
716
715 __timeInterval = None
717 __timeInterval = None
716
718
717 n = None
719 n = None
718
720
719 def __init__(self):
721 def __init__(self):
720
722
721 Operation.__init__(self)
723 Operation.__init__(self)
722
724
723 def setup(self, n=None, timeInterval=None, overlapping=False):
725 def setup(self, n=None, timeInterval=None, overlapping=False):
724 """
726 """
725 Set the parameters of the integration class.
727 Set the parameters of the integration class.
726
728
727 Inputs:
729 Inputs:
728
730
729 n : Number of coherent integrations
731 n : Number of coherent integrations
730 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
732 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
731 overlapping :
733 overlapping :
732
734
733 """
735 """
734
736
735 self.__initime = None
737 self.__initime = None
736 self.__lastdatatime = 0
738 self.__lastdatatime = 0
737
739
738 self.__buffer_spc = 0
740 self.__buffer_spc = 0
739 self.__buffer_cspc = 0
741 self.__buffer_cspc = 0
740 self.__buffer_dc = 0
742 self.__buffer_dc = 0
741
743
742 self.__profIndex = 0
744 self.__profIndex = 0
743 self.__dataReady = False
745 self.__dataReady = False
744 self.__byTime = False
746 self.__byTime = False
745
747
746 if n is None and timeInterval is None:
748 if n is None and timeInterval is None:
747 raise ValueError("n or timeInterval should be specified ...")
749 raise ValueError("n or timeInterval should be specified ...")
748
750
749 if n is not None:
751 if n is not None:
750 self.n = int(n)
752 self.n = int(n)
751 else:
753 else:
752
754
753 self.__integrationtime = int(timeInterval)
755 self.__integrationtime = int(timeInterval)
754 self.n = None
756 self.n = None
755 self.__byTime = True
757 self.__byTime = True
756
758
757 def putData(self, data_spc, data_cspc, data_dc):
759 def putData(self, data_spc, data_cspc, data_dc):
758 """
760 """
759 Add a profile to the __buffer_spc and increase in one the __profileIndex
761 Add a profile to the __buffer_spc and increase in one the __profileIndex
760
762
761 """
763 """
762
764
763 self.__buffer_spc += data_spc
765 self.__buffer_spc += data_spc
764
766
765 if data_cspc is None:
767 if data_cspc is None:
766 self.__buffer_cspc = None
768 self.__buffer_cspc = None
767 else:
769 else:
768 self.__buffer_cspc += data_cspc
770 self.__buffer_cspc += data_cspc
769
771
770 if data_dc is None:
772 if data_dc is None:
771 self.__buffer_dc = None
773 self.__buffer_dc = None
772 else:
774 else:
773 self.__buffer_dc += data_dc
775 self.__buffer_dc += data_dc
774
776
775 self.__profIndex += 1
777 self.__profIndex += 1
776
778
777 return
779 return
778
780
779 def pushData(self):
781 def pushData(self):
780 """
782 """
781 Return the sum of the last profiles and the profiles used in the sum.
783 Return the sum of the last profiles and the profiles used in the sum.
782
784
783 Affected:
785 Affected:
784
786
785 self.__profileIndex
787 self.__profileIndex
786
788
787 """
789 """
788
790
789 data_spc = self.__buffer_spc
791 data_spc = self.__buffer_spc
790 data_cspc = self.__buffer_cspc
792 data_cspc = self.__buffer_cspc
791 data_dc = self.__buffer_dc
793 data_dc = self.__buffer_dc
792 n = self.__profIndex
794 n = self.__profIndex
793
795
794 self.__buffer_spc = 0
796 self.__buffer_spc = 0
795 self.__buffer_cspc = 0
797 self.__buffer_cspc = 0
796 self.__buffer_dc = 0
798 self.__buffer_dc = 0
797 self.__profIndex = 0
799 self.__profIndex = 0
798
800
799 return data_spc, data_cspc, data_dc, n
801 return data_spc, data_cspc, data_dc, n
800
802
801 def byProfiles(self, *args):
803 def byProfiles(self, *args):
802
804
803 self.__dataReady = False
805 self.__dataReady = False
804 avgdata_spc = None
806 avgdata_spc = None
805 avgdata_cspc = None
807 avgdata_cspc = None
806 avgdata_dc = None
808 avgdata_dc = None
807
809
808 self.putData(*args)
810 self.putData(*args)
809
811
810 if self.__profIndex == self.n:
812 if self.__profIndex == self.n:
811
813
812 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
814 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
813 self.n = n
815 self.n = n
814 self.__dataReady = True
816 self.__dataReady = True
815
817
816 return avgdata_spc, avgdata_cspc, avgdata_dc
818 return avgdata_spc, avgdata_cspc, avgdata_dc
817
819
818 def byTime(self, datatime, *args):
820 def byTime(self, datatime, *args):
819
821
820 self.__dataReady = False
822 self.__dataReady = False
821 avgdata_spc = None
823 avgdata_spc = None
822 avgdata_cspc = None
824 avgdata_cspc = None
823 avgdata_dc = None
825 avgdata_dc = None
824
826
825 self.putData(*args)
827 self.putData(*args)
826
828
827 if (datatime - self.__initime) >= self.__integrationtime:
829 if (datatime - self.__initime) >= self.__integrationtime:
828 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
830 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
829 self.n = n
831 self.n = n
830 self.__dataReady = True
832 self.__dataReady = True
831
833
832 return avgdata_spc, avgdata_cspc, avgdata_dc
834 return avgdata_spc, avgdata_cspc, avgdata_dc
833
835
834 def integrate(self, datatime, *args):
836 def integrate(self, datatime, *args):
835
837
836 if self.__profIndex == 0:
838 if self.__profIndex == 0:
837 self.__initime = datatime
839 self.__initime = datatime
838
840
839 if self.__byTime:
841 if self.__byTime:
840 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
842 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
841 datatime, *args)
843 datatime, *args)
842 else:
844 else:
843 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
845 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
844
846
845 if not self.__dataReady:
847 if not self.__dataReady:
846 return None, None, None, None
848 return None, None, None, None
847
849
848 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
850 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
849
851
850 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
852 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
851 if n == 1:
853 if n == 1:
852 return dataOut
854 return dataOut
853
855
854 dataOut.flagNoData = True
856 dataOut.flagNoData = True
855
857
856 if not self.isConfig:
858 if not self.isConfig:
857 self.setup(n, timeInterval, overlapping)
859 self.setup(n, timeInterval, overlapping)
858 self.isConfig = True
860 self.isConfig = True
859
861
860 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
862 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
861 dataOut.data_spc,
863 dataOut.data_spc,
862 dataOut.data_cspc,
864 dataOut.data_cspc,
863 dataOut.data_dc)
865 dataOut.data_dc)
864
866
865 if self.__dataReady:
867 if self.__dataReady:
866
868
867 dataOut.data_spc = avgdata_spc
869 dataOut.data_spc = avgdata_spc
868 dataOut.data_cspc = avgdata_cspc
870 dataOut.data_cspc = avgdata_cspc
869 dataOut.data_dc = avgdata_dc
871 dataOut.data_dc = avgdata_dc
870 dataOut.nIncohInt *= self.n
872 dataOut.nIncohInt *= self.n
871 dataOut.utctime = avgdatatime
873 dataOut.utctime = avgdatatime
872 dataOut.flagNoData = False
874 dataOut.flagNoData = False
873
875
874 return dataOut No newline at end of file
876 return dataOut
General Comments 0
You need to be logged in to leave comments. Login now