##// END OF EJS Templates
schain
RSS
Remove numpy.complex
Juan C. Espinoza -
r1595:fed0fad1169f
parent child
Show More
Side by Side Unified
Context file:
r1595:fed0fad1169f Loading diff...:

The requested changes are too big and content was truncated. Show full diff

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