#!/usr/bin/python
# Copyright (C) 2012 Alex Nitz, Tito Dal Canton
#
# 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 2 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.
"""Utilities to read PSDs from files.
"""
import logging
import numpy
import scipy.interpolate
from pycbc.types import FrequencySeries
logger = logging.getLogger('pycbc.psd.read')
[docs]
def from_numpy_arrays(freq_data, noise_data, length, delta_f, low_freq_cutoff):
"""Interpolate n PSD (as two 1-dimensional arrays of frequency and data)
to the desired length, delta_f and low frequency cutoff.
Parameters
----------
freq_data : array
Array of frequencies.
noise_data : array
PSD values corresponding to frequencies in freq_arr.
length : int
Length of the frequency series in samples.
delta_f : float
Frequency resolution of the frequency series in Herz.
low_freq_cutoff : float
Frequencies below this value are set to zero.
Returns
-------
psd : FrequencySeries
The generated frequency series.
"""
# Only include points above the low frequency cutoff
if freq_data[0] > low_freq_cutoff:
raise ValueError('Lowest frequency in input data '
' is higher than requested low-frequency cutoff ' + str(low_freq_cutoff))
kmin = int(low_freq_cutoff / delta_f)
flow = kmin * delta_f
data_start = (0 if freq_data[0]==low_freq_cutoff else numpy.searchsorted(freq_data, flow) - 1)
data_start = max(0, data_start)
# If the cutoff is exactly in the file, start there
if freq_data[data_start+1] == low_freq_cutoff:
data_start += 1
freq_data = freq_data[data_start:]
noise_data = noise_data[data_start:]
if (length - 1) * delta_f > freq_data[-1]:
logger.warning('Requested number of samples exceeds the highest '
'available frequency in the input data, '
'will use max available frequency instead. '
'(requested %f Hz, available %f Hz)',
(length - 1) * delta_f, freq_data[-1])
length = int(freq_data[-1]/delta_f + 1)
flog = numpy.log(freq_data)
slog = numpy.log(noise_data)
psd_interp = scipy.interpolate.interp1d(
flog, slog, fill_value=(slog[0], slog[-1]), bounds_error=False)
psd = numpy.zeros(length, dtype=numpy.float64)
vals = numpy.log(numpy.arange(kmin, length) * delta_f)
psd[kmin:] = numpy.exp(psd_interp(vals))
return FrequencySeries(psd, delta_f=delta_f)
[docs]
def from_txt(filename, length, delta_f, low_freq_cutoff, is_asd_file=True):
"""Read an ASCII file containing one-sided ASD or PSD data and generate
a frequency series with the corresponding PSD. The ASD or PSD data is
interpolated in order to match the desired resolution of the
generated frequency series.
Parameters
----------
filename : string
Path to a two-column ASCII file. The first column must contain
the frequency (positive frequencies only) and the second column
must contain the amplitude density OR power spectral density.
length : int
Length of the frequency series in samples.
delta_f : float
Frequency resolution of the frequency series in Herz.
low_freq_cutoff : float
Frequencies below this value are set to zero.
is_asd_file : Boolean
If false assume that the second column holds power spectral density.
If true assume that the second column holds amplitude spectral density.
Default: True
Returns
-------
psd : FrequencySeries
The generated frequency series.
Raises
------
ValueError
If the ASCII file contains negative, infinite or NaN frequencies
or amplitude densities.
"""
file_data = numpy.loadtxt(filename)
if (file_data < 0).any() or \
numpy.logical_not(numpy.isfinite(file_data)).any():
raise ValueError('Invalid data in ' + filename)
freq_data = file_data[:, 0]
noise_data = file_data[:, 1]
if is_asd_file:
noise_data = noise_data ** 2
return from_numpy_arrays(freq_data, noise_data, length, delta_f,
low_freq_cutoff)
[docs]
def from_xml(filename, length, delta_f, low_freq_cutoff, ifo_string=None,
root_name='psd'):
"""Read an ASCII file containing one-sided ASD or PSD data and generate
a frequency series with the corresponding PSD. The ASD or PSD data is
interpolated in order to match the desired resolution of the
generated frequency series.
Parameters
----------
filename : string
Path to a two-column ASCII file. The first column must contain
the frequency (positive frequencies only) and the second column
must contain the amplitude density OR power spectral density.
length : int
Length of the frequency series in samples.
delta_f : float
Frequency resolution of the frequency series in Herz.
low_freq_cutoff : float
Frequencies below this value are set to zero.
ifo_string : string
Use the PSD in the file's PSD dictionary with this ifo string.
If not given and only one PSD present in the file return that, if not
given and multiple (or zero) PSDs present an exception will be raised.
root_name : string (default='psd')
If given use this as the root name for the PSD XML file. If this means
nothing to you, then it is probably safe to ignore this option.
Returns
-------
psd : FrequencySeries
The generated frequency series.
"""
import lal.series
from ligo.lw import utils as ligolw_utils
with open(filename, 'rb') as fp:
ct_handler = lal.series.PSDContentHandler
xml_doc = ligolw_utils.load_fileobj(fp, compress='auto',
contenthandler=ct_handler)
psd_dict = lal.series.read_psd_xmldoc(xml_doc, root_name=root_name)
if ifo_string is not None:
psd_freq_series = psd_dict[ifo_string]
elif len(psd_dict.keys()) == 1:
psd_freq_series = psd_dict[tuple(psd_dict.keys())[0]]
else:
err_msg = "No ifo string given and input XML file contains not "
err_msg += "exactly one PSD. Specify which PSD you want to use."
raise ValueError(err_msg)
noise_data = psd_freq_series.data.data[:]
freq_data = numpy.arange(len(noise_data)) * psd_freq_series.deltaF
return from_numpy_arrays(freq_data, noise_data, length, delta_f,
low_freq_cutoff)