1 tplarson 1.3 /*
2 this module creates power spectra out of slices of timeseries
3 based on jtsfiddle, but no detrending or gapfilling
4 */
|
5 tplarson 1.1
6 /*
7 * Detrending/gapfilling/differencing module
8 * ts_fiddle_rml upgrades ts_fiddle
9 */
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <strings.h>
14 #include <time.h>
15 #include <sys/time.h>
16 #include <sys/resource.h>
17 #include <math.h>
18 #include <assert.h>
19 #include <fftw3.h>
20
21 #include "jsoc_main.h"
|
22 kehcheng 1.7 #include "fitsio.h"
|
23 tplarson 1.1
24 #define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0]))
25 #define kNOTSPECIFIED "not specified"
26
|
27 tplarson 1.5 char *module_name = "jtsslice";
|
28 tplarson 1.12 char *cvsinfo_jtsslice = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsslice.c,v 1.11 2012/09/24 20:06:23 tplarson Exp $";
|
29 tplarson 1.1
30 /* Command line arguments: */
31 ModuleArgs_t module_args[] =
32 {
33 {ARG_STRING, "in", "", "input data records"},
34 {ARG_STRING, "out", "", "output data series"},
35 {ARG_STRING, "segin" , kNOTSPECIFIED, "name of input segment if not using segment 0"},
36 {ARG_STRING, "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
37 {ARG_STRING, "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
38 {ARG_STRING, "srclink", "SRCDATA", "name of link to source data"},
39 {ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
40 {ARG_INT, "VERB", "1", "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", NULL},
41 {ARG_STRING, "GAPFILE", "", "record or file containing window function"},
42
43 {ARG_INT, "NDT", "-1", "", ""},
|
44 tplarson 1.3 {ARG_STRING, "TTOTAL", NULL, "total length of time in output"},
45
|
46 tplarson 1.1 {ARG_INT, "IFIRST", "0", "", ""},
|
47 tplarson 1.3 // {ARG_INT, "FLIPM", "0", "", ""},
48
|
49 tplarson 1.1 {ARG_INT, "TSOUT", "0", "", ""},
50 {ARG_INT, "FFTOUT", "0", "", ""},
51 {ARG_INT, "FFT1OUT", "0", "", ""},
52 {ARG_INT, "POWOUT", "0", "", ""},
53 {ARG_INT, "MAVGOUT", "0", "", ""},
54 {ARG_INT, "NAVG", "0", "", ""},
55 {ARG_INT, "NSKIP", "0", "", ""},
56 {ARG_INT, "NOFF", "0", "", ""},
57 {ARG_INT, "IMIN", "0", "", ""},
58 {ARG_INT, "IMAX", "-1", "", ""},
59 {ARG_END},
60 };
61
62 #include "saveparm.c"
63 #include "timing.c"
|
64 tplarson 1.3 #include "set_history.c"
|
65 tplarson 1.1
|
66 tplarson 1.2 #define DIE(code) { fprintf(stderr,"jtsslice died with error code %d\n",(code)); return 0;}
|
67 tplarson 1.1
68 /* global variables holding the values of the command line variables. */
|
69 tplarson 1.3 static char *gapfile;
70 static int lmode, n_sampled_out, flipm;
71 static float tsample;
72 static int ifirst, tsout,fftout, fft1out, powout, mavgout;
|
73 tplarson 1.1 static int navg, nskip, noff, imin, imax;
74
75
76 /* Prototypes for local functions. */
77 static double splitting(int l, int m);
78
79 /* Prototypes for external functions */
80 extern void cffti_(int *, float *);
81 extern void cfftb_(int *, float *, float *);
82 static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out,
83 float tsample_in, float *tsample_out,
84 int *num_sections, int *sections);
85 static void extract_data(int n_in, int *gaps, float *data_in, float *data_out);
86
87 /************ Here begins the main module **************/
88 int DoIt(void)
89 {
90 int i;
91 int start_index[2], counts[2], intervals[2];
92 int m, n_in, n_sampled_in;
93 TIME epoch, step, start, stop;
94 tplarson 1.1 int gapsize, istart, istop, data_dim, detrend_order;
95 int npow, mshift;
96 float tsamplex, df1, c;
97 int *agaps;
|
98 tplarson 1.3 int *gaps;
|
99 tplarson 1.1 float *msum, *pow;
100 float *data, *out_data, *in_data;
101 fftwf_plan plan;
102 int num_sections, *sections;
|
103 tplarson 1.2 char tstartstr[100];
|
104 tplarson 1.1
105 int fstatus = 0;
106 fitsfile *fitsptr;
107 long fdims[1];
108 long fpixel[]={1};
109 int *anynul=0;
110
111 int length[2], startind[2], endind[2];
112 float *datarr;
113 int status=DRMS_SUCCESS;
114 int newstat=0;
|
115 tplarson 1.3 char *ttotal;
116 double nseconds;
|
117 tplarson 1.1
118 DRMS_Segment_t *segin = NULL;
119 DRMS_Segment_t *segout = NULL;
120 DRMS_Array_t *inarr = NULL;
121 DRMS_Array_t *outarr = NULL;
122 DRMS_Record_t *inrec = NULL;
123 DRMS_Record_t *outrec = NULL;
124 DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
125 DRMS_RecLifetime_t lifetime;
126 long long histrecnum = -1;
127 DRMS_Segment_t *gapseg = NULL;
128 DRMS_Array_t *gaparr = NULL;
129
130 HIterator_t outKey_list;
131 DRMS_Keyword_t *outKey;
132 TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
133
134 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
135 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
136
|
137 tplarson 1.6 double tstart, tstartsave, tstop, tstopsave, cadence, tstartout, tstopout;
|
138 tplarson 1.1 double wt0, wt1, wt2, wt3, wt;
139 double ut0, ut1, ut2, ut3, ut;
140 double st0, st1, st2, st3, st;
141 double ct0, ct1, ct2, ct3, ct;
142
143 wt0=getwalltime();
144 ct0=getcputime(&ut0, &st0);
145
146 /* Read input parameters. */
|
147 tplarson 1.2 gapfile = (char *)cmdparams_save_str (&cmdparams, "GAPFILE", &newstat);
|
148 tplarson 1.1 n_sampled_out = cmdparams_save_int (&cmdparams, "NDT", &newstat);
149 ifirst = cmdparams_save_int (&cmdparams, "IFIRST", &newstat);
|
150 tplarson 1.3 // flipm = cmdparams_save_int (&cmdparams, "FLIPM", &newstat);
|
151 tplarson 1.1 tsout = cmdparams_save_int (&cmdparams, "TSOUT", &newstat);
152 fftout = cmdparams_save_int (&cmdparams, "FFTOUT", &newstat);
153 fft1out = cmdparams_save_int (&cmdparams, "FFT1OUT", &newstat);
154 powout = cmdparams_save_int (&cmdparams, "POWOUT", &newstat);
155 mavgout = cmdparams_save_int (&cmdparams, "MAVGOUT", &newstat);
156 navg = cmdparams_save_int (&cmdparams, "NAVG", &newstat);
157 nskip = cmdparams_save_int (&cmdparams, "NSKIP", &newstat);
158 noff = cmdparams_save_int (&cmdparams, "NOFF", &newstat);
159 imin = cmdparams_save_int (&cmdparams, "IMIN", &newstat);
160 imax = cmdparams_save_int (&cmdparams, "IMAX", &newstat);
161
|
162 tplarson 1.2 char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
163 char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
164 char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
165 char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
|
166 tplarson 1.1 int seginflag = strcmp(kNOTSPECIFIED, segnamein);
167 int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
|
168 tplarson 1.2 char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
169 char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
|
170 tplarson 1.1 int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
171 int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
172 if (permflag)
173 lifetime = DRMS_PERMANENT;
174 else
175 lifetime = DRMS_TRANSIENT;
176
|
177 tplarson 1.3 ttotal=(char *)cmdparams_save_str(&cmdparams, "TTOTAL", &newstat);
178 status=drms_names_parseduration(&ttotal, &nseconds, 1);
|
179 tplarson 1.9 if (status != DRMS_SUCCESS)
|
180 tplarson 1.2 {
|
181 tplarson 1.3 fprintf(stderr, "ERROR: problem parsing TTOTAL, = %s\n", ttotal);
182 return 1;
|
183 tplarson 1.2 }
184
185 if (newstat)
186 {
187 fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
188 cpsave_decode_error(newstat);
189 return 1;
190 }
191 else if (savestrlen != strlen(savestr))
192 {
193 fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
194 return 1;
195 }
196
|
197 tplarson 1.1
198 /**** Sanity check of input parameters. ****/
199
200 /* Default is to just output the timeseries. */
|
201 tplarson 1.3 if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0))
|
202 tplarson 1.1 tsout=1;
203 /* Only one type of output allowed at a time */
|
204 tplarson 1.3 if ((tsout+fftout+fft1out+powout+mavgout) !=1)
|
205 tplarson 1.1 {
206 fprintf(stderr, "ERROR: only one type of output allowed at a time\n");
|
207 tplarson 1.2 return 1;
|
208 tplarson 1.1 }
209 if (navg <= 0)
210 navg=1;
211 if (nskip <= 0)
212 nskip=1;
213 if ((navg != 1) && (nskip != 1))
214 {
215 fprintf(stderr, "ERROR: one of navg and nskip must equal 1\n");
|
216 tplarson 1.2 return 1;
|
217 tplarson 1.1 }
218 if (noff < 0 || noff >= nskip)
219 {
220 fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");
|
221 tplarson 1.2 return 1;
|
222 tplarson 1.1 }
|
223 tplarson 1.3
224
225 DRMS_Record_t *tempoutrec = drms_create_record(drms_env,
226 outseries,
227 DRMS_TRANSIENT,
228 &status);
229
230 if (status != DRMS_SUCCESS)
|
231 tplarson 1.1 {
|
232 tplarson 1.3 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
|
233 tplarson 1.2 return 1;
|
234 tplarson 1.1 }
235
|
236 tplarson 1.3 DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
237 DRMS_Link_t *srclink = hcon_lookup_lower(&tempoutrec->links, srclinkname);
238
|
239 tplarson 1.12 // char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsslice.c,v 1.11 2012/09/24 20:06:23 tplarson Exp $");
|
240 tplarson 1.3 if (histlink != NULL)
241 {
|
242 tplarson 1.12 histrecnum=set_history(histlink);
|
243 tplarson 1.3 if (histrecnum < 0)
244 {
245 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
246 return 1;
247 }
248 }
249 else
250 {
251 fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
252 }
|
253 tplarson 1.1
254 // these must be present in the output dataseries and variable, not links or constants
|
255 tplarson 1.5 char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"};
|
256 tplarson 1.1
257 int itest;
258 for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
259 {
260 DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
|
261 tplarson 1.9 if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1)
|
262 tplarson 1.1 {
263 fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
264 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
|
265 tplarson 1.2 return 1;
|
266 tplarson 1.1 }
267 }
268
|
269 tplarson 1.6 cadence=drms_getkey_float(tempoutrec, "T_STEP", &status);
270 double chunksecs = n_sampled_out*cadence;
|
271 tplarson 1.3 if (fmod(nseconds,chunksecs) != 0.0)
272 {
273 fprintf(stderr, "ERROR: input parameters seem incompatible (chunksecs must divide nseconds): nseconds = %f, chunksecs = %f\n", nseconds, chunksecs);
274 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
275 return 1;
276 }
277 int ntimechunks = nseconds/chunksecs;
|
278 tplarson 1.6 tsample=cadence;
|
279 tplarson 1.3
280 int mflipout = drms_getkey_int(tempoutrec, "MFLIPPED", &status);
281 if (status != DRMS_SUCCESS) mflipout=0;
282
|
283 tplarson 1.1 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
284
285
286 if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))
287 {
288 DRMS_RecordSet_t *gaprecset = drms_open_records(drms_env, gapfile, &status);
289 if (status != DRMS_SUCCESS || gaprecset == NULL)
290 {
291 fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status);
|
292 tplarson 1.2 return 1;
|
293 tplarson 1.1 }
294 if (gaprecset->n == 0)
295 {
296 fprintf(stderr,"ERROR: gapfile recordset contains no records\n");
|
297 tplarson 1.2 return 1;
|
298 tplarson 1.1 }
299 if (gaprecset->n > 1)
300 {
301 fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");
|
302 tplarson 1.2 return 1;
|
303 tplarson 1.1 }
304 gapseg = drms_segment_lookupnum(gaprecset->records[0], 0);
305 if (gapseg != NULL)
306 gaparr = drms_segment_read(gapseg, DRMS_TYPE_INT, &status);
307 if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL)
308 {
309 fprintf(stderr, "ERROR: problem reading gaps segment: status = %d\n", status);
|
310 tplarson 1.2 return 1;
|
311 tplarson 1.1 }
312 agaps=(int *)(gaparr->data);
313 gapsize=gaparr->axis[0];
314 }
315 else
316 {
317 fits_get_img_size(fitsptr, 1, fdims, &fstatus);
318 gapsize=fdims[0];
319 agaps=(int *)(malloc(gapsize*sizeof(int)));
320 fits_read_pix(fitsptr, TINT, fpixel, gapsize, NULL, agaps, anynul, &fstatus);
321 fits_close_file(fitsptr, &fstatus);
322 if (fstatus)
323 {
324 fprintf(stderr, "ERROR: ");
325 fits_report_error(stderr, fstatus);
|
326 tplarson 1.2 return 1;
|
327 tplarson 1.1 }
328 }
329
|
330 tplarson 1.9 if (verbflag)
331 printf("gapfile read, gapsize = %d\n",gapsize);
|
332 tplarson 1.1
333 data_dim=gapsize;
334 if (n_sampled_out>gapsize)
335 data_dim=n_sampled_out;
336
|
337 tplarson 1.3 num_sections=1;
338 sections = malloc(2*sizeof(int));
339 sections[0] = 0;
340 sections[1] = data_dim-1;
|
341 tplarson 1.1
342 /* allocate temporary storage */
343 gaps=(int *)(malloc(gapsize*sizeof(int)));
344 data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));
345
|
346 tplarson 1.2 DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status);
|
347 tplarson 1.1
|
348 tplarson 1.2 if (status != DRMS_SUCCESS || inrecset == NULL)
|
349 tplarson 1.1 {
350 fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);
|
351 tplarson 1.2 return 1;
|
352 tplarson 1.1 }
353
|
354 tplarson 1.2 if (inrecset->n == 0)
|
355 tplarson 1.1 {
|
356 tplarson 1.4 fprintf(stderr, "ERROR: input recordset contains no records\n");
|
357 tplarson 1.2 return 1;
|
358 tplarson 1.1 }
359
|
360 tplarson 1.2 inrec=inrecset->records[0];
|
361 tplarson 1.1 char *inchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "T_STEP"};
362
363 for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
364 {
365 DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
|
366 tplarson 1.9 if (inkeytest == NULL)
|
367 tplarson 1.1 {
368 fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
|
369 tplarson 1.2 drms_close_records(inrecset, DRMS_FREE_RECORD);
370 return 1;
|
371 tplarson 1.1 }
372 }
373
|
374 tplarson 1.6 tstartsave=drms_getkey_time(inrec, "T_START", NULL); //must be constant across all input records unless GAPFILE is implemented differently
375
|
376 tplarson 1.3 int mflipin = drms_getkey_int(inrec, "MFLIPPED", &status);
377 if (status != DRMS_SUCCESS) mflipin=0;
378
379 if (mflipout != mflipin)
380 flipm=1;
381 else
382 flipm=0;
383
|
384 tplarson 1.1 // changed n_in to gapsize in following call
385 if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections))
386 {
387 fprintf(stderr, "ERROR: problem in extract_gaps\n");
388 return 0;
389 }
390
391 if (n_sampled_out == -1)
392 n_sampled_out = n_sampled_in;
393 df1 = tsamplex*n_sampled_out;
394
|
395 tplarson 1.3 if (fftout || fft1out || powout || mavgout)
|
396 tplarson 1.1 plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data,
397 (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE);
398
399 /* Set sum to zero before starting */
400 if (mavgout)
401 msum = (float *)calloc(n_sampled_out/2,sizeof(float));
402
|
403 tplarson 1.5 if (fft1out || powout || mavgout)
|
404 tplarson 1.1 {
405 /* Set first and last output index if not set as a parameter. */
406 if (imax < imin)
407 {
408 imin=0;
409 imax=n_sampled_out/2-1;
410 }
411 npow=imax-imin+1;
412 pow=(float *)(malloc(n_sampled_out*sizeof(float)));
413 }
414
415
|
416 tplarson 1.5 status=drms_stage_records(inrecset, 1, 0);
417 if (status != DRMS_SUCCESS)
418 {
419 fprintf(stderr, "ERROR: drms_stage_records returned status = %d\n", status);
420 return 1;
421 }
|
422 tplarson 1.1
|
423 tplarson 1.5 int nrecsout = ntimechunks * inrecset->n;
424 DRMS_RecordSet_t *outrecset = drms_create_records(drms_env, nrecsout, outseries, DRMS_PERMANENT, &status);
425 if (status != DRMS_SUCCESS || outrecset == NULL)
|
426 tplarson 1.1 {
|
427 tplarson 1.5 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
428 return 1;
|
429 tplarson 1.1 }
430
|
431 tplarson 1.9 int gapoffset, gapoffset0 = ifirst;
|
432 tplarson 1.8 double tstart0 = tstartsave+ifirst*tsample;
|
433 tplarson 1.6 int irec;
|
434 tplarson 1.5 for (irec=0;irec<inrecset->n;irec++)
|
435 tplarson 1.3 {
436
|
437 tplarson 1.5 inrec = inrecset->records[irec];
|
438 tplarson 1.1
|
439 tplarson 1.5 int lmin=drms_getkey_int(inrec, "LMIN", NULL);
440 int lmax=drms_getkey_int(inrec, "LMAX", NULL);
441 if (lmin != lmax)
|
442 tplarson 1.3 {
|
443 tplarson 1.5 fprintf(stderr, "ERROR: lmin != lmax not yet supported.\n");
|
444 tplarson 1.3 return 0;
445 }
|
446 tplarson 1.5 lmode=lmin;
|
447 tplarson 1.3
|
448 tplarson 1.9 if (verbflag)
449 printf("starting l=%d\n",lmode);
450
|
451 tplarson 1.6 tstart=drms_getkey_time(inrec, "T_START", NULL);
452 tstop =drms_getkey_time(inrec, "T_STOP", NULL);
453 cadence=drms_getkey_float(inrec, "T_STEP", NULL);
454 if (tstart != tstartsave)
455 {
456 // to lift this restriction must be able to specify multiple gap files/records as input
457 sprint_time(tstartstr, tstart, "TAI", 0);
458 fprintf(stderr, "ERROR: all input records must have same T_START, record %d has %s\n", irec, tstartstr);
459 return 0;
460 }
|
461 tplarson 1.1
|
462 tplarson 1.6 n_in=(tstop-tstart)/cadence;
|
463 tplarson 1.5 if (n_in != gapsize)
|
464 tplarson 1.3 {
|
465 tplarson 1.5 fprintf(stderr, "ERROR: gaps seem incompatible with time-series, gapsize=%d, n_in=%d\n", gapsize, n_in);
466 return 0;
|
467 tplarson 1.3 }
|
468 tplarson 1.1
|
469 tplarson 1.5 /* Create output data set. */
470 if (tsout || fftout)
471 {
472 length[0]=2*n_sampled_out;
473 length[1]=lmode+1;
474 }
475 else
476 {
477 if (fft1out)
478 length[0]=2*npow;
479 else
480 length[0]=npow;
481 if (powout || fft1out)
482 length[1]=2*lmode+1;
483 else
484 length[1]=1;
485 }
|
486 tplarson 1.1
|
487 tplarson 1.5 startind[1]=0;
488 endind[1]=lmode;
489
490 if (seginflag)
491 segin = drms_segment_lookup(inrec, segnamein);
|
492 tplarson 1.3 else
|
493 tplarson 1.5 segin = drms_segment_lookupnum(inrec, 0);
494 if (segin == NULL)
|
495 tplarson 1.3 {
|
496 tplarson 1.5 fprintf(stderr, "ERROR: problem reading input segment: recnum = %ld\n", inrec->recnum);
497 return 0;
|
498 tplarson 1.3 }
|
499 tplarson 1.5
|
500 tplarson 1.6 int iset;
|
501 tplarson 1.5 for (iset=0; iset < ntimechunks; iset++)
502 {
503 tstartout=tstart0 + iset * chunksecs;
504 tstopout=tstartout+chunksecs;
505 sprint_time(tstartstr, tstartout, "TAI", 0);
|
506 tplarson 1.9 gapoffset = gapoffset0 + iset*n_sampled_out;
507
508 if (verbflag>1)
509 {
510 printf(" starting tstart = %s, ",tstartstr);
|
511 tplarson 1.11 if (irec == 0)
512 {
513 int ig, gtotal=0;
514 for (ig=0;ig<n_sampled_out;ig++)
515 gtotal+=gaps[gapoffset+ig];
516 float dutycycle = (float)gtotal/n_sampled_out;
517 printf("dutycycle = %f\n", dutycycle);
518 }
|
519 tplarson 1.9 }
|
520 tplarson 1.5
|
521 tplarson 1.9 // set ifirst to gapoffset for reading slice
|
522 tplarson 1.6 ifirst=gapoffset;
|
523 tplarson 1.9 startind[0]=2*(ifirst);
524 endind[0]=2*(ifirst + n_sampled_out) - 1;
|
525 tplarson 1.6 // set ifirst=0 because data is read starting at ifirst by drms_segment_readslice
526 ifirst=0;
|
527 tplarson 1.5
528 inarr = drms_segment_readslice(segin, usetype, startind, endind, &status);
529 if (status != DRMS_SUCCESS || inarr == NULL)
530 {
531 fprintf(stderr, "ERROR: problem reading input segment: recnum = %ld, status = %d\n", inrec->recnum, status);
532 return 0;
533 }
534 datarr=(float *)(inarr->data);
535
536 /*
537 outrec = drms_create_record(drms_env, outseries, lifetime, &status);
538 if (status != DRMS_SUCCESS || outrec==NULL)
539 {
540 fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status);
541 return 0;
542 }
543 */
|
544 tplarson 1.6 // outrec = outrecset->records[ntimechunks*irec+iset];
545 outrec = outrecset->records[iset*inrecset->n + irec];
|
546 tplarson 1.5
|
547 tplarson 1.9 if (histlink != NULL)
|
548 tplarson 1.5 drms_setlink_static(outrec, histlinkname, histrecnum);
|
549 tplarson 1.9 if (srclink != NULL)
|
550 tplarson 1.5 drms_setlink_static(outrec, srclinkname, inrec->recnum);
551
552 if (segoutflag)
553 segout = drms_segment_lookup(outrec, segnameout);
554 else
555 segout = drms_segment_lookupnum(outrec, 0);
556 outarr = drms_array_create(usetype, 2, length, NULL, &status);
557 if (status != DRMS_SUCCESS || outarr == NULL || segout == NULL)
558 {
559 fprintf(stderr,"ERROR: problem creating output array or segment: length = [%d, %d], status = %d", length[0], length[1], status);
560 return 0;
561 }
562 out_data = outarr->data;
|
563 tplarson 1.1
564 /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/
|
565 tplarson 1.5 for (m=0; m<=lmode; m++)
566 {
|
567 tplarson 1.1
|
568 tplarson 1.9 if (verbflag>2)
569 printf(" starting m=%d\n",m);
570
|
571 tplarson 1.5 in_data=datarr + m*2*n_sampled_out;
|
572 tplarson 1.1
|
573 tplarson 1.5 /* Extracts data copy */
574 extract_data(n_sampled_out, gaps+gapoffset, in_data, data);
|
575 tplarson 1.1
|
576 tplarson 1.5 /* pad with zeros if the last point output (n_sampled_out+ifirst)
577 is after the last data points read in (nread). */
578 memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));
579 if ((ifirst+n_sampled_out) >= n_sampled_in)
580 memset(&data[2*n_sampled_in], 0, 2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));
|
581 tplarson 1.1
|
582 tplarson 1.5 if (fftout || fft1out || powout || mavgout)
583 fftwf_execute(plan);
|
584 tplarson 1.1
|
585 tplarson 1.5 if (tsout || fftout)
586 memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float));
|
587 tplarson 1.1
|
588 tplarson 1.5 else if (fft1out)
|
589 tplarson 1.3 {
|
590 tplarson 1.5 /* Do positive m */
591 memcpy(&out_data[2*(lmode+m)*npow], &data[2*imin], 2*npow*sizeof(float));
592 if (m > 0)
|
593 tplarson 1.3 {
|
594 tplarson 1.5 /* Do negative m */
595 for (i=0; i<npow; i++)
596 {
597 /* Conjugate for negative m */
598 out_data[2*((lmode-m)*npow+i)] = data[2*(n_sampled_out-imin-i)];
599 out_data[2*((lmode-m)*npow+i)+1] = -data[2*(n_sampled_out-imin-i)+1];
600 }
|
601 tplarson 1.3 }
602 }
603 else
604 {
|
605 tplarson 1.5 /* Compute power spectrum from complex FFT. */
606 for (i = 0; i < n_sampled_out; i++)
607 pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1];
608 if (powout)
609 {
610 /* Do positive m */
611 memcpy(&out_data[(lmode+m)*npow], &pow[imin], npow*sizeof(float));
612 if (m > 0)
613 {
614 /* Do negative m */
615 if (imin)
616 out_data[(lmode-m)*npow] = pow[n_sampled_out-imin];
617 else
618 out_data[(lmode-m)*npow] = pow[0];
619 for (i=1; i<npow;i++)
620 out_data[(lmode-m)*npow+i] = pow[n_sampled_out-imin-i];
621 }
622 }
623 else
|
624 tplarson 1.3 {
|
625 tplarson 1.5 /* m-averaged output */
626 /* Sum all freqs, not only those in output. Should be fixed. */
627 /* Maybe should allow for wrapping around Nyquist. */
628 /* Do positive m */
629 mshift = floor(df1*splitting(lmode,m)+0.5f);
630 if (mshift >= 0)
631 for (i=0; i<n_sampled_out/2-mshift; i++)
632 msum[i] += pow[mshift+i];
|
633 tplarson 1.3 else
|
634 tplarson 1.5 for (i=0; i<n_sampled_out/2+mshift; i++)
635 msum[i-mshift] += pow[i];
636 if (m > 0)
637 {
638 /* Do negative m */
639 mshift = floor(df1*splitting(lmode,-m)+0.5f);
640 if (mshift >=0)
641 for (i=1; i<n_sampled_out/2-mshift;i++)
642 msum[i] += pow[n_sampled_out-mshift-i];
643 else
644 for (i=1; i<n_sampled_out/2+mshift; i++)
645 msum[i-mshift] += pow[n_sampled_out-i];
646 }
|
647 tplarson 1.3 }
648 }
649
|
650 tplarson 1.5 } /* End of m loop */
651
652 drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
653 drms_setkey_int(outrec, "QUALITY", 0);
654 drms_setkey_time(outrec, "T_START", tstartout);
655 drms_setkey_time(outrec, "T_STOP", tstopout);
656 drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2);
657 drms_setkey_float(outrec, "T_STEP", tsamplex);
658 drms_setkey_int(outrec, "NDT", n_sampled_out);
659
660 tnow = (double)time(NULL);
661 tnow += UNIX_epoch;
662 drms_setkey_time(outrec, "DATE", tnow);
|
663 tplarson 1.3
|
664 tplarson 1.5 if (mavgout)
665 {
666 c=1.0f/(2*lmode+1);
667 for (i=0; i<npow; i++)
668 out_data[i] = msum[imin+i] * c;
669 }
|
670 tplarson 1.3
|
671 tplarson 1.9 outarr->bzero=segout->bzero;
672 outarr->bscale=segout->bscale;
|
673 tplarson 1.5 status=drms_segment_write(segout, outarr, 0);
674 if (status != DRMS_SUCCESS)
675 {
676 fprintf(stderr, "ERROR: problem writing output segment: status = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", status, tstartstr, lmode, histrecnum);
677 return 0;
678 }
|
679 tplarson 1.8 drms_free_array(outarr);
680 drms_free_array(inarr);
|
681 tplarson 1.5 // drms_close_record(outrec, DRMS_INSERT_RECORD);
|
682 tplarson 1.1
|
683 tplarson 1.3 }
|
684 tplarson 1.1
|
685 tplarson 1.3 }
|
686 tplarson 1.1
|
687 tplarson 1.5 drms_close_records(inrecset, DRMS_FREE_RECORD);
688 drms_close_records(outrecset, DRMS_INSERT_RECORD);
|
689 tplarson 1.2
|
690 tplarson 1.3 if (fftout || fft1out || powout || mavgout)
|
691 tplarson 1.1 fftwf_destroy_plan(plan);
692
693 wt=getwalltime();
694 ct=getcputime(&ut, &st);
695 if (verbflag)
696 {
697 fprintf(stdout, "total time spent: %.2f ms wall time, %.2f ms cpu time\n",
698 wt-wt0, ct-ct0);
699 }
700
|
701 tplarson 1.3 printf("module %s successful completion\n", cmdparams.argv[0]);
|
702 tplarson 1.1 return 0;
703 }
704
705
706 static double splitting(int l, int m)
707 {
708 double a1x[] = {406.0,407.0,408.0,410.5,412.0,411.5,409.0,406.0,406.0};
709 double a3x[] = {10.0,10.0,15.0,19.5,21.0,21.3,21.3,21.3,21.8};
710 double a5x[] = {0.0,0.0,0.0,-1.5,-2.5,-3.5,-4.0,-4.0,-4.5};
711 /* f0 is the frequency used for generating the above a-coeff */
712 double f0 = 1500.0;
713 /* fcol is the frequency for which to optimize the collaps */
714 double fcol = 800.0;
715
716 double ll,x,x3,x5,a1,a2,a3,a5,frac,lx;
717 int ix;
718
719 if (l == 0)
720 ll = 1;
721 else
722 ll = sqrt(l*(l+1.));
723 tplarson 1.1 x = m/ll;
724 x3 = x*x*x;
725 x5 = x3*x*x;
726 lx = 5*log10(l*f0/fcol)-4;
727 ix = floor(lx);
728 frac = lx-ix;
729 if (lx <= 0) {
730 ix = 0;
731 frac = 0.0;
732 }
733 if (lx >=8) {
734 ix = 7;
735 frac = 1.0;
736 }
737 a1 = (1.0-frac)*a1x[ix]+frac*a1x[ix+1];
738 a2 = 0.0;
739 a3 = (1.0-frac)*a3x[ix]+frac*a3x[ix+1];
740 a5 = (1.0-frac)*a5x[ix]+frac*a5x[ix+1];
741
742 return ll*(a1*x+a2*(1.5*x*x-0.5)+a3*(2.5*x3-1.5*x)+a5*(63./8.*x5-70./8.*x3+15./8.*x))/1e9;
743 }
744 tplarson 1.1
745
746
747 /* Extract the data samples actually used by skipping or averaging
748 data. Replace missing data that are not marked as gaps with zero. */
749 void extract_data(int n_in, int *gaps, float *data_in, float *data_out)
750 {
751 int i,j, nmiss = 0;
|
752 tplarson 1.2 // the following check will be skipped in production executables since they will be built with NDEBUG defined.
753 // it doesn't matter because the check is already done in DoIt().
|
754 tplarson 1.1 assert(nskip==1 || navg==1);
755 if ((navg == 1) && (nskip == 1))
756 {
757 for (i=0; i<n_in; i++)
758 {
759 if (gaps[i]==0 )
760 {
761 data_out[2*i] = 0.0f;
762 data_out[2*i+1] = 0.0f;
763 }
764 else
765 {
766 if (isnan(data_in[2*i]) || isnan(data_in[2*i+1]))
767 {
768 data_out[2*i] = 0.0f;
769 data_out[2*i+1] = 0.0f;
770 gaps[i] = 0;
771 nmiss++;
772 }
773 else
774 {
775 tplarson 1.1 data_out[2*i] = data_in[2*i];
776 data_out[2*i+1] = flipm ? -data_in[2*i+1] : data_in[2*i+1];
777 }
778 }
779 }
780 }
781 else if (nskip != 1)
782 {
783 for (i = 0; i<n_in; i++)
784 {
785 if (gaps[i]==0 )
786 {
787 data_out[2*i] = 0.0f;
788 data_out[2*i+1] = 0.0f;
789 }
790 else
791 {
792 int ix = nskip*i+noff;
793 if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
794 {
795 data_out[2*i] = 0.0f;
796 tplarson 1.1 data_out[2*i+1] = 0.0f;
797 gaps[i] = 0;
798 nmiss++;
799 }
800 else
801 {
802 data_out[2*i] = data_in[2*ix];
803 data_out[2*i+1] = flipm ? -data_in[2*ix+1] : data_in[2*ix+1];
804 }
805 }
806 }
807 }
808 else if (navg != 1)
809 {
810 for (i = 0; i<n_in; i++)
811 {
812 if (gaps[i]==0 )
813 {
814 data_out[2*i] = 0.0f;
815 data_out[2*i+1] = 0.0f;
816 }
817 tplarson 1.1 else
818 {
819 float avgr = 0.0f;
820 float avgi = 0.0f;
821 for (j=0; j<navg; j++)
822 {
823 int ix = navg*i+j;
824 if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
825 {
826 data_out[2*i] = 0.0f;
827 data_out[2*i+1] = 0.0f;
828 gaps[i] = 0;
829 nmiss++;
830 break;
831 }
832 else
833 {
834 avgr += data_in[2*ix];
835 avgi += data_in[2*ix+1];
836 }
837 }
838 tplarson 1.1 if (j == navg)
839 {
840 data_out[2*i] = avgr/navg;
841 data_out[2*i+1] = avgi/navg;
842 }
843 }
844 }
845 }
846 }
847
848 /* Extract the windows function for the samples actually used after
849 subsampling or averaging. Modify the list of continous sections
850 accordingly, and make sure we do not average across section boundaries. */
851 static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out,
852 float tsample_in, float *tsample_out,
853 int *num_sections, int *sections)
854 {
|
855 tplarson 1.2 // the following check will be skipped in production executables since they will be built with NDEBUG defined.
856 // it doesn't matter because the check is already done in DoIt().
|
857 tplarson 1.1 assert(nskip==1 || navg==1);
858 int i,j,sect, start,stop;
859
860
861 if (*num_sections<1)
862 {
863 fprintf(stderr,"Apparantly no sections of data available.");
864 return 1;
865 }
866 /* Mask out data that in between sections if this hasn't already been
867 done in gapfile. */
868 for (i=0; i<sections[0] && i<n_in; i++)
869 gaps_in[i] = 0;
870 for (sect=1; sect<(*num_sections); sect++)
871 {
872 for (i=sections[2*sect-1]+1; i<sections[2*sect] && i<n_in; i++)
873 gaps_in[i] = 0;
874 }
875 for (i=sections[2*(*num_sections-1)+1]+1; i<n_in; i++)
876 gaps_in[i] = 0;
877
878 tplarson 1.1
879 if ((navg == 1) && (nskip == 1))
880 {
881 *n_out = n_in;
882 *tsample_out = tsample_in;
883 for (i=0; i<*n_out; i++)
884 gaps_out[i] = gaps_in[i];
885 }
886 else if (nskip != 1)
887 {
888 *n_out = n_in/nskip;
889 *tsample_out = nskip*tsample_in;
890 for (i=0; i<*n_out; i++)
891 {
892 gaps_out[i] = gaps_in[nskip*i+noff];
893 }
894 for (i=0; i<2*(*num_sections); i+=2)
895 {
896 start = (int)ceilf(((float)(sections[i]-noff))/nskip);
897 stop = (int)floorf(((float)(sections[i+1]-noff))/nskip);
898 if (start <= stop)
899 tplarson 1.1 {
900 sections[i] = start;
901 sections[i+1] = stop;
902 }
903 else
904 {
905 /* This section was skipped entirely. */
906 for (j=i; j< 2*(*num_sections-1); j+=2)
907 {
908 sections[j] = sections[j+2];
909 sections[j+1] = sections[j+3];
910 }
911 *num_sections -= 1;
912 }
913 }
914 }
915 else if (navg != 1)
916 {
917 sect = 0;
918 *n_out = n_in/navg;
919 *tsample_out = tsample*navg;
920 tplarson 1.1 for (i=0; i<*n_out; i++)
921 {
922 int igx = 1;
923 while (sect < *num_sections &&
924 !(navg*i >= sections[2*sect] && navg*i >= sections[2*sect+1]))
925 sect++;
926
927 /* Don't allow averaging of data from different sections. */
928 if (navg*i >= sections[2*sect] && (navg*i+navg-1)<=sections[2*sect+1])
929 {
930 for (j=0; j<navg; j++)
931 igx *= gaps_in[navg*i+j];
932 gaps_out[i] = igx;
933 }
934 else
935 gaps_out[i] = 0;
936 }
937 for (i=0; i<2*(*num_sections); i+=2)
938 {
939 start = (int)ceilf(((float)sections[i])/navg);
940 stop = (int)floorf(((float)sections[i+1])/navg);
941 tplarson 1.1 if (start <= stop)
942 {
943 sections[i] = start;
944 sections[i+1] = stop;
945 }
946 else
947 {
948 /* This section was skipped entirely. */
949 for (j=i; j< 2*(*num_sections-1); j+=2)
950 {
951 sections[j] = sections[j+2];
952 sections[j+1] = sections[j+3];
953 }
954 *num_sections -= 1;
955 }
956 }
957 }
958 return 0;
959 }
960
|
Return to
jtsslice.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps