##// END OF EJS Templates
schain
RSS
prueba de con ancho espectral y adicion de snr
avaldez -
r1307:7a6e3ce2f2d5
parent child
Show More
Side by Side Unified
Context file:
r1307:7a6e3ce2f2d5 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,526 +1,512
1 1 import numpy,math,random,time
2 2 #---------------1 Heredamos JRODatareader
3 3 from schainpy.model.io.jroIO_base import *
4 4 #---------------2 Heredamos las propiedades de ProcessingUnit
5 5 from schainpy.model.proc.jroproc_base import ProcessingUnit,Operation,MPDecorator
6 6 #---------------3 Importaremos las clases BascicHeader, SystemHeader, RadarControlHeader, ProcessingHeader
7 7 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader,SystemHeader,RadarControllerHeader, ProcessingHeader
8 8 #---------------4 Importaremos el objeto Voltge
9 9 from schainpy.model.data.jrodata import Voltage
10 10
11 11 class SimulatorReader(JRODataReader, ProcessingUnit):
12 12 incIntFactor = 1
13 13 nFFTPoints = 0
14 14 FixPP_IncInt = 1
15 15 FixRCP_IPP = 1000
16 16 FixPP_CohInt = 1
17 17 Tau_0 = 250
18 18 AcqH0_0 = 70
19 19 H0 = AcqH0_0
20 20 AcqDH_0 = 1.25
21 21 DH0 = AcqDH_0
22 22 Bauds = 32
23 23 BaudWidth = None
24 24 FixRCP_TXA = 40
25 25 FixRCP_TXB = 70
26 26 fAngle = 2.0*math.pi*(1/16)
27 27 DC_level = 500
28 28 stdev = 8
29 29 Num_Codes = 2
30 30 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
31 31 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
32 32 #Dyn_snCode = numpy.array([Num_Codes,Bauds])
33 33 Dyn_snCode = None
34 34 Samples = 200
35 35 channels = 2
36 36 pulses = None
37 37 Reference = None
38 38 pulse_size = None
39 39 prof_gen = None
40 40 Fdoppler = 100
41 41 Hdoppler = 36
42 42 Adoppler = 300
43 43 frequency = 9345
44 44 nTotalReadFiles = 1000
45 45
46 46 def __init__(self):
47 47 """
48 48 Inicializador de la clases SimulatorReader para
49 49 generar datos de voltage simulados.
50 50 Input:
51 51 dataOut: Objeto de la clase Voltage.
52 52 Este Objeto sera utilizado apra almacenar
53 53 un perfil de datos cada vez qe se haga
54 54 un requerimiento (getData)
55 55 """
56 56 ProcessingUnit.__init__(self)
57 57 print(" [ START ] init - Metodo Simulator Reader")
58 58
59 59 self.isConfig = False
60 60 self.basicHeaderObj = BasicHeader(LOCALTIME)
61 61 self.systemHeaderObj = SystemHeader()
62 62 self.radarControllerHeaderObj = RadarControllerHeader()
63 63 self.processingHeaderObj = ProcessingHeader()
64 64 self.profileIndex = 2**32-1
65 65 self.dataOut = Voltage()
66 66 #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1])
67 67 code0 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,1,1,1,-1,1,1,-1,1,-1,-1,-1,1,1,1,-1,1])
68 68 #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0])
69 69 code1 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,-1,-1,1,-1,-1,1,-1,1,1,1,-1,-1,-1,1,-1])
70 70 #self.Dyn_snCode = numpy.array([code0,code1])
71 71 self.Dyn_snCode = None
72 72
73 73 def set_kwargs(self, **kwargs):
74 74 for key, value in kwargs.items():
75 75 setattr(self, key, value)
76 76
77 77 def __hasNotDataInBuffer(self):
78 78
79 79 if self.profileIndex >= self.processingHeaderObj.profilesPerBlock* self.nTxs:
80 80 if self.nReadBlocks>0:
81 81 tmp = self.dataOut.utctime
82 82 tmp_utc = int(self.dataOut.utctime)
83 83 tmp_milisecond = int((tmp-tmp_utc)*1000)
84 84 self.basicHeaderObj.utc = tmp_utc
85 85 self.basicHeaderObj.miliSecond= tmp_milisecond
86 86 return 1
87 87 return 0
88 88
89 89 def setNextFile(self):
90 90 """Set the next file to be readed open it and parse de file header"""
91 91
92 92 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
93 93 self.nReadFiles=self.nReadFiles+1
94 94 if self.nReadFiles > self.nTotalReadFiles:
95 95 self.flagNoMoreFiles=1
96 96 raise schainpy.admin.SchainWarning('No more files to read')
97 97
98 98 print('------------------- [Opening file] ------------------------------',self.nReadFiles)
99 99 self.nReadBlocks = 0
100 100 #if self.nReadBlocks==0:
101 101 # self.readFirstHeader()
102 102
103 103 def __setNewBlock(self):
104 104 self.setNextFile()
105 105 if self.flagIsNewFile:
106 106 return 1
107 107
108 108 def readNextBlock(self):
109 109 while True:
110 110 self.__setNewBlock()
111 111 if not(self.readBlock()):
112 112 return 0
113 113 self.getBasicHeader()
114 114 break
115 115 if self.verbose:
116 116 print("[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
117 117 self.processingHeaderObj.dataBlocksPerFile,
118 118 self.dataOut.datatime.ctime()) )
119 119 return 1
120 120
121 121 def getFirstHeader(self):
122 122 self.getBasicHeader()
123 123 self.dataOut.processingHeaderObj = self.processingHeaderObj.copy()
124 124 self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
125 125 self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
126 126 self.dataOut.dtype = self.dtype
127 127
128 128 self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
129 129 self.dataOut.heightList = numpy.arange(self.processingHeaderObj.nHeights) * self.processingHeaderObj.deltaHeight + self.processingHeaderObj.firstHeight
130 130 self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels))
131 131 self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
132 132 # asumo q la data no esta decodificada
133 133 self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode
134 134 # asumo q la data no esta sin flip
135 135 self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip
136 136 self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft
137 137 self.dataOut.frequency = self.frequency
138 138
139 139 def getBasicHeader(self):
140 140 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
141 141 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
142 142
143 143 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
144 144 self.dataOut.timeZone = self.basicHeaderObj.timeZone
145 145 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
146 146 self.dataOut.errorCount = self.basicHeaderObj.errorCount
147 147 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
148 148 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
149 149
150 150 def readFirstHeader(self):
151 151
152 152 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
153 153 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
154 154 if datatype == 0:
155 155 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
156 156 elif datatype == 1:
157 157 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
158 158 elif datatype == 2:
159 159 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
160 160 elif datatype == 3:
161 161 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
162 162 elif datatype == 4:
163 163 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
164 164 elif datatype == 5:
165 165 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
166 166 else:
167 167 raise ValueError('Data type was not defined')
168 168
169 169 self.dtype = datatype_str
170 170
171 171
172 172 def set_RCH(self, expType=2, nTx=1,ipp=None, txA=0, txB=0,
173 173 nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
174 174 numTaus=0, line6Function=0, line5Function=0, fClock=None,
175 175 prePulseBefore=0, prePulseAfter=0,
176 176 codeType=0, nCode=0, nBaud=0, code=None,
177 177 flip1=0, flip2=0,Taus=0):
178 178 self.radarControllerHeaderObj.expType = expType
179 179 self.radarControllerHeaderObj.nTx = nTx
180 180 self.radarControllerHeaderObj.ipp = float(ipp)
181 181 self.radarControllerHeaderObj.txA = float(txA)
182 182 self.radarControllerHeaderObj.txB = float(txB)
183 183 self.radarControllerHeaderObj.rangeIpp = b'A\n'#ipp
184 184 self.radarControllerHeaderObj.rangeTxA = b''
185 185 self.radarControllerHeaderObj.rangeTxB = b''
186 186
187 187 self.radarControllerHeaderObj.nHeights = int(nHeights)
188 188 self.radarControllerHeaderObj.firstHeight = numpy.array([firstHeight])
189 189 self.radarControllerHeaderObj.deltaHeight = numpy.array([deltaHeight])
190 190 self.radarControllerHeaderObj.samplesWin = numpy.array([nHeights])
191 191
192 192
193 193 self.radarControllerHeaderObj.nWindows = nWindows
194 194 self.radarControllerHeaderObj.numTaus = numTaus
195 195 self.radarControllerHeaderObj.codeType = codeType
196 196 self.radarControllerHeaderObj.line6Function = line6Function
197 197 self.radarControllerHeaderObj.line5Function = line5Function
198 198 #self.radarControllerHeaderObj.fClock = fClock
199 199 self.radarControllerHeaderObj.prePulseBefore= prePulseBefore
200 200 self.radarControllerHeaderObj.prePulseAfter = prePulseAfter
201 201
202 202 self.radarControllerHeaderObj.flip1 = flip1
203 203 self.radarControllerHeaderObj.flip2 = flip2
204 204
205 205 self.radarControllerHeaderObj.code_size = 0
206 206 if self.radarControllerHeaderObj.codeType != 0:
207 207 self.radarControllerHeaderObj.nCode = nCode
208 208 self.radarControllerHeaderObj.nBaud = nBaud
209 209 self.radarControllerHeaderObj.code = code
210 210 self.radarControllerHeaderObj.code_size = int(numpy.ceil(nBaud / 32.)) * nCode * 4
211 211
212 212 if fClock is None and deltaHeight is not None:
213 213 self.fClock = 0.15 / (deltaHeight * 1e-6)
214 214 self.radarControllerHeaderObj.fClock = self.fClock
215 215 if numTaus==0:
216 216 self.radarControllerHeaderObj.Taus = numpy.array(0,'<f4')
217 217 else:
218 218 self.radarControllerHeaderObj.Taus = numpy.array(Taus,'<f4')
219 219
220 220 def set_PH(self, dtype=0, blockSize=0, profilesPerBlock=0,
221 221 dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0,
222 222 nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0,
223 223 deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
224 224 code=0, nBaud=None, shif_fft=False, flag_dc=False,
225 225 flag_cspc=False, flag_decode=False, flag_deflip=False):
226 226
227 227 self.processingHeaderObj.dtype = dtype
228 228 self.processingHeaderObj.profilesPerBlock = profilesPerBlock
229 229 self.processingHeaderObj.dataBlocksPerFile = dataBlocksPerFile
230 230 self.processingHeaderObj.nWindows = nWindows
231 231 self.processingHeaderObj.processFlags = processFlags
232 232 self.processingHeaderObj.nCohInt = nCohInt
233 233 self.processingHeaderObj.nIncohInt = nIncohInt
234 234 self.processingHeaderObj.totalSpectra = totalSpectra
235 235
236 236 self.processingHeaderObj.nHeights = int(nHeights)
237 237 self.processingHeaderObj.firstHeight = firstHeight#numpy.array([firstHeight])#firstHeight
238 238 self.processingHeaderObj.deltaHeight = deltaHeight#numpy.array([deltaHeight])#deltaHeight
239 239 self.processingHeaderObj.samplesWin = nHeights#numpy.array([nHeights])#nHeights
240 240
241 241 def set_BH(self, utc = 0, miliSecond = 0, timeZone = 0):
242 242 self.basicHeaderObj.utc = utc
243 243 self.basicHeaderObj.miliSecond = miliSecond
244 244 self.basicHeaderObj.timeZone = timeZone
245 245
246 246 def set_SH(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=32):
247 247 #self.systemHeaderObj.size = size
248 248 self.systemHeaderObj.nSamples = nSamples
249 249 self.systemHeaderObj.nProfiles = nProfiles
250 250 self.systemHeaderObj.nChannels = nChannels
251 251 self.systemHeaderObj.adcResolution = adcResolution
252 252 self.systemHeaderObj.pciDioBusWidth = pciDioBusWidth
253 253
254 254 def init_acquisition(self):
255 255
256 256 if self.nFFTPoints != 0:
257 257 self.incIntFactor = m_nProfilesperBlock/self.nFFTPoints
258 258 if (self.FixPP_IncInt > self.incIntFactor):
259 259 self.incIntFactor = self.FixPP_IncInt/ self.incIntFactor
260 260 elif(self.FixPP_IncInt< self.incIntFactor):
261 261 print("False alert...")
262 262
263 263 ProfilesperBlock = self.processingHeaderObj.profilesPerBlock
264 264
265 265 self.timeperblock =int(((self.FixRCP_IPP
266 266 *ProfilesperBlock
267 267 *self.FixPP_CohInt
268 268 *self.incIntFactor)
269 269 /150.0)
270 270 *0.9
271 271 +0.5)
272 272 # para cada canal
273 273 self.profiles = ProfilesperBlock*self.FixPP_CohInt
274 274 self.profiles = ProfilesperBlock
275 275 self.Reference = int((self.Tau_0-self.AcqH0_0)/(self.AcqDH_0)+0.5)
276 276 self.BaudWidth = int((self.FixRCP_TXA/self.AcqDH_0)/self.Bauds + 0.5 )
277 277
278 278 if (self.BaudWidth==0):
279 279 self.BaudWidth=1
280 280
281 281 def init_pulse(self,Num_Codes=Num_Codes,Bauds=Bauds,BaudWidth=BaudWidth,Dyn_snCode=Dyn_snCode):
282 282
283 283 Num_Codes = Num_Codes
284 284 Bauds = Bauds
285 285 BaudWidth = BaudWidth
286 286 Dyn_snCode = Dyn_snCode
287 287
288 288 if Dyn_snCode:
289 289 print("EXISTE")
290 290 else:
291 291 print("No existe")
292 292
293 293 if Dyn_snCode: # if Bauds:
294 294 pulses = list(range(0,Num_Codes))
295 295 num_codes = Num_Codes
296 296 for i in range(num_codes):
297 297 pulse_size = Bauds*BaudWidth
298 298 pulses[i] = numpy.zeros(pulse_size)
299 299 for j in range(Bauds):
300 300 for k in range(BaudWidth):
301 301 pulses[i][j*BaudWidth+k] = int(Dyn_snCode[i][j]*600)
302 302 else:
303 303 print("sin code")
304 304 pulses = list(range(1))
305 305 if self.AcqDH_0>0.149:
306 306 pulse_size = int(self.FixRCP_TXB/0.15+0.5)
307 307 else:
308 308 pulse_size = int((self.FixRCP_TXB/self.AcqDH_0)+0.5) #0.0375
309 309 pulses[0] = numpy.ones(pulse_size)
310 310 pulses = 600*pulses[0]
311 311
312 312 return pulses,pulse_size
313 313
314 314 def jro_GenerateBlockOfData(self,Samples=Samples,DC_level= DC_level,stdev=stdev,
315 315 Reference= Reference,pulses= pulses,
316 316 Num_Codes= Num_Codes,pulse_size=pulse_size,
317 317 prof_gen= prof_gen,H0 = H0,DH0=DH0,
318 318 Adoppler=Adoppler,Fdoppler= Fdoppler,Hdoppler=Hdoppler):
319 319 Samples = Samples
320 320 DC_level = DC_level
321 321 stdev = stdev
322 322 m_nR = Reference
323 323 pulses = pulses
324 324 num_codes = Num_Codes
325 325 ps = pulse_size
326 326 prof_gen = prof_gen
327 327 channels = self.channels
328 328 H0 = H0
329 329 DH0 = DH0
330 330 ippSec = self.radarControllerHeaderObj.ippSeconds
331 331 Fdoppler = self.Fdoppler
332 332 Hdoppler = self.Hdoppler
333 333 Adoppler = self.Adoppler
334 334
335 335 self.datablock = numpy.zeros([channels,prof_gen,Samples],dtype= numpy.complex64)
336 336 for i in range(channels):
337 337 for k in range(prof_gen):
338 338 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
339 339 Noise_r = numpy.random.normal(DC_level,stdev,Samples)
340 340 Noise_i = numpy.random.normal(DC_level,stdev,Samples)
341 341 Noise = numpy.zeros(Samples,dtype=complex)
342 342 Noise.real = Noise_r
343 343 Noise.imag = Noise_i
344 344 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·PULSOSΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
345 345 Pulso = numpy.zeros(pulse_size,dtype=complex)
346 346 Pulso.real = pulses[k%num_codes]
347 347 Pulso.imag = pulses[k%num_codes]
348 348 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· PULSES+NOISEΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·
349 349 InBuffer = numpy.zeros(Samples,dtype=complex)
350 350 InBuffer[m_nR:m_nR+ps] = Pulso
351 351 InBuffer = InBuffer+Noise
352 352 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· ANGLE Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
353 353 InBuffer.real[m_nR:m_nR+ps] = InBuffer.real[m_nR:m_nR+ps]*(math.cos( self.fAngle)*5)
354 354 InBuffer.imag[m_nR:m_nR+ps] = InBuffer.imag[m_nR:m_nR+ps]*(math.sin( self.fAngle)*5)
355 355 InBuffer=InBuffer
356 356 self.datablock[i][k]= InBuffer
357 #plot_cts(InBuffer,H0=H0,DH0=DH0)
358 #wave_fft(x=InBuffer,plot_show=True)
359 #time.sleep(1)
357
360 358 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·DOPPLER SIGNAL...............................................
361 359 time_vec = numpy.linspace(0,(prof_gen-1)*ippSec,int(prof_gen))+self.nReadBlocks*ippSec*prof_gen+(self.nReadFiles-1)*ippSec*prof_gen
362 360 fd = Fdoppler #+(600.0/120)*self.nReadBlocks
363 361 d_signal = Adoppler*numpy.array(numpy.exp(1.0j*2.0*math.pi*fd*time_vec),dtype=numpy.complex64)
364 362 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·SeΓ±al con ancho espectralΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
365 #specw_sig = numpy.zeros(int(prof_gen),dtype=complex)
366 #specw_sig.real[200:200+100] = 1*numpy.ones(100)
367 #specw_sig.imag[200:200+100] = 1*numpy.ones(100)
368 # w=2
369 specw_sig = numpy.linspace(-149,150,300)
370 w = 3
371 A = 10
372 specw_sig = specw_sig/w
373 specw_sig = numpy.sinc(specw_sig)
374 specw_sig = A*numpy.array(specw_sig,dtype=numpy.complex64)
363 #specw_sig = numpy.linspace(-149,150,300)
364 #w = 8
365 #A = 20
366 #specw_sig = specw_sig/w
367 #specw_sig = numpy.sinc(specw_sig)
368 #specw_sig = A*numpy.array(specw_sig,dtype=numpy.complex64)
375 369 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DATABLOCK + DOPPLERΒ·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
376 370 HD=int(Hdoppler/self.AcqDH_0)
377 371 for i in range(12):
378 372 self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ d_signal# RESULT
379
380 HD=int(Hdoppler/self.AcqDH_0)
381 HD=int(HD/2)
382 for i in range(12):
383 self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ specw_sig*d_signal# RESULT
384
385 '''
386 a= numpy.zeros(10)
387 for i in range(10):
388 a[i]=i+self.nReadBlocks+20
389 for i in a:
390 self.datablock[0,:,int(i)]=self.datablock[0,:,int(i)]+ d_signal # RESULT
391 '''
373 #Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β· DATABLOCK + DOPPLER*Sinc(x)Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·Β·
374 #HD=int(Hdoppler/self.AcqDH_0)
375 #HD=int(HD/2)
376 #for i in range(12):
377 # self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ specw_sig*d_signal# RESULT
392 378
393 379 def readBlock(self):
394 380
395 381 self.jro_GenerateBlockOfData(Samples= self.samples,DC_level=self.DC_level,
396 382 stdev=self.stdev,Reference= self.Reference,
397 383 pulses = self.pulses,Num_Codes=self.Num_Codes,
398 384 pulse_size=self.pulse_size,prof_gen=self.profiles,
399 385 H0=self.H0,DH0=self.DH0)
400 386
401 387 self.profileIndex = 0
402 388 self.flagIsNewFile = 0
403 389 self.flagIsNewBlock = 1
404 390 self.nTotalBlocks += 1
405 391 self.nReadBlocks += 1
406 392
407 393 return 1
408 394
409 395
410 396 def getData(self):
411 397 if self.flagNoMoreFiles:
412 398 self.dataOut.flagNodata = True
413 399 return 0
414 400 self.flagDiscontinuousBlock = 0
415 401 self.flagIsNewBlock = 0
416 402 if self.__hasNotDataInBuffer(): # aqui es verdad
417 403 if not(self.readNextBlock()): # return 1 y por eso el if not salta a getBasic Header
418 404 return 0
419 405 self.getFirstHeader() # atributo
420 406
421 407 if not self.getByBlock:
422 408 self.dataOut.flagDataAsBlock = False
423 409 self.dataOut.data = self.datablock[:, self.profileIndex, :]
424 410 self.dataOut.profileIndex = self.profileIndex
425 411 self.profileIndex += 1
426 412 else:
427 413 pass
428 414 self.dataOut.flagNoData = False
429 415 self.getBasicHeader()
430 416 self.dataOut.realtime = self.online
431 417 return self.dataOut.data
432 418
433 419
434 420 def setup(self,frequency=49.92e6,incIntFactor= 1, nFFTPoints = 0, FixPP_IncInt=1,FixRCP_IPP=1000,
435 421 FixPP_CohInt= 1,Tau_0= 250,AcqH0_0 = 70 ,AcqDH_0=1.25, Bauds= 32,
436 422 FixRCP_TXA = 40, FixRCP_TXB = 50, fAngle = 2.0*math.pi*(1/16),DC_level= 50,
437 423 stdev= 8,Num_Codes = 1 , Dyn_snCode = None, samples=200,
438 424 channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,nTotalReadFiles=10000,
439 425 **kwargs):
440 426
441 427 self.set_kwargs(**kwargs)
442 428 self.nReadBlocks = 0
443 429 self.nReadFiles = 1
444 430 print('------------------- [Opening file: ] ------------------------------',self.nReadFiles)
445 431
446 432 tmp = time.time()
447 433 tmp_utc = int(tmp)
448 434 tmp_milisecond = int((tmp-tmp_utc)*1000)
449 435 print(" SETUP -basicHeaderObj.utc",datetime.datetime.utcfromtimestamp(tmp))
450 436 if Dyn_snCode is None:
451 437 Num_Codes=1
452 438 Bauds =1
453 439
454 440
455 441
456 442 self.set_BH(utc= tmp_utc,miliSecond= tmp_milisecond,timeZone=300 )
457 443 self.set_RCH( expType=0, nTx=150,ipp=FixRCP_IPP, txA=FixRCP_TXA, txB= FixRCP_TXB,
458 444 nWindows=1 , nHeights=samples, firstHeight=AcqH0_0, deltaHeight=AcqDH_0,
459 445 numTaus=1, line6Function=0, line5Function=0, fClock=None,
460 446 prePulseBefore=0, prePulseAfter=0,
461 447 codeType=0, nCode=Num_Codes, nBaud=32, code=Dyn_snCode,
462 448 flip1=0, flip2=0,Taus=Tau_0)
463 449
464 450 self.set_PH(dtype=0, blockSize=0, profilesPerBlock=300,
465 451 dataBlocksPerFile=120, nWindows=1, processFlags=numpy.array([1024]), nCohInt=1,
466 452 nIncohInt=1, totalSpectra=0, nHeights=samples, firstHeight=AcqH0_0,
467 453 deltaHeight=AcqDH_0, samplesWin=samples, spectraComb=0, nCode=0,
468 454 code=0, nBaud=None, shif_fft=False, flag_dc=False,
469 455 flag_cspc=False, flag_decode=False, flag_deflip=False)
470 456
471 457 self.set_SH(nSamples=samples, nProfiles=300, nChannels=channels)
472 458
473 459 self.readFirstHeader()
474 460
475 461 self.frequency = frequency
476 462 self.incIntFactor = incIntFactor
477 463 self.nFFTPoints = nFFTPoints
478 464 self.FixPP_IncInt = FixPP_IncInt
479 465 self.FixRCP_IPP = FixRCP_IPP
480 466 self.FixPP_CohInt = FixPP_CohInt
481 467 self.Tau_0 = Tau_0
482 468 self.AcqH0_0 = AcqH0_0
483 469 self.H0 = AcqH0_0
484 470 self.AcqDH_0 = AcqDH_0
485 471 self.DH0 = AcqDH_0
486 472 self.Bauds = Bauds
487 473 self.FixRCP_TXA = FixRCP_TXA
488 474 self.FixRCP_TXB = FixRCP_TXB
489 475 self.fAngle = fAngle
490 476 self.DC_level = DC_level
491 477 self.stdev = stdev
492 478 self.Num_Codes = Num_Codes
493 479 self.Dyn_snCode = Dyn_snCode
494 480 self.samples = samples
495 481 self.channels = channels
496 482 self.profiles = None
497 483 self.m_nReference = None
498 484 self.Baudwidth = None
499 485 self.Fdoppler = Fdoppler
500 486 self.Hdoppler = Hdoppler
501 487 self.Adoppler = Adoppler
502 488 self.nTotalReadFiles = int(nTotalReadFiles)
503 489
504 490 print("IPP ", self.FixRCP_IPP)
505 491 print("Tau_0 ",self.Tau_0)
506 492 print("AcqH0_0",self.AcqH0_0)
507 493 print("samples,window ",self.samples)
508 494 print("AcqDH_0",AcqDH_0)
509 495 print("FixRCP_TXA",self.FixRCP_TXA)
510 496 print("FixRCP_TXB",self.FixRCP_TXB)
511 497 print("Dyn_snCode",Dyn_snCode)
512 498 print("Fdoppler", Fdoppler)
513 499 print("Hdoppler",Hdoppler)
514 500 print("Vdopplermax",Fdoppler*(3.0e8/self.frequency)/2.0)
515 501 print("nTotalReadFiles", nTotalReadFiles)
516 502
517 503 self.init_acquisition()
518 504 self.pulses,self.pulse_size=self.init_pulse(Num_Codes=self.Num_Codes,Bauds=self.Bauds,BaudWidth=self.BaudWidth,Dyn_snCode=Dyn_snCode)
519 505 print(" [ END ] - SETUP metodo")
520 506 return
521 507
522 508 def run(self,**kwargs): # metodo propio
523 509 if not(self.isConfig):
524 510 self.setup(**kwargs)
525 511 self.isConfig = True
526 512 self.getData()
Show file before | Show file after | Hide comments
@@ -1,1581 +1,1587
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 167 if minHei == None:
168 168 minHei = self.dataOut.heightList[0]
169 169
170 170 if maxHei == None:
171 171 maxHei = self.dataOut.heightList[-1]
172 172
173 173 if (minHei < self.dataOut.heightList[0]):
174 174 minHei = self.dataOut.heightList[0]
175 175
176 176 if (maxHei > self.dataOut.heightList[-1]):
177 177 maxHei = self.dataOut.heightList[-1]
178 178
179 179 minIndex = 0
180 180 maxIndex = 0
181 181 heights = self.dataOut.heightList
182 182
183 183 inda = numpy.where(heights >= minHei)
184 184 indb = numpy.where(heights <= maxHei)
185 185
186 186 try:
187 187 minIndex = inda[0][0]
188 188 except:
189 189 minIndex = 0
190 190
191 191 try:
192 192 maxIndex = indb[0][-1]
193 193 except:
194 194 maxIndex = len(heights)
195 195
196 196 self.selectHeightsByIndex(minIndex, maxIndex)
197 197
198 198 return self.dataOut
199 199
200 200 def selectHeightsByIndex(self, minIndex, maxIndex):
201 201 """
202 202 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
203 203 minIndex <= index <= maxIndex
204 204
205 205 Input:
206 206 minIndex : valor de indice minimo de altura a considerar
207 207 maxIndex : valor de indice maximo de altura a considerar
208 208
209 209 Affected:
210 210 self.dataOut.data
211 211 self.dataOut.heightList
212 212
213 213 Return:
214 214 1 si el metodo se ejecuto con exito caso contrario devuelve 0
215 215 """
216 216
217 217 if self.dataOut.type == 'Voltage':
218 218 if (minIndex < 0) or (minIndex > maxIndex):
219 219 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
220 220
221 221 if (maxIndex >= self.dataOut.nHeights):
222 222 maxIndex = self.dataOut.nHeights
223 223
224 224 #voltage
225 225 if self.dataOut.flagDataAsBlock:
226 226 """
227 227 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
228 228 """
229 229 data = self.dataOut.data[:,:, minIndex:maxIndex]
230 230 else:
231 231 data = self.dataOut.data[:, minIndex:maxIndex]
232 232
233 233 # firstHeight = self.dataOut.heightList[minIndex]
234 234
235 235 self.dataOut.data = data
236 236 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
237 237
238 238 if self.dataOut.nHeights <= 1:
239 239 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
240 240 elif self.dataOut.type == 'Spectra':
241 241 if (minIndex < 0) or (minIndex > maxIndex):
242 242 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
243 243 minIndex, maxIndex))
244 244
245 245 if (maxIndex >= self.dataOut.nHeights):
246 246 maxIndex = self.dataOut.nHeights - 1
247 247
248 248 # Spectra
249 249 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
250 250
251 251 data_cspc = None
252 252 if self.dataOut.data_cspc is not None:
253 253 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
254 254
255 255 data_dc = None
256 256 if self.dataOut.data_dc is not None:
257 257 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
258 258
259 259 self.dataOut.data_spc = data_spc
260 260 self.dataOut.data_cspc = data_cspc
261 261 self.dataOut.data_dc = data_dc
262 262
263 263 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
264 264
265 265 return 1
266 266
267 267
268 268 class filterByHeights(Operation):
269 269
270 270 def run(self, dataOut, window):
271 271
272 272 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
273 273
274 274 if window == None:
275 275 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
276 276
277 277 newdelta = deltaHeight * window
278 278 r = dataOut.nHeights % window
279 279 newheights = (dataOut.nHeights-r)/window
280 280
281 281 if newheights <= 1:
282 282 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
283 283
284 284 if dataOut.flagDataAsBlock:
285 285 """
286 286 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
287 287 """
288 288 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
289 289 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
290 290 buffer = numpy.sum(buffer,3)
291 291
292 292 else:
293 293 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
294 294 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
295 295 buffer = numpy.sum(buffer,2)
296 296
297 297 dataOut.data = buffer
298 298 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
299 299 dataOut.windowOfFilter = window
300 300
301 301 return dataOut
302 302
303 303
304 304 class setH0(Operation):
305 305
306 306 def run(self, dataOut, h0, deltaHeight = None):
307 307
308 308 if not deltaHeight:
309 309 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
310 310
311 311 nHeights = dataOut.nHeights
312 312
313 313 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
314 314
315 315 dataOut.heightList = newHeiRange
316 316
317 317 return dataOut
318 318
319 319
320 320 class deFlip(Operation):
321 321
322 322 def run(self, dataOut, channelList = []):
323 323
324 324 data = dataOut.data.copy()
325 325
326 326 if dataOut.flagDataAsBlock:
327 327 flip = self.flip
328 328 profileList = list(range(dataOut.nProfiles))
329 329
330 330 if not channelList:
331 331 for thisProfile in profileList:
332 332 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
333 333 flip *= -1.0
334 334 else:
335 335 for thisChannel in channelList:
336 336 if thisChannel not in dataOut.channelList:
337 337 continue
338 338
339 339 for thisProfile in profileList:
340 340 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
341 341 flip *= -1.0
342 342
343 343 self.flip = flip
344 344
345 345 else:
346 346 if not channelList:
347 347 data[:,:] = data[:,:]*self.flip
348 348 else:
349 349 for thisChannel in channelList:
350 350 if thisChannel not in dataOut.channelList:
351 351 continue
352 352
353 353 data[thisChannel,:] = data[thisChannel,:]*self.flip
354 354
355 355 self.flip *= -1.
356 356
357 357 dataOut.data = data
358 358
359 359 return dataOut
360 360
361 361
362 362 class setAttribute(Operation):
363 363 '''
364 364 Set an arbitrary attribute(s) to dataOut
365 365 '''
366 366
367 367 def __init__(self):
368 368
369 369 Operation.__init__(self)
370 370 self._ready = False
371 371
372 372 def run(self, dataOut, **kwargs):
373 373
374 374 for key, value in kwargs.items():
375 375 setattr(dataOut, key, value)
376 376
377 377 return dataOut
378 378
379 379
380 380 class interpolateHeights(Operation):
381 381
382 382 def run(self, dataOut, topLim, botLim):
383 383 #69 al 72 para julia
384 384 #82-84 para meteoros
385 385 if len(numpy.shape(dataOut.data))==2:
386 386 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
387 387 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
388 388 #dataOut.data[:,botLim:limSup+1] = sampInterp
389 389 dataOut.data[:,botLim:topLim+1] = sampInterp
390 390 else:
391 391 nHeights = dataOut.data.shape[2]
392 392 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
393 393 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
394 394 f = interpolate.interp1d(x, y, axis = 2)
395 395 xnew = numpy.arange(botLim,topLim+1)
396 396 ynew = f(xnew)
397 397 dataOut.data[:,:,botLim:topLim+1] = ynew
398 398
399 399 return dataOut
400 400
401 401
402 402 class CohInt(Operation):
403 403
404 404 isConfig = False
405 405 __profIndex = 0
406 406 __byTime = False
407 407 __initime = None
408 408 __lastdatatime = None
409 409 __integrationtime = None
410 410 __buffer = None
411 411 __bufferStride = []
412 412 __dataReady = False
413 413 __profIndexStride = 0
414 414 __dataToPutStride = False
415 415 n = None
416 416
417 417 def __init__(self, **kwargs):
418 418
419 419 Operation.__init__(self, **kwargs)
420 420
421 421 # self.isConfig = False
422 422
423 423 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
424 424 """
425 425 Set the parameters of the integration class.
426 426
427 427 Inputs:
428 428
429 429 n : Number of coherent integrations
430 430 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
431 431 overlapping :
432 432 """
433 433
434 434 self.__initime = None
435 435 self.__lastdatatime = 0
436 436 self.__buffer = None
437 437 self.__dataReady = False
438 438 self.byblock = byblock
439 439 self.stride = stride
440 440
441 441 if n == None and timeInterval == None:
442 442 raise ValueError("n or timeInterval should be specified ...")
443 443
444 444 if n != None:
445 445 self.n = n
446 446 self.__byTime = False
447 447 else:
448 448 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
449 449 self.n = 9999
450 450 self.__byTime = True
451 451
452 452 if overlapping:
453 453 self.__withOverlapping = True
454 454 self.__buffer = None
455 455 else:
456 456 self.__withOverlapping = False
457 457 self.__buffer = 0
458 458
459 459 self.__profIndex = 0
460 460
461 461 def putData(self, data):
462 462
463 463 """
464 464 Add a profile to the __buffer and increase in one the __profileIndex
465 465
466 466 """
467 467
468 468 if not self.__withOverlapping:
469 469 self.__buffer += data.copy()
470 470 self.__profIndex += 1
471 471 return
472 472
473 473 #Overlapping data
474 474 nChannels, nHeis = data.shape
475 475 data = numpy.reshape(data, (1, nChannels, nHeis))
476 476
477 477 #If the buffer is empty then it takes the data value
478 478 if self.__buffer is None:
479 479 self.__buffer = data
480 480 self.__profIndex += 1
481 481 return
482 482
483 483 #If the buffer length is lower than n then stakcing the data value
484 484 if self.__profIndex < self.n:
485 485 self.__buffer = numpy.vstack((self.__buffer, data))
486 486 self.__profIndex += 1
487 487 return
488 488
489 489 #If the buffer length is equal to n then replacing the last buffer value with the data value
490 490 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
491 491 self.__buffer[self.n-1] = data
492 492 self.__profIndex = self.n
493 493 return
494 494
495 495
496 496 def pushData(self):
497 497 """
498 498 Return the sum of the last profiles and the profiles used in the sum.
499 499
500 500 Affected:
501 501
502 502 self.__profileIndex
503 503
504 504 """
505 505
506 506 if not self.__withOverlapping:
507 507 data = self.__buffer
508 508 n = self.__profIndex
509 509
510 510 self.__buffer = 0
511 511 self.__profIndex = 0
512 512
513 513 return data, n
514 514
515 515 #Integration with Overlapping
516 516 data = numpy.sum(self.__buffer, axis=0)
517 517 # print data
518 518 # raise
519 519 n = self.__profIndex
520 520
521 521 return data, n
522 522
523 523 def byProfiles(self, data):
524 524
525 525 self.__dataReady = False
526 526 avgdata = None
527 527 # n = None
528 528 # print data
529 529 # raise
530 530 self.putData(data)
531 531
532 532 if self.__profIndex == self.n:
533 533 avgdata, n = self.pushData()
534 534 self.__dataReady = True
535 535
536 536 return avgdata
537 537
538 538 def byTime(self, data, datatime):
539 539
540 540 self.__dataReady = False
541 541 avgdata = None
542 542 n = None
543 543
544 544 self.putData(data)
545 545
546 546 if (datatime - self.__initime) >= self.__integrationtime:
547 547 avgdata, n = self.pushData()
548 548 self.n = n
549 549 self.__dataReady = True
550 550
551 551 return avgdata
552 552
553 553 def integrateByStride(self, data, datatime):
554 554 # print data
555 555 if self.__profIndex == 0:
556 556 self.__buffer = [[data.copy(), datatime]]
557 557 else:
558 558 self.__buffer.append([data.copy(),datatime])
559 559 self.__profIndex += 1
560 560 self.__dataReady = False
561 561
562 562 if self.__profIndex == self.n * self.stride :
563 563 self.__dataToPutStride = True
564 564 self.__profIndexStride = 0
565 565 self.__profIndex = 0
566 566 self.__bufferStride = []
567 567 for i in range(self.stride):
568 568 current = self.__buffer[i::self.stride]
569 569 data = numpy.sum([t[0] for t in current], axis=0)
570 570 avgdatatime = numpy.average([t[1] for t in current])
571 571 # print data
572 572 self.__bufferStride.append((data, avgdatatime))
573 573
574 574 if self.__dataToPutStride:
575 575 self.__dataReady = True
576 576 self.__profIndexStride += 1
577 577 if self.__profIndexStride == self.stride:
578 578 self.__dataToPutStride = False
579 579 # print self.__bufferStride[self.__profIndexStride - 1]
580 580 # raise
581 581 return self.__bufferStride[self.__profIndexStride - 1]
582 582
583 583
584 584 return None, None
585 585
586 586 def integrate(self, data, datatime=None):
587 587
588 588 if self.__initime == None:
589 589 self.__initime = datatime
590 590
591 591 if self.__byTime:
592 592 avgdata = self.byTime(data, datatime)
593 593 else:
594 594 avgdata = self.byProfiles(data)
595 595
596 596
597 597 self.__lastdatatime = datatime
598 598
599 599 if avgdata is None:
600 600 return None, None
601 601
602 602 avgdatatime = self.__initime
603 603
604 604 deltatime = datatime - self.__lastdatatime
605 605
606 606 if not self.__withOverlapping:
607 607 self.__initime = datatime
608 608 else:
609 609 self.__initime += deltatime
610 610
611 611 return avgdata, avgdatatime
612 612
613 613 def integrateByBlock(self, dataOut):
614 614
615 615 times = int(dataOut.data.shape[1]/self.n)
616 616 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
617 617
618 618 id_min = 0
619 619 id_max = self.n
620 620
621 621 for i in range(times):
622 622 junk = dataOut.data[:,id_min:id_max,:]
623 623 avgdata[:,i,:] = junk.sum(axis=1)
624 624 id_min += self.n
625 625 id_max += self.n
626 626
627 627 timeInterval = dataOut.ippSeconds*self.n
628 628 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
629 629 self.__dataReady = True
630 630 return avgdata, avgdatatime
631 631
632 632 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
633 633
634 634 if not self.isConfig:
635 635 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
636 636 self.isConfig = True
637 637
638 638 if dataOut.flagDataAsBlock:
639 639 """
640 640 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
641 641 """
642 642 avgdata, avgdatatime = self.integrateByBlock(dataOut)
643 643 dataOut.nProfiles /= self.n
644 644 else:
645 645 if stride is None:
646 646 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
647 647 else:
648 648 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
649 649
650 650
651 651 # dataOut.timeInterval *= n
652 652 dataOut.flagNoData = True
653 653
654 654 if self.__dataReady:
655 655 dataOut.data = avgdata
656 656 dataOut.nCohInt *= self.n
657 657 dataOut.utctime = avgdatatime
658 658 # print avgdata, avgdatatime
659 659 # raise
660 660 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
661 661 dataOut.flagNoData = False
662 662 return dataOut
663 663
664 664 class Decoder(Operation):
665 665
666 666 isConfig = False
667 667 __profIndex = 0
668 668
669 669 code = None
670 670
671 671 nCode = None
672 672 nBaud = None
673 673
674 674 def __init__(self, **kwargs):
675 675
676 676 Operation.__init__(self, **kwargs)
677 677
678 678 self.times = None
679 679 self.osamp = None
680 680 # self.__setValues = False
681 681 self.isConfig = False
682 682 self.setupReq = False
683 683 def setup(self, code, osamp, dataOut):
684 684
685 685 self.__profIndex = 0
686 686
687 687 self.code = code
688 688
689 689 self.nCode = len(code)
690 690 self.nBaud = len(code[0])
691 691
692 692 if (osamp != None) and (osamp >1):
693 693 self.osamp = osamp
694 694 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
695 695 self.nBaud = self.nBaud*self.osamp
696 696
697 697 self.__nChannels = dataOut.nChannels
698 698 self.__nProfiles = dataOut.nProfiles
699 699 self.__nHeis = dataOut.nHeights
700 700
701 701 if self.__nHeis < self.nBaud:
702 702 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
703 703
704 704 #Frequency
705 705 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
706 706
707 707 __codeBuffer[:,0:self.nBaud] = self.code
708 708
709 709 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
710 710
711 711 if dataOut.flagDataAsBlock:
712 712
713 713 self.ndatadec = self.__nHeis #- self.nBaud + 1
714 714
715 715 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
716 716
717 717 else:
718 718
719 719 #Time
720 720 self.ndatadec = self.__nHeis #- self.nBaud + 1
721 721
722 722 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
723 723
724 724 def __convolutionInFreq(self, data):
725 725
726 726 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
727 727
728 728 fft_data = numpy.fft.fft(data, axis=1)
729 729
730 730 conv = fft_data*fft_code
731 731
732 732 data = numpy.fft.ifft(conv,axis=1)
733 733
734 734 return data
735 735
736 736 def __convolutionInFreqOpt(self, data):
737 737
738 738 raise NotImplementedError
739 739
740 740 def __convolutionInTime(self, data):
741 741
742 742 code = self.code[self.__profIndex]
743 743 for i in range(self.__nChannels):
744 744 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
745 745
746 746 return self.datadecTime
747 747
748 748 def __convolutionByBlockInTime(self, data):
749 749
750 750 repetitions = int(self.__nProfiles / self.nCode)
751 751 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
752 752 junk = junk.flatten()
753 753 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
754 754 profilesList = range(self.__nProfiles)
755 755
756 756 for i in range(self.__nChannels):
757 757 for j in profilesList:
758 758 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
759 759 return self.datadecTime
760 760
761 761 def __convolutionByBlockInFreq(self, data):
762 762
763 763 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
764 764
765 765
766 766 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
767 767
768 768 fft_data = numpy.fft.fft(data, axis=2)
769 769
770 770 conv = fft_data*fft_code
771 771
772 772 data = numpy.fft.ifft(conv,axis=2)
773 773
774 774 return data
775 775
776 776
777 777 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
778 778
779 779 if dataOut.flagDecodeData:
780 780 print("This data is already decoded, recoding again ...")
781 781
782 782 if not self.isConfig:
783 783
784 784 if code is None:
785 785 if dataOut.code is None:
786 786 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
787 787
788 788 code = dataOut.code
789 789 else:
790 790 code = numpy.array(code).reshape(nCode,nBaud)
791 791 self.setup(code, osamp, dataOut)
792 792
793 793 self.isConfig = True
794 794
795 795 if mode == 3:
796 796 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
797 797
798 798 if times != None:
799 799 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
800 800
801 801 if self.code is None:
802 802 print("Fail decoding: Code is not defined.")
803 803 return
804 804
805 805 self.__nProfiles = dataOut.nProfiles
806 806 datadec = None
807 807
808 808 if mode == 3:
809 809 mode = 0
810 810
811 811 if dataOut.flagDataAsBlock:
812 812 """
813 813 Decoding when data have been read as block,
814 814 """
815 815
816 816 if mode == 0:
817 817 datadec = self.__convolutionByBlockInTime(dataOut.data)
818 818 if mode == 1:
819 819 datadec = self.__convolutionByBlockInFreq(dataOut.data)
820 820 else:
821 821 """
822 822 Decoding when data have been read profile by profile
823 823 """
824 824 if mode == 0:
825 825 datadec = self.__convolutionInTime(dataOut.data)
826 826
827 827 if mode == 1:
828 828 datadec = self.__convolutionInFreq(dataOut.data)
829 829
830 830 if mode == 2:
831 831 datadec = self.__convolutionInFreqOpt(dataOut.data)
832 832
833 833 if datadec is None:
834 834 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
835 835
836 836 dataOut.code = self.code
837 837 dataOut.nCode = self.nCode
838 838 dataOut.nBaud = self.nBaud
839 839
840 840 dataOut.data = datadec
841 841
842 842 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
843 843
844 844 dataOut.flagDecodeData = True #asumo q la data esta decodificada
845 845
846 846 if self.__profIndex == self.nCode-1:
847 847 self.__profIndex = 0
848 848 return dataOut
849 849
850 850 self.__profIndex += 1
851 851
852 852 return dataOut
853 853 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
854 854
855 855
856 856 class ProfileConcat(Operation):
857 857
858 858 isConfig = False
859 859 buffer = None
860 860
861 861 def __init__(self, **kwargs):
862 862
863 863 Operation.__init__(self, **kwargs)
864 864 self.profileIndex = 0
865 865
866 866 def reset(self):
867 867 self.buffer = numpy.zeros_like(self.buffer)
868 868 self.start_index = 0
869 869 self.times = 1
870 870
871 871 def setup(self, data, m, n=1):
872 872 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
873 873 self.nHeights = data.shape[1]#.nHeights
874 874 self.start_index = 0
875 875 self.times = 1
876 876
877 877 def concat(self, data):
878 878
879 879 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
880 880 self.start_index = self.start_index + self.nHeights
881 881
882 882 def run(self, dataOut, m):
883 883 dataOut.flagNoData = True
884 884
885 885 if not self.isConfig:
886 886 self.setup(dataOut.data, m, 1)
887 887 self.isConfig = True
888 888
889 889 if dataOut.flagDataAsBlock:
890 890 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
891 891
892 892 else:
893 893 self.concat(dataOut.data)
894 894 self.times += 1
895 895 if self.times > m:
896 896 dataOut.data = self.buffer
897 897 self.reset()
898 898 dataOut.flagNoData = False
899 899 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
900 900 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
901 901 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
902 902 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
903 903 dataOut.ippSeconds *= m
904 904 return dataOut
905 905
906 906 class ProfileSelector(Operation):
907 907
908 908 profileIndex = None
909 909 # Tamanho total de los perfiles
910 910 nProfiles = None
911 911
912 912 def __init__(self, **kwargs):
913 913
914 914 Operation.__init__(self, **kwargs)
915 915 self.profileIndex = 0
916 916
917 917 def incProfileIndex(self):
918 918
919 919 self.profileIndex += 1
920 920
921 921 if self.profileIndex >= self.nProfiles:
922 922 self.profileIndex = 0
923 923
924 924 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
925 925
926 926 if profileIndex < minIndex:
927 927 return False
928 928
929 929 if profileIndex > maxIndex:
930 930 return False
931 931
932 932 return True
933 933
934 934 def isThisProfileInList(self, profileIndex, profileList):
935 935
936 936 if profileIndex not in profileList:
937 937 return False
938 938
939 939 return True
940 940
941 941 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
942 942
943 943 """
944 944 ProfileSelector:
945 945
946 946 Inputs:
947 947 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
948 948
949 949 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
950 950
951 951 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
952 952
953 953 """
954 954
955 955 if rangeList is not None:
956 956 if type(rangeList[0]) not in (tuple, list):
957 957 rangeList = [rangeList]
958 958
959 959 dataOut.flagNoData = True
960 960
961 961 if dataOut.flagDataAsBlock:
962 962 """
963 963 data dimension = [nChannels, nProfiles, nHeis]
964 964 """
965 965 if profileList != None:
966 966 dataOut.data = dataOut.data[:,profileList,:]
967 967
968 968 if profileRangeList != None:
969 969 minIndex = profileRangeList[0]
970 970 maxIndex = profileRangeList[1]
971 971 profileList = list(range(minIndex, maxIndex+1))
972 972
973 973 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
974 974
975 975 if rangeList != None:
976 976
977 977 profileList = []
978 978
979 979 for thisRange in rangeList:
980 980 minIndex = thisRange[0]
981 981 maxIndex = thisRange[1]
982 982
983 983 profileList.extend(list(range(minIndex, maxIndex+1)))
984 984
985 985 dataOut.data = dataOut.data[:,profileList,:]
986 986
987 987 dataOut.nProfiles = len(profileList)
988 988 dataOut.profileIndex = dataOut.nProfiles - 1
989 989 dataOut.flagNoData = False
990 990
991 991 return dataOut
992 992
993 993 """
994 994 data dimension = [nChannels, nHeis]
995 995 """
996 996
997 997 if profileList != None:
998 998
999 999 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1000 1000
1001 1001 self.nProfiles = len(profileList)
1002 1002 dataOut.nProfiles = self.nProfiles
1003 1003 dataOut.profileIndex = self.profileIndex
1004 1004 dataOut.flagNoData = False
1005 1005
1006 1006 self.incProfileIndex()
1007 1007 return dataOut
1008 1008
1009 1009 if profileRangeList != None:
1010 1010
1011 1011 minIndex = profileRangeList[0]
1012 1012 maxIndex = profileRangeList[1]
1013 1013
1014 1014 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1015 1015
1016 1016 self.nProfiles = maxIndex - minIndex + 1
1017 1017 dataOut.nProfiles = self.nProfiles
1018 1018 dataOut.profileIndex = self.profileIndex
1019 1019 dataOut.flagNoData = False
1020 1020
1021 1021 self.incProfileIndex()
1022 1022 return dataOut
1023 1023
1024 1024 if rangeList != None:
1025 1025
1026 1026 nProfiles = 0
1027 1027
1028 1028 for thisRange in rangeList:
1029 1029 minIndex = thisRange[0]
1030 1030 maxIndex = thisRange[1]
1031 1031
1032 1032 nProfiles += maxIndex - minIndex + 1
1033 1033
1034 1034 for thisRange in rangeList:
1035 1035
1036 1036 minIndex = thisRange[0]
1037 1037 maxIndex = thisRange[1]
1038 1038
1039 1039 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1040 1040
1041 1041 self.nProfiles = nProfiles
1042 1042 dataOut.nProfiles = self.nProfiles
1043 1043 dataOut.profileIndex = self.profileIndex
1044 1044 dataOut.flagNoData = False
1045 1045
1046 1046 self.incProfileIndex()
1047 1047
1048 1048 break
1049 1049
1050 1050 return dataOut
1051 1051
1052 1052
1053 1053 if beam != None: #beam is only for AMISR data
1054 1054 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1055 1055 dataOut.flagNoData = False
1056 1056 dataOut.profileIndex = self.profileIndex
1057 1057
1058 1058 self.incProfileIndex()
1059 1059
1060 1060 return dataOut
1061 1061
1062 1062 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1063 1063
1064 1064
1065 1065 class Reshaper(Operation):
1066 1066
1067 1067 def __init__(self, **kwargs):
1068 1068
1069 1069 Operation.__init__(self, **kwargs)
1070 1070
1071 1071 self.__buffer = None
1072 1072 self.__nitems = 0
1073 1073
1074 1074 def __appendProfile(self, dataOut, nTxs):
1075 1075
1076 1076 if self.__buffer is None:
1077 1077 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1078 1078 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1079 1079
1080 1080 ini = dataOut.nHeights * self.__nitems
1081 1081 end = ini + dataOut.nHeights
1082 1082
1083 1083 self.__buffer[:, ini:end] = dataOut.data
1084 1084
1085 1085 self.__nitems += 1
1086 1086
1087 1087 return int(self.__nitems*nTxs)
1088 1088
1089 1089 def __getBuffer(self):
1090 1090
1091 1091 if self.__nitems == int(1./self.__nTxs):
1092 1092
1093 1093 self.__nitems = 0
1094 1094
1095 1095 return self.__buffer.copy()
1096 1096
1097 1097 return None
1098 1098
1099 1099 def __checkInputs(self, dataOut, shape, nTxs):
1100 1100
1101 1101 if shape is None and nTxs is None:
1102 1102 raise ValueError("Reshaper: shape of factor should be defined")
1103 1103
1104 1104 if nTxs:
1105 1105 if nTxs < 0:
1106 1106 raise ValueError("nTxs should be greater than 0")
1107 1107
1108 1108 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1109 1109 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1110 1110
1111 1111 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1112 1112
1113 1113 return shape, nTxs
1114 1114
1115 1115 if len(shape) != 2 and len(shape) != 3:
1116 1116 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1117 1117
1118 1118 if len(shape) == 2:
1119 1119 shape_tuple = [dataOut.nChannels]
1120 1120 shape_tuple.extend(shape)
1121 1121 else:
1122 1122 shape_tuple = list(shape)
1123 1123
1124 1124 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1125 1125
1126 1126 return shape_tuple, nTxs
1127 1127
1128 1128 def run(self, dataOut, shape=None, nTxs=None):
1129 1129
1130 1130 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1131 1131
1132 1132 dataOut.flagNoData = True
1133 1133 profileIndex = None
1134 1134
1135 1135 if dataOut.flagDataAsBlock:
1136 1136
1137 1137 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1138 1138 dataOut.flagNoData = False
1139 1139
1140 1140 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1141 1141
1142 1142 else:
1143 1143
1144 1144 if self.__nTxs < 1:
1145 1145
1146 1146 self.__appendProfile(dataOut, self.__nTxs)
1147 1147 new_data = self.__getBuffer()
1148 1148
1149 1149 if new_data is not None:
1150 1150 dataOut.data = new_data
1151 1151 dataOut.flagNoData = False
1152 1152
1153 1153 profileIndex = dataOut.profileIndex*nTxs
1154 1154
1155 1155 else:
1156 1156 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1157 1157
1158 1158 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1159 1159
1160 1160 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1161 1161
1162 1162 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1163 1163
1164 1164 dataOut.profileIndex = profileIndex
1165 1165
1166 1166 dataOut.ippSeconds /= self.__nTxs
1167 1167
1168 1168 return dataOut
1169 1169
1170 1170 class SplitProfiles(Operation):
1171 1171
1172 1172 def __init__(self, **kwargs):
1173 1173
1174 1174 Operation.__init__(self, **kwargs)
1175 1175
1176 1176 def run(self, dataOut, n):
1177 1177
1178 1178 dataOut.flagNoData = True
1179 1179 profileIndex = None
1180 1180
1181 1181 if dataOut.flagDataAsBlock:
1182 1182
1183 1183 #nchannels, nprofiles, nsamples
1184 1184 shape = dataOut.data.shape
1185 1185
1186 1186 if shape[2] % n != 0:
1187 1187 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1188 1188
1189 1189 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1190 1190
1191 1191 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1192 1192 dataOut.flagNoData = False
1193 1193
1194 1194 profileIndex = int(dataOut.nProfiles/n) - 1
1195 1195
1196 1196 else:
1197 1197
1198 1198 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1199 1199
1200 1200 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1201 1201
1202 1202 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1203 1203
1204 1204 dataOut.nProfiles = int(dataOut.nProfiles*n)
1205 1205
1206 1206 dataOut.profileIndex = profileIndex
1207 1207
1208 1208 dataOut.ippSeconds /= n
1209 1209
1210 1210 return dataOut
1211 1211
1212 1212 class CombineProfiles(Operation):
1213 1213 def __init__(self, **kwargs):
1214 1214
1215 1215 Operation.__init__(self, **kwargs)
1216 1216
1217 1217 self.__remData = None
1218 1218 self.__profileIndex = 0
1219 1219
1220 1220 def run(self, dataOut, n):
1221 1221
1222 1222 dataOut.flagNoData = True
1223 1223 profileIndex = None
1224 1224
1225 1225 if dataOut.flagDataAsBlock:
1226 1226
1227 1227 #nchannels, nprofiles, nsamples
1228 1228 shape = dataOut.data.shape
1229 1229 new_shape = shape[0], shape[1]/n, shape[2]*n
1230 1230
1231 1231 if shape[1] % n != 0:
1232 1232 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1233 1233
1234 1234 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1235 1235 dataOut.flagNoData = False
1236 1236
1237 1237 profileIndex = int(dataOut.nProfiles*n) - 1
1238 1238
1239 1239 else:
1240 1240
1241 1241 #nchannels, nsamples
1242 1242 if self.__remData is None:
1243 1243 newData = dataOut.data
1244 1244 else:
1245 1245 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1246 1246
1247 1247 self.__profileIndex += 1
1248 1248
1249 1249 if self.__profileIndex < n:
1250 1250 self.__remData = newData
1251 1251 #continue
1252 1252 return
1253 1253
1254 1254 self.__profileIndex = 0
1255 1255 self.__remData = None
1256 1256
1257 1257 dataOut.data = newData
1258 1258 dataOut.flagNoData = False
1259 1259
1260 1260 profileIndex = dataOut.profileIndex/n
1261 1261
1262 1262
1263 1263 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1264 1264
1265 1265 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1266 1266
1267 1267 dataOut.nProfiles = int(dataOut.nProfiles/n)
1268 1268
1269 1269 dataOut.profileIndex = profileIndex
1270 1270
1271 1271 dataOut.ippSeconds *= n
1272 1272
1273 1273 return dataOut
1274 1274
1275 1275 class PulsePairVoltage(Operation):
1276 1276 '''
1277 1277 Function PulsePair(Signal Power, Velocity)
1278 1278 The real component of Lag[0] provides Intensity Information
1279 1279 The imag component of Lag[1] Phase provides Velocity Information
1280 1280
1281 1281 Configuration Parameters:
1282 1282 nPRF = Number of Several PRF
1283 1283 theta = Degree Azimuth angel Boundaries
1284 1284
1285 1285 Input:
1286 1286 self.dataOut
1287 1287 lag[N]
1288 1288 Affected:
1289 1289 self.dataOut.spc
1290 1290 '''
1291 1291 isConfig = False
1292 1292 __profIndex = 0
1293 1293 __initime = None
1294 1294 __lastdatatime = None
1295 1295 __buffer = None
1296 1296 noise = None
1297 1297 __dataReady = False
1298 1298 n = None
1299 1299 __nch = 0
1300 1300 __nHeis = 0
1301 1301 removeDC = False
1302 1302 ipp = None
1303 1303 lambda_ = 0
1304 1304
1305 1305 def __init__(self,**kwargs):
1306 1306 Operation.__init__(self,**kwargs)
1307 1307
1308 1308 def setup(self, dataOut, n = None, removeDC=False):
1309 1309 '''
1310 1310 n= Numero de PRF's de entrada
1311 1311 '''
1312 1312 self.__initime = None
1313 1313 self.__lastdatatime = 0
1314 1314 self.__dataReady = False
1315 1315 self.__buffer = 0
1316 1316 self.__profIndex = 0
1317 1317 self.noise = None
1318 1318 self.__nch = dataOut.nChannels
1319 1319 self.__nHeis = dataOut.nHeights
1320 1320 self.removeDC = removeDC
1321 1321 self.lambda_ = 3.0e8/(9345.0e6)
1322 1322 self.ippSec = dataOut.ippSeconds
1323 1323 self.nCohInt = dataOut.nCohInt
1324 1324 print("IPPseconds",dataOut.ippSeconds)
1325 1325
1326 1326 print("ELVALOR DE n es:", n)
1327 1327 if n == None:
1328 1328 raise ValueError("n should be specified.")
1329 1329
1330 1330 if n != None:
1331 1331 if n<2:
1332 1332 raise ValueError("n should be greater than 2")
1333 1333
1334 1334 self.n = n
1335 1335 self.__nProf = n
1336 1336
1337 1337 self.__buffer = numpy.zeros((dataOut.nChannels,
1338 1338 n,
1339 1339 dataOut.nHeights),
1340 1340 dtype='complex')
1341 self.noise = numpy.zeros([self.__nch,self.__nHeis])
1342 for i in range(self.__nch):
1343 self.noise[i]=dataOut.getNoise(channel=i)
1341 #self.noise = numpy.zeros([self.__nch,self.__nHeis])
1342 #for i in range(self.__nch):
1343 # self.noise[i]=dataOut.getNoise(channel=i)
1344 1344
1345 1345 def putData(self,data):
1346 1346 '''
1347 1347 Add a profile to he __buffer and increase in one the __profiel Index
1348 1348 '''
1349 1349 self.__buffer[:,self.__profIndex,:]= data
1350 1350 self.__profIndex += 1
1351 1351 return
1352 1352
1353 def pushData(self):
1353 def pushData(self,dataOut):
1354 1354 '''
1355 1355 Return the PULSEPAIR and the profiles used in the operation
1356 1356 Affected : self.__profileIndex
1357 1357 '''
1358
1359 1358 if self.removeDC==True:
1360 1359 mean = numpy.mean(self.__buffer,1)
1361 1360 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1362 1361 dc= numpy.tile(tmp,[1,self.__nProf,1])
1363 1362 self.__buffer = self.__buffer - dc
1364 1363
1365 1364 lag_0 = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1)
1366 1365 data_intensity = lag_0/(self.n*self.nCohInt)#*self.nCohInt)
1367 1366
1368 1367 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1369 1368 lag_1 = numpy.sum(pair1,1)
1370 1369 #angle = numpy.angle(numpy.sum(pair1,1))*180/(math.pi)
1371 1370 data_velocity = (-1.0*self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(lag_1)#self.ippSec*self.nCohInt
1372 1371
1372 self.noise = numpy.zeros([self.__nch,self.__nHeis])
1373 for i in range(self.__nch):
1374 self.noise[i]=dataOut.getNoise(channel=i)
1375
1373 1376 lag_0 = lag_0.real/(self.n)
1374 1377 lag_1 = lag_1/(self.n-1)
1375 1378 R1 = numpy.abs(lag_1)
1376 1379 S = (lag_0-self.noise)
1377 #k = R1/S
1378 #k = 1-k
1379 #k =numpy.absolute(k)
1380 #k =numpy.sqrt(k)
1380
1381 data_snrPP = S/self.noise
1382 data_snrPP = numpy.where(data_snrPP<0,1,data_snrPP)
1383
1381 1384 L = S/R1
1382 #print("L",L[0])
1383 1385 L = numpy.where(L<0,1,L)
1384 1386 L = numpy.log(L)
1387
1385 1388 tmp = numpy.sqrt(numpy.absolute(L))
1389
1386 1390 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*tmp*numpy.sign(L)
1387 1391 #data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*k
1388 1392 n = self.__profIndex
1389 1393
1390 1394 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1391 1395 self.__profIndex = 0
1392 return data_intensity,data_velocity,data_specwidth,n
1396 return data_intensity,data_velocity,data_snrPP,data_specwidth,n
1393 1397
1394 def pulsePairbyProfiles(self,data):
1398 def pulsePairbyProfiles(self,dataOut):
1395 1399
1396 1400 self.__dataReady = False
1397 1401 data_intensity = None
1398 1402 data_velocity = None
1399 1403 data_specwidth = None
1400 self.putData(data)
1404 data_snrPP = None
1405 self.putData(data=dataOut.data)
1401 1406 if self.__profIndex == self.n:
1402 1407 #self.noise = numpy.zeros([self.__nch,self.__nHeis])
1403 1408 #for i in range(self.__nch):
1404 1409 # self.noise[i]=data.getNoise(channel=i)
1405 1410 #print(self.noise.shape)
1406 data_intensity, data_velocity,data_specwidth, n = self.pushData()
1411 data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1407 1412 self.__dataReady = True
1408 1413
1409 return data_intensity, data_velocity,data_specwidth
1414 return data_intensity, data_velocity,data_snrPP,data_specwidth
1410 1415
1411 def pulsePairOp(self, data, datatime= None):
1416 def pulsePairOp(self, dataOut, datatime= None):
1412 1417
1413 1418 if self.__initime == None:
1414 1419 self.__initime = datatime
1415
1416 data_intensity, data_velocity,data_specwidth = self.pulsePairbyProfiles(data)
1420 #print("hola")
1421 data_intensity, data_velocity,data_snrPP,data_specwidth = self.pulsePairbyProfiles(dataOut)
1417 1422 self.__lastdatatime = datatime
1418 1423
1419 1424 if data_intensity is None:
1420 return None, None,None, None
1425 return None, None,None,None,None
1421 1426
1422 1427 avgdatatime = self.__initime
1423 1428 deltatime = datatime - self.__lastdatatime
1424 1429 self.__initime = datatime
1425 1430
1426 return data_intensity, data_velocity,data_specwidth,avgdatatime
1431 return data_intensity, data_velocity,data_snrPP,data_specwidth,avgdatatime
1427 1432
1428 1433 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1429 1434
1430 1435 if not self.isConfig:
1431 1436 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1432 1437 self.isConfig = True
1433 data_intensity, data_velocity,data_specwidth, avgdatatime = self.pulsePairOp(dataOut.data, dataOut.utctime)
1438 data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1434 1439 dataOut.flagNoData = True
1435 1440
1436 1441 if self.__dataReady:
1437 1442 dataOut.nCohInt *= self.n
1438 1443 dataOut.data_intensity = data_intensity #valor para intensidad
1439 1444 dataOut.data_velocity = data_velocity #valor para velocidad
1445 dataOut.data_snrPP = data_snrPP # valor para snr
1440 1446 dataOut.data_specwidth = data_specwidth
1441 1447 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1442 1448 dataOut.utctime = avgdatatime
1443 1449 dataOut.flagNoData = False
1444 1450 return dataOut
1445 1451
1446 1452
1447 1453 # import collections
1448 1454 # from scipy.stats import mode
1449 1455 #
1450 1456 # class Synchronize(Operation):
1451 1457 #
1452 1458 # isConfig = False
1453 1459 # __profIndex = 0
1454 1460 #
1455 1461 # def __init__(self, **kwargs):
1456 1462 #
1457 1463 # Operation.__init__(self, **kwargs)
1458 1464 # # self.isConfig = False
1459 1465 # self.__powBuffer = None
1460 1466 # self.__startIndex = 0
1461 1467 # self.__pulseFound = False
1462 1468 #
1463 1469 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1464 1470 #
1465 1471 # #Read data
1466 1472 #
1467 1473 # powerdB = dataOut.getPower(channel = channel)
1468 1474 # noisedB = dataOut.getNoise(channel = channel)[0]
1469 1475 #
1470 1476 # self.__powBuffer.extend(powerdB.flatten())
1471 1477 #
1472 1478 # dataArray = numpy.array(self.__powBuffer)
1473 1479 #
1474 1480 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1475 1481 #
1476 1482 # maxValue = numpy.nanmax(filteredPower)
1477 1483 #
1478 1484 # if maxValue < noisedB + 10:
1479 1485 # #No se encuentra ningun pulso de transmision
1480 1486 # return None
1481 1487 #
1482 1488 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1483 1489 #
1484 1490 # if len(maxValuesIndex) < 2:
1485 1491 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1486 1492 # return None
1487 1493 #
1488 1494 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1489 1495 #
1490 1496 # #Seleccionar solo valores con un espaciamiento de nSamples
1491 1497 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1492 1498 #
1493 1499 # if len(pulseIndex) < 2:
1494 1500 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1495 1501 # return None
1496 1502 #
1497 1503 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1498 1504 #
1499 1505 # #remover senales que se distancien menos de 10 unidades o muestras
1500 1506 # #(No deberian existir IPP menor a 10 unidades)
1501 1507 #
1502 1508 # realIndex = numpy.where(spacing > 10 )[0]
1503 1509 #
1504 1510 # if len(realIndex) < 2:
1505 1511 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1506 1512 # return None
1507 1513 #
1508 1514 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1509 1515 # realPulseIndex = pulseIndex[realIndex]
1510 1516 #
1511 1517 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1512 1518 #
1513 1519 # print "IPP = %d samples" %period
1514 1520 #
1515 1521 # self.__newNSamples = dataOut.nHeights #int(period)
1516 1522 # self.__startIndex = int(realPulseIndex[0])
1517 1523 #
1518 1524 # return 1
1519 1525 #
1520 1526 #
1521 1527 # def setup(self, nSamples, nChannels, buffer_size = 4):
1522 1528 #
1523 1529 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1524 1530 # maxlen = buffer_size*nSamples)
1525 1531 #
1526 1532 # bufferList = []
1527 1533 #
1528 1534 # for i in range(nChannels):
1529 1535 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1530 1536 # maxlen = buffer_size*nSamples)
1531 1537 #
1532 1538 # bufferList.append(bufferByChannel)
1533 1539 #
1534 1540 # self.__nSamples = nSamples
1535 1541 # self.__nChannels = nChannels
1536 1542 # self.__bufferList = bufferList
1537 1543 #
1538 1544 # def run(self, dataOut, channel = 0):
1539 1545 #
1540 1546 # if not self.isConfig:
1541 1547 # nSamples = dataOut.nHeights
1542 1548 # nChannels = dataOut.nChannels
1543 1549 # self.setup(nSamples, nChannels)
1544 1550 # self.isConfig = True
1545 1551 #
1546 1552 # #Append new data to internal buffer
1547 1553 # for thisChannel in range(self.__nChannels):
1548 1554 # bufferByChannel = self.__bufferList[thisChannel]
1549 1555 # bufferByChannel.extend(dataOut.data[thisChannel])
1550 1556 #
1551 1557 # if self.__pulseFound:
1552 1558 # self.__startIndex -= self.__nSamples
1553 1559 #
1554 1560 # #Finding Tx Pulse
1555 1561 # if not self.__pulseFound:
1556 1562 # indexFound = self.__findTxPulse(dataOut, channel)
1557 1563 #
1558 1564 # if indexFound == None:
1559 1565 # dataOut.flagNoData = True
1560 1566 # return
1561 1567 #
1562 1568 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1563 1569 # self.__pulseFound = True
1564 1570 # self.__startIndex = indexFound
1565 1571 #
1566 1572 # #If pulse was found ...
1567 1573 # for thisChannel in range(self.__nChannels):
1568 1574 # bufferByChannel = self.__bufferList[thisChannel]
1569 1575 # #print self.__startIndex
1570 1576 # x = numpy.array(bufferByChannel)
1571 1577 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1572 1578 #
1573 1579 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1574 1580 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1575 1581 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1576 1582 #
1577 1583 # dataOut.data = self.__arrayBuffer
1578 1584 #
1579 1585 # self.__startIndex += self.__newNSamples
1580 1586 #
1581 1587 # return
General Comments 0
You need to be logged in to leave comments. Login now