Loading [MathJax]/extensions/TeX/AMSsymbols.js
LAL 7.7.0.1-3a66518
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
TimeSeriesInterpTest.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014,2015 Kipp Cannon
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19
20
21#include <math.h>
22#include <stdio.h>
23
24
25#include <lal/Date.h>
26#include <lal/LALDatatypes.h>
27#include <lal/TimeSeries.h>
28#include <lal/TimeSeriesInterp.h>
29#include <lal/Units.h>
30
31
33
34
35static REAL8TimeSeries *new_series(double deltaT, unsigned length, double fill)
36{
37 REAL8TimeSeries *new = XLALCreateREAL8TimeSeries("blah", &gps_zero, 0.0, deltaT, &lalDimensionlessUnit, length);
38 unsigned i;
39 for(i = 0 ; i < new->data->length; i++)
40 new->data->data[i] = fill;
41 return new;
42}
43
44
46{
47 return XLALCutREAL8TimeSeries(src, 0, src->data->length);
48}
49
50
51static LIGOTimeGPS t_i(const REAL8TimeSeries *s, unsigned i)
52{
53 LIGOTimeGPS t = s->epoch;
54 XLALGPSAdd(&t, i * s->deltaT);
55 return t;
56}
57
58
59static void add_sine(REAL8TimeSeries *s, LIGOTimeGPS epoch, double ampl, double freq)
60{
61 double t0 = XLALGPSDiff(&s->epoch, &epoch);
62 unsigned i;
63
64 for(i = 0; i < s->data->length; i++) {
65 double t = t0 + i * s->deltaT;
66 s->data->data[i] = ampl * sin(LAL_TWOPI * freq * t);
67 }
68}
69
70
71static void evaluate(REAL8TimeSeries *dst, LALREAL8TimeSeriesInterp *interp, int bounds_check)
72{
73 unsigned i;
74
75 for(i = 0; i < dst->data->length; i++) {
76 LIGOTimeGPS t = t_i(dst, i);
77 dst->data->data[i] = XLALREAL8TimeSeriesInterpEval(interp, &t, bounds_check);
78 }
79}
80
81
83{
84 REAL8TimeSeries *result = copy_series(s1);
85 unsigned i;
86
87 /* assumes series are compatible */
88
89 for(i = 0; i < s1->data->length; i++)
90 result->data->data[i] -= s0->data->data[i];
91
92 return result;
93}
94
95
96static double RMS(const REAL8TimeSeries *s)
97{
98 double rms = 0.;
99 unsigned i;
100
101 for(i = 0; i < s->data->length; i++)
102 rms += s->data->data[i] * s->data->data[i];
103
104 return sqrt(rms / s->data->length);
105}
106
107
108static void minmax(const REAL8TimeSeries *s, double *min, double *max)
109{
110 unsigned i;
111
112 *min = *max = s->data->data[0];
113 for(i = 1; i < s->data->length; i++) {
114 if(s->data->data[i] < *min)
115 *min = s->data->data[i];
116 if(s->data->data[i] > *max)
117 *max = s->data->data[i];
118 }
119}
120
121
122static void check_result(const REAL8TimeSeries *model, const REAL8TimeSeries *result, double rms_bound, double residual_min, double residual_max)
123{
124 REAL8TimeSeries *err = error(model, result);
125 double rms, min, max;
126 rms = RMS(err);
127 minmax(err, &min, &max);
128
129 fprintf(stderr, "error vector: RMS=%g, min=%g, max=%g\n", rms, min, max);
131 if(rms > rms_bound || min < residual_min || max > residual_max) {
132 fprintf(stderr, "error vector larger than allowed\n");
133 exit(1);
134 }
135}
136
137
138int main(void)
139{
140 REAL8TimeSeries *src, *dst, *mdl;
141 LALREAL8TimeSeriesInterp *interp;
142 double f;
143
144 /*
145 * single frequency. include many cycles to approximate source
146 * data with infinite extent. interpolate 1 cycle at high
147 * resolution and compare to model.
148 */
149
150 f = 4000.; /* close to 50% of Nyquist */
151
152 src = new_series(1.0 / 16384, 1024 * 1024, 0.0);
153 add_sine(src, src->epoch, 1.0, f);
154
155 mdl = new_series(1. / 1e8, round(1. / f * 1e8), 0.0);
156 XLALGPSAdd(&mdl->epoch, src->data->length * src->deltaT * .4);
157 dst = copy_series(mdl);
158
159 fprintf(stderr, "interpolating unit amplitude %g kHz sine function sampled at %g Hz to %g MHz\n", f / 1000., 1.0 / src->deltaT, 1.0 / dst->deltaT / 1e6);
160
161 add_sine(mdl, src->epoch, 1.0, f);
162
163 /* in the interpolator's current incarnation, a kernel length this
164 * size at a sample rate of 16384 Hz leads to an internal "no-op
165 * threshold" of about 0.2 ns --- if two samples have different GPS
166 * times then they get their own interpolating kernels, i.e. we
167 * defeat the kernel caching mechanism so that we can watch the
168 * behaviour of the kernel sample-by-sample */
169 interp = XLALREAL8TimeSeriesInterpCreate(src, 65535, NULL, NULL);
170 evaluate(dst, interp, 1);
172
173#if 0
174 {
175 unsigned i;
176 FILE *output = fopen("input.txt", "w");
177 for(i = 0; i < src->data->length; i++) {
178 LIGOTimeGPS t = t_i(src, i);
179 fprintf(output, "%u.%09u %.16g\n", t.gpsSeconds, t.gpsNanoSeconds, src->data->data[i]);
180 }
181 fclose(output);
182 output = fopen("output.txt", "w");
183 for(i = 0; i < mdl->data->length; i++) {
184 LIGOTimeGPS t = t_i(mdl, i);
185 fprintf(output, "%u.%09u %.16g %.16g\n", t.gpsSeconds, t.gpsNanoSeconds, mdl->data->data[i], dst->data->data[i]);
186 }
187 fclose(output);
188 }
189#endif
190
191 check_result(mdl, dst, 1.3e-10, -3.6e-10, +3.6e-10);
192
196
197 /*
198 * higher single frequency, shorter filter, to check for phase
199 * error. interpolate three cycles this time.
200 */
201
202 f = 6500.; /* about 80% of Nyquist */
203
204 src = new_series(1.0 / 16384, 256, 0.0);
205 add_sine(src, src->epoch, 1.0, f);
206
207 mdl = new_series(1. / 1e8, round(3. / f * 1e8), 0.0);
208 XLALGPSAdd(&mdl->epoch, src->data->length * src->deltaT * .4);
209 dst = copy_series(mdl);
210
211 fprintf(stderr, "interpolating unit amplitude %g kHz sine function sampled at %g Hz to %g MHz\n", f / 1000., 1.0 / src->deltaT, 1.0 / dst->deltaT / 1e6);
212
213 add_sine(mdl, src->epoch, 1.0, f);
214
215 interp = XLALREAL8TimeSeriesInterpCreate(src, 9, NULL, NULL);
216 evaluate(dst, interp, 1);
218
219#if 0
220 {
221 unsigned i;
222 FILE *output = fopen("input.txt", "w");
223 for(i = 0; i < src->data->length; i++) {
224 LIGOTimeGPS t = t_i(src, i);
225 fprintf(output, "%u.%09u %.16g\n", t.gpsSeconds, t.gpsNanoSeconds, src->data->data[i]);
226 }
227 fclose(output);
228 output = fopen("output.txt", "w");
229 for(i = 0; i < mdl->data->length; i++) {
230 LIGOTimeGPS t = t_i(mdl, i);
231 fprintf(output, "%u.%09u %.16g %.16g\n", t.gpsSeconds, t.gpsNanoSeconds, mdl->data->data[i], dst->data->data[i]);
232 }
233 fclose(output);
234 }
235#endif
236
237 check_result(mdl, dst, 0.03, -0.078, +0.083);
238
242
243 /*
244 * test behaviour in last sample. allocate series 1 sample longer
245 * than we need it to be so we can control the value of the data
246 * beyond the end of the array, and make sure the interpolator
247 * doesn't return that value.
248 */
249
250 src = new_series(1.0 / 16384, 257, 2.0);/* all 2s */
251 src->data->length--; /* fake the length */
252 src->data->data[src->data->length] = 1; /* place a 1 beyond the end */
253
254 interp = XLALREAL8TimeSeriesInterpCreate(src, 9, NULL, NULL);
255 {
256 LIGOTimeGPS t = src->epoch;
257 double result;
258 XLALGPSAdd(&t, src->data->length * src->deltaT);
259 /* evalute at time of sample beyond end. this should fail */
260 fprintf(stderr, "checking for out-of-bounds failure ...\n");
261 result = XLALREAL8TimeSeriesInterpEval(interp, &t, 1);
262 if(!XLAL_IS_REAL8_FAIL_NAN(result)) {
263 fprintf(stderr, "error: interpolator failed to report error beyond end of array\n");
264 exit(1);
265 } else
266 fprintf(stderr, "... passed\n");
267 /* these should work and return 0., not 1. and not 2 */
268 result = XLALREAL8TimeSeriesInterpEval(interp, &t, 0);
269 if(result != 0.) {
270 fprintf(stderr, "error: interpolator failed in final sample (expected %.16g got %.16g)\n", 0., result);
271 exit(1);
272 }
273 XLALGPSAdd(&t, -1e-9);
274 result = XLALREAL8TimeSeriesInterpEval(interp, &t, 1);
275 if(result != 0.) {
276 fprintf(stderr, "error: interpolator failed in final sample (expected %.16g got %.16g)\n", 0., result);
277 exit(1);
278 }
279 }
281
283
284 /*
285 * success
286 */
287
288 exit(0);
289}
#define fprintf
LALREAL8TimeSeriesInterp * XLALREAL8TimeSeriesInterpCreate(const REAL8TimeSeries *series, int kernel_length, void(*kernel)(double *, int, double, void *), void *kernel_data)
Create a new REAL8TimeSeries interpolator associated with the given REAL8TimeSeries object.
void XLALREAL8TimeSeriesInterpDestroy(LALREAL8TimeSeriesInterp *interp)
Free a LALREAL8TimeSeriesInterp object.
REAL8 XLALREAL8TimeSeriesInterpEval(LALREAL8TimeSeriesInterp *interp, const LIGOTimeGPS *t, int bounds_check)
Evaluate a LALREAL8TimeSeriesInterp at the LIGOTimeGPS t.
static void minmax(const REAL8TimeSeries *s, double *min, double *max)
static REAL8TimeSeries * error(const REAL8TimeSeries *s1, const REAL8TimeSeries *s0)
static REAL8TimeSeries * new_series(double deltaT, unsigned length, double fill)
static void check_result(const REAL8TimeSeries *model, const REAL8TimeSeries *result, double rms_bound, double residual_min, double residual_max)
int main(void)
static LIGOTimeGPS gps_zero
static LIGOTimeGPS t_i(const REAL8TimeSeries *s, unsigned i)
static void evaluate(REAL8TimeSeries *dst, LALREAL8TimeSeriesInterp *interp, int bounds_check)
static void add_sine(REAL8TimeSeries *s, LIGOTimeGPS epoch, double ampl, double freq)
static REAL8TimeSeries * copy_series(const REAL8TimeSeries *src)
static double RMS(const REAL8TimeSeries *s)
static double f(double theta, double y, double xi)
Definition: XLALMarcumQ.c:258
#define LAL_TWOPI
2*pi is circumference of a circle divided by its radius
Definition: LALConstants.h:180
#define LIGOTIMEGPSZERO
Zero-initializer for LIGOTimeGPS structs.
Definition: LALDatatypes.h:464
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
REAL8TimeSeries * XLALCutREAL8TimeSeries(const REAL8TimeSeries *series, size_t first, size_t length)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalDimensionlessUnit
dimensionless units
Definition: UnitDefs.c:156
#define XLAL_IS_REAL8_FAIL_NAN(val)
Tests if val is a XLAL REAL8 failure NaN.
Definition: XLALError.h:393
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
Adds a double to a GPS time.
Definition: XLALTime.c:131
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
Difference between two GPS times as double.
Definition: XLALTime.c:157
Epoch relative to GPS epoch, see LIGOTimeGPS type for more details.
Definition: LALDatatypes.h:458
INT4 gpsSeconds
Seconds since 0h UTC 6 Jan 1980.
Definition: LALDatatypes.h:459
INT4 gpsNanoSeconds
Residual nanoseconds.
Definition: LALDatatypes.h:460
Time series of REAL8 data, see DATATYPE-TimeSeries types for more details.
Definition: LALDatatypes.h:580
REAL8Sequence * data
The sequence of sampled data.
Definition: LALDatatypes.h:586
REAL8 deltaT
The time step between samples of the time series in seconds.
Definition: LALDatatypes.h:583
LIGOTimeGPS epoch
The start time of the time series.
Definition: LALDatatypes.h:582
REAL8 * data
Pointer to the data array.
Definition: LALDatatypes.h:159
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:158
enum @1 epoch
void output(int gps_sec, int output_type)
Definition: tconvert.c:440