forked from scottransom/presto
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmaximize_rzw.c
191 lines (176 loc) · 7.52 KB
/
maximize_rzw.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
#include "accel.h"
double max_rzw_arr(fcomplex * data, long numdata, double rin, double zin,
double win, double *rout, double *zout,
double *wout, rderivs * derivs)
/* Return the Fourier frequency, f-dot, and fdotdot that */
/* maximizes the power. */
/* This isn't a true optimization, but just maximizing the */
/* oversampled F/Fdot/Fdotdot volume. */
{
float locpow = get_localpower3d(data, numdata, rin, zin, win);
float maxpow = 0, pow, powargr, powargi;
int extra = 10, startbin = (int)(rin) - 1;
int numr, numz, numw, nextbin, fftlen, numbetween;
// The factor beyond ACCEL_DR, ACCEL_DZ, ACCEL_DW will interpolate
int interpfac = 4, wind = 0, zind = 0, rind = 0;
fcomplex ***vol, amp;
// Use a more conservative length kernel
int kern_half_width = w_resp_halfwidth(fabs(zin)+2, fabs(win)+20, LOWACC);
numr = 3 * interpfac * ACCEL_RDR;
numz = numw = 2 * interpfac + 1;
numbetween = interpfac * ACCEL_RDR;
fftlen = next2_to_n(numbetween * (2 * kern_half_width + extra));
// printf("fftlen = %d\n", fftlen);
vol = corr_rzw_vol(data, numdata, numbetween, startbin,
zin-ACCEL_DZ, zin+ACCEL_DZ, numz,
win-ACCEL_DW, win+ACCEL_DW, numw,
fftlen, LOWACC, &nextbin);
{
int ii, jj, kk;
for (ii = 0; ii < numw; ii++) {
for (jj = 0; jj < numz; jj++) {
for (kk = 0; kk < numr; kk++) {
amp = vol[ii][jj][kk];
pow = POWER(amp.r, amp.i);
if (pow > maxpow) {
maxpow = pow;
wind = ii;
zind = jj;
rind = kk;
}
}
}
}
}
*rout = startbin + rind * ACCEL_DR / (double) interpfac;
*zout = zin-ACCEL_DZ + zind * ACCEL_DZ / (double) interpfac;
*wout = win-ACCEL_DW + wind * ACCEL_DW / (double) interpfac;
vect_free(vol[0][0]);
vect_free(vol[0]);
vect_free(vol);
get_derivs3d(data, numdata, *rout, *zout, *wout, locpow, derivs);
return maxpow;
}
double max_rzw_file(FILE * fftfile, double rin, double zin, double win,
double *rout, double *zout, double *wout, rderivs * derivs)
/* Return the Fourier frequency, f-dot, and fdotdot that */
/* maximizes the power of the candidate in 'fftfile'. */
{
double maxpow, rin_int, rin_frac;
int filedatalen, extra = 10;
long startbin;
fcomplex *filedata;
// Use a more conservative length kernel
int kern_half_width = w_resp_halfwidth(fabs(zin)+2, fabs(win)+20, LOWACC);
rin_frac = modf(rin, &rin_int);
filedatalen = 2 * kern_half_width + extra;
startbin = (long) rin_int - filedatalen / 2;
filedata = read_fcomplex_file(fftfile, startbin, filedatalen);
maxpow = max_rzw_arr(filedata, filedatalen, rin_frac + filedatalen / 2,
zin, win, rout, zout, wout, derivs);
*rout += startbin;
vect_free(filedata);
return maxpow;
}
void max_rzw_arr_harmonics(fcomplex data[], long numdata,
int num_harmonics,
double rin, double zin, double win,
double *rout, double *zout, double *wout,
rderivs derivs[], double powers[])
/* Return the Fourier frequency, f-dot, and f-dotdot that */
/* maximizes the *summed* power of the multi-harmonic candidate */
{
float maxpow = 0, pow, powargr, powargi, locpow;
long lobin, hhlobin = 0;
int hind, ii, jj, kk, extra = 10, nextbin, fftlen;
int interpfac = 4; // Factor beyond ACCEL_D[RZW] we will interpolate
int numwid = 2; // The number of full peak widths we will search
int numbetween = interpfac * ACCEL_RDR;
int numr = 3 * numbetween; // Search 3 full Fourier bins around high harm
int numz = 2 * numwid * interpfac + 1;
int numw = numz;
int rind = 0, zind = 0, wind = 0;
double dr = 1.0 / (double) numbetween;
// The summed power spectrum, initialized to zeros
float ***powsum = gen_f3Darr(numw, numz, numr);
for (ii = 0; ii < numw * numz * numr; ii++)
powsum[0][0][ii] = 0.0;
for (hind = 0; hind < num_harmonics; hind++) {
int n = num_harmonics - hind; // harmonic number, starting from highest
double rh = rin * n, zh = zin * n, wh = win * n;
// Use a more conservative length kernel
int kern_half_width = w_resp_halfwidth(fabs(zh)+2, fabs(wh)+20, LOWACC);
fcomplex ***vol, amp;
double rh_int, rh_frac, hfrac = n / (double) num_harmonics;
locpow = get_localpower3d(data, numdata, rh, zh, wh);
rh_frac = modf(rh, &rh_int);
// Will do 1+ bins below and 1+ bins above rin
lobin = (long) rh_int - 1;
if (hind==0) hhlobin = lobin;
fftlen = next2_to_n(numbetween * (2 * kern_half_width + extra));
// Create the RZW volume for the harmonic.
// Note that we are computing the z and w values in exact harmonic
// ratios. But the r values are on a power-of-two grid.
vol = corr_rzw_vol(data, numdata, numbetween, lobin,
zh-numwid*hfrac*ACCEL_DZ,
zh+numwid*hfrac*ACCEL_DZ, numz,
wh-numwid*hfrac*ACCEL_DW,
wh+numwid*hfrac*ACCEL_DW, numw,
fftlen, LOWACC, &nextbin);
// Calculate and sum the powers
for (ii = 0; ii < numw; ii++) {
for (jj = 0; jj < numz; jj++) {
for (kk = 0; kk < numr; kk++) {
rind = (long) round(((hhlobin + kk * dr) * hfrac
- lobin) * numbetween);
amp = vol[ii][jj][rind];
powsum[ii][jj][kk] += POWER(amp.r, amp.i) / locpow;
}
}
}
vect_free(vol[0][0]);
vect_free(vol[0]);
vect_free(vol);
}
// Now search the power sums for the highest value
rind = zind = wind = 0;
for (ii = 0; ii < numw; ii++) {
for (jj = 0; jj < numz; jj++) {
for (kk = 0; kk < numr; kk++) {
pow = powsum[ii][jj][kk];
if (pow > maxpow) {
maxpow = pow;
wind = ii;
zind = jj;
rind = kk;
}
}
}
}
// Calculate the proper r, z, and w peak values
{
double nh = num_harmonics;
*rout = (hhlobin + rind * dr) / nh;
*zout = ((zin * nh) - numwid * ACCEL_DZ +
zind * ACCEL_DZ / (double) interpfac) / nh;
*wout = ((win * nh) - numwid * ACCEL_DW +
wind * ACCEL_DW / (double) interpfac) / nh;
}
// Now calculate the derivatives at the peak
for (ii = 0; ii < num_harmonics; ii++) {
int hh = ii + 1;
double rr = *rout * hh, zz = *zout * hh, ww = *wout * hh;
locpow = get_localpower3d(data, numdata, rr, zz, ww);
get_derivs3d(data, numdata, rr, zz, ww, locpow, &(derivs[ii]));
powers[ii] = derivs[ii].pow;
}
/*
printf("numr = %d numz = %d numw = %d\n", numr, numz, numw);
printf("rind = %d zind = %d wind = %d\n", rind, zind, wind);
printf("rin = %f zin = %f win = %f\n", rin , zin , win);
printf("rout = %f zout = %f wout = %f\n", *rout, *zout, *wout);
*/
vect_free(powsum[0][0]);
vect_free(powsum[0]);
vect_free(powsum);
}