diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 8634c4a38..4502d321d 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -19,6 +19,12 @@ All notable changes to GNSS-SDR will be documented in this file. - Fix setting of the signal source gain if the AGC is enabled when using the AD9361 front-end. +### Improvements in Maintainability: + +- Rewritten Viterbi decoder for Galileo navigation messages: encapsulated in a + class instead of being implemented as free inline functions. This improves + memory management. + ### Improvements in Usability: - Added a new monitor to extract the decoded data bits of the navigation diff --git a/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.cc b/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.cc index d2a92bea6..d93d690ec 100644 --- a/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.cc +++ b/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.cc @@ -3,6 +3,7 @@ * \brief Implementation of a Galileo unified INAV and FNAV message demodulator * block * \author Javier Arribas 2018. jarribas(at)cttc.es + * \author Carles Fernandez, 2021. cfernandez(at)cttc.es * * * ----------------------------------------------------------------------------- @@ -18,32 +19,30 @@ #include "galileo_telemetry_decoder_gs.h" -#include "Galileo_E1.h" // for GALILEO_E1_CODE_PERIOD_MS -#include "Galileo_E5a.h" // for GALILEO_E5A_CODE_PERIOD_MS -#include "Galileo_E5b.h" // for GALILEO_E5B_CODE_PERIOD_MS -#include "Galileo_E6.h" // for GALILEO_E6_CODE_PERIOD_MS -#include "convolutional.h" -#include "display.h" +#include "Galileo_E1.h" // for GALILEO_E1_CODE_PERIOD_MS +#include "Galileo_E5a.h" // for GALILEO_E5A_CODE_PERIOD_MS +#include "Galileo_E5b.h" // for GALILEO_E5B_CODE_PERIOD_MS +#include "Galileo_E6.h" // for GALILEO_E6_CODE_PERIOD_MS +#include "display.h" // for colours in terminal: TEXT_BLUE, TEXT_RESET, ... #include "galileo_almanac_helper.h" // for Galileo_Almanac_Helper #include "galileo_ephemeris.h" // for Galileo_Ephemeris #include "galileo_has_page.h" // For Galileo_HAS_page #include "galileo_iono.h" // for Galileo_Iono #include "galileo_utc_model.h" // for Galileo_Utc_Model #include "gnss_sdr_make_unique.h" // for std::make_unique in C++11 -#include "gnss_synchro.h" -#include "tlm_crc_stats.h" -#include "tlm_utils.h" -#include -#include -#include // for make_any -#include // for mp -#include -#include // for fmod -#include // for size_t -#include // for abs -#include // for exception -#include // for cout -#include // for make_shared +#include "gnss_synchro.h" // for Gnss_Synchro +#include "tlm_crc_stats.h" // for Tlm_CRC_Stats +#include "tlm_utils.h" // for save_tlm_matfile, tlm_remove_file +#include "viterbi_decoder.h" // for Viterbi_Decoder +#include // for LOG, DLOG +#include // for gr::io_signature::make +#include // for pmt::make_any +#include // for pmt::mp +#include // for std::array +#include // for std::fmod, std::abs +#include // for size_t +#include // for std::exception +#include // for std::cout #define CRC_ERROR_LIMIT 6 @@ -86,9 +85,13 @@ galileo_telemetry_decoder_gs::galileo_telemetry_decoder_gs( d_remove_dat = conf.remove_dat; d_satellite = Gnss_Satellite(satellite.get_system(), satellite.get_PRN()); d_frame_type = frame_type; - d_nn = 2; - d_KK = 7; - d_mm = d_KK - 1; + + // Viterbi decoder vars + const int32_t nn = 2; // Coding rate 1/n + const int32_t KK = 7; // Constraint Length + const std::array g_encoder{{121, 91}}; // Polynomial G1 and G2 + d_mm = KK - 1; + DLOG(INFO) << "Initializing GALILEO UNIFIED TELEMETRY DECODER"; d_dump_crc_stats = conf.dump_crc_stats; @@ -116,7 +119,7 @@ galileo_telemetry_decoder_gs::galileo_telemetry_decoder_gs( d_preamble_samples = std::vector(d_samples_per_preamble); d_frame_length_symbols = GALILEO_INAV_PAGE_PART_SYMBOLS - GALILEO_INAV_PREAMBLE_LENGTH_BITS; d_codelength = static_cast(d_frame_length_symbols); - d_datalength = (d_codelength / d_nn) - d_mm; + d_datalength = (d_codelength / nn) - d_mm; d_max_symbols_without_valid_frame = GALILEO_INAV_PAGE_SYMBOLS * 30; // rise alarm 60 seconds without valid tlm if (conf.enable_reed_solomon == true) { @@ -136,7 +139,7 @@ galileo_telemetry_decoder_gs::galileo_telemetry_decoder_gs( d_preamble_samples = std::vector(d_samples_per_preamble); d_frame_length_symbols = GALILEO_FNAV_SYMBOLS_PER_PAGE - GALILEO_FNAV_PREAMBLE_LENGTH_BITS; d_codelength = static_cast(d_frame_length_symbols); - d_datalength = (d_codelength / d_nn) - d_mm; + d_datalength = (d_codelength / nn) - d_mm; d_max_symbols_without_valid_frame = GALILEO_FNAV_SYMBOLS_PER_PAGE * 5; // rise alarm 100 seconds without valid tlm break; } @@ -150,7 +153,7 @@ galileo_telemetry_decoder_gs::galileo_telemetry_decoder_gs( d_preamble_samples = std::vector(d_samples_per_preamble); d_frame_length_symbols = GALILEO_CNAV_SYMBOLS_PER_PAGE - GALILEO_CNAV_PREAMBLE_LENGTH_BITS; d_codelength = static_cast(d_frame_length_symbols); - d_datalength = (d_codelength / d_nn) - d_mm; + d_datalength = (d_codelength / nn) - d_mm; d_max_symbols_without_valid_frame = GALILEO_CNAV_SYMBOLS_PER_PAGE * 60; break; } @@ -230,21 +233,20 @@ galileo_telemetry_decoder_gs::galileo_telemetry_decoder_gs( d_cnav_dummy_page = false; d_print_cnav_page = true; + d_inav_nav.init_PRN(d_satellite.get_PRN()); + d_first_eph_sent = false; + // vars for Viterbi decoder const int32_t max_states = 1U << static_cast(d_mm); // 2^d_mm - std::array g_encoder{{121, 91}}; // Polynomial G1 and G2 - d_out0 = std::vector(max_states); d_out1 = std::vector(max_states); d_state0 = std::vector(max_states); d_state1 = std::vector(max_states); - d_inav_nav.init_PRN(d_satellite.get_PRN()); - d_first_eph_sent = false; - - // create appropriate transition matrices - nsc_transit(d_out0.data(), d_state0.data(), 0, g_encoder.data(), d_KK, d_nn); - nsc_transit(d_out1.data(), d_state1.data(), 1, g_encoder.data(), d_KK, d_nn); + // create Vitrebi decoder and appropriate transition vectors + d_viterbi = std::make_unique(KK, nn, d_datalength, g_encoder); + d_viterbi->nsc_transit(d_out0, d_state0, 0); + d_viterbi->nsc_transit(d_out1, d_state1, 1); } @@ -285,13 +287,6 @@ galileo_telemetry_decoder_gs::~galileo_telemetry_decoder_gs() } -void galileo_telemetry_decoder_gs::viterbi_decoder(float *page_part_symbols, int32_t *page_part_bits) -{ - Viterbi(page_part_bits, d_out0.data(), d_state0.data(), d_out1.data(), d_state1.data(), - page_part_symbols, d_KK, d_nn, d_datalength); -} - - void galileo_telemetry_decoder_gs::deinterleaver(int32_t rows, int32_t cols, const float *in, float *out) { for (int32_t r = 0; r < rows; r++) @@ -307,8 +302,8 @@ void galileo_telemetry_decoder_gs::deinterleaver(int32_t rows, int32_t cols, con void galileo_telemetry_decoder_gs::decode_INAV_word(float *page_part_symbols, int32_t frame_length) { // 1. De-interleave - std::vector page_part_symbols_deint(frame_length); - deinterleaver(GALILEO_INAV_INTERLEAVER_ROWS, GALILEO_INAV_INTERLEAVER_COLS, page_part_symbols, page_part_symbols_deint.data()); + std::vector page_part_symbols_soft_value(frame_length); + deinterleaver(GALILEO_INAV_INTERLEAVER_ROWS, GALILEO_INAV_INTERLEAVER_COLS, page_part_symbols, page_part_symbols_soft_value.data()); // 2. Viterbi decoder // 2.1 Take into account the NOT gate in G2 polynomial (Galileo ICD Figure 13, FEC encoder) @@ -317,12 +312,12 @@ void galileo_telemetry_decoder_gs::decode_INAV_word(float *page_part_symbols, in { if ((i + 1) % 2 == 0) { - page_part_symbols_deint[i] = -page_part_symbols_deint[i]; + page_part_symbols_soft_value[i] = -page_part_symbols_soft_value[i]; } } const int32_t decoded_length = frame_length / 2; std::vector page_part_bits(decoded_length); - viterbi_decoder(page_part_symbols_deint.data(), page_part_bits.data()); + d_viterbi->decode(page_part_bits, d_out0, d_state0, d_out1, d_state1, page_part_symbols_soft_value); // 3. Call the Galileo page decoder std::string page_String; @@ -459,8 +454,8 @@ void galileo_telemetry_decoder_gs::decode_INAV_word(float *page_part_symbols, in void galileo_telemetry_decoder_gs::decode_FNAV_word(float *page_symbols, int32_t frame_length) { // 1. De-interleave - std::vector page_symbols_deint(frame_length); - deinterleaver(GALILEO_FNAV_INTERLEAVER_ROWS, GALILEO_FNAV_INTERLEAVER_COLS, page_symbols, page_symbols_deint.data()); + std::vector page_symbols_soft_value(frame_length); + deinterleaver(GALILEO_FNAV_INTERLEAVER_ROWS, GALILEO_FNAV_INTERLEAVER_COLS, page_symbols, page_symbols_soft_value.data()); // 2. Viterbi decoder // 2.1 Take into account the NOT gate in G2 polynomial (Galileo ICD Figure 13, FEC encoder) @@ -469,13 +464,13 @@ void galileo_telemetry_decoder_gs::decode_FNAV_word(float *page_symbols, int32_t { if ((i + 1) % 2 == 0) { - page_symbols_deint[i] = -page_symbols_deint[i]; + page_symbols_soft_value[i] = -page_symbols_soft_value[i]; } } const int32_t decoded_length = frame_length / 2; std::vector page_bits(decoded_length); - viterbi_decoder(page_symbols_deint.data(), page_bits.data()); + d_viterbi->decode(page_bits, d_out0, d_state0, d_out1, d_state1, page_symbols_soft_value); // 3. Call the Galileo page decoder std::string page_String; @@ -533,8 +528,8 @@ void galileo_telemetry_decoder_gs::decode_FNAV_word(float *page_symbols, int32_t void galileo_telemetry_decoder_gs::decode_CNAV_word(float *page_symbols, int32_t page_length) { // 1. De-interleave - std::vector page_symbols_deint(page_length); - deinterleaver(GALILEO_CNAV_INTERLEAVER_ROWS, GALILEO_CNAV_INTERLEAVER_COLS, page_symbols, page_symbols_deint.data()); + std::vector page_symbols_soft_value(page_length); + deinterleaver(GALILEO_CNAV_INTERLEAVER_ROWS, GALILEO_CNAV_INTERLEAVER_COLS, page_symbols, page_symbols_soft_value.data()); // 2. Viterbi decoder // 2.1 Take into account the NOT gate in G2 polynomial (Galileo ICD Figure 13, FEC encoder) @@ -543,12 +538,12 @@ void galileo_telemetry_decoder_gs::decode_CNAV_word(float *page_symbols, int32_t { if ((i + 1) % 2 == 0) { - page_symbols_deint[i] = -page_symbols_deint[i]; + page_symbols_soft_value[i] = -page_symbols_soft_value[i]; } } const int32_t decoded_length = page_length / 2; std::vector page_bits(decoded_length); - viterbi_decoder(page_symbols_deint.data(), page_bits.data()); + d_viterbi->decode(page_bits, d_out0, d_state0, d_out1, d_state1, page_symbols_soft_value); // 3. Call the Galileo page decoder std::string page_String; @@ -622,6 +617,13 @@ void galileo_telemetry_decoder_gs::reset() d_last_valid_preamble = d_sample_counter; d_sent_tlm_failed_msg = false; d_stat = 0; + const int32_t max_states = 1U << static_cast(d_mm); // 2^d_mm + d_out0 = std::vector(max_states); + d_out1 = std::vector(max_states); + d_state0 = std::vector(max_states); + d_state1 = std::vector(max_states); + d_viterbi->nsc_transit(d_out0, d_state0, 0); + d_viterbi->nsc_transit(d_out1, d_state1, 1); DLOG(INFO) << "Telemetry decoder reset for satellite " << d_satellite; } @@ -732,7 +734,7 @@ int galileo_telemetry_decoder_gs::general_work(int noutput_items __attribute__(( corr_value += d_preamble_samples[i]; } } - if (abs(corr_value) >= d_samples_per_preamble) + if (std::abs(corr_value) >= d_samples_per_preamble) { d_preamble_index = d_sample_counter; // record the preamble sample stamp LOG(INFO) << "Preamble detection for Galileo satellite " << this->d_satellite << " in channel " << this->d_channel; @@ -760,11 +762,11 @@ int galileo_telemetry_decoder_gs::general_work(int noutput_items __attribute__(( corr_value += d_preamble_samples[i]; } } - if (abs(corr_value) >= d_samples_per_preamble) + if (std::abs(corr_value) >= d_samples_per_preamble) { // check preamble separation const auto preamble_diff = static_cast(d_sample_counter - d_preamble_index); - if (abs(preamble_diff - d_preamble_period_symbols) == 0) + if (std::abs(preamble_diff - d_preamble_period_symbols) == 0) { // try to decode frame DLOG(INFO) << "Starting page decoder for Galileo satellite " << this->d_satellite; diff --git a/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.h b/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.h index 729c11126..da09f1d29 100644 --- a/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.h +++ b/src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.h @@ -3,6 +3,7 @@ * \brief Implementation of a Galileo unified INAV and FNAV message demodulator * block * \author Javier Arribas 2018. jarribas(at)cttc.es + * \author Carles Fernandez, 2021. cfernandez(at)cttc.es * * * ----------------------------------------------------------------------------- @@ -20,30 +21,30 @@ #ifndef GNSS_SDR_GALILEO_TELEMETRY_DECODER_GS_H #define GNSS_SDR_GALILEO_TELEMETRY_DECODER_GS_H -#include "galileo_cnav_message.h" -#include "galileo_fnav_message.h" -#include "galileo_inav_message.h" -#include "gnss_block_interface.h" -#include "gnss_satellite.h" -#include "nav_message_packet.h" -#include "tlm_conf.h" -#include "tlm_crc_stats.h" -#include -#include // for block -#include // for gr_vector_const_void_star -#include -#include -#include // for std::unique_ptr -#include -#include +#include "galileo_cnav_message.h" // for Galileo_Cnav_Message +#include "galileo_fnav_message.h" // for Galileo_Fnav_Message +#include "galileo_inav_message.h" // for Galileo_Inav_Message +#include "gnss_block_interface.h" // for gnss_shared_ptr (adapts smart pointer type to GNU Radio version) +#include "gnss_satellite.h" // for Gnss_Satellite +#include "nav_message_packet.h" // for Nav_Message_Packet +#include "tlm_conf.h" // for Tlm_Conf +#include // for boost::circular_buffer +#include // for block +#include // for gr_vector_const_void_star +#include // for int32_t, uint32_t +#include // for std::ofstream +#include // for std::unique_ptr +#include // for std::string +#include // for std::vector /** \addtogroup Telemetry_Decoder * \{ */ /** \addtogroup Telemetry_Decoder_gnuradio_blocks * \{ */ - -class galileo_telemetry_decoder_gs; +class Viterbi_Decoder; // forward declaration +class Tlm_CRC_Stats; // forward declaration +class galileo_telemetry_decoder_gs; // forward declaration using galileo_telemetry_decoder_gs_sptr = gnss_shared_ptr; @@ -79,16 +80,13 @@ private: galileo_telemetry_decoder_gs(const Gnss_Satellite &satellite, const Tlm_Conf &conf, int frame_type); - int32_t d_nn; // Coding rate 1/n - int32_t d_KK; // Constraint Length - - void viterbi_decoder(float *page_part_symbols, int32_t *page_part_bits); void deinterleaver(int32_t rows, int32_t cols, const float *in, float *out); void decode_INAV_word(float *page_part_symbols, int32_t frame_length); void decode_FNAV_word(float *page_symbols, int32_t frame_length); void decode_CNAV_word(float *page_symbols, int32_t page_length); - // vars for Viterbi decoder + // data members for Viterbi decoder + std::unique_ptr d_viterbi; std::vector d_preamble_samples; std::vector d_page_part_symbols; std::vector d_out0; diff --git a/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.cc b/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.cc index fdde049eb..d803eee0d 100644 --- a/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.cc +++ b/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.cc @@ -16,7 +16,7 @@ #include "sbas_l1_telemetry_decoder_gs.h" #include "gnss_synchro.h" -#include "viterbi_decoder.h" +#include "viterbi_decoder_sbas.h" #include #include #include // for mp @@ -179,8 +179,8 @@ sbas_l1_telemetry_decoder_gs::Symbol_Aligner_And_Decoder::Symbol_Aligner_And_Dec const int32_t nn = 2; std::array g_encoder{121, 91}; - d_vd1 = std::make_shared(g_encoder.data(), d_KK, nn); - d_vd2 = std::make_shared(g_encoder.data(), d_KK, nn); + d_vd1 = std::make_shared(g_encoder.data(), d_KK, nn); + d_vd2 = std::make_shared(g_encoder.data(), d_KK, nn); d_past_symbol = 0; } diff --git a/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.h b/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.h index 3e12b9eeb..3fdea55fe 100644 --- a/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.h +++ b/src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.h @@ -37,7 +37,7 @@ * \{ */ -class Viterbi_Decoder; +class Viterbi_Decoder_Sbas; class sbas_l1_telemetry_decoder_gs; @@ -128,8 +128,8 @@ private: private: int32_t d_KK; - std::shared_ptr d_vd1; - std::shared_ptr d_vd2; + std::shared_ptr d_vd1; + std::shared_ptr d_vd2; double d_past_symbol; } d_symbol_aligner_and_decoder; diff --git a/src/algorithms/telemetry_decoder/libs/CMakeLists.txt b/src/algorithms/telemetry_decoder/libs/CMakeLists.txt index da1e1d33e..4bf84f6d7 100644 --- a/src/algorithms/telemetry_decoder/libs/CMakeLists.txt +++ b/src/algorithms/telemetry_decoder/libs/CMakeLists.txt @@ -9,17 +9,19 @@ add_subdirectory(libswiftcnav) set(TELEMETRY_DECODER_LIB_SOURCES tlm_conf.cc + tlm_crc_stats.cc tlm_utils.cc viterbi_decoder.cc - tlm_crc_stats.cc + viterbi_decoder_sbas.cc ) set(TELEMETRY_DECODER_LIB_HEADERS + tlm_conf.h - viterbi_decoder.h - convolutional.h - tlm_utils.h tlm_crc_stats.h + tlm_utils.h + viterbi_decoder.h + viterbi_decoder_sbas.h ) list(SORT TELEMETRY_DECODER_LIB_HEADERS) @@ -42,9 +44,8 @@ else() endif() target_link_libraries(telemetry_decoder_libs - PUBLIC - Volkgnsssdr::volkgnsssdr PRIVATE + Volkgnsssdr::volkgnsssdr algorithms_libs Gflags::gflags Glog::glog diff --git a/src/algorithms/telemetry_decoder/libs/convolutional.h b/src/algorithms/telemetry_decoder/libs/convolutional.h deleted file mode 100644 index 2dae42adb..000000000 --- a/src/algorithms/telemetry_decoder/libs/convolutional.h +++ /dev/null @@ -1,278 +0,0 @@ -/*! - * \file convolutional.h - * \brief General functions used to implement convolutional encoding. - * \author Matthew C. Valenti, 2006-2008. - * \author C. Fernandez-Prades, 2019. - * - * ----------------------------------------------------------------------------- - * - * GNSS-SDR is a Global Navigation Satellite System software-defined receiver. - * This file is part of GNSS-SDR. - * - * Copyright (C) 2006-2008 Matthew C. Valenti - * Copyright (C) 2019 C. Fernandez-Prades - * SPDX-License-Identifier: GPL-3.0-or-later - * - * This file is a derived work of the original file, which had this note: - * - * The functions in this file are part of the Iterative Solutions - * Coded Modulation Library. The Iterative Solutions Coded Modulation - * Library is free software; you can redistribute it and/or modify it - * under the terms of the GNU Lesser General Public License as published - * by the Free Software Foundation; either version 2.1 of the License, - * or (at your option) any later version. - * - * This library 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 - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - */ - -#ifndef GNSS_SDR_CONVOLUTIONAL_H -#define GNSS_SDR_CONVOLUTIONAL_H - -#include -#include - -/** \addtogroup Telemetry_Decoder - * \{ */ -/** \addtogroup Telemetry_Decoder_libs - * Utilities for the decoding of GNSS navigation messages. - * \{ */ - - -/* define constants used throughout the library */ -const float MAXLOG = 1e7; /* Define infinity */ - - -/*! - * \brief Determines if a symbol has odd (1) or even (0) parity - * Output parameters: - * \return (returned int): The symbol's parity = 1 for odd and 0 for even - * - * \param[in] symbol The integer-valued symbol - * \param[in] length The highest bit position in the symbol - * - * This function is used by nsc_enc_bit(), rsc_enc_bit(), and rsc_tail() - */ -inline int parity_counter(int symbol, int length) -{ - int counter; - int temp_parity = 0; - - for (counter = 0; counter < length; counter++) - { - temp_parity = temp_parity ^ (symbol & 1); - symbol = symbol >> 1; - } - return (temp_parity); -} - - -/*! - * \brief Convolutionally encodes a single bit using a rate 1/n encoder. - * Takes in one input bit at a time, and produces a n-bit output. - * - * \param[in] input The input data bit (i.e. a 0 or 1). - * \param[in] state_in The starting state of the encoder (an int from 0 to 2^m-1). - * \param[in] g[] An n-element vector containing the code generators in binary form. - * \param[in] KK The constraint length of the convolutional code. - * \param[out] output_p[] An n-element vector containing the encoded bits. - * \param[out] state_out_p[] An integer containing the final state of the encoder - * (i.e. the state after encoding this bit) - * - * This function is used by nsc_transit() - */ -inline int nsc_enc_bit(int state_out_p[], - int input, - int state_in, - const int g[], - int KK, - int nn) -{ - /* declare variables */ - int state, i; - int out_ = 0; - - /* create a word made up of state and new input */ - state = (input << (KK - 1)) ^ state_in; - - /* AND the word with the generators */ - for (i = 0; i < nn; i++) - { - /* update output symbol */ - out_ = (out_ << 1) + parity_counter(state & g[i], KK); - } - - /* shift the state to make the new state */ - state_out_p[0] = state >> 1; - return (out_); -} - - -/*! - * \brief Function that creates the transit and output vectors - */ -inline void nsc_transit(int output_p[], - int trans_p[], - int input, - int g[], - int KK, - int nn) -{ - int nextstate[1]; - int state, states; - states = (1 << (KK - 1)); /* The number of states: 2^mm */ - - /* Determine the output and next state for each possible starting state */ - for (state = 0; state < states; state++) - { - output_p[state] = nsc_enc_bit(nextstate, input, state, g, KK, nn); - trans_p[state] = nextstate[0]; - } - return; -} - - -/*! - * \brief Computes the branch metric used for decoding. - * \return (returned float) The metric between the hypothetical symbol and the received vector - * \param[in] rec_array The received vector, of length nn - * \param[in] symbol The hypothetical symbol - * \param[in] nn The length of the received vector - * - */ -inline float Gamma(const float rec_array[], - int symbol, - int nn) -{ - float rm = 0; - int i; - int mask = 1; - - for (i = 0; i < nn; i++) - { - if (symbol & mask) - { - rm += rec_array[nn - i - 1]; - } - mask = mask << 1; - } - return (rm); -} - - -/*! - * \brief Uses the Viterbi algorithm to perform hard-decision decoding of a convolutional code. - * \param[in] out0[] The output bits for each state if input is a 0. - * \param[in] state0[] The next state if input is a 0. - * \param[in] out1[] The output bits for each state if input is a 1. - * \param[in] state1[] The next state if input is a 1. - * \param[in] r[] The received signal in LLR-form. For BPSK, must be in form r = 2*a*y/(sigma^2). - * \param[in] KK The constraint length of the convolutional code. - * \param[in] LL The number of data bits. - * \param[out] output_u_int[] Hard decisions on the data bits - * - */ -inline void Viterbi(int output_u_int[], - const int out0[], - const int state0[], - const int out1[], - const int state1[], - const float input_c[], - int KK, - int nn, - int LL) -{ - int i, t, state, mm, states; - int number_symbols; - uint32_t max_index; - float metric; - float max_val; - - /* some derived constants */ - mm = KK - 1; - states = 1 << mm; /* 2^mm */ - number_symbols = 1 << nn; /* 2^nn */ - - std::vector prev_section(states, -MAXLOG); - std::vector next_section(states, -MAXLOG); - std::vector prev_bit(states * (LL + mm), 0); - std::vector prev_state(states * (LL + mm), 0); - std::vector rec_array(nn); - std::vector metric_c(number_symbols); - - prev_section[0] = 0.0; /* start in all-zeros state */ - - /* go through trellis */ - for (t = 0; t < LL + mm; t++) - { - rec_array.assign(input_c + nn * t, input_c + nn * t + (nn - 1)); - - /* precompute all possible branch metrics */ - for (i = 0; i < number_symbols; i++) - { - metric_c[i] = Gamma(rec_array.data(), i, nn); - } - - /* step through all states */ - for (state = 0; state < states; state++) - { - /* hypothesis: info bit is a zero */ - metric = prev_section[state] + metric_c[out0[state]]; - - /* store new metric if more than metric in storage */ - if (metric > next_section[state0[state]]) - { - next_section[state0[state]] = metric; - prev_state[t * states + state0[state]] = state; - prev_bit[t * states + state0[state]] = 0; - } - - /* hypothesis: info bit is a one */ - metric = prev_section[state] + metric_c[out1[state]]; - - /* store new metric if more than metric in storage */ - if (metric > next_section[state1[state]]) - { - next_section[state1[state]] = metric; - prev_state[t * states + state1[state]] = state; - prev_bit[t * states + state1[state]] = 1; - } - } - - /* normalize */ - volk_gnsssdr_32f_index_max_32u(&max_index, next_section.data(), states); - max_val = next_section[max_index]; - - for (state = 0; state < states; state++) - { - prev_section[state] = next_section[state] - max_val; - next_section[state] = -MAXLOG; - } - } - - /* trace-back operation */ - state = 0; - - /* tail, no need to output */ - for (t = LL + mm - 1; t >= LL; t--) - { - state = prev_state[t * states + state]; - } - - for (t = LL - 1; t >= 0; t--) - { - output_u_int[t] = prev_bit[t * states + state]; - state = prev_state[t * states + state]; - } -} - - -/** \} */ -/** \} */ -#endif // GNSS_SDR_CONVOLUTIONAL_H diff --git a/src/algorithms/telemetry_decoder/libs/viterbi_decoder.cc b/src/algorithms/telemetry_decoder/libs/viterbi_decoder.cc index 93d8d1fd6..77fdf987c 100644 --- a/src/algorithms/telemetry_decoder/libs/viterbi_decoder.cc +++ b/src/algorithms/telemetry_decoder/libs/viterbi_decoder.cc @@ -1,569 +1,189 @@ /*! * \file viterbi_decoder.cc - * \brief Implementation of a Viterbi decoder class based on the Iterative Solutions - * Coded Modulation Library by Matthew C. Valenti - * \author Daniel Fehr 2013. daniel.co(at)bluewin.ch + * \brief Class that implements a Viterbi decoder + * \author Carles Fernandez, 2021. cfernandez(at)cttc.es * * ----------------------------------------------------------------------------- * * GNSS-SDR is a Global Navigation Satellite System software-defined receiver. * This file is part of GNSS-SDR. * - * Copyright (C) 2010-2020 (see AUTHORS file for a list of contributors) + * Copyright (C) 2010-2021 (see AUTHORS file for a list of contributors) * SPDX-License-Identifier: GPL-3.0-or-later * * ----------------------------------------------------------------------------- */ #include "viterbi_decoder.h" -#include -#include // for fill_n -#include // for operator<<, basic_ostream, char_traits +#include // for volk_gnsssdr_32f_index_max_32u +#include // for std::copy -// logging -#define EVENT 2 // logs important events which don't occur every block -#define FLOW 3 // logs the function calls of block processing functions -#define BLOCK 4 // once per block -#define SAMPLE 5 // about one log entry per sample -#define LMORE 6 // many entries per sample / very specific stuff - - -const float MAXLOG = 1e7; /* Define infinity */ - -Viterbi_Decoder::Viterbi_Decoder(const int g_encoder[], const int KK, const int nn) +Viterbi_Decoder::Viterbi_Decoder(int32_t KK, int32_t nn, int32_t LL, const std::array& g) { - d_nn = nn; // Coding rate 1/n - d_KK = KK; // Constraint Length - - // derived code properties + d_KK = KK; + d_nn = nn; + d_LL = LL; d_mm = d_KK - 1; - d_states = static_cast(1U << d_mm); /* 2^mm */ - d_number_symbols = static_cast(1U << d_nn); /* 2^nn */ - - /* create appropriate transition matrices (trellis) */ - d_out0.reserve(d_states); - d_out1.reserve(d_states); - d_state0.reserve(d_states); - d_state1.reserve(d_states); - - nsc_transit(d_out0.data(), d_state0.data(), 0, g_encoder, d_KK, d_nn); - nsc_transit(d_out1.data(), d_state1.data(), 1, g_encoder, d_KK, d_nn); - - // initialise trellis state - d_trellis_state_is_initialised = false; - Viterbi_Decoder::init_trellis_state(); + d_states = 1 << d_mm; // 2^d_mm + d_number_symbols = 1 << d_nn; // 2^d_nn + d_prev_section = std::vector(d_states, -d_MAXLOG); + d_next_section = std::vector(d_states, -d_MAXLOG); + d_rec_array = std::vector(d_nn); + d_metric_c = std::vector(d_number_symbols); + d_prev_bit = std::vector(d_states * (d_LL + d_mm)); + d_prev_state = std::vector(d_states * (d_LL + d_mm)); + d_g = g; } -void Viterbi_Decoder::reset() +void Viterbi_Decoder::decode(std::vector& output_u_int, + const std::vector& out0, + const std::vector& state0, + const std::vector& out1, + const std::vector& state1, + const std::vector& input_c) { - init_trellis_state(); -} - - -/* Function decode_block() - Description: Uses the Viterbi algorithm to perform hard-decision decoding of a convolutional code. - Input parameters: - r[] The received signal in LLR-form. For BPSK, must be in form r = 2*a*y/(sigma^2). - LL The number of data bits to be decoded (doesn't include the mm zero-tail-bits) - Output parameters: - output_u_int[] Hard decisions on the data bits (without the mm zero-tail-bits) - */ -float Viterbi_Decoder::decode_block(const double input_c[], int output_u_int[], const int LL) -{ - VLOG(FLOW) << "decode_block(): LL=" << LL; - - // init - init_trellis_state(); - // do add compare select - do_acs(input_c, LL + d_mm); - // tail, no need to output -> traceback, but don't decode - const int state = do_traceback(d_mm); - // traceback and decode - const int decoding_length_mismatch = do_tb_and_decode(d_mm, LL, state, output_u_int, d_indicator_metric); - - VLOG(FLOW) << "decoding length mismatch: " << decoding_length_mismatch; - - return d_indicator_metric; -} - - -float Viterbi_Decoder::decode_continuous(const double sym[], - const int traceback_depth, - int bits[], - const int nbits_requested, - int& nbits_decoded) -{ - VLOG(FLOW) << "decode_continuous(): nbits_requested=" << nbits_requested; - - // do add compare select - do_acs(sym, nbits_requested); - // the ML sequence in the newest part of the trellis can not be decoded - // since it depends on the future values -> traceback, but don't decode - const int state = do_traceback(traceback_depth); - // traceback and decode - const int decoding_length_mismatch = do_tb_and_decode(traceback_depth, nbits_requested, state, bits, d_indicator_metric); - nbits_decoded = nbits_requested + decoding_length_mismatch; - - VLOG(FLOW) << "decoding length mismatch (continuous decoding): " << decoding_length_mismatch; - - return d_indicator_metric; -} - - -void Viterbi_Decoder::init_trellis_state() -{ - int state; - // if trellis state has been initialised, free old state memory - if (d_trellis_state_is_initialised) - { - // init trellis state - d_pm_t.clear(); - d_rec_array.clear(); - d_metric_c.clear(); - } - - // reserve new trellis state memory - d_pm_t.reserve(d_states); - d_trellis_paths = std::deque(); - d_rec_array.reserve(d_nn); - d_metric_c.reserve(d_number_symbols); - d_trellis_state_is_initialised = true; - - /* initialize trellis */ - for (state = 0; state < d_states; state++) - { - d_pm_t[state] = -MAXLOG; - // d_pm_t_next[state] = -MAXLOG; - } - d_pm_t[0] = 0; /* start in all-zeros state */ - - d_indicator_metric = 0; -} - - -int Viterbi_Decoder::do_acs(const double sym[], int nbits) -{ - int t; - int i; - int state_at_t; + int32_t i, t, state; + uint32_t max_index; float metric; float max_val; - std::vector pm_t_next(d_states); - /* t: - * - state: state at t - * - d_prev_section[state_at_t]: path metric at t for state state_at_t - * - d_out0[state_at_t]: sent symbols for a data bit 0 if state is state_at_t at time t - * - */ + d_prev_section[0] = 0.0; // start in all-zeros state - for (state_at_t = 0; state_at_t < d_states; state_at_t++) + // go through trellis + for (t = 0; t < d_LL + d_mm; t++) { - pm_t_next[state_at_t] = -MAXLOG; - } + std::copy(input_c.begin() + d_nn * t, input_c.begin() + d_nn * t + (d_nn - 1), d_rec_array.begin()); - /* go through trellis */ - for (t = 0; t < nbits; t++) - { - /* Temporarily store the received symbols current decoding step */ - for (i = 0; i < d_nn; i++) - { - d_rec_array[i] = static_cast(sym[d_nn * t + i]); - } - - /* precompute all possible branch metrics */ + // precompute all possible branch metrics for (i = 0; i < d_number_symbols; i++) { - d_metric_c[i] = gamma(d_rec_array.data(), i, d_nn); - VLOG(LMORE) << "metric for (tx_sym=" << i << "|ry_sym=(" << d_rec_array[0] << ", " << d_rec_array[1] << ") = " << d_metric_c[i]; + d_metric_c[i] = this->Gamma(i); } - // find the survivor branches leading the trellis states at t+1 - Prev next_trellis_states(d_states, t + 1); - /* step through all states */ - for (state_at_t = 0; state_at_t < d_states; state_at_t++) + // step through all states + for (state = 0; state < d_states; state++) { - const int next_state_if_0 = d_state0[state_at_t]; - const int next_state_if_1 = d_state1[state_at_t]; + // hypothesis: info bit is a zero + metric = d_prev_section[state] + d_metric_c[out0[state]]; - /* hypothesis: info bit is a zero */ - const float bm_0 = d_metric_c[d_out0[state_at_t]]; - metric = d_pm_t[state_at_t] + bm_0; // path metric + zerobranch metric - - /* store new metric if more than metric in storage */ - if (metric > pm_t_next[next_state_if_0]) + // store new metric if more than metric in storage + if (metric > d_next_section[state0[state]]) { - pm_t_next[next_state_if_0] = metric; - next_trellis_states.set_current_state_as_ancestor_of_next_state(next_state_if_0, state_at_t); - next_trellis_states.set_decoded_bit_for_next_state(next_state_if_0, 0); - next_trellis_states.set_survivor_branch_metric_of_next_state(next_state_if_0, bm_0); + d_next_section[state0[state]] = metric; + d_prev_state[t * d_states + state0[state]] = state; + d_prev_bit[t * d_states + state0[state]] = 0; } - /* hypothesis: info bit is a one */ - const float bm_1 = d_metric_c[d_out1[state_at_t]]; - metric = d_pm_t[state_at_t] + bm_1; // path metric + onebranch metric + // hypothesis: info bit is a one + metric = d_prev_section[state] + d_metric_c[out1[state]]; - /* store new metric if more than metric in storage */ - if (metric > pm_t_next[next_state_if_1]) + // store new metric if more than metric in storage + if (metric > d_next_section[state1[state]]) { - pm_t_next[next_state_if_1] = metric; - next_trellis_states.set_current_state_as_ancestor_of_next_state(next_state_if_1, state_at_t); - next_trellis_states.set_decoded_bit_for_next_state(next_state_if_1, 1); - next_trellis_states.set_survivor_branch_metric_of_next_state(next_state_if_1, bm_1); + d_next_section[state1[state]] = metric; + d_prev_state[t * d_states + state1[state]] = state; + d_prev_bit[t * d_states + state1[state]] = 1; } } - d_trellis_paths.push_front(next_trellis_states); + // normalize + volk_gnsssdr_32f_index_max_32u(&max_index, d_next_section.data(), d_states); + max_val = d_next_section[max_index]; - /* normalize -> afterwards, the largest metric value is always 0 */ - // max_val = 0; - max_val = -MAXLOG; - for (state_at_t = 0; state_at_t < d_states; state_at_t++) + for (state = 0; state < d_states; state++) { - if (pm_t_next[state_at_t] > max_val) - { - max_val = pm_t_next[state_at_t]; - } - } - VLOG(LMORE) << "max_val at t=" << t << ": " << max_val; - for (state_at_t = 0; state_at_t < d_states; state_at_t++) - { - d_pm_t[state_at_t] = pm_t_next[state_at_t] - max_val; - pm_t_next[state_at_t] = -MAXLOG; + d_prev_section[state] = d_next_section[state] - max_val; + d_next_section[state] = -d_MAXLOG; } } - return t; + // trace-back operation + state = 0; + + // tail, no need to output + for (t = d_LL + d_mm - 1; t >= d_LL; t--) + { + state = d_prev_state[t * d_states + state]; + } + + for (t = d_LL - 1; t >= 0; t--) + { + output_u_int[t] = d_prev_bit[t * d_states + state]; + state = d_prev_state[t * d_states + state]; + } } -int Viterbi_Decoder::do_traceback(size_t traceback_length) +void Viterbi_Decoder::nsc_transit(std::vector& output_p, + std::vector& trans_p, + int32_t input) const { - // traceback_length is in bits - int state; - std::deque::iterator it; + int32_t nextstate; + int32_t state, states; + states = (1 << (d_KK - 1)); // The number of states: 2^d_mm - VLOG(FLOW) << "do_traceback(): traceback_length=" << traceback_length << '\n'; - - if (d_trellis_paths.size() < traceback_length) + // Determine the output and next state for each possible starting state + for (state = 0; state < states; state++) { - traceback_length = d_trellis_paths.size(); + output_p[state] = nsc_enc_bit(&nextstate, input, state); + trans_p[state] = nextstate; } - - state = 0; // maybe start not at state 0, but at state with best metric - for (it = d_trellis_paths.begin(); it < d_trellis_paths.begin() + traceback_length; ++it) - { - state = it->get_anchestor_state_of_current_state(state); - } - return state; + return; } -int Viterbi_Decoder::do_tb_and_decode(int traceback_length, int requested_decoding_length, int state, int output_u_int[], float& indicator_metric) -{ - int n_of_branches_for_indicator_metric = 500; - std::deque::iterator it; - int n_im = 0; - - VLOG(FLOW) << "do_tb_and_decode(): requested_decoding_length=" << requested_decoding_length; - // decode only decode_length bits -> overstep newer bits which are too much - const int decoding_length_mismatch = static_cast(d_trellis_paths.size()) - (traceback_length + requested_decoding_length); - VLOG(BLOCK) << "decoding_length_mismatch=" << decoding_length_mismatch; - const int overstep_length = decoding_length_mismatch >= 0 ? decoding_length_mismatch : 0; - VLOG(BLOCK) << "overstep_length=" << overstep_length; - - for (it = d_trellis_paths.begin() + traceback_length; - it < d_trellis_paths.begin() + traceback_length + overstep_length; ++it) - { - state = it->get_anchestor_state_of_current_state(state); - } - int t_out = static_cast(d_trellis_paths.end() - (d_trellis_paths.begin() + traceback_length + overstep_length) - 1); // requested_decoding_length-1; - indicator_metric = 0; - for (it = d_trellis_paths.begin() + traceback_length + overstep_length; it < d_trellis_paths.end(); ++it) - { - if (it - (d_trellis_paths.begin() + traceback_length + overstep_length) < n_of_branches_for_indicator_metric) - { - n_im++; - indicator_metric += it->get_metric_of_current_state(state); - VLOG(SAMPLE) << "@t=" << it->get_t() << " b=" << it->get_bit_of_current_state(state) << " sm=" << indicator_metric << " d=" << it->get_metric_of_current_state(state); - } - output_u_int[t_out] = it->get_bit_of_current_state(state); - state = it->get_anchestor_state_of_current_state(state); - t_out--; - } - if (n_im > 0) - { - indicator_metric /= static_cast(n_im); - } - - VLOG(BLOCK) << "indicator metric: " << indicator_metric; - // remove old states - if (d_trellis_paths.begin() + traceback_length + overstep_length <= d_trellis_paths.end()) - { - d_trellis_paths.erase(d_trellis_paths.begin() + traceback_length + overstep_length, d_trellis_paths.end()); - } - return decoding_length_mismatch; -} - - -/* function Gamma() - - Description: Computes the branch metric used for decoding. - - Output parameters: - (returned float) The metric between the hypothetical symbol and the recevieved vector - - Input parameters: - rec_array The received vector, of length nn - symbol The hypothetical symbol - nn The length of the received vector - - This function is used by siso() */ -float Viterbi_Decoder::gamma(const float rec_array[], int symbol, int nn) +float Viterbi_Decoder::Gamma(int32_t symbol) const { float rm = 0; - int i; - unsigned int mask = 1U; - float txsym; + int32_t i; + int32_t mask = 1; - for (i = 0; i < nn; i++) + for (i = 0; i < d_nn; i++) { - // if (symbol & mask) rm += rec_array[nn - i - 1]; - txsym = symbol & mask ? 1 : -1; - rm += txsym * rec_array[nn - i - 1]; - mask = mask << 1U; + if (symbol & mask) + { + rm += d_rec_array[d_nn - i - 1]; + } + mask = mask << 1; } - // rm = rm > 50 ? rm : -1000; - return rm; } -/* function that creates the transit and output vectors */ -void Viterbi_Decoder::nsc_transit(int output_p[], int trans_p[], int input, const int g[], - int KK, int nn) +int32_t Viterbi_Decoder::parity_counter(int32_t symbol, int32_t length) const { - int nextstate[1]; - int state; - const int states = static_cast(1U << (KK - 1)); /* The number of states: 2^mm */ + int32_t counter; + int32_t temp_parity = 0; - /* Determine the output and next state for each possible starting state */ - for (state = 0; state < states; state++) - { - output_p[state] = nsc_enc_bit(nextstate, input, state, g, KK, nn); - trans_p[state] = nextstate[0]; - } -} - - -/* Function nsc_enc_bit() - - Description: Convolutionally encodes a single bit using a rate 1/n encoder. - Takes in one input bit at a time, and produces a n-bit output. - - Input parameters: - input The input data bit (i.e. a 0 or 1). - state_in The starting state of the encoder (an int from 0 to 2^m-1). - g[] An n-element vector containing the code generators in binary form. - KK The constraint length of the convolutional code. - nn number of symbols bits per input bits (rate 1/nn) - - Output parameters: - output_p[] An n-element vector containing the encoded bits. - state_out_p[] An integer containing the final state of the encoder - (i.e. the state after encoding this bit) - - This function is used by rsc_encode(), nsc_transit(), rsc_transit(), and nsc_transit() */ -int Viterbi_Decoder::nsc_enc_bit(int state_out_p[], int input, int state_in, - const int g[], int KK, int nn) -{ - /* declare variables */ - int state; - int i; - int out = 0; - - /* create a word made up of state and new input */ - state = (input << (KK - 1)) ^ state_in; - - /* AND the word with the generators */ - for (i = 0; i < nn; i++) - { - /* update output symbol */ - out = (out << 1) + parity_counter(state & g[i], KK); - } - - /* shift the state to make the new state */ - state_out_p[0] = state >> 1; - return (out); -} - - -/* function parity_counter() - - Description: Determines if a symbol has odd (1) or even (0) parity - - Output parameters: - (returned int): The symbol's parity = 1 for odd and 0 for even - - Input parameters: - symbol: The integer-valued symbol - length: The highest bit position in the symbol - - This function is used by nsc_enc_bit(), rsc_enc_bit(), and rsc_tail() */ -int Viterbi_Decoder::parity_counter(int symbol, int length) -{ - int counter; - unsigned int temp_parity = 0; for (counter = 0; counter < length; counter++) { - temp_parity = temp_parity ^ (symbol & 1U); - symbol = symbol >> 1U; + temp_parity = temp_parity ^ (symbol & 1); + symbol = symbol >> 1; } - return static_cast(temp_parity); + return temp_parity; } -// prev helper class -Viterbi_Decoder::Prev::Prev(int states, int t) +int32_t Viterbi_Decoder::nsc_enc_bit(int32_t* state_out_p, + int32_t input, + int32_t state_in) const { - this->t = t; - num_states = states; - state.reserve(num_states); - v_bit.reserve(num_states); - v_metric.reserve(num_states); - refcount = 1; - std::fill_n(state.begin(), num_states, 0); - std::fill_n(v_bit.begin(), num_states, 0); - std::fill_n(v_metric.begin(), num_states, 0.0F); -} + // declare variables + int32_t state, i; + int32_t out_ = 0; + // create a word made up of state and new input + state = (input << (d_KK - 1)) ^ state_in; -// copy constructor -Viterbi_Decoder::Prev::Prev(const Prev& prev) -{ - refcount = prev.refcount; - refcount++; - t = prev.t; - state = prev.state; - num_states = prev.num_states; - v_bit = prev.v_bit; - v_metric = prev.v_metric; - VLOG(LMORE) << "Prev(" - << "?" - << ", " << t << ")" - << " copy, new refcount = " << refcount; -} - - -// assignment constructor -Viterbi_Decoder::Prev& Viterbi_Decoder::Prev::operator=(const Prev& other) -{ - // check for self-assignment - if (&other == this) + // AND the word with the generators + for (i = 0; i < d_nn; i++) { - return *this; + // update output symbol + out_ = (out_ << 1) + parity_counter(state & d_g[i], d_KK); } - // handle old resources - if (refcount != 1) - { // this object is not anymore using them - refcount--; - } - - // increase ref counter for this resource set - refcount = other.refcount; - refcount++; - - // take over resources - t = other.t; - state = other.state; - v_bit = other.v_bit; - v_metric = other.v_metric; - - VLOG(LMORE) << "Prev(" - << "?" - << ", " << t << ")" - << " assignment, new refcount = " << refcount; - return *this; -} - - -Viterbi_Decoder::Prev::~Prev() -{ - if (refcount != 1) - { - refcount--; - VLOG(LMORE) << "~Prev(" - << "?" - << ", " << t << ")" - << " destructor after copy, new refcount = " << refcount; - } -} - - -int Viterbi_Decoder::Prev::get_anchestor_state_of_current_state(int current_state) const -{ - // std::cout << "get prev state: for state " << current_state << " at time " << t << ", the prev state at time " << t - 1 << " is " << state[current_state] << '\n'; - if (num_states > current_state) - { - return state[current_state]; - } - // std::cout << "alarm " << "num_states=" << num_states << " current_state=" << current_state << '\n'; - // return state[current_state]; - return 0; -} - - -int Viterbi_Decoder::Prev::get_bit_of_current_state(int current_state) const -{ - // std::cout << "get prev bit : for state " << current_state << " at time " << t << ", the send bit is " << bit[current_state] << '\n'; - if (num_states > current_state) - { - return v_bit[current_state]; - } - return 0; -} - - -float Viterbi_Decoder::Prev::get_metric_of_current_state(int current_state) const -{ - if (num_states > current_state) - { - return v_metric[current_state]; - } - return 0; -} - - -int Viterbi_Decoder::Prev::get_t() const -{ - return t; -} - - -void Viterbi_Decoder::Prev::set_current_state_as_ancestor_of_next_state(int next_state, int current_state) -{ - if (num_states > next_state) - { - state[next_state] = current_state; - } -} - - -void Viterbi_Decoder::Prev::set_decoded_bit_for_next_state(int next_state, int bit) -{ - if (num_states > next_state) - { - this->v_bit[next_state] = bit; - } -} - - -void Viterbi_Decoder::Prev::set_survivor_branch_metric_of_next_state(int next_state, float metric) -{ - if (num_states > next_state) - { - this->v_metric[next_state] = metric; - } + // shift the state to make the new state + *state_out_p = state >> 1; + return out_; } diff --git a/src/algorithms/telemetry_decoder/libs/viterbi_decoder.h b/src/algorithms/telemetry_decoder/libs/viterbi_decoder.h index a25c0b7b0..ca8a7b896 100644 --- a/src/algorithms/telemetry_decoder/libs/viterbi_decoder.h +++ b/src/algorithms/telemetry_decoder/libs/viterbi_decoder.h @@ -1,15 +1,14 @@ /*! * \file viterbi_decoder.h - * \brief Interface of a Viterbi decoder class based on the Iterative Solutions - * Coded Modulation Library by Matthew C. Valenti - * \author Daniel Fehr 2013. daniel.co(at)bluewin.ch + * \brief Class that implements a Viterbi decoder + * \author Carles Fernandez, 2021. cfernandez(at)cttc.es * * ----------------------------------------------------------------------------- * * GNSS-SDR is a Global Navigation Satellite System software-defined receiver. * This file is part of GNSS-SDR. * - * Copyright (C) 2010-2020 (see AUTHORS file for a list of contributors) + * Copyright (C) 2010-2021 (see AUTHORS file for a list of contributors) * SPDX-License-Identifier: GPL-3.0-or-later * * ----------------------------------------------------------------------------- @@ -18,106 +17,110 @@ #ifndef GNSS_SDR_VITERBI_DECODER_H #define GNSS_SDR_VITERBI_DECODER_H -#include // for size_t -#include +#include +#include #include /** \addtogroup Telemetry_Decoder * \{ */ -/** \addtogroup Telemetry_Decoder_libs telemetry_decoder_libs +/** \addtogroup Telemetry_Decoder_libs + * Utilities for the decoding of GNSS navigation messages. * \{ */ -/*! - * \brief Class that implements a Viterbi decoder - */ class Viterbi_Decoder { public: - Viterbi_Decoder(const int g_encoder[], const int KK, const int nn); - ~Viterbi_Decoder() = default; - void reset(); + /*! + * \brief Constructor of a Viterbi decoder + * \param[in] KK Constraint Length + * \param[in] nn Coding rate 1/n + * \param[in] LL Data length + * \param[in] g Polynomial G1 and G2 + */ + Viterbi_Decoder(int32_t KK, int32_t nn, int32_t LL, const std::array& g); /*! * \brief Uses the Viterbi algorithm to perform hard-decision decoding of a convolutional code. + * \param[out] output_u_int Hard decisions on the data bits + * \param[in] out0 The output bits for each state if input is a 0. + * \param[in] state0 The next state if input is a 0. + * \param[in] out1 The output bits for each state if input is a 1. + * \param[in] state1 The next state if input is a 1. + * \param[in] input_c The received signal in LLR-form. For BPSK, must be in form r = 2*a*y/(sigma^2). * - * \param[in] input_c[] The received signal in LLR-form. For BPSK, must be in form r = 2*a*y/(sigma^2). - * \param[in] LL The number of data bits to be decoded (does not include the mm zero-tail-bits) - * - * \return output_u_int[] Hard decisions on the data bits (without the mm zero-tail-bits) */ - float decode_block(const double input_c[], int* output_u_int, const int LL); + void decode(std::vector& output_u_int, + const std::vector& out0, + const std::vector& state0, + const std::vector& out1, + const std::vector& state1, + const std::vector& input_c); - float decode_continuous(const double sym[], const int traceback_depth, int bits[], - const int nbits_requested, int& nbits_decoded); + /*! + * \brief Function that creates the transit and output vectors + */ + void nsc_transit(std::vector& output_p, + std::vector& trans_p, + int32_t input) const; private: - class Prev - { - public: - int num_states; - Prev(int states, int t); - Prev(const Prev& prev); - Prev& operator=(const Prev& other); - ~Prev(); + /* + * Computes the branch metric used for decoding. + * \return (returned float) The metric between the hypothetical symbol and the received vector + * \param[in] symbol The hypothetical symbol + * + */ + float Gamma(int32_t symbol) const; - int get_anchestor_state_of_current_state(int current_state) const; - int get_bit_of_current_state(int current_state) const; - float get_metric_of_current_state(int current_state) const; - int get_t() const; - void set_current_state_as_ancestor_of_next_state(int next_state, int current_state); - void set_decoded_bit_for_next_state(int next_state, int bit); - void set_survivor_branch_metric_of_next_state(int next_state, float metric); + /* + * Determines if a symbol has odd (1) or even (0) parity + * Output parameters: + * \return (returned int): The symbol's parity = 1 for odd and 0 for even + * + * \param[in] symbol The integer-valued symbol + * \param[in] length The highest bit position in the symbol + * + * This function is used by nsc_enc_bit(), rsc_enc_bit(), and rsc_tail() + */ + int32_t parity_counter(int32_t symbol, int32_t length) const; - private: - std::vector v_metric; - std::vector state; - std::vector v_bit; - int t; - int refcount; - }; + /* + * Convolutionally encodes a single bit using a rate 1/n encoder. + * Takes in one input bit at a time, and produces a n-bit output. + * + * \param[in] input The input data bit (i.e. a 0 or 1). + * \param[in] state_in The starting state of the encoder (an int from 0 to 2^m-1). + * \param[in] g An n-element vector containing the code generators in binary form. + * \param[out] output_p[] An n-element vector containing the encoded bits. + * \param[out] state_out_p[] An integer containing the final state of the encoder + * (i.e. the state after encoding this bit) + * + * This function is used by nsc_transit() + */ + int32_t nsc_enc_bit(int32_t* state_out_p, + int32_t input, + int32_t state_in) const; - // operations on the trellis (change decoder state) - void init_trellis_state(); - int do_acs(const double sym[], int nbits); - int do_traceback(std::size_t traceback_length); - int do_tb_and_decode(int traceback_length, int requested_decoding_length, int state, int output_u_int[], float& indicator_metric); + std::vector d_prev_section{}; + std::vector d_next_section{}; - // branch metric function - float gamma(const float rec_array[], int symbol, int nn); + std::vector d_rec_array{}; + std::vector d_metric_c{}; + std::vector d_prev_bit{}; + std::vector d_prev_state{}; + std::array d_g{}; - // trellis generation - void nsc_transit(int output_p[], int trans_p[], int input, const int g[], int KK, int nn); - int nsc_enc_bit(int state_out_p[], int input, int state_in, const int g[], int KK, int nn); - int parity_counter(int symbol, int length); + float d_MAXLOG = 1e7; // Define infinity + int32_t d_KK{}; + int32_t d_nn{}; + int32_t d_LL{}; - // trellis state - std::deque d_trellis_paths; - std::vector d_pm_t; - std::vector d_metric_c; /* Set of all possible branch metrics */ - std::vector d_rec_array; /* Received values for one trellis section */ - - // trellis definition - std::vector d_out0; - std::vector d_state0; - std::vector d_out1; - std::vector d_state1; - - // measures - float d_indicator_metric; - - // code properties - int d_KK; - int d_nn; - - // derived code properties - int d_mm; - int d_states; - int d_number_symbols; - bool d_trellis_state_is_initialised; + int32_t d_mm{}; + int32_t d_states{}; + int32_t d_number_symbols{}; }; - /** \} */ /** \} */ #endif // GNSS_SDR_VITERBI_DECODER_H diff --git a/src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.cc b/src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.cc new file mode 100644 index 000000000..3b2d2065f --- /dev/null +++ b/src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.cc @@ -0,0 +1,569 @@ +/*! + * \file viterbi_decoder_sbas.cc + * \brief Implementation of a Viterbi decoder class based on the Iterative Solutions + * Coded Modulation Library by Matthew C. Valenti + * \author Daniel Fehr 2013. daniel.co(at)bluewin.ch + * + * ----------------------------------------------------------------------------- + * + * GNSS-SDR is a Global Navigation Satellite System software-defined receiver. + * This file is part of GNSS-SDR. + * + * Copyright (C) 2010-2020 (see AUTHORS file for a list of contributors) + * SPDX-License-Identifier: GPL-3.0-or-later + * + * ----------------------------------------------------------------------------- + */ + +#include "viterbi_decoder_sbas.h" +#include +#include // for fill_n +#include // for operator<<, basic_ostream, char_traits + +// logging +#define EVENT 2 // logs important events which don't occur every block +#define FLOW 3 // logs the function calls of block processing functions +#define BLOCK 4 // once per block +#define SAMPLE 5 // about one log entry per sample +#define LMORE 6 // many entries per sample / very specific stuff + + +const float MAXLOG = 1e7; /* Define infinity */ + +Viterbi_Decoder_Sbas::Viterbi_Decoder_Sbas(const int g_encoder[], const int KK, const int nn) +{ + d_nn = nn; // Coding rate 1/n + d_KK = KK; // Constraint Length + + // derived code properties + d_mm = d_KK - 1; + d_states = static_cast(1U << d_mm); /* 2^mm */ + d_number_symbols = static_cast(1U << d_nn); /* 2^nn */ + + /* create appropriate transition matrices (trellis) */ + d_out0.reserve(d_states); + d_out1.reserve(d_states); + d_state0.reserve(d_states); + d_state1.reserve(d_states); + + nsc_transit(d_out0.data(), d_state0.data(), 0, g_encoder, d_KK, d_nn); + nsc_transit(d_out1.data(), d_state1.data(), 1, g_encoder, d_KK, d_nn); + + // initialise trellis state + d_trellis_state_is_initialised = false; + Viterbi_Decoder_Sbas::init_trellis_state(); +} + + +void Viterbi_Decoder_Sbas::reset() +{ + init_trellis_state(); +} + + +/* Function decode_block() + Description: Uses the Viterbi algorithm to perform hard-decision decoding of a convolutional code. + Input parameters: + r[] The received signal in LLR-form. For BPSK, must be in form r = 2*a*y/(sigma^2). + LL The number of data bits to be decoded (doesn't include the mm zero-tail-bits) + Output parameters: + output_u_int[] Hard decisions on the data bits (without the mm zero-tail-bits) + */ +float Viterbi_Decoder_Sbas::decode_block(const double input_c[], int output_u_int[], const int LL) +{ + VLOG(FLOW) << "decode_block(): LL=" << LL; + + // init + init_trellis_state(); + // do add compare select + do_acs(input_c, LL + d_mm); + // tail, no need to output -> traceback, but don't decode + const int state = do_traceback(d_mm); + // traceback and decode + const int decoding_length_mismatch = do_tb_and_decode(d_mm, LL, state, output_u_int, d_indicator_metric); + + VLOG(FLOW) << "decoding length mismatch: " << decoding_length_mismatch; + + return d_indicator_metric; +} + + +float Viterbi_Decoder_Sbas::decode_continuous(const double sym[], + const int traceback_depth, + int bits[], + const int nbits_requested, + int& nbits_decoded) +{ + VLOG(FLOW) << "decode_continuous(): nbits_requested=" << nbits_requested; + + // do add compare select + do_acs(sym, nbits_requested); + // the ML sequence in the newest part of the trellis can not be decoded + // since it depends on the future values -> traceback, but don't decode + const int state = do_traceback(traceback_depth); + // traceback and decode + const int decoding_length_mismatch = do_tb_and_decode(traceback_depth, nbits_requested, state, bits, d_indicator_metric); + nbits_decoded = nbits_requested + decoding_length_mismatch; + + VLOG(FLOW) << "decoding length mismatch (continuous decoding): " << decoding_length_mismatch; + + return d_indicator_metric; +} + + +void Viterbi_Decoder_Sbas::init_trellis_state() +{ + int state; + // if trellis state has been initialised, free old state memory + if (d_trellis_state_is_initialised) + { + // init trellis state + d_pm_t.clear(); + d_rec_array.clear(); + d_metric_c.clear(); + } + + // reserve new trellis state memory + d_pm_t.reserve(d_states); + d_trellis_paths = std::deque(); + d_rec_array.reserve(d_nn); + d_metric_c.reserve(d_number_symbols); + d_trellis_state_is_initialised = true; + + /* initialize trellis */ + for (state = 0; state < d_states; state++) + { + d_pm_t[state] = -MAXLOG; + // d_pm_t_next[state] = -MAXLOG; + } + d_pm_t[0] = 0; /* start in all-zeros state */ + + d_indicator_metric = 0; +} + + +int Viterbi_Decoder_Sbas::do_acs(const double sym[], int nbits) +{ + int t; + int i; + int state_at_t; + float metric; + float max_val; + std::vector pm_t_next(d_states); + + /* t: + * - state: state at t + * - d_prev_section[state_at_t]: path metric at t for state state_at_t + * - d_out0[state_at_t]: sent symbols for a data bit 0 if state is state_at_t at time t + * + */ + + for (state_at_t = 0; state_at_t < d_states; state_at_t++) + { + pm_t_next[state_at_t] = -MAXLOG; + } + + /* go through trellis */ + for (t = 0; t < nbits; t++) + { + /* Temporarily store the received symbols current decoding step */ + for (i = 0; i < d_nn; i++) + { + d_rec_array[i] = static_cast(sym[d_nn * t + i]); + } + + /* precompute all possible branch metrics */ + for (i = 0; i < d_number_symbols; i++) + { + d_metric_c[i] = gamma(d_rec_array.data(), i, d_nn); + VLOG(LMORE) << "metric for (tx_sym=" << i << "|ry_sym=(" << d_rec_array[0] << ", " << d_rec_array[1] << ") = " << d_metric_c[i]; + } + + // find the survivor branches leading the trellis states at t+1 + Prev next_trellis_states(d_states, t + 1); + /* step through all states */ + for (state_at_t = 0; state_at_t < d_states; state_at_t++) + { + const int next_state_if_0 = d_state0[state_at_t]; + const int next_state_if_1 = d_state1[state_at_t]; + + /* hypothesis: info bit is a zero */ + const float bm_0 = d_metric_c[d_out0[state_at_t]]; + metric = d_pm_t[state_at_t] + bm_0; // path metric + zerobranch metric + + /* store new metric if more than metric in storage */ + if (metric > pm_t_next[next_state_if_0]) + { + pm_t_next[next_state_if_0] = metric; + next_trellis_states.set_current_state_as_ancestor_of_next_state(next_state_if_0, state_at_t); + next_trellis_states.set_decoded_bit_for_next_state(next_state_if_0, 0); + next_trellis_states.set_survivor_branch_metric_of_next_state(next_state_if_0, bm_0); + } + + /* hypothesis: info bit is a one */ + const float bm_1 = d_metric_c[d_out1[state_at_t]]; + metric = d_pm_t[state_at_t] + bm_1; // path metric + onebranch metric + + /* store new metric if more than metric in storage */ + if (metric > pm_t_next[next_state_if_1]) + { + pm_t_next[next_state_if_1] = metric; + next_trellis_states.set_current_state_as_ancestor_of_next_state(next_state_if_1, state_at_t); + next_trellis_states.set_decoded_bit_for_next_state(next_state_if_1, 1); + next_trellis_states.set_survivor_branch_metric_of_next_state(next_state_if_1, bm_1); + } + } + + d_trellis_paths.push_front(next_trellis_states); + + /* normalize -> afterwards, the largest metric value is always 0 */ + // max_val = 0; + max_val = -MAXLOG; + for (state_at_t = 0; state_at_t < d_states; state_at_t++) + { + if (pm_t_next[state_at_t] > max_val) + { + max_val = pm_t_next[state_at_t]; + } + } + VLOG(LMORE) << "max_val at t=" << t << ": " << max_val; + for (state_at_t = 0; state_at_t < d_states; state_at_t++) + { + d_pm_t[state_at_t] = pm_t_next[state_at_t] - max_val; + pm_t_next[state_at_t] = -MAXLOG; + } + } + + return t; +} + + +int Viterbi_Decoder_Sbas::do_traceback(size_t traceback_length) +{ + // traceback_length is in bits + int state; + std::deque::iterator it; + + VLOG(FLOW) << "do_traceback(): traceback_length=" << traceback_length << '\n'; + + if (d_trellis_paths.size() < traceback_length) + { + traceback_length = d_trellis_paths.size(); + } + + state = 0; // maybe start not at state 0, but at state with best metric + for (it = d_trellis_paths.begin(); it < d_trellis_paths.begin() + traceback_length; ++it) + { + state = it->get_anchestor_state_of_current_state(state); + } + return state; +} + + +int Viterbi_Decoder_Sbas::do_tb_and_decode(int traceback_length, int requested_decoding_length, int state, int output_u_int[], float& indicator_metric) +{ + int n_of_branches_for_indicator_metric = 500; + std::deque::iterator it; + int n_im = 0; + + VLOG(FLOW) << "do_tb_and_decode(): requested_decoding_length=" << requested_decoding_length; + // decode only decode_length bits -> overstep newer bits which are too much + const int decoding_length_mismatch = static_cast(d_trellis_paths.size()) - (traceback_length + requested_decoding_length); + VLOG(BLOCK) << "decoding_length_mismatch=" << decoding_length_mismatch; + const int overstep_length = decoding_length_mismatch >= 0 ? decoding_length_mismatch : 0; + VLOG(BLOCK) << "overstep_length=" << overstep_length; + + for (it = d_trellis_paths.begin() + traceback_length; + it < d_trellis_paths.begin() + traceback_length + overstep_length; ++it) + { + state = it->get_anchestor_state_of_current_state(state); + } + int t_out = static_cast(d_trellis_paths.end() - (d_trellis_paths.begin() + traceback_length + overstep_length) - 1); // requested_decoding_length-1; + indicator_metric = 0; + for (it = d_trellis_paths.begin() + traceback_length + overstep_length; it < d_trellis_paths.end(); ++it) + { + if (it - (d_trellis_paths.begin() + traceback_length + overstep_length) < n_of_branches_for_indicator_metric) + { + n_im++; + indicator_metric += it->get_metric_of_current_state(state); + VLOG(SAMPLE) << "@t=" << it->get_t() << " b=" << it->get_bit_of_current_state(state) << " sm=" << indicator_metric << " d=" << it->get_metric_of_current_state(state); + } + output_u_int[t_out] = it->get_bit_of_current_state(state); + state = it->get_anchestor_state_of_current_state(state); + t_out--; + } + if (n_im > 0) + { + indicator_metric /= static_cast(n_im); + } + + VLOG(BLOCK) << "indicator metric: " << indicator_metric; + // remove old states + if (d_trellis_paths.begin() + traceback_length + overstep_length <= d_trellis_paths.end()) + { + d_trellis_paths.erase(d_trellis_paths.begin() + traceback_length + overstep_length, d_trellis_paths.end()); + } + return decoding_length_mismatch; +} + + +/* function Gamma() + + Description: Computes the branch metric used for decoding. + + Output parameters: + (returned float) The metric between the hypothetical symbol and the recevieved vector + + Input parameters: + rec_array The received vector, of length nn + symbol The hypothetical symbol + nn The length of the received vector + + This function is used by siso() */ +float Viterbi_Decoder_Sbas::gamma(const float rec_array[], int symbol, int nn) +{ + float rm = 0; + int i; + unsigned int mask = 1U; + float txsym; + + for (i = 0; i < nn; i++) + { + // if (symbol & mask) rm += rec_array[nn - i - 1]; + txsym = symbol & mask ? 1 : -1; + rm += txsym * rec_array[nn - i - 1]; + mask = mask << 1U; + } + // rm = rm > 50 ? rm : -1000; + + return rm; +} + + +/* function that creates the transit and output vectors */ +void Viterbi_Decoder_Sbas::nsc_transit(int output_p[], int trans_p[], int input, const int g[], + int KK, int nn) +{ + int nextstate[1]; + int state; + const int states = static_cast(1U << (KK - 1)); /* The number of states: 2^mm */ + + /* Determine the output and next state for each possible starting state */ + for (state = 0; state < states; state++) + { + output_p[state] = nsc_enc_bit(nextstate, input, state, g, KK, nn); + trans_p[state] = nextstate[0]; + } +} + + +/* Function nsc_enc_bit() + + Description: Convolutionally encodes a single bit using a rate 1/n encoder. + Takes in one input bit at a time, and produces a n-bit output. + + Input parameters: + input The input data bit (i.e. a 0 or 1). + state_in The starting state of the encoder (an int from 0 to 2^m-1). + g[] An n-element vector containing the code generators in binary form. + KK The constraint length of the convolutional code. + nn number of symbols bits per input bits (rate 1/nn) + + Output parameters: + output_p[] An n-element vector containing the encoded bits. + state_out_p[] An integer containing the final state of the encoder + (i.e. the state after encoding this bit) + + This function is used by rsc_encode(), nsc_transit(), rsc_transit(), and nsc_transit() */ +int Viterbi_Decoder_Sbas::nsc_enc_bit(int state_out_p[], int input, int state_in, + const int g[], int KK, int nn) +{ + /* declare variables */ + int state; + int i; + int out = 0; + + /* create a word made up of state and new input */ + state = (input << (KK - 1)) ^ state_in; + + /* AND the word with the generators */ + for (i = 0; i < nn; i++) + { + /* update output symbol */ + out = (out << 1) + parity_counter(state & g[i], KK); + } + + /* shift the state to make the new state */ + state_out_p[0] = state >> 1; + return (out); +} + + +/* function parity_counter() + + Description: Determines if a symbol has odd (1) or even (0) parity + + Output parameters: + (returned int): The symbol's parity = 1 for odd and 0 for even + + Input parameters: + symbol: The integer-valued symbol + length: The highest bit position in the symbol + + This function is used by nsc_enc_bit(), rsc_enc_bit(), and rsc_tail() */ +int Viterbi_Decoder_Sbas::parity_counter(int symbol, int length) +{ + int counter; + unsigned int temp_parity = 0; + for (counter = 0; counter < length; counter++) + { + temp_parity = temp_parity ^ (symbol & 1U); + symbol = symbol >> 1U; + } + return static_cast(temp_parity); +} + + +// prev helper class +Viterbi_Decoder_Sbas::Prev::Prev(int states, int t) +{ + this->t = t; + num_states = states; + state.reserve(num_states); + v_bit.reserve(num_states); + v_metric.reserve(num_states); + refcount = 1; + std::fill_n(state.begin(), num_states, 0); + std::fill_n(v_bit.begin(), num_states, 0); + std::fill_n(v_metric.begin(), num_states, 0.0F); +} + + +// copy constructor +Viterbi_Decoder_Sbas::Prev::Prev(const Prev& prev) +{ + refcount = prev.refcount; + refcount++; + t = prev.t; + state = prev.state; + num_states = prev.num_states; + v_bit = prev.v_bit; + v_metric = prev.v_metric; + VLOG(LMORE) << "Prev(" + << "?" + << ", " << t << ")" + << " copy, new refcount = " << refcount; +} + + +// assignment constructor +Viterbi_Decoder_Sbas::Prev& Viterbi_Decoder_Sbas::Prev::operator=(const Prev& other) +{ + // check for self-assignment + if (&other == this) + { + return *this; + } + + // handle old resources + if (refcount != 1) + { // this object is not anymore using them + refcount--; + } + + // increase ref counter for this resource set + refcount = other.refcount; + refcount++; + + // take over resources + t = other.t; + state = other.state; + v_bit = other.v_bit; + v_metric = other.v_metric; + + VLOG(LMORE) << "Prev(" + << "?" + << ", " << t << ")" + << " assignment, new refcount = " << refcount; + return *this; +} + + +Viterbi_Decoder_Sbas::Prev::~Prev() +{ + if (refcount != 1) + { + refcount--; + VLOG(LMORE) << "~Prev(" + << "?" + << ", " << t << ")" + << " destructor after copy, new refcount = " << refcount; + } +} + + +int Viterbi_Decoder_Sbas::Prev::get_anchestor_state_of_current_state(int current_state) const +{ + // std::cout << "get prev state: for state " << current_state << " at time " << t << ", the prev state at time " << t - 1 << " is " << state[current_state] << '\n'; + if (num_states > current_state) + { + return state[current_state]; + } + // std::cout << "alarm " << "num_states=" << num_states << " current_state=" << current_state << '\n'; + // return state[current_state]; + return 0; +} + + +int Viterbi_Decoder_Sbas::Prev::get_bit_of_current_state(int current_state) const +{ + // std::cout << "get prev bit : for state " << current_state << " at time " << t << ", the send bit is " << bit[current_state] << '\n'; + if (num_states > current_state) + { + return v_bit[current_state]; + } + return 0; +} + + +float Viterbi_Decoder_Sbas::Prev::get_metric_of_current_state(int current_state) const +{ + if (num_states > current_state) + { + return v_metric[current_state]; + } + return 0; +} + + +int Viterbi_Decoder_Sbas::Prev::get_t() const +{ + return t; +} + + +void Viterbi_Decoder_Sbas::Prev::set_current_state_as_ancestor_of_next_state(int next_state, int current_state) +{ + if (num_states > next_state) + { + state[next_state] = current_state; + } +} + + +void Viterbi_Decoder_Sbas::Prev::set_decoded_bit_for_next_state(int next_state, int bit) +{ + if (num_states > next_state) + { + this->v_bit[next_state] = bit; + } +} + + +void Viterbi_Decoder_Sbas::Prev::set_survivor_branch_metric_of_next_state(int next_state, float metric) +{ + if (num_states > next_state) + { + this->v_metric[next_state] = metric; + } +} diff --git a/src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.h b/src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.h new file mode 100644 index 000000000..40a56c752 --- /dev/null +++ b/src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.h @@ -0,0 +1,123 @@ +/*! + * \file viterbi_decoder_sbas.h + * \brief Interface of a Viterbi decoder class based on the Iterative Solutions + * Coded Modulation Library by Matthew C. Valenti + * \author Daniel Fehr 2013. daniel.co(at)bluewin.ch + * + * ----------------------------------------------------------------------------- + * + * GNSS-SDR is a Global Navigation Satellite System software-defined receiver. + * This file is part of GNSS-SDR. + * + * Copyright (C) 2010-2020 (see AUTHORS file for a list of contributors) + * SPDX-License-Identifier: GPL-3.0-or-later + * + * ----------------------------------------------------------------------------- + */ + +#ifndef GNSS_SDR_VITERBI_DECODER_SBAS_H +#define GNSS_SDR_VITERBI_DECODER_SBAS_H + +#include // for size_t +#include +#include + +/** \addtogroup Telemetry_Decoder + * \{ */ +/** \addtogroup Telemetry_Decoder_libs telemetry_decoder_libs + * \{ */ + + +/*! + * \brief Class that implements a Viterbi decoder + */ +class Viterbi_Decoder_Sbas +{ +public: + Viterbi_Decoder_Sbas(const int g_encoder[], const int KK, const int nn); + ~Viterbi_Decoder_Sbas() = default; + void reset(); + + /*! + * \brief Uses the Viterbi algorithm to perform hard-decision decoding of a convolutional code. + * + * \param[in] input_c[] The received signal in LLR-form. For BPSK, must be in form r = 2*a*y/(sigma^2). + * \param[in] LL The number of data bits to be decoded (does not include the mm zero-tail-bits) + * + * \return output_u_int[] Hard decisions on the data bits (without the mm zero-tail-bits) + */ + float decode_block(const double input_c[], int* output_u_int, const int LL); + + float decode_continuous(const double sym[], const int traceback_depth, int bits[], + const int nbits_requested, int& nbits_decoded); + +private: + class Prev + { + public: + int num_states; + Prev(int states, int t); + Prev(const Prev& prev); + Prev& operator=(const Prev& other); + ~Prev(); + + int get_anchestor_state_of_current_state(int current_state) const; + int get_bit_of_current_state(int current_state) const; + float get_metric_of_current_state(int current_state) const; + int get_t() const; + void set_current_state_as_ancestor_of_next_state(int next_state, int current_state); + void set_decoded_bit_for_next_state(int next_state, int bit); + void set_survivor_branch_metric_of_next_state(int next_state, float metric); + + private: + std::vector v_metric; + std::vector state; + std::vector v_bit; + int t; + int refcount; + }; + + // operations on the trellis (change decoder state) + void init_trellis_state(); + int do_acs(const double sym[], int nbits); + int do_traceback(std::size_t traceback_length); + int do_tb_and_decode(int traceback_length, int requested_decoding_length, int state, int output_u_int[], float& indicator_metric); + + // branch metric function + float gamma(const float rec_array[], int symbol, int nn); + + // trellis generation + void nsc_transit(int output_p[], int trans_p[], int input, const int g[], int KK, int nn); + int nsc_enc_bit(int state_out_p[], int input, int state_in, const int g[], int KK, int nn); + int parity_counter(int symbol, int length); + + // trellis state + std::deque d_trellis_paths; + std::vector d_pm_t; + std::vector d_metric_c; /* Set of all possible branch metrics */ + std::vector d_rec_array; /* Received values for one trellis section */ + + // trellis definition + std::vector d_out0; + std::vector d_state0; + std::vector d_out1; + std::vector d_state1; + + // measures + float d_indicator_metric; + + // code properties + int d_KK; + int d_nn; + + // derived code properties + int d_mm; + int d_states; + int d_number_symbols; + bool d_trellis_state_is_initialised; +}; + + +/** \} */ +/** \} */ +#endif // GNSS_SDR_VITERBI_DECODER_SBAS_H diff --git a/src/tests/unit-tests/signal-processing-blocks/telemetry_decoder/galileo_fnav_inav_decoder_test.cc b/src/tests/unit-tests/signal-processing-blocks/telemetry_decoder/galileo_fnav_inav_decoder_test.cc index 9c3f05e23..9166a6b59 100644 --- a/src/tests/unit-tests/signal-processing-blocks/telemetry_decoder/galileo_fnav_inav_decoder_test.cc +++ b/src/tests/unit-tests/signal-processing-blocks/telemetry_decoder/galileo_fnav_inav_decoder_test.cc @@ -16,13 +16,15 @@ * ----------------------------------------------------------------------------- */ -#include "convolutional.h" #include "galileo_fnav_message.h" #include "galileo_inav_message.h" +#include "viterbi_decoder.h" #include -#include +#include // for copy +#include #include #include +#include // for std::back_inserter #include #include @@ -30,21 +32,39 @@ class Galileo_FNAV_INAV_test : public ::testing::Test { public: + Galileo_FNAV_INAV_test() + { + // vars for Viterbi decoder + viterbi_inav->nsc_transit(i_out0, i_state0, 0); + viterbi_inav->nsc_transit(i_out1, i_state1, 1); + viterbi_fnav->nsc_transit(f_out0, f_state0, 0); + viterbi_fnav->nsc_transit(f_out1, f_state1, 1); + flag_even_word_arrived = 0; + } + + ~Galileo_FNAV_INAV_test() = default; + Galileo_Inav_Message INAV_decoder; Galileo_Fnav_Message FNAV_decoder; - // vars for Viterbi decoder - int32_t *out0, *out1, *state0, *state1; - int32_t g_encoder[2]; const int32_t nn = 2; // Coding rate 1/n const int32_t KK = 7; // Constraint Length int32_t mm = KK - 1; int32_t flag_even_word_arrived; - void viterbi_decoder(float *page_part_symbols, int32_t *page_part_bits, int32_t _datalength) - { - Viterbi(page_part_bits, out0, state0, out1, state1, - page_part_symbols, KK, nn, _datalength); - } + const std::array g_encoder{{121, 91}}; + std::shared_ptr viterbi_fnav = std::make_shared(KK, nn, ((488 / nn) - mm), g_encoder); + std::shared_ptr viterbi_inav = std::make_shared(KK, nn, ((240 / nn) - mm), g_encoder); + int32_t max_states = 1 << mm; // 2^mm + + std::vector i_out0 = std::vector(max_states); + std::vector i_out1 = std::vector(max_states); + std::vector i_state0 = std::vector(max_states); + std::vector i_state1 = std::vector(max_states); + + std::vector f_out0 = std::vector(max_states); + std::vector f_out1 = std::vector(max_states); + std::vector f_state0 = std::vector(max_states); + std::vector f_state1 = std::vector(max_states); void deinterleaver(int32_t rows, int32_t cols, const float *in, float *out) { @@ -57,12 +77,12 @@ public: } } - bool decode_INAV_word(float *page_part_symbols, int32_t frame_length) { // 1. De-interleave - auto *page_part_symbols_deint = static_cast(volk_gnsssdr_malloc(frame_length * sizeof(float), volk_gnsssdr_get_alignment())); - deinterleaver(GALILEO_INAV_INTERLEAVER_ROWS, GALILEO_INAV_INTERLEAVER_COLS, page_part_symbols, page_part_symbols_deint); + std::vector page_part_symbols_deint = std::vector(frame_length / 2); + std::copy(&page_part_symbols[0], &page_part_symbols[frame_length / 2], std::back_inserter(page_part_symbols_deint)); + deinterleaver(GALILEO_INAV_INTERLEAVER_ROWS, GALILEO_INAV_INTERLEAVER_COLS, page_part_symbols, page_part_symbols_deint.data()); // 2. Viterbi decoder // 2.1 Take into account the NOT gate in G2 polynomial (Galileo ICD Figure 13, FEC encoder) @@ -75,12 +95,8 @@ public: } } - auto *page_part_bits = static_cast(volk_gnsssdr_malloc((frame_length / 2) * sizeof(int32_t), volk_gnsssdr_get_alignment())); - - const int32_t CodeLength = 240; - int32_t DataLength = (CodeLength / nn) - mm; - viterbi_decoder(page_part_symbols_deint, page_part_bits, DataLength); - volk_gnsssdr_free(page_part_symbols_deint); + std::vector page_part_bits = std::vector(frame_length / 2); + viterbi_inav->decode(page_part_bits, i_out0, i_state0, i_out1, i_state1, page_part_symbols_deint); // 3. Call the Galileo page decoder std::string page_String; @@ -103,8 +119,6 @@ public: INAV_decoder.split_page(page_String, flag_even_word_arrived); if (INAV_decoder.get_flag_CRC_test() == true) { - std::cout << "Galileo E1 INAV PAGE CRC correct \n"; - // std::cout << "Galileo E1 CRC correct on channel " << d_channel << " from satellite " << d_satellite << '\n'; crc_ok = true; } flag_even_word_arrived = 0; @@ -115,15 +129,15 @@ public: INAV_decoder.split_page(page_String.c_str(), flag_even_word_arrived); flag_even_word_arrived = 1; } - volk_gnsssdr_free(page_part_bits); return crc_ok; } bool decode_FNAV_word(float *page_symbols, int32_t frame_length) { // 1. De-interleave - auto *page_symbols_deint = static_cast(volk_gnsssdr_malloc(frame_length * sizeof(float), volk_gnsssdr_get_alignment())); - deinterleaver(GALILEO_FNAV_INTERLEAVER_ROWS, GALILEO_FNAV_INTERLEAVER_COLS, page_symbols, page_symbols_deint); + std::vector page_symbols_deint = std::vector(frame_length); + std::copy(&page_symbols[0], &page_symbols[frame_length / 2], std::back_inserter(page_symbols_deint)); + deinterleaver(GALILEO_FNAV_INTERLEAVER_ROWS, GALILEO_FNAV_INTERLEAVER_COLS, page_symbols, page_symbols_deint.data()); // 2. Viterbi decoder // 2.1 Take into account the NOT gate in G2 polynomial (Galileo ICD Figure 13, FEC encoder) @@ -135,13 +149,9 @@ public: page_symbols_deint[i] = -page_symbols_deint[i]; } } - auto *page_bits = static_cast(volk_gnsssdr_malloc((frame_length / 2) * sizeof(int32_t), volk_gnsssdr_get_alignment())); - const int32_t CodeLength = 488; - int32_t DataLength = (CodeLength / nn) - mm; - viterbi_decoder(page_symbols_deint, page_bits, DataLength); - - volk_gnsssdr_free(page_symbols_deint); + std::vector page_bits = std::vector(frame_length); + viterbi_fnav->decode(page_bits, f_out0, f_state0, f_out1, f_state1, page_symbols_deint); // 3. Call the Galileo page decoder std::string page_String; @@ -156,43 +166,18 @@ public: page_String.push_back('0'); } } - volk_gnsssdr_free(page_bits); // DECODE COMPLETE WORD (even + odd) and TEST CRC FNAV_decoder.split_page(page_String); if (FNAV_decoder.get_flag_CRC_test() == true) { - std::cout << "Galileo E5a FNAV PAGE CRC correct \n"; return true; } return false; } - - Galileo_FNAV_INAV_test() - { - // vars for Viterbi decoder - int32_t max_states = 1 << mm; // 2^mm - g_encoder[0] = 121; // Polynomial G1 - g_encoder[1] = 91; // Polynomial G2 - out0 = static_cast(volk_gnsssdr_malloc(max_states * sizeof(int32_t), volk_gnsssdr_get_alignment())); - out1 = static_cast(volk_gnsssdr_malloc(max_states * sizeof(int32_t), volk_gnsssdr_get_alignment())); - state0 = static_cast(volk_gnsssdr_malloc(max_states * sizeof(int32_t), volk_gnsssdr_get_alignment())); - state1 = static_cast(volk_gnsssdr_malloc(max_states * sizeof(int32_t), volk_gnsssdr_get_alignment())); - // create appropriate transition matrices - nsc_transit(out0, state0, 0, g_encoder, KK, nn); - nsc_transit(out1, state1, 1, g_encoder, KK, nn); - flag_even_word_arrived = 0; - } - - ~Galileo_FNAV_INAV_test() - { - volk_gnsssdr_free(out0); - volk_gnsssdr_free(out1); - volk_gnsssdr_free(state0); - volk_gnsssdr_free(state1); - } }; + TEST_F(Galileo_FNAV_INAV_test, ValidationOfResults) { std::chrono::time_point start, end; @@ -266,5 +251,5 @@ TEST_F(Galileo_FNAV_INAV_test, ValidationOfResults) }) << "Exception during INAV frame decoding"; end = std::chrono::system_clock::now(); elapsed_seconds = end - start; - std::cout << "Galileo FNAV/INAV Test completed in " << elapsed_seconds.count() * 1e6 << " microseconds\n"; + std::cout << "Galileo INAV/FNAV CRC and Viterbi decoder test completed in " << elapsed_seconds.count() * 1e6 << " microseconds\n"; }