Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Hrms-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 "Hrms-TDT-Reference.h"
25 // uncomment to disable assert()
26 // #define NDEBUG
27 #include <cassert>
28 #include <cstring>
29 
30 namespace ska {
31 namespace cheetah {
32 namespace hrms {
33 namespace test {
34 
35 HrmsTDTReference::HrmsTDTReference(Harmonics harmonics)
36  :_harmonics(harmonics)
37 {
38 }
39 
40 void HrmsTDTReference::process(tdrt::FrequencySeries<float>& h_input,
41  std::vector<std::shared_ptr<tdrt::FrequencySeries<float>>>& h_output)
42 {
43 
44  // std::cout << "_harmonics " << _harmonics << std::endl;
45  // std::cout << "h_output.size() " << h_output.size() << std::endl;
46 
47  unsigned n_harmonics= _harmonics+1;
48  assert(h_output.size() >= n_harmonics);
49 
50  size_t length = h_input.get_length();
51  float delta_f = h_input.get_delta_f();
52 
53  for (std::vector<std::shared_ptr<tdrt::FrequencySeries<float>>>::const_iterator it = h_output.begin();
54  it!=h_output.end(); ++it)
55  {
56  assert ((*it)->get_length() == length);
57  assert ((*it)->get_delta_f() == delta_f);
58  }
59 
60  // Reference section doing summing like sigprog sumhrm.f: -----------------------------------
61  // (omitting to just copy the first harmonic sum)
62 
63  const float Pmax = 9.9999; // largest pulsar period to find
64  size_t nf1 = 2 * length / (h_input.get_delta_f() * Pmax); // var from sigprog sumhrm.f
65 
66  //std::cout << "Pmax " << Pmax << std::endl;
67  //std::cout << "nf1 " << nf1 << std::endl;
68 
69  for (unsigned fold= 0; fold <= _harmonics; ++fold)
70  {
71  float fval = 1 << fold; // 1,2,4,8,16 ...
72  size_t NFS = fval*nf1 - fval/2; // var from sigprog sumhrm.f
73  float xdiv = 1.0 / fval; // var from sigprog sumhrm.f
74 
75  //std::cout << "summing fold " << fold << " with fval " << fval << std::endl;
76  //std::cout << "NFS " << NFS << std::endl;
77  //std::cout << "xdiv " << xdiv << std::endl;
78 
79  for (size_t n = NFS; n < length; ++n)
80  {
81  h_output[fold]->get_data()[n] = 0.0; // zero init output array
82  long long last_idx = -1;
83 
84  for (unsigned m=0; m < fval; ++m) // summing over 1,2,4,8,16 ...
85  {
86  float stretch_idx_float = n * m * xdiv + 0.5;
87  long long stretch_idx_int = (size_t) stretch_idx_float;
88  //std::cout << "stretch index " << stretch_idx_float <<" "<< stretch_idx_int << std::endl;
89  if ( stretch_idx_int > 1 && stretch_idx_int != last_idx)
90  {
91  h_output[fold]->get_data()[n] += h_input.get_data()[stretch_idx_int]; // summing
92  }
93  last_idx= stretch_idx_int;
94  }
95 
96  }
97  }
98 }
99 
100 void HrmsTDTReference::harmonics_train(size_t period, tdrt::FrequencySeries<float>& h_harmonics_train_out)
101 {
102  assert(h_harmonics_train_out.get_length() > period * 10); // accept at least 10 harmonics
103  assert(period > 32); // accept at least 32 bins between harmonics
104 
105  float* fsHarmonics = h_harmonics_train_out.get_data();
106 
107  for(size_t i = 0; i < h_harmonics_train_out.get_length(); ++i)
108  {
109  // only 16 harmonics, and DC is removed
110  if(i % period == 0 && i <= 16 * period && i > 10)
111  fsHarmonics[i]=1.0;
112  else
113  fsHarmonics[i]=0.0;
114  }
115 }
116 
117 } // test
118 } // namespace hrms
119 } // namespace cheetah
120 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35