24 #include "Tdrt-TDT-Reference.h" 36 TdrtTDTReference::TdrtTDTReference()
41 TimeSeries<float>& h_output,
44 assert(h_input.get_length() == h_output.get_length());
45 assert(h_input.get_delta_t() == h_output.get_delta_t());
47 long long numTimeSamps= h_input.get_length();;
48 double sampTime= h_input.get_delta_t();
49 double pulsarTime= 0.0;
51 double earthTime= 0.0;
53 const double c = 299792458.0;
55 float obsLen = sampTime * numTimeSamps;
56 float tau0 = sampTime/(1 + acceleration * obsLen/(2 * c));
62 float* input = h_input.get_data();
63 float* output = h_output.get_data();
64 for (nsamp=0; nsamp < numTimeSamps; nsamp++)
67 earthTime = nsamp * sampTime;
69 timeThru = pulsarTime;
72 pulsarTime = pulsarTime + tau0 * (1 + acceleration * earthTime/c);
75 samp_out = (int) round(timeThru/sampTime);
80 if (samp_out >= 0 && samp_out < numTimeSamps)
83 output[samp_out] = input[nsamp];
94 void TdrtTDTReference::const_accel_pulsar(
double acceleration,
float period,
95 TimeSeries<float>& h_output)
97 const double c = 299792458.0;
98 const double tsamp = h_output.get_delta_t();
100 float f0= 1.0 / period;
101 float fdot= acceleration / c * f0;
102 float Nrot, Nrot_next;
104 float* tsout = h_output.get_data();
109 for(
size_t i= 0; i < h_output.get_length(); i++)
111 Nrot = tsamp * i * f0 + 0.5 * fdot * tsamp *tsamp * i * i;
112 Nrot_next = tsamp * (i+1) * f0 + 0.5 * fdot * tsamp * tsamp * (i+1) * (i+1);
116 if( ((
int) floor(Nrot_next) - (int) floor(Nrot)) == 1)
131 void TdrtTDTReference::sine_calibration_signal(
float period, TimeSeries<float>& h_output)
133 const float tsamp = h_output.get_delta_t();
134 const float pi = 3.14159265;
135 float* tsout = h_output.get_data();
137 for(
size_t i= 0; i < h_output.get_length(); i++)
139 tsout[i]= sin( 2.0*pi*i * tsamp / period );
143 double TdrtTDTReference::resampled_period(
double acceleration,
double period,
double tsamp,
int length)
145 double c = 299792458.0;
148 new_period = period/(1+(acceleration*tsamp*length)/(2*c));
150 std::cout <<
"Period after resampling: " << new_period << std::endl;
155 float TdrtTDTReference::tdrt_correlation(TimeSeries<float>& h_tsA, TimeSeries<float>& h_tsB)
157 assert(h_tsA.get_length() == h_tsB.get_length());
158 assert(h_tsA.get_delta_t() == h_tsB.get_delta_t());
160 size_t length = h_tsA.get_length();
161 float* tsA = h_tsA.get_data();
162 float* tsB = h_tsB.get_data();
165 double correlation = 0.0;
167 for(
size_t i=0; i<length; ++i)
169 correlation += tsA[i] * tsB[i];
172 if((tsA[i] != 0) && (found==0))
177 if((tsA[i] == 0) && (found==1))
187 return correlation/npul;
190 float TdrtTDTReference::tdrt_bin_diff(TimeSeries<float>& h_tsA, TimeSeries<float>& h_tsB)
192 assert(h_tsA.get_length() == h_tsB.get_length());
193 assert(h_tsA.get_delta_t() == h_tsB.get_delta_t());
195 size_t length = h_tsA.get_length();
196 float* tsA = h_tsA.get_data();
197 float* tsB = h_tsB.get_data();
210 float* Alist = (
float*)malloc(
sizeof(
float)*length/2);
211 float* Blist = (
float*)malloc(
sizeof(
float)*length/2);
213 for(i=0; i<length; ++i)
221 else if((tsA[i] == 0) && (counting == 1))
223 Alist[npulA]=TOA_A/amp;
231 for(i=0; i<length; ++i)
239 else if((tsB[i] == 0) && (counting == 1))
241 Blist[npulB]=TOA_B/amp;
249 if(npulA<npulB) npul=npulA;
252 for(i=0; i<npul; ++i)
254 TOA_diff+=fabs(Alist[i]-Blist[i]);
264 return TOA_diff/npul;
267 std::tuple<float, size_t> TdrtTDTReference::tdrt_recover_amp(TimeSeries<float>& h_ts)
270 size_t length = h_ts.get_length();
271 float* ts = h_ts.get_data();
278 for(
size_t i=0; i<length; ++i)
288 else if((ts[i] == 0) && (counting==1))
300 return std::make_tuple(aTOA, npul);
303 void TdrtTDTReference::tdrt_bin_spacing(TimeSeries<float>& h_ts)
305 size_t length = h_ts.get_length();
306 float* ts = h_ts.get_data();
308 int* bin_list=(
int *) malloc(length);
309 int* bin_diff=(
int *) malloc(length);
311 unsigned int npulses=0;
315 for(ii=0; ii<length; ii++)
319 bin_list[npulses]=ii;
325 for(ii=0; ii<npulses-1; ii++)
326 bin_diff[ii]=bin_list[ii+1]-bin_list[ii];
333 for(ii=0; ii<npulses-1; ii++)
335 if(bin_diff[ii+1]<bin_diff[ii])
338 bin_diff[ii]=bin_diff[ii+1];
339 bin_diff[ii+1]=holder;
348 for(ii=0; ii<npulses-1; ii++)
353 printf(
"Mean: %lf\n", mean);
357 for(ii=0; ii<npulses-1; ii++)
359 variance+=(bin_diff[ii]-mean)*(bin_diff[ii]-mean);
361 variance/=(npulses-1);
362 variance=sqrt(variance);
363 printf(
"Variance: %lf\n", variance);
366 bin_diff[npulses]=0.0;
367 for(ii=0; ii<npulses-1; ii++)
369 if(bin_diff[ii+1] != bin_diff[ii])
371 printf(
"%d: %d\n", bin_diff[ii], count);
381 double TdrtTDTReference::fold_pulsar(TimeSeries<float>& h_ts,
double period,
int adjacent)
383 int int_period = (int) round(period);
384 int bin_to_write = 0;
386 int length = h_ts.get_length();
390 double pulse_sum = 0.0;
392 float* ts = h_ts.get_data();
393 double* profile = (
double *) malloc(
sizeof(
double)*int_period);
395 for(ii=0; ii<int_period; ii++) profile[ii]=0.0;
397 std::cout <<
"Folding at a period of " << int_period <<
" bins" << std::endl;
400 for(ii=0; ii<length; ii++)
402 bin_to_write=ii%int_period;
403 profile[bin_to_write]+=ts[ii];
404 if(bin_to_write == int_period-1) num_pulses++;
407 std::cout <<
"Found " << num_pulses <<
" pulses" << std::endl;
409 for(ii=0; ii<int_period; ii++) printf(
"%lf ", profile[ii]);
413 max_bin=0; max_val=profile[0];
414 for(ii=0; ii<int_period; ii++)
416 if(profile[ii]>max_val)
424 pulse_sum=profile[max_bin];
426 for(ii=1; ii<adjacent+1; ii++)
428 if(max_bin+ii<=int_period-ii)
429 pulse_sum+=profile[max_bin+ii];
431 pulse_sum+=profile[ii-1];
434 pulse_sum+=profile[max_bin-ii];
436 pulse_sum+=profile[int_period-ii];
439 std::cout <<
"WARNING: This function only works for input periods that are integers" << std::endl;
440 std::cout <<
"Folded pulse sum: " << pulse_sum << std::endl;
444 return pulse_sum/num_pulses;
Some limits and constants for FLDO.
void process(TimeSeries< float > &h_input, TimeSeries< float > &h_output, float acceleration)