# Copyright (C) 2017 Hunter A. Gabbard, Andrew Lundgren,
# Duncan Macleod, Alex Nitz
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
"""
This module retrives a timeseries and then calculates
the q-transform of that time series
"""
import numpy
from numpy import ceil, log, exp
from pycbc.types.timeseries import FrequencySeries, TimeSeries
from pycbc.fft import ifft
from pycbc.types import zeros
[docs]def qplane(qplane_tile_dict, fseries, return_complex=False):
"""Performs q-transform on each tile for each q-plane and selects
tile with the maximum energy. Q-transform can then
be interpolated to a desired frequency and time resolution.
Parameters
----------
qplane_tile_dict:
Dictionary containing a list of q-tile tupples for each q-plane
fseries: 'pycbc FrequencySeries'
frequency-series data set
return_complex: {False, bool}
Return the raw complex series instead of the normalized power.
Returns
-------
q : float
The q of the maximum q plane
times : numpy.ndarray
The time that the qtransform is sampled.
freqs : numpy.ndarray
The frequencies that the qtransform is samled.
qplane : numpy.ndarray (2d)
The two dimensional interpolated qtransform of this time series.
"""
# store q-transforms for each q in a dict
qplanes = {}
max_energy, max_key = None, None
for i, q in enumerate(qplane_tile_dict):
energies = []
for f0 in qplane_tile_dict[q]:
energy = qseries(fseries, q, f0, return_complex=return_complex)
menergy = abs(energy).max()
energies.append(energy)
if i == 0 or menergy > max_energy:
max_energy = menergy
max_key = q
qplanes[q] = energies
# record q-transform output for peak q
plane = qplanes[max_key]
frequencies = qplane_tile_dict[max_key]
times = plane[0].sample_times.numpy()
plane = numpy.array([v.numpy() for v in plane])
return max_key, times, frequencies, numpy.array(plane)
[docs]def qtiling(fseries, qrange, frange, mismatch=0.2):
"""Iterable constructor of QTile tuples
Parameters
----------
fseries: 'pycbc FrequencySeries'
frequency-series data set
qrange:
upper and lower bounds of q range
frange:
upper and lower bounds of frequency range
mismatch:
percentage of desired fractional mismatch
Returns
-------
qplane_tile_dict: 'dict'
dictionary containing Q-tile tuples for a set of Q-planes
"""
qplane_tile_dict = {}
qs = list(_iter_qs(qrange, deltam_f(mismatch)))
for q in qs:
qtilefreq = _iter_frequencies(q, frange, mismatch, fseries.duration)
qplane_tile_dict[q] = numpy.array(list(qtilefreq))
return qplane_tile_dict
[docs]def deltam_f(mismatch):
"""Fractional mismatch between neighbouring tiles
Parameters
----------
mismatch: 'float'
percentage of desired fractional mismatch
Returns
-------
:type: 'float'
"""
return 2 * (mismatch / 3.) ** (1/2.)
def _iter_qs(qrange, deltam):
"""Iterate over the Q values
Parameters
----------
qrange:
upper and lower bounds of q range
deltam:
Fractional mismatch between neighbouring tiles
Returns
-------
Q-value:
Q value for Q-tile
"""
# work out how many Qs we need
cumum = log(float(qrange[1]) / qrange[0]) / 2**(1/2.)
nplanes = int(max(ceil(cumum / deltam), 1))
dq = cumum / nplanes
for i in range(nplanes):
yield qrange[0] * exp(2**(1/2.) * dq * (i + .5))
return
def _iter_frequencies(q, frange, mismatch, dur):
"""Iterate over the frequencies of this 'QPlane'
Parameters
----------
q:
q value
frange: 'list'
upper and lower bounds of frequency range
mismatch:
percentage of desired fractional mismatch
dur:
duration of timeseries in seconds
Returns
-------
frequencies:
Q-Tile frequency
"""
# work out how many frequencies we need
minf, maxf = frange
fcum_mismatch = log(float(maxf) / minf) * (2 + q**2)**(1/2.) / 2.
nfreq = int(max(1, ceil(fcum_mismatch / deltam_f(mismatch))))
fstep = fcum_mismatch / nfreq
fstepmin = 1. / dur
# for each frequency, yield a QTile
for i in range(nfreq):
yield (float(minf) *
exp(2 / (2 + q**2)**(1/2.) * (i + .5) * fstep) //
fstepmin * fstepmin)
return
[docs]def qseries(fseries, Q, f0, return_complex=False):
"""Calculate the energy 'TimeSeries' for the given fseries
Parameters
----------
fseries: 'pycbc FrequencySeries'
frequency-series data set
Q:
q value
f0:
central frequency
return_complex: {False, bool}
Return the raw complex series instead of the normalized power.
Returns
-------
energy: '~pycbc.types.TimeSeries'
A 'TimeSeries' of the normalized energy from the Q-transform of
this tile against the data.
"""
# normalize and generate bi-square window
qprime = Q / 11**(1/2.)
norm = numpy.sqrt(315. * qprime / (128. * f0))
window_size = 2 * int(f0 / qprime * fseries.duration) + 1
xfrequencies = numpy.linspace(-1., 1., window_size)
start = int((f0 - (f0 / qprime)) * fseries.duration)
end = int(start + window_size)
center = (start + end) // 2
windowed = fseries[start:end] * (1 - xfrequencies ** 2) ** 2 * norm
tlen = (len(fseries)-1) * 2
windowed.resize(tlen)
windowed.roll(-center)
# calculate the time series for this q -value
windowed = FrequencySeries(windowed, delta_f=fseries.delta_f,
epoch=fseries.start_time)
ctseries = TimeSeries(zeros(tlen, dtype=numpy.complex128),
delta_t=fseries.delta_t)
ifft(windowed, ctseries)
if return_complex:
return ctseries
else:
energy = ctseries.squared_norm()
medianenergy = numpy.median(energy.numpy())
return energy / float(medianenergy)