##// END OF EJS Templates
schain
RSS
Update metod deFlip y dopplerFlip
Alexander Valdez -
r1667:7984c78c901f
parent child
Show More
Side by Side Unified
Context file:
r1667:7984c78c901f Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,898 +1,898
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 # repositorio
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
158
159 if nVoltProfiles == nProfiles:
159 if nVoltProfiles == nProfiles:
160 self.buffer = self.dataIn.data.copy()
160 self.buffer = self.dataIn.data.copy()
161 self.profIndex = nVoltProfiles
161 self.profIndex = nVoltProfiles
162
162
163 elif nVoltProfiles < nProfiles:
163 elif nVoltProfiles < nProfiles:
164
164
165 if self.profIndex == 0:
165 if self.profIndex == 0:
166 self.id_min = 0
166 self.id_min = 0
167 self.id_max = nVoltProfiles
167 self.id_max = nVoltProfiles
168
168
169 self.buffer[:, self.id_min:self.id_max,
169 self.buffer[:, self.id_min:self.id_max,
170 :] = self.dataIn.data
170 :] = self.dataIn.data
171 self.profIndex += nVoltProfiles
171 self.profIndex += nVoltProfiles
172 self.id_min += nVoltProfiles
172 self.id_min += nVoltProfiles
173 self.id_max += nVoltProfiles
173 self.id_max += nVoltProfiles
174 else:
174 else:
175 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
175 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
176 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
176 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
177 self.dataOut.flagNoData = True
177 self.dataOut.flagNoData = True
178 else:
178 else:
179 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
179 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
180 self.profIndex += 1
180 self.profIndex += 1
181
181
182 if self.firstdatatime == None:
182 if self.firstdatatime == None:
183 self.firstdatatime = self.dataIn.utctime
183 self.firstdatatime = self.dataIn.utctime
184
184
185 if self.profIndex == nProfiles:
185 if self.profIndex == nProfiles:
186 self.__updateSpecFromVoltage()
186 self.__updateSpecFromVoltage()
187 if pairsList == None:
187 if pairsList == None:
188 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
188 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
189 else:
189 else:
190 self.dataOut.pairsList = pairsList
190 self.dataOut.pairsList = pairsList
191 self.__getFft()
191 self.__getFft()
192 self.dataOut.flagNoData = False
192 self.dataOut.flagNoData = False
193 self.firstdatatime = None
193 self.firstdatatime = None
194 self.profIndex = 0
194 self.profIndex = 0
195 else:
195 else:
196 raise ValueError("The type of input object '%s' is not valid".format(
196 raise ValueError("The type of input object '%s' is not valid".format(
197 self.dataIn.type))
197 self.dataIn.type))
198
198
199 def __selectPairs(self, pairsList):
199 def __selectPairs(self, pairsList):
200
200
201 if not pairsList:
201 if not pairsList:
202 return
202 return
203
203
204 pairs = []
204 pairs = []
205 pairsIndex = []
205 pairsIndex = []
206
206
207 for pair in pairsList:
207 for pair in pairsList:
208 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
208 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
209 continue
209 continue
210 pairs.append(pair)
210 pairs.append(pair)
211 pairsIndex.append(pairs.index(pair))
211 pairsIndex.append(pairs.index(pair))
212
212
213 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
213 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
214 self.dataOut.pairsList = pairs
214 self.dataOut.pairsList = pairs
215
215
216 return
216 return
217
217
218 def selectFFTs(self, minFFT, maxFFT ):
218 def selectFFTs(self, minFFT, maxFFT ):
219 """
219 """
220 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
220 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
221 minFFT<= FFT <= maxFFT
221 minFFT<= FFT <= maxFFT
222 """
222 """
223
223
224 if (minFFT > maxFFT):
224 if (minFFT > maxFFT):
225 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
225 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
226
226
227 if (minFFT < self.dataOut.getFreqRange()[0]):
227 if (minFFT < self.dataOut.getFreqRange()[0]):
228 minFFT = self.dataOut.getFreqRange()[0]
228 minFFT = self.dataOut.getFreqRange()[0]
229
229
230 if (maxFFT > self.dataOut.getFreqRange()[-1]):
230 if (maxFFT > self.dataOut.getFreqRange()[-1]):
231 maxFFT = self.dataOut.getFreqRange()[-1]
231 maxFFT = self.dataOut.getFreqRange()[-1]
232
232
233 minIndex = 0
233 minIndex = 0
234 maxIndex = 0
234 maxIndex = 0
235 FFTs = self.dataOut.getFreqRange()
235 FFTs = self.dataOut.getFreqRange()
236
236
237 inda = numpy.where(FFTs >= minFFT)
237 inda = numpy.where(FFTs >= minFFT)
238 indb = numpy.where(FFTs <= maxFFT)
238 indb = numpy.where(FFTs <= maxFFT)
239
239
240 try:
240 try:
241 minIndex = inda[0][0]
241 minIndex = inda[0][0]
242 except:
242 except:
243 minIndex = 0
243 minIndex = 0
244
244
245 try:
245 try:
246 maxIndex = indb[0][-1]
246 maxIndex = indb[0][-1]
247 except:
247 except:
248 maxIndex = len(FFTs)
248 maxIndex = len(FFTs)
249
249
250 self.selectFFTsByIndex(minIndex, maxIndex)
250 self.selectFFTsByIndex(minIndex, maxIndex)
251
251
252 return 1
252 return 1
253
253
254 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
254 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
255 newheis = numpy.where(
255 newheis = numpy.where(
256 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
256 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
257
257
258 if hei_ref != None:
258 if hei_ref != None:
259 newheis = numpy.where(self.dataOut.heightList > hei_ref)
259 newheis = numpy.where(self.dataOut.heightList > hei_ref)
260
260
261 minIndex = min(newheis[0])
261 minIndex = min(newheis[0])
262 maxIndex = max(newheis[0])
262 maxIndex = max(newheis[0])
263 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
263 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
264 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
264 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
265
265
266 # determina indices
266 # determina indices
267 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
267 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
268 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
268 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
269 avg_dB = 10 * \
269 avg_dB = 10 * \
270 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
270 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
271 beacon_dB = numpy.sort(avg_dB)[-nheis:]
271 beacon_dB = numpy.sort(avg_dB)[-nheis:]
272 beacon_heiIndexList = []
272 beacon_heiIndexList = []
273 for val in avg_dB.tolist():
273 for val in avg_dB.tolist():
274 if val >= beacon_dB[0]:
274 if val >= beacon_dB[0]:
275 beacon_heiIndexList.append(avg_dB.tolist().index(val))
275 beacon_heiIndexList.append(avg_dB.tolist().index(val))
276
276
277 #data_spc = data_spc[:,:,beacon_heiIndexList]
277 #data_spc = data_spc[:,:,beacon_heiIndexList]
278 data_cspc = None
278 data_cspc = None
279 if self.dataOut.data_cspc is not None:
279 if self.dataOut.data_cspc is not None:
280 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
280 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
281 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
281 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
282
282
283 data_dc = None
283 data_dc = None
284 if self.dataOut.data_dc is not None:
284 if self.dataOut.data_dc is not None:
285 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
285 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
286 #data_dc = data_dc[:,beacon_heiIndexList]
286 #data_dc = data_dc[:,beacon_heiIndexList]
287
287
288 self.dataOut.data_spc = data_spc
288 self.dataOut.data_spc = data_spc
289 self.dataOut.data_cspc = data_cspc
289 self.dataOut.data_cspc = data_cspc
290 self.dataOut.data_dc = data_dc
290 self.dataOut.data_dc = data_dc
291 self.dataOut.heightList = heightList
291 self.dataOut.heightList = heightList
292 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
292 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
293
293
294 return 1
294 return 1
295
295
296 def selectFFTsByIndex(self, minIndex, maxIndex):
296 def selectFFTsByIndex(self, minIndex, maxIndex):
297 """
297 """
298
298
299 """
299 """
300
300
301 if (minIndex < 0) or (minIndex > maxIndex):
301 if (minIndex < 0) or (minIndex > maxIndex):
302 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
302 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
303
303
304 if (maxIndex >= self.dataOut.nProfiles):
304 if (maxIndex >= self.dataOut.nProfiles):
305 maxIndex = self.dataOut.nProfiles-1
305 maxIndex = self.dataOut.nProfiles-1
306
306
307 #Spectra
307 #Spectra
308 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
308 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
309
309
310 data_cspc = None
310 data_cspc = None
311 if self.dataOut.data_cspc is not None:
311 if self.dataOut.data_cspc is not None:
312 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
312 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
313
313
314 data_dc = None
314 data_dc = None
315 if self.dataOut.data_dc is not None:
315 if self.dataOut.data_dc is not None:
316 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
316 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
317
317
318 self.dataOut.data_spc = data_spc
318 self.dataOut.data_spc = data_spc
319 self.dataOut.data_cspc = data_cspc
319 self.dataOut.data_cspc = data_cspc
320 self.dataOut.data_dc = data_dc
320 self.dataOut.data_dc = data_dc
321
321
322 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
322 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
323 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
323 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
324 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
324 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
325
325
326 return 1
326 return 1
327
327
328 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
328 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
329 # validacion de rango
329 # validacion de rango
330 if minHei == None:
330 if minHei == None:
331 minHei = self.dataOut.heightList[0]
331 minHei = self.dataOut.heightList[0]
332
332
333 if maxHei == None:
333 if maxHei == None:
334 maxHei = self.dataOut.heightList[-1]
334 maxHei = self.dataOut.heightList[-1]
335
335
336 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
336 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
337 print('minHei: %.2f is out of the heights range' % (minHei))
337 print('minHei: %.2f is out of the heights range' % (minHei))
338 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
338 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
339 minHei = self.dataOut.heightList[0]
339 minHei = self.dataOut.heightList[0]
340
340
341 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
341 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
342 print('maxHei: %.2f is out of the heights range' % (maxHei))
342 print('maxHei: %.2f is out of the heights range' % (maxHei))
343 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
343 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
344 maxHei = self.dataOut.heightList[-1]
344 maxHei = self.dataOut.heightList[-1]
345
345
346 # validacion de velocidades
346 # validacion de velocidades
347 velrange = self.dataOut.getVelRange(1)
347 velrange = self.dataOut.getVelRange(1)
348
348
349 if minVel == None:
349 if minVel == None:
350 minVel = velrange[0]
350 minVel = velrange[0]
351
351
352 if maxVel == None:
352 if maxVel == None:
353 maxVel = velrange[-1]
353 maxVel = velrange[-1]
354
354
355 if (minVel < velrange[0]) or (minVel > maxVel):
355 if (minVel < velrange[0]) or (minVel > maxVel):
356 print('minVel: %.2f is out of the velocity range' % (minVel))
356 print('minVel: %.2f is out of the velocity range' % (minVel))
357 print('minVel is setting to %.2f' % (velrange[0]))
357 print('minVel is setting to %.2f' % (velrange[0]))
358 minVel = velrange[0]
358 minVel = velrange[0]
359
359
360 if (maxVel > velrange[-1]) or (maxVel < minVel):
360 if (maxVel > velrange[-1]) or (maxVel < minVel):
361 print('maxVel: %.2f is out of the velocity range' % (maxVel))
361 print('maxVel: %.2f is out of the velocity range' % (maxVel))
362 print('maxVel is setting to %.2f' % (velrange[-1]))
362 print('maxVel is setting to %.2f' % (velrange[-1]))
363 maxVel = velrange[-1]
363 maxVel = velrange[-1]
364
364
365 # seleccion de indices para rango
365 # seleccion de indices para rango
366 minIndex = 0
366 minIndex = 0
367 maxIndex = 0
367 maxIndex = 0
368 heights = self.dataOut.heightList
368 heights = self.dataOut.heightList
369
369
370 inda = numpy.where(heights >= minHei)
370 inda = numpy.where(heights >= minHei)
371 indb = numpy.where(heights <= maxHei)
371 indb = numpy.where(heights <= maxHei)
372
372
373 try:
373 try:
374 minIndex = inda[0][0]
374 minIndex = inda[0][0]
375 except:
375 except:
376 minIndex = 0
376 minIndex = 0
377
377
378 try:
378 try:
379 maxIndex = indb[0][-1]
379 maxIndex = indb[0][-1]
380 except:
380 except:
381 maxIndex = len(heights)
381 maxIndex = len(heights)
382
382
383 if (minIndex < 0) or (minIndex > maxIndex):
383 if (minIndex < 0) or (minIndex > maxIndex):
384 raise ValueError("some value in (%d,%d) is not valid" % (
384 raise ValueError("some value in (%d,%d) is not valid" % (
385 minIndex, maxIndex))
385 minIndex, maxIndex))
386
386
387 if (maxIndex >= self.dataOut.nHeights):
387 if (maxIndex >= self.dataOut.nHeights):
388 maxIndex = self.dataOut.nHeights - 1
388 maxIndex = self.dataOut.nHeights - 1
389
389
390 # seleccion de indices para velocidades
390 # seleccion de indices para velocidades
391 indminvel = numpy.where(velrange >= minVel)
391 indminvel = numpy.where(velrange >= minVel)
392 indmaxvel = numpy.where(velrange <= maxVel)
392 indmaxvel = numpy.where(velrange <= maxVel)
393 try:
393 try:
394 minIndexVel = indminvel[0][0]
394 minIndexVel = indminvel[0][0]
395 except:
395 except:
396 minIndexVel = 0
396 minIndexVel = 0
397
397
398 try:
398 try:
399 maxIndexVel = indmaxvel[0][-1]
399 maxIndexVel = indmaxvel[0][-1]
400 except:
400 except:
401 maxIndexVel = len(velrange)
401 maxIndexVel = len(velrange)
402
402
403 # seleccion del espectro
403 # seleccion del espectro
404 data_spc = self.dataOut.data_spc[:,
404 data_spc = self.dataOut.data_spc[:,
405 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
405 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
406 # estimacion de ruido
406 # estimacion de ruido
407 noise = numpy.zeros(self.dataOut.nChannels)
407 noise = numpy.zeros(self.dataOut.nChannels)
408
408
409 for channel in range(self.dataOut.nChannels):
409 for channel in range(self.dataOut.nChannels):
410 daux = data_spc[channel, :, :]
410 daux = data_spc[channel, :, :]
411 sortdata = numpy.sort(daux, axis=None)
411 sortdata = numpy.sort(daux, axis=None)
412 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
412 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
413
413
414 self.dataOut.noise_estimation = noise.copy()
414 self.dataOut.noise_estimation = noise.copy()
415
415
416 return 1
416 return 1
417
417
418 class removeDC(Operation):
418 class removeDC(Operation):
419
419
420 def run(self, dataOut, mode=2):
420 def run(self, dataOut, mode=2):
421 self.dataOut = dataOut
421 self.dataOut = dataOut
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 num_chan = jspectra.shape[0]
425 num_chan = jspectra.shape[0]
426 num_hei = jspectra.shape[2]
426 num_hei = jspectra.shape[2]
427
427
428 if jcspectra is not None:
428 if jcspectra is not None:
429 jcspectraExist = True
429 jcspectraExist = True
430 num_pairs = jcspectra.shape[0]
430 num_pairs = jcspectra.shape[0]
431 else:
431 else:
432 jcspectraExist = False
432 jcspectraExist = False
433
433
434 freq_dc = int(jspectra.shape[1] / 2)
434 freq_dc = int(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 ind_vel = ind_vel.astype(int)
436 ind_vel = ind_vel.astype(int)
437
437
438 if ind_vel[0] < 0:
438 if ind_vel[0] < 0:
439 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
439 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
440
440
441 if mode == 1:
441 if mode == 1:
442 jspectra[:, freq_dc, :] = (
442 jspectra[:, freq_dc, :] = (
443 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
443 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
444
444
445 if jcspectraExist:
445 if jcspectraExist:
446 jcspectra[:, freq_dc, :] = (
446 jcspectra[:, freq_dc, :] = (
447 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
447 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()] = (
468 jspectra[ich, freq_dc, junkid.nonzero()] = (
469 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
469 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
470
470
471 if jcspectraExist:
471 if jcspectraExist:
472 for ip in range(num_pairs):
472 for ip in range(num_pairs):
473 yy = jcspectra[ip, ind_vel, :]
473 yy = jcspectra[ip, ind_vel, :]
474 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
474 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
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 self.dataOut
479 return self.dataOut
480
480
481 class removeInterference(Operation):
481 class removeInterference(Operation):
482
482
483 def removeInterference2(self):
483 def removeInterference2(self):
484
484
485 cspc = self.dataOut.data_cspc
485 cspc = self.dataOut.data_cspc
486 spc = self.dataOut.data_spc
486 spc = self.dataOut.data_spc
487 Heights = numpy.arange(cspc.shape[2])
487 Heights = numpy.arange(cspc.shape[2])
488 realCspc = numpy.abs(cspc)
488 realCspc = numpy.abs(cspc)
489
489
490 for i in range(cspc.shape[0]):
490 for i in range(cspc.shape[0]):
491 LinePower= numpy.sum(realCspc[i], axis=0)
491 LinePower= numpy.sum(realCspc[i], axis=0)
492 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
492 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
493 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
493 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
494 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
494 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
495 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
495 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
496 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
496 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
497
497
498
498
499 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
499 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
500 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
500 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
501 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
501 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
502 cspc[i,InterferenceRange,:] = numpy.NaN
502 cspc[i,InterferenceRange,:] = numpy.NaN
503
503
504 self.dataOut.data_cspc = cspc
504 self.dataOut.data_cspc = cspc
505
505
506 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
506 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
507
507
508 jspectra = self.dataOut.data_spc
508 jspectra = self.dataOut.data_spc
509 jcspectra = self.dataOut.data_cspc
509 jcspectra = self.dataOut.data_cspc
510 jnoise = self.dataOut.getNoise()
510 jnoise = self.dataOut.getNoise()
511 num_incoh = self.dataOut.nIncohInt
511 num_incoh = self.dataOut.nIncohInt
512
512
513 num_channel = jspectra.shape[0]
513 num_channel = jspectra.shape[0]
514 num_prof = jspectra.shape[1]
514 num_prof = jspectra.shape[1]
515 num_hei = jspectra.shape[2]
515 num_hei = jspectra.shape[2]
516
516
517 # hei_interf
517 # hei_interf
518 if hei_interf is None:
518 if hei_interf is None:
519 count_hei = int(num_hei / 2)
519 count_hei = int(num_hei / 2)
520 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
520 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
521 hei_interf = numpy.asarray(hei_interf)[0]
521 hei_interf = numpy.asarray(hei_interf)[0]
522 # nhei_interf
522 # nhei_interf
523 if (nhei_interf == None):
523 if (nhei_interf == None):
524 nhei_interf = 5
524 nhei_interf = 5
525 if (nhei_interf < 1):
525 if (nhei_interf < 1):
526 nhei_interf = 1
526 nhei_interf = 1
527 if (nhei_interf > count_hei):
527 if (nhei_interf > count_hei):
528 nhei_interf = count_hei
528 nhei_interf = count_hei
529 if (offhei_interf == None):
529 if (offhei_interf == None):
530 offhei_interf = 0
530 offhei_interf = 0
531
531
532 ind_hei = list(range(num_hei))
532 ind_hei = list(range(num_hei))
533 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
533 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
534 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
534 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
535 mask_prof = numpy.asarray(list(range(num_prof)))
535 mask_prof = numpy.asarray(list(range(num_prof)))
536 num_mask_prof = mask_prof.size
536 num_mask_prof = mask_prof.size
537 comp_mask_prof = [0, num_prof / 2]
537 comp_mask_prof = [0, num_prof / 2]
538
538
539 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
539 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
540 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
540 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
541 jnoise = numpy.nan
541 jnoise = numpy.nan
542 noise_exist = jnoise[0] < numpy.Inf
542 noise_exist = jnoise[0] < numpy.Inf
543
543
544 # Subrutina de Remocion de la Interferencia
544 # Subrutina de Remocion de la Interferencia
545 for ich in range(num_channel):
545 for ich in range(num_channel):
546 # Se ordena los espectros segun su potencia (menor a mayor)
546 # Se ordena los espectros segun su potencia (menor a mayor)
547 power = jspectra[ich, mask_prof, :]
547 power = jspectra[ich, mask_prof, :]
548 power = power[:, hei_interf]
548 power = power[:, hei_interf]
549 power = power.sum(axis=0)
549 power = power.sum(axis=0)
550 psort = power.ravel().argsort()
550 psort = power.ravel().argsort()
551
551
552 # Se estima la interferencia promedio en los Espectros de Potencia empleando
552 # Se estima la interferencia promedio en los Espectros de Potencia empleando
553 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
553 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
554 offhei_interf, nhei_interf + offhei_interf))]]]
554 offhei_interf, nhei_interf + offhei_interf))]]]
555
555
556 if noise_exist:
556 if noise_exist:
557 # tmp_noise = jnoise[ich] / num_prof
557 # tmp_noise = jnoise[ich] / num_prof
558 tmp_noise = jnoise[ich]
558 tmp_noise = jnoise[ich]
559 junkspc_interf = junkspc_interf - tmp_noise
559 junkspc_interf = junkspc_interf - tmp_noise
560 #junkspc_interf[:,comp_mask_prof] = 0
560 #junkspc_interf[:,comp_mask_prof] = 0
561
561
562 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
562 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
563 jspc_interf = jspc_interf.transpose()
563 jspc_interf = jspc_interf.transpose()
564 # Calculando el espectro de interferencia promedio
564 # Calculando el espectro de interferencia promedio
565 noiseid = numpy.where(
565 noiseid = numpy.where(
566 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
566 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
567 noiseid = noiseid[0]
567 noiseid = noiseid[0]
568 cnoiseid = noiseid.size
568 cnoiseid = noiseid.size
569 interfid = numpy.where(
569 interfid = numpy.where(
570 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
570 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
571 interfid = interfid[0]
571 interfid = interfid[0]
572 cinterfid = interfid.size
572 cinterfid = interfid.size
573
573
574 if (cnoiseid > 0):
574 if (cnoiseid > 0):
575 jspc_interf[noiseid] = 0
575 jspc_interf[noiseid] = 0
576
576
577 # Expandiendo los perfiles a limpiar
577 # Expandiendo los perfiles a limpiar
578 if (cinterfid > 0):
578 if (cinterfid > 0):
579 new_interfid = (
579 new_interfid = (
580 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
580 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
581 new_interfid = numpy.asarray(new_interfid)
581 new_interfid = numpy.asarray(new_interfid)
582 new_interfid = {x for x in new_interfid}
582 new_interfid = {x for x in new_interfid}
583 new_interfid = numpy.array(list(new_interfid))
583 new_interfid = numpy.array(list(new_interfid))
584 new_cinterfid = new_interfid.size
584 new_cinterfid = new_interfid.size
585 else:
585 else:
586 new_cinterfid = 0
586 new_cinterfid = 0
587
587
588 for ip in range(new_cinterfid):
588 for ip in range(new_cinterfid):
589 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
589 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
590 jspc_interf[new_interfid[ip]
590 jspc_interf[new_interfid[ip]
591 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
591 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
592
592
593 jspectra[ich, :, ind_hei] = jspectra[ich, :,
593 jspectra[ich, :, ind_hei] = jspectra[ich, :,
594 ind_hei] - jspc_interf # Corregir indices
594 ind_hei] - jspc_interf # Corregir indices
595
595
596 # Removiendo la interferencia del punto de mayor interferencia
596 # Removiendo la interferencia del punto de mayor interferencia
597 ListAux = jspc_interf[mask_prof].tolist()
597 ListAux = jspc_interf[mask_prof].tolist()
598 maxid = ListAux.index(max(ListAux))
598 maxid = ListAux.index(max(ListAux))
599
599
600 if cinterfid > 0:
600 if cinterfid > 0:
601 for ip in range(cinterfid * (interf == 2) - 1):
601 for ip in range(cinterfid * (interf == 2) - 1):
602 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
602 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
603 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
603 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
604 cind = len(ind)
604 cind = len(ind)
605
605
606 if (cind > 0):
606 if (cind > 0):
607 jspectra[ich, interfid[ip], ind] = tmp_noise * \
607 jspectra[ich, interfid[ip], ind] = tmp_noise * \
608 (1 + (numpy.random.uniform(cind) - 0.5) /
608 (1 + (numpy.random.uniform(cind) - 0.5) /
609 numpy.sqrt(num_incoh))
609 numpy.sqrt(num_incoh))
610
610
611 ind = numpy.array([-2, -1, 1, 2])
611 ind = numpy.array([-2, -1, 1, 2])
612 xx = numpy.zeros([4, 4])
612 xx = numpy.zeros([4, 4])
613
613
614 for id1 in range(4):
614 for id1 in range(4):
615 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
615 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
616
616
617 xx_inv = numpy.linalg.inv(xx)
617 xx_inv = numpy.linalg.inv(xx)
618 xx = xx_inv[:, 0]
618 xx = xx_inv[:, 0]
619 ind = (ind + maxid + num_mask_prof) % num_mask_prof
619 ind = (ind + maxid + num_mask_prof) % num_mask_prof
620 yy = jspectra[ich, mask_prof[ind], :]
620 yy = jspectra[ich, mask_prof[ind], :]
621 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
621 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
622 yy.transpose(), xx)
622 yy.transpose(), xx)
623
623
624 indAux = (jspectra[ich, :, :] < tmp_noise *
624 indAux = (jspectra[ich, :, :] < tmp_noise *
625 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
625 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
626 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
626 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
627 (1 - 1 / numpy.sqrt(num_incoh))
627 (1 - 1 / numpy.sqrt(num_incoh))
628
628
629 # Remocion de Interferencia en el Cross Spectra
629 # Remocion de Interferencia en el Cross Spectra
630 if jcspectra is None:
630 if jcspectra is None:
631 return jspectra, jcspectra
631 return jspectra, jcspectra
632 num_pairs = int(jcspectra.size / (num_prof * num_hei))
632 num_pairs = int(jcspectra.size / (num_prof * num_hei))
633 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
633 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
634
634
635 for ip in range(num_pairs):
635 for ip in range(num_pairs):
636
636
637 #-------------------------------------------
637 #-------------------------------------------
638
638
639 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
639 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
640 cspower = cspower[:, hei_interf]
640 cspower = cspower[:, hei_interf]
641 cspower = cspower.sum(axis=0)
641 cspower = cspower.sum(axis=0)
642
642
643 cspsort = cspower.ravel().argsort()
643 cspsort = cspower.ravel().argsort()
644 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
644 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
645 offhei_interf, nhei_interf + offhei_interf))]]]
645 offhei_interf, nhei_interf + offhei_interf))]]]
646 junkcspc_interf = junkcspc_interf.transpose()
646 junkcspc_interf = junkcspc_interf.transpose()
647 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
647 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
648
648
649 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
649 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
650
650
651 median_real = int(numpy.median(numpy.real(
651 median_real = int(numpy.median(numpy.real(
652 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
652 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
653 median_imag = int(numpy.median(numpy.imag(
653 median_imag = int(numpy.median(numpy.imag(
654 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
654 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
655 comp_mask_prof = [int(e) for e in comp_mask_prof]
655 comp_mask_prof = [int(e) for e in comp_mask_prof]
656 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
656 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
657 median_real, median_imag)
657 median_real, median_imag)
658
658
659 for iprof in range(num_prof):
659 for iprof in range(num_prof):
660 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
660 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
661 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
661 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
662
662
663 # Removiendo la Interferencia
663 # Removiendo la Interferencia
664 jcspectra[ip, :, ind_hei] = jcspectra[ip,
664 jcspectra[ip, :, ind_hei] = jcspectra[ip,
665 :, ind_hei] - jcspc_interf
665 :, ind_hei] - jcspc_interf
666
666
667 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
667 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
668 maxid = ListAux.index(max(ListAux))
668 maxid = ListAux.index(max(ListAux))
669
669
670 ind = numpy.array([-2, -1, 1, 2])
670 ind = numpy.array([-2, -1, 1, 2])
671 xx = numpy.zeros([4, 4])
671 xx = numpy.zeros([4, 4])
672
672
673 for id1 in range(4):
673 for id1 in range(4):
674 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
674 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
675
675
676 xx_inv = numpy.linalg.inv(xx)
676 xx_inv = numpy.linalg.inv(xx)
677 xx = xx_inv[:, 0]
677 xx = xx_inv[:, 0]
678
678
679 ind = (ind + maxid + num_mask_prof) % num_mask_prof
679 ind = (ind + maxid + num_mask_prof) % num_mask_prof
680 yy = jcspectra[ip, mask_prof[ind], :]
680 yy = jcspectra[ip, mask_prof[ind], :]
681 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
681 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
682
682
683 # Guardar Resultados
683 # Guardar Resultados
684 self.dataOut.data_spc = jspectra
684 self.dataOut.data_spc = jspectra
685 self.dataOut.data_cspc = jcspectra
685 self.dataOut.data_cspc = jcspectra
686
686
687 return 1
687 return 1
688
688
689 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
689 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
690
690
691 self.dataOut = dataOut
691 self.dataOut = dataOut
692
692
693 if mode == 1:
693 if mode == 1:
694 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
694 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
695 elif mode == 2:
695 elif mode == 2:
696 self.removeInterference2()
696 self.removeInterference2()
697
697
698 return self.dataOut
698 return self.dataOut
699
699
700
700
701 class IncohInt(Operation):
701 class IncohInt(Operation):
702
702
703 __profIndex = 0
703 __profIndex = 0
704 __withOverapping = False
704 __withOverapping = False
705
705
706 __byTime = False
706 __byTime = False
707 __initime = None
707 __initime = None
708 __lastdatatime = None
708 __lastdatatime = None
709 __integrationtime = None
709 __integrationtime = None
710
710
711 __buffer_spc = None
711 __buffer_spc = None
712 __buffer_cspc = None
712 __buffer_cspc = None
713 __buffer_dc = None
713 __buffer_dc = None
714
714
715 __dataReady = False
715 __dataReady = False
716
716
717 __timeInterval = None
717 __timeInterval = None
718
718
719 n = None
719 n = None
720
720
721 def __init__(self):
721 def __init__(self):
722
722
723 Operation.__init__(self)
723 Operation.__init__(self)
724
724
725 def setup(self, n=None, timeInterval=None, overlapping=False):
725 def setup(self, n=None, timeInterval=None, overlapping=False):
726 """
726 """
727 Set the parameters of the integration class.
727 Set the parameters of the integration class.
728
728
729 Inputs:
729 Inputs:
730
730
731 n : Number of coherent integrations
731 n : Number of coherent integrations
732 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
732 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
733 overlapping :
733 overlapping :
734
734
735 """
735 """
736
736
737 self.__initime = None
737 self.__initime = None
738 self.__lastdatatime = 0
738 self.__lastdatatime = 0
739
739
740 self.__buffer_spc = 0
740 self.__buffer_spc = 0
741 self.__buffer_cspc = 0
741 self.__buffer_cspc = 0
742 self.__buffer_dc = 0
742 self.__buffer_dc = 0
743
743
744 self.__profIndex = 0
744 self.__profIndex = 0
745 self.__dataReady = False
745 self.__dataReady = False
746 self.__byTime = False
746 self.__byTime = False
747
747
748 if n is None and timeInterval is None:
748 if n is None and timeInterval is None:
749 raise ValueError("n or timeInterval should be specified ...")
749 raise ValueError("n or timeInterval should be specified ...")
750
750
751 if n is not None:
751 if n is not None:
752 self.n = int(n)
752 self.n = int(n)
753 else:
753 else:
754
754
755 self.__integrationtime = int(timeInterval)
755 self.__integrationtime = int(timeInterval)
756 self.n = None
756 self.n = None
757 self.__byTime = True
757 self.__byTime = True
758
758
759 def putData(self, data_spc, data_cspc, data_dc):
759 def putData(self, data_spc, data_cspc, data_dc):
760 """
760 """
761 Add a profile to the __buffer_spc and increase in one the __profileIndex
761 Add a profile to the __buffer_spc and increase in one the __profileIndex
762
762
763 """
763 """
764
764
765 self.__buffer_spc += data_spc
765 self.__buffer_spc += data_spc
766
766
767 if data_cspc is None:
767 if data_cspc is None:
768 self.__buffer_cspc = None
768 self.__buffer_cspc = None
769 else:
769 else:
770 self.__buffer_cspc += data_cspc
770 self.__buffer_cspc += data_cspc
771
771
772 if data_dc is None:
772 if data_dc is None:
773 self.__buffer_dc = None
773 self.__buffer_dc = None
774 else:
774 else:
775 self.__buffer_dc += data_dc
775 self.__buffer_dc += data_dc
776
776
777 self.__profIndex += 1
777 self.__profIndex += 1
778
778
779 return
779 return
780
780
781 def pushData(self):
781 def pushData(self):
782 """
782 """
783 Return the sum of the last profiles and the profiles used in the sum.
783 Return the sum of the last profiles and the profiles used in the sum.
784
784
785 Affected:
785 Affected:
786
786
787 self.__profileIndex
787 self.__profileIndex
788
788
789 """
789 """
790
790
791 data_spc = self.__buffer_spc
791 data_spc = self.__buffer_spc
792 data_cspc = self.__buffer_cspc
792 data_cspc = self.__buffer_cspc
793 data_dc = self.__buffer_dc
793 data_dc = self.__buffer_dc
794 n = self.__profIndex
794 n = self.__profIndex
795
795
796 self.__buffer_spc = 0
796 self.__buffer_spc = 0
797 self.__buffer_cspc = 0
797 self.__buffer_cspc = 0
798 self.__buffer_dc = 0
798 self.__buffer_dc = 0
799 self.__profIndex = 0
799 self.__profIndex = 0
800
800
801 return data_spc, data_cspc, data_dc, n
801 return data_spc, data_cspc, data_dc, n
802
802
803 def byProfiles(self, *args):
803 def byProfiles(self, *args):
804
804
805 self.__dataReady = False
805 self.__dataReady = False
806 avgdata_spc = None
806 avgdata_spc = None
807 avgdata_cspc = None
807 avgdata_cspc = None
808 avgdata_dc = None
808 avgdata_dc = None
809
809
810 self.putData(*args)
810 self.putData(*args)
811
811
812 if self.__profIndex == self.n:
812 if self.__profIndex == self.n:
813
813
814 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
814 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
815 self.n = n
815 self.n = n
816 self.__dataReady = True
816 self.__dataReady = True
817
817
818 return avgdata_spc, avgdata_cspc, avgdata_dc
818 return avgdata_spc, avgdata_cspc, avgdata_dc
819
819
820 def byTime(self, datatime, *args):
820 def byTime(self, datatime, *args):
821
821
822 self.__dataReady = False
822 self.__dataReady = False
823 avgdata_spc = None
823 avgdata_spc = None
824 avgdata_cspc = None
824 avgdata_cspc = None
825 avgdata_dc = None
825 avgdata_dc = None
826
826
827 self.putData(*args)
827 self.putData(*args)
828
828
829 if (datatime - self.__initime) >= self.__integrationtime:
829 if (datatime - self.__initime) >= self.__integrationtime:
830 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
830 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
831 self.n = n
831 self.n = n
832 self.__dataReady = True
832 self.__dataReady = True
833
833
834 return avgdata_spc, avgdata_cspc, avgdata_dc
834 return avgdata_spc, avgdata_cspc, avgdata_dc
835
835
836 def integrate(self, datatime, *args):
836 def integrate(self, datatime, *args):
837
837
838 if self.__profIndex == 0:
838 if self.__profIndex == 0:
839 self.__initime = datatime
839 self.__initime = datatime
840
840
841 if self.__byTime:
841 if self.__byTime:
842 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
842 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
843 datatime, *args)
843 datatime, *args)
844 else:
844 else:
845 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
845 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
846
846
847 if not self.__dataReady:
847 if not self.__dataReady:
848 return None, None, None, None
848 return None, None, None, None
849
849
850 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
850 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
851
851
852 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
852 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
853 if n == 1:
853 if n == 1:
854 return dataOut
854 return dataOut
855
855
856 dataOut.flagNoData = True
856 dataOut.flagNoData = True
857
857
858 if not self.isConfig:
858 if not self.isConfig:
859 self.setup(n, timeInterval, overlapping)
859 self.setup(n, timeInterval, overlapping)
860 self.isConfig = True
860 self.isConfig = True
861
861
862 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
862 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
863 dataOut.data_spc,
863 dataOut.data_spc,
864 dataOut.data_cspc,
864 dataOut.data_cspc,
865 dataOut.data_dc)
865 dataOut.data_dc)
866
866
867 if self.__dataReady:
867 if self.__dataReady:
868
868
869 dataOut.data_spc = avgdata_spc
869 dataOut.data_spc = avgdata_spc
870 dataOut.data_cspc = avgdata_cspc
870 dataOut.data_cspc = avgdata_cspc
871 dataOut.data_dc = avgdata_dc
871 dataOut.data_dc = avgdata_dc
872 dataOut.nIncohInt *= self.n
872 dataOut.nIncohInt *= self.n
873 dataOut.utctime = avgdatatime
873 dataOut.utctime = avgdatatime
874 dataOut.flagNoData = False
874 dataOut.flagNoData = False
875
875
876 return dataOut
876 return dataOut
877
877
878 class dopplerFlip(Operation):
878 class dopplerFlip(Operation):
879
879
880 def run(self, dataOut):
880 def run(self, dataOut):
881 # arreglo 1: (num_chan, num_profiles, num_heights)
881 # arreglo 1: (num_chan, num_profiles, num_heights)
882 self.dataOut = dataOut
882 self.dataOut = dataOut
883 # JULIA-oblicua, indice 2
883 # JULIA-oblicua, indice 2
884 # arreglo 2: (num_profiles, num_heights)
884 # arreglo 2: (num_profiles, num_heights)
885 jspectra = self.dataOut.data_spc[2]
885 jspectra = self.dataOut.data_spc[2]
886 jspectra_tmp = numpy.zeros(jspectra.shape)
886 jspectra_tmp = numpy.zeros(jspectra.shape)
887 num_profiles = jspectra.shape[0]
887 num_profiles = jspectra.shape[0]
888 freq_dc = int(num_profiles / 2)
888 freq_dc = int(num_profiles / 2)
889 # Flip con for
889 # Flip con for
890 for j in range(num_profiles):
890 for j in range(num_profiles):
891 jspectra_tmp[num_profiles-j-1]= jspectra[j]
891 jspectra_tmp[num_profiles-j-1]= jspectra[j]
892 # Intercambio perfil de DC con perfil inmediato anterior
892 # Intercambio perfil de DC con perfil inmediato anterior
893 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
893 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
894 jspectra_tmp[freq_dc]= jspectra[freq_dc]
894 jspectra_tmp[freq_dc]= jspectra[freq_dc]
895 # canal modificado es re-escrito en el arreglo de canales
895 # canal modificado es re-escrito en el arreglo de canales
896 self.dataOut.data_spc[2] = jspectra_tmp
896 self.dataOut.data_spc[2] = jspectra_tmp
897
897
898 return self.dataOut No newline at end of file
898 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now