-
Notifications
You must be signed in to change notification settings - Fork 0
/
presto2einstein.c
324 lines (255 loc) · 11.8 KB
/
presto2einstein.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
/* presto2einstein.c - convert a presto dedispersed time series and metadata file into an */
/* Einstein@Home .binary work unit - siemion@berkeley.edu - 2012 */
/* */
/* usage: presto2einstein <filestem> (optionally <filestem>.dat or <filestem>.inf) */
/* reads <filestem>.dat and <filestem>.inf - creates <filestem>.binary */
/* */
/* presto input is 32 bit floating point data, output binary is 8 bit signed, mean 0, */
/* quantized at +/- 6 sigma */
/* */
/* Various funcionality has been appropriated from Einstein@Home (Radio Pulsar Edition), */
/* license and attribution below: */
/* */
/***************************************************************************
* Copyright (C) 2008 by Benjamin Knispel, Holger Pletsch *
* benjamin.knispel[AT]aei.mpg.de *
* Copyright (C) 2009, 2010 by Oliver Bock *
* oliver.bock[AT]aei.mpg.de *
* Copyright (C) 2009, 2010 by Heinz-Bernd Eggenstein *
* * * *
* Einstein@Home 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, version 2 of the License. *
* *
* Einstein@Home 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 Einstein@Home. If not, see <http://www.gnu.org/licenses/>. *
* *
***************************************************************************/
#include <stdio.h>
#include <math.h>
#include <inttypes.h>
#include <stdlib.h>
#include <string.h>
// useful definitions
#define FN_LENGTH 256 // length of the filenames
#define N_BINS_SS 40 // number of bins in the screensaver powerspectrum
#define MICROSEC 1.0e-6 // conversion factor from seconds to microseconds
#define MEGAHZ 1.0e6 // conversion factor from Hz to MHz
//header struct from data files
struct data_header
{
double tsample; // Sample time in us
double tobs; // Observation time in s
double timestamp; // Time stamp (MJD)
double fcenter; // Center freq in MHz
double fchan; // Channel band in kHz
double RA; // right ascencsion (J2000)
double DEC; // declination (J2000)
double gal_l; // galactic longitude
double gal_b; // galactic latitude
double AZstart; // AZ at start
double ZAstart; // ZA at start
double ASTstart; // AST at start
double LSTstart; // LST at start
double DM; // trial dispersion measure value pc cm^{-3}
double scale; // scale factor for compressed data. If compressed data value = 1, uncompressed value is 1/scale
uint32_t filesize; // filesize in bytes -> refers to the file generated by filterbank, not the file in which this header sits
uint32_t datasize; // datasize in bytes -> refers to the file generated by filterbank, not the file in which this header sits
uint32_t nsamples; // number of samples
uint16_t smprec; // Number of samples/record
uint16_t nchan; // Number of channels/record
uint16_t nifs; // Nifs = number of polarizations -> for PALFA channel 1 contains sum of two polarizations
uint16_t lagformat; // lagformat
uint16_t sum; // sum
uint16_t level; // level
char name[FN_LENGTH]; // name of the object
char originalfile[FN_LENGTH]; // Original WAPP file
char proj_id[FN_LENGTH]; // project ID
char observers[FN_LENGTH]; // observers
} __attribute__((__packed__));
typedef struct data_header Data_Header;
unsigned char quantize(float d, float min, float max);
int file_exists(char *filename);
void comp_stats(float *mean, float *stddev, float *vec, long int veclen);
int main (int argc, char *argv[]) {
char tmpfilename[255];
char tmpstr1[128];
char tmpstr2[128];
double tmpdbl1=0;
double tmpdbl2=0;
double tmpdbl3=0;
double tmpsign=0;
FILE *fp;
Data_Header einstein_header;
float *timeseries_data;
signed char *quantized_timeseries;
long int i;
/* initialize defaults */
einstein_header.tsample = 1; // Sample time in us
einstein_header.tobs = 1; // Observation time in s
einstein_header.timestamp = 55000; // Time stamp (MJD)
einstein_header.fcenter = 1400; // Center freq in MHz
einstein_header.fchan = 1; // Channel band in kHz
einstein_header.RA = 0.0; // right ascencsion (J2000)
einstein_header.DEC = 0.0; // declination (J2000)
einstein_header.gal_l = 0.0; // galactic longitude
einstein_header.gal_b = 0.0; // galactic latitude
einstein_header.AZstart = 0.0; // AZ at start
einstein_header.ZAstart = 0.0; // ZA at start
einstein_header.ASTstart = 0.0; // AST at start
einstein_header.LSTstart = 0.0; // LST at start
einstein_header.filesize = 1024000;
einstein_header.datasize = 1024000;
einstein_header.nsamples = 1024000;
einstein_header.smprec = 1024;
einstein_header.nchan = 1024;
einstein_header.nifs = 1;
einstein_header.lagformat = 0;
einstein_header.sum = 1;
einstein_header.level = 1;
einstein_header.scale = 100;
sprintf(einstein_header.name, "unknown");
sprintf(einstein_header.originalfile, "unknown");
sprintf(einstein_header.proj_id, "unknown");
sprintf(einstein_header.observers, "unknown");
/* clip .dat or .inf if it was appended to the file stem */
if(strstr(argv[1], ".dat") != NULL) memset(argv[1] + strlen(argv[1]) - 4, 0x0, 4);
if(strstr(argv[1], ".inf") != NULL) memset(argv[1] + strlen(argv[1]) - 4, 0x0, 4);
sprintf(tmpfilename, "%s.inf", argv[1]);
if(file_exists(tmpfilename)) {
fp = fopen(tmpfilename, "rb");
} else {
fprintf(stderr, "%s doesn't exist!\n", tmpfilename);
exit(0);
}
fscanf(fp, " Data file name without suffix = %s", einstein_header.originalfile);
fscanf(fp, " Telescope used = %s", tmpstr1);
fscanf(fp, " Instrument used = %s", tmpstr2);
sprintf(einstein_header.proj_id, "%s_%s", tmpstr1, tmpstr2);
fscanf(fp, " Object being observed = %s", einstein_header.name);
fscanf(fp, " J2000 Right Ascension (hh:mm:ss.ssss) = %s", tmpstr1);
sscanf(tmpstr1, "%lf:%lf:%lf", &tmpdbl1, &tmpdbl2, &tmpdbl3);
einstein_header.RA = tmpdbl1 + tmpdbl2/60.0 + tmpdbl3/3600.0;
fscanf(fp, " J2000 Declination (dd:mm:ss.ssss) = %s", tmpstr1);
sscanf(tmpstr1, "%lf:%lf:%lf", &tmpdbl1, &tmpdbl2, &tmpdbl3);
if(tmpdbl1 < 0) {
tmpsign = -1;
} else {
tmpsign = 1;
}
einstein_header.DEC = tmpsign * (fabs(tmpdbl1) + tmpdbl2/60.0 + tmpdbl3/(60.0 * 360.0));
fscanf(fp, " Data observed by = %s", einstein_header.observers);
fscanf(fp, " Epoch of observation (MJD) = %lf", &(einstein_header.timestamp));
fscanf(fp, " Barycentered? (1=yes, 0=no) = %lf", &tmpdbl1);
fscanf(fp, " Number of bins in the time series = %lf", &tmpdbl1);
fscanf(fp, " Width of each time series bin (sec) = %lf", &tmpdbl2);
einstein_header.tsample = tmpdbl2 * 1000000.0;
einstein_header.tobs = tmpdbl1 * tmpdbl2;
einstein_header.nsamples = (uint32_t) tmpdbl1;
fscanf(fp, " Any breaks in the data? (1=yes, 0=no) = %lf", &tmpdbl1);
fscanf(fp, " Type of observation (EM band) = %s", tmpstr1);
fscanf(fp, " Beam diameter (arcsec) = %lf", &tmpdbl1);
fscanf(fp, " Dispersion measure (cm-3 pc) = %lf", &(einstein_header.DM));
fscanf(fp, " Central freq of low channel (Mhz) = %lf", &(einstein_header.fcenter));
fscanf(fp, " Total bandwidth (Mhz) = %lf", &tmpdbl1);
fscanf(fp, " Number of channels = %" SCNu16, &(einstein_header.nchan));
fscanf(fp, " Channel bandwidth (Mhz) = %lf", &(einstein_header.fchan));
fscanf(fp, " Data analyzed by = %s", tmpstr1);
printf("Header info: \n");
printf("Original File: %s\n", einstein_header.originalfile);
printf("Project ID: %s\n", einstein_header.proj_id);
printf("Object Name: %s\n", einstein_header.name);
printf("RA: %g\n", einstein_header.RA);
printf("Dec: %g\n", einstein_header.DEC);
printf("Observer: %s\n", einstein_header.observers);
printf("Time Stamp: %15.15g\n", einstein_header.timestamp);
printf("Sample Time: %15.15g\n", einstein_header.tsample);
printf("Obs. Length: %15.15g\n", einstein_header.tobs);
printf("NChan: %u\n", einstein_header.nchan);
printf("NSamp: %u\n", einstein_header.nsamples);
//printf("%g %g %g\n", tmpdbl1, tmpdbl2, tmpdbl3);
fclose(fp);
sprintf(tmpfilename, "%s.dat", argv[1]);
if(file_exists(tmpfilename)) {
fp = fopen(tmpfilename, "rb");
} else {
fprintf(stderr, "%s doesn't exist!\n", tmpfilename);
exit(0);
}
timeseries_data = malloc(sizeof(float) * einstein_header.nsamples);
quantized_timeseries = malloc(sizeof(char) * einstein_header.nsamples);
fread(timeseries_data, sizeof(float), einstein_header.nsamples, fp);
fclose(fp);
/* calculate min/max */
float min;
float max;
/* this works ok for well behaved data */
/*
min = timeseries_data[0];
max = timeseries_data[0];
for(i = 1; i < einstein_header.nsamples; i++) {
if(timeseries_data[i] < min) min = timeseries_data[i];
if(timeseries_data[i] > max) max = timeseries_data[i];
}
*/
float mean=0;
float stddev=0;
/* this method is better */
comp_stats(&mean, &stddev, timeseries_data, einstein_header.nsamples);
min = mean - 6 * stddev;
max = mean + 6 * stddev;
printf("%f %f\n", min, max);
/* quantize returns a unsigned char centered at 128, e@h wants a signed char */
for(i = 0; i < einstein_header.nsamples; i++) {
//printf("%u\n", quantize(timeseries_data[i], min, max));
quantized_timeseries[i] = (signed char) ((float) quantize(timeseries_data[i], min, max) - 128.0);
//if(quantized_timeseries[i] < 0) printf("%d\n", quantized_timeseries[i]);
//usleep(50000);
//printf("%d %f\n", quantized_timeseries[i], timeseries_data[i]);
}
free(timeseries_data);
sprintf(tmpfilename, "%s.binary", argv[1]);
fp = fopen(tmpfilename, "wb");
if(fp == NULL) {
fprintf(stderr, "couldn't open %s for writing\n", tmpfilename);
exit(0);
}
fwrite(&einstein_header, sizeof(einstein_header), 1, fp);
fwrite(quantized_timeseries, sizeof(signed char), einstein_header.nsamples, fp);
fclose(fp);
free(quantized_timeseries);
return 0;
}
unsigned char quantize(float d, float min, float max)
{
if(d > max) d = max;
if(d < min) d = min;
return (unsigned char)(((d - min) / (max-min)) * 256.0);
}
int file_exists(char *filename)
{
if ((fopen(filename,"rb"))==NULL) {
return(0);
} else {
return(1);
}
}
void comp_stats(float *mean, float *stddev, float *vec, long int veclen){
//compute mean and stddev of floating point vector vec
long int i,j,k;
double tmean = 0;
double tstddev = 0;
for(i=0;i<veclen;i++) {
tmean = tmean + (double) vec[i];
tstddev = tstddev + ((double) vec[i] * vec[i]);
}
tstddev = pow((tstddev - ((tmean * tmean)/veclen))/(veclen - 1), 0.5);
tmean = tmean / (veclen);
*mean = (float) tmean;
*stddev = (float) tstddev;
}