Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Pwft-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 "Pwft-TDT-Reference.h"
25 #include <iostream>
26 // uncomment to disable assert()
27 // #define NDEBUG
28 #include <cassert>
29 #include <cstring>
30 
31 namespace ska {
32 namespace cheetah {
33 namespace pwft {
34 namespace test {
35 
36 PwftTDTReference::PwftTDTReference(size_t length)
37  : _length(length)
38 {
39  _fft_in = (float*) malloc(sizeof(float) * _length);
40  _fft_out = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * (_length/2+1));
41  _plan = fftwf_plan_dft_r2c_1d(_length, _fft_in, _fft_out, FFTW_MEASURE);
42 }
43 
44 PwftTDTReference::~PwftTDTReference()
45 {
46  fftwf_destroy_plan(_plan);
47  fftwf_free(_fft_in);
48  free(_fft_out);
49 }
50 
51 void PwftTDTReference::process(tdrt::TimeSeries<float>& h_input,
52  tdrt::FrequencySeries<float>& h_output)
53 {
54  assert(h_input.get_length() == _length);
55  assert(h_output.get_length() == _length);
56  assert( 1.0 / h_input.get_delta_t() == h_output.get_delta_f());
57 
58  std::memcpy(_fft_in, h_input.get_data(), sizeof(float) * _length);
59 
60  fftwf_execute(_plan);
61 
62  float* tmp_fs= _fft_in;
63  for(size_t i=0; i < _length/2 + 1; ++i)
64  {
65  __complex__ float tmp;
66  __real__ tmp= _fft_out[i][0];
67  __imag__ tmp= _fft_out[i][1];
68  tmp_fs[i] = cabsf(tmp);
69  }
70 
71  float* tsOut = h_output.get_data();
72  float bprev, bnext, b0;
73 
74  for(size_t i=0; i < _length/2 + 1; ++i)
75  {
76  if(i==0)
77  {
78  bnext=0.5*powf(fabsf(tmp_fs[i]+tmp_fs[i+1]), 2);
79  tsOut[i]=fmaxf(tmp_fs[i]*tmp_fs[i], bnext);
80  }
81  else if(i==_length/2)
82  {
83  bprev=0.5*powf(fabsf(tmp_fs[i]+tmp_fs[i-1]), 2);
84  tsOut[i]=fmaxf(tmp_fs[i]*tmp_fs[i], bprev);
85  }
86  else
87  {
88  bprev=0.5*powf(fabsf(tmp_fs[i]+tmp_fs[i-1]), 2);
89  bnext=0.5*powf(fabsf(tmp_fs[i]+tmp_fs[i+1]), 2);
90  b0=fmaxf(bprev, bnext);
91  tsOut[i]=fmax(b0, tmp_fs[i]*tmp_fs[i]);
92  }
93  }
94 }
95 
96 void PwftTDTReference::sine_signal(float period, tdrt::TimeSeries<float>& h_output)
97 {
98  const float tsamp = h_output.get_delta_t();
99  const float pi = 3.14159265;
100  float* tsout = h_output.get_data();
101 
102  for(size_t i= 0; i < h_output.get_length(); i++)
103  {
104  tsout[i]= sin( 2.0*pi*i * tsamp / period );
105  }
106 }
107 
108 float PwftTDTReference::pwft_correlation(tdrt::FrequencySeries<float>& h_fsA, tdrt::FrequencySeries<float>& h_fsB)
109 {
110  assert(h_fsA.get_length() == h_fsB.get_length());
111  assert(h_fsA.get_delta_f() == h_fsB.get_delta_f());
112 
113  size_t length = h_fsA.get_length();
114  float* psA = h_fsA.get_data();
115  float* psB = h_fsB.get_data();
116 
117  double correlation=0.0;
118  double sumA=0.0;
119 
120  for(size_t i=0; i<length; ++i)
121  {
122  correlation+=psA[i]*psB[i];
123  sumA+=(psA[i]*psA[i]);
124  }
125  //std::cout << "Correlation: " << correlation << std::endl;
126  //std::cout << "SumA: " << sumA << std::endl;
127  //std::cout << "Normalized correlation: " << correlation/sumA << std::endl;
128 
129  if(sumA != 0.0)
130  {
131  correlation/=sumA;
132  }
133  else
134  {
135  correlation=0.0;
136  }
137  return correlation;
138 }
139 
140 } // namespace test
141 } // namespace pwft
142 } // namespace cheetah
143 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35