mirrors/gnss-sdr
1
0
Fork 0
mirror of https://github.com/gnss-sdr/gnss-sdr synced 2025-11-18 16:15:21 +00:00
Code Issues Releases Wiki Activity

Merge branch 'next' of https://github.com/gnss-sdr/gnss-sdr into glamountain-kf

Browse Source
This commit is contained in:
Carles Fernandez
2018-07-20 12:59:46 +02:00
parent 486ac195db 23a036f58f
commit 5a908fcf7d
63 changed files with 1948 additions and 1122 deletions
Download Patch File Download Diff File Expand all files Collapse all files

363
src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition.cc
View File

@@ -52,8 +52,8 @@ pcps_acquisition_sptr pcps_make_acquisition(const Acq_Conf& conf_)
pcps_acquisition::pcps_acquisition(const Acq_Conf& conf_) : gr::block("pcps_acquisition",
gr::io_signature::make(1, 1, conf_.it_size * conf_.sampled_ms * conf_.samples_per_ms * (conf_.bit_transition_flag ? 2 : 1)),
gr::io_signature::make(0, 0, conf_.it_size * conf_.sampled_ms * conf_.samples_per_ms * (conf_.bit_transition_flag ? 2 : 1)))
gr::io_signature::make(1, 1, conf_.it_size * std::floor(conf_.sampled_ms * conf_.samples_per_ms) * (conf_.bit_transition_flag ? 2 : 1)),
gr::io_signature::make(0, 0, conf_.it_size))
{
this->message_port_register_out(pmt::mp("events"));
@@ -63,8 +63,17 @@ pcps_acquisition::pcps_acquisition(const Acq_Conf& conf_) : gr::block("pcps_acqu
d_positive_acq = 0;
d_state = 0;
d_old_freq = 0;
d_well_count = 0;
d_fft_size = acq_parameters.sampled_ms * acq_parameters.samples_per_ms;
d_num_noncoherent_integrations_counter = 0;
d_consumed_samples = acq_parameters.sampled_ms * acq_parameters.samples_per_ms * (acq_parameters.bit_transition_flag ? 2 : 1);
if (acq_parameters.sampled_ms == acq_parameters.ms_per_code)
{
d_fft_size = d_consumed_samples;
}
else
{
d_fft_size = d_consumed_samples * 2;
}
//d_fft_size = next power of two? ////
d_mag = 0;
d_input_power = 0.0;
d_num_doppler_bins = 0;
@@ -94,12 +103,14 @@ pcps_acquisition::pcps_acquisition(const Acq_Conf& conf_) : gr::block("pcps_acqu
// size of the input buffer and padding the code with zeros.
if (acq_parameters.bit_transition_flag)
{
d_fft_size *= 2;
acq_parameters.max_dwells = 1; //Activation of acq_parameters.bit_transition_flag invalidates the value of acq_parameters.max_dwells
d_fft_size = d_consumed_samples * 2;
acq_parameters.max_dwells = 1; // Activation of acq_parameters.bit_transition_flag invalidates the value of acq_parameters.max_dwells
}
d_tmp_buffer = static_cast<float*>(volk_gnsssdr_malloc(d_fft_size * sizeof(float), volk_gnsssdr_get_alignment()));
d_fft_codes = static_cast<gr_complex*>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
d_magnitude = static_cast<float*>(volk_gnsssdr_malloc(d_fft_size * sizeof(float), volk_gnsssdr_get_alignment()));
d_input_signal = static_cast<gr_complex*>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
// Direct FFT
d_fft_if = new gr::fft::fft_complex(d_fft_size, true);
@@ -110,11 +121,12 @@ pcps_acquisition::pcps_acquisition(const Acq_Conf& conf_) : gr::block("pcps_acqu
d_gnss_synchro = 0;
d_grid_doppler_wipeoffs = nullptr;
d_grid_doppler_wipeoffs_step_two = nullptr;
d_magnitude_grid = nullptr;
d_worker_active = false;
d_data_buffer = static_cast<gr_complex*>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
d_data_buffer = static_cast<gr_complex*>(volk_gnsssdr_malloc(d_consumed_samples * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
if (d_cshort)
{
d_data_buffer_sc = static_cast<lv_16sc_t*>(volk_gnsssdr_malloc(d_fft_size * sizeof(lv_16sc_t), volk_gnsssdr_get_alignment()));
d_data_buffer_sc = static_cast<lv_16sc_t*>(volk_gnsssdr_malloc(d_consumed_samples * sizeof(lv_16sc_t), volk_gnsssdr_get_alignment()));
}
else
{
@@ -124,6 +136,16 @@ pcps_acquisition::pcps_acquisition(const Acq_Conf& conf_) : gr::block("pcps_acqu
d_step_two = false;
d_dump_number = 0;
d_dump_channel = acq_parameters.dump_channel;
d_samplesPerChip = acq_parameters.samples_per_chip;
// todo: CFAR statistic not available for non-coherent integration
if (acq_parameters.max_dwells == 1)
{
d_use_CFAR_algorithm_flag = acq_parameters.use_CFAR_algorithm_flag;
}
else
{
d_use_CFAR_algorithm_flag = false;
}
}
@@ -134,8 +156,10 @@ pcps_acquisition::~pcps_acquisition()
for (unsigned int i = 0; i < d_num_doppler_bins; i++)
{
volk_gnsssdr_free(d_grid_doppler_wipeoffs[i]);
volk_gnsssdr_free(d_magnitude_grid[i]);
}
delete[] d_grid_doppler_wipeoffs;
delete[] d_magnitude_grid;
}
if (acq_parameters.make_2_steps)
{
@@ -147,6 +171,8 @@ pcps_acquisition::~pcps_acquisition()
}
volk_gnsssdr_free(d_fft_codes);
volk_gnsssdr_free(d_magnitude);
volk_gnsssdr_free(d_tmp_buffer);
volk_gnsssdr_free(d_input_signal);
delete d_ifft;
delete d_fft_if;
volk_gnsssdr_free(d_data_buffer);
@@ -179,7 +205,15 @@ void pcps_acquisition::set_local_code(std::complex<float>* code)
}
else
{
memcpy(d_fft_if->get_inbuf(), code, sizeof(gr_complex) * d_fft_size);
if (acq_parameters.sampled_ms == acq_parameters.ms_per_code)
{
memcpy(d_fft_if->get_inbuf(), code, sizeof(gr_complex) * d_consumed_samples);
}
else
{
std::fill_n(d_fft_if->get_inbuf(), d_fft_size - d_consumed_samples, gr_complex(0.0, 0.0));
memcpy(d_fft_if->get_inbuf() + d_consumed_samples, code, sizeof(gr_complex) * d_consumed_samples);
}
}
d_fft_if->execute(); // We need the FFT of local code
@@ -244,12 +278,20 @@ void pcps_acquisition::init()
d_grid_doppler_wipeoffs_step_two[doppler_index] = static_cast<gr_complex*>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
}
}
d_magnitude_grid = new float*[d_num_doppler_bins];
for (unsigned int doppler_index = 0; doppler_index < d_num_doppler_bins; doppler_index++)
{
d_grid_doppler_wipeoffs[doppler_index] = static_cast<gr_complex*>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
d_magnitude_grid[doppler_index] = static_cast<float*>(volk_gnsssdr_malloc(d_fft_size * sizeof(float), volk_gnsssdr_get_alignment()));
for (unsigned k = 0; k < d_fft_size; k++)
{
d_magnitude_grid[doppler_index][k] = 0.0;
}
int doppler = -static_cast<int>(acq_parameters.doppler_max) + d_doppler_step * doppler_index;
update_local_carrier(d_grid_doppler_wipeoffs[doppler_index], d_fft_size, d_old_freq + doppler);
}
d_worker_active = false;
if (acq_parameters.dump)
@@ -269,6 +311,7 @@ void pcps_acquisition::update_grid_doppler_wipeoffs()
}
}
void pcps_acquisition::update_grid_doppler_wipeoffs_step2()
{
for (unsigned int doppler_index = 0; doppler_index < acq_parameters.num_doppler_bins_step2; doppler_index++)
@@ -278,6 +321,7 @@ void pcps_acquisition::update_grid_doppler_wipeoffs_step2()
}
}
void pcps_acquisition::set_state(int state)
{
gr::thread::scoped_lock lock(d_setlock); // require mutex with work function called by the scheduler
@@ -288,7 +332,6 @@ void pcps_acquisition::set_state(int state)
d_gnss_synchro->Acq_doppler_hz = 0.0;
d_gnss_synchro->Acq_doppler_step = 0;
d_gnss_synchro->Acq_samplestamp_samples = 0;
d_well_count = 0;
d_mag = 0.0;
d_input_power = 0.0;
d_test_statistics = 0.0;
@@ -419,117 +462,202 @@ void pcps_acquisition::dump_results(int effective_fft_size)
}
float pcps_acquisition::max_to_input_power_statistic(uint32_t& indext, int& doppler, float input_power, unsigned int num_doppler_bins, int doppler_max, int doppler_step)
{
float grid_maximum = 0.0;
unsigned int index_doppler = 0;
uint32_t tmp_intex_t = 0;
uint32_t index_time = 0;
float fft_normalization_factor = static_cast<float>(d_fft_size) * static_cast<float>(d_fft_size);
// Find the correlation peak and the carrier frequency
for (unsigned int i = 0; i < num_doppler_bins; i++)
{
volk_gnsssdr_32f_index_max_32u(&tmp_intex_t, d_magnitude_grid[i], d_fft_size);
if (d_magnitude_grid[i][tmp_intex_t] > grid_maximum)
{
grid_maximum = d_magnitude_grid[i][tmp_intex_t];
index_doppler = i;
index_time = tmp_intex_t;
}
}
indext = index_time;
if (!d_step_two)
{
doppler = -static_cast<int>(doppler_max) + doppler_step * static_cast<int>(index_doppler);
}
else
{
doppler = static_cast<int>(d_doppler_center_step_two + (index_doppler - (acq_parameters.num_doppler_bins_step2 / 2.0) * acq_parameters.doppler_step2));
}
float magt = grid_maximum / (fft_normalization_factor * fft_normalization_factor);
return magt / input_power;
}
float pcps_acquisition::first_vs_second_peak_statistic(uint32_t& indext, int& doppler, unsigned int num_doppler_bins, int doppler_max, int doppler_step)
{
// Look for correlation peaks in the results
// Find the highest peak and compare it to the second highest peak
// The second peak is chosen not closer than 1 chip to the highest peak
float firstPeak = 0.0;
unsigned int index_doppler = 0;
uint32_t tmp_intex_t = 0;
uint32_t index_time = 0;
// Find the correlation peak and the carrier frequency
for (unsigned int i = 0; i < num_doppler_bins; i++)
{
volk_gnsssdr_32f_index_max_32u(&tmp_intex_t, d_magnitude_grid[i], d_fft_size);
if (d_magnitude_grid[i][tmp_intex_t] > firstPeak)
{
firstPeak = d_magnitude_grid[i][tmp_intex_t];
index_doppler = i;
index_time = tmp_intex_t;
}
}
indext = index_time;
if (!d_step_two)
{
doppler = -static_cast<int>(doppler_max) + doppler_step * static_cast<int>(index_doppler);
}
else
{
doppler = static_cast<int>(d_doppler_center_step_two + (index_doppler - (acq_parameters.num_doppler_bins_step2 / 2.0) * acq_parameters.doppler_step2));
}
// Find 1 chip wide code phase exclude range around the peak
int32_t excludeRangeIndex1 = index_time - d_samplesPerChip;
int32_t excludeRangeIndex2 = index_time + d_samplesPerChip;
// Correct code phase exclude range if the range includes array boundaries
if (excludeRangeIndex1 < 0)
{
excludeRangeIndex1 = d_fft_size + excludeRangeIndex1;
}
else if (excludeRangeIndex2 >= static_cast<int>(d_fft_size))
{
excludeRangeIndex2 = excludeRangeIndex2 - d_fft_size;
}
int32_t idx = excludeRangeIndex1;
memcpy(d_tmp_buffer, d_magnitude_grid[index_doppler], d_fft_size);
do
{
d_tmp_buffer[idx] = 0.0;
idx++;
if (idx == static_cast<int>(d_fft_size)) idx = 0;
}
while (idx != excludeRangeIndex2);
// Find the second highest correlation peak in the same freq. bin ---
volk_gnsssdr_32f_index_max_32u(&tmp_intex_t, d_tmp_buffer, d_fft_size);
float secondPeak = d_tmp_buffer[tmp_intex_t];
// Compute the test statistics and compare to the threshold
return firstPeak / secondPeak;
}
void pcps_acquisition::acquisition_core(unsigned long int samp_count)
{
gr::thread::scoped_lock lk(d_setlock);
// initialize acquisition algorithm
int doppler = 0;
uint32_t indext = 0;
float magt = 0.0;
const gr_complex* in = d_data_buffer; // Get the input samples pointer
int effective_fft_size = (acq_parameters.bit_transition_flag ? d_fft_size / 2 : d_fft_size);
if (d_cshort)
{
volk_gnsssdr_16ic_convert_32fc(d_data_buffer, d_data_buffer_sc, d_fft_size);
volk_gnsssdr_16ic_convert_32fc(d_data_buffer, d_data_buffer_sc, d_consumed_samples);
}
float fft_normalization_factor = static_cast<float>(d_fft_size) * static_cast<float>(d_fft_size);
memcpy(d_input_signal, d_data_buffer, d_consumed_samples * sizeof(gr_complex));
if (d_fft_size > d_consumed_samples)
{
for (unsigned int i = d_consumed_samples; i < d_fft_size; i++)
{
d_input_signal[i] = gr_complex(0.0, 0.0);
}
}
const gr_complex* in = d_input_signal; // Get the input samples pointer
d_input_power = 0.0;
d_mag = 0.0;
d_well_count++;
d_num_noncoherent_integrations_counter++;
DLOG(INFO) << "Channel: " << d_channel
<< " , doing acquisition of satellite: " << d_gnss_synchro->System << " " << d_gnss_synchro->PRN
<< " ,sample stamp: " << samp_count << ", threshold: "
<< d_threshold << ", doppler_max: " << acq_parameters.doppler_max
<< ", doppler_step: " << d_doppler_step
<< ", use_CFAR_algorithm_flag: " << (acq_parameters.use_CFAR_algorithm_flag ? "true" : "false");
<< ", use_CFAR_algorithm_flag: " << (d_use_CFAR_algorithm_flag ? "true" : "false");
lk.unlock();
if (acq_parameters.use_CFAR_algorithm_flag)
if (d_use_CFAR_algorithm_flag or acq_parameters.bit_transition_flag)
{
// 1- (optional) Compute the input signal power estimation
volk_32fc_magnitude_squared_32f(d_magnitude, in, d_fft_size);
volk_32f_accumulator_s32f(&d_input_power, d_magnitude, d_fft_size);
// Compute the input signal power estimation
volk_32fc_magnitude_squared_32f(d_tmp_buffer, in, d_fft_size);
volk_32f_accumulator_s32f(&d_input_power, d_tmp_buffer, d_fft_size);
d_input_power /= static_cast<float>(d_fft_size);
}
// 2- Doppler frequency search loop
// Doppler frequency grid loop
if (!d_step_two)
{
for (unsigned int doppler_index = 0; doppler_index < d_num_doppler_bins; doppler_index++)
{
// doppler search steps
int doppler = -static_cast<int>(acq_parameters.doppler_max) + d_doppler_step * doppler_index;
// Remove Doppler
volk_32fc_x2_multiply_32fc(d_fft_if->get_inbuf(), in, d_grid_doppler_wipeoffs[doppler_index], d_fft_size);
// 3- Perform the FFT-based convolution (parallel time search)
// Perform the FFT-based convolution (parallel time search)
// Compute the FFT of the carrier wiped--off incoming signal
d_fft_if->execute();
// Multiply carrier wiped--off, Fourier transformed incoming signal
// with the local FFT'd code reference using SIMD operations with VOLK library
// Multiply carrier wiped--off, Fourier transformed incoming signal with the local FFT'd code reference
volk_32fc_x2_multiply_32fc(d_ifft->get_inbuf(), d_fft_if->get_outbuf(), d_fft_codes, d_fft_size);
// compute the inverse FFT
// Compute the inverse FFT
d_ifft->execute();
// Search maximum
// Compute squared magnitude (and accumulate in case of non-coherent integration)
size_t offset = (acq_parameters.bit_transition_flag ? effective_fft_size : 0);
volk_32fc_magnitude_squared_32f(d_magnitude, d_ifft->get_outbuf() + offset, effective_fft_size);
volk_gnsssdr_32f_index_max_32u(&indext, d_magnitude, effective_fft_size);
magt = d_magnitude[indext];
if (acq_parameters.use_CFAR_algorithm_flag)
if (d_num_noncoherent_integrations_counter == 1)
{
// Normalize the maximum value to correct the scale factor introduced by FFTW
magt = d_magnitude[indext] / (fft_normalization_factor * fft_normalization_factor);
volk_32fc_magnitude_squared_32f(d_magnitude_grid[doppler_index], d_ifft->get_outbuf() + offset, effective_fft_size);
}
// 4- record the maximum peak and the associated synchronization parameters
if (d_mag < magt)
else
{
d_mag = magt;
if (!acq_parameters.use_CFAR_algorithm_flag)
{
// Search grid noise floor approximation for this doppler line
volk_32f_accumulator_s32f(&d_input_power, d_magnitude, effective_fft_size);
d_input_power = (d_input_power - d_mag) / (effective_fft_size - 1);
}
// In case that acq_parameters.bit_transition_flag = true, we compare the potentially
// new maximum test statistics (d_mag/d_input_power) with the value in
// d_test_statistics. When the second dwell is being processed, the value
// of d_mag/d_input_power could be lower than d_test_statistics (i.e,
// the maximum test statistics in the previous dwell is greater than
// current d_mag/d_input_power). Note that d_test_statistics is not
// restarted between consecutive dwells in multidwell operation.
if (d_test_statistics < (d_mag / d_input_power) or !acq_parameters.bit_transition_flag)
{
d_gnss_synchro->Acq_delay_samples = static_cast<double>(indext % acq_parameters.samples_per_code);
d_gnss_synchro->Acq_doppler_hz = static_cast<double>(doppler);
d_gnss_synchro->Acq_samplestamp_samples = samp_count;
// 5- Compute the test statistics and compare to the threshold
//d_test_statistics = 2 * d_fft_size * d_mag / d_input_power;
d_test_statistics = d_mag / d_input_power;
}
volk_32fc_magnitude_squared_32f(d_tmp_buffer, d_ifft->get_outbuf() + offset, effective_fft_size);
volk_32f_x2_add_32f(d_magnitude_grid[doppler_index], d_magnitude_grid[doppler_index], d_tmp_buffer, effective_fft_size);
}
// Record results to file if required
if (acq_parameters.dump and d_channel == d_dump_channel)
{
memcpy(grid_.colptr(doppler_index), d_magnitude, sizeof(float) * effective_fft_size);
memcpy(grid_.colptr(doppler_index), d_magnitude_grid[doppler_index], sizeof(float) * effective_fft_size);
}
}
// Compute the test statistic
if (d_use_CFAR_algorithm_flag)
{
d_test_statistics = max_to_input_power_statistic(indext, doppler, d_input_power, d_num_doppler_bins, acq_parameters.doppler_max, d_doppler_step);
}
else
{
d_test_statistics = first_vs_second_peak_statistic(indext, doppler, d_num_doppler_bins, acq_parameters.doppler_max, d_doppler_step);
}
d_gnss_synchro->Acq_delay_samples = static_cast<double>(std::fmod(static_cast<float>(indext), acq_parameters.samples_per_code));
d_gnss_synchro->Acq_doppler_hz = static_cast<double>(doppler);
d_gnss_synchro->Acq_samplestamp_samples = samp_count;
}
else
{
for (unsigned int doppler_index = 0; doppler_index < acq_parameters.num_doppler_bins_step2; doppler_index++)
{
// doppler search steps
float doppler = d_doppler_center_step_two + (static_cast<float>(doppler_index) - static_cast<float>(acq_parameters.num_doppler_bins_step2) / 2.0) * acq_parameters.doppler_step2;
volk_32fc_x2_multiply_32fc(d_fft_if->get_inbuf(), in, d_grid_doppler_wipeoffs_step_two[doppler_index], d_fft_size);
// 3- Perform the FFT-based convolution (parallel time search)
@@ -543,55 +671,31 @@ void pcps_acquisition::acquisition_core(unsigned long int samp_count)
// compute the inverse FFT
d_ifft->execute();
// Search maximum
size_t offset = (acq_parameters.bit_transition_flag ? effective_fft_size : 0);
volk_32fc_magnitude_squared_32f(d_magnitude, d_ifft->get_outbuf() + offset, effective_fft_size);
volk_gnsssdr_32f_index_max_32u(&indext, d_magnitude, effective_fft_size);
magt = d_magnitude[indext];
if (acq_parameters.use_CFAR_algorithm_flag)
if (d_num_noncoherent_integrations_counter == 1)
{
// Normalize the maximum value to correct the scale factor introduced by FFTW
magt = d_magnitude[indext] / (fft_normalization_factor * fft_normalization_factor);
volk_32fc_magnitude_squared_32f(d_magnitude_grid[doppler_index], d_ifft->get_outbuf() + offset, effective_fft_size);
}
// 4- record the maximum peak and the associated synchronization parameters
if (d_mag < magt)
else
{
d_mag = magt;
if (!acq_parameters.use_CFAR_algorithm_flag)
{
// Search grid noise floor approximation for this doppler line
volk_32f_accumulator_s32f(&d_input_power, d_magnitude, effective_fft_size);
d_input_power = (d_input_power - d_mag) / (effective_fft_size - 1);
}
// In case that acq_parameters.bit_transition_flag = true, we compare the potentially
// new maximum test statistics (d_mag/d_input_power) with the value in
// d_test_statistics. When the second dwell is being processed, the value
// of d_mag/d_input_power could be lower than d_test_statistics (i.e,
// the maximum test statistics in the previous dwell is greater than
// current d_mag/d_input_power). Note that d_test_statistics is not
// restarted between consecutive dwells in multidwell operation.
if (d_test_statistics < (d_mag / d_input_power) or !acq_parameters.bit_transition_flag)
{
d_gnss_synchro->Acq_delay_samples = static_cast<double>(indext % acq_parameters.samples_per_code);
d_gnss_synchro->Acq_doppler_hz = static_cast<double>(doppler);
d_gnss_synchro->Acq_samplestamp_samples = samp_count;
// 5- Compute the test statistics and compare to the threshold
//d_test_statistics = 2 * d_fft_size * d_mag / d_input_power;
d_test_statistics = d_mag / d_input_power;
}
}
// Record results to file if required
if (acq_parameters.dump and d_channel == d_dump_channel)
{
memcpy(grid_.colptr(doppler_index), d_magnitude, sizeof(float) * effective_fft_size);
volk_32fc_magnitude_squared_32f(d_tmp_buffer, d_ifft->get_outbuf() + offset, effective_fft_size);
volk_32f_x2_add_32f(d_magnitude_grid[doppler_index], d_magnitude_grid[doppler_index], d_tmp_buffer, effective_fft_size);
}
}
// Compute the test statistic
if (d_use_CFAR_algorithm_flag)
{
d_test_statistics = max_to_input_power_statistic(indext, doppler, d_input_power, acq_parameters.num_doppler_bins_step2, static_cast<int>(d_doppler_center_step_two - (static_cast<float>(acq_parameters.num_doppler_bins_step2) / 2.0) * acq_parameters.doppler_step2), acq_parameters.doppler_step2);
}
else
{
d_test_statistics = first_vs_second_peak_statistic(indext, doppler, acq_parameters.num_doppler_bins_step2, static_cast<int>(d_doppler_center_step_two - (static_cast<float>(acq_parameters.num_doppler_bins_step2) / 2.0) * acq_parameters.doppler_step2), acq_parameters.doppler_step2);
}
d_gnss_synchro->Acq_delay_samples = static_cast<double>(std::fmod(static_cast<float>(indext), acq_parameters.samples_per_code));
d_gnss_synchro->Acq_doppler_hz = static_cast<double>(doppler);
d_gnss_synchro->Acq_samplestamp_samples = samp_count;
}
lk.lock();
if (!acq_parameters.bit_transition_flag)
{
@@ -618,12 +722,13 @@ void pcps_acquisition::acquisition_core(unsigned long int samp_count)
d_state = 0; // Positive acquisition
}
}
else if (d_well_count == acq_parameters.max_dwells)
if (d_num_noncoherent_integrations_counter == acq_parameters.max_dwells)
{
if (d_state != 0) send_negative_acquisition();
d_state = 0;
d_active = false;
d_step_two = false;
send_negative_acquisition();
}
}
else
@@ -659,10 +764,24 @@ void pcps_acquisition::acquisition_core(unsigned long int samp_count)
}
}
d_worker_active = false;
// Record results to file if required
if (acq_parameters.dump and d_channel == d_dump_channel)
if ((d_num_noncoherent_integrations_counter == acq_parameters.max_dwells) or (d_positive_acq == 1))
{
pcps_acquisition::dump_results(effective_fft_size);
// Record results to file if required
if (acq_parameters.dump and d_channel == d_dump_channel)
{
pcps_acquisition::dump_results(effective_fft_size);
}
d_num_noncoherent_integrations_counter = 0;
d_positive_acq = 0;
// Reset grid
for (unsigned int i = 0; i < d_num_doppler_bins; i++)
{
for (unsigned k = 0; k < d_fft_size; k++)
{
d_magnitude_grid[i][k] = 0.0;
}
}
}
}
@@ -681,13 +800,12 @@ int pcps_acquisition::general_work(int noutput_items __attribute__((unused)),
* 5. Compute the test statistics and compare to the threshold
* 6. Declare positive or negative acquisition using a message port
*/
gr::thread::scoped_lock lk(d_setlock);
if (!d_active or d_worker_active)
{
if (!acq_parameters.blocking_on_standby)
{
d_sample_counter += d_fft_size * ninput_items[0];
d_sample_counter += d_consumed_samples * ninput_items[0];
consume_each(ninput_items[0]);
}
if (d_step_two)
@@ -708,14 +826,13 @@ int pcps_acquisition::general_work(int noutput_items __attribute__((unused)),
d_gnss_synchro->Acq_delay_samples = 0.0;
d_gnss_synchro->Acq_doppler_hz = 0.0;
d_gnss_synchro->Acq_samplestamp_samples = 0;
d_well_count = 0;
d_mag = 0.0;
d_input_power = 0.0;
d_test_statistics = 0.0;
d_state = 1;
if (!acq_parameters.blocking_on_standby)
{
d_sample_counter += d_fft_size * ninput_items[0]; // sample counter
d_sample_counter += d_consumed_samples * ninput_items[0]; // sample counter
consume_each(ninput_items[0]);
}
break;
@@ -726,11 +843,11 @@ int pcps_acquisition::general_work(int noutput_items __attribute__((unused)),
// Copy the data to the core and let it know that new data is available
if (d_cshort)
{
memcpy(d_data_buffer_sc, input_items[0], d_fft_size * sizeof(lv_16sc_t));
memcpy(d_data_buffer_sc, input_items[0], d_consumed_samples * sizeof(lv_16sc_t));
}
else
{
memcpy(d_data_buffer, input_items[0], d_fft_size * sizeof(gr_complex));
memcpy(d_data_buffer, input_items[0], d_consumed_samples * sizeof(gr_complex));
}
if (acq_parameters.blocking)
{
@@ -742,7 +859,7 @@ int pcps_acquisition::general_work(int noutput_items __attribute__((unused)),
gr::thread::thread d_worker(&pcps_acquisition::acquisition_core, this, d_sample_counter);
d_worker_active = true;
}
d_sample_counter += d_fft_size;
d_sample_counter += d_consumed_samples;
consume_each(1);
break;
}

13
src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition.h
View File

@@ -95,24 +95,33 @@ private:
void dump_results(int effective_fft_size);
float first_vs_second_peak_statistic(uint32_t& indext, int& doppler, unsigned int num_doppler_bins, int doppler_max, int doppler_step);
float max_to_input_power_statistic(uint32_t& indext, int& doppler, float input_power, unsigned int num_doppler_bins, int doppler_max, int doppler_step);
Acq_Conf acq_parameters;
bool d_active;
bool d_worker_active;
bool d_cshort;
bool d_step_two;
bool d_use_CFAR_algorithm_flag;
int d_positive_acq;
float d_threshold;
float d_mag;
float d_input_power;
float d_test_statistics;
float* d_magnitude;
float** d_magnitude_grid;
float* d_tmp_buffer;
gr_complex* d_input_signal;
uint32_t d_samplesPerChip;
long d_old_freq;
int d_state;
unsigned int d_channel;
unsigned int d_doppler_step;
float d_doppler_center_step_two;
unsigned int d_well_count;
unsigned int d_num_noncoherent_integrations_counter;
unsigned int d_fft_size;
unsigned int d_consumed_samples;
unsigned int d_num_doppler_bins;
unsigned long int d_sample_counter;
gr_complex** d_grid_doppler_wipeoffs;
@@ -150,7 +159,7 @@ public:
}
/*!
* \brief Initializes acquisition algorithm.
* \brief Initializes acquisition algorithm and reserves memory.
*/
void init();

408
src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc
View File

@@ -40,46 +40,41 @@
#include <volk_gnsssdr/volk_gnsssdr.h>
#include <algorithm> // std::rotate, std::fill_n
#include <sstream>
#include <matio.h>
using google::LogMessage;
pcps_acquisition_fine_doppler_cc_sptr pcps_make_acquisition_fine_doppler_cc(
int max_dwells, unsigned int sampled_ms, int doppler_max, int doppler_min,
long fs_in, int samples_per_ms, bool dump,
std::string dump_filename)
pcps_acquisition_fine_doppler_cc_sptr pcps_make_acquisition_fine_doppler_cc(const Acq_Conf &conf_)
{
return pcps_acquisition_fine_doppler_cc_sptr(
new pcps_acquisition_fine_doppler_cc(max_dwells, sampled_ms, doppler_max, doppler_min,
fs_in, samples_per_ms, dump, dump_filename));
new pcps_acquisition_fine_doppler_cc(conf_));
}
pcps_acquisition_fine_doppler_cc::pcps_acquisition_fine_doppler_cc(
int max_dwells, unsigned int sampled_ms, int doppler_max, int doppler_min,
long fs_in, int samples_per_ms, bool dump,
std::string dump_filename) : gr::block("pcps_acquisition_fine_doppler_cc",
gr::io_signature::make(1, 1, sizeof(gr_complex)),
gr::io_signature::make(0, 0, sizeof(gr_complex)))
pcps_acquisition_fine_doppler_cc::pcps_acquisition_fine_doppler_cc(const Acq_Conf &conf_)
: gr::block("pcps_acquisition_fine_doppler_cc",
gr::io_signature::make(1, 1, sizeof(gr_complex)),
gr::io_signature::make(0, 0, sizeof(gr_complex)))
{
this->message_port_register_out(pmt::mp("events"));
acq_parameters = conf_;
d_sample_counter = 0; // SAMPLE COUNTER
d_active = false;
d_fs_in = fs_in;
d_samples_per_ms = samples_per_ms;
d_sampled_ms = sampled_ms;
d_config_doppler_max = doppler_max;
d_config_doppler_min = doppler_min;
d_fft_size = d_sampled_ms * d_samples_per_ms;
d_fs_in = conf_.fs_in;
d_samples_per_ms = conf_.samples_per_ms;
d_config_doppler_max = conf_.doppler_max;
d_fft_size = d_samples_per_ms;
// HS Acquisition
d_max_dwells = max_dwells;
d_max_dwells = conf_.max_dwells;
d_gnuradio_forecast_samples = d_fft_size;
d_input_power = 0.0;
d_state = 0;
d_carrier = static_cast<gr_complex *>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
d_fft_codes = static_cast<gr_complex *>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
d_magnitude = static_cast<float *>(volk_gnsssdr_malloc(d_fft_size * sizeof(float), volk_gnsssdr_get_alignment()));
d_10_ms_buffer = static_cast<gr_complex *>(volk_gnsssdr_malloc(50 * d_samples_per_ms * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
// Direct FFT
d_fft_if = new gr::fft::fft_complex(d_fft_size, true);
@@ -87,10 +82,10 @@ pcps_acquisition_fine_doppler_cc::pcps_acquisition_fine_doppler_cc(
d_ifft = new gr::fft::fft_complex(d_fft_size, false);
// For dumping samples into a file
d_dump = dump;
d_dump_filename = dump_filename;
d_dump = conf_.dump;
d_dump_filename = conf_.dump_filename;
d_doppler_resolution = 0;
d_n_samples_in_buffer = 0;
d_threshold = 0;
d_num_doppler_points = 0;
d_doppler_step = 0;
@@ -102,21 +97,46 @@ pcps_acquisition_fine_doppler_cc::pcps_acquisition_fine_doppler_cc(
d_test_statistics = 0;
d_well_count = 0;
d_channel = 0;
d_positive_acq = 0;
d_dump_number = 0;
d_dump_channel = 0; // this implementation can only produce dumps in channel 0
//todo: migrate config parameters to the unified acquisition config class
}
// Finds next power of two
// for n. If n itself is a
// power of two then returns n
unsigned int pcps_acquisition_fine_doppler_cc::nextPowerOf2(unsigned int n)
{
n--;
n |= n >> 1;
n |= n >> 2;
n |= n >> 4;
n |= n >> 8;
n |= n >> 16;
n++;
return n;
}
void pcps_acquisition_fine_doppler_cc::set_doppler_step(unsigned int doppler_step)
{
d_doppler_step = doppler_step;
// Create the search grid array
d_num_doppler_points = floor(std::abs(d_config_doppler_max - d_config_doppler_min) / d_doppler_step);
d_num_doppler_points = floor(std::abs(2 * d_config_doppler_max) / d_doppler_step);
d_grid_data = new float *[d_num_doppler_points];
for (int i = 0; i < d_num_doppler_points; i++)
{
d_grid_data[i] = static_cast<float *>(volk_gnsssdr_malloc(d_fft_size * sizeof(float), volk_gnsssdr_get_alignment()));
}
if (d_dump)
{
grid_ = arma::fmat(d_fft_size, d_num_doppler_points, arma::fill::zeros);
}
update_carrier_wipeoff();
}
@@ -138,12 +158,9 @@ pcps_acquisition_fine_doppler_cc::~pcps_acquisition_fine_doppler_cc()
volk_gnsssdr_free(d_carrier);
volk_gnsssdr_free(d_fft_codes);
volk_gnsssdr_free(d_magnitude);
volk_gnsssdr_free(d_10_ms_buffer);
delete d_ifft;
delete d_fft_if;
if (d_dump)
{
d_dump_file.close();
}
free_grid_memory();
}
@@ -167,7 +184,6 @@ void pcps_acquisition_fine_doppler_cc::init()
d_gnss_synchro->Acq_delay_samples = 0.0;
d_gnss_synchro->Acq_doppler_hz = 0.0;
d_gnss_synchro->Acq_samplestamp_samples = 0;
d_input_power = 0.0;
d_state = 0;
}
@@ -187,6 +203,7 @@ void pcps_acquisition_fine_doppler_cc::reset_grid()
d_well_count = 0;
for (int i = 0; i < d_num_doppler_points; i++)
{
//todo: use memset here
for (unsigned int j = 0; j < d_fft_size; j++)
{
d_grid_data[i][j] = 0.0;
@@ -203,7 +220,7 @@ void pcps_acquisition_fine_doppler_cc::update_carrier_wipeoff()
d_grid_doppler_wipeoffs = new gr_complex *[d_num_doppler_points];
for (int doppler_index = 0; doppler_index < d_num_doppler_points; doppler_index++)
{
doppler_hz = d_config_doppler_min + d_doppler_step * doppler_index;
doppler_hz = d_doppler_step * doppler_index - d_config_doppler_max;
// doppler search steps
// compute the carrier doppler wipe-off signal and store it
phase_step_rad = static_cast<float>(GPS_TWO_PI) * doppler_hz / static_cast<float>(d_fs_in);
@@ -215,52 +232,70 @@ void pcps_acquisition_fine_doppler_cc::update_carrier_wipeoff()
}
double pcps_acquisition_fine_doppler_cc::search_maximum()
double pcps_acquisition_fine_doppler_cc::compute_CAF()
{
float magt = 0.0;
float fft_normalization_factor;
float firstPeak = 0.0;
int index_doppler = 0;
uint32_t tmp_intex_t = 0;
uint32_t index_time = 0;
// Look for correlation peaks in the results ==============================
// Find the highest peak and compare it to the second highest peak
// The second peak is chosen not closer than 1 chip to the highest peak
//--- Find the correlation peak and the carrier frequency --------------
for (int i = 0; i < d_num_doppler_points; i++)
{
volk_gnsssdr_32f_index_max_32u(&tmp_intex_t, d_grid_data[i], d_fft_size);
if (d_grid_data[i][tmp_intex_t] > magt)
if (d_grid_data[i][tmp_intex_t] > firstPeak)
{
magt = d_grid_data[i][tmp_intex_t];
//std::cout<<magt<<std::endl;
firstPeak = d_grid_data[i][tmp_intex_t];
index_doppler = i;
index_time = tmp_intex_t;
}
// Record results to file if required
if (d_dump and d_channel == d_dump_channel)
{
memcpy(grid_.colptr(i), d_grid_data[i], sizeof(float) * d_fft_size);
}
}
// Normalize the maximum value to correct the scale factor introduced by FFTW
fft_normalization_factor = static_cast<float>(d_fft_size) * static_cast<float>(d_fft_size);
magt = magt / (fft_normalization_factor * fft_normalization_factor);
// -- - Find 1 chip wide code phase exclude range around the peak
uint32_t samplesPerChip = ceil(GPS_L1_CA_CHIP_PERIOD * static_cast<float>(this->d_fs_in));
int32_t excludeRangeIndex1 = index_time - samplesPerChip;
int32_t excludeRangeIndex2 = index_time + samplesPerChip;
// -- - Correct code phase exclude range if the range includes array boundaries
if (excludeRangeIndex1 < 0)
{
excludeRangeIndex1 = d_fft_size + excludeRangeIndex1;
}
else if (excludeRangeIndex2 >= static_cast<int>(d_fft_size))
{
excludeRangeIndex2 = excludeRangeIndex2 - d_fft_size;
}
int32_t idx = excludeRangeIndex1;
do
{
d_grid_data[index_doppler][idx] = 0.0;
idx++;
if (idx == static_cast<int>(d_fft_size)) idx = 0;
}
while (idx != excludeRangeIndex2);
//--- Find the second highest correlation peak in the same freq. bin ---
volk_gnsssdr_32f_index_max_32u(&tmp_intex_t, d_grid_data[index_doppler], d_fft_size);
float secondPeak = d_grid_data[index_doppler][tmp_intex_t];
// 5- Compute the test statistics and compare to the threshold
d_test_statistics = magt / (d_input_power * std::sqrt(d_well_count));
d_test_statistics = firstPeak / secondPeak;
// 4- record the maximum peak and the associated synchronization parameters
d_gnss_synchro->Acq_delay_samples = static_cast<double>(index_time);
d_gnss_synchro->Acq_doppler_hz = static_cast<double>(index_doppler * d_doppler_step + d_config_doppler_min);
d_gnss_synchro->Acq_doppler_hz = static_cast<double>(index_doppler * d_doppler_step - d_config_doppler_max);
d_gnss_synchro->Acq_samplestamp_samples = d_sample_counter;
// Record results to file if required
if (d_dump)
{
std::stringstream filename;
std::streamsize n = 2 * sizeof(float) * (d_fft_size); // complex file write
filename.str("");
filename << "../data/test_statistics_" << d_gnss_synchro->System
<< "_" << d_gnss_synchro->Signal << "_sat_"
<< d_gnss_synchro->PRN << "_doppler_" << d_gnss_synchro->Acq_doppler_hz << ".dat";
d_dump_file.open(filename.str().c_str(), std::ios::out | std::ios::binary);
d_dump_file.write(reinterpret_cast<char *>(d_grid_data[index_doppler]), n); //write directly |abs(x)|^2 in this Doppler bin?
d_dump_file.close();
}
return d_test_statistics;
}
@@ -296,6 +331,7 @@ int pcps_acquisition_fine_doppler_cc::compute_and_accumulate_grid(gr_vector_cons
// doppler search steps
// Perform the carrier wipe-off
volk_32fc_x2_multiply_32fc(d_fft_if->get_inbuf(), in, d_grid_doppler_wipeoffs[doppler_index], d_fft_size);
// 3- Perform the FFT-based convolution (parallel time search)
// Compute the FFT of the carrier wiped--off incoming signal
d_fft_if->execute();
@@ -308,29 +344,38 @@ int pcps_acquisition_fine_doppler_cc::compute_and_accumulate_grid(gr_vector_cons
d_ifft->execute();
// save the grid matrix delay file
volk_32fc_magnitude_squared_32f(p_tmp_vector, d_ifft->get_outbuf(), d_fft_size);
const float *old_vector = d_grid_data[doppler_index];
volk_32f_x2_add_32f(d_grid_data[doppler_index], old_vector, p_tmp_vector, d_fft_size);
//accumulate grid values
volk_32f_x2_add_32f(d_grid_data[doppler_index], d_grid_data[doppler_index], p_tmp_vector, d_fft_size);
}
volk_gnsssdr_free(p_tmp_vector);
return d_fft_size;
//debug
// std::cout << "iff=[";
// for (int n = 0; n < d_fft_size; n++)
// {
// std::cout << std::real(d_ifft->get_outbuf()[n]) << "+" << std::imag(d_ifft->get_outbuf()[n]) << "i,";
// }
// std::cout << "]\n";
// getchar();
}
int pcps_acquisition_fine_doppler_cc::estimate_Doppler(gr_vector_const_void_star &input_items)
int pcps_acquisition_fine_doppler_cc::estimate_Doppler()
{
// Direct FFT
int zero_padding_factor = 2;
int fft_size_extended = d_fft_size * zero_padding_factor;
int zero_padding_factor = 8;
int prn_replicas = 10;
int signal_samples = prn_replicas * d_fft_size;
//int fft_size_extended = nextPowerOf2(signal_samples * zero_padding_factor);
int fft_size_extended = signal_samples * zero_padding_factor;
gr::fft::fft_complex *fft_operator = new gr::fft::fft_complex(fft_size_extended, true);
//zero padding the entire vector
std::fill_n(fft_operator->get_inbuf(), fft_size_extended, gr_complex(0.0, 0.0));
//1. generate local code aligned with the acquisition code phase estimation
gr_complex *code_replica = static_cast<gr_complex *>(volk_gnsssdr_malloc(d_fft_size * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
gr_complex *code_replica = static_cast<gr_complex *>(volk_gnsssdr_malloc(signal_samples * sizeof(gr_complex), volk_gnsssdr_get_alignment()));
gps_l1_ca_code_gen_complex_sampled(code_replica, d_gnss_synchro->PRN, d_fs_in, 0);
@@ -342,10 +387,12 @@ int pcps_acquisition_fine_doppler_cc::estimate_Doppler(gr_vector_const_void_star
std::rotate(code_replica, code_replica + (d_fft_size - shift_index), code_replica + d_fft_size - 1);
}
for (int n = 0; n < prn_replicas - 1; n++)
{
memcpy(&code_replica[(n + 1) * d_fft_size], code_replica, d_fft_size * sizeof(gr_complex));
}
//2. Perform code wipe-off
const gr_complex *in = reinterpret_cast<const gr_complex *>(input_items[0]); //Get the input samples pointer
volk_32fc_x2_multiply_32fc(fft_operator->get_inbuf(), in, code_replica, d_fft_size);
volk_32fc_x2_multiply_32fc(fft_operator->get_inbuf(), d_10_ms_buffer, code_replica, signal_samples);
// 3. Perform the FFT (zero padded!)
fft_operator->execute();
@@ -360,8 +407,8 @@ int pcps_acquisition_fine_doppler_cc::estimate_Doppler(gr_vector_const_void_star
//case even
int counter = 0;
float *fftFreqBins = new float[fft_size_extended];
float fftFreqBins[fft_size_extended];
std::fill_n(fftFreqBins, fft_size_extended, 0.0);
for (int k = 0; k < (fft_size_extended / 2); k++)
@@ -372,7 +419,7 @@ int pcps_acquisition_fine_doppler_cc::estimate_Doppler(gr_vector_const_void_star
for (int k = fft_size_extended / 2; k > 0; k--)
{
fftFreqBins[counter] = ((-static_cast<float>(d_fs_in) / 2) * static_cast<float>(k)) / (static_cast<float>(fft_size_extended) / 2.0);
fftFreqBins[counter] = ((-static_cast<float>(d_fs_in) / 2.0) * static_cast<float>(k)) / (static_cast<float>(fft_size_extended) / 2.0);
counter++;
}
@@ -380,50 +427,56 @@ int pcps_acquisition_fine_doppler_cc::estimate_Doppler(gr_vector_const_void_star
if (std::abs(fftFreqBins[tmp_index_freq] - d_gnss_synchro->Acq_doppler_hz) < 1000)
{
d_gnss_synchro->Acq_doppler_hz = static_cast<double>(fftFreqBins[tmp_index_freq]);
//std::cout<<"FFT maximum present at "<<fftFreqBins[tmp_index_freq]<<" [Hz]"<<std::endl;
//std::cout << "FFT maximum present at " << fftFreqBins[tmp_index_freq] << " [Hz]" << std::endl;
}
else
{
DLOG(INFO) << "Abs(Grid Doppler - FFT Doppler)=" << std::abs(fftFreqBins[tmp_index_freq] - d_gnss_synchro->Acq_doppler_hz);
DLOG(INFO) << "Error estimating fine frequency Doppler";
//debug log
//
// std::cout<<"FFT maximum present at "<<fftFreqBins[tmp_index_freq]<<" [Hz]"<<std::endl;
// std::stringstream filename;
// std::streamsize n = sizeof(gr_complex) * (d_fft_size);
//
// filename.str("");
// filename << "../data/code_prn_" << d_gnss_synchro->PRN << ".dat";
// d_dump_file.open(filename.str().c_str(), std::ios::out
// | std::ios::binary);
// d_dump_file.write(reinterpret_cast<char*>(code_replica), n); //write directly |abs(x)|^2 in this Doppler bin?
// d_dump_file.close();
//
// filename.str("");
// filename << "../data/signal_prn_" << d_gnss_synchro->PRN << ".dat";
// d_dump_file.open(filename.str().c_str(), std::ios::out
// | std::ios::binary);
// d_dump_file.write(reinterpret_cast<char*>(in), n); //write directly |abs(x)|^2 in this Doppler bin?
// d_dump_file.close();
//
//
// n = sizeof(float) * (fft_size_extended);
// filename.str("");
// filename << "../data/fft_prn_" << d_gnss_synchro->PRN << ".dat";
// d_dump_file.open(filename.str().c_str(), std::ios::out
// | std::ios::binary);
// d_dump_file.write(reinterpret_cast<char*>(p_tmp_vector), n); //write directly |abs(x)|^2 in this Doppler bin?
// d_dump_file.close();
}
// free memory!!
delete fft_operator;
volk_gnsssdr_free(code_replica);
volk_gnsssdr_free(p_tmp_vector);
delete[] fftFreqBins;
return d_fft_size;
}
// Called by gnuradio to enable drivers, etc for i/o devices.
bool pcps_acquisition_fine_doppler_cc::start()
{
d_sample_counter = 0;
return true;
}
void pcps_acquisition_fine_doppler_cc::set_state(int state)
{
//gr::thread::scoped_lock lock(d_setlock); // require mutex with work function called by the scheduler
d_state = state;
if (d_state == 1)
{
d_gnss_synchro->Acq_delay_samples = 0.0;
d_gnss_synchro->Acq_doppler_hz = 0.0;
d_gnss_synchro->Acq_samplestamp_samples = 0;
d_well_count = 0;
d_test_statistics = 0.0;
d_active = true;
reset_grid();
}
else if (d_state == 0)
{
}
else
{
LOG(ERROR) << "State can only be set to 0 or 1";
}
}
int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
gr_vector_int &ninput_items __attribute__((unused)), gr_vector_const_void_star &input_items,
gr_vector_void_star &output_items __attribute__((unused)))
@@ -442,29 +495,36 @@ int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
* S5. Negative_Acq: Send message and stop acq -> S0
*/
int samples_remaining;
switch (d_state)
{
case 0: // S0. StandBy
//DLOG(INFO) <<"S0"<<std::endl;
if (d_active == true)
{
reset_grid();
d_n_samples_in_buffer = 0;
d_state = 1;
}
if (!acq_parameters.blocking_on_standby)
{
d_sample_counter += d_fft_size; // sample counter
consume_each(d_fft_size);
}
break;
case 1: // S1. ComputeGrid
//DLOG(INFO) <<"S1"<<std::endl;
compute_and_accumulate_grid(input_items);
memcpy(&d_10_ms_buffer[d_n_samples_in_buffer], reinterpret_cast<const gr_complex *>(input_items[0]), d_fft_size * sizeof(gr_complex));
d_n_samples_in_buffer += d_fft_size;
d_well_count++;
if (d_well_count >= d_max_dwells)
{
d_state = 2;
}
d_sample_counter += d_fft_size; // sample counter
consume_each(d_fft_size);
break;
case 2: // Compute test statistics and decide
//DLOG(INFO) <<"S2"<<std::endl;
d_input_power = estimate_input_power(input_items);
d_test_statistics = search_maximum();
d_test_statistics = compute_CAF();
if (d_test_statistics > d_threshold)
{
d_state = 3; //perform fine doppler estimation
@@ -472,16 +532,34 @@ int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
else
{
d_state = 5; //negative acquisition
d_n_samples_in_buffer = 0;
}
break;
case 3: // Fine doppler estimation
//DLOG(INFO) <<"S3"<<std::endl;
DLOG(INFO) << "Performing fine Doppler estimation";
estimate_Doppler(input_items); //disabled in repo
d_state = 4;
samples_remaining = 10 * d_samples_per_ms - d_n_samples_in_buffer;
if (samples_remaining > noutput_items)
{
memcpy(&d_10_ms_buffer[d_n_samples_in_buffer], reinterpret_cast<const gr_complex *>(input_items[0]), noutput_items * sizeof(gr_complex));
d_n_samples_in_buffer += noutput_items;
d_sample_counter += noutput_items; // sample counter
consume_each(noutput_items);
}
else
{
if (samples_remaining > 0)
{
memcpy(&d_10_ms_buffer[d_n_samples_in_buffer], reinterpret_cast<const gr_complex *>(input_items[0]), samples_remaining * sizeof(gr_complex));
d_sample_counter += samples_remaining; // sample counter
consume_each(samples_remaining);
}
estimate_Doppler(); //disabled in repo
d_n_samples_in_buffer = 0;
d_state = 4;
}
break;
case 4: // Positive_Acq
//DLOG(INFO) <<"S4"<<std::endl;
DLOG(INFO) << "positive acquisition";
DLOG(INFO) << "satellite " << d_gnss_synchro->System << " " << d_gnss_synchro->PRN;
DLOG(INFO) << "sample_stamp " << d_sample_counter;
@@ -489,15 +567,23 @@ int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
DLOG(INFO) << "test statistics threshold " << d_threshold;
DLOG(INFO) << "code phase " << d_gnss_synchro->Acq_delay_samples;
DLOG(INFO) << "doppler " << d_gnss_synchro->Acq_doppler_hz;
DLOG(INFO) << "input signal power " << d_input_power;
d_positive_acq = 1;
d_active = false;
// Record results to file if required
if (d_dump and d_channel == d_dump_channel)
{
dump_results(d_fft_size);
}
// Send message to channel port //0=STOP_CHANNEL 1=ACQ_SUCCEES 2=ACQ_FAIL
this->message_port_pub(pmt::mp("events"), pmt::from_long(1));
d_state = 0;
if (!acq_parameters.blocking_on_standby)
{
d_sample_counter += noutput_items; // sample counter
consume_each(noutput_items);
}
break;
case 5: // Negative_Acq
//DLOG(INFO) <<"S5"<<std::endl;
DLOG(INFO) << "negative acquisition";
DLOG(INFO) << "satellite " << d_gnss_synchro->System << " " << d_gnss_synchro->PRN;
DLOG(INFO) << "sample_stamp " << d_sample_counter;
@@ -505,20 +591,108 @@ int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
DLOG(INFO) << "test statistics threshold " << d_threshold;
DLOG(INFO) << "code phase " << d_gnss_synchro->Acq_delay_samples;
DLOG(INFO) << "doppler " << d_gnss_synchro->Acq_doppler_hz;
DLOG(INFO) << "input signal power " << d_input_power;
d_positive_acq = 0;
d_active = false;
// Record results to file if required
if (d_dump and d_channel == d_dump_channel)
{
dump_results(d_fft_size);
}
// Send message to channel port //0=STOP_CHANNEL 1=ACQ_SUCCEES 2=ACQ_FAIL
this->message_port_pub(pmt::mp("events"), pmt::from_long(2));
d_state = 0;
if (!acq_parameters.blocking_on_standby)
{
d_sample_counter += noutput_items; // sample counter
consume_each(noutput_items);
}
break;
default:
d_state = 0;
if (!acq_parameters.blocking_on_standby)
{
d_sample_counter += noutput_items; // sample counter
consume_each(noutput_items);
}
break;
}
//DLOG(INFO)<<"d_sample_counter="<<d_sample_counter<<std::endl;
d_sample_counter += d_fft_size; // sample counter
consume_each(d_fft_size);
return noutput_items;
return 0;
}
void pcps_acquisition_fine_doppler_cc::dump_results(int effective_fft_size)
{
d_dump_number++;
std::string filename = d_dump_filename;
filename.append("_");
filename.append(1, d_gnss_synchro->System);
filename.append("_");
filename.append(1, d_gnss_synchro->Signal[0]);
filename.append(1, d_gnss_synchro->Signal[1]);
filename.append("_ch_");
filename.append(std::to_string(d_channel));
filename.append("_");
filename.append(std::to_string(d_dump_number));
filename.append("_sat_");
filename.append(std::to_string(d_gnss_synchro->PRN));
filename.append(".mat");
mat_t *matfp = Mat_CreateVer(filename.c_str(), NULL, MAT_FT_MAT73);
if (matfp == NULL)
{
std::cout << "Unable to create or open Acquisition dump file" << std::endl;
d_dump = false;
}
else
{
size_t dims[2] = {static_cast<size_t>(effective_fft_size), static_cast<size_t>(d_num_doppler_points)};
matvar_t *matvar = Mat_VarCreate("acq_grid", MAT_C_SINGLE, MAT_T_SINGLE, 2, dims, grid_.memptr(), 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
dims[0] = static_cast<size_t>(1);
dims[1] = static_cast<size_t>(1);
matvar = Mat_VarCreate("doppler_max", MAT_C_UINT32, MAT_T_UINT32, 1, dims, &d_config_doppler_max, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
matvar = Mat_VarCreate("doppler_step", MAT_C_UINT32, MAT_T_UINT32, 1, dims, &d_doppler_step, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
matvar = Mat_VarCreate("d_positive_acq", MAT_C_INT32, MAT_T_INT32, 1, dims, &d_positive_acq, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
float aux = static_cast<float>(d_gnss_synchro->Acq_doppler_hz);
matvar = Mat_VarCreate("acq_doppler_hz", MAT_C_SINGLE, MAT_T_SINGLE, 1, dims, &aux, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
aux = static_cast<float>(d_gnss_synchro->Acq_delay_samples);
matvar = Mat_VarCreate("acq_delay_samples", MAT_C_SINGLE, MAT_T_SINGLE, 1, dims, &aux, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
matvar = Mat_VarCreate("test_statistic", MAT_C_SINGLE, MAT_T_SINGLE, 1, dims, &d_test_statistics, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
matvar = Mat_VarCreate("threshold", MAT_C_SINGLE, MAT_T_SINGLE, 1, dims, &d_threshold, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
aux = 0.0;
matvar = Mat_VarCreate("input_power", MAT_C_SINGLE, MAT_T_SINGLE, 1, dims, &aux, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
matvar = Mat_VarCreate("sample_counter", MAT_C_UINT64, MAT_T_UINT64, 1, dims, &d_sample_counter, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
matvar = Mat_VarCreate("PRN", MAT_C_UINT32, MAT_T_UINT32, 1, dims, &d_gnss_synchro->PRN, 0);
Mat_VarWrite(matfp, matvar, MAT_COMPRESSION_ZLIB); // or MAT_COMPRESSION_NONE
Mat_VarFree(matvar);
Mat_Close(matfp);
}
}

61
src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.h
View File

@@ -1,8 +1,9 @@
/*!
* \file pcps_acquisition_fine_doppler_acquisition_cc.h
* \brief This class implements a Parallel Code Phase Search Acquisition with multi-dwells and fine Doppler estimation
* for GPS L1 C/A signal
*
* Acquisition strategy (Kay Borre book + CFAR threshold).
* Acquisition strategy (Kay Borre book).
* <ol>
* <li> Compute the input signal power estimation
* <li> Doppler serial search loop
@@ -49,6 +50,8 @@
#define GNSS_SDR_PCPS_ACQUISITION_FINE_DOPPLER_CC_H_
#include "gnss_synchro.h"
#include "acq_conf.h"
#include <armadillo>
#include <gnuradio/block.h>
#include <gnuradio/gr_complex.h>
#include <gnuradio/fft/fft.h>
@@ -61,59 +64,44 @@ typedef boost::shared_ptr<pcps_acquisition_fine_doppler_cc>
pcps_acquisition_fine_doppler_cc_sptr;
pcps_acquisition_fine_doppler_cc_sptr
pcps_make_acquisition_fine_doppler_cc(int max_dwells, unsigned int sampled_ms,
int doppler_max, int doppler_min, long fs_in, int samples_per_ms,
bool dump, std::string dump_filename);
pcps_make_acquisition_fine_doppler_cc(const Acq_Conf& conf_);
/*!
* \brief This class implements a Parallel Code Phase Search Acquisition.
*
* Check \ref Navitec2012 "An Open Source Galileo E1 Software Receiver",
* Algorithm 1, for a pseudocode description of this implementation.
*/
class pcps_acquisition_fine_doppler_cc : public gr::block
{
private:
friend pcps_acquisition_fine_doppler_cc_sptr
pcps_make_acquisition_fine_doppler_cc(int max_dwells, unsigned int sampled_ms,
int doppler_max, int doppler_min, long fs_in,
int samples_per_ms, bool dump,
std::string dump_filename);
pcps_acquisition_fine_doppler_cc(int max_dwells, unsigned int sampled_ms,
int doppler_max, int doppler_min, long fs_in,
int samples_per_ms, bool dump,
std::string dump_filename);
void calculate_magnitudes(gr_complex* fft_begin, int doppler_shift,
int doppler_offset);
pcps_make_acquisition_fine_doppler_cc(const Acq_Conf& conf_);
pcps_acquisition_fine_doppler_cc(const Acq_Conf& conf_);
int compute_and_accumulate_grid(gr_vector_const_void_star& input_items);
int estimate_Doppler(gr_vector_const_void_star& input_items);
int estimate_Doppler();
float estimate_input_power(gr_vector_const_void_star& input_items);
double search_maximum();
double compute_CAF();
void reset_grid();
void update_carrier_wipeoff();
void free_grid_memory();
bool start();
Acq_Conf acq_parameters;
long d_fs_in;
int d_samples_per_ms;
int d_max_dwells;
unsigned int d_doppler_resolution;
int d_gnuradio_forecast_samples;
float d_threshold;
std::string d_satellite_str;
int d_config_doppler_max;
int d_config_doppler_min;
int d_num_doppler_points;
int d_doppler_step;
unsigned int d_sampled_ms;
unsigned int d_fft_size;
unsigned long int d_sample_counter;
gr_complex* d_carrier;
gr_complex* d_fft_codes;
gr_complex* d_10_ms_buffer;
float* d_magnitude;
float** d_grid_data;
@@ -124,17 +112,22 @@ private:
Gnss_Synchro* d_gnss_synchro;
unsigned int d_code_phase;
float d_doppler_freq;
float d_input_power;
float d_test_statistics;
std::ofstream d_dump_file;
int d_positive_acq;
int d_state;
bool d_active;
int d_well_count;
int d_n_samples_in_buffer;
bool d_dump;
unsigned int d_channel;
std::string d_dump_filename;
arma::fmat grid_;
long int d_dump_number;
unsigned int d_dump_channel;
public:
/*!
* \brief Default destructor.
@@ -187,6 +180,7 @@ public:
inline void set_channel(unsigned int channel)
{
d_channel = channel;
d_dump_channel = d_channel;
}
/*!
@@ -214,6 +208,13 @@ public:
*/
void set_doppler_step(unsigned int doppler_step);
/*!
* \brief If set to 1, ensures that acquisition starts at the
* first available sample.
* \param state - int=1 forces start of acquisition
*/
void set_state(int state);
/*!
* \brief Parallel Code Phase Search Acquisition signal processing.
*/
@@ -222,6 +223,14 @@ public:
gr_vector_void_star& output_items);
void forecast(int noutput_items, gr_vector_int& ninput_items_required);
/*!
* \brief Obtains the next power of 2 greater or equal to the input parameter
* \param n - Integer value to obtain the next power of 2.
*/
unsigned int nextPowerOf2(unsigned int n);
void dump_results(int effective_fft_size);
};
#endif /* pcps_acquisition_fine_doppler_cc*/