pylal.seriesutils

1 # Copyright (C) 2012 Duncan M. Macleod 2 # 3 # This program is free software; you can redistribute it and/or modify it 4 # under the terms of the GNU General Public License as published by the 5 # Free Software Foundation; either version 3 of the License, or (at your 6 # option) any later version. 7 # 8 # This program is distributed in the hope that it will be useful, but 9 # WITHOUT ANY WARRANTY; without even the implied warranty of 10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 # Public License for more details. 12 # 13 # You should have received a copy of the GNU General Public License along 14 # with this program; if not, write to the Free Software Foundation, Inc., 15 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 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 # Preamble 23 # ============================================================================= 24 25 from __future__ import division 26 import os,re,tempfile,warnings 27 28 # import swig LALSuite packages 29 from lal import lal 30 from lalframe import lalframe 31 32 # import gluecache 33 from glue import lal as _gluecache 34 35 # type code dict 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 # Data conversion 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 # Frame reading helpers 89 # ============================================================================= 90

91 -def fromNDS(chname, start, duration, server="nds.ligo.caltech.edu",\ 92 port=31200):

93 """ 94 Read data from NDS into a REAL?TimeSeries 95 """ 96 97 import nds 98 99 connection = nds.daq(server, port) 100 data = connection.fetch(start, start+duration, chname)[0] 101 epoch = lal.LIGOTimeGPS(start) 102 deltaT = duration/data.size 103 return fromarray(data, name=chname, epoch=epoch, deltaT=deltaT)

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 # open frame 131 framefile = os.path.abspath(framefile) 132 stream = lalframe.FrOpen('', framefile) 133 134 # return 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 # read cache from file object into gluecache for the time being 165 if hasattr(cache, 'readline'): 166 cache = gluecache.Cache.fromfile(cache) 167 168 # read cache from glue.lal.Cache, or lalframe.FrCache, or file 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 # open cache to stream 176 stream = lalframe.FrCacheOpen(cache) 177 178 # return 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 # set mode 209 if verbose: mode = lalframe.FR_STREAM_VERBOSE_MODE 210 else: mode = lalframe.FR_STREAM_DEFAULT_MODE 211 lalframe.FrSetMode(mode, stream) 212 213 # set time 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 # get series type 222 frdatatype = lalframe.FrStreamGetTimeSeriesType(chname, stream) 223 if datatype == -1: 224 datatype = frdatatype 225 226 TYPESTR = _typestr[datatype] 227 228 # get data 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 # return 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

373 -def duplicate(series):

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 # Average spectrum 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 # check data 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 # generate window 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 # generate FFT plan 459 if not plan: 460 func = getattr(lal, "CreateForward%sFFTPlan" % TYPESTR) 461 plan = func(seglen, 1) 462 463 # generate output spectrum 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 # calculate medianmean spectrum 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 # dereference 483 if destroywindow: 484 del window 485 if destroyplan: 486 del plan 487 488 # return 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 # get number of segments 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 # generate window 541 if not window: 542 func = getattr(lal, "CreateKaiser%sWindow" % TYPESTR) 543 window = func(seglen, 24) 544 545 # generate FFT plan 546 if not plan: 547 func = getattr(lal, "CreateForward%sFFTPlan" % TYPESTR) 548 plan = func(seglen, 1) 549 550 # set up return object 551 func = getattr(lal, "Create%sVectorSequence" % TYPESTR) 552 out = func(numseg, seglen//2+1) 553 554 # loop over steps 555 cut = getattr(lal, "Cut%sTimeSeries" % TYPESTR) 556 for i in range(numseg): 557 # extract portion of time series 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 # Get data type from series 568 # ============================================================================= 569

570 -def typecode(t):

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 # Helper functions 581 # ============================================================================= 582

583 -def gluecache_to_FrCache(cache):

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

Source Code for Module pylal.seriesutils