forked from scottransom/psrfits2psrfits
-
Notifications
You must be signed in to change notification settings - Fork 1
/
testrsc.c
137 lines (116 loc) · 3.01 KB
/
testrsc.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
/******
* testrsc.c SC 20091123
* Driver to test rescaling
* gcc rescale.c testrsc.c -o testrsc -lm
******/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "rescale.h"
#define VERBOSE 1
#define CASE 3
#define NDATA 15625
float randg(void);
float dataval[NDATA], outval[NDATA], outval2[NDATA];
int runtest(void); // Only uses globals
/******************************/
int main(int argc, char *argv[])
{
int i, j, Nval;
#if(CASE==1)
/***** Fake Gaussian data *****/
for (i = 0; i < NDATA; i++) {
dataval[i] = (int) (32000.0 + randg() * 10240.0) % 65536;
/* Poor random values, but who cares? */
}
if (runtest() != 0) {
printf("Error!?\n");
return (-1);
}
#endif
#if(CASE==2)
/***** Fake Gaussian data, bad bias *****/
for (i = 0; i < NDATA; i++) {
dataval[i] = roundf(1000.0 + randg() * 3000.0);
if (dataval[i] < 0.0) {
dataval[i] = rand() % 100;
}
if (dataval[i] > 65536.0) {
dataval[i] = (65536 - rand() % 100);
}
}
if (runtest() != 0) {
printf("Error!?\n");
return (-1);
}
#endif
#if(CASE==3)
/***** Long tail case *****/
for (i = 0; i < NDATA; i++) {
j = rand() % 256;
dataval[i] = (j * j) % 65536;
}
if (runtest() != 0) {
printf("Error!?\n");
return (-1);
}
#endif
#if(CASE==4)
/***** Bad range case *****/
for (i = 0; i < NDATA; i++) {
dataval[i] = ((int) (randg() * 1000.0 + 20000 + rand() % 20000)) % 65536;
}
if (runtest() != 0) {
printf("Error!?\n");
return (-1);
}
#endif
return (0);
}
/******************************/
int runtest(void)
{
int i;
float offset, scale;
// Obtain scale parameters:
if (rescale(dataval, NDATA, &offset, &scale) != 0) {
printf("Rescale routine failed!\n");
return (-1);
}
#if(VERBOSE)
printf("# Calculated offset and scale: %.3f, %.3f\n", offset, scale);
#endif
// Generate rescaled data:
for (i = 0; i < NDATA; i++) {
// outval reflects results of rescale()
outval[i] = roundf((dataval[i] - offset) / scale);
if (outval[i] < 0.0) {
outval[i] = 0.0;
}
// outval2 reflects basicscale(), always 2048.0. 4096.0:
outval2[i] = roundf((dataval[i] - 2048.0) / 4096.0);
}
// Now print it out for plotting:
for (i = 0; i < NDATA; i++) {
printf("%8.1f %8.1f %8.1f\n", dataval[i], outval2[i], outval[i]);
}
return (0);
}
/******************************/
float randg(void)
{
float x1, x2, w;
/* WARNING: Uses system rand()! */
do {
x1 = 2.0 * ((float) rand() / RAND_MAX) - 1.0;
x2 = 2.0 * ((float) rand() / RAND_MAX) - 1.0;
w = x1 * x1 + x2 * x2;
} while (w >= 1.0);
w = sqrt((-2.0 * log(w)) / w);
// x1*w and x2*w are two Gaussian random variables
// (mu=0,sigma=1)
// Return one:
return (x1 * w);
}
/******************************/