diff options
| author | Johnathan Corgan <johnathan@corganlabs.com> | 2012-11-12 20:36:05 -0800 |
|---|---|---|
| committer | Johnathan Corgan <johnathan@corganlabs.com> | 2012-11-12 20:36:05 -0800 |
| commit | ce1213e615023fa5c43f75ba2d2ec54649025b8f (patch) | |
| tree | 8567615201922dc4f1196a965161f4bba0e24dd7 /gr-fec/lib | |
| parent | 5b568abe194f027860670317915c27c174fab52e (diff) | |
fec: add new top-level component for FEC blocks
Removed viterbi and reed-solomon from gnuradio-core
Diffstat (limited to 'gr-fec/lib')
31 files changed, 2481 insertions, 0 deletions
diff --git a/gr-fec/lib/CMakeLists.txt b/gr-fec/lib/CMakeLists.txt new file mode 100644 index 0000000000..422acb3607 --- /dev/null +++ b/gr-fec/lib/CMakeLists.txt @@ -0,0 +1,56 @@ +# Copyright 2012 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# 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. + +# These are convenience libraries of 3rd party code. +# Change to test for distro provided packages +GR_INCLUDE_SUBDIRECTORY(reed-solomon) +GR_INCLUDE_SUBDIRECTORY(viterbi) + +######################################################################## +# Setup the include and linker paths +######################################################################## +include_directories( + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_BINARY_DIR} + ${GNURADIO_CORE_INCLUDE_DIRS} + ${GR_FEC_INCLUDE_DIRS} + ${Boost_INCLUDE_DIRS} +) + +link_directories(${FEC_LIBRARY_DIRS}) +link_directories(${Boost_LIBRARY_DIRS}) + +######################################################################## +# Setup library +######################################################################## +list(APPEND gnuradio_fec_sources + decode_ccsds_27_fb_impl.cc + encode_ccsds_27_bb_impl.cc +) + +list(APPEND gnuradio_fec_libs + gnuradio-core + ${Boost_LIBRARIES} + ${FEC_LIBRARIES} +) + +add_library(gnuradio-fec SHARED ${gnuradio_fec_sources}) +target_link_libraries(gnuradio-fec ${gnuradio_fec_libs}) +GR_LIBRARY_FOO(gnuradio-fec RUNTIME_COMPONENT "fec_runtime" DEVEL_COMPONENT "fec_devel") + diff --git a/gr-fec/lib/decode_ccsds_27_fb_impl.cc b/gr-fec/lib/decode_ccsds_27_fb_impl.cc new file mode 100644 index 0000000000..87bdc46402 --- /dev/null +++ b/gr-fec/lib/decode_ccsds_27_fb_impl.cc @@ -0,0 +1,89 @@ +/* -*- c++ -*- */ +/* + * Copyright 2012 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "decode_ccsds_27_fb_impl.h" +#include <gr_io_signature.h> + +namespace gr { + namespace fec { + + decode_ccsds_27_fb::sptr decode_ccsds_27_fb::make() + { + return gnuradio::get_initial_sptr(new decode_ccsds_27_fb_impl()); + } + + decode_ccsds_27_fb_impl::decode_ccsds_27_fb_impl() + : gr_sync_decimator("decode_ccsds_27_fb", + gr_make_io_signature (1, 1, sizeof(float)), + gr_make_io_signature (1, 1, sizeof(char)), + 2*8) // Rate 1/2 code, unpacked to packed conversion + { + float RATE = 0.5; + float ebn0 = 12.0; + float esn0 = RATE*pow(10.0, ebn0/10.0); + + gen_met(d_mettab, 100, esn0, 0.0, 256); + viterbi_chunks_init(d_state0); + } + + int + decode_ccsds_27_fb_impl::work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) + { + const float *in = (const float *)input_items[0]; + unsigned char *out = (unsigned char *)output_items[0]; + + for (int i = 0; i < noutput_items*16; i++) { + // Translate and clip [-1.0..1.0] to [28..228] + float sample = in[i]*100.0+128.0; + if (sample > 255.0) + sample = 255.0; + else if (sample < 0.0) + sample = 0.0; + unsigned char sym = (unsigned char)(floor(sample)); + + d_viterbi_in[d_count % 4] = sym; + if ((d_count % 4) == 3) { + // Every fourth symbol, perform butterfly operation + viterbi_butterfly2(d_viterbi_in, d_mettab, d_state0, d_state1); + + // Every sixteenth symbol, read out a byte + if (d_count % 16 == 11) { + // long metric = + viterbi_get_output(d_state0, out++); + // printf("%li\n", *(out-1), metric); + } + } + + d_count++; + } + + return noutput_items; + } + + } /* namespace fec */ +}/* namespace gr */ diff --git a/gr-fec/lib/decode_ccsds_27_fb_impl.h b/gr-fec/lib/decode_ccsds_27_fb_impl.h new file mode 100644 index 0000000000..c02e284e70 --- /dev/null +++ b/gr-fec/lib/decode_ccsds_27_fb_impl.h @@ -0,0 +1,57 @@ +/* -*- c++ -*- */ +/* + * Copyright 2012 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +#ifndef INCLUDED_FEC_DECODE_CCSDS_27_FB_IMPL_H +#define INCLUDED_FEC_DECODE_CCSDS_27_FB_IMPL_H + +#include <fec/decode_ccsds_27_fb.h> + +extern "C" { +#include "viterbi/viterbi.h" +} + +namespace gr { + namespace fec { + + class FEC_API decode_ccsds_27_fb_impl : public decode_ccsds_27_fb + { + private: + // Viterbi state + int d_mettab[2][256]; + struct viterbi_state d_state0[64]; + struct viterbi_state d_state1[64]; + unsigned char d_viterbi_in[16]; + + int d_count; + + public: + decode_ccsds_27_fb_impl(); + + int work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + }; + + } /* namespace fec */ +} /* namespace gr */ + +#endif /* INCLUDED_FEC_DECODE_CCSDS_27_FB_IMPL_H */ diff --git a/gr-fec/lib/encode_ccsds_27_bb_impl.cc b/gr-fec/lib/encode_ccsds_27_bb_impl.cc new file mode 100644 index 0000000000..28ee450ef4 --- /dev/null +++ b/gr-fec/lib/encode_ccsds_27_bb_impl.cc @@ -0,0 +1,65 @@ +/* -*- c++ -*- */ +/* + * Copyright 2012 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "encode_ccsds_27_bb_impl.h" +#include <gr_io_signature.h> + +extern "C" { +#include "viterbi/viterbi.h" +} + +namespace gr { + namespace fec { + + encode_ccsds_27_bb::sptr encode_ccsds_27_bb::make() + { + return gnuradio::get_initial_sptr(new encode_ccsds_27_bb_impl()); + } + + encode_ccsds_27_bb_impl::encode_ccsds_27_bb_impl() + : gr_sync_interpolator("encode_ccsds_27_bb", + gr_make_io_signature (1, 1, sizeof(char)), + gr_make_io_signature (1, 1, sizeof(char)), + 16) // Rate 1/2 code, packed to unpacked conversion + { + d_encstate = 0; + } + + int + encode_ccsds_27_bb_impl::work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) + { + unsigned char *in = (unsigned char *)input_items[0]; + unsigned char *out = (unsigned char *)output_items[0]; + + d_encstate = encode(out, in, noutput_items/16, d_encstate); + + return noutput_items; + } + + } /* namespace fec */ +}/* namespace gr */ diff --git a/gr-fec/lib/encode_ccsds_27_bb_impl.h b/gr-fec/lib/encode_ccsds_27_bb_impl.h new file mode 100644 index 0000000000..140f8c2ae4 --- /dev/null +++ b/gr-fec/lib/encode_ccsds_27_bb_impl.h @@ -0,0 +1,47 @@ +/* -*- c++ -*- */ +/* + * Copyright 2012 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +#ifndef INCLUDED_FEC_ENCODE_CCSDS_27_BB_IMPL_H +#define INCLUDED_FEC_ENCODE_CCSDS_27_BB_IMPL_H + +#include <fec/encode_ccsds_27_bb.h> + +namespace gr { + namespace fec { + + class FEC_API encode_ccsds_27_bb_impl : public encode_ccsds_27_bb + { + private: + unsigned char d_encstate; + + public: + encode_ccsds_27_bb_impl(); + + int work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + }; + + } /* namespace fec */ +} /* namespace gr */ + +#endif /* INCLUDED_FEC_ENCODE_CCSDS_27_BB_IMPL_H */ diff --git a/gr-fec/lib/reed-solomon/CMakeLists.txt b/gr-fec/lib/reed-solomon/CMakeLists.txt new file mode 100644 index 0000000000..0858a670c5 --- /dev/null +++ b/gr-fec/lib/reed-solomon/CMakeLists.txt @@ -0,0 +1,62 @@ +# Copyright 2010-2012 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# 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. + +######################################################################## +# This file included, use CMake directory variables +######################################################################## +#MSVC workaround: we cant have dynamically sized arrays. +#So ifdef a max array bounds that is larger than NN and NROOTS +#Its a bit of a hack, but if you look at the code, its so full of ifdefs, +#and lacks optimization where it should be pre-allocating these arrays. +if(MSVC) + set_source_files_properties( + ${CMAKE_CURRENT_SOURCE_DIR}/exercise.c + ${CMAKE_CURRENT_SOURCE_DIR}/decode_rs.c + PROPERTIES COMPILE_DEFINITIONS "MAX_ARRAY=256;" + ) +endif(MSVC) + +set(gr_fec_rs_sources + ${CMAKE_CURRENT_SOURCE_DIR}/encode_rs.c + ${CMAKE_CURRENT_SOURCE_DIR}/decode_rs.c + ${CMAKE_CURRENT_SOURCE_DIR}/init_rs.c +) + +######################################################################## +# Setup sources and includes +######################################################################## +list(APPEND gnuradio_fec_sources ${gr_fec_rs_sources}) + +#install( +# FILES ${CMAKE_CURRENT_SOURCE_DIR}/rs.h +# DESTINATION ${GR_INCLUDE_DIR}/gnuradio/fec +# COMPONENT "fec_devel" +#) + +######################################################################## +# Register unit tests +######################################################################## +if(ENABLE_TESTING) +add_executable(gr_fec_rstest + ${gr_fec_rs_sources} + ${CMAKE_CURRENT_SOURCE_DIR}/rstest.c + ${CMAKE_CURRENT_SOURCE_DIR}/exercise.c +) +add_test(gr-fec-reed-solomon-test gr_fec_rstest) +endif(ENABLE_TESTING) diff --git a/gr-fec/lib/reed-solomon/Makefile.in.karn b/gr-fec/lib/reed-solomon/Makefile.in.karn new file mode 100644 index 0000000000..8550b41581 --- /dev/null +++ b/gr-fec/lib/reed-solomon/Makefile.in.karn @@ -0,0 +1,99 @@ +# Copyright 2002 Phil Karn, KA9Q +# May be used under the terms of the GNU General Public License (GPL) +# @configure_input@ +srcdir = @srcdir@ +prefix = @prefix@ +exec_prefix=@exec_prefix@ +VPATH = @srcdir@ +CC=@CC@ + +CFLAGS=@CFLAGS@ @ARCH_OPTION@ -Wall + +LIB= encode_rs_char.o encode_rs_int.o encode_rs_8.o \ + decode_rs_char.o decode_rs_int.o decode_rs_8.o \ + init_rs_char.o init_rs_int.o ccsds_tab.o \ + encode_rs_ccsds.o decode_rs_ccsds.o ccsds_tal.o + +all: librs.a librs.so.@SO_VERSION@ + +test: rstest + ./rstest + +rstest: rstest.o exercise_int.o exercise_char.o exercise_8.o exercise_ccsds.o \ + librs.a + gcc -g -o $@ $^ + +install: all + install -D -m 644 -p librs.a librs.so.@SO_VERSION@ @libdir@ + (cd @libdir@;ln -f -s librs.so.@SO_VERSION@ librs.so) + ldconfig + install -m 644 -p rs.h @includedir@ + install -m 644 rs.3 @mandir@/man3 + +librs.a: $(LIB) + ar rv $@ $^ + +librs.so.@SO_VERSION@: librs.a + gcc -shared -Xlinker -soname=librs.so.@SO_NAME@ -o $@ -Wl,-whole-archive $^ -Wl,-no-whole-archive -lc + +encode_rs_char.o: encode_rs.c + gcc $(CFLAGS) -c -o $@ $^ + +encode_rs_int.o: encode_rs.c + gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^ + +encode_rs_8.o: encode_rs.c + gcc -DFIXED=1 $(CFLAGS) -c -o $@ $^ + +decode_rs_char.o: decode_rs.c + gcc $(CFLAGS) -c -o $@ $^ + +decode_rs_int.o: decode_rs.c + gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^ + +decode_rs_8.o: decode_rs.c + gcc -DFIXED=1 $(CFLAGS) -c -o $@ $^ + +init_rs_char.o: init_rs.c + gcc $(CFLAGS) -c -o $@ $^ + +init_rs_int.o: init_rs.c + gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^ + +ccsds_tab.o: ccsds_tab.c + +ccsds_tab.c: gen_ccsds + ./gen_ccsds > ccsds_tab.c + +gen_ccsds: gen_ccsds.o init_rs_char.o + gcc -o $@ $^ + +gen_ccsds.o: gen_ccsds.c + gcc $(CFLAGS) -c -o $@ $^ + +ccsds_tal.o: ccsds_tal.c + +ccsds_tal.c: gen_ccsds_tal + ./gen_ccsds_tal > ccsds_tal.c + +exercise_char.o: exercise.c + gcc $(CFLAGS) -c -o $@ $^ + +exercise_int.o: exercise.c + gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^ + +exercise_8.o: exercise.c + gcc -DFIXED=1 $(CFLAGS) -c -o $@ $^ + +exercise_ccsds.o: exercise.c + gcc -DCCSDS=1 $(CFLAGS) -c -o $@ $^ + + +clean: + rm -f *.o *.a ccsds_tab.c ccsds_tal.c gen_ccsds gen_ccsds_tal \ + rstest librs.so.@SO_VERSION@ + +distclean: clean + rm -f config.log config.cache config.status config.h makefile + + diff --git a/gr-fec/lib/reed-solomon/README b/gr-fec/lib/reed-solomon/README new file mode 100644 index 0000000000..5c867638ef --- /dev/null +++ b/gr-fec/lib/reed-solomon/README @@ -0,0 +1,2 @@ +This code is from http://people.qualcomm.com/karn/code/fec +It is based on reed-soloman-3.1.1 (1 Jan 2002). diff --git a/gr-fec/lib/reed-solomon/README.karn b/gr-fec/lib/reed-solomon/README.karn new file mode 100644 index 0000000000..f30644ffea --- /dev/null +++ b/gr-fec/lib/reed-solomon/README.karn @@ -0,0 +1,22 @@ +This package implements a general purpose Reed-Solomon encoding and decoding +facility. See the rs.3 man page for details. + +To install, simply do the following after extracting this tarball into +an empty directory: + +./configure +make +make install + +The command "make test" runs a battery of encode/decode tests using a +variety of RS codes using random data and random errors. Each test +should pass with an "OK"; if any fail, please let me know so I can +track down the problem. + +Phil Karn (karn@ka9q.net) 1 Jan 2002 + +Copyright 2002, Phil Karn, KA9Q +This software may be used under the terms of the GNU General Public License (GPL). + + + diff --git a/gr-fec/lib/reed-solomon/ccsds.h b/gr-fec/lib/reed-solomon/ccsds.h new file mode 100644 index 0000000000..0f2bde6186 --- /dev/null +++ b/gr-fec/lib/reed-solomon/ccsds.h @@ -0,0 +1 @@ +extern unsigned char Taltab[],Tal1tab[]; diff --git a/gr-fec/lib/reed-solomon/char.h b/gr-fec/lib/reed-solomon/char.h new file mode 100644 index 0000000000..bc4fd8a1ba --- /dev/null +++ b/gr-fec/lib/reed-solomon/char.h @@ -0,0 +1,57 @@ +/* Include file to configure the RS codec for character symbols + * + * Copyright 2002, Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ + +#define DTYPE unsigned char + +#include <fec/api.h> + +/* Reed-Solomon codec control block */ +struct rs { + unsigned int mm; /* Bits per symbol */ + unsigned int nn; /* Symbols per block (= (1<<mm)-1) */ + unsigned char *alpha_to; /* log lookup table */ + unsigned char *index_of; /* Antilog lookup table */ + unsigned char *genpoly; /* Generator polynomial */ + unsigned int nroots; /* Number of generator roots = number of parity symbols */ + unsigned char fcr; /* First consecutive root, index form */ + unsigned char prim; /* Primitive element, index form */ + unsigned char iprim; /* prim-th root of 1, index form */ +}; + +static inline unsigned int modnn(struct rs *rs, unsigned int x){ + while (x >= rs->nn) { + x -= rs->nn; + x = (x >> rs->mm) + (x & rs->nn); + } + return x; +} +#define MODNN(x) modnn(rs,x) + +#define MM (rs->mm) +#define NN (rs->nn) +#define ALPHA_TO (rs->alpha_to) +#define INDEX_OF (rs->index_of) +#define GENPOLY (rs->genpoly) +#define NROOTS (rs->nroots) +#define FCR (rs->fcr) +#define PRIM (rs->prim) +#define IPRIM (rs->iprim) +#define A0 (NN) + +#define ENCODE_RS encode_rs_char +#define DECODE_RS decode_rs_char +#define INIT_RS init_rs_char +#define FREE_RS free_rs_char + +FEC_API void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity); +FEC_API int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras); +FEC_API void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned int fcr, + unsigned int prim,unsigned int nroots); +FEC_API void FREE_RS(void *p); + + + + diff --git a/gr-fec/lib/reed-solomon/decode_rs.c b/gr-fec/lib/reed-solomon/decode_rs.c new file mode 100644 index 0000000000..f9438cf266 --- /dev/null +++ b/gr-fec/lib/reed-solomon/decode_rs.c @@ -0,0 +1,270 @@ +/* Reed-Solomon decoder + * Copyright 2002 Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ + +#ifdef DEBUG +#include <stdio.h> +#endif + +#include <string.h> + +#define NULL ((void *)0) +#define min(a,b) ((a) < (b) ? (a) : (b)) + +#ifdef FIXED +#include "fixed.h" +#elif defined(BIGSYM) +#include "int.h" +#else +#include "char.h" +#endif + +int DECODE_RS( +#ifndef FIXED +void *p, +#endif +DTYPE *data, int *eras_pos, int no_eras){ + +#ifndef FIXED + struct rs *rs = (struct rs *)p; +#endif + int deg_lambda, el, deg_omega; + int i, j, r, k; +#ifdef MAX_ARRAY + DTYPE u,q,tmp,num1,num2,den,discr_r; + DTYPE lambda[MAX_ARRAY], s[MAX_ARRAY]; /* Err+Eras Locator poly + * and syndrome poly */ + DTYPE b[MAX_ARRAY], t[MAX_ARRAY], omega[MAX_ARRAY]; + DTYPE root[MAX_ARRAY], reg[MAX_ARRAY], loc[MAX_ARRAY]; +#else + DTYPE u,q,tmp,num1,num2,den,discr_r; + DTYPE lambda[NROOTS+1], s[NROOTS]; /* Err+Eras Locator poly + * and syndrome poly */ + DTYPE b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1]; + DTYPE root[NROOTS], reg[NROOTS+1], loc[NROOTS]; +#endif + int syn_error, count; + + /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */ + for(i=0;(unsigned int)i<NROOTS;i++) + s[i] = data[0]; + + for(j=1;(unsigned int)j<NN;j++){ + for(i=0;(unsigned int)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;(unsigned int)i<NROOTS;i++){ + syn_error |= s[i]; + s[i] = INDEX_OF[s[i]]; + } + + if (!syn_error) { + /* if syndrome is zero, data[] is a codeword and there are no + * errors to correct. So return data[] unmodified + */ + count = 0; + goto finish; + } + memset(&lambda[1],0,NROOTS*sizeof(lambda[0])); + lambda[0] = 1; + + if (no_eras > 0) { + /* Init lambda to be the erasure locator polynomial */ + lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))]; + for (i = 1; i < no_eras; i++) { + u = MODNN(PRIM*(NN-1-eras_pos[i])); + for (j = i+1; j > 0; j--) { + tmp = INDEX_OF[lambda[j - 1]]; + if(tmp != A0) + lambda[j] ^= ALPHA_TO[MODNN(u + tmp)]; + } + } + +#if DEBUG >= 1 + /* Test code that verifies the erasure locator polynomial just constructed + Needed only for decoder debugging. */ + + /* find roots of the erasure location polynomial */ + for(i=1;i<=no_eras;i++) + reg[i] = INDEX_OF[lambda[i]]; + + count = 0; + for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) { + q = 1; + for (j = 1; j <= no_eras; j++) + if (reg[j] != A0) { + reg[j] = MODNN(reg[j] + j); + q ^= ALPHA_TO[reg[j]]; + } + if (q != 0) + continue; + /* store root and error location number indices */ + root[count] = i; + loc[count] = k; + count++; + } + if (count != no_eras) { + printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras); + count = -1; + goto finish; + } +#if DEBUG >= 2 + printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n"); + for (i = 0; i < count; i++) + printf("%d ", loc[i]); + printf("\n"); +#endif +#endif + } + for(i=0;(unsigned int)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 ((unsigned int)(++r) <= NROOTS) { /* r is the step number */ + /* Compute discrepancy at the r-th step in poly-form */ + discr_r = 0; + for (i = 0; i < r; i++){ + if ((lambda[i] != 0) && (s[r-i-1] != A0)) { + discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])]; + } + } + discr_r = INDEX_OF[discr_r]; /* Index form */ + if (discr_r == A0) { + /* 2 lines below: B(x) <-- x*B(x) */ + memmove(&b[1],b,NROOTS*sizeof(b[0])); + b[0] = A0; + } else { + /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */ + t[0] = lambda[0]; + for (i = 0 ; (unsigned int)i < NROOTS; i++) { + if(b[i] != A0) + t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])]; + else + t[i+1] = lambda[i+1]; + } + if (2 * el <= r + no_eras - 1) { + el = r + no_eras - el; + /* + * 2 lines below: B(x) <-- inv(discr_r) * + * lambda(x) + */ + for (i = 0; (unsigned int)i <= NROOTS; i++) + b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN); + } else { + /* 2 lines below: B(x) <-- x*B(x) */ + memmove(&b[1],b,NROOTS*sizeof(b[0])); + b[0] = A0; + } + memcpy(lambda,t,(NROOTS+1)*sizeof(t[0])); + } + } + + /* Convert lambda to index form and compute deg(lambda(x)) */ + deg_lambda = 0; + for(i=0;(unsigned int)i<NROOTS+1;i++){ + lambda[i] = INDEX_OF[lambda[i]]; + if(lambda[i] != A0) + deg_lambda = i; + } + /* Find roots of the error+erasure locator polynomial by Chien search */ + memcpy(®[1],&lambda[1],NROOTS*sizeof(reg[0])); + count = 0; /* Number of roots of lambda(x) */ + for (i = 1,k=IPRIM-1; (unsigned int)i <= NN; i++,k = MODNN(k+IPRIM)) { + q = 1; /* lambda[0] is always 0 */ + for (j = deg_lambda; j > 0; j--){ + if (reg[j] != A0) { + reg[j] = MODNN(reg[j] + j); + q ^= ALPHA_TO[reg[j]]; + } + } + if (q != 0) + continue; /* Not a root */ + /* store root (index-form) and error location number */ +#if DEBUG>=2 + printf("count %d root %d loc %d\n",count,i,k); +#endif + root[count] = i; + loc[count] = k; + /* If we've already found max possible roots, + * abort the search to save time + */ + if(++count == deg_lambda) + break; + } + if (deg_lambda != count) { + /* + * deg(lambda) unequal to number of roots => uncorrectable + * error detected + */ + count = -1; + goto finish; + } + /* + * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo + * x**NROOTS). in index form. Also find deg(omega). + */ + deg_omega = 0; + for (i = 0; (unsigned int)i < NROOTS;i++){ + tmp = 0; + j = (deg_lambda < i) ? deg_lambda : i; + for(;j >= 0; j--){ + if ((s[i - j] != A0) && (lambda[j] != A0)) + tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])]; + } + if(tmp != 0) + deg_omega = i; + omega[i] = INDEX_OF[tmp]; + } + omega[NROOTS] = A0; + + /* + * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 = + * inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form + */ + for (j = count-1; j >=0; j--) { + num1 = 0; + for (i = deg_omega; i >= 0; i--) { + if (omega[i] != A0) + num1 ^= ALPHA_TO[MODNN(omega[i] + i * root[j])]; + } + num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)]; + den = 0; + + /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */ + for (i = (int)min((unsigned int)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 (den == 0) { +#if DEBUG >= 1 + printf("\n ERROR: denominator = 0\n"); +#endif + count = -1; + goto finish; + } + /* Apply error to data */ + if (num1 != 0) { + data[loc[j]] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])]; + } + } + finish: + if(eras_pos != NULL){ + for(i=0;i<count;i++) + eras_pos[i] = loc[i]; + } + return count; +} diff --git a/gr-fec/lib/reed-solomon/decode_rs_ccsds.c b/gr-fec/lib/reed-solomon/decode_rs_ccsds.c new file mode 100644 index 0000000000..2543d3a640 --- /dev/null +++ b/gr-fec/lib/reed-solomon/decode_rs_ccsds.c @@ -0,0 +1,27 @@ +/* This function wraps around the fixed 8-bit decoder, performing the + * basis transformations necessary to meet the CCSDS standard + * + * Copyright 2002, Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#define FIXED 1 +#include "fixed.h" +#include "ccsds.h" + +int decode_rs_ccsds(unsigned char *data,int *eras_pos,int no_eras){ + int i,r; + unsigned char cdata[NN]; + + /* Convert data from dual basis to conventional */ + for(i=0;i<NN;i++) + cdata[i] = Tal1tab[data[i]]; + + r = decode_rs_8(cdata,eras_pos,no_eras); + + if(r > 0){ + /* Convert from conventional to dual basis */ + for(i=0;i<NN;i++) + data[i] = Taltab[cdata[i]]; + } + return r; +} diff --git a/gr-fec/lib/reed-solomon/encode_rs.c b/gr-fec/lib/reed-solomon/encode_rs.c new file mode 100644 index 0000000000..cd31f32c6e --- /dev/null +++ b/gr-fec/lib/reed-solomon/encode_rs.c @@ -0,0 +1,47 @@ +/* Reed-Solomon encoder + * Copyright 2002, Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#include <string.h> + +#ifdef FIXED +#include "fixed.h" +#elif defined(BIGSYM) +#include "int.h" +#else +#include "char.h" +#endif + +void ENCODE_RS( +#ifndef FIXED +void *p, +#endif +DTYPE *data, DTYPE *bb){ +#ifndef FIXED + struct rs *rs = (struct rs *)p; +#endif + unsigned int i, j; + DTYPE feedback; + + memset(bb,0,NROOTS*sizeof(DTYPE)); + + for(i=0;i<NN-NROOTS;i++){ + feedback = INDEX_OF[data[i] ^ bb[0]]; + if(feedback != A0){ /* feedback term is non-zero */ +#ifdef UNNORMALIZED + /* This line is unnecessary when GENPOLY[NROOTS] is unity, as it must + * always be for the polynomials constructed by init_rs() + */ + feedback = MODNN(NN - GENPOLY[NROOTS] + feedback); +#endif + for(j=1;j<NROOTS;j++) + bb[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS-j])]; + } + /* Shift */ + memmove(&bb[0],&bb[1],sizeof(DTYPE)*(NROOTS-1)); + if(feedback != A0) + bb[NROOTS-1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])]; + else + bb[NROOTS-1] = 0; + } +} diff --git a/gr-fec/lib/reed-solomon/encode_rs_ccsds.c b/gr-fec/lib/reed-solomon/encode_rs_ccsds.c new file mode 100644 index 0000000000..a748b34689 --- /dev/null +++ b/gr-fec/lib/reed-solomon/encode_rs_ccsds.c @@ -0,0 +1,24 @@ +/* This function wraps around the fixed 8-bit encoder, performing the + * basis transformations necessary to meet the CCSDS standard + * + * Copyright 2002, Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#define FIXED +#include "fixed.h" +#include "ccsds.h" + +void encode_rs_ccsds(unsigned char *data,unsigned char *parity){ + int i; + unsigned char cdata[NN-NROOTS]; + + /* Convert data from dual basis to conventional */ + for(i=0;i<NN-NROOTS;i++) + cdata[i] = Tal1tab[data[i]]; + + encode_rs_8(cdata,parity); + + /* Convert parity from conventional to dual basis */ + for(i=0;i<NN-NROOTS;i++) + parity[i] = Taltab[parity[i]]; +} diff --git a/gr-fec/lib/reed-solomon/exercise.c b/gr-fec/lib/reed-solomon/exercise.c new file mode 100644 index 0000000000..de33a6bff3 --- /dev/null +++ b/gr-fec/lib/reed-solomon/exercise.c @@ -0,0 +1,134 @@ +/* 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 General Public License (GPL) + */ +#define FLAG_ERASURE 1 /* Randomly flag 50% of errors as erasures */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#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 +#if MAX_ARRAY + DTYPE block[MAX_ARRAY],tblock[MAX_ARRAY]; + unsigned int i; + int errors; + int errlocs[MAX_ARRAY]; + int derrlocs[MAX_ARRAY]; +#else + DTYPE block[NN],tblock[NN]; + unsigned int i; + int errors; + int errlocs[NN]; + int derrlocs[NROOTS]; +#endif + 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;(unsigned int)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]); +#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<(unsigned int)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); +#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<(unsigned int)derrors;i++){ + if(errlocs[derrlocs[i]] == 0){ + PRINTPARM + printf(" decoder indicates error in location %d without error\n",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; +} diff --git a/gr-fec/lib/reed-solomon/fixed.h b/gr-fec/lib/reed-solomon/fixed.h new file mode 100644 index 0000000000..ee4ae2e714 --- /dev/null +++ b/gr-fec/lib/reed-solomon/fixed.h @@ -0,0 +1,40 @@ +/* Configure the RS codec with fixed parameters for CCSDS standard + * (255,223) code over GF(256). Note: the conventional basis is still + * used; the dual-basis mappings are performed in [en|de]code_rs_ccsds.c + * + * Copyright 2002 Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#define DTYPE unsigned char + +#include <fec/api.h> + +static inline int mod255(int x){ + while (x >= 255) { + x -= 255; + x = (x >> 8) + (x & 255); + } + return x; +} +#define MODNN(x) mod255(x) + +extern unsigned char CCSDS_alpha_to[]; +extern unsigned char CCSDS_index_of[]; +extern unsigned char CCSDS_poly[]; + +#define MM 8 +#define NN 255 +#define ALPHA_TO CCSDS_alpha_to +#define INDEX_OF CCSDS_index_of +#define GENPOLY CCSDS_poly +#define NROOTS 32 +#define FCR 112 +#define PRIM 11 +#define IPRIM 116 +#define A0 (NN) + +#define ENCODE_RS encode_rs_8 +#define DECODE_RS decode_rs_8 + +FEC_API void ENCODE_RS(DTYPE *data,DTYPE *parity); +FEC_API int DECODE_RS(DTYPE *data, int *eras_pos, int no_eras);
\ No newline at end of file diff --git a/gr-fec/lib/reed-solomon/gen_ccsds.c b/gr-fec/lib/reed-solomon/gen_ccsds.c new file mode 100644 index 0000000000..1e4e4f5363 --- /dev/null +++ b/gr-fec/lib/reed-solomon/gen_ccsds.c @@ -0,0 +1,34 @@ +/* Generate tables for CCSDS code + * Copyright 2002 Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#include <stdio.h> +#include "char.h" + +int main(){ + struct rs *rs; + int i; + + rs = init_rs_char(8,0x187,112,11,32); /* CCSDS standard */ + printf("unsigned char CCSDS_alpha_to[] = {"); + for(i=0;i<256;i++){ + if((i % 16) == 0) + printf("\n"); + printf("0x%02x,",rs->alpha_to[i]); + } + printf("\n};\n\nunsigned char CCSDS_index_of[] = {"); + for(i=0;i<256;i++){ + if((i % 16) == 0) + printf("\n"); + printf("%3d,",rs->index_of[i]); + } + printf("\n};\n\nunsigned char CCSDS_poly[] = {"); + for(i=0;i<33;i++){ + if((i % 16) == 0) + printf("\n"); + + printf("%3d,",rs->genpoly[i]); + } + printf("\n};\n"); + exit(0); +} diff --git a/gr-fec/lib/reed-solomon/gen_ccsds_tal.c b/gr-fec/lib/reed-solomon/gen_ccsds_tal.c new file mode 100644 index 0000000000..9dde18917b --- /dev/null +++ b/gr-fec/lib/reed-solomon/gen_ccsds_tal.c @@ -0,0 +1,50 @@ +/* Conversion lookup tables from conventional alpha to Berlekamp's + * dual-basis representation. Used in the CCSDS version only. + * taltab[] -- convert conventional to dual basis + * tal1tab[] -- convert dual basis to conventional + + * Note: the actual RS encoder/decoder works with the conventional basis. + * So data is converted from dual to conventional basis before either + * encoding or decoding and then converted back. + * + * Copyright 2002 Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#include <stdio.h> +unsigned char Taltab[256],Tal1tab[256]; + +static unsigned char tal[] = { 0x8d, 0xef, 0xec, 0x86, 0xfa, 0x99, 0xaf, 0x7b }; + +/* Generate conversion lookup tables between conventional alpha representation + * (@**7, @**6, ...@**0) + * and Berlekamp's dual basis representation + * (l0, l1, ...l7) + */ +int main(){ + int i,j,k; + + 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); + } + Tal1tab[Taltab[i]] = i; + } + printf("unsigned char Taltab[] = {\n"); + for(i=0;i<256;i++){ + if((i % 16) == 0) + printf("\n"); + printf("0x%02x,",Taltab[i]); + } + printf("\n};\n\nunsigned char Tal1tab[] = {"); + for(i=0;i<256;i++){ + if((i % 16) == 0) + printf("\n"); + printf("0x%02x,",Tal1tab[i]); + } + printf("\n};\n"); + exit(0); +} + diff --git a/gr-fec/lib/reed-solomon/init_rs.c b/gr-fec/lib/reed-solomon/init_rs.c new file mode 100644 index 0000000000..4ec77cd723 --- /dev/null +++ b/gr-fec/lib/reed-solomon/init_rs.c @@ -0,0 +1,129 @@ +/* Initialize a RS codec + * + * Copyright 2002 Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#include <stdlib.h> + +#ifdef CCSDS +#include "ccsds.h" +#elif defined(BIGSYM) +#include "int.h" +#else +#include "char.h" +#endif + +#define NULL ((void *)0) + +void FREE_RS(void *p){ + struct rs *rs = (struct rs *)p; + + free(rs->alpha_to); + free(rs->index_of); + free(rs->genpoly); + free(rs); +} + +/* Initialize a Reed-Solomon codec + * symsize = symbol size, bits (1-8) + * gfpoly = Field generator polynomial coefficients + * fcr = first root of RS code generator polynomial, index form + * prim = primitive element to generate polynomial roots + * nroots = RS code generator polynomial degree (number of roots) + */ +void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned fcr,unsigned prim, + unsigned int nroots){ + struct rs *rs; + int sr,root,iprim; + unsigned int i, j; + + if(symsize > 8*sizeof(DTYPE)) + return NULL; /* Need version with ints rather than chars */ + + if(fcr >= (1u<<symsize)) + return NULL; + if(prim == 0 || prim >= (1u<<symsize)) + return NULL; + if(nroots >= (1u<<symsize)) + return NULL; /* Can't have more roots than symbol values! */ + + rs = (struct rs *)calloc(1,sizeof(struct rs)); + rs->mm = symsize; + rs->nn = (1<<symsize)-1; + + rs->alpha_to = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1)); + if(rs->alpha_to == NULL){ + free(rs); + return NULL; + } + rs->index_of = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1)); + if(rs->index_of == NULL){ + free(rs->alpha_to); + free(rs); + return NULL; + } + + /* Generate Galois field lookup tables */ + rs->index_of[0] = A0; /* log(zero) = -inf */ + rs->alpha_to[A0] = 0; /* alpha**-inf = 0 */ + sr = 1; + for(i=0;i<rs->nn;i++){ + rs->index_of[sr] = i; + rs->alpha_to[i] = sr; + sr <<= 1; + if(sr & (1<<symsize)) + sr ^= gfpoly; + sr &= rs->nn; + } + if(sr != 1){ + /* field generator polynomial is not primitive! */ + free(rs->alpha_to); + free(rs->index_of); + free(rs); + return NULL; + } + + /* Form RS code generator polynomial from its roots */ + rs->genpoly = (DTYPE *)malloc(sizeof(DTYPE)*(nroots+1)); + if(rs->genpoly == NULL){ + free(rs->alpha_to); + free(rs->index_of); + free(rs); + return NULL; + } + rs->fcr = fcr; + rs->prim = prim; + rs->nroots = nroots; + + /* Find prim-th root of 1, used in decoding */ + for(iprim=1;(iprim % prim) != 0;iprim += rs->nn) + ; + rs->iprim = iprim / prim; + + rs->genpoly[0] = 1; + for (i = 0,root=fcr*prim; i < nroots; i++,root += prim) { + rs->genpoly[i+1] = 1; + + /* Multiply rs->genpoly[] by @**(root + x) */ + for (j = i; j > 0; j--){ + if (rs->genpoly[j] != 0) + rs->genpoly[j] = rs->genpoly[j-1] ^ rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[j]] + root)]; + else + rs->genpoly[j] = rs->genpoly[j-1]; + } + /* rs->genpoly[0] can never be zero */ + rs->genpoly[0] = rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[0]] + root)]; + } + /* convert rs->genpoly[] to index form for quicker encoding */ + for (i = 0; i <= nroots; i++) + rs->genpoly[i] = rs->index_of[rs->genpoly[i]]; + +#if 0 + printf ("genpoly:\n"); + for (i = nroots; i >= 0; i--){ + printf (" %3d*X^%d\n", rs->alpha_to[rs->genpoly[i]], i); + } +#endif + + return rs; +} diff --git a/gr-fec/lib/reed-solomon/int.h b/gr-fec/lib/reed-solomon/int.h new file mode 100644 index 0000000000..6bb00fd102 --- /dev/null +++ b/gr-fec/lib/reed-solomon/int.h @@ -0,0 +1,55 @@ +/* Include file to configure the RS codec for integer symbols + * + * Copyright 2002, Phil Karn, KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ +#define DTYPE int + +#include <fec/api.h> + +/* Reed-Solomon codec control block */ +struct FEC_API rs { + unsigned int mm; /* Bits per symbol */ + unsigned int nn; /* Symbols per block (= (1<<mm)-1) */ + int *alpha_to; /* log lookup table */ + int *index_of; /* Antilog lookup table */ + int *genpoly; /* Generator polynomial */ + unsigned int nroots; /* Number of generator roots = number of parity symbols */ + unsigned int fcr; /* First consecutive root, index form */ + unsigned int prim; /* Primitive element, index form */ + unsigned int iprim; /* prim-th root of 1, index form */ +}; + +static inline int modnn(struct rs *rs,int x){ + while (x >= rs->nn) { + x -= rs->nn; + x = (x >> rs->mm) + (x & rs->nn); + } + return x; +} +#define MODNN(x) modnn(rs,x) + +#define MM (rs->mm) +#define NN (rs->nn) +#define ALPHA_TO (rs->alpha_to) +#define INDEX_OF (rs->index_of) +#define GENPOLY (rs->genpoly) +#define NROOTS (rs->nroots) +#define FCR (rs->fcr) +#define PRIM (rs->prim) +#define IPRIM (rs->iprim) +#define A0 (NN) + +#define ENCODE_RS encode_rs_int +#define DECODE_RS decode_rs_int +#define INIT_RS init_rs_int +#define FREE_RS free_rs_int + +FEC_API void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity); +FEC_API int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras); +void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned int fcr, + unsigned int prim,unsigned int nroots); +FEC_API void FREE_RS(void *p); + + + diff --git a/gr-fec/lib/reed-solomon/rs.3 b/gr-fec/lib/reed-solomon/rs.3 new file mode 100644 index 0000000000..c3953ce57a --- /dev/null +++ b/gr-fec/lib/reed-solomon/rs.3 @@ -0,0 +1,170 @@ +.TH REED-SOLOMON 3 +.SH NAME +init_rs_int, encode_rs_int, decode_rs_int, free_rs_int, +init_rs_char, encode_rs_char, decode_rs_char, free_rs_char, +encode_rs_8, decode_rs_8, encode_rs_ccsds, decode_rs_ccsds +.SH SYNOPSIS +.nf +.ft B +#include "rs.h" + +void *init_rs_int(unsigned int symsize,unsigned int gfpoly,unsigned fcr, +unsigned prim,unsigned int nroots); +void encode_rs_int(void *rs,int *data,int *parity); +int decode_rs_int(void *rs,int *data,int *eras_pos,int no_eras); +void free_rs_int(void *rs); + +void *init_rs_char(unsigned int symsize,unsigned int gfpoly,unsigned fcr, +unsigned prim,unsigned int nroots); +void encode_rs_char(void *rs,unsigned char *data,unsigned char *parity); +int decode_rs_char(void *rs,unsigned char *data,int *eras_pos,int no_eras); +void free_rs_char(void *rs); + +void encode_rs_8(unsigned char *data,unsigned char *parity); +int decode_rs_8(unsigned char *data,int *eras_pos,int no_eras); + +void encode_rs_ccsds(unsigned char *data,unsigned char *parity); +int decode_rs_ccsds(unsigned char *data,int *eras_pos,int no_eras); + +unsigned char Taltab[256]; +unsigned char Tal1tab[256]; + +.fi + +.SH DESCRIPTION +These functions implement Reed-Solomon error control encoding and +decoding. For optimal performance in a variety of applications, three +sets of functions are supplied. To access these functions, add "-lrs" +to your linker command line. + +The functions with names ending in "_int" handle data in integer arrays, +permitting arbitrarily large codewords limited only by machine +resources. + +The functions with names ending in "_char" take unsigned char arrays and can +handle codes with symbols of 8 bits or less (i.e., with codewords of +255 symbols or less). + +\fBencode_rs_8\fR and \fBdecode_rs_8\fR implement a specific +(255,223) code with 8-bit symbols specified by the CCSDS: +a field generator of 1 + X + X^2 + X^7 + X^8 and a code +generator with first consecutive root = 112 and a primitive element of +11. These functions use the conventional +polynomial form, \fBnot\fR the dual-basis specified in +the CCSDS standard, to represent symbols. + +For full CCSDS compatibility, \fBencode_rs_ccsds\fR and +\fBdecode_rs_ccsds\fR are provided. These functions use two lookup +tables, \fBTaltab\fR to convert from conventional to dual-basis, and +\fBTal1tab\fR to perform the inverse mapping from dual-basis to +conventional form, before and after calls to \fBencode_rs_8\fR +and \fBdecode_rs_8\fR. + +The _8 and _ccsds functions do not require initialization. +To use the general purpose RS encoder or decoder (i.e., +the _char or _int versions), the user must first +call \fBinit_rs_int\fR or \fBinit_rs_char\fR as appropriate. The +arguments are as follows: + +\fBsymsize\fR gives the symbol size in bits, up to 8 for \fBinit_rs_char\fR +or 32 for \fBinit_rs_int\fR on a machine with 32-bit ints (though such a +huge code would exhaust memory limits on a 32-bit machine). The resulting +Reed-Solomon code word will have 2^\fBsymsize\fR - 1 symbols, +each containing \fBsymsize\fR bits. + +\fBgfpoly\fR gives the extended Galois field generator polynomial coefficients, +with the 0th coefficient in the low order bit. The polynomial +\fImust\fR be primitive; if not, the call will fail and NULL will be +returned. + +\fBfcr\fR gives, in index form, the first consecutive root of the +Reed Solomon code generator polynomial. + +\fBprim\fR gives, in index form, the primitive element in the Galois field +used to generate the Reed Solomon code generator polynomial. + +\fBnroots\fR gives the number of roots in the Reed Solomon code +generator polynomial. This equals the number of parity symbols +per code block. + +The resulting Reed-Solomon code has parameters (N,K), where +N = 2^\fBsymsize\fR-1 and K = N-\fBnroots\fR. + +The \fBencode_rs_char\fR and \fBencode_rs_int\fR functions accept +the pointer returned by \fBinit_rs_char\fR or +\fBinit_rs_int\fR, respectively, to +encode a block of data using the specified code. +The input data array is expected to +contain K symbols (of \fBsymsize\fR bits each, right justified +in each char or int) and \fBnroots\fR parity symbols will be placed +into the \fBparity\fR array, right justified. + +The \fBdecode_rs_char\fR and \fBdecode_rs_int\fR functions correct +the errors in a Reed-Solomon codeword up to the capability of the code. +An optional list of "erased" symbol indices may be given in the \fBeras_pos\fR +array to assist the decoder; this parameter may be NULL if no erasures +are given. The number of erased symbols must be given in the \fBno_eras\fR +parameter. + +To maximize performance, the encode and decode functions perform no +"sanity checking" of their inputs. Decoder failure may result if +\fBeras_pos\fR contains duplicate entries, and both encoder and +decoder will fail if an input symbol exceeds its allowable range. +(Symbol range overflow cannot occur with the _8 or _ccsds functions, +or with the _char functions when 8-bit symbols are specified.) + +The decoder corrects the symbols "in place", returning the number +of symbols in error. If the codeword is uncorrectable, -1 is returned +and the data block is unchanged. If \fBeras_pos\fR is non-null, it is +used to return a list of corrected symbol positions, in no particular +order. This means that the +array passed through this parameter \fImust\fR have at least \fBnroots\fR +elements to prevent a possible buffer overflow. + +The \fBfree_rs_int\fR and \fBfree_rs_char\fR functions free the internal +space allocated by the \fBinit_rs_int\fR and \fBinit_rs_char\fR functions, +respecitively. + +The functions \fBencode_rs_8\fR and \fBdecode_rs_8\fR do not have +corresponding \fBinit\fR and \fBfree\fR, nor do they take the +\fBrs\fR argument accepted by the other functions as their parameters +are statically compiled. These functions implement a code +equivalent to calling + +\fBinit_rs_char\fR(8,0x187,112,11,32); + +and using the resulting pointer with \fBencode_rs_char\fR and +\fBdecode_rs_char\fR. + +.SH RETURN VALUES +\fBinit_rs_int\fR and \fBinit_rs_char\fR return a pointer to an internal +control structure that must be passed to the corresponding encode, decode +and free functions. These functions return NULL on error. + +The decode functions return a count of corrected +symbols, or -1 if the block was uncorrectible. + +.SH AUTHOR +Phil Karn, KA9Q (karn@ka9q.net), based heavily on earlier work by Robert +Morelos-Zaragoza (rober@spectra.eng.hawaii.edu) and Hari Thirumoorthy +(harit@spectra.eng.hawaii.edu). + +.SH COPYRIGHT +Copyright 2002, Phil Karn, KA9Q. May be used under the terms of the +GNU General Public License (GPL). + +.SH SEE ALSO +CCSDS 101.0-B-5: Telemetry Channel Coding. +http://www.ccsds.org/documents/pdf/CCSDS-101.0-B-5.pdf + +.SH NOTE +CCSDS chose the "dual basis" symbol representation because it +simplified the implementation of a Reed-Solomon encoder in dedicated +hardware. However, this approach holds no advantages for a software +implementation on a general purpose computer, so use of the dual basis +is recommended only if compatibility with the CCSDS standard is needed, +e.g., to decode data from an existing spacecraft using the CCSDS +standard. If you just want a fast (255,223) RS codec without needing +to interoperate with a CCSDS standard code, use \fBencode_rs_8\fR +and \fBdecode_rs_8\fR. + diff --git a/gr-fec/lib/reed-solomon/rs.h b/gr-fec/lib/reed-solomon/rs.h new file mode 100644 index 0000000000..be7ed51b42 --- /dev/null +++ b/gr-fec/lib/reed-solomon/rs.h @@ -0,0 +1,31 @@ +#include <fec/api.h> +/* User include file for the Reed-Solomon codec + * Copyright 2002, Phil Karn KA9Q + * May be used under the terms of the GNU General Public License (GPL) + */ + +/* General purpose RS codec, 8-bit symbols */ +FEC_API void encode_rs_char(void *rs,unsigned char *data,unsigned char *parity); +FEC_API int decode_rs_char(void *rs,unsigned char *data,int *eras_pos, + int no_eras); +FEC_API void *init_rs_char(unsigned int symsize,unsigned int gfpoly, + unsigned int fcr,unsigned int prim,unsigned int nroots); +FEC_API void free_rs_char(void *rs); + +/* General purpose RS codec, integer symbols */ +FEC_API void encode_rs_int(void *rs,int *data,int *parity); +FEC_API int decode_rs_int(void *rs,int *data,int *eras_pos,int no_eras); +FEC_API void *init_rs_int(unsigned int symsize,unsigned int gfpoly,unsigned int fcr, + unsigned int prim,unsigned int nroots); +FEC_API void free_rs_int(void *rs); + +/* CCSDS standard (255,223) RS codec with conventional (*not* dual-basis) + * symbol representation + */ +FEC_API void encode_rs_8(unsigned char *data,unsigned char *parity); +FEC_API int decode_rs_8(unsigned char *data,int *eras_pos,int no_eras); + +/* Tables to map from conventional->dual (Taltab) and + * dual->conventional (Tal1tab) bases + */ +extern unsigned char Taltab[],Tal1tab[]; diff --git a/gr-fec/lib/reed-solomon/rstest.c b/gr-fec/lib/reed-solomon/rstest.c new file mode 100644 index 0000000000..d8fc5448a7 --- /dev/null +++ b/gr-fec/lib/reed-solomon/rstest.c @@ -0,0 +1,117 @@ +/* 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 General Public License (GPL) + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdio.h> +#include <stdlib.h> +#include <time.h> +#include "rs.h" + +int exercise_char(void *,int); + +#ifdef ALL_VERSIONS +int exercise_int(void *,int); +int exercise_8(int); +int exercise_ccsds(int); +#endif + +struct { + 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 */ +#ifdef ALL_VESIONS + {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 }, +#endif + {0, 0, 0, 0, 0}, +}; + +int main(){ + void *handle; + int errs,terrs; + int i; + + terrs = 0; + srandom(time(NULL)); + +#ifdef ALL_VERSIONS + printf("Testing fixed (255,223) RS codec..."); + fflush(stdout); + errs = exercise_8(10); + terrs += errs; + if(errs == 0){ + printf("OK\n"); + } + printf("Testing CCSDS standard (255,223) RS codec..."); + fflush(stdout); + errs = exercise_ccsds(10); + terrs += errs; + if(errs == 0){ + printf("OK\n"); + } +#endif + + 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) RS codec...",nn,kk); + fflush(stdout); + if(Tab[i].symsize <= 8){ + if((handle = init_rs_char(Tab[i].symsize,Tab[i].genpoly,Tab[i].fcs,Tab[i].prim,Tab[i].nroots)) == NULL){ + printf("init_rs_char failed!\n"); + continue; + } + errs = exercise_char(handle,Tab[i].ntrials); + } else { +#ifdef ALL_VERSIONS + if((handle = init_rs_int(Tab[i].symsize,Tab[i].genpoly,Tab[i].fcs,Tab[i].prim,Tab[i].nroots)) == NULL){ + printf("init_rs_int failed!\n"); + continue; + } + errs = exercise_int(handle,Tab[i].ntrials); +#else + printf ("init_rs_init support is not enabled\n"); + exit (1); +#endif + + } + terrs += errs; + if(errs == 0){ + printf("OK\n"); + } + free_rs_char(handle); + } + if(terrs == 0) + printf("All codec tests passed!\n"); + + exit(0); +} + + diff --git a/gr-fec/lib/viterbi/CMakeLists.txt b/gr-fec/lib/viterbi/CMakeLists.txt new file mode 100644 index 0000000000..6605556a84 --- /dev/null +++ b/gr-fec/lib/viterbi/CMakeLists.txt @@ -0,0 +1,63 @@ +# Copyright 2010-2011 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# 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. + +######################################################################## +# This file included, use CMake directory variables +######################################################################## + +set(viterbi_sources + ${CMAKE_CURRENT_SOURCE_DIR}/metrics.c + ${CMAKE_CURRENT_SOURCE_DIR}/tab.c + ${CMAKE_CURRENT_SOURCE_DIR}/viterbi.c +) + +######################################################################## +# define missing erf function with C linkage (hack for metrics.c) +######################################################################## +if(MSVC) +file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/boost_math_erf.cc " +#include <boost/math/special_functions/erf.hpp> +extern \"C\" double erf(double x){ + return boost::math::erf(x); +} +") +list(APPEND viterbi_sources ${CMAKE_CURRENT_BINARY_DIR}/boost_math_erf.cc) +endif(MSVC) + +######################################################################## +# Append gnuradio-core library sources +######################################################################## +list(APPEND gnuradio_fec_sources ${viterbi_sources}) + +######################################################################## +# Install runtime headers +######################################################################## +#install( +# FILES ${CMAKE_CURRENT_SOURCE_DIR}/viterbi.h +# DESTINATION ${GR_INCLUDE_DIR}/gnuradio/fec +# COMPONENT "fec_devel" +#) + +######################################################################## +# Create some text executables (not registered tests) +# Its not much to build so the sources are just re-listed, +# rather than create a new library just for these two apps. +######################################################################## +#ADD_EXECUTABLE(viterbi_encode ${CMAKE_CURRENT_SOURCE_DIR}/encode.cc ${viterbi_sources}) +#ADD_EXECUTABLE(viterbi_decode ${CMAKE_CURRENT_SOURCE_DIR}/decode.cc ${viterbi_sources}) diff --git a/gr-fec/lib/viterbi/decode.cc b/gr-fec/lib/viterbi/decode.cc new file mode 100644 index 0000000000..368e697134 --- /dev/null +++ b/gr-fec/lib/viterbi/decode.cc @@ -0,0 +1,88 @@ +/* -*- c++ -*- */ +/* + * Copyright 2008 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +/* + * This is a minimal example demonstrating how to call the Viterbi decoder + * in continuous streaming mode. It accepts data on stdin and writes to + * stdout. + * + */ + +extern "C" { +#include "viterbi.h" +} + +#include <cstdio> +#include <cmath> + +#define MAXCHUNKSIZE 4096 +#define MAXENCSIZE MAXCHUNKSIZE*16 + +int main() +{ + unsigned char data[MAXCHUNKSIZE]; + signed char syms[MAXENCSIZE]; + int count = 0; + + // Initialize metric table + int mettab[2][256]; + int amp = 100; + float RATE=0.5; + float ebn0 = 12.0; + float esn0 = RATE*pow(10.0, ebn0/10); + gen_met(mettab, amp, esn0, 0.0, 4); + + // Initialize decoder state + struct viterbi_state state0[64]; + struct viterbi_state state1[64]; + unsigned char viterbi_in[16]; + viterbi_chunks_init(state0); + + while (!feof(stdin)) { + unsigned int n = fread(syms, 1, MAXENCSIZE, stdin); + unsigned char *out = data; + + for (unsigned int i = 0; i < n; i++) { + + // FIXME: This implements hard decoding by slicing the input stream + unsigned char sym = syms[i] > 0 ? -amp : amp; + + // Write the symbol to the decoder input + viterbi_in[count % 4] = sym; + + // Every four symbols, perform the butterfly2 operation + if ((count % 4) == 3) { + viterbi_butterfly2(viterbi_in, mettab, state0, state1); + + // Every sixteen symbols, perform the readback operation + if ((count > 64) && (count % 16) == 11) { + viterbi_get_output(state0, out); + fwrite(out++, 1, 1, stdout); + } + } + + count++; + } + } + + return 0; +} diff --git a/gr-fec/lib/viterbi/encode.cc b/gr-fec/lib/viterbi/encode.cc new file mode 100644 index 0000000000..83a85fcacb --- /dev/null +++ b/gr-fec/lib/viterbi/encode.cc @@ -0,0 +1,54 @@ +/* -*- c++ -*- */ +/* + * Copyright 2008 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +/* + * This is a minimal example demonstrating how to call the ECC encoder + * in continuous streaming mode. It accepts data on stdin and writes to + * stdout. + * + * FIXME: This does not flush the final bits out of the encoder. + * + */ + +extern "C" { +#include "viterbi.h" +} + +#include <cstdio> + +#define MAXCHUNKSIZE 4096 +#define MAXENCSIZE MAXCHUNKSIZE*16 + +int main() +{ + unsigned char encoder_state = 0; + unsigned char data[MAXCHUNKSIZE]; + unsigned char syms[MAXENCSIZE]; + + while (!feof(stdin)) { + unsigned int n = fread(data, 1, MAXCHUNKSIZE, stdin); + encoder_state = encode(syms, data, n, encoder_state); + fwrite(syms, 1, n*16, stdout); + } + + return 0; +} diff --git a/gr-fec/lib/viterbi/metrics.c b/gr-fec/lib/viterbi/metrics.c new file mode 100644 index 0000000000..0d91c301ff --- /dev/null +++ b/gr-fec/lib/viterbi/metrics.c @@ -0,0 +1,126 @@ +/* + * Copyright 1995 Phil Karn, KA9Q + * Copyright 2008 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +/* + * Generate metric tables for a soft-decision convolutional decoder + * assuming gaussian noise on a PSK channel. + * + * Works from "first principles" by evaluating the normal probability + * function and then computing the log-likelihood function + * for every possible received symbol value + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +/* Symbols are offset-binary, with 128 corresponding to an erased (no + * information) symbol + */ +#define OFFSET 128 + +#include <stdlib.h> +#include <math.h> + +//declare erf in case it was missing in math.h and provided for by the build system +extern double erf(double x); + +/* Normal function integrated from -Inf to x. Range: 0-1 */ +#define normal(x) (0.5 + 0.5*erf((x)/M_SQRT2)) + +/* Logarithm base 2 */ +#define gr_log2(x) (log(x)*M_LOG2E) + +/* Generate log-likelihood metrics for 8-bit soft quantized channel + * assuming AWGN and BPSK + */ +void +gen_met(int mettab[2][256], /* Metric table, [sent sym][rx symbol] */ + int amp, /* Signal amplitude, units */ + double esn0, /* Es/N0 ratio in dB */ + double bias, /* Metric bias; 0 for viterbi, rate for sequential */ + int scale) /* Scale factor */ +{ + double noise; + int s,bit; + double metrics[2][256]; + double p0,p1; + + /* Es/N0 as power ratio */ + esn0 = pow(10.,esn0/10); + + noise = 0.5/esn0; /* only half the noise for BPSK */ + noise = sqrt(noise); /* noise/signal Voltage ratio */ + + /* Zero is a special value, since this sample includes all + * lower samples that were clipped to this value, i.e., it + * takes the whole lower tail of the curve + */ + p1 = normal(((0-OFFSET+0.5)/amp - 1)/noise); /* P(s|1) */ + + /* Prob of this value occurring for a 0-bit */ /* P(s|0) */ + p0 = normal(((0-OFFSET+0.5)/amp + 1)/noise); + metrics[0][0] = gr_log2(2*p0/(p1+p0)) - bias; + metrics[1][0] = gr_log2(2*p1/(p1+p0)) - bias; + + for(s=1;s<255;s++){ + /* P(s|1), prob of receiving s given 1 transmitted */ + p1 = normal(((s-OFFSET+0.5)/amp - 1)/noise) - + normal(((s-OFFSET-0.5)/amp - 1)/noise); + + /* P(s|0), prob of receiving s given 0 transmitted */ + p0 = normal(((s-OFFSET+0.5)/amp + 1)/noise) - + normal(((s-OFFSET-0.5)/amp + 1)/noise); + +#ifdef notdef + printf("P(%d|1) = %lg, P(%d|0) = %lg\n",s,p1,s,p0); +#endif + metrics[0][s] = gr_log2(2*p0/(p1+p0)) - bias; + metrics[1][s] = gr_log2(2*p1/(p1+p0)) - bias; + } + /* 255 is also a special value */ + /* P(s|1) */ + p1 = 1 - normal(((255-OFFSET-0.5)/amp - 1)/noise); + /* P(s|0) */ + p0 = 1 - normal(((255-OFFSET-0.5)/amp + 1)/noise); + + metrics[0][255] = gr_log2(2*p0/(p1+p0)) - bias; + metrics[1][255] = gr_log2(2*p1/(p1+p0)) - bias; +#ifdef notdef + /* The probability of a raw symbol error is the probability + * that a 1-bit would be received as a sample with value + * 0-128. This is the offset normal curve integrated from -Inf to 0. + */ + printf("symbol Pe = %lg\n",normal(-1/noise)); +#endif + for(bit=0;bit<2;bit++){ + for(s=0;s<256;s++){ + /* Scale and round to nearest integer */ + mettab[bit][s] = floor(metrics[bit][s] * scale + 0.5); +#ifdef notdef + printf("metrics[%d][%d] = %lg, mettab = %d\n", + bit,s,metrics[bit][s],mettab[bit][s]); +#endif + } + } +} diff --git a/gr-fec/lib/viterbi/tab.c b/gr-fec/lib/viterbi/tab.c new file mode 100644 index 0000000000..1c135acfee --- /dev/null +++ b/gr-fec/lib/viterbi/tab.c @@ -0,0 +1,57 @@ +/* + * Copyright 1995 Phil Karn, KA9Q + * Copyright 2008 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +/* 8-bit parity lookup table, generated by partab.c */ +unsigned char Partab[] = { + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, + 0, 1, 1, 0, 1, 0, 0, 1, + 0, 1, 1, 0, 1, 0, 0, 1, + 1, 0, 0, 1, 0, 1, 1, 0, +}; diff --git a/gr-fec/lib/viterbi/viterbi.c b/gr-fec/lib/viterbi/viterbi.c new file mode 100644 index 0000000000..fc88866035 --- /dev/null +++ b/gr-fec/lib/viterbi/viterbi.c @@ -0,0 +1,355 @@ +/* + * Copyright 1995 Phil Karn, KA9Q + * Copyright 2008 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +/* + * Viterbi decoder for K=7 rate=1/2 convolutional code + * Some modifications from original Karn code by Matt Ettus + */ + +#include "viterbi.h" + +/* The two generator polynomials for the NASA Standard K=7 code. + * Since these polynomials are known to be optimal for this constraint + * length there is not much point in changing them. But if you do, you + * will have to regenerate the BUTTERFLY macro calls in viterbi() + */ +#define POLYA 0x6d +#define POLYB 0x4f + +/* The basic Viterbi decoder operation, called a "butterfly" + * operation because of the way it looks on a trellis diagram. Each + * butterfly involves an Add-Compare-Select (ACS) operation on the two nodes + * where the 0 and 1 paths from the current node merge at the next step of + * the trellis. + * + * The code polynomials are assumed to have 1's on both ends. Given a + * function encode_state() that returns the two symbols for a given + * encoder state in the low two bits, such a code will have the following + * identities for even 'n' < 64: + * + * encode_state(n) = encode_state(n+65) + * encode_state(n+1) = encode_state(n+64) = (3 ^ encode_state(n)) + * + * Any convolutional code you would actually want to use will have + * these properties, so these assumptions aren't too limiting. + * + * Doing this as a macro lets the compiler evaluate at compile time the + * many expressions that depend on the loop index and encoder state and + * emit them as immediate arguments. + * This makes an enormous difference on register-starved machines such + * as the Intel x86 family where evaluating these expressions at runtime + * would spill over into memory. + */ +#define BUTTERFLY(i,sym) { \ + int m0,m1;\ +\ + /* ACS for 0 branch */\ + m0 = state[i].metric + mets[sym]; /* 2*i */\ + m1 = state[i+32].metric + mets[3^sym]; /* 2*i + 64 */\ + if(m0 > m1){\ + next[2*i].metric = m0;\ + next[2*i].path = state[i].path << 1;\ + } else {\ + next[2*i].metric = m1;\ + next[2*i].path = (state[i+32].path << 1)|1;\ + }\ + /* ACS for 1 branch */\ + m0 = state[i].metric + mets[3^sym]; /* 2*i + 1 */\ + m1 = state[i+32].metric + mets[sym]; /* 2*i + 65 */\ + if(m0 > m1){\ + next[2*i+1].metric = m0;\ + next[2*i+1].path = state[i].path << 1;\ + } else {\ + next[2*i+1].metric = m1;\ + next[2*i+1].path = (state[i+32].path << 1)|1;\ + }\ +} + +extern unsigned char Partab[]; /* Parity lookup table */ + +/* Convolutionally encode data into binary symbols */ +unsigned char +encode(unsigned char *symbols, + unsigned char *data, + unsigned int nbytes, + unsigned char encstate) +{ + int i; + + while(nbytes-- != 0){ + for(i=7;i>=0;i--){ + encstate = (encstate << 1) | ((*data >> i) & 1); + *symbols++ = Partab[encstate & POLYA]; + *symbols++ = Partab[encstate & POLYB]; + } + data++; + } + + return encstate; +} + +/* Viterbi decoder */ +int +viterbi(unsigned long *metric, /* Final path metric (returned value) */ + unsigned char *data, /* Decoded output data */ + unsigned char *symbols, /* Raw deinterleaved input symbols */ + unsigned int nbits, /* Number of output bits */ + int mettab[2][256] /* Metric table, [sent sym][rx symbol] */ + ){ + unsigned int bitcnt = 0; + int mets[4]; + long bestmetric; + int beststate,i; + struct viterbi_state state0[64],state1[64],*state,*next; + + state = state0; + next = state1; + + /* Initialize starting metrics to prefer 0 state */ + state[0].metric = 0; + for(i=1;i<64;i++) + state[i].metric = -999999; + state[0].path = 0; + + for(bitcnt = 0;bitcnt < nbits;bitcnt++){ + /* Read input symbol pair and compute all possible branch + * metrics + */ + mets[0] = mettab[0][symbols[0]] + mettab[0][symbols[1]]; + mets[1] = mettab[0][symbols[0]] + mettab[1][symbols[1]]; + mets[2] = mettab[1][symbols[0]] + mettab[0][symbols[1]]; + mets[3] = mettab[1][symbols[0]] + mettab[1][symbols[1]]; + symbols += 2; + + /* These macro calls were generated by genbut.c */ + BUTTERFLY(0,0); + BUTTERFLY(1,1); + BUTTERFLY(2,3); + BUTTERFLY(3,2); + BUTTERFLY(4,3); + BUTTERFLY(5,2); + BUTTERFLY(6,0); + BUTTERFLY(7,1); + BUTTERFLY(8,0); + BUTTERFLY(9,1); + BUTTERFLY(10,3); + BUTTERFLY(11,2); + BUTTERFLY(12,3); + BUTTERFLY(13,2); + BUTTERFLY(14,0); + BUTTERFLY(15,1); + BUTTERFLY(16,2); + BUTTERFLY(17,3); + BUTTERFLY(18,1); + BUTTERFLY(19,0); + BUTTERFLY(20,1); + BUTTERFLY(21,0); + BUTTERFLY(22,2); + BUTTERFLY(23,3); + BUTTERFLY(24,2); + BUTTERFLY(25,3); + BUTTERFLY(26,1); + BUTTERFLY(27,0); + BUTTERFLY(28,1); + BUTTERFLY(29,0); + BUTTERFLY(30,2); + BUTTERFLY(31,3); + + /* Swap current and next states */ + if(bitcnt & 1){ + state = state0; + next = state1; + } else { + state = state1; + next = state0; + } + // ETTUS + //if(bitcnt > nbits-7){ + /* In tail, poison non-zero nodes */ + //for(i=1;i<64;i += 2) + // state[i].metric = -9999999; + //} + /* Produce output every 8 bits once path memory is full */ + if((bitcnt % 8) == 5 && bitcnt > 32){ + /* Find current best path */ + bestmetric = state[0].metric; + beststate = 0; + for(i=1;i<64;i++){ + if(state[i].metric > bestmetric){ + bestmetric = state[i].metric; + beststate = i; + } + } +#ifdef notdef + printf("metrics[%d] = %d state = %lx\n",beststate, + state[beststate].metric,state[beststate].path); +#endif + *data++ = state[beststate].path >> 24; + } + + } + /* Output remaining bits from 0 state */ + // ETTUS Find best state instead + bestmetric = state[0].metric; + beststate = 0; + for(i=1;i<64;i++){ + if(state[i].metric > bestmetric){ + bestmetric = state[i].metric; + beststate = i; + } + } + if((i = bitcnt % 8) != 6) + state[beststate].path <<= 6-i; + + *data++ = state[beststate].path >> 24; + *data++ = state[beststate].path >> 16; + *data++ = state[beststate].path >> 8; + *data = state[beststate].path; + //printf ("BS = %d\tBSM = %d\tM0 = %d\n",beststate,state[beststate].metric,state[0].metric); + *metric = state[beststate].metric; + return 0; +} + + +void +viterbi_chunks_init(struct viterbi_state* state) { + // Initialize starting metrics to prefer 0 state + int i; + state[0].metric = 0; + state[0].path = 0; + for(i=1;i<64;i++) + state[i].metric = -999999; +} + +void +viterbi_butterfly8(unsigned char *symbols, int mettab[2][256], struct viterbi_state *state0, struct viterbi_state *state1) +{ + unsigned int bitcnt; + int mets[4]; + + struct viterbi_state *state, *next; + state = state0; + next = state1; + // Operate on 16 symbols (8 bits) at a time + for(bitcnt = 0;bitcnt < 8;bitcnt++){ + // Read input symbol pair and compute all possible branch metrics + mets[0] = mettab[0][symbols[0]] + mettab[0][symbols[1]]; + mets[1] = mettab[0][symbols[0]] + mettab[1][symbols[1]]; + mets[2] = mettab[1][symbols[0]] + mettab[0][symbols[1]]; + mets[3] = mettab[1][symbols[0]] + mettab[1][symbols[1]]; + symbols += 2; + + // These macro calls were generated by genbut.c + BUTTERFLY(0,0);BUTTERFLY(1,1);BUTTERFLY(2,3);BUTTERFLY(3,2); + BUTTERFLY(4,3);BUTTERFLY(5,2);BUTTERFLY(6,0);BUTTERFLY(7,1); + BUTTERFLY(8,0);BUTTERFLY(9,1);BUTTERFLY(10,3);BUTTERFLY(11,2); + BUTTERFLY(12,3);BUTTERFLY(13,2);BUTTERFLY(14,0);BUTTERFLY(15,1); + BUTTERFLY(16,2);BUTTERFLY(17,3);BUTTERFLY(18,1);BUTTERFLY(19,0); + BUTTERFLY(20,1);BUTTERFLY(21,0);BUTTERFLY(22,2);BUTTERFLY(23,3); + BUTTERFLY(24,2);BUTTERFLY(25,3);BUTTERFLY(26,1);BUTTERFLY(27,0); + BUTTERFLY(28,1);BUTTERFLY(29,0);BUTTERFLY(30,2);BUTTERFLY(31,3); + + // Swap current and next states + if(bitcnt & 1){ + state = state0; + next = state1; + } else { + state = state1; + next = state0; + } + } +} + +void +viterbi_butterfly2(unsigned char *symbols, int mettab[2][256], struct viterbi_state *state0, struct viterbi_state *state1) +{ + //unsigned int bitcnt; + int mets[4]; + + struct viterbi_state *state, *next; + state = state0; + next = state1; + // Operate on 4 symbols (2 bits) at a time + + // Read input symbol pair and compute all possible branch metrics + mets[0] = mettab[0][symbols[0]] + mettab[0][symbols[1]]; + mets[1] = mettab[0][symbols[0]] + mettab[1][symbols[1]]; + mets[2] = mettab[1][symbols[0]] + mettab[0][symbols[1]]; + mets[3] = mettab[1][symbols[0]] + mettab[1][symbols[1]]; + + // These macro calls were generated by genbut.c + BUTTERFLY(0,0);BUTTERFLY(1,1);BUTTERFLY(2,3);BUTTERFLY(3,2); + BUTTERFLY(4,3);BUTTERFLY(5,2);BUTTERFLY(6,0);BUTTERFLY(7,1); + BUTTERFLY(8,0);BUTTERFLY(9,1);BUTTERFLY(10,3);BUTTERFLY(11,2); + BUTTERFLY(12,3);BUTTERFLY(13,2);BUTTERFLY(14,0);BUTTERFLY(15,1); + BUTTERFLY(16,2);BUTTERFLY(17,3);BUTTERFLY(18,1);BUTTERFLY(19,0); + BUTTERFLY(20,1);BUTTERFLY(21,0);BUTTERFLY(22,2);BUTTERFLY(23,3); + BUTTERFLY(24,2);BUTTERFLY(25,3);BUTTERFLY(26,1);BUTTERFLY(27,0); + BUTTERFLY(28,1);BUTTERFLY(29,0);BUTTERFLY(30,2);BUTTERFLY(31,3); + + state = state1; + next = state0; + + // Read input symbol pair and compute all possible branch metrics + mets[0] = mettab[0][symbols[2]] + mettab[0][symbols[3]]; + mets[1] = mettab[0][symbols[2]] + mettab[1][symbols[3]]; + mets[2] = mettab[1][symbols[2]] + mettab[0][symbols[3]]; + mets[3] = mettab[1][symbols[2]] + mettab[1][symbols[3]]; + + // These macro calls were generated by genbut.c + BUTTERFLY(0,0);BUTTERFLY(1,1);BUTTERFLY(2,3);BUTTERFLY(3,2); + BUTTERFLY(4,3);BUTTERFLY(5,2);BUTTERFLY(6,0);BUTTERFLY(7,1); + BUTTERFLY(8,0);BUTTERFLY(9,1);BUTTERFLY(10,3);BUTTERFLY(11,2); + BUTTERFLY(12,3);BUTTERFLY(13,2);BUTTERFLY(14,0);BUTTERFLY(15,1); + BUTTERFLY(16,2);BUTTERFLY(17,3);BUTTERFLY(18,1);BUTTERFLY(19,0); + BUTTERFLY(20,1);BUTTERFLY(21,0);BUTTERFLY(22,2);BUTTERFLY(23,3); + BUTTERFLY(24,2);BUTTERFLY(25,3);BUTTERFLY(26,1);BUTTERFLY(27,0); + BUTTERFLY(28,1);BUTTERFLY(29,0);BUTTERFLY(30,2);BUTTERFLY(31,3); +} + +unsigned char +viterbi_get_output(struct viterbi_state *state, unsigned char *outbuf) { + // Produce output every 8 bits once path memory is full + // if((bitcnt % 8) == 5 && bitcnt > 32) { + + // Find current best path + unsigned int i,beststate; + int bestmetric; + + bestmetric = state[0].metric; + beststate = 0; + for(i=1;i<64;i++) + if(state[i].metric > bestmetric) { + bestmetric = state[i].metric; + beststate = i; + } + *outbuf = state[beststate].path >> 24; + return bestmetric; +} + + +//printf ("BS = %d\tBSM = %d\tM0 = %d\n",beststate,state[beststate].metric,state[0].metric); +// In tail, poison non-zero nodes +//if(bits_out > packet_size-7) +// for(i=1;i<64;i += 2) +// state[i].metric = -9999999; + diff --git a/gr-fec/lib/viterbi/viterbi.h b/gr-fec/lib/viterbi/viterbi.h new file mode 100644 index 0000000000..97e09d4a8d --- /dev/null +++ b/gr-fec/lib/viterbi/viterbi.h @@ -0,0 +1,53 @@ +/* + * Copyright 2008 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * 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. + */ + +/* The path memory for each state is 32 bits. This is slightly shorter + * than we'd like for K=7, especially since we chain back every 8 bits. + * But it fits so nicely into a 32-bit machine word... + */ + +#include <fec/api.h> + +struct viterbi_state { + unsigned long path; /* Decoded path to this state */ + long metric; /* Cumulative metric to this state */ +}; + +FEC_API +int gen_met(int mettab[2][256], /* Metric table */ + int amp, /* Signal amplitude */ + double esn0, /* Es/N0 ratio in dB */ + double bias, /* Metric bias */ + int scale); /* Scale factor */ + +FEC_API unsigned char +encode(unsigned char *symbols, unsigned char *data, + unsigned int nbytes,unsigned char encstate); + +FEC_API void +viterbi_chunks_init(struct viterbi_state* state); + + FEC_API void +viterbi_butterfly2(unsigned char *symbols, int mettab[2][256], + struct viterbi_state *state0, struct viterbi_state *state1); + +FEC_API unsigned char +viterbi_get_output(struct viterbi_state *state, unsigned char *outbuf); |