1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 """
18 This module provides a bunch of user-friendly wrappers to the SWIG-bound REAL8TimeSeries and REAL8FrequencySeries objects and their associated functions.
19 """
20
21
22
23
24
25 from __future__ import division
26 import os,re,tempfile,warnings
27
28
29 from lal import lal
30 from lalframe import lalframe
31
32
33 from glue import lal as _gluecache
34
35
36 _typestr = {\
37 lal.I2_TYPE_CODE: 'INT2',\
38 lal.I4_TYPE_CODE: 'INT4',\
39 lal.I8_TYPE_CODE: 'INT8',\
40 lal.U2_TYPE_CODE: 'UINT2',\
41 lal.U4_TYPE_CODE: 'UINT4',\
42 lal.U8_TYPE_CODE: 'UINT8',\
43 lal.S_TYPE_CODE: 'REAL4',\
44 lal.D_TYPE_CODE: 'REAL8',\
45 lal.C_TYPE_CODE: 'COMPLEX8',\
46 lal.Z_TYPE_CODE: 'COMPLEX16',\
47 }
48
49
50
51
52
53 -def fromarray(array, name="", epoch=lal.LIGOTimeGPS(), f0=0, deltaT=1,\
54 sampleUnits=lal.DimensionlessUnit, frequencyseries=False):
55 """
56 Convert numpy.array to REAL8TimeSeries. Use frequencyseries to return
57 REAL8FrequencySeries.
58
59 Arguments:
60
61 array : numpy.array
62 input data array
63
64 Keyword arguments:
65
66 name : str
67 name of data source
68 epoch : lal.LIGOTimeGPS
69 GPS start time for array
70 deltaT : float
71 sampling time for data (or frequency spacing for FrequencySeries)
72 f0 : float
73 lower frequency limit for data
74 sampleUnits : lal.Unit
75 amplitude unit for array
76 """
77
78 if frequencyseries:
79 series = lal.CreateREAL8FrequencySeries(name, epoch, f0, deltaT,\
80 sampleUnits, array.size)
81 else:
82 series = lal.CreateREAL8TimeSeries(name, epoch, f0, deltaT,\
83 sampleUnits, array.size)
84 series.data.data = array.astype(float)
85 return series
86
87
88
89
90
91 -def fromNDS(chname, start, duration, server="nds.ligo.caltech.edu",\
92 port=31200):
104
105 -def fromframefile(framefile, chname, start=-1, duration=1,\
106 datatype=-1, verbose=False):
107 """
108 Read data from the GWF format framefile into a REAL?TimeSeries.
109 Restrict data with the gpsstart and duration parameters.
110
111 Arguments:
112
113 framefile : str
114 path to input frame file
115 chname : str
116 name of channel to read
117
118 Keyword arguments:
119
120 start : lal.LIGOTimeGPS
121 GPS start time of requested data
122 duration : float
123 length of requested data series in seconds
124 datatype : int
125 LAL enum for requested datatype, -1 == read from frame metadata
126 verbose : [ True | False ]
127 verbose output
128 """
129
130
131 framefile = os.path.abspath(framefile)
132 stream = lalframe.FrOpen('', framefile)
133
134
135 series = fromFrStream(stream, chname, start=start, duration=duration,\
136 datatype=datatype, verbose=verbose)
137 return series
138
139 -def fromlalcache(cache, chname, start=-1, duration=1, datatype=-1,\
140 verbose=False):
141 """
142 Read data from a cache object into a REAL?TimeSeries.
143 Restrict data with the gpsstart and duration parameters.
144
145 Arguments:
146
147 cache : [ filepath | file | glue.lal.Cache | lalframe.FrCache ]
148 object to be read to lalframe.FrCache
149 chname : str
150 name of channel to read
151
152 Keyword arguments:
153
154 start : lal.LIGOTimeGPS
155 GPS start time of requested data
156 duration : float
157 length of requested data series in seconds
158 datatype : int
159 LAL enum for requested datatype, -1 == read from frame metadata
160 verbose : [ True | False ]
161 verbose output
162 """
163
164
165 if hasattr(cache, 'readline'):
166 cache = gluecache.Cache.fromfile(cache)
167
168
169 if cache and isinstance(cache, _gluecache.Cache):
170 cache.sort(key=lambda e: e.segment[0])
171 cache = gluecache_to_FrCache(cache)
172 elif isinstance(cache, str):
173 cache = lal.CacheImport(cache)
174
175
176 stream = lalframe.FrCacheOpen(cache)
177
178
179 series = fromFrStream(stream, chname, start=start, duration=duration,\
180 datatype=datatype, verbose=verbose)
181 return series
182
183 -def fromFrStream(stream, chname, start=-1, duration=1, datatype=-1,\
184 verbose=False):
185 """
186 Read data from the lalframe.LALFrStream object stream into a REAL?TimeSeries.
187 Restrict data with the gpsstart and duration parameters.
188
189 Arguments:
190
191 stream : lalframe.FrStream
192 frame stream to read
193 chname : str
194 name of channel to read
195
196 Keyword arguments:
197
198 start : lal.LIGOTimeGPS
199 GPS start time of requested data
200 duration : float
201 length of requested data series in seconds
202 datatype : int
203 LAL enum for requested datatype, -1 == read from frame metadata
204 verbose : [ True | False ]
205 verbose output
206 """
207
208
209 if verbose: mode = lalframe.FR_STREAM_VERBOSE_MODE
210 else: mode = lalframe.FR_STREAM_DEFAULT_MODE
211 lalframe.FrSetMode(mode, stream)
212
213
214 if int(start) == -1:
215 start = stream.epoch
216 else:
217 start = lal.LIGOTimeGPS(float(start))
218 duration = float(duration)
219 lalframe.FrSeek(start, stream)
220
221
222 frdatatype = lalframe.FrStreamGetTimeSeriesType(chname, stream)
223 if datatype == -1:
224 datatype = frdatatype
225
226 TYPESTR = _typestr[datatype]
227
228
229 if frdatatype == datatype:
230 func = getattr(lalframe, 'FrStreamRead%sTimeSeries' % TYPESTR)
231 series = func(stream, chname, start, duration, 0)
232 else:
233 dblseries = lalframe.FrStreamInputREAL8TimeSeries(stream, chname,
234 start, duration, 0)
235 func = getattr(lal, 'Create%sTimeSeries' % TYPESTR)
236 series = func(dblseries.name, dblseries.epoch, dblseries.f0,\
237 dblseries.deltaT, dblseries.sampleUnits,\
238 dblseries.data.length)
239 series.data.data = dblseries.data.data.astype(type(series.data.data[0]))
240 del dblseries
241
242
243 return series
244
245 -def resample(series, rate, inplace=True):
246 """
247 Resample a REAL?TimeSeries to the new rate (Hertz). By default resampling
248 is done in-place and a copy of the series is returned, but if requested
249 the original series is duplicated before resampling and so is unaffected.
250
251 Arguments:
252
253 series : swiglal.REAL?TimeSeries
254 input timeseries to resample
255 rate : float
256 sampling rate (Hertz) to resample to
257
258 Keyword arguments:
259
260 inplace : [ True | False ]
261 perform resampling inplace on input series, default: True.
262 If False, input series is duplicated and so is left in the original
263 state.
264 """
265
266 datatype = typecode(type(series))
267 TYPESTR = _typestr[datatype]
268
269 func = getattr(lal, 'Resample%sTimeSeries' % TYPESTR)
270 if inplace:
271 func(series, 1/rate)
272 return series
273 else:
274 series2 = duplicate(series)
275 func(series2, 1/rate)
276 return series2
277
278 -def highpass(series, frequency, amplitude=0.9, filtorder=8, inplace=True):
279 """
280 Apply Butterworth high pass filter to REAL?TimeSeries.
281
282 Arguments:
283
284 series : swiglal.REAL?TimeSeries
285 input series to highpass
286 frequency : float
287 frequency below which to attenuate
288
289 Keyword arguments:
290
291 amplitude : float
292 desired amplitude response of filter
293 filtorder : int
294 desired order of filter
295 inplace : [ True | False ]
296 perform resampling inplace on input series, default: True.
297 If False, input series is duplicated and so is left in the original
298 state.
299 """
300
301 datatype = typecode(type(series))
302 TYPESTR = _typestr[datatype]
303
304 func = getattr(lal, 'HighPass%sTimeSeries' % TYPESTR)
305 if inplace:
306 func(series, frequency, amplitude, filtorder)
307 return series
308 else:
309 series2 = duplicate(series)
310 func(series2, frequency, amplitude, filtorder)
311 return series2
312
313 -def lowpass(series, frequency, amplitude=0.9, filtorder=8, inplace=True):
314 """
315 Apply Butterworth low pass filter to REAL?TimeSeries.
316
317 Arguments:
318
319 series : swiglal.REAL?TimeSeries
320 input series to lowpass
321 frequency : float
322 frequency above which to attenuate
323
324 Keyword arguments:
325
326 amplitude : float
327 desired amplitude response of filter
328 filtorder : int
329 desired order of filter
330 inplace : [ True | False ]
331 perform resampling inplace on input series, default: True.
332 If False, input series is duplicated and so is left in the original
333 state.
334 """
335
336 datatype = typecode(type(series))
337 TYPESTR = _typestr[datatype]
338
339 func = getattr(lal, 'LowPass%sTimeSeries' % TYPESTR)
340 if inplace:
341 func(series, frequency, amplitude, filtorder)
342 return series
343 else:
344 series2 = duplicate(series)
345 func(series2, frequency, amplitude, filtorder)
346 return series2
347
348 -def bandpass(series, fmin, fmax, amplitude=0.9, filtorder=8, inplace=True):
349 """
350 Apply Butterworth band pass filter to REAL?TimeSeries.
351
352 Arguments:
353
354 series : swiglal.REAL?TimeSeries
355 input series to lowpass
356 fmin : float
357 frequency below which to attenuate
358 fmax : float
359 frequency above which to attenuate
360
361 Keyword arguments:
362
363 amplitude : float
364 desired amplitude response of filter
365 filtorder : int
366 desired order of filter
367 """
368
369 series = highpass(series, fmin, amplitude, filtorder, inplace)
370 series = lowpass(series, fmax, amplitude, filtorder, inplace)
371 return series
372
374 """
375 Duplicate a TimeSeries or FrequencySeries.
376
377 Arguments:
378
379 series : [ TimeSeries | FrequencySeries ]
380 input series to duplicate
381 """
382 seriestype = str(type(series)).split("'")[1]
383 datatype = typecode(seriestype)
384 TYPESTR = _typestr[datatype]
385 func = getattr(lal, 'Create%s' % seriestype)
386 if seriestype.endswith("FrequencySeries"):
387 out = func(series.name, series.epoch, series.f0, series.deltaF,\
388 series.sampleUnits, series.data.length)
389 elif seriestype.endswith("TimeSeries"):
390 out = func(series.name, series.epoch, series.f0, series.deltaT,\
391 series.sampleUnits, series.data.length)
392 out.data.data = series.data.data
393 return out
394
395
396
397
398
399 -def compute_average_spectrum(series, seglen, stride, window=None, plan=None,\
400 average='medianmean', unit=lal.StrainUnit):
401 """
402 Calculate the average (power) spectrum of the given REAL?TimeSeries
403
404 Arguments:
405
406 seglen : int
407 length of FFT (in samples)
408 stride : int
409 gap between FFTs (in samples)
410
411 Keyword arguments:
412
413 window : lal.REAL8Window
414 swiglal window
415 plan : lal.REAL8FFTPlan
416 plan for FFT
417 spectrum : [ 'median' | 'medianmean' | 'welch' ]
418 averaging method for spectrum, default: 'medianmean'
419 unit : lal.Unit
420 LAL unit for data
421 """
422
423 datatype = typecode(type(series))
424 TYPESTR = _typestr[datatype]
425
426 seglen = int(seglen)
427 stride = int(stride)
428
429
430 reclen = series.data.length
431 numseg = 1 + int((reclen - seglen)/stride)
432 worklen = (numseg - 1)*stride + seglen
433 if worklen != reclen:
434 warnings.warn("Data array is not long enough to be covered completely, "
435 "the trailing %d samples will not be used for spectrum "
436 "calculation" % (reclen-worklen))
437 series = duplicate(series)
438 func = getattr(lal, 'Resize%sTimeSeries' % TYPESTR)
439 func(series, 0, worklen)
440 reclen = series.data.length
441 numseg = 1 + int((reclen - seglen)/stride)
442 if numseg % 2 == 1:
443 warnings.warn("Data array is long enough for an odd number of FFT "
444 "averages. The last one will be discarded for the "
445 "median-mean spectrum.")
446 worklen = reclen - (seglen-stride)
447 series = duplicate(series)
448 func = getattr(lal, 'Resize%sTimeSeries' % TYPESTR)
449 series = func(series, 0, worklen)
450 reclen = series.data.length
451
452
453 destroywindow = not window
454 destroyplan = not plan
455 if not window:
456 func = getattr(lal, "CreateKaiser%sWindow" % TYPESTR)
457 window = func(seglen, 24)
458
459 if not plan:
460 func = getattr(lal, "CreateForward%sFFTPlan" % TYPESTR)
461 plan = func(seglen, 1)
462
463
464 f0 = (1/series.deltaT) * (1/seglen)
465 deltaF = (1/seglen)
466 func = getattr(lal, "Create%sFrequencySeries" % TYPESTR)
467 spectrum = func(series.name, series.epoch, f0, deltaF, lal.StrainUnit,\
468 seglen//2+1)
469
470
471 if re.match('median?mean\Z', average, re.I):
472 func = getattr(lal, "%sAverageSpectrumMedianMean" % TYPESTR)
473 elif re.match('median\Z', average, re.I):
474 func = getattr(lal, "%sAverageSpectrumMedian" % TYPESTR)
475 elif re.match('welch\Z', average, re.I):
476 func = getattr(lal, "%sAverageSpectrumWelch" % TYPESTR)
477 else:
478 raise NotImplementedError("Sorry, only 'median' and 'medianmean' "+\
479 " and 'welch' average methods are available.")
480 func(spectrum, series, seglen, stride, window, plan)
481
482
483 if destroywindow:
484 del window
485 if destroyplan:
486 del plan
487
488
489 return spectrum
490
491 -def compute_average_spectrogram(series, step, seglen, stride, window=None,\
492 plan=None, average='medianmean',\
493 unit=lal.StrainUnit):
494 """
495 Compute the average (power) spectrogram of the given REAL?TimeSeries by
496 stacking together average spectra for each timestep.
497
498 Arguments:
499
500 series : REAL8TimeSeries
501 input data
502 step : int
503 length of single average spectrum (in samples)
504 seglen : int
505 length of FFT (in samples)
506 stride : int
507 gap between FFTs (in samples)
508
509 Keyword arguments:
510
511 window : lal.REAL8Window
512 swiglal window
513 plan : lal.REAL8FFTPlan
514 plan for FFT
515 spectrum : [ 'median' | 'medianmean' | 'welch' ]
516 averaging method for spectrum, default: 'medianmean'
517 unit : lal.Unit
518 LAL unit for data
519 """
520
521 datatype = typecode(type(series))
522 TYPESTR = _typestr[datatype]
523
524 step = int(step)
525 seglen = int(seglen)
526 stride = int(stride)
527
528
529 duration = series.data.length
530 numseg = int(duration//step)
531 if numseg == 0:
532 raise ValueError("Data array is too short to compute even a single "
533 "average.")
534 if duration % step != 0:
535 warnings.warn("data is not the right size for complete coverage in %d "\
536 "point steps. %d steps will be computed and the "\
537 "remaining %d samples will be discarded."\
538 % (step, numseg, duration % step))
539
540
541 if not window:
542 func = getattr(lal, "CreateKaiser%sWindow" % TYPESTR)
543 window = func(seglen, 24)
544
545
546 if not plan:
547 func = getattr(lal, "CreateForward%sFFTPlan" % TYPESTR)
548 plan = func(seglen, 1)
549
550
551 func = getattr(lal, "Create%sVectorSequence" % TYPESTR)
552 out = func(numseg, seglen//2+1)
553
554
555 cut = getattr(lal, "Cut%sTimeSeries" % TYPESTR)
556 for i in range(numseg):
557
558 stepseries = cut(series, i*step, step)
559 spectrum = compute_average_spectrum(stepseries, seglen, stride,\
560 window=window, plan=plan,\
561 average=average)
562 out.data[i,:] = spectrum.data.data
563
564 return out, spectrum.deltaF, spectrum.f0
565
566
567
568
569
571 """
572 Return LAL type code for this type string.
573 """
574 for datatype,key in _typestr.iteritems():
575 if re.search(key, str(t)):
576 return datatype
577 raise TypeError("Cannot find type code for %s" % str(t))
578
579
580
581
582
584 """
585 Convert a glue.lal.Cache object to a lalframe.FrCache object.
586 Writes cache to temporary file and reads to FrCache.
587
588 Arguments:
589
590 cache : glue.lal.Cache
591 LAL cache object from GLUE to convert
592 """
593 with tempfile.NamedTemporaryFile(delete=False) as t:
594 cache = cache
595 for e in cache:
596 e.segment = type(e.segment)(int(e.segment[0]), int(e.segment[1]))
597 cache.tofile(t)
598 frcache = lal.CacheImport(t.name)
599 os.remove(t.name)
600 return frcache
601