import Reed-Solomon LGPL implementation for use on packet radio.

Encodes 223 bytes in 255 bytes, allowing upto 16 errors.
This commit is contained in:
gardners 2013-08-30 16:46:03 +09:30 committed by Jeremy Lakeman
parent 91fcf8b1b3
commit c28fdcf558
15 changed files with 1170 additions and 1 deletions

2
fec-3.0.1/README Normal file
View File

@ -0,0 +1,2 @@
Code from http://www.ka9q.net/code/fec/
Used under the LGPL.

43
fec-3.0.1/ccsds_tables.c Normal file
View File

@ -0,0 +1,43 @@
char CCSDS_alpha_to[] = {
0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80,0x87,0x89,0x95,0xad,0xdd,0x3d,0x7a,0xf4,
0x6f,0xde,0x3b,0x76,0xec,0x5f,0xbe,0xfb,0x71,0xe2,0x43,0x86,0x8b,0x91,0xa5,0xcd,
0x1d,0x3a,0x74,0xe8,0x57,0xae,0xdb,0x31,0x62,0xc4,0x0f,0x1e,0x3c,0x78,0xf0,0x67,
0xce,0x1b,0x36,0x6c,0xd8,0x37,0x6e,0xdc,0x3f,0x7e,0xfc,0x7f,0xfe,0x7b,0xf6,0x6b,
0xd6,0x2b,0x56,0xac,0xdf,0x39,0x72,0xe4,0x4f,0x9e,0xbb,0xf1,0x65,0xca,0x13,0x26,
0x4c,0x98,0xb7,0xe9,0x55,0xaa,0xd3,0x21,0x42,0x84,0x8f,0x99,0xb5,0xed,0x5d,0xba,
0xf3,0x61,0xc2,0x03,0x06,0x0c,0x18,0x30,0x60,0xc0,0x07,0x0e,0x1c,0x38,0x70,0xe0,
0x47,0x8e,0x9b,0xb1,0xe5,0x4d,0x9a,0xb3,0xe1,0x45,0x8a,0x93,0xa1,0xc5,0x0d,0x1a,
0x34,0x68,0xd0,0x27,0x4e,0x9c,0xbf,0xf9,0x75,0xea,0x53,0xa6,0xcb,0x11,0x22,0x44,
0x88,0x97,0xa9,0xd5,0x2d,0x5a,0xb4,0xef,0x59,0xb2,0xe3,0x41,0x82,0x83,0x81,0x85,
0x8d,0x9d,0xbd,0xfd,0x7d,0xfa,0x73,0xe6,0x4b,0x96,0xab,0xd1,0x25,0x4a,0x94,0xaf,
0xd9,0x35,0x6a,0xd4,0x2f,0x5e,0xbc,0xff,0x79,0xf2,0x63,0xc6,0x0b,0x16,0x2c,0x58,
0xb0,0xe7,0x49,0x92,0xa3,0xc1,0x05,0x0a,0x14,0x28,0x50,0xa0,0xc7,0x09,0x12,0x24,
0x48,0x90,0xa7,0xc9,0x15,0x2a,0x54,0xa8,0xd7,0x29,0x52,0xa4,0xcf,0x19,0x32,0x64,
0xc8,0x17,0x2e,0x5c,0xb8,0xf7,0x69,0xd2,0x23,0x46,0x8c,0x9f,0xb9,0xf5,0x6d,0xda,
0x33,0x66,0xcc,0x1f,0x3e,0x7c,0xf8,0x77,0xee,0x5b,0xb6,0xeb,0x51,0xa2,0xc3,0x00,
};
char CCSDS_index_of[] = {
255, 0, 1, 99, 2,198,100,106, 3,205,199,188,101,126,107, 42,
4,141,206, 78,200,212,189,225,102,221,127, 49,108, 32, 43,243,
5, 87,142,232,207,172, 79,131,201,217,213, 65,190,148,226,180,
103, 39,222,240,128,177, 50, 53,109, 69, 33, 18, 44, 13,244, 56,
6,155, 88, 26,143,121,233,112,208,194,173,168, 80,117,132, 72,
202,252,218,138,214, 84, 66, 36,191,152,149,249,227, 94,181, 21,
104, 97, 40,186,223, 76,241, 47,129,230,178, 63, 51,238, 54, 16,
110, 24, 70,166, 34,136, 19,247, 45,184, 14, 61,245,164, 57, 59,
7,158,156,157, 89,159, 27, 8,144, 9,122, 28,234,160,113, 90,
209, 29,195,123,174, 10,169,145, 81, 91,118,114,133,161, 73,235,
203,124,253,196,219, 30,139,210,215,146, 85,170, 67, 11, 37,175,
192,115,153,119,150, 92,250, 82,228,236, 95, 74,182,162, 22,134,
105,197, 98,254, 41,125,187,204,224,211, 77,140,242, 31, 48,220,
130,171,231, 86,179,147, 64,216, 52,176,239, 38, 55, 12, 17, 68,
111,120, 25,154, 71,116,167,193, 35, 83,137,251, 20, 93,248,151,
46, 75,185, 96, 15,237, 62,229,246,135,165, 23, 58,163, 60,183,
};
char CCSDS_poly[] = {
0,249, 59, 66, 4, 43,126,251, 97, 30, 3,213, 50, 66,170, 5,
24, 5,170, 66, 50,213, 3, 30, 97,251,126, 43, 4, 66, 59,249,
0,
};

24
fec-3.0.1/char.h Normal file
View File

@ -0,0 +1,24 @@
/* Stuff specific to the 8-bit symbol version of the general purpose RS codecs
*
* Copyright 2003, Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
typedef unsigned char data_t;
#define MODNN(x) modnn(rs,x)
#define MM (rs->mm)
#define NN (rs->nn)
#define ALPHA_TO (rs->alpha_to)
#define INDEX_OF (rs->index_of)
#define GENPOLY (rs->genpoly)
#define NROOTS (rs->nroots)
#define FCR (rs->fcr)
#define PRIM (rs->prim)
#define IPRIM (rs->iprim)
#define PAD (rs->pad)
#define A0 (NN)

298
fec-3.0.1/decode_rs.h Normal file
View File

@ -0,0 +1,298 @@
/* The guts of the Reed-Solomon decoder, meant to be #included
* into a function body with the following typedefs, macros and variables supplied
* according to the code parameters:
* data_t - a typedef for the data symbol
* data_t data[] - array of NN data and parity symbols to be corrected in place
* retval - an integer lvalue into which the decoder's return code is written
* NROOTS - the number of roots in the RS code generator polynomial,
* which is the same as the number of parity symbols in a block.
Integer variable or literal.
* NN - the total number of symbols in a RS block. Integer variable or literal.
* PAD - the number of pad symbols in a block. Integer variable or literal.
* ALPHA_TO - The address of an array of NN elements to convert Galois field
* elements in index (log) form to polynomial form. Read only.
* INDEX_OF - The address of an array of NN elements to convert Galois field
* elements in polynomial form to index (log) form. Read only.
* MODNN - a function to reduce its argument modulo NN. May be inline or a macro.
* FCR - An integer literal or variable specifying the first consecutive root of the
* Reed-Solomon generator polynomial. Integer variable or literal.
* PRIM - The primitive root of the generator poly. Integer variable or literal.
* DEBUG - If set to 1 or more, do various internal consistency checking. Leave this
* undefined for production code
* The memset(), memmove(), and memcpy() functions are used. The appropriate header
* file declaring these functions (usually <string.h>) must be included by the calling
* program.
*/
#if !defined(NROOTS)
#error "NROOTS not defined"
#endif
#if !defined(NN)
#error "NN not defined"
#endif
#if !defined(PAD)
#error "PAD not defined"
#endif
#if !defined(ALPHA_TO)
#error "ALPHA_TO not defined"
#endif
#if !defined(INDEX_OF)
#error "INDEX_OF not defined"
#endif
#if !defined(MODNN)
#error "MODNN not defined"
#endif
#if !defined(FCR)
#error "FCR not defined"
#endif
#if !defined(PRIM)
#error "PRIM not defined"
#endif
#if !defined(NULL)
#define NULL ((void *)0)
#endif
#undef MIN
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#undef A0
#define A0 (NN)
{
int deg_lambda, el, deg_omega;
int i, j, r,k;
data_t u,q,tmp,num1,num2,den,discr_r;
data_t lambda[NROOTS+1], s[NROOTS]; /* Err+Eras Locator poly
* and syndrome poly */
data_t b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1];
data_t root[NROOTS], reg[NROOTS+1], loc[NROOTS];
int syn_error, count;
/* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
for(i=0;i<NROOTS;i++)
s[i] = data[0];
for(j=1;j<NN-PAD;j++){
for(i=0;i<NROOTS;i++){
if(s[i] == 0){
s[i] = data[j];
} else {
s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR+i)*PRIM)];
}
}
}
/* Convert syndromes to index form, checking for nonzero condition */
syn_error = 0;
for(i=0;i<NROOTS;i++){
syn_error |= s[i];
s[i] = INDEX_OF[s[i]];
}
if (!syn_error) {
/* if syndrome is zero, data[] is a codeword and there are no
* errors to correct. So return data[] unmodified
*/
count = 0;
goto finish;
}
memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
lambda[0] = 1;
if (no_eras > 0) {
/* Init lambda to be the erasure locator polynomial */
lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
for (i = 1; i < no_eras; i++) {
u = MODNN(PRIM*(NN-1-eras_pos[i]));
for (j = i+1; j > 0; j--) {
tmp = INDEX_OF[lambda[j - 1]];
if(tmp != A0)
lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
}
}
#if DEBUG >= 1
/* Test code that verifies the erasure locator polynomial just constructed
Needed only for decoder debugging. */
/* find roots of the erasure location polynomial */
for(i=1;i<=no_eras;i++)
reg[i] = INDEX_OF[lambda[i]];
count = 0;
for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
q = 1;
for (j = 1; j <= no_eras; j++)
if (reg[j] != A0) {
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
}
if (q != 0)
continue;
/* store root and error location number indices */
root[count] = i;
loc[count] = k;
count++;
}
if (count != no_eras) {
printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras);
count = -1;
goto finish;
}
#if DEBUG >= 2
printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
for (i = 0; i < count; i++)
printf("%d ", loc[i]);
printf("\n");
#endif
#endif
}
for(i=0;i<NROOTS+1;i++)
b[i] = INDEX_OF[lambda[i]];
/*
* Begin Berlekamp-Massey algorithm to determine error+erasure
* locator polynomial
*/
r = no_eras;
el = no_eras;
while (++r <= NROOTS) { /* r is the step number */
/* Compute discrepancy at the r-th step in poly-form */
discr_r = 0;
for (i = 0; i < r; i++){
if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])];
}
}
discr_r = INDEX_OF[discr_r]; /* Index form */
if (discr_r == A0) {
/* 2 lines below: B(x) <-- x*B(x) */
memmove(&b[1],b,NROOTS*sizeof(b[0]));
b[0] = A0;
} else {
/* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
t[0] = lambda[0];
for (i = 0 ; i < NROOTS; i++) {
if(b[i] != A0)
t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
else
t[i+1] = lambda[i+1];
}
if (2 * el <= r + no_eras - 1) {
el = r + no_eras - el;
/*
* 2 lines below: B(x) <-- inv(discr_r) *
* lambda(x)
*/
for (i = 0; i <= NROOTS; i++)
b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
} else {
/* 2 lines below: B(x) <-- x*B(x) */
memmove(&b[1],b,NROOTS*sizeof(b[0]));
b[0] = A0;
}
memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
}
}
/* Convert lambda to index form and compute deg(lambda(x)) */
deg_lambda = 0;
for(i=0;i<NROOTS+1;i++){
lambda[i] = INDEX_OF[lambda[i]];
if(lambda[i] != A0)
deg_lambda = i;
}
/* Find roots of the error+erasure locator polynomial by Chien search */
memcpy(&reg[1],&lambda[1],NROOTS*sizeof(reg[0]));
count = 0; /* Number of roots of lambda(x) */
for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
q = 1; /* lambda[0] is always 0 */
for (j = deg_lambda; j > 0; j--){
if (reg[j] != A0) {
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
}
}
if (q != 0)
continue; /* Not a root */
/* store root (index-form) and error location number */
#if DEBUG>=2
printf("count %d root %d loc %d\n",count,i,k);
#endif
root[count] = i;
loc[count] = k;
/* If we've already found max possible roots,
* abort the search to save time
*/
if(++count == deg_lambda)
break;
}
if (deg_lambda != count) {
/*
* deg(lambda) unequal to number of roots => uncorrectable
* error detected
*/
count = -1;
goto finish;
}
/*
* Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
* x**NROOTS). in index form. Also find deg(omega).
*/
deg_omega = deg_lambda-1;
for (i = 0; i <= deg_omega;i++){
tmp = 0;
for(j=i;j >= 0; j--){
if ((s[i - j] != A0) && (lambda[j] != A0))
tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
}
omega[i] = INDEX_OF[tmp];
}
/*
* Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
* inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
*/
for (j = count-1; j >=0; j--) {
num1 = 0;
for (i = deg_omega; i >= 0; i--) {
if (omega[i] != A0)
num1 ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
}
num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
den = 0;
/* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
for (i = MIN(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
if(lambda[i+1] != A0)
den ^= ALPHA_TO[MODNN(lambda[i+1] + i * root[j])];
}
#if DEBUG >= 1
if (den == 0) {
printf("\n ERROR: denominator = 0\n");
count = -1;
goto finish;
}
#endif
/* Apply error to data */
if (num1 != 0 && loc[j] >= PAD) {
data[loc[j]-PAD] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
}
}
finish:
if(eras_pos != NULL){
for(i=0;i<count;i++)
eras_pos[i] = loc[i];
}
retval = count;
}

24
fec-3.0.1/decode_rs_8.c Normal file
View File

@ -0,0 +1,24 @@
/* General purpose Reed-Solomon decoder for 8-bit symbols or less
* Copyright 2003 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#ifdef DEBUG
#include <stdio.h>
#endif
#include <string.h>
#include "fixed.h"
int decode_rs_8(data_t *data, int *eras_pos, int no_eras, int pad){
int retval;
if(pad < 0 || pad > 222){
return -1;
}
#include "decode_rs.h"
return retval;
}

58
fec-3.0.1/encode_rs.h Normal file
View File

@ -0,0 +1,58 @@
/* The guts of the Reed-Solomon encoder, meant to be #included
* into a function body with the following typedefs, macros and variables supplied
* according to the code parameters:
* data_t - a typedef for the data symbol
* data_t data[] - array of NN-NROOTS-PAD and type data_t to be encoded
* data_t parity[] - an array of NROOTS and type data_t to be written with parity symbols
* NROOTS - the number of roots in the RS code generator polynomial,
* which is the same as the number of parity symbols in a block.
Integer variable or literal.
*
* NN - the total number of symbols in a RS block. Integer variable or literal.
* PAD - the number of pad symbols in a block. Integer variable or literal.
* ALPHA_TO - The address of an array of NN elements to convert Galois field
* elements in index (log) form to polynomial form. Read only.
* INDEX_OF - The address of an array of NN elements to convert Galois field
* elements in polynomial form to index (log) form. Read only.
* MODNN - a function to reduce its argument modulo NN. May be inline or a macro.
* GENPOLY - an array of NROOTS+1 elements containing the generator polynomial in index form
* The memset() and memmove() functions are used. The appropriate header
* file declaring these functions (usually <string.h>) must be included by the calling
* program.
* Copyright 2004, Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#undef A0
#define A0 (NN) /* Special reserved value encoding zero in index form */
{
int i, j;
data_t feedback;
memset(parity,0,NROOTS*sizeof(data_t));
for(i=0;i<NN-NROOTS-PAD;i++){
feedback = INDEX_OF[data[i] ^ parity[0]];
if(feedback != A0){ /* feedback term is non-zero */
#ifdef UNNORMALIZED
/* This line is unnecessary when GENPOLY[NROOTS] is unity, as it must
* always be for the polynomials constructed by init_rs()
*/
feedback = MODNN(NN - GENPOLY[NROOTS] + feedback);
#endif
for(j=1;j<NROOTS;j++)
parity[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS-j])];
}
/* Shift */
memmove(&parity[0],&parity[1],sizeof(data_t)*(NROOTS-1));
if(feedback != A0)
parity[NROOTS-1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])];
else
parity[NROOTS-1] = 0;
}
}

81
fec-3.0.1/encode_rs_8.c Normal file
View File

@ -0,0 +1,81 @@
/* Reed-Solomon encoder
* Copyright 2004, Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#include <string.h>
#include "fixed.h"
#ifdef __VEC__
#include <sys/sysctl.h>
#endif
static enum {UNKNOWN=0,MMX,SSE,SSE2,ALTIVEC,PORT} cpu_mode;
static void encode_rs_8_c(data_t *data, data_t *parity,int pad);
#if __vec__
static void encode_rs_8_av(data_t *data, data_t *parity,int pad);
#endif
void encode_rs_8(data_t *data, data_t *parity,int pad){
if(cpu_mode == UNKNOWN){
cpu_mode = PORT;
}
switch(cpu_mode){
#if __vec__
case ALTIVEC:
encode_rs_8_av(data,parity,pad);
return;
#endif
#if __i386__
case MMX:
case SSE:
case SSE2:
#endif
default:
encode_rs_8_c(data,parity,pad);
return;
}
}
#if __vec__ /* PowerPC G4/G5 Altivec instructions are available */
static vector unsigned char reverse = (vector unsigned char)(0,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1);
static vector unsigned char shift_right = (vector unsigned char)(15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30);
/* Lookup table for feedback multiplications
* These are the low half of the coefficients. Since the generator polynomial is
* palindromic, we form the other half by reversing this one
*/
extern static union { vector unsigned char v; unsigned char c[16]; } table[256];
static void encode_rs_8_av(data_t *data, data_t *parity,int pad){
union { vector unsigned char v[2]; unsigned char c[32]; } shift_register;
int i;
shift_register.v[0] = (vector unsigned char)(0);
shift_register.v[1] = (vector unsigned char)(0);
for(i=0;i<NN-NROOTS-pad;i++){
vector unsigned char feedback0,feedback1;
unsigned char f;
f = data[i] ^ shift_register.c[31];
feedback1 = table[f].v;
feedback0 = vec_perm(feedback1,feedback1,reverse);
/* Shift right one byte */
shift_register.v[1] = vec_perm(shift_register.v[0],shift_register.v[1],shift_right) ^ feedback1;
shift_register.v[0] = vec_sro(shift_register.v[0],(vector unsigned char)(8)) ^ feedback0;
shift_register.c[0] = f;
}
for(i=0;i<NROOTS;i++)
parity[NROOTS-i-1] = shift_register.c[i];
}
#endif
/* Portable C version */
static void encode_rs_8_c(data_t *data, data_t *parity,int pad){
#include "encode_rs.h"
}

347
fec-3.0.1/fec.h Normal file
View File

@ -0,0 +1,347 @@
/* User include file for libfec
* Copyright 2004, Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#ifndef _FEC_H_
#define _FEC_H_
/* r=1/2 k=7 convolutional encoder polynomials
* The NASA-DSN convention is to use V27POLYA inverted, then V27POLYB
* The CCSDS/NASA-GSFC convention is to use V27POLYB, then V27POLYA inverted
*/
#define V27POLYA 0x6d
#define V27POLYB 0x4f
void *create_viterbi27(int len);
void set_viterbi27_polynomial(int polys[2]);
int init_viterbi27(void *vp,int starting_state);
int update_viterbi27_blk(void *vp,unsigned char sym[],int npairs);
int chainback_viterbi27(void *vp, unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi27(void *vp);
#ifdef __VEC__
void *create_viterbi27_av(int len);
void set_viterbi27_polynomial_av(int polys[2]);
int init_viterbi27_av(void *p,int starting_state);
int chainback_viterbi27_av(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi27_av(void *p);
int update_viterbi27_blk_av(void *p,unsigned char *syms,int nbits);
#endif
#ifdef __i386__
void *create_viterbi27_mmx(int len);
void set_viterbi27_polynomial_mmx(int polys[2]);
int init_viterbi27_mmx(void *p,int starting_state);
int chainback_viterbi27_mmx(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi27_mmx(void *p);
int update_viterbi27_blk_mmx(void *p,unsigned char *syms,int nbits);
void *create_viterbi27_sse(int len);
void set_viterbi27_polynomial_sse(int polys[2]);
int init_viterbi27_sse(void *p,int starting_state);
int chainback_viterbi27_sse(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi27_sse(void *p);
int update_viterbi27_blk_sse(void *p,unsigned char *syms,int nbits);
void *create_viterbi27_sse2(int len);
void set_viterbi27_polynomial_sse2(int polys[2]);
int init_viterbi27_sse2(void *p,int starting_state);
int chainback_viterbi27_sse2(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi27_sse2(void *p);
int update_viterbi27_blk_sse2(void *p,unsigned char *syms,int nbits);
#endif
void *create_viterbi27_port(int len);
void set_viterbi27_polynomial_port(int polys[2]);
int init_viterbi27_port(void *p,int starting_state);
int chainback_viterbi27_port(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi27_port(void *p);
int update_viterbi27_blk_port(void *p,unsigned char *syms,int nbits);
/* r=1/2 k=9 convolutional encoder polynomials */
#define V29POLYA 0x1af
#define V29POLYB 0x11d
void *create_viterbi29(int len);
void set_viterbi29_polynomial(int polys[2]);
int init_viterbi29(void *vp,int starting_state);
int update_viterbi29_blk(void *vp,unsigned char syms[],int nbits);
int chainback_viterbi29(void *vp, unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi29(void *vp);
#ifdef __VEC__
void *create_viterbi29_av(int len);
void set_viterbi29_polynomial_av(int polys[2]);
int init_viterbi29_av(void *p,int starting_state);
int chainback_viterbi29_av(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi29_av(void *p);
int update_viterbi29_blk_av(void *p,unsigned char *syms,int nbits);
#endif
#ifdef __i386__
void *create_viterbi29_mmx(int len);
void set_viterbi29_polynomial_mmx(int polys[2]);
int init_viterbi29_mmx(void *p,int starting_state);
int chainback_viterbi29_mmx(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi29_mmx(void *p);
int update_viterbi29_blk_mmx(void *p,unsigned char *syms,int nbits);
void *create_viterbi29_sse(int len);
void set_viterbi29_polynomial_sse(int polys[2]);
int init_viterbi29_sse(void *p,int starting_state);
int chainback_viterbi29_sse(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi29_sse(void *p);
int update_viterbi29_blk_sse(void *p,unsigned char *syms,int nbits);
void *create_viterbi29_sse2(int len);
void set_viterbi29_polynomial_sse2(int polys[2]);
int init_viterbi29_sse2(void *p,int starting_state);
int chainback_viterbi29_sse2(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi29_sse2(void *p);
int update_viterbi29_blk_sse2(void *p,unsigned char *syms,int nbits);
#endif
void *create_viterbi29_port(int len);
void set_viterbi29_polynomial_port(int polys[2]);
int init_viterbi29_port(void *p,int starting_state);
int chainback_viterbi29_port(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi29_port(void *p);
int update_viterbi29_blk_port(void *p,unsigned char *syms,int nbits);
/* r=1/3 k=9 convolutional encoder polynomials */
#define V39POLYA 0x1ed
#define V39POLYB 0x19b
#define V39POLYC 0x127
void *create_viterbi39(int len);
void set_viterbi39_polynomial(int polys[3]);
int init_viterbi39(void *vp,int starting_state);
int update_viterbi39_blk(void *vp,unsigned char syms[],int nbits);
int chainback_viterbi39(void *vp, unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi39(void *vp);
#ifdef __VEC__
void *create_viterbi39_av(int len);
void set_viterbi39_polynomial_av(int polys[3]);
int init_viterbi39_av(void *p,int starting_state);
int chainback_viterbi39_av(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi39_av(void *p);
int update_viterbi39_blk_av(void *p,unsigned char *syms,int nbits);
#endif
#ifdef __i386__
void *create_viterbi39_mmx(int len);
void set_viterbi39_polynomial_mmx(int polys[3]);
int init_viterbi39_mmx(void *p,int starting_state);
int chainback_viterbi39_mmx(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi39_mmx(void *p);
int update_viterbi39_blk_mmx(void *p,unsigned char *syms,int nbits);
void *create_viterbi39_sse(int len);
void set_viterbi39_polynomial_sse(int polys[3]);
int init_viterbi39_sse(void *p,int starting_state);
int chainback_viterbi39_sse(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi39_sse(void *p);
int update_viterbi39_blk_sse(void *p,unsigned char *syms,int nbits);
void *create_viterbi39_sse2(int len);
void set_viterbi39_polynomial_sse2(int polys[3]);
int init_viterbi39_sse2(void *p,int starting_state);
int chainback_viterbi39_sse2(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi39_sse2(void *p);
int update_viterbi39_blk_sse2(void *p,unsigned char *syms,int nbits);
#endif
void *create_viterbi39_port(int len);
void set_viterbi39_polynomial_port(int polys[3]);
int init_viterbi39_port(void *p,int starting_state);
int chainback_viterbi39_port(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi39_port(void *p);
int update_viterbi39_blk_port(void *p,unsigned char *syms,int nbits);
/* r=1/6 k=15 Cassini convolutional encoder polynomials without symbol inversion
* dfree = 56
* These bits may be left-right flipped from some textbook representations;
* here I have the bits entering the shift register from the right (low) end
*
* Some other spacecraft use the same code, but with the polynomials in a different order.
* E.g., Mars Pathfinder and STEREO swap POLYC and POLYD. All use alternate symbol inversion,
* so use set_viterbi615_polynomial() as appropriate.
*/
#define V615POLYA 042631
#define V615POLYB 047245
#define V615POLYC 056507
#define V615POLYD 073363
#define V615POLYE 077267
#define V615POLYF 064537
void *create_viterbi615(int len);
void set_viterbi615_polynomial(int polys[6]);
int init_viterbi615(void *vp,int starting_state);
int update_viterbi615_blk(void *vp,unsigned char *syms,int nbits);
int chainback_viterbi615(void *vp, unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi615(void *vp);
#ifdef __VEC__
void *create_viterbi615_av(int len);
void set_viterbi615_polynomial_av(int polys[6]);
int init_viterbi615_av(void *p,int starting_state);
int chainback_viterbi615_av(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi615_av(void *p);
int update_viterbi615_blk_av(void *p,unsigned char *syms,int nbits);
#endif
#ifdef __i386__
void *create_viterbi615_mmx(int len);
void set_viterbi615_polynomial_mmx(int polys[6]);
int init_viterbi615_mmx(void *p,int starting_state);
int chainback_viterbi615_mmx(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi615_mmx(void *p);
int update_viterbi615_blk_mmx(void *p,unsigned char *syms,int nbits);
void *create_viterbi615_sse(int len);
void set_viterbi615_polynomial_sse(int polys[6]);
int init_viterbi615_sse(void *p,int starting_state);
int chainback_viterbi615_sse(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi615_sse(void *p);
int update_viterbi615_blk_sse(void *p,unsigned char *syms,int nbits);
void *create_viterbi615_sse2(int len);
void set_viterbi615_polynomial_sse2(int polys[6]);
int init_viterbi615_sse2(void *p,int starting_state);
int chainback_viterbi615_sse2(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi615_sse2(void *p);
int update_viterbi615_blk_sse2(void *p,unsigned char *syms,int nbits);
#endif
void *create_viterbi615_port(int len);
void set_viterbi615_polynomial_port(int polys[6]);
int init_viterbi615_port(void *p,int starting_state);
int chainback_viterbi615_port(void *p,unsigned char *data,unsigned int nbits,unsigned int endstate);
void delete_viterbi615_port(void *p);
int update_viterbi615_blk_port(void *p,unsigned char *syms,int nbits);
/* General purpose RS codec, 8-bit symbols */
void encode_rs_char(void *rs,unsigned char *data,unsigned char *parity);
int decode_rs_char(void *rs,unsigned char *data,int *eras_pos,
int no_eras);
void *init_rs_char(int symsize,int gfpoly,
int fcr,int prim,int nroots,
int pad);
void free_rs_char(void *rs);
/* General purpose RS codec, integer symbols */
void encode_rs_int(void *rs,int *data,int *parity);
int decode_rs_int(void *rs,int *data,int *eras_pos,int no_eras);
void *init_rs_int(int symsize,int gfpoly,int fcr,
int prim,int nroots,int pad);
void free_rs_int(void *rs);
/* CCSDS standard (255,223) RS codec with conventional (*not* dual-basis)
* symbol representation
*/
void encode_rs_8(unsigned char *data,unsigned char *parity,int pad);
int decode_rs_8(unsigned char *data,int *eras_pos,int no_eras,int pad);
/* CCSDS standard (255,223) RS codec with dual-basis symbol representation */
void encode_rs_ccsds(unsigned char *data,unsigned char *parity,int pad);
int decode_rs_ccsds(unsigned char *data,int *eras_pos,int no_eras,int pad);
/* Tables to map from conventional->dual (Taltab) and
* dual->conventional (Tal1tab) bases
*/
extern unsigned char Taltab[],Tal1tab[];
/* CPU SIMD instruction set available */
extern enum cpu_mode {UNKNOWN=0,PORT,MMX,SSE,SSE2,ALTIVEC} Cpu_mode;
void find_cpu_mode(void); /* Call this once at startup to set Cpu_mode */
/* Determine parity of argument: 1 = odd, 0 = even */
#ifdef __i386__
static inline int parityb(unsigned char x){
__asm__ __volatile__ ("test %1,%1;setpo %0" : "=g"(x) : "r" (x));
return x;
}
#else
void partab_init();
static inline int parityb(unsigned char x){
extern unsigned char Partab[256];
extern int P_init;
if(!P_init){
partab_init();
}
return Partab[x];
}
#endif
static inline int parity(int x){
/* Fold down to one byte */
x ^= (x >> 16);
x ^= (x >> 8);
return parityb(x);
}
/* Useful utilities for simulation */
double normal_rand(double mean, double std_dev);
unsigned char addnoise(int sym,double amp,double gain,double offset,int clip);
extern int Bitcnt[];
/* Dot product functions */
void *initdp(signed short coeffs[],int len);
void freedp(void *dp);
long dotprod(void *dp,signed short a[]);
void *initdp_port(signed short coeffs[],int len);
void freedp_port(void *dp);
long dotprod_port(void *dp,signed short a[]);
#ifdef __i386__
void *initdp_mmx(signed short coeffs[],int len);
void freedp_mmx(void *dp);
long dotprod_mmx(void *dp,signed short a[]);
void *initdp_sse(signed short coeffs[],int len);
void freedp_sse(void *dp);
long dotprod_sse(void *dp,signed short a[]);
void *initdp_sse2(signed short coeffs[],int len);
void freedp_sse2(void *dp);
long dotprod_sse2(void *dp,signed short a[]);
#endif
#ifdef __VEC__
void *initdp_av(signed short coeffs[],int len);
void freedp_av(void *dp);
long dotprod_av(void *dp,signed short a[]);
#endif
/* Sum of squares - accepts signed shorts, produces unsigned long long */
unsigned long long sumsq(signed short *in,int cnt);
unsigned long long sumsq_port(signed short *in,int cnt);
#ifdef __i386__
unsigned long long sumsq_mmx(signed short *in,int cnt);
unsigned long long sumsq_sse(signed short *in,int cnt);
unsigned long long sumsq_sse2(signed short *in,int cnt);
#endif
#ifdef __VEC__
unsigned long long sumsq_av(signed short *in,int cnt);
#endif
/* Low-level data structures and routines */
int cpu_features(void);
#endif /* _FEC_H_ */

33
fec-3.0.1/fixed.h Normal file
View File

@ -0,0 +1,33 @@
/* Stuff specific to the CCSDS (255,223) RS codec
* (255,223) code over GF(256). Note: the conventional basis is still
* used; the dual-basis mappings are performed in [en|de]code_rs_ccsds.c
*
* Copyright 2003 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
typedef unsigned char data_t;
static inline int mod255(int x){
while (x >= 255) {
x -= 255;
x = (x >> 8) + (x & 255);
}
return x;
}
#define MODNN(x) mod255(x)
extern data_t CCSDS_alpha_to[];
extern data_t CCSDS_index_of[];
extern data_t CCSDS_poly[];
#define MM 8
#define NN 255
#define ALPHA_TO CCSDS_alpha_to
#define INDEX_OF CCSDS_index_of
#define GENPOLY CCSDS_poly
#define NROOTS 32
#define FCR 112
#define PRIM 11
#define IPRIM 116
#define PAD pad

39
fec-3.0.1/gen_ccsds.c Normal file
View File

@ -0,0 +1,39 @@
/* Generate tables for CCSDS code
* Copyright 2002 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "char.h"
#include "rs-common.h"
#include "fec.h"
int main(){
struct rs *rs;
int i;
rs = init_rs_char(8,0x187,112,11,32,0); /* CCSDS standard */
assert(rs != NULL);
printf("char CCSDS_alpha_to[] = {");
for(i=0;i<256;i++){
if((i % 16) == 0)
printf("\n");
printf("0x%02x,",rs->alpha_to[i]);
}
printf("\n};\n\nchar CCSDS_index_of[] = {");
for(i=0;i<256;i++){
if((i % 16) == 0)
printf("\n");
printf("%3d,",rs->index_of[i]);
}
printf("\n};\n\nchar CCSDS_poly[] = {");
for(i=0;i<33;i++){
if((i % 16) == 0)
printf("\n");
printf("%3d,",rs->genpoly[i]);
}
printf("\n};\n");
exit(0);
}

106
fec-3.0.1/init_rs.h Normal file
View File

@ -0,0 +1,106 @@
/* Common code for intializing a Reed-Solomon control block (char or int symbols)
* Copyright 2004 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#undef NULL
#define NULL ((void *)0)
{
int i, j, sr,root,iprim;
rs = NULL;
/* Check parameter ranges */
if(symsize < 0 || symsize > 8*sizeof(data_t)){
goto done;
}
if(fcr < 0 || fcr >= (1<<symsize))
goto done;
if(prim <= 0 || prim >= (1<<symsize))
goto done;
if(nroots < 0 || nroots >= (1<<symsize))
goto done; /* Can't have more roots than symbol values! */
if(pad < 0 || pad >= ((1<<symsize) -1 - nroots))
goto done; /* Too much padding */
rs = (struct rs *)calloc(1,sizeof(struct rs));
if(rs == NULL)
goto done;
rs->mm = symsize;
rs->nn = (1<<symsize)-1;
rs->pad = pad;
rs->alpha_to = (data_t *)malloc(sizeof(data_t)*(rs->nn+1));
if(rs->alpha_to == NULL){
free(rs);
rs = NULL;
goto done;
}
rs->index_of = (data_t *)malloc(sizeof(data_t)*(rs->nn+1));
if(rs->index_of == NULL){
free(rs->alpha_to);
free(rs);
rs = NULL;
goto done;
}
/* Generate Galois field lookup tables */
rs->index_of[0] = A0; /* log(zero) = -inf */
rs->alpha_to[A0] = 0; /* alpha**-inf = 0 */
sr = 1;
for(i=0;i<rs->nn;i++){
rs->index_of[sr] = i;
rs->alpha_to[i] = sr;
sr <<= 1;
if(sr & (1<<symsize))
sr ^= gfpoly;
sr &= rs->nn;
}
if(sr != 1){
/* field generator polynomial is not primitive! */
free(rs->alpha_to);
free(rs->index_of);
free(rs);
rs = NULL;
goto done;
}
/* Form RS code generator polynomial from its roots */
rs->genpoly = (data_t *)malloc(sizeof(data_t)*(nroots+1));
if(rs->genpoly == NULL){
free(rs->alpha_to);
free(rs->index_of);
free(rs);
rs = NULL;
goto done;
}
rs->fcr = fcr;
rs->prim = prim;
rs->nroots = nroots;
/* Find prim-th root of 1, used in decoding */
for(iprim=1;(iprim % prim) != 0;iprim += rs->nn)
;
rs->iprim = iprim / prim;
rs->genpoly[0] = 1;
for (i = 0,root=fcr*prim; i < nroots; i++,root += prim) {
rs->genpoly[i+1] = 1;
/* Multiply rs->genpoly[] by @**(root + x) */
for (j = i; j > 0; j--){
if (rs->genpoly[j] != 0)
rs->genpoly[j] = rs->genpoly[j-1] ^ rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[j]] + root)];
else
rs->genpoly[j] = rs->genpoly[j-1];
}
/* rs->genpoly[0] can never be zero */
rs->genpoly[0] = rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[0]] + root)];
}
/* convert rs->genpoly[] to index form for quicker encoding */
for (i = 0; i <= nroots; i++)
rs->genpoly[i] = rs->index_of[rs->genpoly[i]];
done:;
}

35
fec-3.0.1/init_rs_char.c Normal file
View File

@ -0,0 +1,35 @@
/* Initialize a RS codec
*
* Copyright 2002 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#include <stdlib.h>
#include "char.h"
#include "rs-common.h"
void free_rs_char(void *p){
struct rs *rs = (struct rs *)p;
free(rs->alpha_to);
free(rs->index_of);
free(rs->genpoly);
free(rs);
}
/* Initialize a Reed-Solomon codec
* symsize = symbol size, bits
* gfpoly = Field generator polynomial coefficients
* fcr = first root of RS code generator polynomial, index form
* prim = primitive element to generate polynomial roots
* nroots = RS code generator polynomial degree (number of roots)
* pad = padding bytes at front of shortened block
*/
void *init_rs_char(int symsize,int gfpoly,int fcr,int prim,
int nroots,int pad){
struct rs *rs;
#include "init_rs.h"
return rs;
}

49
fec-3.0.1/main.c Normal file
View File

@ -0,0 +1,49 @@
#include <stdio.h>
#include "fixed.h"
void encode_rs_8(data_t *data, data_t *parity,int pad);
int decode_rs_8(data_t *data, int *eras_pos, int no_eras, int pad);
int dump(char *name,unsigned char *addr,int len)
{
int i,j;
fprintf(stderr,"Dump of %s\n",name);
for(i=0;i<len;i+=16)
{
fprintf(stderr," %04x :",i);
for(j=0;j<16&&(i+j)<len;j++) fprintf(stderr," %02x",addr[i+j]);
for(;j<16;j++) fprintf(stderr," ");
fprintf(stderr," ");
for(j=0;j<16&&(i+j)<len;j++) fprintf(stderr,"%c",addr[i+j]>=' '&&addr[i+j]<0x7f?addr[i+j]:'.');
fprintf(stderr,"\n");
}
return 0;
}
int main(int argc,char **argv)
{
unsigned char in[255];
unsigned char out[255];
srandom(getpid());
int i;
for(i=0;i<223;i++) in[i]=i;
encode_rs_8(&in[0],&in[223],0);
bcopy(in,out,255);
dump("data with parity",out,255);
for(i=0;i<16;i++) out[random()%255]^=0xff;
dump("data with errors added",out,255);
decode_rs_8(out,NULL,0,0);
dump("data after error correction",out,223);
return 0;
}

26
fec-3.0.1/rs-common.h Normal file
View File

@ -0,0 +1,26 @@
/* Stuff common to all the general-purpose Reed-Solomon codecs
* Copyright 2004 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
/* Reed-Solomon codec control block */
struct rs {
int mm; /* Bits per symbol */
int nn; /* Symbols per block (= (1<<mm)-1) */
data_t *alpha_to; /* log lookup table */
data_t *index_of; /* Antilog lookup table */
data_t *genpoly; /* Generator polynomial */
int nroots; /* Number of generator roots = number of parity symbols */
int fcr; /* First consecutive root, index form */
int prim; /* Primitive element, index form */
int iprim; /* prim-th root of 1, index form */
int pad; /* Padding bytes in shortened block */
};
static inline int modnn(struct rs *rs,int x){
while (x >= rs->nn) {
x -= rs->nn;
x = (x >> rs->mm) + (x & rs->nn);
}
return x;
}

View File

@ -69,4 +69,8 @@ SERVAL_SOURCES = \
$(SERVAL_BASE)strlcpy.c \
$(SERVAL_BASE)vomp.c \
$(SERVAL_BASE)vomp_console.c \
$(SERVAL_BASE)xprintf.c
$(SERVAL_BASE)xprintf.c \
$(SERVAL_BASE)fec-3.0.1/ccsds_tables.c \
$(SERVAL_BASE)fec-3.0.1/decode_rs_8.c \
$(SERVAL_BASE)fec-3.0.1/encode_rs_8.c \
$(SERVAL_BASE)fec-3.0.1/init_rs_char.c