gstlal_psd_xml_from_asd_txt

#!/usr/bin/env python
#
# Copyright (C) 2011  Chad Hanna
#
# 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.

import numpy
from scipy.interpolate import interp1d
import sys
from gstlal import reference_psd
from pylal.xlal import datatypes
from pylal.xlal.datatypes import lalunit
from optparse import OptionParser

## @file
# gstlal_psd_xml_from_asd_txt
#
#
# A program to turn an ASCII PSD or ASD into an xml file suitable for consumption by various gstlal programs; see gstlal_psd_xml_from_asd_txt for details
#
# 
# ### Usage
#
# -# an early V1 input ASD from <a href=http://www.lsc-group.phys.uwm.edu/cgit/gstlal/plain/gstlal/share/v1_early_asd.txt>here</a>.
#
#       $ gstlal_psd_xml_from_asd_txt --instrument=V1 --output v1psd.xml.gz v1_early_asd.txt
#
#
# -# other things possible, please add some!
#
# ### Command line options
#
#   + `--instrument` [string]: Instrument, e.g. H1
#   + `--output` [filename]: Set the name of the LIGO light-weight XML file to output
#   + `--df` [float] (Hz): Set the frequency resolution to interpolate to, default = 0.25
#   + `--type` [ASD|PSD]: Input is ASD or PSD, default ASD
#
# specify the ASD or PSD text file as a file argument at the end
#
#
# ### Notes
#
# - The PSD will be interpolated with with linear interpolation if needed.  The user is highly encouraged to plot the result to see if it is satisfactory.  See, e.g., gstlal_plot_psd
#
#
# ### Review status
#
# | Reviewers               | git hash                  | date      |
# | ----------------------------------- | ---------------------             | ------------- |
# | Florent, Duncan Me., Jolien, Kipp, Chad |abedb7c19ad3d5348c000bc5f9435d5d2e810e42   | 2014-04-29    |
#
#

parser = OptionParser(description = __doc__)
parser.add_option("--instrument", help="instrument, e.g. H1")
parser.add_option("--output", metavar = "filename", help = "Set the name of the LIGO light-weight XML file to output")
parser.add_option("--df", metavar = "float", type = "float", default = 0.25, help = "set the frequency resolution to interpolate to, default = 0.25")
parser.add_option("--type", metavar = "ASD|PSD", default = "ASD", help = "input is ASD or PSD, default ASD")

options, filenames = parser.parse_args()

data = numpy.loadtxt(filenames[0], comments = '#')
f = data[:,0]
if options.type == "ASD":
    psd = data[:,1]**2
elif options.type == "PSD":
    psd = data[:,1]
else:
    raise ValueError("--type must be ASD or PSD")

#FIXME hack to pad the series since it doesn't start at 0 :(
f_padded = numpy.concatenate((numpy.arange(0, f[0], options.df), f))
psd_padded = numpy.concatenate((numpy.ones(len(f_padded) - len(f)) * psd[0], psd))

uniformf = numpy.arange(0, f_padded.max(), options.df)
psdinterp = interp1d(f_padded, psd_padded)

psd = psdinterp(uniformf)

psdseries = datatypes.real8frequencyseries.REAL8FrequencySeries(deltaF=options.df, data=psd, sampleUnits=lalunit.LALUnit("s strain^2"), f0=0)
reference_psd.write_psd(options.output, {options.instrument: psdseries})