Commit 4237bb05 authored by no_author's avatar no_author
Browse files

This commit was manufactured by cvs2svn to create tag 'v6r26'.

git-svn-id: https://svn.ego-gw.it/svn/advsw/Fr/tags/v6r26@22390 47e613b6-93f3-4244-b775-e06b1756610f
parent 64524fde
version: v6r03 date: 2003 02 26 16:11:38 user: swmgr reason: Import of v6r03
version: v6r06 date: 2003 05 20 18:01:56 user: swmgr reason: Frame Library, import of v6r06
version: v6r07 date: 2003 06 27 09:23:24 user: swmgr reason: import of v6r07
version: v6r07p1 currently modified by masserot previous: v6r07
version: v6r08 date: 2003 09 01 09:58:33 user: swmgr reason: See end of FrDoc for ChangeLog, fixed SPR #22 previous: v6r07
version: v6r09 date: 2003 10 02 10:36:25 user: swmgr reason: Fixed SPR #36 previous: v6r08
version: v6r10 date: 2003 12 15 10:21:06 user: swmgr reason: See end of FrDoc for ChangeLog, update to root v3.10.01_1, added saveadc and saveproc on octave
version: v6r11 date: 2004 01 23 15:09:03 user: swmgr reason: See end of FrDoc for ChangeLog, update to root v3.10.02
version: v6r12 date: 2004 05 28 18:05:47 user: swmgr reason: Used for VCS-2.0, see end of FrDoc for changes
version: v6r14 date: 2004 08 09 12:07:10 user: swmgr reason: See end of FrDoc
version: v6r15 date: 2004 09 14 09:52:54 user: swmgr reason: Bug fix - FrVectConcat and octave interface fix. See User manual for detail. previous: v6r14
version: v6r16 date: 2005 03 10 17:09:20 user: swmgr reason: See User manual for Bug fixes and New features
version: v6r17 date: 2005 03 10 17:11:21 user: swmgr reason: Allow the non alignement on GPS for multiple output files, See User manual for detail
version: v6r19 date: 2005 05 16 18:30:10 user: swmgr reason: See end of FrDoc
Frame Library Software Terms and Conditions
The authors hereby grant permission to use, copy, and distribute this
software and its documentation for any purpose, provided that existing
copyright notices are retained in all copies and that this notice is
included verbatim in any distributions. Additionally, the authors grant
permission to modify this software and its documentation for any purpose,
provided that such modifications are not distributed without the explicit
consent of the authors and that existing copyright notices are retained in
all copies. Users of the software are asked to feed back problems, benefits,
and/or suggestions about the authors.
Support for this software - fixing of bugs, incorporation of new features -
is done on a best effort basis. All bug fixes and enhancements will be made
available under the same terms and conditions as the original software,
IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR
DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF,
EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,
INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT. THIS SOFTWARE IS
PROVIDED ON AN "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO
OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
MODIFICATIONS.
Python frgetvect
----------------
This module is intended to allow access to GWF frame files in Python. It
interfaces the FrameL library.
Written by Nick Fotopoulos (nvf@ligo.mit.edu)
Installation instructions
-------------------------
1. Compile and install FrameL. Make sure that the libraries are in your
LD_LIBRARY_PATH.
2. Place the pyfrgetvect source files in a subdirectory of the main FrameL
directory. e.g., v6r20/pyfrgetvect
3. python setup.py build
4. sudo python setup.py install
Due to sudo's environment handling, in some environments, I have had to use
an alternate installation step:
4'. sudo su
5'. python setup.py install
6'. <Ctrl>-d
Usage
-----
The input arguments are:
1) file name (single file only)
2) channel name (it could be an ADC, SIM or PROC channel)
3) (optional) starting GPS time (default = first frame in the file)
4) (optional) vector length in second (default = 1 second)
5) (optional) debug level (default = 0 = no output)
Returned data (in a tuple):
1) ADC or SIM or PROC data stored as double
2) x axis values relative to the first data point.
This is usual time but it could be frequency in the case
of a frequency series
3) frequency values in the case of time series
4) GPS starting time (in seconds)
5) starting time as a string
6) vector unitX as a string
Tested configurations
---------------------
Tested on Debian GNU/Linux (x86), Red Hat Linux (x86,x86_64) and Mac OS X (PPC).
Code is written for handling complex data, but that has not been tested.
/*
Nick Fotopoulos (nvf@ligo.mit.edu)
Python port of frgetvect (based on Matlab frgetvect by Benoit Mours)
$Id: frgetvect.c,v 1.1.1.5 2009-03-17 08:38:50 swmgr Exp $
Requires: numpy, FrameL
To do: Check inputs, allow multiple files, raise appropriate exceptions,
test on complex-valued frames, make more Pythonic (use keyword
arguments, etc)
*/
const char docstring[] = "Nick Fotopoulos (nvf@mit.edu)\n"
"Python port of frgetvect (based on Matlab frgetvect)\n"
"\n"
"The input arguments are:\n"
" 1) file name\n"
" 2) channel name (it could be an ADC, SIM or PROC channel)\n"
" 3) (optional) starting GPS time (default = first frame in the file)\n"
" 4) (optional) vector length in second (default = 1 second)\n"
" 5) (optional) debug level (default = 0 = no output)\n"
"\n"
"Returned data (in a tuple):\n"
" 1) ADC or SIM or PROC data stored as double\n"
" 2) x axis values relative to the first data point. This is usually\n"
" time, but it could be frequency in the case of a frequency series\n"
" 3) frequency values, in the case of time series\n"
" 4) GPS starting time (in seconds)\n"
" 5) starting time as a string\n"
" 6) vector unitX as a string\n"
" 7) vector unitY as a string\n";
#include "Python.h"
#include "numpy/arrayobject.h"
#include "FrameL.h"
static PyObject *frgetvect(PyObject *self, PyObject *args)
{
import_array();
int nrhs = 5; // I don't actually use this at present.
PyArrayObject *out1, *out2, *out3;
double out4;
char *out5, *out6, *out7;
struct FrFile *iFile;
struct FrVect *vect;
long debugLvl, i, nData, utc;
double *data, dt, df, t0, duration;
char *fileName, *vectName, *tempstr;
npy_intp shape[1];
int ok;
/*--------------- unpack arguments --------------------*/
/* Give bad default arguments as a way of counting how many
arguments were passed to this function. */
t0 = -1.;
duration = -1.;
debugLvl = -1;
/* The | in the format string indicates the next arguments are
optional. They are simply not assigned anything. */
ok = PyArg_ParseTuple(args, "ss|ddi", &fileName, &vectName,
&t0, &duration, &debugLvl);
// Handle defaults
if (debugLvl < 0)
{
debugLvl = 0;
nrhs = 4;
}
if (duration < 0.)
{
duration = 1.;
nrhs = 3;
}
if (t0 < 0.)
{
t0 = 0.;
nrhs = 2;
}
FrLibSetLvl(debugLvl);
/*-------------- open file --------------------------*/
if (debugLvl > 0)
{
const char *msg = "Opening %s for channel %s (t0=%.2f, duration=%.2f).\n";
printf(msg, fileName, vectName, t0, duration);
}
iFile = FrFileINew(fileName);
if (iFile == NULL)
{
printf("%s\n", FrErrorGetHistory());
return;
}
if (debugLvl > 0)
{
const char *msg = "Opened %s for channel %s!\n";
printf(msg, fileName, vectName);
}
/*-------------- get vector --------------------------*/
vect = FrFileIGetV(iFile, vectName, t0, duration);
if(debugLvl > 0) FrVectDump(vect, stdout, debugLvl);
if(vect == NULL)
{
printf("In file %s, vector not found: %s\n",fileName,vectName);
FrFileIEnd(iFile);
return;
}
if (debugLvl > 0)
{
const char *msg = "Extracted channel %s successfully!\n";
printf(msg, vectName);
}
nData = vect->nData;
dt = vect->dx[0];
if (0.0 == dt)
{
printf("dt==0\n");
df = 0.0;
}
else
{
df = 1.0/(dt*nData);
}
/*-------- copy data ------*/
shape[0] = (npy_intp)nData;
if (vect->type == FR_VECT_8C || vect->type == FR_VECT_16C) {
out1 = (PyArrayObject *)PyArray_SimpleNew(1, shape, PyArray_CDOUBLE);
} else {
out1 = (PyArrayObject *)PyArray_SimpleNew(1, shape, PyArray_DOUBLE);
}
data = (double *)PyArray_DATA(out1);
if (data==NULL) {
printf("Unable to allocate space for data.\n");
return;
}
if(vect->type == FR_VECT_2S)
{for(i=0; i<nData; i++) {data[i] = vect->dataS[i];}}
else if(vect->type == FR_VECT_4S)
{for(i=0; i<nData; i++) {data[i] = vect->dataI[i];}}
else if(vect->type == FR_VECT_8S)
{for(i=0; i<nData; i++) {data[i] = vect->dataL[i];}}
else if(vect->type == FR_VECT_1U)
{for(i=0; i<nData; i++) {data[i] = vect->dataU[i];}}
else if(vect->type == FR_VECT_2U)
{for(i=0; i<nData; i++) {data[i] = vect->dataUS[i];}}
else if(vect->type == FR_VECT_4U)
{for(i=0; i<nData; i++) {data[i] = vect->dataUI[i];}}
else if(vect->type == FR_VECT_8U)
{for(i=0; i<nData; i++) {data[i] = vect->dataUL[i];}}
else if(vect->type == FR_VECT_4R)
{for(i=0; i<nData; i++) {data[i] = vect->dataF[i];}}
else if(vect->type == FR_VECT_8R)
{for(i=0; i<nData; i++) {data[i] = vect->dataD[i];}}
// Note the 2*nData in the for loop for complex types
else if(vect->type == FR_VECT_8C)
{for(i=0; i<2*nData; i++) {data[i] = vect->dataF[i];}}
else if(vect->type == FR_VECT_16C)
{for(i=0; i<2*nData; i++) {data[i] = vect->dataD[i];}}
else
{
printf("No numpy type for this channel");
FrVectFree(vect);
FrFileIEnd(iFile);
return;
}
/*------------- fill time and frequency array --------*/
// output2 = x-axis values relative to first data point.
// Usually time, but could be frequency in the case of a frequency
// series
out2 = (PyArrayObject *)PyArray_SimpleNew(1, shape, PyArray_DOUBLE);
data = (double *)PyArray_DATA(out2);
for(i=0; i<nData; i++) {
data[i] = vect->startX[0]+(double)i*dt;
}
// output3 = frequency values in the case of time series
shape[0] = nData/2;
out3 = (PyArrayObject *)PyArray_SimpleNew(1, shape, PyArray_DOUBLE);
data = (double *)PyArray_DATA(out2);
for (i=0;i<nData/2;i++) {
data[i] = (double)i*df;
}
// output4 = gps start time
out4 = (double)vect->GTime;
// output5 = gps start time as a string
utc = vect->GTime - vect->ULeapS + FRGPSTAI;
tempstr = (char *)malloc(200*sizeof(char));
sprintf(tempstr,"Starting GPS time:%.1f UTC=%s",
vect->GTime,FrStrGTime(utc));
out5 = (PyStringObject *)PyString_FromString(tempstr);
free(tempstr);
// output6 = unitX as a string
out6 = (PyStringObject *)PyString_FromString(vect->unitX[0]);
// output7 = unitY as a string
out7 = (PyStringObject *)PyString_FromString(vect->unitY);
/*------------- clean up -----------------------------*/
FrVectFree(vect);
FrFileIEnd(iFile);
return Py_BuildValue("(NNNdNNN)",out1,out2,out3,out4,out5,out6,out7);
}
static PyMethodDef frgetvectMethods[] = {
{"frgetvect", frgetvect, METH_VARARGS, docstring},
{NULL, NULL, 0, NULL}
};
PyMODINIT_FUNC initfrgetvect(void){
(void) Py_InitModule("frgetvect", frgetvectMethods);
}
#!/usr/bin/env python
import frgetvect
# test general call
#output = frgetvect.frgetvect("../data/test.dat","Adc1",750517171,10.,1)
# test call with default values
output = frgetvect.frgetvect("../data/test.dat","Adc1")
print "len(output[0]) = %d" % len(output[0])
print "len(output[1]) = %d" % len(output[1])
print "len(output[2]) = %d" % len(output[2])
print "output[3:] = " + str(output[3:])
print output[0][:5] #should print the first 5 elements of the y-axis.
include $(CMTROOT)/src/Makefile.header
include $(CMTROOT)/src/constituents.make
package Fr
use VirgoPolicy v2r2
macro Fr_linkopts " -L$(FRROOT)/$(Fr_tag) -lFrame -lm "
macro+ Fr_linkopts "" \
Linux " -Wl,-rpath,$(FRROOT)/$(Fr_tag)"
alias FrDump "${FRROOT}/${Fr_tag}/FrDump.exe"
alias FrCopy "${FRROOT}/${Fr_tag}/FrCopy.exe"
alias FrCheck "${FRROOT}/${Fr_tag}/FrCheck.exe"
private
macro+ cflags " -O "
macro utilities_constituents "FrDump FrCopy FrCheck"
macro examples_constituents "FrCompress FrCopyFile FrCopyFrame FrFileDump \
FrFull FrMark FrMultiR FrMultiTOC FrMultiW FrOnline FrReshape FrSpeed FrStat"
macro constituents "Frame $(utilities_constituents) "
macro+ constituents "" \
OSF1 "$(examples_constituents)"
public
# The constituents
# The standard library
library Frame -no_prototypes FrameL.c FrIO.c FrFilter.c -s=zlib *.c
# The standard utilities
application FrDump FrDump.c
application FrCopy FrCopy.c
application FrCheck FrCheck.c
# The examples and test program
application FrCompress exampleCompress.c
application FrCopyFile exampleCopyFile.c
application FrCopyFrame exampleCopyFrame.c
application FrFileDump exampleFileDump.c
application FrFull exampleFull.c
application FrMark exampleMark.c
application FrMultiR exampleMultiR.c
application FrMultiTOC exampleMultiTOC.c
application FrMultiW exampleMultiW.c
application FrOnline exampleOnline.c
application FrReshape exampleReshape.c
application FrSpeed exampleSpeed.c
application FrStat exampleStat.c
File added
%---- define file name and channel name
%
file = '../data/test.dat';
channel = 'fastAdc1';
%
% ---first extarct data from frames -------------
%
[a,t,f,t0,t0s,c,u,more] = frextract(file, channel,2,5);
%
%---------- plot time serie --------------------
%
subplot(2,1,1)
plot(t,a)
ylabel(channel)
xlabel('time [s]')
title(t0s)
%
%------ compute and plot FFT --------------------
%
b = fft(a);
m = abs(b(1:length(b)/2));
subplot(2,1,2)
semilogy(f,m)
ylabel('power')
xlabel('frequency [Hz]')
title(['FFT for ',channel])
%
%---- write an audio file (listen it with Netscape) ----
%
auwrite(a/max(abs(a)),16384.,'audio.au')
%---- define file name and channel name
%
file = '../data/test.dat';
channel = 'fastAdc1';
%
% ---first extract data from frames -------------
%
[a,t,f,t0,t0s,c,u,more] = frextract(file,channel,2,3);
%
%---------- plot time serie --------------------
%
subplot(2,1,1)
plot(t,a)
ylabel(channel)
xlabel('time [s]')
title(t0s)
%
%------ compute and plot FFT --------------------
%
b = fft(a);
m = abs(b(1:length(b)/2));
subplot(2,1,2)
loglog(f,m)
ylabel('power')
xlabel('frequency [Hz]')
title(['FFT for ',channel])
%---- define file name and channel name
%
file = '../data/test.dat';
channel = 'fastProc';
%
% ---first extract data from frames (sequential read) ---
%
[a,t,f,t0,t0s,c,u,more] = frextract(file,channel,2,3);
%
subplot(3,1,1)
plot(t,a)
ylabel(channel)
xlabel('time [s]')
title(t0s)
%
%----- then do a random access to extract the same data --
%
[a,t,f,t0,t0s,c,u] = frgetvect(file,channel,t0,3.);
%
subplot(3,1,2)
plot(t,a)
ylabel(channel)
xlabel('time [s]')
title(t0s)
%
%----- then do a random access to extract part of the data --
%
[a,t,f,t0,t0s,c,u] = frgetvect(file,channel,t0+.6,0.6);
%
subplot(3,1,3)
plot(t,a)
ylabel(channel)
xlabel('time [s]')
title(t0s)
%---- define file name and channel name
%
file = '../data/test.dat';
channel = 'fastAdc1';
%
% ---first get data from frames -------------
%
[a,t,f,t0,t0s,c,u] = frgetvect(file,channel,750517172.,4.);
[a,t,f,t0,t0s,c,u] = frgetvect(file,channel,0.,4.,1);
%
%---------- plot time serie --------------------
%
subplot(2,1,1)
plot(t,a)
ylabel(channel)
xlabel('time [s]')
title(t0s)
%
%------ compute and plot FFT --------------------
%
b = fft(a);
m = abs(b(1:length(b)/2));
subplot(2,1,2)
loglog(f,m)
ylabel('power')
xlabel('frequency [Hz]')
title(['FFT for ',channel])
/*---------------------------------------------------------------------*/
/* FrExtract.c by B. Mours LAPP Oct 22, 2002 */
/* && Increase file name length - K. Thorne - Sept 19, 2006 && */
/* */
/* This Matlab mex file extract from a frame file the data for */
/* one ADC channel (this is for Matlab version 5) */
/* */
/* The input arguments are: */
/* 1) file name(s). This could be a single file, a list of file */
/* separated by space or a frame file list (ffl) */
/* 2) ADC or PROC name */
/* 3) (optional) first frame (default=0 first frame in the file) */
/* 4) (optional) number of frames to read (default = Frame) */
/* 5) (optional) debug level (default = 0 = no output) */
/* */
/* Returned matlab data: */
/* 1) ADC or PROC data (time series) */
/* 2) (optional) x axis values relative to the first data point. */
/* This is usual time but it could be frequency in the case */
/* of a frequency series (double) */
/* 3) (optional) frequency values in the case of time series (double)*/
/* 4) (optional) GPS starting time (in second.nanosec) */
/* 5) (optional) starting time as a string */
/* 6) (optional) ADC comment as a string */
/* 7) (optional) ADC unit as a string */
/* 8) (optional) additional info: it is a 9 words vector which */
/* content:crate #, channel#, nBits#, biais, slope, */
/* sampleRate, timeOffset(S.N), fShift, overRange */
/* All this values are stored as double */
/*---------------------------------------------------------------------*/
#include "mex.h"
#include "FrameL.h"
#define MAX_FILE_NAME 1024
#define MAX_VECT_NAME 100
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{FrFile *iFile;
FrameH *frame;
FrAdcData *adc;
FrVect *vect;
FrProcData *proc;
int status, first, nFrames;
long buffSize, debugLvl, i, nBin, nData, index, utc;
double *data, *datai, dt, df, sampleRate;
char fileName[MAX_FILE_NAME];
char adcName[MAX_VECT_NAME];
char *buff, *name, *comment, *units;
/*--------------- check arguments ---------------------*/
if(nrhs<2)
{mexErrMsgTxt("Please provide the file name(s) and ADC or PROC name as argument\n");
return;}
if(!mxIsChar(prhs[0]))
{mexErrMsgTxt("The first arguement is not a filename\n");
return;}
if(!mxIsChar(prhs[1]))
{mexErrMsgTxt("The second arguement is not the adc name \n");
return;}
if(nlhs==0)
{mexErrMsgTxt("Please provide at least one output argument\n");
return;}
if(nrhs > 2)
{data = mxGetPr(prhs[2]);
first = data[0];}
else{first = 0;}
if(nrhs > 3)
{data = mxGetPr(prhs[3]);
nFrames = data[0];}
else{nFrames = 1;}
if(nrhs > 4)
{data = mxGetPr(prhs[4]);
debugLvl = data[0];
FrLibSetLvl(debugLvl);}
else{debugLvl = 0;}