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