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