mirrors
/
gnss-sdr
mirror of https://github.com/gnss-sdr/gnss-sdr
1
0
Fork 0
Code Issues Releases Wiki Activity

Rewriting of the Viterbi decoding 

New class implementing the Viterbi decoder. All memory allocated in the constructor

Old SBAS Viterbi decoder moved to Viterbi_Decoder_Sbas class.
This commit is contained in:
Carles Fernandez 2021-09-25 20:29:51 +02:00
parent 010069b86a
commit e52ccfb893
No known key found for this signature in database
GPG Key ID: 4C583C52B0C3877D
12 changed files with 1022 additions and 993 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
docs/CHANGELOG.md
View File

@ -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

114
src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.cc
View File

@ -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 <glog/logging.h>
#include <gnuradio/io_signature.h>
#include <pmt/pmt.h> // for make_any
#include <pmt/pmt_sugar.h> // for mp
#include <array>
#include <cmath> // for fmod
#include <cstddef> // for size_t
#include <cstdlib> // for abs
#include <exception> // for exception
#include <iostream> // for cout
#include <memory> // 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 <glog/logging.h> // for LOG, DLOG
#include <gnuradio/io_signature.h> // for gr::io_signature::make
#include <pmt/pmt.h> // for pmt::make_any
#include <pmt/pmt_sugar.h> // for pmt::mp
#include <array> // for std::array
#include <cmath> // for std::fmod, std::abs
#include <cstddef> // for size_t
#include <exception> // for std::exception
#include <iostream> // 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<int32_t, 2> 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<int32_t>(d_samples_per_preamble);
d_frame_length_symbols = GALILEO_INAV_PAGE_PART_SYMBOLS - GALILEO_INAV_PREAMBLE_LENGTH_BITS;
d_codelength = static_cast<int32_t>(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<int32_t>(d_samples_per_preamble);
d_frame_length_symbols = GALILEO_FNAV_SYMBOLS_PER_PAGE - GALILEO_FNAV_PREAMBLE_LENGTH_BITS;
d_codelength = static_cast<int32_t>(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<int32_t>(d_samples_per_preamble);
d_frame_length_symbols = GALILEO_CNAV_SYMBOLS_PER_PAGE - GALILEO_CNAV_PREAMBLE_LENGTH_BITS;
d_codelength = static_cast<int32_t>(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<uint32_t>(d_mm); // 2^d_mm
std::array<int32_t, 2> g_encoder{{121, 91}}; // Polynomial G1 and G2
d_out0 = std::vector<int32_t>(max_states);
d_out1 = std::vector<int32_t>(max_states);
d_state0 = std::vector<int32_t>(max_states);
d_state1 = std::vector<int32_t>(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<Viterbi_Decoder>(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<float> 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<float> 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<int32_t> 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<float> page_symbols_deint(frame_length);
deinterleaver(GALILEO_FNAV_INTERLEAVER_ROWS, GALILEO_FNAV_INTERLEAVER_COLS, page_symbols, page_symbols_deint.data());
std::vector<float> 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<int32_t> 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<float> page_symbols_deint(page_length);
deinterleaver(GALILEO_CNAV_INTERLEAVER_ROWS, GALILEO_CNAV_INTERLEAVER_COLS, page_symbols, page_symbols_deint.data());
std::vector<float> 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<int32_t> 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<uint32_t>(d_mm); // 2^d_mm
d_out0 = std::vector<int32_t>(max_states);
d_out1 = std::vector<int32_t>(max_states);
d_state0 = std::vector<int32_t>(max_states);
d_state1 = std::vector<int32_t>(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<int32_t>(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;

44
src/algorithms/telemetry_decoder/gnuradio_blocks/galileo_telemetry_decoder_gs.h
View File

@ -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 <boost/circular_buffer.hpp>
#include <gnuradio/block.h> // for block
#include <gnuradio/types.h> // for gr_vector_const_void_star
#include <cstdint>
#include <fstream>
#include <memory> // for std::unique_ptr
#include <string>
#include <vector>
#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 <boost/circular_buffer.hpp> // for boost::circular_buffer
#include <gnuradio/block.h> // for block
#include <gnuradio/types.h> // for gr_vector_const_void_star
#include <cstdint> // for int32_t, uint32_t
#include <fstream> // for std::ofstream
#include <memory> // for std::unique_ptr
#include <string> // for std::string
#include <vector> // 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<galileo_telemetry_decoder_gs>;
@ -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<Viterbi_Decoder> d_viterbi;
std::vector<int32_t> d_preamble_samples;
std::vector<float> d_page_part_symbols;
std::vector<int32_t> d_out0;

6
src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.cc
View File

@ -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 <glog/logging.h>
#include <gnuradio/io_signature.h>
#include <pmt/pmt_sugar.h> // 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<int32_t, nn> g_encoder{121, 91};
d_vd1 = std::make_shared<Viterbi_Decoder>(g_encoder.data(), d_KK, nn);
d_vd2 = std::make_shared<Viterbi_Decoder>(g_encoder.data(), d_KK, nn);
d_vd1 = std::make_shared<Viterbi_Decoder_Sbas>(g_encoder.data(), d_KK, nn);
d_vd2 = std::make_shared<Viterbi_Decoder_Sbas>(g_encoder.data(), d_KK, nn);
d_past_symbol = 0;
}

6
src/algorithms/telemetry_decoder/gnuradio_blocks/sbas_l1_telemetry_decoder_gs.h
View File

@ -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<Viterbi_Decoder> d_vd1;
std::shared_ptr<Viterbi_Decoder> d_vd2;
std::shared_ptr<Viterbi_Decoder_Sbas> d_vd1;
std::shared_ptr<Viterbi_Decoder_Sbas> d_vd2;
double d_past_symbol;
} d_symbol_aligner_and_decoder;

13
src/algorithms/telemetry_decoder/libs/CMakeLists.txt
View File

@ -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

278
src/algorithms/telemetry_decoder/libs/convolutional.h
View File

@ -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 <volk_gnsssdr/volk_gnsssdr.h>
#include <vector>
/** \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<float> prev_section(states, -MAXLOG);
std::vector<float> next_section(states, -MAXLOG);
std::vector<int> prev_bit(states * (LL + mm), 0);
std::vector<int> prev_state(states * (LL + mm), 0);
std::vector<float> rec_array(nn);
std::vector<float> 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

598
src/algorithms/telemetry_decoder/libs/viterbi_decoder.cc
View File

@ -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 <glog/logging.h>
#include <algorithm> // for fill_n
#include <ostream> // for operator<<, basic_ostream, char_traits
#include <volk_gnsssdr/volk_gnsssdr.h> // for volk_gnsssdr_32f_index_max_32u
#include <algorithm> // 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<int32_t, 2>& 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<int>(1U << d_mm); /* 2^mm */
d_number_symbols = static_cast<int>(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<float>(d_states, -d_MAXLOG);
d_next_section = std::vector<float>(d_states, -d_MAXLOG);
d_rec_array = std::vector<float>(d_nn);
d_metric_c = std::vector<float>(d_number_symbols);
d_prev_bit = std::vector<int32_t>(d_states * (d_LL + d_mm));
d_prev_state = std::vector<int32_t>(d_states * (d_LL + d_mm));
d_g = g;
}
void Viterbi_Decoder::reset()
void Viterbi_Decoder::decode(std::vector<int32_t>& output_u_int,
const std::vector<int32_t>& out0,
const std::vector<int32_t>& state0,
const std::vector<int32_t>& out1,
const std::vector<int32_t>& state1,
const std::vector<float>& 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<Prev>();
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<float> 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<float>(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<int32_t>& output_p,
std::vector<int32_t>& trans_p,
int32_t input) const
{
// traceback_length is in bits
int state;
std::deque<Prev>::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<Prev>::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<int>(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<int>(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<float>(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<int>(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<int>(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_;
}

159
src/algorithms/telemetry_decoder/libs/viterbi_decoder.h
View File

@ -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 <cstddef> // for size_t
#include <deque>
#include <array>
#include <cstdint>
#include <vector>
/** \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<int32_t, 2>& 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<int32_t>& output_u_int,
const std::vector<int32_t>& out0,
const std::vector<int32_t>& state0,
const std::vector<int32_t>& out1,
const std::vector<int32_t>& state1,
const std::vector<float>& 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<int32_t>& output_p,
std::vector<int32_t>& 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<float> v_metric;
std::vector<int> state;
std::vector<int> 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<float> d_prev_section{};
std::vector<float> d_next_section{};
// branch metric function
float gamma(const float rec_array[], int symbol, int nn);
std::vector<float> d_rec_array{};
std::vector<float> d_metric_c{};
std::vector<int32_t> d_prev_bit{};
std::vector<int32_t> d_prev_state{};
std::array<int32_t, 2> 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<Prev> d_trellis_paths;
std::vector<float> d_pm_t;
std::vector<float> d_metric_c; /* Set of all possible branch metrics */
std::vector<float> d_rec_array; /* Received values for one trellis section */
// trellis definition
std::vector<int> d_out0;
std::vector<int> d_state0;
std::vector<int> d_out1;
std::vector<int> 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

569
src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.cc Normal file
View File

@ -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 <glog/logging.h>
#include <algorithm> // for fill_n
#include <ostream> // 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<int>(1U << d_mm); /* 2^mm */
d_number_symbols = static_cast<int>(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<Prev>();
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<float> 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<float>(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<Prev>::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<Prev>::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<int>(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<int>(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<float>(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<int>(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<int>(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;
}
}

123
src/algorithms/telemetry_decoder/libs/viterbi_decoder_sbas.h Normal file
View File

@ -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 <cstddef> // for size_t
#include <deque>
#include <vector>
/** \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<float> v_metric;
std::vector<int> state;
std::vector<int> 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<Prev> d_trellis_paths;
std::vector<float> d_pm_t;
std::vector<float> d_metric_c; /* Set of all possible branch metrics */
std::vector<float> d_rec_array; /* Received values for one trellis section */
// trellis definition
std::vector<int> d_out0;
std::vector<int> d_state0;
std::vector<int> d_out1;
std::vector<int> 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

99
src/tests/unit-tests/signal-processing-blocks/telemetry_decoder/galileo_fnav_inav_decoder_test.cc
View File

@ -16,13 +16,15 @@
* -----------------------------------------------------------------------------
*/
#include "convolutional.h"
#include "galileo_fnav_message.h"
#include "galileo_inav_message.h"
#include "viterbi_decoder.h"
#include <gtest/gtest.h>
#include <volk_gnsssdr/volk_gnsssdr.h>
#include <algorithm> // for copy
#include <array>
#include <chrono>
#include <exception>
#include <iterator> // for std::back_inserter
#include <string>
#include <unistd.h>
@ -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<int32_t, 2> g_encoder{{121, 91}};
std::shared_ptr<Viterbi_Decoder> viterbi_fnav = std::make_shared<Viterbi_Decoder>(KK, nn, ((488 / nn) - mm), g_encoder);
std::shared_ptr<Viterbi_Decoder> viterbi_inav = std::make_shared<Viterbi_Decoder>(KK, nn, ((240 / nn) - mm), g_encoder);
int32_t max_states = 1 << mm; // 2^mm
std::vector<int32_t> i_out0 = std::vector<int32_t>(max_states);
std::vector<int32_t> i_out1 = std::vector<int32_t>(max_states);
std::vector<int32_t> i_state0 = std::vector<int32_t>(max_states);
std::vector<int32_t> i_state1 = std::vector<int32_t>(max_states);
std::vector<int32_t> f_out0 = std::vector<int32_t>(max_states);
std::vector<int32_t> f_out1 = std::vector<int32_t>(max_states);
std::vector<int32_t> f_state0 = std::vector<int32_t>(max_states);
std::vector<int32_t> f_state1 = std::vector<int32_t>(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<float *>(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<float> page_part_symbols_deint = std::vector<float>(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<int32_t *>(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<int32_t> page_part_bits = std::vector<int32_t>(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<float *>(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<float> page_symbols_deint = std::vector<float>(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<int32_t *>(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<int32_t> page_bits = std::vector<int32_t>(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<int32_t *>(volk_gnsssdr_malloc(max_states * sizeof(int32_t), volk_gnsssdr_get_alignment()));
out1 = static_cast<int32_t *>(volk_gnsssdr_malloc(max_states * sizeof(int32_t), volk_gnsssdr_get_alignment()));
state0 = static_cast<int32_t *>(volk_gnsssdr_malloc(max_states * sizeof(int32_t), volk_gnsssdr_get_alignment()));
state1 = static_cast<int32_t *>(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<std::chrono::system_clock> 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";
}