diff --git a/src-core/CMakeLists.txt b/src-core/CMakeLists.txt index 7faf08384..8541877fc 100644 --- a/src-core/CMakeLists.txt +++ b/src-core/CMakeLists.txt @@ -16,6 +16,7 @@ file(GLOB_RECURSE SatDump_core_CPPS *.cpp common/map/maidenhead.c libs/correct/*.c libs/openjp2/*.c + libs/deepspace-turbo/*.c ) find_package(PkgConfig) diff --git a/src-core/common/codings/turbo/ccsds_turbo.cpp b/src-core/common/codings/turbo/ccsds_turbo.cpp new file mode 100644 index 000000000..ba69e3088 --- /dev/null +++ b/src-core/common/codings/turbo/ccsds_turbo.cpp @@ -0,0 +1,200 @@ +#include "ccsds_turbo.h" +#include +#include +#include + +namespace codings +{ + namespace turbo + { + CCSDSTurbo::CCSDSTurbo(turbo_base_t base, turbo_rate_t type) + : d_base(base), d_code_type(type) + { + // d_pack = new blocks::kernel::pack_k_bits(8); + d_info_length = base * 8; + + int p[8] = {31, 37, 43, 47, 53, 59, 61, 67}; + int k1 = 8; + int k2 = base; + + d_pi = (int *)malloc(d_info_length * sizeof *d_pi); + + for (int s = 1; s <= d_info_length; ++s) + { + int m = (s - 1) % 2; + int i = (int)floor((s - 1) / (2 * k2)); + int j = (int)floor((s - 1) / 2) - i * k2; + int t = (19 * i + 1) % (k1 / 2); + int q = t % 8 + 1; + int c = (p[q - 1] * j + 21 * m) % k2; + d_pi[s - 1] = 2 * (t + c * (k1 / 2) + 1) - m - 1; + } + + int N_components_upper; + int N_components_lower; + + d_backward = "0011"; + + switch (d_code_type) + { + case RATE_1_2: + { + N_components_upper = 2; + N_components_lower = 1; + + d_forward_upper[0] = "10011"; // systematic output + d_forward_upper[1] = "11011"; + + d_forward_lower[0] = "11011"; + // need to define puncturing pattern here maybe with a pointer to function + // 110 101 110 101 110 101 + + d_code1 = convcode_initialize((char **)d_forward_upper, (char *)d_backward, N_components_upper); + d_code2 = convcode_initialize((char **)d_forward_lower, (char *)d_backward, N_components_lower); + d_turbo = turbo_initialize(d_code1, d_code2, d_pi, d_info_length); + d_rate = 1.0 / 2.0; + d_encoded_length = d_turbo.encoded_length * 2 / 3; + break; + } + case RATE_1_3: + { + N_components_upper = 2; + N_components_lower = 1; + + d_forward_upper[0] = "10011"; // systematic output + d_forward_upper[1] = "11011"; + + d_forward_lower[0] = "11011"; // no need for puncturing + + d_code1 = convcode_initialize((char **)d_forward_upper, (char *)d_backward, N_components_upper); + d_code2 = convcode_initialize((char **)d_forward_lower, (char *)d_backward, N_components_lower); + d_turbo = turbo_initialize(d_code1, d_code2, d_pi, d_info_length); + d_rate = 1.0 / 3.0; + d_encoded_length = d_turbo.encoded_length; + break; + } + case RATE_1_4: + { + N_components_upper = 3; + N_components_lower = 1; + + d_forward_upper[0] = "10011"; // systematic output + d_forward_upper[1] = "10101"; + d_forward_upper[2] = "11111"; + + d_forward_lower[0] = "11011"; // no need for puncturing + + d_code1 = convcode_initialize((char **)d_forward_upper, (char *)d_backward, N_components_upper); + d_code2 = convcode_initialize((char **)d_forward_lower, (char *)d_backward, N_components_lower); + d_turbo = turbo_initialize(d_code1, d_code2, d_pi, d_info_length); + d_rate = 1.0 / 4.0; + d_encoded_length = d_turbo.encoded_length; + break; + } + case RATE_1_6: + { + N_components_upper = 4; + N_components_lower = 2; + + d_forward_upper[0] = "10011"; // systematic output + d_forward_upper[1] = "11011"; + d_forward_upper[2] = "10101"; + d_forward_upper[3] = "11111"; + + d_forward_lower[0] = "11011"; // no need for puncturing + d_forward_lower[1] = "11111"; + + d_code1 = convcode_initialize((char **)d_forward_upper, (char *)d_backward, N_components_upper); + d_code2 = convcode_initialize((char **)d_forward_lower, (char *)d_backward, N_components_lower); + d_turbo = turbo_initialize(d_code1, d_code2, d_pi, d_info_length); + d_rate = 1.0 / 6.0; + d_encoded_length = d_turbo.encoded_length; + break; + } + } + } + + CCSDSTurbo::~CCSDSTurbo() + { + delete[] d_pi; + } + + void CCSDSTurbo::encode(uint8_t *frame, uint8_t *codeword) + { + int *bits_in = (int *)malloc(d_encoded_length * sizeof(int *)); + for (int i = 0; i < d_info_length / 8; i++) + for (int j = 0; j < 8; j++) + bits_in[i * 8 + j] = (frame[i] & (0x80 >> j)) ? 1 : 0; + + int *encoded = turbo_encode(bits_in, d_turbo); + + uint8_t *encoded_u8 = (uint8_t *)malloc(d_encoded_length * sizeof(uint8_t *)); + if (d_code_type == RATE_1_2) + { + int j = 0; + for (int i = 0; i < d_turbo.encoded_length; i++) + { + if (puncturing(i)) + { + encoded_u8[j] = encoded[i]; + j++; + } + } + } + else + { + for (int i = 0; i < d_encoded_length; i++) + encoded_u8[i] = encoded[i]; + } + + memset(codeword, 0, d_encoded_length / 8); + + for (int i = 0; i < d_encoded_length; i++) + codeword[i / 8] = codeword[i / 8] << 1 | encoded_u8[i]; + } + + void CCSDSTurbo::decode(float *codeword, uint8_t *frame, int iterations) + { + d_turbo.interleaver = d_pi; + + const float *bits_in = codeword; + + double *bits_depunctured = (double *)malloc(sizeof(double) * d_turbo.encoded_length); + if (d_code_type == RATE_1_2) + { + int j = 0; + for (int i = 0; i < d_turbo.encoded_length; i++) + { + if (puncturing(i)) + { + bits_depunctured[i] = bits_in[j]; + j++; + } + else + { + bits_depunctured[i] = 0.0; + } + } + } + else + { + for (int i = 0; i < d_encoded_length; i++) + bits_depunctured[i] = bits_in[i]; + } + + int *decoded = turbo_decode(bits_depunctured, iterations, d_sigma * d_sigma, d_turbo); + + uint8_t *decoded_u8 = frame; + for (int i = 0; i < d_info_length / 8; i++) + { + decoded_u8[i] = 0; + for (int j = 0; j < 8; j++) + decoded_u8[i] |= decoded[i * 8 + j] ? (0x80 >> j) : 0; + } + + free(bits_depunctured); + free(decoded); + // free(decoded_u8); + } + } +} \ No newline at end of file diff --git a/src-core/common/codings/turbo/ccsds_turbo.h b/src-core/common/codings/turbo/ccsds_turbo.h new file mode 100644 index 000000000..7a9f29029 --- /dev/null +++ b/src-core/common/codings/turbo/ccsds_turbo.h @@ -0,0 +1,95 @@ +#pragma once + +extern "C" +{ +#include "libs/deepspace-turbo/libconvcodes.h" +#include "libs/deepspace-turbo/libturbocodes.h" +} +#include + +namespace codings +{ + namespace turbo + { + enum turbo_base_t + { + BASE_223 = 223, + BASE_446 = 446, + BASE_892 = 892, + BASE_1115 = 1115, + }; + + enum turbo_rate_t + { + RATE_1_2, + RATE_1_3, + RATE_1_4, + RATE_1_6, + }; + + /* + CCSDS Turbo Decoder based on : + https://github.com/khawatkom/gr-ccsds-1 + */ + class CCSDSTurbo + { + private: + int d_base; + int d_octets; + turbo_rate_t d_code_type; + + float d_rate; + int d_info_length; + int d_encoded_length; + + float d_sigma = 0.707; + + static const int MAX_COMPONENTS = 4; + + int *d_pi; + const char *d_forward_upper[MAX_COMPONENTS]; + const char *d_forward_lower[MAX_COMPONENTS]; + const char *d_backward; + t_convcode d_code1; + t_convcode d_code2; + t_turbocode d_turbo; + + int puncturing(int k) + { + + int bit_idx = k % 3; + + // bit 0,3,6,... corresponding to systematic output + if (!bit_idx) + return 1; + + // get block index + int block_idx = k / 3; + + // on odd blocks puncture second bit + if (block_idx % 2) + return bit_idx != 1; + + // on even blocks puncture third bit + return bit_idx != 2; + } + + public: + CCSDSTurbo(turbo_base_t base, turbo_rate_t type); + ~CCSDSTurbo(); + + // Get specifics of the current code + int frame_length() { return d_info_length; } + int codeword_length() { return d_encoded_length; } + + // Set Sigma for decoding + float set_sigma(float sigma) { d_sigma = sigma; } + + // Encode a Turbo codeword, takes bytes in, output bytes + void encode(uint8_t *frame, uint8_t *codeword); + + // Decode a Turbo codeword, takes soft-bits floats in, outputs bytes + void decode(float *codeword, uint8_t *frame, int iterations = 10); + }; + } +} \ No newline at end of file diff --git a/src-core/libs/deepspace-turbo/LICENSE.md b/src-core/libs/deepspace-turbo/LICENSE.md new file mode 100644 index 000000000..71267c88e --- /dev/null +++ b/src-core/libs/deepspace-turbo/LICENSE.md @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2017 Gianluca Marcon + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/src-core/libs/deepspace-turbo/libconvcodes.c b/src-core/libs/deepspace-turbo/libconvcodes.c new file mode 100644 index 000000000..21df3e833 --- /dev/null +++ b/src-core/libs/deepspace-turbo/libconvcodes.c @@ -0,0 +1,549 @@ +// +// Created by gianluca on 20/02/17. +// +#include +#include +#include +#include +#include "libconvcodes.h" + +static int get_bit(int num, int position) +{ + return (num >> position) & 1; +} + +static char* state2str(int state, int memory) +{ + char *str_state = malloc(memory + 1);/*{{{*/ + str_state[memory] = '\0'; + + for (int i = 0; i < memory; i++) { + str_state[i] = '0' + get_bit(state, memory - 1 - i); + } + + return str_state;/*}}}*/ +} + +static int convcode_stateupdate(int state, int input, t_convcode code) +{ + int memory = code.memory;/*{{{*/ + + int first_reg = 0; + for (int i = 0; i < memory; i++) + first_reg = (first_reg + code.backward_connections[i]*get_bit(state, memory - 1 - i)) % 2; + + // shift the content of the registers + int new_state = state >> 1; + + // compute the new content of the first register (MSB) + first_reg = (first_reg + input) % 2; + + // switch last bit + new_state ^= (-first_reg ^ new_state) & (1 << (memory - 1)); + + return new_state;/*}}}*/ +} + +static int *convcode_output(int state, int input, t_convcode code) +{ + int *output = calloc(code.components, sizeof(int));/*{{{*/ + int new_state = convcode_stateupdate(state, input, code); + + // get content of first register of the new state + // we have to add it to the feedforward part + int first_reg = get_bit(new_state, code.memory - 1); + + for (int c = 0; c < code.components; c++) { + output[c] = code.forward_connections[c][0]*first_reg; + + for (int i = 0; i < code.memory; i++) + output[c] = (output[c] + code.forward_connections[c][i+1]*get_bit(state, code.memory - 1 - i)) % 2; + } + + return output;/*}}}*/ +} + +t_convcode convcode_initialize(char *forward[], char *backward, int N_components) +{ +/*{{{*/ + // code initialized + t_convcode code; + + code.components = N_components; + + // number of shift registers + int code_memory = strlen(backward); + code.memory = code_memory; + + // initialize connection arrays + int **fwd_con = malloc(N_components * sizeof(int*)); + int *bwd_con = malloc(code_memory*sizeof(int)); + + // convert input strings to arrays + for (int i = 0; i < N_components; i++) { + fwd_con[i] = malloc((code_memory+1) * sizeof(int)); + + int j = 0; + for (; j < code_memory; j++) { + fwd_con[i][j] = forward[i][j] - '0'; + bwd_con[j] = backward[j] - '0'; + } + fwd_con[i][j] = forward[i][j] - '0'; + } + + code.forward_connections = fwd_con; + code.backward_connections = bwd_con; + + int N_states = 2 << (code_memory - 1); + int **neighbors = malloc(N_states * sizeof(int*)); + + // populate lookup table for state-update function + // and create neighbors array + int **next_state = malloc(N_states * sizeof(int*)); + for (int i = 0; i < N_states; i++) { + // initialize to 0 + neighbors[i] = calloc(2, sizeof(int)); + } + + for (int i = 0; i < N_states; i++) { + next_state[i] = malloc(2 * sizeof(int)); + + int updated0 = convcode_stateupdate(i, 0, code); + next_state[i][0] = updated0; + + + // save to neighbords array, use minus sign if input is 0 + // plus sign if input is 1. check whether it's possible to + // write by checking if it's content is zero. + // Exploit the fact that in binary codes a state only + // has two neighbors + if (!neighbors[updated0][0]) + neighbors[updated0][0] = -(i + 1); + else + neighbors[updated0][1] = -(i + 1); + + int updated1 = convcode_stateupdate(i, 1, code); + next_state[i][1] = updated1; + + if (!neighbors[updated1][0]) + neighbors[updated1][0] = i + 1; + else + neighbors[updated1][1] = i + 1; + } + + code.next_state = next_state; + code.neighbors = neighbors; + + + // populate output function lookup table + int ***output; + output = malloc(N_states * sizeof(int**)); + for (int i = 0; i < N_states; i++) { + output[i] = malloc(2*sizeof(int*)); + for (int j = 0; j < 2; j++) + output[i][j] = convcode_output(i, j, code); + } + code.output = output; + + return code;/*}}}*/ +} + +void convcode_clear(t_convcode code) +{ + for (int i = 0; i < code.components; i++) {/*{{{*/ + /* printf("Component %d \t Address %p\n", i, code.forward_connections[i]); */ + free(code.forward_connections[i]); + free(code.next_state[i]); + free(code.neighbors[i]); + } + free(code.output); + free(code.forward_connections); + free(code.backward_connections); + free(code.next_state); + free(code.neighbors);/*}}}*/ +} +/*void convcode_clear(t_convcode *code) +{ + for (int i = 0; i < code->components; i++) + // printf("Component %d \t Address %p\n", i, code.forward_connections[i]); + free(code->forward_connections[i]); + free(code->next_state[i]); + free(code->neighbors[i]); + + + for (int j = 0; j < 2; ++j) { + free(code->output[i][j]); + } + } + + free(code->output); + free(code->forward_connections); + free(code->backward_connections); + free(code->next_state); + free(code->neighbors); +}*/ +int* convcode_encode(int *packet, int packet_length, t_convcode code) +{ + // add support for puncturing patterns?/*{{{*/ + int encoded_length = (packet_length + code.memory) * code.components; + int *encoded_packet = malloc(encoded_length * sizeof *encoded_packet); + + int state = 0; + + for (int i = 0; i < packet_length; i++) + { + int current_bit = packet[i]; + int *output = code.output[state][current_bit]; + state = code.next_state[state][current_bit]; + + for (int c = 0; c < code.components; c++) + { + int out = output[c]; + encoded_packet[code.components * i + c] = output[c]; + } + } + + // add trellis termination + for (int i = packet_length; i < packet_length + code.memory; i++) + { + int input = 0; + + // input is equal to the feedback part in order to inject zeros into the registers + for (int j = 0; j < code.memory; j++) + input = (input + code.backward_connections[j]*get_bit(state, code.memory - 1 - j)) % 2; + + int *output = code.output[state][input]; + state = code.next_state[state][input]; + + for (int c = 0; c < code.components; c++) + encoded_packet[code.components * i + c] = output[c]; + + } + + return encoded_packet;/*}}}*/ +} + +int* convcode_decode(double *received, int length, t_convcode code) +{ + int N_states = 2 << (code.memory - 1);/*{{{*/ + int packet_length = length / code.components - code.memory; + int *decoded_packet = malloc(packet_length * sizeof *decoded_packet); + + // allocate matrix containing survivor sequences and metric vector + double *metric = malloc(N_states * sizeof *metric); + int **data_matrix; + data_matrix = malloc(N_states * sizeof(int*)); + + for (int i = 0; i < N_states; i++ ) + { + data_matrix[i] = malloc((packet_length + code.memory)* sizeof(int)); + metric[i] = 1e6; // should be Infinity + } + + // trellis starts at state 0 + metric[0] = 0; + + double *tmp_metric = malloc(N_states * sizeof *tmp_metric); + double *rho = malloc(code.components * sizeof *rho); + for (int k = 0; k < packet_length + code.memory; k++) { + + // get received symbol + for (int r = 0; r < code.components; r++) + rho[r] = received[k*code.components + r]; + + for (int s = 0; s < N_states; s++) { + + // get neighbors + int nA = abs(code.neighbors[s][0]) - 1; + int uA = (code.neighbors[s][0] > 0); + int nB = abs(code.neighbors[s][1]) - 1; + int uB = (code.neighbors[s][1] > 0); + + int *outA = code.output[nA][uA]; + int *outB = code.output[nB][uB]; + + double costA = 0; + double costB = 0; + for (int i = 0; i < code.components; i++) { + costA += pow(rho[i] - 2*outA[i] + 1, 2); + costB += pow(rho[i] - 2*outB[i] + 1, 2); + } + + costA += metric[nA]; + costB += metric[nB]; + + double minimum_cost = (costA > costB) ? costB : costA; + int idx = minimum_cost == costB; + tmp_metric[s] = minimum_cost; + + data_matrix[s][k] = code.neighbors[s][idx]; + } + + // find minimum + double min_metric = tmp_metric[0]; + for (int s = 0; s < N_states; s++) + min_metric = (min_metric < tmp_metric[s]) ? min_metric : tmp_metric[s]; + + // normalize + for (int s = 0; s < N_states; s++) + metric[s] = tmp_metric[s] - min_metric; + + } + + // backtrack + int state = 0; // trellis is terminated + for (int k = packet_length + code.memory - 1; k >= 0; k--) + { + int input = (data_matrix[state][k] > 0); + state = abs(data_matrix[state][k]) - 1; + + if (k < packet_length) + decoded_packet[k] = input; + } + + // free memory + free(metric); + free(rho); + free(tmp_metric); + + for (int i = 0; i < N_states; i++ ) + free(data_matrix[i]); + free(data_matrix); + + return decoded_packet;/*}}}*/ +} + +void print_neighbors(t_convcode code) +{ + int N_states = 2 << (code.memory - 1);/*{{{*/ + + for (int i = 0; i < 34; i++){ + if (i % 11) + printf("-"); + else + printf("+"); + } + printf("\n"); + printf("|%-10s|%-10s|%-10s|\n", "STATE", "NEIGHBOR", "INPUT"); + for (int i = 0; i < 34; i++){ + if (i % 11) + printf("-"); + else + printf("+"); + } + printf("\n"); + + for (int i = 0; i < N_states; i++) { + int s0 = abs(code.neighbors[i][0])-1; + int s1 = abs(code.neighbors[i][1])-1; + + int u0 = (code.neighbors[i][0] > 0) ? 1 : 0; + int u1 = (code.neighbors[i][1] > 0) ? 1 : 0; + + printf("|%-10s|%-10s|%-10d|\n", state2str(i, code.memory), state2str(s0, code.memory), u0); + printf("|%-10s|%-10s|%-10d|\n", state2str(i, code.memory), state2str(s1, code.memory), u1); + } + for (int i = 0; i < 34; i++){ + if (i % 11) + printf("-"); + else + printf("+"); + } + printf("\n");/*}}}*/ +} + +int *convcode_extrinsic(double *received, double length, double ***a_priori, t_convcode code, double noise_variance, + int decision) +{ + int N_states = 2 << (code.memory - 1);/*{{{*/ + int packet_length = (int) length / code.components - code.memory; + + long int threshold = 1e10; + // copy a priori probabilities on local array + double **app = malloc(2 * sizeof(double*));/*{{{*/ + + for (int i = 0; i < 2; ++i) + app[i] = malloc((packet_length + code.memory) * sizeof *app); + + for (int i = 0; i < packet_length; ++i){ + app[0][i] = (*a_priori)[0][i]; + app[1][i] = (*a_priori)[1][i]; + } + + for (int i = 0; i < code.memory; i++) { + app[0][packet_length + i] = log(0.5); + app[1][packet_length + i] = log(0.5); + } + + /*}}}*/ + + // initialize backward messages + double **backward = malloc(N_states * sizeof(double*));/*{{{*/ + for (int k = 0; k < N_states; ++k) { + backward[k] = malloc((packet_length + code.memory) * sizeof(double)); + backward[k][packet_length + code.memory - 1] = -threshold; + } + + backward[0][packet_length + code.memory - 1] = 0; + + double *rho = malloc(code.components * sizeof *rho); + + for (int i = packet_length + code.memory - 2; i >= 0; i--) { + + for (int j = 0; j < code.components; ++j) + rho[j] = received[code.components*(i+1) + j]; + + for (int s = 0; s < N_states; ++s) { + double B = -threshold; + + for (int u = 0; u < 2; ++u) { + int next = code.next_state[s][u]; + int *out = code.output[s][u]; + + double g = 0; + for (int j = 0; j < code.components; ++j) + g += pow(rho[j]- (2*out[j] - 1), 2); + + B = exp_sum(B, app[u][i+1] + backward[next][i+1] + (-g/(2*noise_variance))); + } + + backward[s][i] = B; + } + + // normalize + double max = backward[0][i]; + for (int s = 0; s < N_states; ++s) + max = backward[s][i] > max ? backward[s][i] : max; + + for (int s = 0; s < N_states; ++s) + backward[s][i] -= max; + }/*}}}*/ + + // initialize forward messages + double **forward = malloc(N_states * sizeof(double*));/*{{{*/ + for (int k = 0; k < N_states; ++k) { + forward[k] = malloc((packet_length + code.memory) * sizeof(double)); + forward[k][0] = -threshold; + } + forward[0][0] = 0; + + + for (int i = 1; i < packet_length + code.memory; ++i) { + + for (int j = 0; j < code.components; ++j) + rho[j] = received[code.components*(i-1) + j]; + + for (int s = 0; s < N_states; ++s) { + + double F = -threshold; + + // pass through each neighbour + int *neigh = code.neighbors[s]; + for (int n = 0; n < 2; ++n) { + int state = abs(neigh[n]) - 1; + int input = neigh[n] > 0; + + int *out = code.output[state][input]; + + double g = 0; + // compute g + for (int j = 0; j < code.components; ++j) + g += pow(rho[j] - (2*out[j] - 1),2); + + F = exp_sum(F, app[input][i-1] + forward[state][i-1] + (-g/(2*noise_variance))); + } + + forward[s][i] = F; + + } + + // normalize + double max = forward[0][i]; + for (int s = 0; s < N_states; ++s) + max = forward[s][i] > max ? forward[s][i] : max; + + for (int s = 0; s < N_states; ++s) + forward[s][i] -= max; + }/*}}}*/ + + // initialize extrinsic messages + double **extrinsic = malloc(2 * sizeof(double*));/*{{{*/ + + for (int k = 0; k < 2; ++k) { + extrinsic[k] = malloc((packet_length * code.memory) * sizeof(double)); + } + + for (int i = 0; i < packet_length + code.memory; ++i) { + for (int j = 0; j < code.components; ++j) + rho[j] = received[code.components*i + j]; + + + for (int u = 0; u < 2; ++u) { + double E = -threshold; + for (int s = 0; s < N_states; ++s) { + + int state = code.next_state[s][u]; + + double g = 0; + + int *out = code.output[s][u]; + for (int j = 0; j < code.components; ++j) + g += pow(rho[j] - (2*out[j] - 1),2); + + double fwd = forward[s][i]; + double bwd = backward[state][i]; + E = exp_sum(E, fwd + bwd + (-g/(2*noise_variance))); + } + + extrinsic[u][i] = E; + } + +// double normalization = log(exp(extrinsic[0][i]) + exp(extrinsic[1][i])); +// extrinsic[0][i] -= normalization; +// extrinsic[1][i] -= normalization; + if (i < packet_length) + { + (*a_priori)[0][i] = extrinsic[0][i]; + (*a_priori)[1][i] = extrinsic[1][i]; + } + }/*}}}*/ + + // decision + int *decoded = NULL; + + if (decision){ + decoded = malloc(packet_length * sizeof(int) ); //sizeof *decoded + for (int i = 0; i < packet_length; ++i) { + double one = app[1][i] + extrinsic[1][i]; + double zero = app[0][i] + extrinsic[0][i]; + decoded[i] = one > zero; + } + } + + // free memory + for (int l = 0; l < N_states; ++l) {/*{{{*/ + free(backward[l]); + free(forward[l]); + } + free(backward); + free(forward); + + for (int i = 0; i < 2; i++) { + free(extrinsic[i]); + free(app[i]); + } + free(extrinsic); + free(app); + free(rho);/*}}}*/ + + return decoded; + + /*}}}*/ +} + +static double exp_sum(double a, double b) +{ + double diff = a-b;/*{{{*/ + return (a > b) ? a : b + log(1 + exp(-diff > 0 ? diff : -diff));/*}}}*//*}}}*/ +} + diff --git a/src-core/libs/deepspace-turbo/libconvcodes.h b/src-core/libs/deepspace-turbo/libconvcodes.h new file mode 100644 index 000000000..cccf94417 --- /dev/null +++ b/src-core/libs/deepspace-turbo/libconvcodes.h @@ -0,0 +1,36 @@ +// +// Created by gianluca on 20/02/17. +// + +#ifndef DEEPSPACE_TURBO_LIBCONVCODES_H +#define DEEPSPACE_TURBO_LIBCONVCODES_H + +typedef struct str_convcode{ + int components; + int memory; + int **forward_connections; + int *backward_connections; + int **next_state; + int **neighbors; + int ***output; +} t_convcode; + +static int get_bit(int num, int position); +static char* state2str(int state, int memory); +static int convcode_stateupdate(int state, int input, t_convcode code); +static int *convcode_output(int state, int input, t_convcode code); + +t_convcode convcode_initialize(char *forward[], char *backward, int N_components); +void convcode_clear(t_convcode code); +int* convcode_encode(int *packet, int packet_length, t_convcode code); +int* convcode_decode(double *received, int length, t_convcode code); + +void print_neighbors(t_convcode code); + +// BCJR decoding +int * convcode_extrinsic(double *received, double length, double ***a_priori, t_convcode code, double noise_variance, + int decision); + +static double exp_sum(double a, double b); + +#endif //DEEPSPACE_TURBO_LIBCONVCODES_H diff --git a/src-core/libs/deepspace-turbo/libturbocodes.c b/src-core/libs/deepspace-turbo/libturbocodes.c new file mode 100644 index 000000000..81633cc7c --- /dev/null +++ b/src-core/libs/deepspace-turbo/libturbocodes.c @@ -0,0 +1,201 @@ +// +// Created by gianluca on 22/02/17. +// + +#include "libturbocodes.h" +//#include "utilities.h" +#include +#include + + +static int *turbo_interleave(int *packet, t_turbocode code) +{ + int *interleaved_packet = malloc(code.packet_length * sizeof(int));// {{{ + for (int j = 0; j < code.packet_length; ++j) { + interleaved_packet[j] = packet[code.interleaver[j]]; + } + + return interleaved_packet;// }}} +} + +static int *turbo_deinterleave(int *packet, t_turbocode code) +{ + int *local = malloc(code.packet_length*sizeof(int));// {{{ + for (int i = 0; i < code.packet_length; ++i) { + local[code.interleaver[i]] = packet[i]; + } + + return local;// }}} +} + +static void message_interleave(double ***messages, t_turbocode code) +{ + // local array// {{{ + double **local = malloc(2*sizeof(double*)); + local[0] = malloc(code.packet_length * sizeof(double)); + local[1] = malloc(code.packet_length * sizeof(double)); + + for (int i = 0; i < code.packet_length; ++i) { + local[0][i] = (*messages)[0][code.interleaver[i]]; + local[1][i] = (*messages)[1][code.interleaver[i]]; + } + + for (int i = 0; i < code.packet_length; ++i) { + (*messages)[0][i] = local[0][i]; + (*messages)[1][i] = local[1][i]; + } + + free(local[0]); + free(local[1]); + free(local);// }}} +} + +static void message_deinterleave(double ***messages, t_turbocode code) +{ + // local array// {{{ + double **local = malloc(2*sizeof(double*)); + local[0] = malloc(code.packet_length * sizeof(double)); + local[1] = malloc(code.packet_length * sizeof(double)); + + for (int i = 0; i < code.packet_length; ++i) { + local[0][code.interleaver[i]] = (*messages)[0][i]; + local[1][code.interleaver[i]] = (*messages)[1][i]; + } + + for (int i = 0; i < code.packet_length; ++i) { + (*messages)[0][i] = local[0][i]; + (*messages)[1][i] = local[1][i]; + } + + free(local[0]); + free(local[1]); + free(local);// }}} +} + + +t_turbocode turbo_initialize(t_convcode upper, t_convcode lower, int *interleaver, int packet_length) +{ + t_turbocode code;/*{{{*/ + code.upper_code = upper; + code.lower_code = lower; + + code.packet_length = packet_length; + code.interleaver = interleaver; + + // compute encoded length + int turbo_length = 0; + turbo_length += upper.components * (code.packet_length + upper.memory); + turbo_length += lower.components * (code.packet_length + lower.memory); + + code.encoded_length = turbo_length; + return code;/*}}}*/ +} + +int *turbo_encode(int *packet, t_turbocode code) +{ + int *interleaved_packet = turbo_interleave(packet, code);/*{{{*/ + + // reference to encoded messages + int **conv_encoded = malloc(2 * sizeof(int*)); + int turbo_length = code.encoded_length; + conv_encoded[0] = convcode_encode(packet, code.packet_length, code.upper_code); + conv_encoded[1] = convcode_encode(interleaved_packet, code.packet_length, code.lower_code); + + int *turbo_encoded = malloc(turbo_length * sizeof *turbo_encoded); + + t_convcode codes[2] = {code.upper_code, code.lower_code}; + // parallel to serial + + int k = 0, c = 0, cw = 0;/*{{{*/ + while (k < turbo_length) { + t_convcode cc = codes[c]; + + // number of components of cc + int comps = cc.components; + + // copy bits from cc output to turbo_encoded + for (int i = 0; i < comps; i++) { + int bit = conv_encoded[c][cw*comps + i]; + turbo_encoded[k++] = bit; + } + + c = (c + 1) % 2; + // when c = 0 the first codeword is complete + cw = !c ? cw+1 : cw; + }/*}}}*/ + + free(conv_encoded[0]); + free(conv_encoded[1]); + free(conv_encoded); + + free(interleaved_packet); + + + return turbo_encoded;/*}}}*/ +} + +int* turbo_decode(double *received, int iterations, double noise_variance, t_turbocode code) +{ + // serial to parallel/*{{{*/ + int *lengths = malloc(2 * sizeof *lengths);/*{{{*/ + double **streams = malloc(2 * sizeof(double*)); + t_convcode codes[2] = {code.upper_code, code.lower_code}; + for (int i = 0; i < 2; i++) { + t_convcode cc = codes[i]; + lengths[i] = cc.components * (code.packet_length + cc.memory); + streams[i] = malloc(lengths[i] * sizeof(double)); + } + + int k = 0, c = 0, cw = 0; + while (k < code.encoded_length) { + t_convcode cc = codes[c]; + + for (int i = 0; i < cc.components; i++) + streams[c][cw*cc.components + i] = received[k++]; + + c = (c + 1) % 2; + cw = !c ? cw + 1 : cw; + }/*}}}*/ + + // initial messages + double **messages = malloc(2 * sizeof(double *)); + for (int i = 0; i < 2; i++) { + messages[i] = malloc(code.packet_length * sizeof(double)); + for (int j = 0; j < code.packet_length; j++) { + messages[i][j] = log(0.5); + } + } + + int *turbo_decoded = NULL; + int *turbo_decoded_1 = NULL; + for (int i = 0; i < iterations; i++) { + + // run BCJR on upper code + turbo_decoded_1 = convcode_extrinsic(streams[0], lengths[0], &messages, code.upper_code, noise_variance, 0); + + // apply interleaver + message_interleave(&messages, code); + + // run BCJR on lower code + turbo_decoded = convcode_extrinsic(streams[1], lengths[1], &messages, code.lower_code, noise_variance, i == (iterations - 1)); + + // deinterleave + message_deinterleave(&messages, code); + } + + int *decoded_deinterleaved = turbo_deinterleave(turbo_decoded, code); + //decoded_deinterleaved = turbo_deinterleave(turbo_decoded, code); + + for (int i = 0; i < 2; i++) + free(streams[i]); + free(streams); + free(turbo_decoded); + free(turbo_decoded_1); + free(lengths); + free(messages[0]); + free(messages[1]); + free(messages); + + //length of the + return decoded_deinterleaved; /*}}}*/ +} diff --git a/src-core/libs/deepspace-turbo/libturbocodes.h b/src-core/libs/deepspace-turbo/libturbocodes.h new file mode 100644 index 000000000..cd7236e21 --- /dev/null +++ b/src-core/libs/deepspace-turbo/libturbocodes.h @@ -0,0 +1,30 @@ +// +// Created by gianluca on 22/02/17. +// + +#ifndef DEEPSPACE_TURBO_LIBTURBOCODES_H +#define DEEPSPACE_TURBO_LIBTURBOCODES_H + +#include "libconvcodes.h" + +typedef struct str_turbocode{ + t_convcode upper_code; + t_convcode lower_code; + + int *interleaver; + int packet_length; + int encoded_length; +} t_turbocode; + +static int *turbo_interleave(int *packet, t_turbocode code); +static int *turbo_deinterleave(int *packet, t_turbocode code); +static void message_interleave(double ***messages, t_turbocode code); +static void message_deinterleave(double ***messages, t_turbocode code); + +t_turbocode turbo_initialize(t_convcode upper, t_convcode lower, int *interleaver, int packet_length); + +int *turbo_encode(int *packet, t_turbocode code); + +int *turbo_decode(double* received, int iterations, double noise_variance, t_turbocode code); +//void turbo_decode(double *received, int iterations, double noise_variance, t_turbocode code, int *decoded_deinterleaved); +#endif //DEEPSPACE_TURBO_LIBTURBOCODES_H diff --git a/src-testing/main.cpp b/src-testing/main.cpp index 820e52312..82410e45f 100644 --- a/src-testing/main.cpp +++ b/src-testing/main.cpp @@ -11,8 +11,45 @@ **********************************************************************/ #include "logger.h" +#include "common/codings/turbo/ccsds_turbo.h" +#include int main(int argc, char *argv[]) { initLogger(); + + codings::turbo::CCSDSTurbo ccsds_turbo(codings::turbo::BASE_223, codings::turbo::RATE_1_2); + + logger->critical("Frame length {:d}", ccsds_turbo.frame_length()); + logger->critical("Codeword length {:d}", ccsds_turbo.codeword_length()); + + uint8_t turbo_frame[223]; + uint8_t turbo_codeword[3576 / 8]; + std::ifstream ts_in("/home/alan/Downloads/sk8.ts"); + std::ofstream encoded_in("/home/alan/encoded_turbo.ts"); + + float turbo_soft_frame[3576]; + + while (!ts_in.eof()) + { + ts_in.read((char *)turbo_frame, 188); + + ccsds_turbo.encode(turbo_frame, turbo_codeword); + + for (int i = 0; i < 3576; i++) + { + uint8_t bit = (turbo_codeword[i / 8] >> (7 - (i % 8))) & 1; + turbo_soft_frame[i] = bit ? 1 : -1; + } + + for (int i = 0; i < 100; i++) + { + int pos = rand() % 3576; + turbo_soft_frame[pos] = -turbo_soft_frame[pos]; + } + + ccsds_turbo.decode(turbo_soft_frame, turbo_frame, 5); + + encoded_in.write((char *)turbo_frame, 188); + } }