Integrate libfec as part of the gr-satnogs

This commit is contained in:
Manolis Surligas 2019-12-10 21:54:06 +02:00
parent 2a03fc66bf
commit cbc733c25e
164 changed files with 6712 additions and 5893 deletions

View File

@ -90,7 +90,8 @@ endif()
########################################################################
find_package(PythonLibs 3)
find_package(Gnuradio "3.8" REQUIRED runtime fft blocks filter analog digital)
find_package(Gnuradio "3.8" REQUIRED
COMPONENTS runtime blocks fft analog filter digital pmt)
include(GrVersion)
include(GrPlatform) #define LIB_SUFFIX
@ -148,50 +149,6 @@ find_package(JsonCpp REQUIRED)
option(INCLUDE_DEBUG_BLOCKS
"Enable/Disable blocks that are used for debugging purposes" ON)
########################################################################
# Find gr-satnogs external build dependencies
########################################################################
########################################################################
# Search for the libfec if it is already installed in the system
# If not, install the internal one.
########################################################################
find_package(Fec)
if(NOT FEC_FOUND)
message(WARNING "libfec is not installed. The internal libfec will be automatically build and install.")
include(ExternalProject)
ExternalProject_Add(FEC_EXTERNAL
SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/libfec
BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/libfec
CMAKE_ARGS "-DCMAKE_C_FLAGS=${CMAKE_C_FLAGS}"
"-DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}"
"-DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX}"
INSTALL_COMMAND ""
)
ExternalProject_Get_Property(FEC_EXTERNAL binary_dir)
add_library(fec SHARED IMPORTED)
set_property(TARGET fec PROPERTY IMPORTED_LOCATION ${install_dir}/libfec.so)
add_dependencies(fec FEC_EXTERNAL)
set(FEC_LIBRARIES "${binary_dir}/libfec.so")
set(FEC_INCLUDE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/libfec")
# Install the header and the library in the standard places
install(FILES
"${FEC_INCLUDE_DIRS}/fec.h"
DESTINATION "include"
)
install(FILES
${FEC_LIBRARIES}
DESTINATION lib${LIB_SUFFIX}
)
else()
add_library(fec INTERFACE)
endif()
########################################################################
# Setup doxygen option
########################################################################

View File

@ -1,25 +0,0 @@
INCLUDE(FindPkgConfig)
PKG_CHECK_MODULES(PC_FEC fec)
FIND_PATH(
FEC_INCLUDE_DIRS
NAMES fec.h
HINTS $ENV{FEC_DIR}/include
${PC_FEC_INCLUDEDIR}
PATHS /usr/local/include
/usr/include
)
FIND_LIBRARY(
FEC_LIBRARIES
NAMES fec
HINTS $ENV{FEC_DIR}/lib
${PC_FEC_LIBDIR}
PATHS /usr/local/lib
/usr/local/lib64
/usr/lib
/usr/lib64
)
INCLUDE(FindPackageHandleStandardArgs)
FIND_PACKAGE_HANDLE_STANDARD_ARGS(FEC DEFAULT_MSG FEC_LIBRARIES FEC_INCLUDE_DIRS)

View File

@ -1,6 +1,7 @@
usr/bin/*
usr/include/*
usr/lib/*/libgnuradio-satnogs.so
usr/lib/*/libgnuradio-satnogs-fec.so
usr/lib/*/cmake/*
usr/lib/python*
usr/share/*

View File

@ -1,2 +1 @@
usr/lib/*/libfec.so
usr/lib/*/lib*.so.*

View File

@ -21,6 +21,8 @@
########################################################################
# Install public header files
########################################################################
add_subdirectory(libfec)
list(APPEND DEBUG_HEADER_FILES
cw_encoder.h
debug_msg_source_raw.h

View File

@ -0,0 +1,26 @@
# Copyright 2011,2012 Free Software Foundation, Inc.
#
# This file was generated by gr_modtool, a tool from the GNU Radio framework
# This file is a part of gr-satnogs
#
# GNU Radio is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3, or (at your option)
# any later version.
#
# GNU Radio 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 General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with GNU Radio; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street,
# Boston, MA 02110-1301, USA.
########################################################################
# Install public header files
########################################################################
install(FILES
fec.h DESTINATION include/satnogs/libfec
)

View File

@ -0,0 +1,419 @@
/* 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_
#include <satnogs/api.h>
#ifdef __cplusplus
extern "C" {
#endif
/* 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
SATNOGS_API void *create_viterbi27(int len);
SATNOGS_API void set_viterbi27_polynomial(int polys[2]);
SATNOGS_API int init_viterbi27(void *vp, int starting_state);
SATNOGS_API int update_viterbi27_blk(void *vp, unsigned char sym[], int npairs);
SATNOGS_API int chainback_viterbi27(void *vp, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi27(void *vp);
#ifdef __VEC__
SATNOGS_API void *create_viterbi27_av(int len);
SATNOGS_API void set_viterbi27_polynomial_av(int polys[2]);
SATNOGS_API int init_viterbi27_av(void *p, int starting_state);
SATNOGS_API int chainback_viterbi27_av(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi27_av(void *p);
SATNOGS_API int update_viterbi27_blk_av(void *p, unsigned char *syms,
int nbits);
#endif
#ifdef __i386__
SATNOGS_API void *create_viterbi27_mmx(int len);
SATNOGS_API void set_viterbi27_polynomial_mmx(int polys[2]);
SATNOGS_API int init_viterbi27_mmx(void *p, int starting_state);
SATNOGS_API int chainback_viterbi27_mmx(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi27_mmx(void *p);
SATNOGS_API int update_viterbi27_blk_mmx(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi27_sse(int len);
SATNOGS_API void set_viterbi27_polynomial_sse(int polys[2]);
SATNOGS_API int init_viterbi27_sse(void *p, int starting_state);
SATNOGS_API int chainback_viterbi27_sse(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi27_sse(void *p);
SATNOGS_API int update_viterbi27_blk_sse(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi27_sse2(int len);
SATNOGS_API void set_viterbi27_polynomial_sse2(int polys[2]);
SATNOGS_API int init_viterbi27_sse2(void *p, int starting_state);
SATNOGS_API int chainback_viterbi27_sse2(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi27_sse2(void *p);
SATNOGS_API int update_viterbi27_blk_sse2(void *p, unsigned char *syms,
int nbits);
#endif
SATNOGS_API void *create_viterbi27_port(int len);
SATNOGS_API void set_viterbi27_polynomial_port(int polys[2]);
SATNOGS_API int init_viterbi27_port(void *p, int starting_state);
SATNOGS_API int chainback_viterbi27_port(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi27_port(void *p);
SATNOGS_API 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
SATNOGS_API void *create_viterbi29(int len);
SATNOGS_API void set_viterbi29_polynomial(int polys[2]);
SATNOGS_API int init_viterbi29(void *vp, int starting_state);
SATNOGS_API int update_viterbi29_blk(void *vp, unsigned char syms[], int nbits);
SATNOGS_API int chainback_viterbi29(void *vp, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi29(void *vp);
#ifdef __VEC__
SATNOGS_API void *create_viterbi29_av(int len);
SATNOGS_API void set_viterbi29_polynomial_av(int polys[2]);
SATNOGS_API int init_viterbi29_av(void *p, int starting_state);
SATNOGS_API int chainback_viterbi29_av(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi29_av(void *p);
SATNOGS_API int update_viterbi29_blk_av(void *p, unsigned char *syms,
int nbits);
#endif
#ifdef __i386__
SATNOGS_API void *create_viterbi29_mmx(int len);
SATNOGS_API void set_viterbi29_polynomial_mmx(int polys[2]);
SATNOGS_API int chainback_viterbi29_mmx(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi29_mmx(void *p);
SATNOGS_API int update_viterbi29_blk_mmx(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi29_sse(int len);
SATNOGS_API void set_viterbi29_polynomial_sse(int polys[2]);
SATNOGS_API int init_viterbi29_sse(void *p, int starting_state);
SATNOGS_API int chainback_viterbi29_sse(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi29_sse(void *p);
SATNOGS_API int update_viterbi29_blk_sse(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi29_sse2(int len);
SATNOGS_API void set_viterbi29_polynomial_sse2(int polys[2]);
SATNOGS_API int init_viterbi29_sse2(void *p, int starting_state);
SATNOGS_API int chainback_viterbi29_sse2(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi29_sse2(void *p);
SATNOGS_API int update_viterbi29_blk_sse2(void *p, unsigned char *syms,
int nbits);
#endif
SATNOGS_API void *create_viterbi29_port(int len);
SATNOGS_API void set_viterbi29_polynomial_port(int polys[2]);
SATNOGS_API int init_viterbi29_port(void *p, int starting_state);
SATNOGS_API int chainback_viterbi29_port(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi29_port(void *p);
SATNOGS_API 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
SATNOGS_API void *create_viterbi39(int len);
SATNOGS_API void set_viterbi39_polynomial(int polys[3]);
SATNOGS_API int init_viterbi39(void *vp, int starting_state);
SATNOGS_API int update_viterbi39_blk(void *vp, unsigned char syms[], int nbits);
SATNOGS_API int chainback_viterbi39(void *vp, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi39(void *vp);
#ifdef __VEC__
SATNOGS_API void *create_viterbi39_av(int len);
SATNOGS_API void set_viterbi39_polynomial_av(int polys[3]);
SATNOGS_API int init_viterbi39_av(void *p, int starting_state);
SATNOGS_API int chainback_viterbi39_av(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi39_av(void *p);
SATNOGS_API int update_viterbi39_blk_av(void *p, unsigned char *syms,
int nbits);
#endif
#ifdef __i386__
SATNOGS_API void *create_viterbi39_mmx(int len);
SATNOGS_API void set_viterbi39_polynomial_mmx(int polys[3]);
SATNOGS_API int init_viterbi39_mmx(void *p, int starting_state);
SATNOGS_API int chainback_viterbi39_mmx(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi39_mmx(void *p);
SATNOGS_API int update_viterbi39_blk_mmx(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi39_sse(int len);
SATNOGS_API void set_viterbi39_polynomial_sse(int polys[3]);
SATNOGS_API int init_viterbi39_sse(void *p, int starting_state);
SATNOGS_API int chainback_viterbi39_sse(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi39_sse(void *p);
SATNOGS_API int update_viterbi39_blk_sse(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi39_sse2(int len);
SATNOGS_API void set_viterbi39_polynomial_sse2(int polys[3]);
SATNOGS_API int init_viterbi39_sse2(void *p, int starting_state);
SATNOGS_API int chainback_viterbi39_sse2(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi39_sse2(void *p);
SATNOGS_API int update_viterbi39_blk_sse2(void *p, unsigned char *syms,
int nbits);
#endif
SATNOGS_API void *create_viterbi39_port(int len);
SATNOGS_API void set_viterbi39_polynomial_port(int polys[3]);
SATNOGS_API int init_viterbi39_port(void *p, int starting_state);
SATNOGS_API int chainback_viterbi39_port(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi39_port(void *p);
SATNOGS_API 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
SATNOGS_API void *create_viterbi615(int len);
SATNOGS_API void set_viterbi615_polynomial(int polys[6]);
SATNOGS_API int init_viterbi615(void *vp, int starting_state);
SATNOGS_API int update_viterbi615_blk(void *vp, unsigned char *syms, int nbits);
SATNOGS_API int chainback_viterbi615(void *vp, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi615(void *vp);
#ifdef __VEC__
SATNOGS_API void *create_viterbi615_av(int len);
SATNOGS_API void set_viterbi615_polynomial_av(int polys[6]);
SATNOGS_API int init_viterbi615_av(void *p, int starting_state);
SATNOGS_API int chainback_viterbi615_av(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi615_av(void *p);
SATNOGS_API int update_viterbi615_blk_av(void *p, unsigned char *syms,
int nbits);
#endif
#ifdef __i386__
SATNOGS_API void *create_viterbi615_mmx(int len);
SATNOGS_API void set_viterbi615_polynomial_mmx(int polys[6]);
SATNOGS_API int init_viterbi615_mmx(void *p, int starting_state);
SATNOGS_API int chainback_viterbi615_mmx(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi615_mmx(void *p);
SATNOGS_API int update_viterbi615_blk_mmx(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi615_sse(int len);
SATNOGS_API void set_viterbi615_polynomial_sse(int polys[6]);
SATNOGS_API int init_viterbi615_sse(void *p, int starting_state);
SATNOGS_API int chainback_viterbi615_sse(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi615_sse(void *p);
SATNOGS_API int update_viterbi615_blk_sse(void *p, unsigned char *syms,
int nbits);
SATNOGS_API void *create_viterbi615_sse2(int len);
SATNOGS_API void set_viterbi615_polynomial_sse2(int polys[6]);
SATNOGS_API int init_viterbi615_sse2(void *p, int starting_state);
SATNOGS_API int chainback_viterbi615_sse2(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi615_sse2(void *p);
SATNOGS_API int update_viterbi615_blk_sse2(void *p, unsigned char *syms,
int nbits);
#endif
SATNOGS_API void *create_viterbi615_port(int len);
SATNOGS_API void set_viterbi615_polynomial_port(int polys[6]);
SATNOGS_API int init_viterbi615_port(void *p, int starting_state);
SATNOGS_API int chainback_viterbi615_port(void *p, unsigned char *data,
unsigned int nbits, unsigned int endstate);
SATNOGS_API void delete_viterbi615_port(void *p);
SATNOGS_API int update_viterbi615_blk_port(void *p, unsigned char *syms,
int nbits);
/* General purpose RS codec, 8-bit symbols */
SATNOGS_API void encode_rs_char(void *rs, unsigned char *data,
unsigned char *parity);
SATNOGS_API int decode_rs_char(void *rs, unsigned char *data, int *eras_pos,
int no_eras);
SATNOGS_API void *init_rs_char(int symsize, int gfpoly,
int fcr, int prim, int nroots,
int pad);
SATNOGS_API void free_rs_char(void *rs);
/* General purpose RS codec, integer symbols */
SATNOGS_API void encode_rs_int(void *rs, int *data, int *parity);
SATNOGS_API int decode_rs_int(void *rs, int *data, int *eras_pos, int no_eras);
SATNOGS_API void *init_rs_int(int symsize, int gfpoly, int fcr,
int prim, int nroots, int pad);
SATNOGS_API void free_rs_int(void *rs);
/* CCSDS standard (255,223) RS codec with conventional (*not* dual-basis)
* symbol representation
*/
SATNOGS_API void encode_rs_8(unsigned char *data, unsigned char *parity,
int pad);
SATNOGS_API 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 */
SATNOGS_API void encode_rs_ccsds(unsigned char *data, unsigned char *parity,
int pad);
SATNOGS_API 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 */
SATNOGS_API extern enum cpu_mode {UNKNOWN = 0, PORT, MMX, SSE, SSE2, ALTIVEC} Cpu_mode;
SATNOGS_API void find_cpu_mode(
void); /* Call this once at startup to set Cpu_mode */
/* Determine parity of argument: 1 = odd, 0 = even */
#if defined(__i386__) || defined(__x86_64__)
static inline int parityb(unsigned char x)
{
__asm__ __volatile__("test %1,%1;setpo %0" : "=q"(x) : "q"(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 */
SATNOGS_API double normal_rand(double mean, double std_dev);
SATNOGS_API unsigned char addnoise(int sym, double amp, double gain,
double offset, int clip);
extern int Bitcnt[];
/* Dot product functions */
SATNOGS_API void *initdp(signed short coeffs[], int len);
SATNOGS_API void freedp(void *dp);
SATNOGS_API long dotprod(void *dp, signed short a[]);
SATNOGS_API void *initdp_port(signed short coeffs[], int len);
SATNOGS_API void freedp_port(void *dp);
SATNOGS_API long dotprod_port(void *dp, signed short a[]);
#ifdef __i386__
SATNOGS_API void *initdp_mmx(signed short coeffs[], int len);
SATNOGS_API void freedp_mmx(void *dp);
SATNOGS_API long dotprod_mmx(void *dp, signed short a[]);
SATNOGS_API void *initdp_sse(signed short coeffs[], int len);
SATNOGS_API void freedp_sse(void *dp);
SATNOGS_API long dotprod_sse(void *dp, signed short a[]);
SATNOGS_API void *initdp_sse2(signed short coeffs[], int len);
SATNOGS_API void freedp_sse2(void *dp);
SATNOGS_API long dotprod_sse2(void *dp, signed short a[]);
#endif
#ifdef __x86_64__
SATNOGS_API void *initdp_sse2(signed short coeffs[], int len);
SATNOGS_API void freedp_sse2(void *dp);
SATNOGS_API long dotprod_sse2(void *dp, signed short a[]);
#endif
#ifdef __VEC__
SATNOGS_API void *initdp_av(signed short coeffs[], int len);
SATNOGS_API void freedp_av(void *dp);
SATNOGS_API long dotprod_av(void *dp, signed short a[]);
#endif
/* Sum of squares - accepts signed shorts, produces unsigned long long */
SATNOGS_API unsigned long long sumsq(signed short *in, int cnt);
SATNOGS_API unsigned long long sumsq_port(signed short *in, int cnt);
#ifdef __i386__
SATNOGS_API unsigned long long sumsq_mmx(signed short *in, int cnt);
SATNOGS_API unsigned long long sumsq_sse(signed short *in, int cnt);
SATNOGS_API unsigned long long sumsq_sse2(signed short *in, int cnt);
#endif
#ifdef __x86_64__
SATNOGS_API unsigned long long sumsq_sse2(signed short *in, int cnt);
#endif
#ifdef __VEC__
SATNOGS_API unsigned long long sumsq_av(signed short *in, int cnt);
#endif
/* Low-level data structures and routines */
SATNOGS_API int cpu_features(void);
#ifdef __cplusplus
}
#endif
#endif /* _FEC_H_ */

View File

@ -23,6 +23,7 @@
#include <cstdint>
#include <cmath>
#include <cstdio>
#include <arpa/inet.h>
namespace gr {
@ -180,7 +181,7 @@ update_crc32(uint32_t crc, const uint8_t *data, size_t len)
0x2D02EF8DL
};
register uint32_t i;
uint32_t i;
for (i = 0; i < len; i++) {
crc = (crc >> 8) ^ crc32_lut[(crc ^ data[i]) & 0xff];
}

View File

@ -23,6 +23,8 @@
########################################################################
include(GrPlatform) #define LIB_SUFFIX
add_subdirectory(libfec)
list(APPEND satnogs_debug_sources
morse_debug_source_impl.cc
debug_msg_source_impl.cc
@ -78,23 +80,26 @@ if(NOT satnogs_sources)
endif(NOT satnogs_sources)
add_library(gnuradio-satnogs SHARED ${satnogs_sources})
add_dependencies(gnuradio-satnogs fec)
add_dependencies(gnuradio-satnogs gnuradio-satnogs-fec)
target_link_libraries(gnuradio-satnogs PUBLIC
target_link_libraries(gnuradio-satnogs
gnuradio-satnogs-fec
gnuradio::gnuradio-runtime
gnuradio::gnuradio-fft
gnuradio::gnuradio-analog
gnuradio::gnuradio-blocks
gnuradio::gnuradio-fft
gnuradio::gnuradio-digital
gnuradio::gnuradio-pmt
${VOLK_LIBRARIES}
${OGGVORBIS_LIBRARIES}
${PNG_LIBRARIES}
${png++_LIBRARIES}
${FEC_LIBRARIES}
${JSONCPP_LIBRARY}
)
target_include_directories(gnuradio-satnogs
PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../include>
PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/>
PUBLIC $<INSTALL_INTERFACE:include>
)
set_target_properties(gnuradio-satnogs PROPERTIES DEFINE_SYMBOL "gnuradio_satnogs_EXPORTS")
@ -128,6 +133,7 @@ include(GrTest)
list(APPEND test_satnogs_sources
qa_golay24.cc
)
# Anything we need to link to for the unit tests go here
list(APPEND GR_TEST_TARGET_DEPS gnuradio-satnogs)

View File

@ -26,10 +26,7 @@
#include <satnogs/amsat_duv_decoder.h>
#include <satnogs/metadata.h>
#include <gnuradio/blocks/count_bits.h>
extern "C" {
#include <fec.h>
}
#include <satnogs/libfec/fec.h>
namespace gr {
namespace satnogs {

View File

@ -28,10 +28,7 @@
#include <satnogs/metadata.h>
#include <satnogs/utils.h>
#include <satnogs/log.h>
extern "C" {
#include <fec.h>
}
#include <satnogs/libfec/fec.h>
namespace gr {
namespace satnogs {

View File

@ -23,9 +23,7 @@
#include "config.h"
#endif
#include <gnuradio/io_signature.h>
#include <satnogs/golay24.h>
#include <satnogs/utils.h>
namespace gr {

View File

@ -1,42 +1,7 @@
########################################################################
# Project setup
########################################################################
cmake_minimum_required(VERSION 2.8)
project(libfec ASM C)
option(BUILD_32BIT_ON_64BIT "Build a 32-bit library on a 64-bit system" OFF)
# Select the release build type by default to get optimization flags
if(NOT CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE "Release")
message(STATUS "Build type not specified: defaulting to release.")
endif(NOT CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE ${CMAKE_BUILD_TYPE} CACHE STRING "")
list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake/Modules)
if(NOT LIB_INSTALL_DIR)
set(LIB_INSTALL_DIR lib)
endif()
########################################################################
# Version information
########################################################################
set(VERSION_INFO_MAJOR 3)
set(VERSION_INFO_MINOR 0)
set(VERSION_INFO_PATCH 0)
if(NOT DEFINED VERSION_INFO_EXTRA)
set(VERSION_INFO_EXTRA "git")
endif()
include(Version)
if(NOT DEFINED VERSION)
#set(VERSION "\"${VERSION_INFO_MAJOR}.${VERSION_INFO_MINOR}.${VERSION_INFO_PATCH}\"")
set(VERSION "\"${VERSION_INFO}\"")
endif()
include_directories(
${PROJECT_SOURCE_DIR}/include
${PROJECT_SOURCE_DIR}/include/satnogs/libfec
)
########################################################################
# Compiler specific setup
@ -125,9 +90,9 @@ if(TARGET_ARCH MATCHES "x64")
sumsq.c
sumsq_port.c
cpu_mode_x86_64.c
##asm
#sse2bfly27-64.s
#sse2bfly29-64.s
##asm
#sse2bfly27-64.s
#sse2bfly29-64.s
)
elseif(TARGET_ARCH MATCHES "x86")
@ -154,24 +119,24 @@ elseif(TARGET_ARCH MATCHES "x86")
sumsq_sse2.c
sumsq_mmx.c
cpu_mode_x86.c
#asm
cpu_features.s
dotprod_mmx_assist.s
dotprod_sse2_assist.s
mmxbfly27.s
mmxbfly29.s
peak_mmx_assist.s
peak_sse2_assist.s
peak_sse_assist.s
peakval_mmx_assist.s
peakval_sse2_assist.s
peakval_sse_assist.s
sse2bfly27.s
sse2bfly29.s
ssebfly27.s
ssebfly29.s
sumsq_mmx_assist.s
sumsq_sse2_assist.s
#asm
cpu_features.s
dotprod_mmx_assist.s
dotprod_sse2_assist.s
mmxbfly27.s
mmxbfly29.s
peak_mmx_assist.s
peak_sse2_assist.s
peak_sse_assist.s
peakval_mmx_assist.s
peakval_sse2_assist.s
peakval_sse_assist.s
sse2bfly27.s
sse2bfly29.s
ssebfly27.s
ssebfly29.s
sumsq_mmx_assist.s
sumsq_sse2_assist.s
)
elseif(TARGET_ARCH MATCHES "ppc|ppc64")
@ -225,47 +190,6 @@ list(APPEND libfec_sources
ccsds_tal.c
)
################################################################################
# Generate pkg-config file
################################################################################
foreach(inc ${LIBFEC_INCLUDE_DIR})
list(APPEND LIBFEC_PC_CFLAGS "-I${inc}")
endforeach()
foreach(lib ${LIBFEC_LIBRARY_DIRS})
list(APPEND LIBFEC_PC_PRIV_LIBS "-L${lib}")
endforeach()
set(LIBFEC_PC_PREFIX ${CMAKE_INSTALL_PREFIX})
set(LIBFEC_PC_EXEC_PREFIX \${prefix})
set(LIBFEC_PC_LIBDIR \${exec_prefix}/${LIB_INSTALL_DIR})
set(LIBFEC_PC_INCLUDEDIR \${prefix}/include)
set(LIBFEC_PC_VERSION ${VERSION})
set(LIBFEC_PC_LIBS "-lfec")
# Use space-delimiter in the .pc file, rather than CMake's semicolon separator
string(REPLACE ";" " " LIBFEC_PC_CFLAGS "${LIBFEC_PC_CFLAGS}")
string(REPLACE ";" " " LIBFEC_PC_LIBS "${LIBFEC_PC_LIBS}")
# Unset these to avoid hard-coded paths in a cross-environment
if(CMAKE_CROSSCOMPILING)
unset(LIBFEC_PC_CFLAGS)
unset(LIBFEC_PC_LIBS)
endif()
configure_file(
${CMAKE_CURRENT_SOURCE_DIR}/libfec.pc.in
${CMAKE_CURRENT_BINARY_DIR}/libfec.pc
@ONLY
)
install(
FILES ${CMAKE_CURRENT_BINARY_DIR}/libfec.pc
DESTINATION ${LIB_INSTALL_DIR}/pkgconfig/
)
########################################################################
# Setup libraries
########################################################################
@ -273,51 +197,41 @@ install(
# generate ccsds_tab.c
add_executable(gen_ccsds gen_ccsds.c init_rs_char.c)
add_custom_command(
OUTPUT ${CMAKE_BINARY_DIR}/ccsds_tab.c
COMMAND ${CMAKE_BINARY_DIR}/gen_ccsds > ccsds_tab.c
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/ccsds_tab.c
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/gen_ccsds > ccsds_tab.c
DEPENDS gen_ccsds
)
# generate ccsds_tal.c
add_executable(gen_ccsds_tal gen_ccsds_tal.c)
add_custom_command(
OUTPUT ${CMAKE_BINARY_DIR}/ccsds_tal.c
COMMAND ${CMAKE_BINARY_DIR}/gen_ccsds_tal > ccsds_tal.c
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/ccsds_tal.c
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/gen_ccsds_tal > ccsds_tal.c
DEPENDS gen_ccsds_tal
)
# libfec
add_library(libfec_shared SHARED ${libfec_sources})
set_target_properties(libfec_shared PROPERTIES OUTPUT_NAME fec)
target_link_libraries(libfec_shared ${M_LIB})
add_library(gnuradio-satnogs-fec SHARED ${libfec_sources})
target_link_libraries(gnuradio-satnogs-fec ${M_LIB})
install(TARGETS libfec_shared
DESTINATION ${LIB_INSTALL_DIR})
install(FILES "${PROJECT_SOURCE_DIR}/fec.h"
DESTINATION include)
target_include_directories(gnuradio-satnogs-fec
PUBLIC
$<INSTALL_INTERFACE:${CMAKE_INSTALL_PREFIX}/include/satnogs/libfec>
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../../include/satnogs/libfec>
PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/
)
if(APPLE)
set_target_properties(gnuradio-satnogs-fec PROPERTIES
INSTALL_NAME_DIR "${CMAKE_INSTALL_PREFIX}/lib"
)
endif(APPLE)
########################################################################
# Create uninstall target
# Install built library files
########################################################################
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_uninstall.cmake.in"
"${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake"
IMMEDIATE @ONLY)
add_custom_target(uninstall
COMMAND ${CMAKE_COMMAND} -P ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake)
include(GrMiscUtils)
GR_LIBRARY_FOO(gnuradio-satnogs-fec)
########################################################################
# Print Summary
########################################################################
message(STATUS "")
message(STATUS "##########################################################")
message(STATUS "## Building for version: ${VERSION}")
message(STATUS "## Target Architecture: ${TARGET_ARCH}")
message(STATUS "## Using install prefix: ${CMAKE_INSTALL_PREFIX}")
message(STATUS "##########################################################")
message(STATUS "")

View File

@ -1,5 +1,5 @@
typedef unsigned char data_t;
extern unsigned char Taltab[],Tal1tab[];
extern unsigned char Taltab[], Tal1tab[];
#define NN 255
#define NROOTS 32

View File

@ -9,7 +9,7 @@ typedef unsigned char data_t;
#define MM (rs->mm)
#define NN (rs->nn)
#define ALPHA_TO (rs->alpha_to)
#define ALPHA_TO (rs->alpha_to)
#define INDEX_OF (rs->index_of)
#define GENPOLY (rs->genpoly)
#define NROOTS (rs->nroots)

View File

@ -8,6 +8,7 @@
enum cpu_mode Cpu_mode;
// Use the portable code for this unknown CPU
void find_cpu_mode(void) {
void find_cpu_mode(void)
{
Cpu_mode = PORT;
}

45
lib/libfec/cpu_mode_ppc.c Normal file
View File

@ -0,0 +1,45 @@
/* Determine CPU support for SIMD on Power PC
* Copyright 2004 Phil Karn, KA9Q
*/
#include <stdio.h>
#include "fec.h"
#ifdef __VEC__
#include <sys/sysctl.h>
#endif
/* Various SIMD instruction set names */
char *Cpu_modes[] = {"Unknown", "Portable C", "x86 Multi Media Extensions (MMX)",
"x86 Streaming SIMD Extensions (SSE)",
"x86 Streaming SIMD Extensions 2 (SSE2)",
"PowerPC G4/G5 Altivec/Velocity Engine"
};
enum cpu_mode Cpu_mode;
void find_cpu_mode(void)
{
if (Cpu_mode != UNKNOWN) {
return;
}
#ifdef __VEC__
{
/* Ask the OS if we have Altivec support */
int selectors[2] = { CTL_HW, HW_VECTORUNIT };
int hasVectorUnit = 0;
size_t length = sizeof(hasVectorUnit);
int error = sysctl(selectors, 2, &hasVectorUnit, &length, NULL, 0);
if (0 == error && hasVectorUnit) {
Cpu_mode = ALTIVEC;
}
else {
Cpu_mode = PORT;
}
}
#else
Cpu_mode = PORT;
#endif
fprintf(stderr, "SIMD CPU detect: %s\n", Cpu_modes[Cpu_mode]);
}

39
lib/libfec/cpu_mode_x86.c Normal file
View File

@ -0,0 +1,39 @@
/* Determine CPU support for SIMD
* Copyright 2004 Phil Karn, KA9Q
*/
#include <stdio.h>
#include "fec.h"
/* Various SIMD instruction set names */
char *Cpu_modes[] = {"Unknown", "Portable C", "x86 Multi Media Extensions (MMX)",
"x86 Streaming SIMD Extensions (SSE)",
"x86 Streaming SIMD Extensions 2 (SSE2)",
"PowerPC G4/G5 Altivec/Velocity Engine"
};
enum cpu_mode Cpu_mode;
void find_cpu_mode(void)
{
int f;
if (Cpu_mode != UNKNOWN) {
return;
}
/* Figure out what kind of CPU we have */
f = cpu_features();
if (f & (1 << 26)) { /* SSE2 is present */
Cpu_mode = SSE2;
}
else if (f & (1 << 25)) { /* SSE is present */
Cpu_mode = SSE;
}
else if (f & (1 << 23)) { /* MMX is present */
Cpu_mode = MMX;
}
else { /* No SIMD at all */
Cpu_mode = PORT;
}
fprintf(stderr, "SIMD CPU detect: %s\n", Cpu_modes[Cpu_mode]);
}

View File

@ -0,0 +1,30 @@
/* Determine CPU support for SIMD
* Copyright 2004 Phil Karn, KA9Q
*
* Modified in 2012 by Matthias P. Braendli, HB9EGM
*/
#include <stdio.h>
#include "fec.h"
/* Various SIMD instruction set names */
char *Cpu_modes[] = {"Unknown", "Portable C", "x86 Multi Media Extensions (MMX)",
"x86 Streaming SIMD Extensions (SSE)",
"x86 Streaming SIMD Extensions 2 (SSE2)",
"PowerPC G4/G5 Altivec/Velocity Engine"
};
enum cpu_mode Cpu_mode;
void find_cpu_mode(void)
{
int f;
if (Cpu_mode != UNKNOWN) {
return;
}
/* According to the wikipedia entry x86-64, all x86-64 processors have SSE2 */
/* The same assumption is also in other source files ! */
Cpu_mode = SSE2;
fprintf(stderr, "CPU: x86-64, using portable C implementation\n");
}

View File

@ -10,7 +10,7 @@
#include <string.h>
#define NULL ((void *)0)
#define min(a,b) ((a) < (b) ? (a) : (b))
#define min(a,b) ((a) < (b) ? (a) : (b))
#ifdef FIXED
#include "fixed.h"
@ -22,43 +22,48 @@
int DECODE_RS(
#ifdef FIXED
data_t *data, int *eras_pos, int no_eras,int pad){
data_t *data, int *eras_pos, int no_eras, int pad)
{
#else
void *p,data_t *data, int *eras_pos, int no_eras){
void *p, data_t *data, int *eras_pos, int no_eras)
{
struct rs *rs = (struct rs *)p;
#endif
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 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;
#ifdef FIXED
/* Check pad parameter for validity */
if(pad < 0 || pad >= NN)
if (pad < 0 || pad >= NN) {
return -1;
}
#endif
/* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
for(i=0;i<NROOTS;i++)
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)];
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++){
for (i = 0; i < NROOTS; i++) {
syn_error |= s[i];
s[i] = INDEX_OF[s[i]];
}
@ -70,136 +75,149 @@ void *p,data_t *data, int *eras_pos, int no_eras){
count = 0;
goto finish;
}
memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
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]))];
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)];
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++)
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)) {
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;
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);
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++)
for (i = 0; i < count; i++) {
printf("%d ", loc[i]);
}
printf("\n");
#endif
#endif
}
for(i=0;i<NROOTS+1;i++)
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 */
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])];
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 */
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]));
memmove(&b[1], b, NROOTS * sizeof(b[0]));
b[0] = A0;
} else {
}
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 (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;
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);
}
}
memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
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++){
for (i = 0; i < NROOTS + 1; i++) {
lambda[i] = INDEX_OF[lambda[i]];
if(lambda[i] != A0)
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)) {
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--){
for (j = deg_lambda; j > 0; j--) {
if (reg[j] != A0) {
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
}
}
if (q != 0)
continue; /* Not a root */
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);
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)
if (++count == deg_lambda) {
break;
}
}
if (deg_lambda != count) {
/*
@ -213,12 +231,13 @@ void *p,data_t *data, int *eras_pos, int no_eras){
* 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++){
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])];
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];
}
@ -227,19 +246,21 @@ void *p,data_t *data, int *eras_pos, int no_eras){
* 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--) {
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])];
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])];
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) {
@ -250,13 +271,15 @@ void *p,data_t *data, int *eras_pos, int no_eras){
#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])];
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++)
finish:
if (eras_pos != NULL) {
for (i = 0; i < count; i++) {
eras_pos[i] = loc[i];
}
}
return count;
}

View File

@ -64,180 +64,205 @@
#endif
#undef MIN
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#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 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++)
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)];
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++){
for (i = 0; i < NROOTS; i++)
{
syn_error |= s[i];
s[i] = INDEX_OF[s[i]];
}
if (!syn_error) {
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]));
memset(&lambda[1], 0, NROOTS * sizeof(lambda[0]));
lambda[0] = 1;
if (no_eras > 0) {
if (no_eras > 0)
{
/* Init lambda to be the erasure locator polynomial */
lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
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)];
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++)
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)) {
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;
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);
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++)
for (i = 0; i < count; i++) {
printf("%d ", loc[i]);
}
printf("\n");
#endif
#endif
}
for(i=0;i<NROOTS+1;i++)
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 */
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])];
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 */
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]));
memmove(&b[1], b, NROOTS * sizeof(b[0]));
b[0] = A0;
} else {
}
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 (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;
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);
}
}
memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
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++){
for (i = 0; i < NROOTS + 1; i++)
{
lambda[i] = INDEX_OF[lambda[i]];
if(lambda[i] != A0)
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)) {
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--){
for (j = deg_lambda; j > 0; j--) {
if (reg[j] != A0) {
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
}
}
if (q != 0)
continue; /* Not a root */
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);
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)
if (++count == deg_lambda) {
break;
}
}
if (deg_lambda != count) {
if (deg_lambda != count)
{
/*
* deg(lambda) unequal to number of roots => uncorrectable
* error detected
@ -249,12 +274,14 @@
* 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++){
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])];
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];
}
@ -263,19 +290,22 @@
* 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--) {
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])];
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])];
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) {
@ -286,13 +316,16 @@
#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])];
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++)
finish:
if (eras_pos != NULL)
{
for (i = 0; i < count; i++) {
eras_pos[i] = loc[i];
}
}
retval = count;
}

View File

@ -7,18 +7,20 @@
#include <stdio.h>
#endif
#include <satnogs/api.h>
#include <string.h>
#include "fixed.h"
int decode_rs_8(data_t *data, int *eras_pos, int no_eras, int pad){
SATNOGS_API int decode_rs_8(data_t *data, int *eras_pos, int no_eras, int pad)
{
int retval;
if(pad < 0 || pad > 222){
if (pad < 0 || pad > 222) {
return -1;
}
#include "decode_rs.h"
return retval;
}

View File

@ -7,20 +7,25 @@
#include "ccsds.h"
#include "fec.h"
int decode_rs_ccsds(data_t *data,int *eras_pos,int no_eras,int pad){
int i,r;
#include <satnogs/api.h>
SATNOGS_API int decode_rs_ccsds(data_t *data, int *eras_pos, int no_eras, int pad)
{
int i, r;
data_t cdata[NN];
/* Convert data from dual basis to conventional */
for(i=0;i<NN-pad;i++)
for (i = 0; i < NN - pad; i++) {
cdata[i] = Tal1tab[data[i]];
}
r = decode_rs_8(cdata,eras_pos,no_eras,pad);
r = decode_rs_8(cdata, eras_pos, no_eras, pad);
if(r > 0){
if (r > 0) {
/* Convert from conventional to dual basis */
for(i=0;i<NN-pad;i++)
for (i = 0; i < NN - pad; i++) {
data[i] = Taltab[cdata[i]];
}
}
return r;
}

View File

@ -7,16 +7,19 @@
#include <stdio.h>
#endif
#include <satnogs/api.h>
#include <string.h>
#include "char.h"
#include "rs-common.h"
int decode_rs_char(void *p, data_t *data, int *eras_pos, int no_eras){
SATNOGS_API int decode_rs_char(void *p, data_t *data, int *eras_pos, int no_eras)
{
int retval;
struct rs *rs = (struct rs *)p;
#include "decode_rs.h"
return retval;
}

View File

@ -6,17 +6,18 @@
#ifdef DEBUG
#include <stdio.h>
#endif
#include <satnogs/api.h>
#include <string.h>
#include "int.h"
#include "rs-common.h"
int decode_rs_int(void *p, data_t *data, int *eras_pos, int no_eras){
SATNOGS_API int decode_rs_int(void *p, data_t *data, int *eras_pos, int no_eras)
{
int retval;
struct rs *rs = (struct rs *)p;
#include "decode_rs.h"
return retval;
}

View File

@ -6,57 +6,59 @@
#include <stdlib.h>
#include "fec.h"
void *initdp_port(signed short coeffs[],int len);
long dotprod_port(void *p,signed short *b);
void *initdp_port(signed short coeffs[], int len);
long dotprod_port(void *p, signed short *b);
void freedp_port(void *p);
#ifdef __i386__
void *initdp_mmx(signed short coeffs[],int len);
void *initdp_sse2(signed short coeffs[],int len);
long dotprod_mmx(void *p,signed short *b);
long dotprod_sse2(void *p,signed short *b);
void *initdp_mmx(signed short coeffs[], int len);
void *initdp_sse2(signed short coeffs[], int len);
long dotprod_mmx(void *p, signed short *b);
long dotprod_sse2(void *p, signed short *b);
void freedp_mmx(void *p);
void freedp_sse2(void *p);
#endif
#ifdef __VEC__
void *initdp_av(signed short coeffs[],int len);
long dotprod_av(void *p,signed short *b);
void *initdp_av(signed short coeffs[], int len);
long dotprod_av(void *p, signed short *b);
void freedp_av(void *p);
#endif
/* Create and return a descriptor for use with the dot product function */
void *initdp(signed short coeffs[],int len){
void *initdp(signed short coeffs[], int len)
{
find_cpu_mode();
switch(Cpu_mode){
switch (Cpu_mode) {
case PORT:
default:
return initdp_port(coeffs,len);
return initdp_port(coeffs, len);
#ifdef __i386__
case MMX:
case SSE:
return initdp_mmx(coeffs,len);
return initdp_mmx(coeffs, len);
case SSE2:
return initdp_sse2(coeffs,len);
return initdp_sse2(coeffs, len);
#endif
#ifdef __x86_64__
case SSE2:
return initdp_port(coeffs,len);
return initdp_port(coeffs, len);
#endif
#ifdef __VEC__
case ALTIVEC:
return initdp_av(coeffs,len);
return initdp_av(coeffs, len);
#endif
}
}
/* Free a dot product descriptor created earlier */
void freedp(void *p){
switch(Cpu_mode){
void freedp(void *p)
{
switch (Cpu_mode) {
case PORT:
default:
return freedp_port(p);
@ -83,27 +85,28 @@ void freedp(void *p){
/* Compute a dot product given a descriptor and an input array
* The length is taken from the descriptor
*/
long dotprod(void *p,signed short a[]){
switch(Cpu_mode){
long dotprod(void *p, signed short a[])
{
switch (Cpu_mode) {
case PORT:
default:
return dotprod_port(p,a);
return dotprod_port(p, a);
#ifdef __i386__
case MMX:
case SSE:
return dotprod_mmx(p,a);
return dotprod_mmx(p, a);
case SSE2:
return dotprod_sse2(p,a);
return dotprod_sse2(p, a);
#endif
#ifdef __x86_64__
case SSE2:
return dotprod_port(p,a);
return dotprod_port(p, a);
#endif
#ifdef __VEC__
case ALTIVEC:
return dotprod_av(p,a);
return dotprod_av(p, a);
#endif
}
}

View File

@ -16,77 +16,86 @@ struct dotprod {
};
/* Create and return a descriptor for use with the dot product function */
void *initdp_av(signed short coeffs[],int len){
void *initdp_av(signed short coeffs[], int len)
{
struct dotprod *dp;
int i,j;
int i, j;
if(len == 0)
if (len == 0) {
return NULL;
}
dp = (struct dotprod *)calloc(1,sizeof(struct dotprod));
dp = (struct dotprod *)calloc(1, sizeof(struct dotprod));
dp->len = len;
/* Make 8 copies of coefficients, one for each data alignment,
* each aligned to 16-byte boundary
*/
for(i=0;i<8;i++){
dp->coeffs[i] = calloc(1+(len+i-1)/8,sizeof(vector signed short));
for(j=0;j<len;j++)
dp->coeffs[i][j+i] = coeffs[j];
for (i = 0; i < 8; i++) {
dp->coeffs[i] = calloc(1 + (len + i - 1) / 8, sizeof(vector signed short));
for (j = 0; j < len; j++) {
dp->coeffs[i][j + i] = coeffs[j];
}
}
return (void *)dp;
}
/* Free a dot product descriptor created earlier */
void freedp_av(void *p){
void freedp_av(void *p)
{
struct dotprod *dp = (struct dotprod *)p;
int i;
for(i=0;i<8;i++)
if(dp->coeffs[i] != NULL)
for (i = 0; i < 8; i++)
if (dp->coeffs[i] != NULL) {
free(dp->coeffs[i]);
}
free(dp);
}
/* Compute a dot product given a descriptor and an input array
* The length is taken from the descriptor
*/
long dotprod_av(void *p,signed short a[]){
long dotprod_av(void *p, signed short a[])
{
struct dotprod *dp = (struct dotprod *)p;
int al;
vector signed short *ar,*d;
vector signed int sums0,sums1,sums2,sums3;
union { vector signed int v; signed int w[4];} s;
vector signed short *ar, *d;
vector signed int sums0, sums1, sums2, sums3;
union {
vector signed int v;
signed int w[4];
} s;
int nblocks;
/* round ar down to beginning of 16-byte block containing 0th element of
* input buffer. Then set d to one of 8 sets of shifted coefficients
*/
ar = (vector signed short *)((int)a & ~15);
al = ((int)a & 15)/sizeof(signed short);
al = ((int)a & 15) / sizeof(signed short);
d = (vector signed short *)dp->coeffs[al];
nblocks = (dp->len+al-1)/8+1;
nblocks = (dp->len + al - 1) / 8 + 1;
/* Sum into four vectors each holding four 32-bit partial sums */
sums3 = sums2 = sums1 = sums0 = (vector signed int)(0);
while(nblocks >= 4){
sums0 = vec_msums(ar[nblocks-1],d[nblocks-1],sums0);
sums1 = vec_msums(ar[nblocks-2],d[nblocks-2],sums1);
sums2 = vec_msums(ar[nblocks-3],d[nblocks-3],sums2);
sums3 = vec_msums(ar[nblocks-4],d[nblocks-4],sums3);
while (nblocks >= 4) {
sums0 = vec_msums(ar[nblocks - 1], d[nblocks - 1], sums0);
sums1 = vec_msums(ar[nblocks - 2], d[nblocks - 2], sums1);
sums2 = vec_msums(ar[nblocks - 3], d[nblocks - 3], sums2);
sums3 = vec_msums(ar[nblocks - 4], d[nblocks - 4], sums3);
nblocks -= 4;
}
sums0 = vec_adds(sums0,sums1);
sums2 = vec_adds(sums2,sums3);
sums0 = vec_adds(sums0,sums2);
while(nblocks-- > 0){
sums0 = vec_msums(ar[nblocks],d[nblocks],sums0);
sums0 = vec_adds(sums0, sums1);
sums2 = vec_adds(sums2, sums3);
sums0 = vec_adds(sums0, sums2);
while (nblocks-- > 0) {
sums0 = vec_msums(ar[nblocks], d[nblocks], sums0);
}
/* Sum 4 partial sums into final result */
s.v = vec_sums(sums0,(vector signed int)(0));
s.v = vec_sums(sums0, (vector signed int)(0));
return s.w[3];
}

View File

@ -16,50 +16,56 @@ struct dotprod {
*/
signed short *coeffs[4];
};
long dotprod_mmx_assist(signed short *a,signed short *b,int cnt);
long dotprod_mmx_assist(signed short *a, signed short *b, int cnt);
/* Create and return a descriptor for use with the dot product function */
void *initdp_mmx(signed short coeffs[],int len){
void *initdp_mmx(signed short coeffs[], int len)
{
struct dotprod *dp;
int i,j;
int i, j;
if(len == 0)
if (len == 0) {
return NULL;
}
dp = (struct dotprod *)calloc(1,sizeof(struct dotprod));
dp = (struct dotprod *)calloc(1, sizeof(struct dotprod));
dp->len = len;
/* Make 4 copies of coefficients, one for each data alignment */
for(i=0;i<4;i++){
dp->coeffs[i] = (signed short *)calloc(1+(len+i-1)/4,
4*sizeof(signed short));
for(j=0;j<len;j++)
dp->coeffs[i][j+i] = coeffs[j];
for (i = 0; i < 4; i++) {
dp->coeffs[i] = (signed short *)calloc(1 + (len + i - 1) / 4,
4 * sizeof(signed short));
for (j = 0; j < len; j++) {
dp->coeffs[i][j + i] = coeffs[j];
}
}
return (void *)dp;
}
/* Free a dot product descriptor created earlier */
void freedp_mmx(void *p){
void freedp_mmx(void *p)
{
struct dotprod *dp = (struct dotprod *)p;
int i;
for(i=0;i<4;i++)
if(dp->coeffs[i] != NULL)
for (i = 0; i < 4; i++)
if (dp->coeffs[i] != NULL) {
free(dp->coeffs[i]);
}
free(dp);
}
/* Compute a dot product given a descriptor and an input array
* The length is taken from the descriptor
*/
long dotprod_mmx(void *p,signed short a[]){
long dotprod_mmx(void *p, signed short a[])
{
struct dotprod *dp = (struct dotprod *)p;
int al;
signed short *ar;
/* Round input data address down to 8 byte boundary
* NB: depending on the alignment of a[], memory
* before a[] will be accessed. The contents don't matter since they'll
@ -68,14 +74,14 @@ long dotprod_mmx(void *p,signed short a[]){
* in the x86 machines is done on much larger boundaries
*/
ar = (signed short *)((int)a & ~7);
/* Choose one of 4 sets of pre-shifted coefficients. al is both the
* index into dp->coeffs[] and the number of 0 words padded onto
* that coefficients array for alignment purposes
*/
al = a - ar;
/* Call assembler routine to do the work, passing number of 4-word blocks */
return dotprod_mmx_assist(ar,dp->coeffs[al],(dp->len+al-1)/4+1);
return dotprod_mmx_assist(ar, dp->coeffs[al], (dp->len + al - 1) / 4 + 1);
}

View File

@ -13,43 +13,49 @@ struct dotprod {
};
/* Create and return a descriptor for use with the dot product function */
void *initdp_port(signed short coeffs[],int len){
void *initdp_port(signed short coeffs[], int len)
{
struct dotprod *dp;
int j;
if(len == 0)
if (len == 0) {
return NULL;
}
dp = (struct dotprod *)calloc(1,sizeof(struct dotprod));
dp = (struct dotprod *)calloc(1, sizeof(struct dotprod));
dp->len = len;
/* Just one copy of the coefficients for the C version */
dp->coeffs = (signed short *)calloc(len,sizeof(signed short));
for(j=0;j<len;j++)
dp->coeffs = (signed short *)calloc(len, sizeof(signed short));
for (j = 0; j < len; j++) {
dp->coeffs[j] = coeffs[j];
}
return (void *)dp;
}
/* Free a dot product descriptor created earlier */
void freedp_port(void *p){
void freedp_port(void *p)
{
struct dotprod *dp = (struct dotprod *)p;
if(dp->coeffs != NULL)
free(dp->coeffs);
if (dp->coeffs != NULL) {
free(dp->coeffs);
}
free(dp);
}
/* Compute a dot product given a descriptor and an input array
* The length is taken from the descriptor
*/
long dotprod_port(void *p,signed short a[]){
long dotprod_port(void *p, signed short a[])
{
struct dotprod *dp = (struct dotprod *)p;
long corr;
int i;
corr = 0;
for(i=0;i<dp->len;i++){
for (i = 0; i < dp->len; i++) {
corr += (long)a[i] * dp->coeffs[i];
}
return corr;

View File

@ -18,55 +18,61 @@ struct dotprod {
signed short *coeffs[8];
};
long dotprod_sse2_assist(signed short *a,signed short *b,int cnt);
long dotprod_sse2_assist(signed short *a, signed short *b, int cnt);
/* Create and return a descriptor for use with the dot product function */
void *initdp_sse2(signed short coeffs[],int len){
void *initdp_sse2(signed short coeffs[], int len)
{
struct dotprod *dp;
int i,j,blksize;
int i, j, blksize;
if(len == 0)
if (len == 0) {
return NULL;
}
dp = (struct dotprod *)calloc(1,sizeof(struct dotprod));
dp = (struct dotprod *)calloc(1, sizeof(struct dotprod));
dp->len = len;
/* Make 8 copies of coefficients, one for each data alignment,
* each aligned to 16-byte boundary
*/
for(i=0;i<8;i++){
blksize = (1+(len+i-1)/8) * 8*sizeof(signed short);
posix_memalign((void **)&dp->coeffs[i],16,blksize);
memset(dp->coeffs[i],0,blksize);
for(j=0;j<len;j++)
dp->coeffs[i][j+i] = coeffs[j];
for (i = 0; i < 8; i++) {
blksize = (1 + (len + i - 1) / 8) * 8 * sizeof(signed short);
posix_memalign((void **)&dp->coeffs[i], 16, blksize);
memset(dp->coeffs[i], 0, blksize);
for (j = 0; j < len; j++) {
dp->coeffs[i][j + i] = coeffs[j];
}
}
return (void *)dp;
}
/* Free a dot product descriptor created earlier */
void freedp_sse2(void *p){
void freedp_sse2(void *p)
{
struct dotprod *dp = (struct dotprod *)p;
int i;
for(i=0;i<8;i++)
if(dp->coeffs[i] != NULL)
for (i = 0; i < 8; i++)
if (dp->coeffs[i] != NULL) {
free(dp->coeffs[i]);
}
free(dp);
}
/* Compute a dot product given a descriptor and an input array
* The length is taken from the descriptor
*/
long dotprod_sse2(void *p,signed short a[]){
long dotprod_sse2(void *p, signed short a[])
{
struct dotprod *dp = (struct dotprod *)p;
int al;
signed short *ar;
ar = (signed short *)((int)a & ~15);
al = a - ar;
/* Call assembler routine to do the work, passing number of 8-word blocks */
return dotprod_sse2_assist(ar,dp->coeffs[al],(dp->len+al-1)/8+1);
return dotprod_sse2_assist(ar, dp->coeffs[al], (dp->len + al - 1) / 8 + 1);
}

View File

@ -12,28 +12,29 @@
#if HAVE_GETOPT_LONG
struct option Options[] = {
{"force-altivec",0,NULL,'a'},
{"force-port",0,NULL,'p'},
{"force-mmx",0,NULL,'m'},
{"force-sse",0,NULL,'s'},
{"force-sse2",0,NULL,'t'},
{"trials",0,NULL,'n'},
{"force-altivec", 0, NULL, 'a'},
{"force-port", 0, NULL, 'p'},
{"force-mmx", 0, NULL, 'm'},
{"force-sse", 0, NULL, 's'},
{"force-sse2", 0, NULL, 't'},
{"trials", 0, NULL, 'n'},
{NULL},
};
#endif
int main(int argc,char *argv[]){
int main(int argc, char *argv[])
{
short coeffs[512];
short input[2048];
int trials=1000,d;
int trials = 1000, d;
int errors = 0;
#if HAVE_GETOPT_LONG
while((d = getopt_long(argc,argv,"apmstn:",Options,NULL)) != EOF){
while ((d = getopt_long(argc, argv, "apmstn:", Options, NULL)) != EOF) {
#else
while((d = getopt(argc,argv,"apmstn:")) != EOF){
while ((d = getopt(argc, argv, "apmstn:")) != EOF) {
#endif
switch(d){
switch (d) {
case 'a':
Cpu_mode = ALTIVEC;
break;
@ -55,45 +56,47 @@ int main(int argc,char *argv[]){
}
}
while(trials--){
while (trials--) {
long port_result;
long simd_result;
int ntaps;
int i;
int csum = 0;
int offset;
void *dp_simd,*dp_port;
void *dp_simd, *dp_port;
/* Generate set of coefficients
* limit sum of absolute values to 32767 to avoid overflow
*/
memset(coeffs,0,sizeof(coeffs));
for(i=0;i<512;i++){
memset(coeffs, 0, sizeof(coeffs));
for (i = 0; i < 512; i++) {
double gv;
gv = normal_rand(0.,100.);
if(csum + fabs(gv) > 32767)
break;
gv = normal_rand(0., 100.);
if (csum + fabs(gv) > 32767) {
break;
}
coeffs[i] = gv;
csum += fabs(gv);
}
ntaps = i;
/* Compare results to portable C version for a bunch of random data buffers and offsets */
dp_simd = initdp(coeffs,ntaps);
dp_port = initdp_port(coeffs,ntaps);
for(i=0;i<2048;i++)
dp_simd = initdp(coeffs, ntaps);
dp_port = initdp_port(coeffs, ntaps);
for (i = 0; i < 2048; i++) {
input[i] = random();
}
offset = random() & 511;
simd_result = dotprod(dp_simd,input+offset);
port_result = dotprod_port(dp_port,input+offset);
if(simd_result != port_result){
simd_result = dotprod(dp_simd, input + offset);
port_result = dotprod_port(dp_port, input + offset);
if (simd_result != port_result) {
errors++;
}
}
printf("dtest: %d errors\n",errors);
printf("dtest: %d errors\n", errors);
exit(0);
}

View File

@ -14,9 +14,11 @@
void ENCODE_RS(
#ifdef FIXED
data_t *data, data_t *bb,int pad){
data_t *data, data_t *bb, int pad)
{
#else
void *p,data_t *data, data_t *bb){
void *p, data_t *data, data_t *bb)
{
struct rs *rs = (struct rs *)p;
#endif
int i, j;
@ -24,29 +26,33 @@ void *p,data_t *data, data_t *bb){
#ifdef FIXED
/* Check pad parameter for validity */
if(pad < 0 || pad >= NN)
if (pad < 0 || pad >= NN) {
return;
}
#endif
memset(bb,0,NROOTS*sizeof(data_t));
memset(bb, 0, NROOTS * sizeof(data_t));
for(i=0;i<NN-NROOTS-PAD;i++){
for (i = 0; i < NN - NROOTS - PAD; i++) {
feedback = INDEX_OF[data[i] ^ bb[0]];
if(feedback != A0){ /* feedback term is non-zero */
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++)
bb[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS-j])];
for (j = 1; j < NROOTS; j++) {
bb[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS - j])];
}
}
/* Shift */
memmove(&bb[0],&bb[1],sizeof(data_t)*(NROOTS-1));
if(feedback != A0)
bb[NROOTS-1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])];
else
bb[NROOTS-1] = 0;
memmove(&bb[0], &bb[1], sizeof(data_t) * (NROOTS - 1));
if (feedback != A0) {
bb[NROOTS - 1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])];
}
else {
bb[NROOTS - 1] = 0;
}
}
}

View File

@ -8,7 +8,7 @@
* 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
@ -34,25 +34,29 @@
int i, j;
data_t feedback;
memset(parity,0,NROOTS*sizeof(data_t));
memset(parity, 0, NROOTS * sizeof(data_t));
for(i=0;i<NN-NROOTS-PAD;i++){
for (i = 0; i < NN - NROOTS - PAD; i++)
{
feedback = INDEX_OF[data[i] ^ parity[0]];
if(feedback != A0){ /* feedback term is non-zero */
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])];
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;
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;
}
}
}

View File

@ -9,29 +9,33 @@
#endif
static enum {UNKNOWN=0,MMX,SSE,SSE2,ALTIVEC,PORT} cpu_mode;
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);
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);
static void encode_rs_8_av(data_t *data, data_t *parity, int pad);
#endif
#if __i386__
int cpu_features(void);
#endif
void encode_rs_8(data_t *data, data_t *parity,int pad){
if(cpu_mode == UNKNOWN){
void encode_rs_8(data_t *data, data_t *parity, int pad)
{
if (cpu_mode == UNKNOWN) {
#ifdef __i386__
int f;
/* Figure out what kind of CPU we have */
f = cpu_features();
if(f & (1<<26)){ /* SSE2 is present */
if (f & (1 << 26)) { /* SSE2 is present */
cpu_mode = SSE2;
} else if(f & (1<<25)){ /* SSE is present */
}
else if (f & (1 << 25)) { /* SSE is present */
cpu_mode = SSE;
} else if(f & (1<<23)){ /* MMX is present */
}
else if (f & (1 << 23)) { /* MMX is present */
cpu_mode = MMX;
} else { /* No SIMD at all */
}
else { /* No SIMD at all */
cpu_mode = PORT;
}
#elif __x86_64__
@ -42,18 +46,20 @@ void encode_rs_8(data_t *data, data_t *parity,int pad){
int hasVectorUnit = 0;
size_t length = sizeof(hasVectorUnit);
int error = sysctl(selectors, 2, &hasVectorUnit, &length, NULL, 0);
if(0 == error && hasVectorUnit)
if (0 == error && hasVectorUnit) {
cpu_mode = ALTIVEC;
else
}
else {
cpu_mode = PORT;
}
#else
cpu_mode = PORT;
#endif
}
switch(cpu_mode){
switch (cpu_mode) {
#if __vec__
case ALTIVEC:
encode_rs_8_av(data,parity,pad);
encode_rs_8_av(data, parity, pad);
return;
#endif
@ -68,49 +74,62 @@ void encode_rs_8(data_t *data, data_t *parity,int pad){
#endif
default:
encode_rs_8_c(data,parity,pad);
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);
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];
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;
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;
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);
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.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];
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){
static void encode_rs_8_c(data_t *data, data_t *parity, int pad)
{
#include "encode_rs.h"

75
lib/libfec/encode_rs_av.c Normal file
View File

@ -0,0 +1,75 @@
/* Fast Reed-Solomon encoder for (255,223) CCSDS code on PowerPC G4/G5 using Altivec instructions
* Copyright 2004, Phil Karn KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#include <stdio.h>
#include <string.h>
#include "fixed.h"
/* Lookup table for feedback multiplications
* These are the low half of the coefficients. Since the generator polynomial is
* palindromic, we form it by reversing these on the fly
*/
static union {
vector unsigned char v;
unsigned char c[16];
} table[256];
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);
extern data_t CCSDS_alpha_to[];
extern data_t CCSDS_index_of[];
extern data_t CCSDS_poly[];
void rs_init_av()
{
int i, j;
/* The PowerPC is big-endian, so the low-order byte of each vector contains the highest order term in the polynomial */
for (j = 0; j < 16; j++) {
table[0].c[j] = 0;
for (i = 1; i < 256; i++) {
table[i].c[16 - j - 1] = CCSDS_alpha_to[MODNN(CCSDS_poly[j + 1] +
CCSDS_index_of[i])];
}
}
#if 0
for (i = 0; i < 256; i++) {
printf("table[%3d] = %3vu\n", i, table[i].v);
}
#endif
}
void encode_rs_av(unsigned char *data, unsigned char *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];
}
}

View File

@ -8,17 +8,20 @@
#include "ccsds.h"
#include "fec.h"
void encode_rs_ccsds(data_t *data,data_t *parity,int pad){
void encode_rs_ccsds(data_t *data, data_t *parity, int pad)
{
int i;
data_t cdata[NN-NROOTS];
data_t cdata[NN - NROOTS];
/* Convert data from dual basis to conventional */
for(i=0;i<NN-NROOTS-pad;i++)
for (i = 0; i < NN - NROOTS - pad; i++) {
cdata[i] = Tal1tab[data[i]];
}
encode_rs_8(cdata,parity,pad);
encode_rs_8(cdata, parity, pad);
/* Convert parity from conventional to dual basis */
for(i=0;i<NROOTS;i++)
for (i = 0; i < NROOTS; i++) {
parity[i] = Taltab[parity[i]];
}
}

View File

@ -7,7 +7,8 @@
#include "char.h"
#include "rs-common.h"
void encode_rs_char(void *p,data_t *data, data_t *parity){
void encode_rs_char(void *p, data_t *data, data_t *parity)
{
struct rs *rs = (struct rs *)p;
#include "encode_rs.h"

View File

@ -7,7 +7,8 @@
#include "int.h"
#include "rs-common.h"
void encode_rs_int(void *p,data_t *data, data_t *parity){
void encode_rs_int(void *p, data_t *data, data_t *parity)
{
struct rs *rs = (struct rs *)p;
#include "encode_rs.h"

128
lib/libfec/exercise.c Normal file
View File

@ -0,0 +1,128 @@
/* Exercise an RS codec a specified number of times using random
* data and error patterns
*
* Copyright 2002 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#define FLAG_ERASURE 1 /* Randomly flag 50% of errors as erasures */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef FIXED
#include "fixed.h"
#define EXERCISE exercise_8
#elif defined(CCSDS)
#include "fixed.h"
#include "ccsds.h"
#define EXERCISE exercise_ccsds
#elif defined(BIGSYM)
#include "int.h"
#define EXERCISE exercise_int
#else
#include "char.h"
#define EXERCISE exercise_char
#endif
#ifdef FIXED
#define PRINTPARM printf("(255,223):");
#elif defined(CCSDS)
#define PRINTPARM printf("CCSDS (255,223):");
#else
#define PRINTPARM printf("(%d,%d):",rs->nn,rs->nn-rs->nroots);
#endif
/* Exercise the RS codec passed as an argument */
int EXERCISE(
#if !defined(CCSDS) && !defined(FIXED)
void *p,
#endif
int trials)
{
#if !defined(CCSDS) && !defined(FIXED)
struct rs *rs = (struct rs *)p;
#endif
data_t block[NN], tblock[NN];
int i;
int errors;
int errlocs[NN];
int derrlocs[NROOTS];
int derrors;
int errval, errloc;
int erasures;
int decoder_errors = 0;
while (trials-- != 0) {
/* Test up to the error correction capacity of the code */
for (errors = 0; errors <= NROOTS / 2; errors++) {
/* Load block with random data and encode */
for (i = 0; i < NN - NROOTS; i++) {
block[i] = random() & NN;
}
#if defined(CCSDS) || defined(FIXED)
ENCODE_RS(&block[0], &block[NN - NROOTS], 0);
#else
ENCODE_RS(rs, &block[0], &block[NN - NROOTS]);
#endif
/* Make temp copy, seed with errors */
memcpy(tblock, block, sizeof(tblock));
memset(errlocs, 0, sizeof(errlocs));
memset(derrlocs, 0, sizeof(derrlocs));
erasures = 0;
for (i = 0; i < errors; i++) {
do {
errval = random() & NN;
}
while (errval == 0); /* Error value must be nonzero */
do {
errloc = random() % NN;
}
while (errlocs[errloc] != 0); /* Must not choose the same location twice */
errlocs[errloc] = 1;
#if FLAG_ERASURE
if (random() & 1) { /* 50-50 chance */
derrlocs[erasures++] = errloc;
}
#endif
tblock[errloc] ^= errval;
}
/* Decode the errored block */
#if defined(CCSDS) || defined(FIXED)
derrors = DECODE_RS(tblock, derrlocs, erasures, 0);
#else
derrors = DECODE_RS(rs, tblock, derrlocs, erasures);
#endif
if (derrors != errors) {
PRINTPARM
printf(" decoder says %d errors, true number is %d\n", derrors, errors);
decoder_errors++;
}
for (i = 0; i < derrors; i++) {
if (errlocs[derrlocs[i]] == 0) {
PRINTPARM
printf(" decoder indicates error in location %d without error\n", derrlocs[i]);
decoder_errors++;
}
}
if (memcmp(tblock, block, sizeof(tblock)) != 0) {
PRINTPARM
printf(" uncorrected errors! output ^ input:");
decoder_errors++;
for (i = 0; i < NN; i++) {
printf(" %02x", tblock[i] ^ block[i]);
}
printf("\n");
}
}
}
return decoder_errors;
}

68
lib/libfec/fec.c Normal file
View File

@ -0,0 +1,68 @@
/* Utility routines for FEC support
* Copyright 2004, Phil Karn, KA9Q
*/
#include <stdio.h>
#include "fec.h"
unsigned char Partab[256];
int P_init;
/* Create 256-entry odd-parity lookup table
* Needed only on non-ia32 machines
*/
void partab_init(void)
{
int i, cnt, ti;
/* Initialize parity lookup table */
for (i = 0; i < 256; i++) {
cnt = 0;
ti = i;
while (ti) {
if (ti & 1) {
cnt++;
}
ti >>= 1;
}
Partab[i] = cnt & 1;
}
P_init = 1;
}
/* Lookup table giving count of 1 bits for integers 0-255 */
int Bitcnt[] = {
0, 1, 1, 2, 1, 2, 2, 3,
1, 2, 2, 3, 2, 3, 3, 4,
1, 2, 2, 3, 2, 3, 3, 4,
2, 3, 3, 4, 3, 4, 4, 5,
1, 2, 2, 3, 2, 3, 3, 4,
2, 3, 3, 4, 3, 4, 4, 5,
2, 3, 3, 4, 3, 4, 4, 5,
3, 4, 4, 5, 4, 5, 5, 6,
1, 2, 2, 3, 2, 3, 3, 4,
2, 3, 3, 4, 3, 4, 4, 5,
2, 3, 3, 4, 3, 4, 4, 5,
3, 4, 4, 5, 4, 5, 5, 6,
2, 3, 3, 4, 3, 4, 4, 5,
3, 4, 4, 5, 4, 5, 5, 6,
3, 4, 4, 5, 4, 5, 5, 6,
4, 5, 5, 6, 5, 6, 6, 7,
1, 2, 2, 3, 2, 3, 3, 4,
2, 3, 3, 4, 3, 4, 4, 5,
2, 3, 3, 4, 3, 4, 4, 5,
3, 4, 4, 5, 4, 5, 5, 6,
2, 3, 3, 4, 3, 4, 4, 5,
3, 4, 4, 5, 4, 5, 5, 6,
3, 4, 4, 5, 4, 5, 5, 6,
4, 5, 5, 6, 5, 6, 6, 7,
2, 3, 3, 4, 3, 4, 4, 5,
3, 4, 4, 5, 4, 5, 5, 6,
3, 4, 4, 5, 4, 5, 5, 6,
4, 5, 5, 6, 5, 6, 6, 7,
3, 4, 4, 5, 4, 5, 5, 6,
4, 5, 5, 6, 5, 6, 6, 7,
4, 5, 5, 6, 5, 6, 6, 7,
5, 6, 6, 7, 6, 7, 7, 8,
};

View File

@ -7,7 +7,8 @@
*/
typedef unsigned char data_t;
static inline int mod255(int x){
static inline int mod255(int x)
{
while (x >= 255) {
x -= 255;
x = (x >> 8) + (x & 255);

View File

@ -9,30 +9,34 @@
#include "rs-common.h"
#include "fec.h"
int main(){
int main()
{
struct rs *rs;
int i;
rs = init_rs_char(8,0x187,112,11,32,0); /* CCSDS standard */
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)
for (i = 0; i < 256; i++) {
if ((i % 16) == 0) {
printf("\n");
printf("0x%02x,",rs->alpha_to[i]);
}
printf("0x%02x,", rs->alpha_to[i]);
}
printf("\n};\n\nchar CCSDS_index_of[] = {");
for(i=0;i<256;i++){
if((i % 16) == 0)
for (i = 0; i < 256; i++) {
if ((i % 16) == 0) {
printf("\n");
printf("%3d,",rs->index_of[i]);
}
printf("%3d,", rs->index_of[i]);
}
printf("\n};\n\nchar CCSDS_poly[] = {");
for(i=0;i<33;i++){
if((i % 16) == 0)
for (i = 0; i < 33; i++) {
if ((i % 16) == 0) {
printf("\n");
}
printf("%3d,",rs->genpoly[i]);
printf("%3d,", rs->genpoly[i]);
}
printf("\n};\n");
exit(0);

View File

@ -14,7 +14,7 @@
#include <stdlib.h>
#define DTYPE unsigned char
DTYPE Taltab[256],Tal1tab[256];
DTYPE Taltab[256], Tal1tab[256];
static DTYPE tal[] = { 0x8d, 0xef, 0xec, 0x86, 0xfa, 0x99, 0xaf, 0x7b };
@ -23,29 +23,33 @@ static DTYPE tal[] = { 0x8d, 0xef, 0xec, 0x86, 0xfa, 0x99, 0xaf, 0x7b };
* and Berlekamp's dual basis representation
* (l0, l1, ...l7)
*/
int main(){
int i,j,k;
int main()
{
int i, j, k;
for(i=0;i<256;i++){/* For each value of input */
for (i = 0; i < 256; i++) { /* For each value of input */
Taltab[i] = 0;
for(j=0;j<8;j++) /* for each column of matrix */
for(k=0;k<8;k++){ /* for each row of matrix */
if(i & (1<<k))
Taltab[i] ^= tal[7-k] & (1<<j);
for (j = 0; j < 8; j++) /* for each column of matrix */
for (k = 0; k < 8; k++) { /* for each row of matrix */
if (i & (1 << k)) {
Taltab[i] ^= tal[7 - k] & (1 << j);
}
}
Tal1tab[Taltab[i]] = i;
}
printf("unsigned char Taltab[] = {\n");
for(i=0;i<256;i++){
if((i % 16) == 0)
for (i = 0; i < 256; i++) {
if ((i % 16) == 0) {
printf("\n");
printf("0x%02x,",Taltab[i]);
}
printf("0x%02x,", Taltab[i]);
}
printf("\n};\n\nunsigned char Tal1tab[] = {");
for(i=0;i<256;i++){
if((i % 16) == 0)
for (i = 0; i < 256; i++) {
if ((i % 16) == 0) {
printf("\n");
printf("0x%02x,",Tal1tab[i]);
}
printf("0x%02x,", Tal1tab[i]);
}
printf("\n};\n");
exit(0);

View File

@ -12,7 +12,8 @@
#include "rs-common.h"
void free_rs(void *p){
void free_rs(void *p)
{
struct rs *rs = (struct rs *)p;
free(rs->alpha_to);
@ -29,8 +30,9 @@ void free_rs(void *p){
* nroots = RS code generator polynomial degree (number of roots)
* pad = padding bytes at front of shortened block
*/
void *init_rs_common(int symsize,int gfpoly,int fcr,int prim,
int nroots,int pad){
void *init_rs_common(int symsize, int gfpoly, int fcr, int prim,
int nroots, int pad)
{
struct rs *rs;
#include "init_rs.h"

View File

@ -6,39 +6,52 @@
#define NULL ((void *)0)
{
int i, j, sr,root,iprim;
int i, j, sr, root, iprim;
rs = NULL;
/* Check parameter ranges */
if(symsize < 0 || symsize > 8*sizeof(data_t)){
if (symsize < 0 || symsize > 8 * sizeof(data_t))
{
goto done;
}
if(fcr < 0 || fcr >= (1<<symsize))
if (fcr < 0 || fcr >= (1 << symsize))
{
goto done;
if(prim <= 0 || prim >= (1<<symsize))
}
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 */
}
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)
rs = (struct rs *)calloc(1, sizeof(struct rs));
if (rs == NULL)
{
goto done;
}
rs->mm = symsize;
rs->nn = (1<<symsize)-1;
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){
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){
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;
@ -49,15 +62,18 @@
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++){
for (i = 0; i < rs->nn; i++)
{
rs->index_of[sr] = i;
rs->alpha_to[i] = sr;
sr <<= 1;
if(sr & (1<<symsize))
if (sr & (1 << symsize)) {
sr ^= gfpoly;
}
sr &= rs->nn;
}
if(sr != 1){
if (sr != 1)
{
/* field generator polynomial is not primitive! */
free(rs->alpha_to);
free(rs->index_of);
@ -67,8 +83,9 @@
}
/* Form RS code generator polynomial from its roots */
rs->genpoly = (data_t *)malloc(sizeof(data_t)*(nroots+1));
if(rs->genpoly == NULL){
rs->genpoly = (data_t *)malloc(sizeof(data_t) * (nroots + 1));
if (rs->genpoly == NULL)
{
free(rs->alpha_to);
free(rs->index_of);
free(rs);
@ -80,27 +97,34 @@
rs->nroots = nroots;
/* Find prim-th root of 1, used in decoding */
for(iprim=1;(iprim % prim) != 0;iprim += rs->nn)
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;
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];
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)];
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:;
}
done:
;
}

View File

@ -8,7 +8,8 @@
#include "char.h"
#include "rs-common.h"
void free_rs_char(void *p){
void free_rs_char(void *p)
{
struct rs *rs = (struct rs *)p;
free(rs->alpha_to);
@ -25,8 +26,9 @@ void free_rs_char(void *p){
* 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){
void *init_rs_char(int symsize, int gfpoly, int fcr, int prim,
int nroots, int pad)
{
struct rs *rs;
#include "init_rs.h"

View File

@ -8,7 +8,8 @@
#include "int.h"
#include "rs-common.h"
void free_rs_int(void *p){
void free_rs_int(void *p)
{
struct rs *rs = (struct rs *)p;
free(rs->alpha_to);
@ -25,8 +26,9 @@ void free_rs_int(void *p){
* nroots = RS code generator polynomial degree (number of roots)
* pad = padding bytes at front of shortened block
*/
void *init_rs_int(int symsize,int gfpoly,int fcr,int prim,
int nroots,int pad){
void *init_rs_int(int symsize, int gfpoly, int fcr, int prim,
int nroots, int pad)
{
struct rs *rs;
#include "init_rs.h"

View File

@ -9,7 +9,7 @@ typedef unsigned int data_t;
#define MM (rs->mm)
#define NN (rs->nn)
#define ALPHA_TO (rs->alpha_to)
#define ALPHA_TO (rs->alpha_to)
#define INDEX_OF (rs->index_of)
#define GENPOLY (rs->genpoly)
#define NROOTS (rs->nroots)

View File

@ -9,30 +9,33 @@
#define NSAMP 200002
#define OFFSET 1
int peakval(signed short *,int);
int peakval_port(signed short *,int);
int peakval(signed short *, int);
int peakval_port(signed short *, int);
int main(){
int i,s;
int result,rresult;
int main()
{
int i, s;
int result, rresult;
signed short samples[NSAMP];
srandom(time(NULL));
for(i=0;i<NSAMP;i++){
for (i = 0; i < NSAMP; i++) {
do {
s = random() & 0x0fff;
} while(s == 0x8000);
}
while (s == 0x8000);
samples[i] = s;
}
samples[5] = 25000;
rresult = peakval_port(&samples[OFFSET],NSAMP-OFFSET);
result = peakval(&samples[OFFSET],NSAMP-OFFSET);
if(result == rresult){
rresult = peakval_port(&samples[OFFSET], NSAMP - OFFSET);
result = peakval(&samples[OFFSET], NSAMP - OFFSET);
if (result == rresult) {
printf("OK\n");
} else {
printf("peak mismatch: %d != %d\n",result,rresult);
}
else {
printf("peak mismatch: %d != %d\n", result, rresult);
}
exit(0);
}

51
lib/libfec/peakval.c Normal file
View File

@ -0,0 +1,51 @@
/* Switch to appropriate version of peakval routine
* Copyright 2004, Phil Karn, KA9Q
*/
#include <stdlib.h>
#include "fec.h"
int peakval_port(signed short *b, int cnt);
#ifdef __i386__
int peakval_mmx(signed short *b, int cnt);
int peakval_sse(signed short *b, int cnt);
int peakval_sse2(signed short *b, int cnt);
#endif
#ifdef __x86_64__
int peakval_sse2(signed short *b, int cnt);
#endif
#ifdef __VEC__
int peakval_av(signed short *b, int cnt);
#endif
int peakval(signed short *b, int cnt)
{
find_cpu_mode();
switch (Cpu_mode) {
case PORT:
default:
return peakval_port(b, cnt);
#ifdef __i386__
case MMX:
return peakval_mmx(b, cnt);
case SSE:
return peakval_sse(b, cnt);
case SSE2:
return peakval_sse2(b, cnt);
#endif
#ifdef __x86_64__
case SSE2:
return peakval_port(b, cnt);
//return peakval_sse2(b,cnt);
#endif
#ifdef __VEC__
case ALTIVEC:
return peakval_av(b, cnt);
#endif
}
}

67
lib/libfec/peakval_av.c Normal file
View File

@ -0,0 +1,67 @@
/* Return the largest absolute value of a vector of signed shorts
* This is the Altivec SIMD version.
* Copyright 2004 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#include "fec.h"
signed short peakval_av(signed short *in, int cnt)
{
vector signed short x;
int pad;
union {
vector signed char cv;
vector signed short hv;
signed short s[8];
signed char c[16];
} s;
vector signed short smallest, largest;
smallest = (vector signed short)(0);
largest = (vector signed short)(0);
if ((pad = (int)in & 15) != 0) {
/* Load unaligned leading word */
x = vec_perm(vec_ld(0, in), (vector signed short)(0), vec_lvsl(0, in));
if (cnt < 8) { /* Shift right to chop stuff beyond end of short block */
s.c[15] = (8 - cnt) << 4;
x = vec_sro(x, s.cv);
}
smallest = vec_min(smallest, x);
largest = vec_max(largest, x);
in += 8 - pad / 2;
cnt -= 8 - pad / 2;
}
/* Everything is now aligned, rip through most of the block */
while (cnt >= 8) {
x = vec_ld(0, in);
smallest = vec_min(smallest, x);
largest = vec_max(largest, x);
in += 8;
cnt -= 8;
}
/* Handle trailing fragment, if any */
if (cnt > 0) {
x = vec_ld(0, in);
s.c[15] = (8 - cnt) << 4;
x = vec_sro(x, s.cv);
smallest = vec_min(smallest, x);
largest = vec_max(largest, x);
}
/* Combine and extract result */
largest = vec_max(largest, vec_abs(smallest));
s.c[15] = 64; /* Shift right four 16-bit words */
largest = vec_max(largest, vec_sro(largest, s.cv));
s.c[15] = 32; /* Shift right two 16-bit words */
largest = vec_max(largest, vec_sro(largest, s.cv));
s.c[15] = 16; /* Shift right one 16-bit word */
largest = vec_max(largest, vec_sro(largest, s.cv));
s.hv = largest;
return s.s[7];
}

View File

@ -4,31 +4,35 @@
#include <stdlib.h>
int peakval_mmx_assist(signed short *,int);
int peakval_mmx_assist(signed short *, int);
int peakval_mmx(signed short *b,int cnt){
int peakval_mmx(signed short *b, int cnt)
{
int peak = 0;
int a;
while(((int)b & 7) != 0 && cnt != 0){
while (((int)b & 7) != 0 && cnt != 0) {
a = abs(*b);
if(a > peak)
if (a > peak) {
peak = a;
}
b++;
cnt--;
}
a = peakval_mmx_assist(b,cnt);
if(a > peak)
a = peakval_mmx_assist(b, cnt);
if (a > peak) {
peak = a;
}
b += cnt & ~3;
cnt &= 3;
while(cnt != 0){
while (cnt != 0) {
a = abs(*b);
if(a > peak)
if (a > peak) {
peak = a;
}
b++;
cnt--;
}
return peak;
}
}

View File

@ -3,14 +3,16 @@
*/
#include <stdlib.h>
#include "fec.h"
int peakval_port(signed short *b,int len){
int peakval_port(signed short *b, int len)
{
int peak = 0;
int a,i;
int a, i;
for(i=0;i<len;i++){
for (i = 0; i < len; i++) {
a = abs(b[i]);
if(a > peak)
if (a > peak) {
peak = a;
}
}
return peak;
}

View File

@ -5,31 +5,35 @@
#include <stdlib.h>
#include "fec.h"
int peakval_sse_assist(signed short *,int);
int peakval_sse_assist(signed short *, int);
int peakval_sse(signed short *b,int cnt){
int peakval_sse(signed short *b, int cnt)
{
int peak = 0;
int a;
while(((int)b & 7) != 0 && cnt != 0){
while (((int)b & 7) != 0 && cnt != 0) {
a = abs(*b);
if(a > peak)
if (a > peak) {
peak = a;
}
b++;
cnt--;
}
a = peakval_sse_assist(b,cnt);
if(a > peak)
a = peakval_sse_assist(b, cnt);
if (a > peak) {
peak = a;
}
b += cnt & ~3;
cnt &= 3;
while(cnt != 0){
while (cnt != 0) {
a = abs(*b);
if(a > peak)
if (a > peak) {
peak = a;
}
b++;
cnt--;
}
return peak;
}
}

View File

@ -4,31 +4,35 @@
#include <stdlib.h>
#include "fec.h"
int peakval_sse2_assist(signed short *,int);
int peakval_sse2_assist(signed short *, int);
int peakval_sse2(signed short *b,int cnt){
int peakval_sse2(signed short *b, int cnt)
{
int peak = 0;
int a;
while(((int)b & 15) != 0 && cnt != 0){
while (((int)b & 15) != 0 && cnt != 0) {
a = abs(*b);
if(a > peak)
if (a > peak) {
peak = a;
}
b++;
cnt--;
}
a = peakval_sse2_assist(b,cnt);
if(a > peak)
a = peakval_sse2_assist(b, cnt);
if (a > peak) {
peak = a;
}
b += cnt & ~7;
cnt &= 7;
while(cnt != 0){
while (cnt != 0) {
a = abs(*b);
if(a > peak)
if (a > peak) {
peak = a;
}
b++;
cnt--;
}
return peak;
}
}

View File

@ -17,7 +17,8 @@ struct rs {
int pad; /* Padding bytes in shortened block */
};
static inline int modnn(struct rs *rs,int x){
static inline int modnn(struct rs *rs, int x)
{
while (x >= rs->nn) {
x -= rs->nn;
x = (x >> rs->mm) + (x & rs->nn);

60
lib/libfec/rs_speedtest.c Normal file
View File

@ -0,0 +1,60 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include "fec.h"
int main()
{
unsigned char block[255];
int i;
void *rs;
struct rusage start, finish;
double extime;
int trials = 10000;
for (i = 0; i < 223; i++) {
block[i] = 0x01;
}
rs = init_rs_char(8, 0x187, 112, 11, 32, 0);
encode_rs_char(rs, block, &block[223]);
getrusage(RUSAGE_SELF, &start);
for (i = 0; i < trials; i++) {
#if 0
block[0] ^= 0xff; /* Introduce an error */
block[2] ^= 0xff; /* Introduce an error */
#endif
decode_rs_char(rs, block, NULL, 0);
}
getrusage(RUSAGE_SELF, &finish);
extime = finish.ru_utime.tv_sec - start.ru_utime.tv_sec + 1e-6 *
(finish.ru_utime.tv_usec - start.ru_utime.tv_usec);
printf("Execution time for %d Reed-Solomon blocks using general decoder: %.2f sec\n",
trials, extime);
printf("decoder speed: %g bits/s\n", trials * 223 * 8 / extime);
encode_rs_8(block, &block[223], 0);
getrusage(RUSAGE_SELF, &start);
for (i = 0; i < trials; i++) {
#if 0
block[0] ^= 0xff; /* Introduce an error */
block[2] ^= 0xff; /* Introduce an error */
#endif
decode_rs_8(block, NULL, 0, 0);
}
getrusage(RUSAGE_SELF, &finish);
extime = finish.ru_utime.tv_sec - start.ru_utime.tv_sec + 1e-6 *
(finish.ru_utime.tv_usec - start.ru_utime.tv_usec);
printf("Execution time for %d Reed-Solomon blocks using CCSDS decoder: %.2f sec\n",
trials, extime);
printf("decoder speed: %g bits/s\n", trials * 223 * 8 / extime);
exit(0);
}

324
lib/libfec/rstest.c Normal file
View File

@ -0,0 +1,324 @@
/* Test the Reed-Solomon codecs
* for various block sizes and with random data and random error patterns
*
* 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 <memory.h>
#include <time.h>
#include "fec.h"
struct etab {
int symsize;
int genpoly;
int fcs;
int prim;
int nroots;
int ntrials;
} Tab[] = {
{2, 0x7, 1, 1, 1, 10 },
{3, 0xb, 1, 1, 2, 10 },
{4, 0x13, 1, 1, 4, 10 },
{5, 0x25, 1, 1, 6, 10 },
{6, 0x43, 1, 1, 8, 10 },
{7, 0x89, 1, 1, 10, 10 },
{8, 0x11d, 1, 1, 32, 10 },
{8, 0x187, 112, 11, 32, 10 }, /* Duplicates CCSDS codec */
{9, 0x211, 1, 1, 32, 10 },
{10, 0x409, 1, 1, 32, 10 },
{11, 0x805, 1, 1, 32, 10 },
{12, 0x1053, 1, 1, 32, 5 },
{13, 0x201b, 1, 1, 32, 2 },
{14, 0x4443, 1, 1, 32, 1 },
{15, 0x8003, 1, 1, 32, 1 },
{16, 0x1100b, 1, 1, 32, 1 },
{0, 0, 0, 0, 0},
};
int exercise_char(struct etab *e);
int exercise_int(struct etab *e);
int exercise_8(void);
int main()
{
int i;
srandom(time(NULL));
printf("Testing fixed CCSDS encoder...\n");
exercise_8();
for (i = 0; Tab[i].symsize != 0; i++) {
int nn, kk;
nn = (1 << Tab[i].symsize) - 1;
kk = nn - Tab[i].nroots;
printf("Testing (%d,%d) code...\n", nn, kk);
if (Tab[i].symsize <= 8) {
exercise_char(&Tab[i]);
}
else {
exercise_int(&Tab[i]);
}
}
exit(0);
}
int exercise_8(void)
{
int nn = 255;
unsigned char block[nn], tblock[nn];
int errlocs[nn], derrlocs[nn];
int i;
int errors;
int derrors, kk;
int errval, errloc;
int erasures;
int decoder_errors = 0;
/* Compute code parameters */
kk = 223;
/* Test up to the error correction capacity of the code */
for (errors = 0; errors <= (nn - kk) / 2; errors++) {
/* Load block with random data and encode */
for (i = 0; i < kk; i++) {
block[i] = random() & nn;
}
memcpy(tblock, block, sizeof(block));
encode_rs_8(block, &block[kk], 0);
/* Make temp copy, seed with errors */
memcpy(tblock, block, sizeof(block));
memset(errlocs, 0, sizeof(errlocs));
memset(derrlocs, 0, sizeof(derrlocs));
erasures = 0;
for (i = 0; i < errors; i++) {
do {
errval = random() & nn;
}
while (errval == 0); /* Error value must be nonzero */
do {
errloc = random() % nn;
}
while (errlocs[errloc] != 0); /* Must not choose the same location twice */
errlocs[errloc] = 1;
#if FLAG_ERASURE
if (random() & 1) { /* 50-50 chance */
derrlocs[erasures++] = errloc;
}
#endif
tblock[errloc] ^= errval;
}
/* Decode the errored block */
derrors = decode_rs_8(tblock, derrlocs, erasures, 0);
if (derrors != errors) {
printf("(%d,%d) decoder says %d errors, true number is %d\n", nn, kk, derrors,
errors);
decoder_errors++;
}
for (i = 0; i < derrors; i++) {
if (errlocs[derrlocs[i]] == 0) {
printf("(%d,%d) decoder indicates error in location %d without error\n", nn, kk,
derrlocs[i]);
decoder_errors++;
}
}
if (memcmp(tblock, block, sizeof(tblock)) != 0) {
printf("(%d,%d) decoder uncorrected errors! output ^ input:", nn, kk);
decoder_errors++;
for (i = 0; i < nn; i++) {
printf(" %02x", tblock[i] ^ block[i]);
}
printf("\n");
}
}
return decoder_errors;
}
int exercise_char(struct etab *e)
{
int nn = (1 << e->symsize) - 1;
unsigned char block[nn], tblock[nn];
int errlocs[nn], derrlocs[nn];
int i;
int errors;
int derrors, kk;
int errval, errloc;
int erasures;
int decoder_errors = 0;
void *rs;
if (e->symsize > 8) {
return -1;
}
/* Compute code parameters */
kk = nn - e->nroots;
rs = init_rs_char(e->symsize, e->genpoly, e->fcs, e->prim, e->nroots, 0);
if (rs == NULL) {
printf("init_rs_char failed!\n");
return -1;
}
/* Test up to the error correction capacity of the code */
for (errors = 0; errors <= e->nroots / 2; errors++) {
/* Load block with random data and encode */
for (i = 0; i < kk; i++) {
block[i] = random() & nn;
}
memcpy(tblock, block, sizeof(block));
encode_rs_char(rs, block, &block[kk]);
/* Make temp copy, seed with errors */
memcpy(tblock, block, sizeof(block));
memset(errlocs, 0, sizeof(errlocs));
memset(derrlocs, 0, sizeof(derrlocs));
erasures = 0;
for (i = 0; i < errors; i++) {
do {
errval = random() & nn;
}
while (errval == 0); /* Error value must be nonzero */
do {
errloc = random() % nn;
}
while (errlocs[errloc] != 0); /* Must not choose the same location twice */
errlocs[errloc] = 1;
#if FLAG_ERASURE
if (random() & 1) { /* 50-50 chance */
derrlocs[erasures++] = errloc;
}
#endif
tblock[errloc] ^= errval;
}
/* Decode the errored block */
derrors = decode_rs_char(rs, tblock, derrlocs, erasures);
if (derrors != errors) {
printf("(%d,%d) decoder says %d errors, true number is %d\n", nn, kk, derrors,
errors);
decoder_errors++;
}
for (i = 0; i < derrors; i++) {
if (errlocs[derrlocs[i]] == 0) {
printf("(%d,%d) decoder indicates error in location %d without error\n", nn, kk,
derrlocs[i]);
decoder_errors++;
}
}
if (memcmp(tblock, block, sizeof(tblock)) != 0) {
printf("(%d,%d) decoder uncorrected errors! output ^ input:", nn, kk);
decoder_errors++;
for (i = 0; i < nn; i++) {
printf(" %02x", tblock[i] ^ block[i]);
}
printf("\n");
}
}
free_rs_char(rs);
return 0;
}
int exercise_int(struct etab *e)
{
int nn = (1 << e->symsize) - 1;
int block[nn], tblock[nn];
int errlocs[nn], derrlocs[nn];
int i;
int errors;
int derrors, kk;
int errval, errloc;
int erasures;
int decoder_errors = 0;
void *rs;
/* Compute code parameters */
kk = nn - e->nroots;
rs = init_rs_int(e->symsize, e->genpoly, e->fcs, e->prim, e->nroots, 0);
if (rs == NULL) {
printf("init_rs_int failed!\n");
return -1;
}
/* Test up to the error correction capacity of the code */
for (errors = 0; errors <= e->nroots / 2; errors++) {
/* Load block with random data and encode */
for (i = 0; i < kk; i++) {
block[i] = random() & nn;
}
memcpy(tblock, block, sizeof(block));
encode_rs_int(rs, block, &block[kk]);
/* Make temp copy, seed with errors */
memcpy(tblock, block, sizeof(block));
memset(errlocs, 0, sizeof(errlocs));
memset(derrlocs, 0, sizeof(derrlocs));
erasures = 0;
for (i = 0; i < errors; i++) {
do {
errval = random() & nn;
}
while (errval == 0); /* Error value must be nonzero */
do {
errloc = random() % nn;
}
while (errlocs[errloc] != 0); /* Must not choose the same location twice */
errlocs[errloc] = 1;
#if FLAG_ERASURE
if (random() & 1) { /* 50-50 chance */
derrlocs[erasures++] = errloc;
}
#endif
tblock[errloc] ^= errval;
}
/* Decode the errored block */
derrors = decode_rs_int(rs, tblock, derrlocs, erasures);
if (derrors != errors) {
printf("(%d,%d) decoder says %d errors, true number is %d\n", nn, kk, derrors,
errors);
decoder_errors++;
}
for (i = 0; i < derrors; i++) {
if (errlocs[derrlocs[i]] == 0) {
printf("(%d,%d) decoder indicates error in location %d without error\n", nn, kk,
derrlocs[i]);
decoder_errors++;
}
}
if (memcmp(tblock, block, sizeof(tblock)) != 0) {
printf("(%d,%d) decoder uncorrected errors! output ^ input:", nn, kk);
decoder_errors++;
for (i = 0; i < nn; i++) {
printf(" %02x", tblock[i] ^ block[i]);
}
printf("\n");
}
}
free_rs_int(rs);
return 0;
}

View File

@ -2,19 +2,19 @@
#include <stdlib.h>
#include "fec.h"
#define MAX_RANDOM 0x7fffffff
#define MAX_RANDOM 0x7fffffff
/* Generate gaussian random double with specified mean and std_dev */
double normal_rand(double mean, double std_dev)
{
double fac,rsq,v1,v2;
double fac, rsq, v1, v2;
static double gset;
static int iset;
if(iset){
if (iset) {
/* Already got one */
iset = 0;
return mean + std_dev*gset;
return mean + std_dev * gset;
}
/* Generate two evenly distributed numbers between -1 and +1
* that are inside the unit circle
@ -22,22 +22,27 @@ double normal_rand(double mean, double std_dev)
do {
v1 = 2.0 * (double)random() / MAX_RANDOM - 1;
v2 = 2.0 * (double)random() / MAX_RANDOM - 1;
rsq = v1*v1 + v2*v2;
} while(rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0*log(rsq)/rsq);
gset = v1*fac;
rsq = v1 * v1 + v2 * v2;
}
while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log(rsq) / rsq);
gset = v1 * fac;
iset++;
return mean + std_dev*v2*fac;
return mean + std_dev * v2 * fac;
}
unsigned char addnoise(int sym,double amp,double gain,double offset,int clip){
unsigned char addnoise(int sym, double amp, double gain, double offset,
int clip)
{
int sample;
sample = offset + gain*normal_rand(sym?amp:-amp,1.0);
sample = offset + gain * normal_rand(sym ? amp : -amp, 1.0);
/* Clip to 8-bit offset range */
if(sample < 0)
if (sample < 0) {
sample = 0;
else if(sample > clip)
}
else if (sample > clip) {
sample = clip;
}
return sample;
}

46
lib/libfec/sqtest.c Normal file
View File

@ -0,0 +1,46 @@
/* Verify correctness of the sum-of-square routines */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
/* These values should trigger leading/trailing array fragment handling */
#define NSAMP 200002
#define OFFSET 1
long long sumsq_wq(signed short *in, int cnt);
long long sumsq_wq_ref(signed short *in, int cnt);
int main()
{
int i;
long long result, rresult;
signed short samples[NSAMP];
srandom(time(NULL));
for (i = 0; i < NSAMP; i++) {
samples[i] = random() & 0xffff;
}
rresult = sumsq_wq(&samples[OFFSET], NSAMP - OFFSET);
result = sumsq_wq(&samples[OFFSET], NSAMP - OFFSET);
if (result == rresult) {
printf("OK\n");
}
else {
printf("sum mismatch: %lld != %lld\n", result, rresult);
}
exit(0);
}
long long sumsq_wq_ref(signed short *in, int cnt)
{
long long sum = 0;
int i;
for (i = 0; i < cnt; i++) {
sum += (long)in[i] * in[i];
}
return sum;
}

51
lib/libfec/sumsq.c Normal file
View File

@ -0,0 +1,51 @@
/* Compute the sum of the squares of a vector of signed shorts
* Copyright 2004 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#include <stdlib.h>
#include "fec.h"
unsigned long long sumsq_port(signed short *, int);
#ifdef __i386__
unsigned long long sumsq_mmx(signed short *, int);
unsigned long long sumsq_sse(signed short *, int);
unsigned long long sumsq_sse2(signed short *, int);
#endif
#ifdef __x86_64__
unsigned long long sumsq_sse2(signed short *, int);
#endif
#ifdef __VEC__
unsigned long long sumsq_av(signed short *, int);
#endif
unsigned long long sumsq(signed short *in, int cnt)
{
switch (Cpu_mode) {
case PORT:
default:
return sumsq_port(in, cnt);
#ifdef __i386__
case SSE:
case MMX:
return sumsq_mmx(in, cnt);
case SSE2:
return sumsq_sse2(in, cnt);
#endif
#ifdef __x86_64__
case SSE2:
return sumsq_port(in, cnt);
//return sumsq_sse2(in,cnt);
#endif
#ifdef __VEC__
case ALTIVEC:
return sumsq_av(in, cnt);
#endif
}
}

84
lib/libfec/sumsq_av.c Normal file
View File

@ -0,0 +1,84 @@
/* Compute the sum of the squares of a vector of signed shorts
* This is the Altivec SIMD version. It's a little hairy because Altivec
* does not do 64-bit operations directly, so we have to accumulate separate
* 32-bit sums and carries
* Copyright 2004 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
#include "fec.h"
unsigned long long sumsq_av(signed short *in, int cnt)
{
long long sum;
vector signed short x;
vector unsigned int sums, carries, s1, s2;
int pad;
union {
vector unsigned char cv;
vector unsigned int iv;
unsigned int w[4];
unsigned char c[16];
} s;
carries = sums = (vector unsigned int)(0);
if ((pad = (int)in & 15) != 0) {
/* Load unaligned leading word */
x = vec_perm(vec_ld(0, in), (vector signed short)(0), vec_lvsl(0, in));
if (cnt < 8) { /* Shift right to chop stuff beyond end of short block */
s.c[15] = (8 - cnt) << 4;
x = vec_sro(x, s.cv);
}
sums = (vector unsigned int)vec_msum(x, x, (vector signed int)(0));
in += 8 - pad / 2;
cnt -= 8 - pad / 2;
}
/* Everything is now aligned, rip through most of the block */
while (cnt >= 8) {
x = vec_ld(0, in);
/* A single vec_msum cannot overflow, but we have to sum it with
* the earlier terms separately to handle the carries
* The cast to unsigned is OK because squares are always positive
*/
s1 = (vector unsigned int)vec_msum(x, x, (vector signed int)(0));
carries = vec_add(carries, vec_addc(sums, s1));
sums = vec_add(sums, s1);
in += 8;
cnt -= 8;
}
/* Handle trailing fragment, if any */
if (cnt > 0) {
x = vec_ld(0, in);
s.c[15] = (8 - cnt) << 4;
x = vec_sro(x, s.cv);
s1 = (vector unsigned int)vec_msum(x, x, (vector signed int)(0));
carries = vec_add(carries, vec_addc(sums, s1));
sums = vec_add(sums, s1);
}
/* Combine 4 sub-sums and carries */
s.c[15] = 64; /* Shift right two 32-bit words */
s1 = vec_sro(sums, s.cv);
s2 = vec_sro(carries, s.cv);
carries = vec_add(carries, vec_addc(sums, s1));
sums = vec_add(sums, s1);
carries = vec_add(carries, s2);
s.c[15] = 32; /* Shift right one 32-bit word */
s1 = vec_sro(sums, s.cv);
s2 = vec_sro(carries, s.cv);
carries = vec_add(carries, vec_addc(sums, s1));
sums = vec_add(sums, s1);
carries = vec_add(carries, s2);
/* Extract sum and carries from right-hand words and combine into result */
s.iv = sums;
sum = s.w[3];
s.iv = carries;
sum += (long long)s.w[3] << 32;
return sum;
}

View File

@ -10,23 +10,24 @@
* May be used under the terms of the GNU Lesser Public License (LGPL)
*/
long long sumsq_mmx_assist(signed short *,int);
long long sumsq_mmx_assist(signed short *, int);
long long sumsq_mmx(signed short *in,int cnt){
long long sumsq_mmx(signed short *in, int cnt)
{
long long sum = 0;
/* Handle stuff before the next 8-byte boundary */
while(((int)in & 7) != 0 && cnt != 0){
while (((int)in & 7) != 0 && cnt != 0) {
sum += (long)in[0] * in[0];
in++;
cnt--;
}
sum += sumsq_mmx_assist(in,cnt);
sum += sumsq_mmx_assist(in, cnt);
in += cnt & ~7;
cnt &= 7;
/* Handle up to 7 words at end */
while(cnt != 0){
while (cnt != 0) {
sum += (long)in[0] * in[0];
in++;
cnt--;

View File

@ -5,11 +5,12 @@
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*/
unsigned long long sumsq_port(signed short *in,int cnt){
unsigned long long sumsq_port(signed short *in, int cnt)
{
long long sum = 0;
int i;
for(i=0;i<cnt;i++){
for (i = 0; i < cnt; i++) {
sum += (int)in[i] * (int)in[i];
}
return sum;

Some files were not shown because too many files have changed in this diff Show More