1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
//
// Copyright 2010 Ettus Research LLC
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
#include <vector>
#include <complex>
#include <cmath>
#include <fstream>
struct result_t{double freq, real_corr, imag_corr, sup;};
/***********************************************************************
* Constants
**********************************************************************/
static const double tau = 6.28318531;
static const double alpha = 0.0001; //very tight iir filter
static const size_t wave_table_len = 8192;
static const size_t num_search_steps = 5;
static const size_t num_search_iters = 7;
static const size_t skip_initial_samps = 20;
/***********************************************************************
* Sinusoid wave table
**********************************************************************/
static inline std::vector<std::complex<float> > gen_table(void){
std::vector<std::complex<float> > wave_table(wave_table_len);
for (size_t i = 0; i < wave_table_len; i++){
wave_table[i] = std::polar<float>(1.0, (tau*i)/wave_table_len);
}
return wave_table;
}
static inline std::complex<float> wave_table_lookup(const size_t index){
static const std::vector<std::complex<float> > wave_table = gen_table();
return wave_table[index % wave_table_len];
}
/***********************************************************************
* Compute power of a tone
**********************************************************************/
static inline double compute_tone_dbrms(
const std::vector<std::complex<float> > &samples,
const double freq //freq is fractional
){
//shift the samples so the tone at freq is down at DC
std::vector<std::complex<double> > shifted(samples.size() - skip_initial_samps);
for (size_t i = 0; i < shifted.size(); i++){
shifted[i] = std::complex<double>(samples[i+skip_initial_samps]) * std::polar<double>(1.0, -freq*tau*i);
}
//filter the samples with a narrow low pass
std::complex<double> iir_output = 0, iir_last = 0;
double output = 0;
for (size_t i = 0; i < shifted.size(); i++){
iir_output = alpha * shifted[i] + (1-alpha)*iir_last;
iir_last = iir_output;
output += std::abs(iir_output);
}
return 20*std::log10(output/shifted.size());
}
/***********************************************************************
* Write a dat file
**********************************************************************/
static inline void write_samples_to_file(
const std::vector<std::complex<float> > &samples, const std::string &file
){
std::ofstream outfile(file.c_str(), std::ofstream::binary);
outfile.write((const char*)&samples.front(), samples.size()*sizeof(std::complex<float>));
outfile.close();
}
|
