Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Tdrt-TDT-Reference.cpp
1 /*
2  * The MIT License (MIT)
3  *
4  * Copyright (c) 2016 The SKA organisation
5  *
6  * Permission is hereby granted, free of charge, to any person obtaining a copy
7  * of this software and associated documentation files (the "Software"), to deal
8  * in the Software without restriction, including without limitation the rights
9  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10  * copies of the Software, and to permit persons to whom the Software is
11  * furnished to do so, subject to the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be included in all
14  * copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
24 #include "Tdrt-TDT-Reference.h"
25 #include <iostream>
26 #include <cmath>
27 // uncomment to disable assert()
28 // #define NDEBUG
29 #include <cassert>
30 
31 namespace ska {
32 namespace cheetah {
33 namespace tdrt {
34 namespace test {
35 
36 TdrtTDTReference::TdrtTDTReference()
37 {
38 }
39 
40 void TdrtTDTReference::process(TimeSeries<float>& h_input,
41  TimeSeries<float>& h_output,
42  float acceleration)
43 {
44  assert(h_input.get_length() == h_output.get_length());
45  assert(h_input.get_delta_t() == h_output.get_delta_t());
46 
47  long long numTimeSamps= h_input.get_length();;
48  double sampTime= h_input.get_delta_t();
49  double pulsarTime= 0.0;
50  double timeThru= 0.0;
51  double earthTime= 0.0;
52 
53  const double c = 299792458.0; // speed of light
54 
55  float obsLen = sampTime * numTimeSamps;
56  float tau0 = sampTime/(1 + acceleration * obsLen/(2 * c));
57  // float earlyWeight, lateWeight;
58 
59  long long nsamp;
60  long long samp_out;
61 
62  float* input = h_input.get_data();
63  float* output = h_output.get_data();
64  for (nsamp=0; nsamp < numTimeSamps; nsamp++)
65  {
66  // Time interval on Earth
67  earthTime = nsamp * sampTime;
68 
69  timeThru = pulsarTime;
70 
71  // Relate the time interval at Earth to that of the pulsar
72  pulsarTime = pulsarTime + tau0 * (1 + acceleration * earthTime/c);
73 
74  // samp_out = (int) floor(timeThru/sampTime);
75  samp_out = (int) round(timeThru/sampTime);
76 
77  // earlyWeight = (timeThru - samp_out * sampTime)/sampTime;
79 
80  if (samp_out >= 0 && samp_out < numTimeSamps)
81  {
82  // output[samp_out] += lateWeight * input[nsamp];
83  output[samp_out] = input[nsamp];
84  }
85  /*
86  if (samp_out >= -1 && samp_out < (numTimeSamps - 1))
87  {
88  output[samp_out+1] = earlyWeight * input[nsamp];
89  }
90  */
91  }
92 }
93 
94 void TdrtTDTReference::const_accel_pulsar(double acceleration, float period,
95  TimeSeries<float>& h_output)
96 {
97  const double c = 299792458.0; // speed of light
98  const double tsamp = h_output.get_delta_t();
99 
100  float f0= 1.0 / period;
101  float fdot= acceleration / c * f0;
102  float Nrot, Nrot_next;
103 
104  float* tsout = h_output.get_data();
105 
106  //Initialize tdum array, calc Nrot and Nrotint
107 
108  // for(size_t i= 0; i < h_output.get_length() - 1; i++)
109  for(size_t i= 0; i < h_output.get_length(); i++)
110  {
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);
113 
114  // float floc;
115  // float binfrac;
116  if( ((int) floor(Nrot_next) - (int) floor(Nrot)) == 1)
117  {
118  // floc = f0 + fdot * tsamp * i;
119  tsout[i]= 1.0;
120  /*binfrac=fabs(Nrot-lround(Nrot))/floc/tsamp;
121  tsout[i]=binfrac;
122  tsout[i+1]=1-binfrac; */
123  }
124  else
125  {
126  tsout[i]=0.0;
127  }
128  }
129 }
130 
131 void TdrtTDTReference::sine_calibration_signal(float period, TimeSeries<float>& h_output)
132 {
133  const float tsamp = h_output.get_delta_t();
134  const float pi = 3.14159265;
135  float* tsout = h_output.get_data();
136 
137  for(size_t i= 0; i < h_output.get_length(); i++)
138  {
139  tsout[i]= sin( 2.0*pi*i * tsamp / period );
140  }
141 }
142 
143 double TdrtTDTReference::resampled_period(double acceleration, double period, double tsamp, int length)
144 {
145  double c = 299792458.0; // speed of light (m/s)
146  double new_period;
147 
148  new_period = period/(1+(acceleration*tsamp*length)/(2*c));
149 
150  std::cout << "Period after resampling: " << new_period << std::endl;
151 
152  return new_period;
153 }
154 
155 float TdrtTDTReference::tdrt_correlation(TimeSeries<float>& h_tsA, TimeSeries<float>& h_tsB)
156 {
157  assert(h_tsA.get_length() == h_tsB.get_length());
158  assert(h_tsA.get_delta_t() == h_tsB.get_delta_t());
159 
160  size_t length = h_tsA.get_length();
161  float* tsA = h_tsA.get_data();
162  float* tsB = h_tsB.get_data();
163  size_t npul = 0;
164  size_t found = 0;
165  double correlation = 0.0;
166 
167  for(size_t i=0; i<length; ++i)
168  {
169  correlation += tsA[i] * tsB[i];
170 
171  /*This code is only needed to count pulses*/
172  if((tsA[i] != 0) && (found==0))
173  {
174  npul+=1;
175  found=1;
176  }
177  if((tsA[i] == 0) && (found==1))
178  {
179  found=0;
180  }
181  }
182 
183  //std::cout << "Correlation: " << correlation << std::endl;
184  //std::cout << "Number of pulses: " << npul << std::endl;
185  //std::cout << "Normalized correlation: " << correlation/npul << std::endl;
186 
187  return correlation/npul;
188 }
189 
190 float TdrtTDTReference::tdrt_bin_diff(TimeSeries<float>& h_tsA, TimeSeries<float>& h_tsB)
191 {
192  assert(h_tsA.get_length() == h_tsB.get_length());
193  assert(h_tsA.get_delta_t() == h_tsB.get_delta_t());
194 
195  size_t length = h_tsA.get_length();
196  float* tsA = h_tsA.get_data();
197  float* tsB = h_tsB.get_data();
198 
199  size_t i = 0;
200  size_t npulA=0;
201  size_t npulB=0;
202  size_t npul=0;
203  size_t counting=0;
204 
205  float TOA_A=0.0;
206  float TOA_B=0.0;
207  float TOA_diff=0.0;
208 
209  float amp=0;
210  float* Alist = (float*)malloc(sizeof(float)*length/2);
211  float* Blist = (float*)malloc(sizeof(float)*length/2);
212 
213  for(i=0; i<length; ++i)
214  {
215  if(tsA[i] > 0)
216  {
217  TOA_A+=tsA[i]*i;
218  amp+=tsA[i];
219  counting=1;
220  }
221  else if((tsA[i] == 0) && (counting == 1))
222  {
223  Alist[npulA]=TOA_A/amp;
224  amp=0;
225  TOA_A=0;
226  npulA+=1;
227  counting=0;
228  }
229  }
230 
231  for(i=0; i<length; ++i)
232  {
233  if(tsB[i] > 0)
234  {
235  TOA_B+=tsB[i]*i;
236  amp+=tsB[i];
237  counting=1;
238  }
239  else if((tsB[i] == 0) && (counting == 1))
240  {
241  Blist[npulB]=TOA_B/amp;
242  amp=0;
243  TOA_B=0;
244  npulB+=1;
245  counting=0;
246  }
247  }
248 
249  if(npulA<npulB) npul=npulA;
250  else npul=npulB;
251 
252  for(i=0; i<npul; ++i)
253  {
254  TOA_diff+=fabs(Alist[i]-Blist[i]);
255  }
256 
257  //std::cout << "Num pulses: " << npul << std::endl;
258  //std::cout << "TOA diff: " << TOA_diff << std::endl;
259  //std::cout << "Norm TOA diff: " << TOA_diff/npul << std::endl;
260 
261  free(Alist);
262  free(Blist);
263 
264  return TOA_diff/npul;
265 }
266 
267 std::tuple<float, size_t> TdrtTDTReference::tdrt_recover_amp(TimeSeries<float>& h_ts)
268 {
269 
270  size_t length = h_ts.get_length();
271  float* ts = h_ts.get_data();
272 
273  size_t counting=0;
274  size_t npul=1;
275 
276  float aTOA = 0.0;
277 
278  for(size_t i=0; i<length; ++i)
279  {
280  if(ts[i] > 0)
281  {
282  if(i==0)
283  npul=0;
284 
285  aTOA+=ts[i];
286  counting=1;
287  }
288  else if((ts[i] == 0) && (counting==1))
289  {
290  npul+=1;
291  counting=0;
292  }
293  }
294 
295  aTOA/=npul;
296 
297  //std::cout << "Number of pulses: " << npul << std::endl;
298  //std::cout << "Amplitude recovered: " << aTOA << std::endl;
299 
300  return std::make_tuple(aTOA, npul);
301 }
302 
303 void TdrtTDTReference::tdrt_bin_spacing(TimeSeries<float>& h_ts)
304 {
305  size_t length = h_ts.get_length();
306  float* ts = h_ts.get_data();
307 
308  int* bin_list=(int *) malloc(length);
309  int* bin_diff=(int *) malloc(length);
310  unsigned int ii=0;
311  unsigned int npulses=0;
312  int holder=0;
313 
314  /*First loop over time series and find locations of pulses*/
315  for(ii=0; ii<length; ii++)
316  {
317  if(ts[ii] != 0.0)
318  {
319  bin_list[npulses]=ii;
320  npulses++;
321  }
322  }
323 
324  /*Calculate differences*/
325  for(ii=0; ii<npulses-1; ii++)
326  bin_diff[ii]=bin_list[ii+1]-bin_list[ii];
327 
328  /*Sort*/
329  int noflip=0;
330  while(noflip == 0)
331  {
332  noflip=1;
333  for(ii=0; ii<npulses-1; ii++)
334  {
335  if(bin_diff[ii+1]<bin_diff[ii])
336  {
337  holder=bin_diff[ii];
338  bin_diff[ii]=bin_diff[ii+1];
339  bin_diff[ii+1]=holder;
340  noflip=0;
341  }
342  }
343  }
344 
345  /*Calc mean */
346  double mean=0.0;
347 
348  for(ii=0; ii<npulses-1; ii++)
349  {
350  mean+=bin_diff[ii];
351  }
352  mean/=(npulses-1);
353  printf("Mean: %lf\n", mean);
354 
355  /* Calc stddev */
356  double variance=0.0;
357  for(ii=0; ii<npulses-1; ii++)
358  {
359  variance+=(bin_diff[ii]-mean)*(bin_diff[ii]-mean);
360  }
361  variance/=(npulses-1);
362  variance=sqrt(variance);
363  printf("Variance: %lf\n", variance);
364 
365  int count=1;
366  bin_diff[npulses]=0.0; /*This is a hack to allow loop to go over edge*/
367  for(ii=0; ii<npulses-1; ii++)
368  {
369  if(bin_diff[ii+1] != bin_diff[ii])
370  {
371  printf("%d: %d\n", bin_diff[ii], count);
372  count=1;
373  }
374  else
375  {
376  count++;
377  }
378  }
379 }
380 
381 double TdrtTDTReference::fold_pulsar(TimeSeries<float>& h_ts, double period, int adjacent)
382 {
383  int int_period = (int) round(period);
384  int bin_to_write = 0;
385  int num_pulses = 0;
386  int length = h_ts.get_length();
387  int ii=0;
388  int max_bin=-1;
389  double max_val=-1.0;
390  double pulse_sum = 0.0;
391 
392  float* ts = h_ts.get_data();
393  double* profile = (double *) malloc(sizeof(double)*int_period);
394 
395  for(ii=0; ii<int_period; ii++) profile[ii]=0.0;
396 
397  std::cout << "Folding at a period of " << int_period << " bins" << std::endl;
398 
399  /* Fold profile */
400  for(ii=0; ii<length; ii++)
401  {
402  bin_to_write=ii%int_period;
403  profile[bin_to_write]+=ts[ii];
404  if(bin_to_write == int_period-1) num_pulses++;
405  }
406 
407  std::cout << "Found " << num_pulses << " pulses" << std::endl;
408 
409  for(ii=0; ii<int_period; ii++) printf("%lf ", profile[ii]);
410  printf("\n");
411 
412  /*Find max bin*/
413  max_bin=0; max_val=profile[0];
414  for(ii=0; ii<int_period; ii++)
415  {
416  if(profile[ii]>max_val)
417  {
418  max_val=profile[ii];
419  max_bin=ii;
420  }
421  }
422 
423  /* Sum profile including adjecent bins set by input param */
424  pulse_sum=profile[max_bin];
425 
426  for(ii=1; ii<adjacent+1; ii++)
427  {
428  if(max_bin+ii<=int_period-ii) /*Check for wrapping*/
429  pulse_sum+=profile[max_bin+ii];
430  else
431  pulse_sum+=profile[ii-1];
432 
433  if(max_bin-ii>=0)
434  pulse_sum+=profile[max_bin-ii];
435  else
436  pulse_sum+=profile[int_period-ii];
437  }
438 
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;
441 
442  free(profile);
443 
444  return pulse_sum/num_pulses;
445 }
446 
447 } // namespace test
448 } // namespace tdrt
449 } // namespace cheetah
450 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35
void process(TimeSeries< float > &h_input, TimeSeries< float > &h_output, float acceleration)