24 #include "Pwft-TDT-Reference.h" 36 PwftTDTReference::PwftTDTReference(
size_t length)
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);
44 PwftTDTReference::~PwftTDTReference()
46 fftwf_destroy_plan(_plan);
51 void PwftTDTReference::process(tdrt::TimeSeries<float>& h_input,
52 tdrt::FrequencySeries<float>& h_output)
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());
58 std::memcpy(_fft_in, h_input.get_data(),
sizeof(float) * _length);
62 float* tmp_fs= _fft_in;
63 for(
size_t i=0; i < _length/2 + 1; ++i)
65 __complex__
float tmp;
66 __real__ tmp= _fft_out[i][0];
67 __imag__ tmp= _fft_out[i][1];
68 tmp_fs[i] = cabsf(tmp);
71 float* tsOut = h_output.get_data();
72 float bprev, bnext, b0;
74 for(
size_t i=0; i < _length/2 + 1; ++i)
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);
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);
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]);
96 void PwftTDTReference::sine_signal(
float period, tdrt::TimeSeries<float>& h_output)
98 const float tsamp = h_output.get_delta_t();
99 const float pi = 3.14159265;
100 float* tsout = h_output.get_data();
102 for(
size_t i= 0; i < h_output.get_length(); i++)
104 tsout[i]= sin( 2.0*pi*i * tsamp / period );
108 float PwftTDTReference::pwft_correlation(tdrt::FrequencySeries<float>& h_fsA, tdrt::FrequencySeries<float>& h_fsB)
110 assert(h_fsA.get_length() == h_fsB.get_length());
111 assert(h_fsA.get_delta_f() == h_fsB.get_delta_f());
113 size_t length = h_fsA.get_length();
114 float* psA = h_fsA.get_data();
115 float* psB = h_fsB.get_data();
117 double correlation=0.0;
120 for(
size_t i=0; i<length; ++i)
122 correlation+=psA[i]*psB[i];
123 sumA+=(psA[i]*psA[i]);
Some limits and constants for FLDO.