Initial commit

CMakeLists.txt

@@ -0,0 +1,90 @@
+cmake_minimum_required(VERSION 3.0)
+project(horus C)
+
+include(GNUInstallDirs)
+mark_as_advanced(CLEAR
+    CMAKE_INSTALL_BINDIR
+    CMAKE_INSTALL_INCLUDEDIR
+    CMAKE_INSTALL_LIBDIR
+)
+
+#
+# Prevent in-source builds
+# If an in-source build is attempted, you will still need to clean up a few
+# files manually.
+#
+set(CMAKE_DISABLE_SOURCE_CHANGES ON)
+set(CMAKE_DISABLE_IN_SOURCE_BUILD ON)
+if("${CMAKE_SOURCE_DIR}" STREQUAL "${CMAKE_BINARY_DIR}")
+  message(FATAL_ERROR "In-source builds in ${CMAKE_BINARY_DIR} are not "
+   "allowed, please remove ./CMakeCache.txt and ./CMakeFiles/, create a "
+   "separate build directory and run cmake from there.")
+endif("${CMAKE_SOURCE_DIR}" STREQUAL "${CMAKE_BINARY_DIR}")
+
+
+# Set project version information. This should probably be done via external
+# file at some point.
+#
+set(HORUS_VERSION_MAJOR 0)
+set(HORUS_VERSION_MINOR 2)
+# Set to patch level if needed, otherwise leave FALSE.
+# Must be positive (non-zero) if set, since 0 == FALSE in CMake.
+set(HORUS_VERSION_PATCH FALSE)
+set(HORUS_VERSION "${HORUS_VERSION_MAJOR}.${HORUS_VERSION_MINOR}")
+# Patch level version bumps should not change API/ABI.
+set(SOVERSION "${HORUS_VERSION_MAJOR}.${HORUS_VERSION_MINOR}")
+if(HORUS_VERSION_PATCH)
+    set(HORUS_VERSION "${HORUS_VERSION}.${HORUS_VERSION_PATCH}")
+endif()
+message(STATUS "Horuslib version: ${HORUS_VERSION}")
+
+# Set default flags
+set(CMAKE_C_FLAGS "-Wall -Wextra -Wno-unused-function -Wno-strict-overflow -O3 -g -I. -MD ${CMAKE_C_FLAGS} -DENABLE_ASSERTIONS")
+
+# Arch specific stuff here
+message(STATUS "Host system arch is: ${CMAKE_SYSTEM_PROCESSOR}")
+
+
+add_subdirectory(src)
+
+# Ctests ----------------------------------------------------------------------
+
+include(CTest)
+enable_testing()
+
+     add_test(NAME test_horus_binary
+             COMMAND sh -c "cd ${CMAKE_CURRENT_BINARY_DIR}/src;
+                            sox ${CMAKE_CURRENT_SOURCE_DIR}/samples/horus_binary_ebno_4.5db.wav -r 48000 -t raw - |
+                            ./horus_demod -m binary - -"
+              )
+             set_tests_properties(test_horus_binary PROPERTIES PASS_REGULAR_EXPRESSION "1C9A9545")
+
+     add_test(NAME test_horus_binary_iq
+             COMMAND sh -c "cd ${CMAKE_CURRENT_BINARY_DIR}/src;
+                            cat ${CMAKE_CURRENT_SOURCE_DIR}/samples/horusb_iq_s16.raw | 
+                            ./horus_demod -q -m binary --fsk_lower=1000 --fsk_upper=20000 - -"
+             )
+             set_tests_properties(test_horus_binary_iq PROPERTIES
+             PASS_REGULAR_EXPRESSION "000900071E2A000000000000000000000000259A6B14")
+    
+     # Wenet - Using Mask estimator (Assuming ~120 kHz tone spacing)
+     add_test(NAME test_wenet_mask
+        COMMAND sh -c "cd ${CMAKE_CURRENT_BINARY_DIR}/src;
+                       cat ${CMAKE_CURRENT_SOURCE_DIR}/samples/wenet_sample.c8 | 
+                       ./fsk_demod --cu8 -s --mask=120000 2 921416 115177 - - |
+                       ./drs232_ldpc - - 2>&1 | strings"
+        )
+        set_tests_properties(test_wenet_mask PROPERTIES
+        PASS_REGULAR_EXPRESSION "packet_errors: 0 PER: 0.000")
+
+     # Using regular frequency estimator, tell freq est to avoid first 100 kHz due to a DC line 
+     add_test(NAME test_wenet_nomask
+        COMMAND sh -c "cd ${CMAKE_CURRENT_BINARY_DIR}/src; 
+                       cat ${CMAKE_CURRENT_SOURCE_DIR}/samples/wenet_sample.c8 | 
+                       ./fsk_demod --cu8 -s --fsk_lower 100000 2 921416 115177 - - |
+                       ./drs232_ldpc - - 2>&1 | strings"
+        )
+        set_tests_properties(test_wenet_nomask PROPERTIES
+        PASS_REGULAR_EXPRESSION "packet_errors: 0 PER: 0.000")
+
+

README.md

@@ -1,2 +1,29 @@
-# horuslib
+# High Altitude Balloon (HAB) telemetry library
-High Altitude Balloon Telemetry Library
+
+## Building
+
+```
+$ git clone https://github.com/drowe67/codec2.git
+$ cd codec2 && mkdir build_linux && cd build_linux && cmake ../ && make
+$ cd ~
+$ git clone https://github.com/drowe67/hab.git
+$ cd hab && mkdir build_linux && cd build_linux
+$ cmake -DCODEC2_BUILD_DIR=~/codec2/build_linux ..
+$ make
+```
+
+## Testing
+
+```
+$ cd hab/build_linux
+$ ctest
+```
+
+## Further Reading
+
+   Here are some links to projects and blog posts that use this code:
+
+   1. [Horus Binary](https://github.com/projecthorus/horusbinary) High Altitude Balloon (HAB) telemetry protocol, 3 second updates, works at 7dB lower SNR that RTTY.
+   1. [Testing HAB Telemetry, Horus binary waveform](http://www.rowetel.com/?p=5906)
+   1. [Wenet](https://github.com/projecthorus/wenet) - high speed SSTV images from balloons at the edge of space
+   1. [Wenet High speed SSTV images](http://www.rowetel.com/?p=5344)

src/CMakeLists.txt

@@ -0,0 +1,64 @@
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+set(horus_srcs
+  fsk.c
+  kiss_fft.c
+  kiss_fftr.c
+  mpdecode_core.c
+  H_256_768_22.c
+  golay23.c
+  phi0.c
+  horus_api.c
+  horus_l2.c
+)
+
+add_library(horus SHARED ${horus_srcs})
+target_link_libraries(horus m)
+set_target_properties(horus PROPERTIES
+    PUBLIC_HEADER horus_api.h
+)
+target_include_directories(horus INTERFACE
+    $<INSTALL_INTERFACE:include/horus>
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+)
+install(TARGETS horus EXPORT horus-config
+    LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+    ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+    RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+    PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/horus
+)
+
+install(EXPORT horus-config
+    DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/horus
+)
+# Export libhab target for import into build trees of other projects.
+export(TARGETS horus
+    FILE ${CMAKE_BINARY_DIR}/horus.cmake
+)
+
+add_executable(fsk_mod fsk_mod.c)
+target_link_libraries(fsk_mod m horus ${CMAKE_REQUIRED_LIBRARIES})
+
+add_executable(fsk_demod fsk_demod.c modem_probe.c octave.c)
+target_link_libraries(fsk_demod m horus ${CMAKE_REQUIRED_LIBRARIES})
+
+add_executable(fsk_get_test_bits fsk_get_test_bits.c)
+target_link_libraries(fsk_get_test_bits)
+
+add_executable(fsk_put_test_bits fsk_put_test_bits.c)
+target_link_libraries(fsk_put_test_bits ${CMAKE_REQUIRED_LIBRARIES})
+
+add_executable(drs232 drs232.c)
+target_link_libraries(drs232 m horus ${CMAKE_REQUIRED_LIBRARIES})
+
+add_executable(drs232_ldpc drs232_ldpc.c)
+target_link_libraries(drs232_ldpc m horus ${CMAKE_REQUIRED_LIBRARIES})
+
+add_definitions(-DINTERLEAVER -DSCRAMBLER -DRUN_TIME_TABLES)
+add_executable(horus_gen_test_bits horus_gen_test_bits.c)
+target_link_libraries(horus_gen_test_bits m horus)
+
+add_definitions(-DHORUS_L2_RX -DINTERLEAVER -DSCRAMBLER -DRUN_TIME_TABLES)
+add_executable(horus_demod horus_demod.c horus_api.c horus_l2.c golay23.c fsk.c kiss_fft.c)
+target_link_libraries(horus_demod m horus ${CMAKE_REQUIRED_LIBRARIES})
+

src/H_128_384_23.h

@@ -0,0 +1,19 @@
+/*
+  FILE....: H_128_384_23.h
+
+  Static arrays for LDPC codec H_128_384_23, generated by ldpc_gen_c_h_file.m.
+*/
+
+#define H_128_384_23_NUMBERPARITYBITS 256
+#define H_128_384_23_MAX_ROW_WEIGHT 3
+#define H_128_384_23_CODELENGTH 384
+#define H_128_384_23_NUMBERROWSHCOLS 128
+#define H_128_384_23_MAX_COL_WEIGHT 5
+#define H_128_384_23_DEC_TYPE 0
+#define H_128_384_23_MAX_ITER 100
+
+extern const uint16_t H_128_384_23_H_rows[];
+extern const uint16_t H_128_384_23_H_cols[];
+extern const float H_128_384_23_input[];
+extern const char H_128_384_23_detected_data[];
+

src/H_256_768_22.h

@@ -0,0 +1,19 @@
+/*
+  FILE....: H_256_768_22.h
+
+  Static arrays for LDPC codec H_256_768_22, generated by ldpc_gen_c_h_file.m.
+*/
+
+#define H_256_768_22_NUMBERPARITYBITS 512
+#define H_256_768_22_MAX_ROW_WEIGHT 2
+#define H_256_768_22_CODELENGTH 768
+#define H_256_768_22_NUMBERROWSHCOLS 256
+#define H_256_768_22_MAX_COL_WEIGHT 4
+#define H_256_768_22_DEC_TYPE 0
+#define H_256_768_22_MAX_ITER 100
+
+extern const uint16_t H_256_768_22_H_rows[];
+extern const uint16_t H_256_768_22_H_cols[];
+extern const float H_256_768_22_input[];
+extern const char H_256_768_22_detected_data[];
+

src/_kiss_fft_guts.h

@@ -0,0 +1,164 @@
+/*
+Copyright (c) 2003-2010, Mark Borgerding
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+    * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/* kiss_fft.h
+   defines kiss_fft_scalar as either short or a float type
+   and defines
+   typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */
+#include "kiss_fft.h"
+#include <limits.h>
+
+#define MAXFACTORS 32
+/* e.g. an fft of length 128 has 4 factors
+ as far as kissfft is concerned
+ 4*4*4*2
+ */
+
+struct kiss_fft_state{
+    int nfft;
+    int inverse;
+    int factors[2*MAXFACTORS];
+    kiss_fft_cpx twiddles[1];
+};
+
+/*
+  Explanation of macros dealing with complex math:
+
+   C_MUL(m,a,b)         : m = a*b
+   C_FIXDIV( c , div )  : if a fixed point impl., c /= div. noop otherwise
+   C_SUB( res, a,b)     : res = a - b
+   C_SUBFROM( res , a)  : res -= a
+   C_ADDTO( res , a)    : res += a
+ * */
+#ifdef FIXED_POINT
+#if (FIXED_POINT==32)
+# define FRACBITS 31
+# define SAMPPROD int64_t
+#define SAMP_MAX 2147483647
+#else
+# define FRACBITS 15
+# define SAMPPROD int32_t
+#define SAMP_MAX 32767
+#endif
+
+#define SAMP_MIN -SAMP_MAX
+
+#if defined(CHECK_OVERFLOW)
+#  define CHECK_OVERFLOW_OP(a,op,b)  \
+	if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
+		fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) );  }
+#endif
+
+
+#   define smul(a,b) ( (SAMPPROD)(a)*(b) )
+#   define sround( x )  (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS )
+
+#   define S_MUL(a,b) sround( smul(a,b) )
+
+#   define C_MUL(m,a,b) \
+      do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
+          (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
+
+#   define DIVSCALAR(x,k) \
+	(x) = sround( smul(  x, SAMP_MAX/k ) )
+
+#   define C_FIXDIV(c,div) \
+	do {    DIVSCALAR( (c).r , div);  \
+		DIVSCALAR( (c).i  , div); }while (0)
+
+#   define C_MULBYSCALAR( c, s ) \
+    do{ (c).r =  sround( smul( (c).r , s ) ) ;\
+        (c).i =  sround( smul( (c).i , s ) ) ; }while(0)
+
+#else  /* not FIXED_POINT*/
+
+#   define S_MUL(a,b) ( (a)*(b) )
+#define C_MUL(m,a,b) \
+    do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
+        (m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
+#   define C_FIXDIV(c,div) /* NOOP */
+#   define C_MULBYSCALAR( c, s ) \
+    do{ (c).r *= (s);\
+        (c).i *= (s); }while(0)
+#endif
+
+#ifndef CHECK_OVERFLOW_OP
+#  define CHECK_OVERFLOW_OP(a,op,b) /* noop */
+#endif
+
+#define  C_ADD( res, a,b)\
+    do { \
+	    CHECK_OVERFLOW_OP((a).r,+,(b).r)\
+	    CHECK_OVERFLOW_OP((a).i,+,(b).i)\
+	    (res).r=(a).r+(b).r;  (res).i=(a).i+(b).i; \
+    }while(0)
+#define  C_SUB( res, a,b)\
+    do { \
+	    CHECK_OVERFLOW_OP((a).r,-,(b).r)\
+	    CHECK_OVERFLOW_OP((a).i,-,(b).i)\
+	    (res).r=(a).r-(b).r;  (res).i=(a).i-(b).i; \
+    }while(0)
+#define C_ADDTO( res , a)\
+    do { \
+	    CHECK_OVERFLOW_OP((res).r,+,(a).r)\
+	    CHECK_OVERFLOW_OP((res).i,+,(a).i)\
+	    (res).r += (a).r;  (res).i += (a).i;\
+    }while(0)
+
+#define C_SUBFROM( res , a)\
+    do {\
+	    CHECK_OVERFLOW_OP((res).r,-,(a).r)\
+	    CHECK_OVERFLOW_OP((res).i,-,(a).i)\
+	    (res).r -= (a).r;  (res).i -= (a).i; \
+    }while(0)
+
+
+#ifdef FIXED_POINT
+#  define KISS_FFT_COS(phase)  floorf(.5+SAMP_MAX * cosf (phase))
+#  define KISS_FFT_SIN(phase)  floorf(.5+SAMP_MAX * sinf (phase))
+#  define HALF_OF(x) ((x)>>1)
+#elif defined(USE_SIMD)
+#  define KISS_FFT_COS(phase) _mm_set1_ps( cosf(phase) )
+#  define KISS_FFT_SIN(phase) _mm_set1_ps( sinf(phase) )
+#  define HALF_OF(x) ((x)*_mm_set1_ps(.5))
+#else
+#  define KISS_FFT_COS(phase) (kiss_fft_scalar) cosf(phase)
+#  define KISS_FFT_SIN(phase) (kiss_fft_scalar) sinf(phase)
+#  define HALF_OF(x) ((x)*.5)
+#endif
+
+#define  kf_cexp(x,phase) \
+	do{ \
+		(x)->r = KISS_FFT_COS(phase);\
+		(x)->i = KISS_FFT_SIN(phase);\
+	}while(0)
+
+
+/* a debugging function */
+#define pcpx(c)\
+    fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) )
+
+
+#ifdef KISS_FFT_USE_ALLOCA
+// define this to allow use of alloca instead of malloc for temporary buffers
+// Temporary buffers are used in two case:
+// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5
+// 2. "in-place" FFTs.  Notice the quotes, since kissfft does not really do an in-place transform.
+#include <alloca.h>
+#define  KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes)
+#define  KISS_FFT_TMP_FREE(ptr)
+#else
+#define  KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes)
+#define  KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr)
+#endif

src/codec2_fdmdv.h

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: codec2_fdmdv.h
+  AUTHOR......: David Rowe
+  DATE CREATED: April 14 2012
+
+  A 1400 bit/s (nominal) Frequency Division Multiplexed Digital Voice
+  (FDMDV) modem.  Used for digital audio over HF SSB. See
+  README_fdmdv.txt for more information, and fdmdv_mod.c and
+  fdmdv_demod.c for example usage.
+
+  The name codec2_fdmdv.h is used to make it unique when "make
+  installed".
+
+  References:
+
+    [1] http://n1su.com/fdmdv/FDMDV_Docs_Rel_1_4b.pdf
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2012 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __FDMDV__
+#define __FDMDV__
+
+#include "comp.h"
+#include "modem_stats.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* set up the calling convention for DLL function import/export for
+   WIN32 cross compiling */
+
+#ifdef __CODEC2_WIN32__
+#ifdef __CODEC2_BUILDING_DLL__
+#define CODEC2_WIN32SUPPORT __declspec(dllexport) __stdcall
+#else
+#define CODEC2_WIN32SUPPORT __declspec(dllimport) __stdcall
+#endif
+#else
+#define CODEC2_WIN32SUPPORT
+#endif
+
+#define FDMDV_NC                      14  /* default number of data carriers                                */
+#define FDMDV_NC_MAX                  20  /* maximum number of data carriers                                */
+#define FDMDV_BITS_PER_FRAME          28  /* 20ms frames, for nominal 1400 bit/s                            */
+#define FDMDV_NOM_SAMPLES_PER_FRAME  160  /* modulator output samples/frame and nominal demod samples/frame */
+                                          /* at 8000 Hz sample rate                                         */
+#define FDMDV_MAX_SAMPLES_PER_FRAME  200  /* max demod samples/frame, use this to allocate storage          */
+#define FDMDV_SCALE                 1000  /* suggested scaling for 16 bit shorts                            */
+#define FDMDV_FCENTRE               1500  /* Centre frequency, Nc/2 carriers below this, Nc/2 carriers above (Hz) */
+
+/* 8 to 48 kHz sample rate conversion */
+
+#define FDMDV_OS                 2                            /* oversampling rate                   */
+#define FDMDV_OS_TAPS_16K       48                            /* number of OS filter taps at 16kHz   */
+#define FDMDV_OS_TAPS_8K        (FDMDV_OS_TAPS_16K/FDMDV_OS)  /* number of OS filter taps at 8kHz    */
+
+/* FDMDV states and stats structures */
+
+struct FDMDV;
+
+struct FDMDV * fdmdv_create(int Nc);
+void           fdmdv_destroy(struct FDMDV *fdmdv_state);
+void           fdmdv_use_old_qpsk_mapping(struct FDMDV *fdmdv_state);
+int            fdmdv_bits_per_frame(struct FDMDV *fdmdv_state);
+float          fdmdv_get_fsep(struct FDMDV *fdmdv_state);
+void           fdmdv_set_fsep(struct FDMDV *fdmdv_state, float fsep);
+
+void           fdmdv_mod(struct FDMDV *fdmdv_state, COMP tx_fdm[], int tx_bits[], int *sync_bit);
+void           fdmdv_demod(struct FDMDV *fdmdv_state, int rx_bits[], int *reliable_sync_bit, COMP rx_fdm[], int *nin);
+
+void           fdmdv_get_test_bits(struct FDMDV *fdmdv_state, int tx_bits[]);
+int            fdmdv_error_pattern_size(struct FDMDV *fdmdv_state);
+void           fdmdv_put_test_bits(struct FDMDV *f, int *sync, short error_pattern[], int *bit_errors, int *ntest_bits, int rx_bits[]);
+
+void           fdmdv_get_demod_stats(struct FDMDV *fdmdv_state, struct MODEM_STATS *stats);
+
+void           fdmdv_8_to_16(float out16k[], float in8k[], int n);
+void           fdmdv_8_to_16_short(short out16k[], short in8k[], int n);
+void           fdmdv_16_to_8(float out8k[], float in16k[], int n);
+void           fdmdv_16_to_8_short(short out8k[], short in16k[], int n);
+
+void           fdmdv_freq_shift(COMP rx_fdm_fcorr[], COMP rx_fdm[], float foff, COMP *foff_phase_rect, int nin);
+
+/* debug/development function(s) */
+
+void fdmdv_dump_osc_mags(struct FDMDV *f);
+void fdmdv_simulate_channel(float *sig_pwr_av, COMP samples[], int nin, float target_snr);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
+

src/comp.h

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: comp.h
+  AUTHOR......: David Rowe
+  DATE CREATED: 24/08/09
+
+  Complex number definition.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __COMP__
+#define __COMP__
+
+/* Complex number */
+
+typedef struct {
+  float real;
+  float imag;
+} COMP;
+
+#endif

src/comp_prim.h

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: comp_prim.h
+  AUTHOR......: David Rowe
+  DATE CREATED: Marh 2015
+
+  Complex number maths primitives.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2015 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __COMP_PRIM__
+#define __COMP_PRIM__
+
+/*---------------------------------------------------------------------------*\
+
+                               FUNCTIONS
+
+\*---------------------------------------------------------------------------*/
+
+inline static COMP cneg(COMP a)
+{
+    COMP res;
+
+    res.real = -a.real;
+    res.imag = -a.imag;
+
+    return res;
+}
+
+inline static COMP cconj(COMP a)
+{
+    COMP res;
+
+    res.real = a.real;
+    res.imag = -a.imag;
+
+    return res;
+}
+
+inline static COMP cmult(COMP a, COMP b)
+{
+    COMP res;
+
+    res.real = a.real*b.real - a.imag*b.imag;
+    res.imag = a.real*b.imag + a.imag*b.real;
+
+    return res;
+}
+
+inline static COMP fcmult(float a, COMP b)
+{
+    COMP res;
+
+    res.real = a*b.real;
+    res.imag = a*b.imag;
+
+    return res;
+}
+
+inline static COMP cadd(COMP a, COMP b)
+{
+    COMP res;
+
+    res.real = a.real + b.real;
+    res.imag = a.imag + b.imag;
+
+    return res;
+}
+
+inline static float cabsolute(COMP a)
+{
+    return sqrtf((a.real * a.real) + (a.imag * a.imag) );
+}
+
+/*
+ * Euler's formula in a new convenient function
+ */
+inline static COMP comp_exp_j(float phi){
+    COMP res;
+    res.real = cosf(phi);
+    res.imag = sinf(phi);
+    return res;
+}
+
+/*
+ * Quick and easy complex 0
+ */
+inline static COMP comp0(){
+    COMP res;
+    res.real = 0;
+    res.imag = 0;
+    return res;
+}
+
+/*
+ * Quick and easy complex subtract
+ */
+inline static COMP csub(COMP a, COMP b){
+    COMP res;
+    res.real = a.real-b.real;
+    res.imag = a.imag-b.imag;
+    return res;
+}
+
+/*
+ * Compare the magnitude of a and b. if |a|>|b|, return true, otw false.
+ * This needs no square roots
+ */
+inline static int comp_mag_gt(COMP a,COMP b){
+    return ((a.real*a.real)+(a.imag*a.imag)) > ((b.real*b.real)+(b.imag*b.imag));
+}
+
+/*
+ * Normalize a complex number's magnitude to 1
+ */
+inline static COMP comp_normalize(COMP a){
+    COMP b;
+    float av = cabsolute(a);
+    b.real = a.real/av;
+    b.imag = a.imag/av;
+    return b;
+}
+
+#endif

src/debug_alloc.h

@@ -0,0 +1,61 @@
+/* debug_alloc.h
+ *
+ * Some macros which can report on malloc results.
+ *
+ * Enable with "-D DEBUG_ALLOC"
+ */
+
+#ifndef DEBUG_ALLOC_H
+#define DEBUG_ALLOC_H
+
+#include <stdio.h>
+
+// Debug calls
+
+#ifdef CORTEX_M4
+extern char * __heap_end;
+register char * sp asm ("sp");
+#endif
+
+  static inline void * DEBUG_MALLOC(const char *func, size_t size) {
+    void *ptr = malloc(size);
+    fprintf(stderr, "MALLOC: %s %p %d", func, ptr, (int)size);
+#ifdef CORTEX_M4
+
+    fprintf(stderr, " : sp %p ", sp);
+#endif
+    if (!ptr) fprintf(stderr, " ** FAILED **");
+    fprintf(stderr, "\n");
+    return(ptr);
+    }
+
+  static inline void * DEBUG_CALLOC(const char *func, size_t nmemb, size_t size) {
+    void *ptr = calloc(nmemb, size);
+    fprintf(stderr, "CALLOC: %s %p %d %d", func, ptr, (int)nmemb, (int)size);
+#ifdef CORTEX_M4
+    fprintf(stderr, " : sp %p ", sp);
+#endif
+    if (!ptr) fprintf(stderr, " ** FAILED **");
+    fprintf(stderr, "\n");
+    return(ptr);
+    }
+  static inline void DEBUG_FREE(const char *func, void *ptr) {
+    free(ptr);
+    fprintf(stderr, "FREE: %s %p\n", func, ptr);
+    }
+
+#ifdef DEBUG_ALLOC
+  #define MALLOC(size) DEBUG_MALLOC(__func__, size)
+  #define CALLOC(nmemb, size) DEBUG_CALLOC(__func__, nmemb, size)
+  #define FREE(ptr) DEBUG_FREE(__func__, ptr)
+#else //DEBUG_ALLOC
+// Default to normal calls
+  #define MALLOC(size) malloc(size)
+
+  #define CALLOC(nmemb, size) calloc(nmemb, size)
+
+  #define FREE(ptr) free(ptr)
+
+#endif //DEBUG_ALLOC
+
+#endif //DEBUG_ALLOC_H

src/drs232.c

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: drs232.c
+  AUTHOR......: David Rowe
+  DATE CREATED: March 2016
+
+  Looks for a unique word in series of bits.  When found, deframes a RS232
+  encoded frame of bytes.  Used for high bit rate Horus SSTV reception.
+
+  Frame format:
+
+    16 bytes 0x55 - 0xabcdef01 UW - 256 bytes of payload - 2 bytes CRC
+
+  Each byte is encoded as a 10 bit RS232 serial word: 
+  
+    0 LSB .... MSB 1
+
+  Building:
+   
+    $ gcc drs232.c -o drs232 -Wall
+
+\*---------------------------------------------------------------------------*/
+
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+
+/* states -----------------------------------------------*/
+
+#define LOOK_FOR_UW    0
+#define COLLECT_PACKET 1
+
+/* packet parameters */
+
+#define UW_BYTES               4
+#define UW_BITS                40
+#define UW_ALLOWED_ERRORS      5
+#define BYTES_PER_PACKET       256
+#define CRC_BYTES              2
+#define BITS_PER_BYTE          10
+#define UNPACKED_PACKET_BYTES  ((UW_BYTES+BYTES_PER_PACKET+CRC_BYTES)*BITS_PER_BYTE)
+
+/* UW pattern we look for, including start/stop bits */
+
+uint8_t uw[] = {
+    /* 0xb                0xa */
+    0, 1, 1, 0, 1, 0, 1, 0, 1, 1,
+    /* 0xd                0xc */
+    0, 1, 0, 1, 1, 0, 0, 1, 1, 1,
+    /* 0xf                0xe */
+    0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
+    /* 0x1                0x0 */
+    0, 1, 0, 0, 0, 0, 0, 0, 0, 1,
+};
+
+ // from http://stackoverflow.com/questions/10564491/function-to-calculate-a-crc16-checksum
+
+unsigned short gen_crc16(unsigned char* data_p, int length){
+    unsigned char x;
+    unsigned short crc = 0xFFFF;
+
+    while (length--){
+        x = crc >> 8 ^ *data_p++;
+        x ^= x>>4;
+        crc = (crc << 8) ^ ((unsigned short)(x << 12)) ^ ((unsigned short)(x <<5)) ^ ((unsigned short)x);
+    }
+    
+    return crc;
+}
+
+
+int main(int argc, char *argv[]) {
+    FILE       *fin, *fout;
+    int         state, next_state, i, j, k, ind, score, bits_read;
+    char        bit;
+    uint8_t     unpacked_packet[UNPACKED_PACKET_BYTES];
+    uint8_t     packet[BYTES_PER_PACKET+CRC_BYTES];
+    uint8_t     abyte;
+    uint16_t    tx_checksum, rx_checksum;
+    int         verbose, packet_errors, packets;
+
+    if (argc < 3) {
+	fprintf(stderr, "usage: drs232 InputOneBitPerChar OutputPackets [-v[v]]\n");
+ 	exit(1);
+    }
+
+    if (strcmp(argv[1], "-")  == 0) fin = stdin;
+    else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
+	fprintf(stderr, "Error opening input file: %s: %s.\n",
+         argv[1], strerror(errno));
+	exit(1);
+    }
+
+    if (strcmp(argv[2], "-") == 0) fout = stdout;
+    else if ( (fout = fopen(argv[2],"wb")) == NULL ) {
+	fprintf(stderr, "Error opening output file: %s: %s.\n",
+         argv[2], strerror(errno));
+	exit(1);
+    }
+
+    verbose = 0;
+    if (argc > 3) {
+        if (strcmp(argv[3], "-v") == 0) {
+            verbose = 1;
+        }
+        if (strcmp(argv[3], "-vv") == 0) {
+            verbose = 2;
+        }
+    }
+
+    state = LOOK_FOR_UW;
+    for(i=0; i<UNPACKED_PACKET_BYTES; i++)
+        unpacked_packet[i] = 0;
+    bits_read = 0; 
+
+    packet_errors = packets = 0;
+
+    while(fread(&bit, sizeof(char), 1, fin) == 1) {
+        bits_read++;
+
+        next_state = state;
+
+        if (state == LOOK_FOR_UW) {
+
+            /* put latest input bit into sliding buffer */
+            
+            for(i=0; i<UW_BITS-1; i++) {
+                unpacked_packet[i] = unpacked_packet[i+1];
+            }
+            unpacked_packet[i] = bit;
+
+            /* lets see if it matches the UW */
+
+            score = 0;
+            for(i=0; i<UW_BITS; i++) {
+                score += (unpacked_packet[i] == uw[i]);
+                /* if (i == BITS_PER_BYTE)
+                    printf(" ");
+                    printf("%1d", unpacked_packet[i]); */
+            }
+            //printf("\n");
+            
+            //fprintf(stderr,"UW score: %d\n", score);
+            if (score > (UW_BITS-UW_ALLOWED_ERRORS)) {
+                if (verbose == 2) {
+                    fprintf(stderr,"UW found!\n");
+                }
+                ind = UW_BITS;
+                next_state = COLLECT_PACKET;
+            }             
+        }
+
+        if (state == COLLECT_PACKET) {
+            unpacked_packet[ind++] = bit;
+ 
+           if (ind == UNPACKED_PACKET_BYTES) {
+
+               /* OK we have enough bits, remove RS232 sync bits and pack */
+
+               for(i=UW_BITS,k=0; i<UNPACKED_PACKET_BYTES; i+=BITS_PER_BYTE,k++) {
+                   //for(j=0; j<BITS_PER_BYTE; j++)
+                   //    printf("%1d", unpacked_packet[i+j]);
+                   //printf("\n");
+                   abyte = 0;
+                   for(j=0; j<8; j++) {
+                       abyte |= unpacked_packet[i+j+1] << j;
+                   }
+                   packet[k] = abyte;
+                   //printf("k:%d 0x%02x\n", k, packet[k]);
+                   //if (k == 4)
+                   //  exit(0);
+               }
+
+               /* output if CRC check is OK */
+
+               rx_checksum = gen_crc16(packet, BYTES_PER_PACKET);
+               tx_checksum = packet[BYTES_PER_PACKET] + (packet[BYTES_PER_PACKET+1] << 8);
+
+               if (verbose == 2) {
+                   fprintf(stderr, "k=%d bytes after UW  tx_checksum: 0x%02x rx_checksum: 0x%02x\n", 
+                           k, tx_checksum, rx_checksum);
+               }
+
+               packets++;
+               if (rx_checksum == tx_checksum) {
+                   fwrite(packet, sizeof(char), BYTES_PER_PACKET, fout);
+               }
+               else
+                   packet_errors++;
+
+               if (verbose) {
+                   fprintf(stderr, "packets: %d packet_errors: %d PER: %4.3f\n", 
+                           packets, packet_errors, 
+                           (float)packet_errors/packets);
+               }
+
+               next_state = LOOK_FOR_UW;
+            }
+
+        }
+        //if (bits_read == (16*10 + UNPACKED_PACKET_BYTES))
+        //    exit(0);
+
+        state = next_state;      
+    }
+
+    fclose(fin);
+    fclose(fout);
+
+    fprintf(stderr, "packets: %d packet_errors: %d PER: %4.3f\n", packets, packet_errors, 
+            (float)packet_errors/packets);
+
+    return 0;
+}

src/drs232_ldpc.c

@@ -0,0 +1,286 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: drs232_ldpc.c
+  AUTHOR......: David Rowe
+  DATE CREATED: Sep 2016
+
+  Looks for a unique word in series of soft decision symbols.  When
+  found, deframes a RS232 encoded frame of soft decision bit, LDPC
+  decodes, and outputs a frame of packed bytes.  Used for high bit
+  rate Horus SSTV reception.
+
+  Frame format:
+
+    16 bytes 0x55 - 0xabcdef01 UW - 256 bytes of payload - 2 bytes CRC - 65 bytes LPDC Parity
+
+  Each byte is encoded as a 10 bit RS232 serial word: 
+  
+    0 LSB .... MSB 1
+
+  Building:
+   
+    $ gcc drs232_ldpc.c mpdecode_core.c -o drs232_ldpc -Wall -lm
+
+\*---------------------------------------------------------------------------*/
+
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+
+#include "mpdecode_core.h"
+
+/* Machine generated consts, H_rows, H_cols, test input/output data to
+   change LDPC code regenerate this file. */
+
+#include "H2064_516_sparse.h"  
+
+/* states -----------------------------------------------*/
+
+#define LOOK_FOR_UW    0
+#define COLLECT_PACKET 1
+
+/* packet parameters */
+
+#define UW_BYTES               4
+#define UW_BITS                40
+#define UW_ALLOWED_ERRORS      5
+#define BYTES_PER_PACKET       256
+#define CRC_BYTES              2
+#define PARITY_BYTES           65
+#define BITS_PER_BYTE          10
+#define UNPACKED_PACKET_BYTES  ((UW_BYTES+BYTES_PER_PACKET+CRC_BYTES)*BITS_PER_BYTE)
+#define SYMBOLS_PER_PACKET     (BYTES_PER_PACKET+CRC_BYTES+PARITY_BYTES)*BITS_PER_BYTE
+
+/* UW pattern we look for, including start/stop bits */
+
+uint8_t uw[] = {
+    /* 0xb                0xa */
+    0, 1, 1, 0, 1, 0, 1, 0, 1, 1,
+    /* 0xd                0xc */
+    0, 1, 0, 1, 1, 0, 0, 1, 1, 1,
+    /* 0xf                0xe */
+    0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
+    /* 0x1                0x0 */
+    0, 1, 0, 0, 0, 0, 0, 0, 0, 1,
+};
+
+
+// from http://stackoverflow.com/questions/10564491/function-to-calculate-a-crc16-checksum
+
+unsigned short gen_crc16(unsigned char* data_p, int length){
+    unsigned char x;
+    unsigned short crc = 0xFFFF;
+
+    while (length--){
+        x = crc >> 8 ^ *data_p++;
+        x ^= x>>4;
+        crc = (crc << 8) ^ ((unsigned short)(x << 12)) ^ ((unsigned short)(x <<5)) ^ ((unsigned short)x);
+    }
+    
+    return crc;
+}
+
+
+int main(int argc, char *argv[]) {
+    FILE       *fin, *fout;
+    int         state, next_state, i, j, k, ind, score, verbose;
+    float       symbol;
+    uint8_t     bit, bit_buffer[UW_BITS];
+    double      symbol_buf[SYMBOLS_PER_PACKET];
+    double      symbol_buf_no_rs232[SYMBOLS_PER_PACKET];
+    float       llr[SYMBOLS_PER_PACKET];
+    uint8_t     unpacked_packet[CODELENGTH];
+    uint8_t     packet[BYTES_PER_PACKET+CRC_BYTES];
+    uint8_t     abyte;
+    uint16_t    tx_checksum, rx_checksum, packet_errors, packets;
+    int         CodeLength, iter, parityCheckCount;
+    struct LDPC ldpc;
+
+    assert(sizeof(uw) == UW_BITS);
+
+    /* LDPC parameters */
+
+    CodeLength = CODELENGTH;                    /* length of entire codeword in bits */
+
+    /* set up LDPC code from include file constants */
+
+    ldpc.max_iter = MAX_ITER;
+    ldpc.dec_type = 0;
+    ldpc.q_scale_factor = 1;
+    ldpc.r_scale_factor = 1;
+    ldpc.CodeLength = CODELENGTH;
+    ldpc.NumberParityBits = NUMBERPARITYBITS;
+    ldpc.NumberRowsHcols = NUMBERROWSHCOLS;
+    ldpc.max_row_weight = MAX_ROW_WEIGHT;
+    ldpc.max_col_weight = MAX_COL_WEIGHT;
+    ldpc.H_rows = H_rows;
+    ldpc.H_cols = H_cols;
+
+    /* process command line ----------------------------------------------*/
+
+    if (argc < 3) {
+	fprintf(stderr, "usage: drs232 InputOneSymbolPerFloat OutputPackets [-v[v]]\n");
+ 	exit(1);
+    }
+
+    if (strcmp(argv[1], "-")  == 0) fin = stdin;
+    else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
+	fprintf(stderr, "Error opening input file: %s: %s.\n",
+         argv[1], strerror(errno));
+	exit(1);
+    }
+
+    if (strcmp(argv[2], "-") == 0) fout = stdout;
+    else if ( (fout = fopen(argv[2],"wb")) == NULL ) {
+	fprintf(stderr, "Error opening output file: %s: %s.\n",
+         argv[2], strerror(errno));
+	exit(1);
+    }
+
+    verbose = 0;
+    if (argc > 3) {
+        if (strcmp(argv[3], "-v") == 0) {
+            verbose = 1;
+        }
+        if (strcmp(argv[3], "-vv") == 0) {
+            verbose = 2;
+        }
+    }
+
+    state = LOOK_FOR_UW;
+    memset(bit_buffer,0,  sizeof(bit_buffer));
+
+    packet_errors = packets = 0;
+
+    while(fread(&symbol, sizeof(float), 1, fin) == 1) {
+
+        /* make hard decision for purpose of UW detection */
+
+        bit = symbol < 0;
+        //printf("symbol; %f bit: %d\n", symbol, bit);
+        next_state = state;
+        if (state == LOOK_FOR_UW) {
+
+            /* put latest input bit into sliding buffer */
+            
+            for(i=0; i<UW_BITS-1; i++) {
+                bit_buffer[i] = bit_buffer[i+1];
+            }
+            bit_buffer[i] = bit;
+
+            /* lets see if it matches the UW */
+
+            score = 0;
+            for(i=0; i<UW_BITS; i++) {
+                score += (bit_buffer[i] == uw[i]);
+                /* if (i == BITS_PER_BYTE)
+                    printf(" ");
+                    printf("%1d", unpacked_packet[i]); */
+            }
+            //printf("\n");
+            
+            //fprintf(stderr,"UW score: %d\n", score);
+            if (score >= (UW_BITS-UW_ALLOWED_ERRORS)) {
+                //fprintf(stderr,"UW found! score: %d\n verbose: %d\n", score, verbose);
+                ind = 0;
+                next_state = COLLECT_PACKET;
+            }             
+        }
+
+        if (state == COLLECT_PACKET) {
+            symbol_buf[ind++] = symbol;
+ 
+            if (ind == SYMBOLS_PER_PACKET) {
+
+               /* OK we have enough bits, remove RS232 sync symbols.
+                  This is set up for bit<->byte ordering as per python
+                  tx code */
+
+               for(i=0,k=0; i<SYMBOLS_PER_PACKET; i+=BITS_PER_BYTE) {
+                   for(j=0; j<8; j++) {
+                       symbol_buf_no_rs232[k+j] = symbol_buf[i+7-j+1];
+                   }
+                   k += 8;
+               }
+
+               /* now LDPC decode */
+
+               sd_to_llr(llr, symbol_buf_no_rs232, CodeLength);
+               iter = run_ldpc_decoder(&ldpc, unpacked_packet, llr, &parityCheckCount);
+
+               /* pack into bytes */
+
+               for(i=0; i<BYTES_PER_PACKET+CRC_BYTES; i++) {
+                   abyte = 0;
+                   for(j=0; j<8; j++)
+                       abyte |= unpacked_packet[8*i+j] << (7-j);
+                   packet[i] = abyte;
+               }
+
+               /* then output if CRC check is OK */
+
+               rx_checksum = gen_crc16(packet, BYTES_PER_PACKET);
+               tx_checksum = packet[BYTES_PER_PACKET] + (packet[BYTES_PER_PACKET+1] << 8);
+
+               if (verbose == 2) {
+                   if (rx_checksum != tx_checksum) {
+                       fprintf(stderr, "tx_checksum: 0x%02x rx_checksum: 0x%02x\n", 
+                               tx_checksum, rx_checksum);
+                   }
+               }
+               
+               packets++;
+               if (rx_checksum == tx_checksum) {
+                   fwrite(packet, sizeof(char), BYTES_PER_PACKET, fout);
+                   fflush(fout);
+               }
+               else
+                   packet_errors++;
+
+               if (verbose) {
+                   fprintf(stderr, "packets: %d packet_errors: %d PER: %4.3f iter: %d\n", 
+                           packets, packet_errors, 
+                           (float)packet_errors/packets, iter);
+               }
+               //exit(0);
+               next_state = LOOK_FOR_UW;
+            }
+
+        }
+        //if (bits_read == (16*10 + UNPACKED_PACKET_BYTES))
+        //    exit(0);
+
+        state = next_state;      
+    }
+
+    fclose(fin);
+    fclose(fout);
+
+    fprintf(stderr, "packets: %d packet_errors: %d PER: %4.3f\n", packets, packet_errors, 
+            (float)packet_errors/packets);
+
+    return 0;
+}
+
+

src/fsk.h

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: fsk.h
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 6 January 2016
+
+  C Implementation of 2FSK/4FSK modulator/demodulator, based on octave/fsk_horus.m
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef __C2FSK_H
+#define __C2FSK_H
+
+#include <stdint.h>
+#include "comp.h"
+#include "kiss_fftr.h"
+#include "modem_stats.h"
+
+#define MODE_2FSK 2
+#define MODE_4FSK 4
+
+#define MODE_M_MAX 4
+
+#define FSK_SCALE 16383
+
+/* default internal parameters */
+#define FSK_DEFAULT_P 8
+#define FSK_DEFAULT_NSYM 50
+
+struct FSK {
+    /*  Static parameters set up by fsk_init */
+    int Ndft;               /* freq offset est fft */
+    int Fs;                 /* sample freq */
+    int N;                  /* processing buffer size */
+    int Rs;                 /* symbol rate */
+    int Ts;                 /* samples per symbol */
+    int Nmem;               /* size of extra mem for timing adj */
+    int P;                  /* oversample rate for timing est/adj */
+    int Nsym;               /* Number of symbols spat out in a processing frame */
+    int Nbits;              /* Number of bits spat out in a processing frame */
+    int f1_tx;              /* f1 for modulator */
+    int fs_tx;              /* Space between TX freqs for modulatosr */
+    int mode;               /* 2FSK or 4FSK */
+    float tc;               /* time constant for smoothing FFTs */
+    int est_min;            /* Minimum frequency for freq. estimator */
+    int est_max;            /* Maximum frequency for freq. estimaotr */
+    int est_space;          /* Minimum frequency spacing for freq. estimator */
+    float* hann_table;	    /* Precomputed or runtime computed hann window table */
+    
+    /*  Parameters used by demod */
+    float* Sf;	            /* Average of magnitude spectrum */
+    COMP phi_c[MODE_M_MAX]; /* phase of each demod local oscillator */
+    COMP *f_dc;             /* down converted samples               */
+    
+    kiss_fft_cfg fft_cfg;   /* Config for KISS FFT, used in freq est */
+    float norm_rx_timing;   /* Normalized RX timing */
+        
+    
+    /*  Parameters used by mod */
+    COMP tx_phase_c;        /* TX phase, but complex */ 
+    
+    /*  Statistics generated by demod */
+    float EbNodB;            /* Estimated EbNo in dB */
+    float f_est[MODE_M_MAX]; /* Estimated frequencies (peak method) */
+    float f2_est[MODE_M_MAX];/* Estimated frequencies (mask method) */
+    int   freq_est_type;     /* which estimator to use              */
+    float ppm;               /* Estimated PPM clock offset */
+    
+    /*  Parameters used by mod/demod and driving code */
+    int nin;                /* Number of samples to feed the next demod cycle */
+    int burst_mode;         /* enables/disables 'burst' mode */
+    int lock_nin;           /* locks nin during testing */
+    
+    /*  modem statistic struct */
+    struct MODEM_STATS *stats;
+    int normalise_eye;      /* enables/disables normalisation of eye diagram */
+};
+
+/*
+ * Create an FSK config/state struct from a set of config parameters
+ * 
+ * int Fs - Sample frequency
+ * int Rs - Symbol rate
+ * int tx_f1 - '0' frequency
+ * int tx_fs - frequency spacing
+ */
+struct FSK * fsk_create(int Fs, int Rs, int M, int tx_f1, int tx_fs);
+
+/*
+ * Create an FSK config/state struct from a set of config parameters
+ * 
+ * int Fs - Sample frequency
+ * int Rs - Symbol rate
+ * int tx_f1 - '0' frequency
+ * int tx_fs - frequency spacing
+ */
+struct FSK * fsk_create_hbr(int Fs, int Rs, int M, int P, int Nsym, int tx_f1, int tx_fs);
+
+/*
+ * Set the minimum and maximum frequencies at which the freq. estimator can find tones
+ */
+void fsk_set_freq_est_limits(struct FSK *fsk,int fmin, int fmax);
+
+/* 
+ * Clear the estimator states
+ */
+void fsk_clear_estimators(struct FSK *fsk);
+
+/*
+ * Fills MODEM_STATS struct with demod statistics
+ */
+void fsk_get_demod_stats(struct FSK *fsk, struct MODEM_STATS *stats);
+
+/*
+ * Destroy an FSK state struct and free it's memory
+ * 
+ * struct FSK *fsk - FSK config/state struct to be destroyed
+ */
+void fsk_destroy(struct FSK *fsk);
+
+/*
+ * Modulates Nsym bits into N samples
+ * 
+ * struct FSK *fsk - FSK config/state struct, set up by fsk_create
+ * float fsk_out[] - Buffer for N samples of modulated FSK
+ * uint8_t tx_bits[] - Buffer containing Nbits unpacked bits
+ */
+void fsk_mod(struct FSK *fsk, float fsk_out[], uint8_t tx_bits[]);
+
+/*
+ * Modulates Nsym bits into N samples
+ * 
+ * struct FSK *fsk - FSK config/state struct, set up by fsk_create
+ * float fsk_out[] - Buffer for N samples of "voltage" used to modulate an external VCO
+ * uint8_t tx_bits[] - Buffer containing Nbits unpacked bits
+ */
+void fsk_mod_ext_vco(struct FSK *fsk, float vco_out[], uint8_t tx_bits[]);
+
+/*
+ * Modulates Nsym bits into N complex samples
+ * 
+ * struct FSK *fsk - FSK config/state struct, set up by fsk_create
+ * comp fsk_out[] - Buffer for N samples of modulated FSK
+ * uint8_t tx_bits[] - Buffer containing Nbits unpacked bits
+ */
+void fsk_mod_c(struct FSK *fsk, COMP fsk_out[], uint8_t tx_bits[]);
+
+
+/*
+ * Returns the number of samples needed for the next fsk_demod() cycle
+ *
+ * struct FSK *fsk - FSK config/state struct, set up by fsk_create
+ * returns - number of samples to be fed into fsk_demod next cycle 
+ */
+uint32_t fsk_nin(struct FSK *fsk);
+
+
+/*
+ * Demodulate some number of FSK samples. The number of samples to be 
+ *  demodulated can be found by calling fsk_nin().
+ * 
+ * struct FSK *fsk - FSK config/state struct, set up by fsk_create
+ * uint8_t rx_bits[] - Buffer for Nbits unpacked bits to be written
+ * float fsk_in[] - nin samples of modulated FSK
+ */
+void fsk_demod(struct FSK *fsk, uint8_t rx_bits[],COMP fsk_in[]);
+
+/*
+ * Demodulate some number of FSK samples. The number of samples to be 
+ *  demodulated can be found by calling fsk_nin().
+ * 
+ * struct FSK *fsk - FSK config/state struct, set up by fsk_create
+ * float rx_bits[] - Buffer for Nbits soft decision bits to be written
+ * float fsk_in[] - nin samples of modualted FSK
+ */
+void fsk_demod_sd(struct FSK *fsk, float rx_bits[],COMP fsk_in[]);
+
+/* enables/disables normalisation of eye diagram samples */
+  
+void fsk_stats_normalise_eye(struct FSK *fsk, int normalise_enable);
+
+/* Set the FSK modem into burst demod mode */
+
+void fsk_enable_burst_mode(struct FSK *fsk);
+
+/* Set freq est algorithm 0: peak 1:mask */
+void fsk_set_freq_est_alg(struct FSK *fsk, int est_type);
+
+#endif

src/fsk_demod.c

@@ -0,0 +1,445 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: fsk_demod.c
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 8 January 2016
+
+  C test driver for fsk_demod in fsk.c. Reads in a stream of 32 bit cpu endian
+  floats and writes out the detected bits
+   
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#define TEST_FRAME_SIZE 100  /* must match fsk_get_test_bits.c */
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <getopt.h>
+#include <time.h>
+#include <signal.h>
+#include <unistd.h>
+
+#include "fsk.h"
+#include "codec2_fdmdv.h"
+#include "modem_stats.h"
+
+/* cleanly exit when we get a SIGTERM */
+
+void sig_handler(int signo)
+{
+    if (signo == SIGTERM) {
+        exit(0);
+    }
+}
+
+int main(int argc,char *argv[]){
+    struct FSK *fsk;
+    struct MODEM_STATS stats;
+    int Fs,Rs,M,P,stats_ctr,stats_loop;
+    float loop_time;
+    int enable_stats = 0;
+    FILE *fin,*fout;
+    uint8_t *bitbuf = NULL;
+    int16_t *rawbuf;
+    COMP *modbuf;
+    float *sdbuf = NULL;
+    int i,j,Ndft;
+    int soft_dec_mode = 0;
+    stats_loop = 0;
+    int complex_input = 1, bytes_per_sample = 2;
+    int stats_rate = 8;
+    int testframe_mode = 0;
+    P = 8;  /* default */
+    M = 0;
+    int fsk_lower = 0;
+    int fsk_upper = 0;
+    int user_fsk_lower = 0;
+    int user_fsk_upper = 0;
+    int nsym = FSK_DEFAULT_NSYM;
+    int mask = 0;
+    int tx_tone_separation = 100;
+    
+    int o = 0;
+    int opt_idx = 0;
+    while( o != -1 ){
+        static struct option long_opts[] = {
+            {"help",      no_argument,        0, 'h'},
+            {"conv",      required_argument,  0, 'p'},
+            {"cs16",      no_argument,        0, 'c'},
+            {"cu8",       no_argument,        0, 'd'},
+            {"fsk_lower", required_argument,  0, 'b'},
+            {"fsk_upper", required_argument,  0, 'u'},
+            {"stats",     optional_argument,  0, 't'},
+            {"soft-dec",  no_argument,        0, 's'},
+            {"testframes",no_argument,        0, 'f'},
+            {"nsym",      required_argument,  0, 'n'},
+            {"mask",      required_argument,  0, 'm'},
+            {0, 0, 0, 0}
+        };
+        
+        o = getopt_long(argc,argv,"fhlp:cdt::sb:u:m",long_opts,&opt_idx);
+        
+        switch(o){
+        case 'c':
+            complex_input = 2;
+            bytes_per_sample = 2;
+            break;
+        case 'd':
+            complex_input = 2;
+            bytes_per_sample = 1;
+            break;
+        case 'f':
+            testframe_mode = 1;
+            break;
+        case 't':
+            enable_stats = 1;
+            if(optarg != NULL){
+                stats_rate = atoi(optarg);
+                if(stats_rate == 0){
+                    stats_rate = 8;
+                }
+            }
+            break;
+        case 's':
+            soft_dec_mode = 1;
+            break;
+        case 'p':
+            P = atoi(optarg);
+            break;
+        case 'b':
+            if (optarg != NULL) {
+                fsk_lower = atoi(optarg);
+                user_fsk_lower = 1;
+            }
+            break;
+        case 'u':
+            if (optarg != NULL){
+                fsk_upper = atoi(optarg);
+                 user_fsk_upper = 1;
+            }
+            break;
+        case 'n':
+            if (optarg != NULL) {
+                nsym = atoi(optarg);
+            }
+            break;
+        case 'm':
+            mask = 1;
+            tx_tone_separation = atoi(optarg);
+            break;
+        case 'h':
+        case '?':
+            goto helpmsg;
+            break;
+        }
+    }
+    int dx = optind;
+    
+    if( (argc - dx) < 5){
+        fprintf(stderr, "Too few arguments\n");
+        goto helpmsg;
+    }
+    
+    if( (argc - dx) > 5) {
+        fprintf(stderr, "Too many arguments\n");
+    helpmsg:
+        fprintf(stderr,"usage: %s [options] (2|4) SampleRate SymbolRate InputModemRawFile OutputFile\n",argv[0]);
+        fprintf(stderr," -c --cs16          The raw input file will be in complex signed 16 bit format.\n");
+        fprintf(stderr," -d --cu8           The raw input file will be in complex unsigned 8 bit format.\n");
+        fprintf(stderr,"                    If neither -c nor -d are used, the input should be in signed 16 bit format.\n");
+        fprintf(stderr," -f --testframes    Testframe mode, prints stats to stderr when a testframe is detected, if -t (JSON) \n");
+        fprintf(stderr,"                    is enabled stats will be in JSON format\n");
+        fprintf(stderr," -t[r] --stats=[r]  Print out modem statistics to stderr in JSON.\n");
+        fprintf(stderr,"                    r, if provided, sets the number of modem frames between statistic printouts.\n");
+        fprintf(stderr," -s --soft-dec      The output file will be in a soft-decision format, with one 32-bit float per bit.\n");
+        fprintf(stderr,"                    If -s is not used, the output will be in a 1 byte-per-bit format.\n");
+        fprintf(stderr," -p P               The demod internals operate at a rate of Fs/P, default %d\n", FSK_DEFAULT_P);
+        fprintf(stderr,"                    P must be divisible by the symbol rate. Smaller P values will result in faster\n");
+        fprintf(stderr,"                    processing but lower demodulation performance. Default %d\n", FSK_DEFAULT_P);
+        fprintf(stderr," --fsk_lower freq   lower limit of freq estimator (default 0 for real input, -Fs/2  for complex input)\n");
+        fprintf(stderr," --fsk_upper freq   upper limit of freq estimator (default Fs/2)\n");
+        fprintf(stderr," --nsym Nsym        number of symbols used for estimators. Default %d\n", FSK_DEFAULT_NSYM);
+        fprintf(stderr," --mask TxFreqSpace Use \"mask\" freq estimator (default is \"peak\" estimator)\n");
+        exit(1);
+    }
+    
+    /* Extract parameters */
+    M = atoi(argv[dx]);
+    Fs = atoi(argv[dx + 1]);
+    Rs = atoi(argv[dx + 2]);
+    
+    if( (M!=2) && (M!=4) ){
+        fprintf(stderr,"Mode %d is not valid. Mode must be 2 or 4.\n",M);
+        goto helpmsg;
+    }
+    
+    /* Open files */
+    if(strcmp(argv[dx + 3],"-")==0){
+        fin = stdin;
+    }else{
+        fin = fopen(argv[dx + 3],"r");
+    }
+    
+    if(strcmp(argv[dx + 4],"-")==0){
+        fout = stdout;
+    }else{
+        fout = fopen(argv[dx + 4],"w");
+    }
+
+    /* set up FSK */
+    #define UNUSED 1000
+    fsk = fsk_create_hbr(Fs,Rs,M,P,nsym,UNUSED,tx_tone_separation);
+
+    /* set freq estimator limits */
+    if (!user_fsk_lower) {
+        if (complex_input == 1)
+            fsk_lower = 0;
+        else
+            fsk_lower = -Fs/2;
+    }
+    if (!user_fsk_upper) {
+        fsk_upper = Fs/2;
+    }
+    fprintf(stderr,"Setting estimator limits to %d to %d Hz.\n", fsk_lower, fsk_upper);
+    fsk_set_freq_est_limits(fsk,fsk_lower,fsk_upper);
+
+    fsk_set_freq_est_alg(fsk, mask);
+    
+    if(fin==NULL || fout==NULL || fsk==NULL){
+        fprintf(stderr,"Couldn't open files\n");
+        exit(1);
+    }
+
+    /* set up testframe mode */
+         
+    int      testframecnt, bitcnt, biterr, testframe_detected;
+    uint8_t *bitbuf_tx = NULL, *bitbuf_rx = NULL;
+    if (testframe_mode) {
+        bitbuf_tx = (uint8_t*)malloc(sizeof(uint8_t)*TEST_FRAME_SIZE); assert(bitbuf_tx != NULL);
+        bitbuf_rx = (uint8_t*)malloc(sizeof(uint8_t)*TEST_FRAME_SIZE); assert(bitbuf_rx != NULL);
+    
+        /* Generate known tx frame from known seed */
+        
+        srand(158324);
+        for(i=0; i<TEST_FRAME_SIZE; i++){
+            bitbuf_tx[i] = rand()&0x1;
+            bitbuf_rx[i] = 0;
+        }
+
+        testframecnt = 0;
+        bitcnt = 0;
+        biterr = 0;
+    }
+    
+    if(enable_stats){
+        loop_time = ((float)fsk_nin(fsk))/((float)Fs);
+        stats_loop = (int)(1/(stats_rate*loop_time));
+        stats_ctr = 0;
+    }
+    
+    /* allocate buffers for processing */
+    if(soft_dec_mode){
+        sdbuf = (float*)malloc(sizeof(float)*fsk->Nbits); assert(sdbuf != NULL);
+    }else{
+        bitbuf = (uint8_t*)malloc(sizeof(uint8_t)*fsk->Nbits); assert(bitbuf != NULL);
+    }
+    rawbuf = (int16_t*)malloc(bytes_per_sample*(fsk->N+fsk->Ts*2)*complex_input);
+    modbuf = (COMP*)malloc(sizeof(COMP)*(fsk->N+fsk->Ts*2));
+
+    /* set up signal handler so we can terminate gracefully */
+    
+    if (signal(SIGTERM, sig_handler) == SIG_ERR) {
+        printf("\ncan't catch SIGTERM\n");
+    }
+    
+    /* Demodulate! */
+    
+    while( fread(rawbuf,bytes_per_sample*complex_input,fsk_nin(fsk),fin) == fsk_nin(fsk) ){
+        /* convert input to a buffer of floats.  Note scaling isn't really necessary for FSK */
+
+        if (complex_input == 1) {
+            /* S16 real input */
+            for(i=0;i<fsk_nin(fsk);i++){
+                modbuf[i].real = ((float)rawbuf[i])/FDMDV_SCALE;
+                modbuf[i].imag = 0.0;
+            }
+        }
+        else {
+            if (bytes_per_sample == 1) {
+                /* U8 complex */
+                uint8_t *rawbuf_u8 = (uint8_t*)rawbuf;
+                for(i=0;i<fsk_nin(fsk);i++){
+                    modbuf[i].real = ((float)rawbuf_u8[2*i]-127.0)/128.0;
+                    modbuf[i].imag = ((float)rawbuf_u8[2*i+1]-127.0)/128.0;
+                }
+            }
+            else {
+                /* S16 complex */
+                for(i=0;i<fsk_nin(fsk);i++){
+                    modbuf[i].real = ((float)rawbuf[2*i])/FDMDV_SCALE;
+                    modbuf[i].imag = ((float)rawbuf[2*i+1]/FDMDV_SCALE);
+                }
+            }            
+        }
+
+        if(soft_dec_mode){
+            fsk_demod_sd(fsk,sdbuf,modbuf);
+        }else{
+            fsk_demod(fsk,bitbuf,modbuf);
+        }
+        
+        testframe_detected = 0;
+        if (testframe_mode) {
+            /* attempt to find a testframe and update stats */
+            /* update silding window of input bits */
+
+            int errs;
+            for(j=0; j<fsk->Nbits; j++) {
+                for(i=0; i<TEST_FRAME_SIZE-1; i++) {
+                    bitbuf_rx[i] = bitbuf_rx[i+1];
+                }
+                if (soft_dec_mode == 1) {
+                    bitbuf_rx[TEST_FRAME_SIZE-1] = sdbuf[j] < 0.0;
+                }
+                else {
+                    bitbuf_rx[TEST_FRAME_SIZE-1] = bitbuf[j];
+                }
+
+                /* compare to know tx frame.  If they are time aligned, there
+                   will be a fairly low bit error rate */
+
+                errs = 0;
+                for(i=0; i<TEST_FRAME_SIZE; i++) {
+                    if (bitbuf_rx[i] != bitbuf_tx[i]) {
+                        errs++;
+                    }
+                }
+                
+                if (errs < 0.1*TEST_FRAME_SIZE) {
+                    /* OK, we have a valid test frame sync, so lets count errors */
+                    testframe_detected = 1;
+                    testframecnt++;
+                    bitcnt += TEST_FRAME_SIZE;
+                    biterr += errs;
+                    if (enable_stats == 0) {
+                        fprintf(stderr,"errs: %d FSK BER %f, bits tested %d, bit errors %d\n",
+                                errs, ((float)biterr/(float)bitcnt),bitcnt,biterr);
+                    }
+                }
+            }
+        } /* if (testframe_mode) ... */
+        
+        if (enable_stats) {
+            if ((stats_ctr < 0) || testframe_detected) {
+                fsk_get_demod_stats(fsk,&stats);
+
+                /* Print standard 2FSK stats */
+
+                fprintf(stderr,"{");
+                time_t seconds  = time(NULL);
+
+                fprintf(stderr,"\"secs\": %ld, \"EbNodB\": %5.1f, \"ppm\": %4d,",seconds, stats.snr_est, (int)fsk->ppm);
+                float *f_est;
+                if (fsk->freq_est_type)
+                    f_est = fsk->f2_est;
+                else
+                    f_est = fsk->f_est;
+                fprintf(stderr," \"f1_est\":%.1f, \"f2_est\":%.1f",f_est[0],f_est[1]);
+
+                /* Print 4FSK stats if in 4FSK mode */
+
+                if(fsk->mode == 4){
+                    fprintf(stderr,", \"f3_est\":%.1f, \"f4_est\":%.1f",f_est[2],f_est[3]);
+                }
+	    
+                if (testframe_mode == 0) {
+                    /* Print the eye diagram */
+
+                    fprintf(stderr,",\t\"eye_diagram\":[");                 
+                    for(i=0;i<stats.neyetr;i++){
+                        fprintf(stderr,"[");
+                        for(j=0;j<stats.neyesamp;j++){
+                            fprintf(stderr,"%f ",stats.rx_eye[i][j]);
+                            if(j<stats.neyesamp-1) fprintf(stderr,",");
+                        }
+                        fprintf(stderr,"]");
+                        if(i<stats.neyetr-1) fprintf(stderr,",");
+                    }
+                    fprintf(stderr,"],");
+	    
+                    /* Print a sample of the FFT from the freq estimator */
+                    fprintf(stderr,"\"samp_fft\":[");
+                    Ndft = fsk->Ndft/2;
+                    for(i=0; i<Ndft; i++){
+                        fprintf(stderr,"%f ",(fsk->Sf)[i]);
+                        if(i<Ndft-1) fprintf(stderr,",");
+                    }
+                    fprintf(stderr,"]");
+                }
+                
+                if (testframe_mode) {
+                    fprintf(stderr,", \"frames\":%d, \"bits\":%d, \"errs\":%d",testframecnt,bitcnt,biterr);
+                }
+                
+                fprintf(stderr,"}\n");                
+
+                if (stats_ctr < 0) {
+                    stats_ctr = stats_loop;
+                }
+            }
+            if (testframe_mode == 0) {
+                stats_ctr--;
+            }
+        }
+
+        if(soft_dec_mode){
+            fwrite(sdbuf,sizeof(float),fsk->Nbits,fout);
+        }else{
+            fwrite(bitbuf,sizeof(uint8_t),fsk->Nbits,fout);
+        }
+
+        if(fin == stdin || fout == stdin){
+            fflush(fin);
+            fflush(fout);
+        }
+    } /* while(fread ...... */
+
+    if (testframe_mode) {
+        free(bitbuf_tx);
+        free(bitbuf_rx);
+    }
+    
+    if(soft_dec_mode){
+        free(sdbuf);
+    }else{
+        free(bitbuf);
+    }
+    
+    free(rawbuf);
+    free(modbuf);
+    
+    fclose(fin);
+    fclose(fout);
+    fsk_destroy(fsk);
+
+    return 0;
+}
+

src/fsk_get_test_bits.c

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: fsk_get_test_bits.c
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: January 2016
+
+  Generates a pseudorandom sequence of bits for testing of fsk_mod and fsk_demod
+
+\*---------------------------------------------------------------------------*/
+
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <stdio.h>
+#include <string.h>
+#include "fsk.h"
+
+#define TEST_FRAME_SIZE 100  /* arbitrary chice, repeats after this
+                                many bits, sets frame size for rx
+                                processing */
+
+int main(int argc,char *argv[]){
+    int bitcnt, framecnt;
+    int framesize = TEST_FRAME_SIZE;
+    int i;
+    FILE *fout;
+    uint8_t *bitbuf;
+    
+    if(argc < 3){
+        fprintf(stderr,"usage: %s OutputBitsOnePerByte numBits [framesize]\n",argv[0]);
+        exit(1);
+    }
+
+    if (argc == 4){
+        framesize = atoi(argv[3]);
+        fprintf(stderr, "Using custom frame size of %d bits\n", framesize);
+    }
+    
+    /* Extract parameters */
+    bitcnt = atoi(argv[2]);
+    framecnt = bitcnt/framesize;
+    if (framecnt == 0) {
+        fprintf(stderr,"Need a minimum of %d bits\n", framesize);
+        exit(1);
+    }
+
+    if(strcmp(argv[1],"-")==0){
+        fout = stdout;
+    }else{
+        fout = fopen(argv[1],"w");
+    }
+    
+    if(fout==NULL){
+        fprintf(stderr,"Couldn't open output file: %s\n", argv[1]);
+        goto cleanup;
+    }
+    
+    /* allocate buffers for processing */
+    bitbuf = (uint8_t*)alloca(sizeof(uint8_t)*framesize);
+    
+    /* Generate buffer of test frame bits from known seed */
+    srand(158324);
+    for(i=0; i<framesize; i++){
+	bitbuf[i] = rand()&0x1;
+    }
+        
+    /* Output test frames */
+    srand(158324);
+    for(i=0; i<framecnt; i++){
+	fwrite(bitbuf,sizeof(uint8_t),framesize,fout);
+	if(fout == stdout){
+	    fflush(fout);
+	}
+    }
+    
+ cleanup:
+    fclose(fout);
+
+    return 0;
+}

src/fsk_mod.c

@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: fsk_mod.c
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 8 January 2016
+
+  C test driver for fsk_mod in fsk.c. Reads in a set of bits to modulate
+   from a file, passed as a parameter, and writes modulated output to
+   another file
+   
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdio.h>
+#include <string.h>
+#include <getopt.h>
+#include "fsk.h"
+#include "codec2_fdmdv.h"
+
+int main(int argc,char *argv[]){
+    struct FSK *fsk;
+    int Fs,Rs,f1,fs,M;
+    int i;
+    int p, user_p = 0;
+    FILE *fin,*fout;
+    uint8_t *bitbuf;
+    int16_t *rawbuf;
+    float *modbuf;
+    
+    char usage[] = "usage: %s [-p P] Mode SampleFreq SymbolFreq TxFreq1 TxFreqSpace InputOneBitPerCharFile OutputModRawFile\n";
+
+    int opt;
+    while ((opt = getopt(argc, argv, "p:")) != -1) {
+        switch (opt) {
+        case 'p':
+            p = atoi(optarg);
+            user_p = 1;
+            break;
+        default:
+            fprintf(stderr, usage, argv[0]);
+            exit(1);
+        }
+    }
+
+    if (argc<8){
+        fprintf(stderr, usage, argv[0]);
+        exit(1);
+    }
+    
+    /* Extract parameters */
+    M = atoi(argv[optind++]);
+    Fs = atoi(argv[optind++]);
+    Rs = atoi(argv[optind++]);
+    f1 = atoi(argv[optind++]);
+    fs = atoi(argv[optind++]);
+    
+    if(strcmp(argv[optind],"-")==0){
+        fin = stdin;
+    }else{
+        fin = fopen(argv[optind],"r");
+    }
+    optind++;
+    
+    if(strcmp(argv[optind],"-")==0){
+        fout = stdout;
+    }else{
+        fout = fopen(argv[optind],"w");
+    }
+
+    /* p is not actually used for the modulator, but we need to set it for fsk_create() to be happy */    
+    if (!user_p)
+        p = Fs/Rs;
+    
+    /* set up FSK */
+    fsk = fsk_create_hbr(Fs,Rs,M,p,FSK_DEFAULT_NSYM,f1,fs);
+    
+    if(fin==NULL || fout==NULL || fsk==NULL){
+        fprintf(stderr,"Couldn't open test vector files\n");
+        goto cleanup;
+    }
+    
+    
+    /* allocate buffers for processing */
+    bitbuf = (uint8_t*)malloc(sizeof(uint8_t)*fsk->Nbits);
+    rawbuf = (int16_t*)malloc(sizeof(int16_t)*fsk->N);
+    modbuf = (float*)malloc(sizeof(float)*fsk->N);
+    
+    /* Modulate! */
+    while( fread(bitbuf,sizeof(uint8_t),fsk->Nbits,fin) == fsk->Nbits ){
+        fsk_mod(fsk,modbuf,bitbuf);
+        for(i=0; i<fsk->N; i++){
+			rawbuf[i] = (int16_t)(modbuf[i]*(float)FDMDV_SCALE);
+		}
+        fwrite(rawbuf,sizeof(int16_t),fsk->N,fout);
+        
+		if(fin == stdin || fout == stdin){
+			fflush(fin);
+			fflush(fout);
+		}
+    }
+    free(bitbuf);
+    free(rawbuf);
+    free(modbuf);
+    
+    cleanup:
+    fclose(fin);
+    fclose(fout);
+    fsk_destroy(fsk);
+    exit(0);
+}

src/fsk_put_test_bits.c

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: fsk_get_test_bits.c
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: January 2016
+
+  Generates a pseudorandom sequence of bits for testing of fsk_mod and
+  fsk_demod.
+
+\*---------------------------------------------------------------------------*/
+
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdio.h>
+#include <string.h>
+#include <getopt.h>
+#include "fsk.h"
+
+#define TEST_FRAME_SIZE 100  /* must match fsk_get_test_bits.c */
+
+#define VALID_PACKET_BER_THRESH 0.1
+
+int main(int argc,char *argv[]){
+    int bitcnt,biterr,i,errs,packetcnt;
+    int framesize = TEST_FRAME_SIZE;
+    float valid_packet_ber_thresh = VALID_PACKET_BER_THRESH;
+    int packet_pass_thresh = 0;
+    float ber_pass_thresh = 0;
+    FILE *fin;
+    uint8_t *bitbuf_tx, *bitbuf_rx, abit;
+    int verbose = 1;
+    
+    char usage[] = "usage: %s [-f frameSizeBits] [-t VaildFrameBERThreshold] [-b BERPass] [-p numPacketsPass] InputOneBitPerByte\n";
+
+    int opt;
+    while ((opt = getopt(argc, argv, "f:b:p:hqt:")) != -1) {
+        switch (opt) {
+        case 't':
+            valid_packet_ber_thresh = atof(optarg);
+            break;
+          case 'b':
+            ber_pass_thresh = atof(optarg);
+            break;
+        case 'p':
+            packet_pass_thresh = atoi(optarg);
+            break;
+        case 'f':
+            framesize = atoi(optarg);
+            break;
+        case 'q':
+            verbose = 0;
+            break;
+        case 'h':
+        default:
+            fprintf(stderr, usage, argv[0]);
+            exit(1);
+        }
+    }
+    if (argc == 1) {
+        fprintf(stderr, usage, argv[0]);
+        exit(1);
+    }
+    char *fname = argv[optind++];
+    if ((strcmp(fname,"-")==0) || (argc<2)){
+        fin = stdin;
+    } else {
+        fin = fopen(fname,"r");
+    }
+    
+    if(fin==NULL){
+        fprintf(stderr,"Couldn't open input file: %s\n", argv[1]);
+        exit(1);
+    }
+
+    /* allocate buffers for processing */
+    bitbuf_tx = (uint8_t*)alloca(sizeof(uint8_t)*framesize);
+    bitbuf_rx = (uint8_t*)alloca(sizeof(uint8_t)*framesize);
+    
+    /* Generate known tx frame from known seed */
+    srand(158324);
+    for(i=0; i<framesize; i++){
+	bitbuf_tx[i] = rand()&0x1;
+	bitbuf_rx[i] = 0;
+    }
+    
+    bitcnt = 0; biterr = 0; packetcnt = 0;
+    float ber = 0.5;
+    
+    while(fread(&abit,sizeof(uint8_t),1,fin)>0){
+
+        /* update sliding window of input bits */
+
+        for(i=0; i<framesize-1; i++) {
+            bitbuf_rx[i] = bitbuf_rx[i+1];
+        }
+        bitbuf_rx[framesize-1] = abit;
+
+        /* compare to know tx frame.  If they are time aligned, there
+           will be a fairly low bit error rate */
+
+        errs = 0;
+        for(i=0; i<framesize; i++) {
+            if (bitbuf_rx[i] != bitbuf_tx[i]) {
+                errs++;
+            }
+        }
+        if (errs < valid_packet_ber_thresh * framesize) {
+            /* OK, we have a valid packet, so lets count errors */
+            packetcnt++;
+            bitcnt += framesize;
+            biterr += errs;
+            ber =  (float)biterr/(float)bitcnt;
+            if (verbose) {
+                fprintf(stderr,"[%04d] BER %5.3f, bits tested %6d, bit errors %6d errs: %4d \n",
+                        packetcnt, ber, bitcnt, biterr, errs);
+            }
+        }
+    }
+ 
+    fclose(fin);
+
+    fprintf(stderr,"[%04d] BER %5.3f, bits tested %6d, bit errors %4d\n", packetcnt, ber, bitcnt, biterr);
+    if ((packetcnt >= packet_pass_thresh) && (ber <= ber_pass_thresh)) {
+        fprintf(stderr,"PASS\n");
+        return 0;
+    }
+    else {
+        fprintf(stderr,"FAIL\n");
+        return 1;
+    }
+}

src/golay23.c

@@ -0,0 +1,311 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: golay23.c
+  AUTHOR......: Tomas Härdin & David Rowe
+  DATE CREATED: 3 March 2013
+
+  To test:
+
+     src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_UNITTEST                   && ./golay23
+     src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_UNITTEST -DRUN_TIME_TABLES && ./golay23
+     src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_UNITTEST -DNO_TABLES       && ./golay23
+
+  To generate tables:
+     src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_MAKETABLES                 && ./golay23
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 Tomas Härdin & David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not,see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+
+#ifdef GOLAY23_MAKETABLES
+#define RUN_TIME_TABLES
+#endif
+
+#ifndef NO_TABLES
+#ifdef RUN_TIME_TABLES
+int static encoding_table[4096];
+int static decoding_table[2048];
+static int inited = 0;
+#else
+//default is to use precomputed tables
+#include "golayenctable.h"
+#include "golaydectable.h"
+#endif
+#endif
+
+//since we want to avoid bit-reversing inside syndrome() we bit-reverse the polynomial instead
+#define GOLAY_POLYNOMIAL    0xC75   //AE3 reversed
+
+int golay23_syndrome(int c) {
+    //could probably be done slightly smarter, but works
+    int x;
+    for (x = 11; x >= 0; x--) {
+        if (c & ((1<<11) << x)) {
+            c ^= GOLAY_POLYNOMIAL << x;
+        }
+    }
+    return c;
+}
+
+#ifdef __GNUC__
+#define popcount __builtin_popcount
+#elif defined(__MSC_VER)
+#include <intrin.h>
+#define popcount __popcnt
+#else
+static int popcount(unsigned int c) {
+    int ret = 0;
+    while (c) {
+        if (c & 1) {
+            ret++;
+        }
+        c >>= 1;
+    }
+    return ret;
+}
+#endif
+
+#if defined(NO_TABLES) || defined(RUN_TIME_TABLES)
+static int golay23_encode_no_tables(int c) {
+    c <<= 11;
+    return golay23_syndrome(c) | c;
+}
+#endif
+
+#ifdef NO_TABLES
+static int unrotate(unsigned int c, int x) {
+    return ((c << x) & 0x7FFFFF) | (c >> (23 - x));
+}
+
+static int golay23_decode_no_tables(int c) {
+    //TODO: optimize?
+    int x;
+    c = unrotate(c, 12);
+
+    for (x = 0; x < 23; x++) {
+        int t;
+        int s = golay23_syndrome(c);
+
+        if (popcount(s) <= 3) {
+            return unrotate(c ^ s, x) & 0xFFF;
+        }
+
+        for (t = 0; t < 23; t++) {
+            int c2 = c ^ (1 << t);
+            int s = golay23_syndrome(c2);
+
+            if (popcount(s) <= 2) {
+                return unrotate(c2 ^ s, x) & 0xFFF;
+            }
+        }
+
+        //rotate
+        c = (c >> 1) | ((c & 1) << 22);
+    }
+
+    //shouldn't reach here..
+    assert("Something is wrong with golay23_decode_no_tables()..");
+    return c & 0xFFF;
+}
+#endif
+
+void golay23_init(void) {
+#ifdef RUN_TIME_TABLES
+    int x, y, z;
+    inited = 1;
+    for (x = 0; x < 4096; x++) {
+        encoding_table[x] = golay23_encode_no_tables(x);
+    }
+
+    decoding_table[0] = 0;
+    //1-bit errors
+    for (x = 0; x < 23; x++) {
+        int d = 1<<x;
+        decoding_table[golay23_syndrome(d)] = d;
+    }
+    //2-bit errors
+    for (x = 0; x < 22; x++) {
+        for (y = x+1; y < 23; y++) {
+            int d = (1<<x) | (1<<y);
+            decoding_table[golay23_syndrome(d)] = d;
+        }
+    }
+    //3-bit errors
+    for (x = 0; x < 21; x++) {
+        for (y = x+1; y < 22; y++) {
+            for (z = y+1; z < 23; z++) {
+                int d = (1<<x) | (1<<y) | (1<<z);
+                decoding_table[golay23_syndrome(d)] = d;
+            }
+        }
+    }
+#endif
+}
+
+int  golay23_encode(int c) {
+    assert(c >= 0 && c <= 0xFFF);
+#ifdef RUN_TIME_TABLES
+    assert(inited);
+#endif
+
+#ifdef NO_TABLES
+    return golay23_encode_no_tables(c);
+#else
+    return encoding_table[c];
+#endif
+}
+
+int  golay23_decode(int c) {
+    assert(c >= 0 && c <= 0x7FFFFF);
+#ifdef RUN_TIME_TABLES
+    assert(inited);
+#endif
+
+#ifdef NO_TABLES
+    //duplicate old golay23_decode()'s shift
+    return unrotate(golay23_decode_no_tables(c), 11);
+#else
+    //message is shifted 11 places left in the return value
+    return c ^ decoding_table[golay23_syndrome(c)];
+#endif
+}
+
+int  golay23_count_errors(int recd_codeword, int corrected_codeword) {
+    return popcount(recd_codeword ^ corrected_codeword);
+}
+
+/**
+ * Table generation and testing code below
+ */
+
+#ifdef GOLAY23_MAKETABLES
+#include <stdio.h>
+
+int main() {
+    int x;
+    //generate and dump
+    golay23_init();
+
+    FILE *enc = fopen("golayenctable.h", "w");
+    FILE *dec = fopen("golaydectable.h", "w");
+
+    fprintf(enc, "/* Generated by golay23.c -DGOLAY23_MAKETABLE */\n\
+\n\
+const int static encoding_table[]={\n");
+    for (x = 0; x < 4096; x++) {
+        fprintf(enc, x < 4095 ? "  0x%x,\n" : "  0x%x\n", encoding_table[x]);
+    }
+    fprintf(enc, "};\n");
+
+    fprintf(dec, "/* Generated by golay23.c -DGOLAY23_MAKETABLE */\n\
+\n\
+const int static decoding_table[]={\n");
+    for (x = 0; x < 2048; x++) {
+        fprintf(dec, x < 2047 ? "  0x%x,\n" : "  0x%x\n", decoding_table[x]);
+    }
+    fprintf(dec, "};\n");
+
+    fclose(enc);
+    fclose(dec);
+
+    return 0;
+}
+
+#elif defined(GOLAY23_UNITTEST)
+#include <stdio.h>
+#include <stdlib.h>
+#include <memory.h>
+
+int main() {
+    int c;
+
+    golay23_init();
+
+    //keep track of whether every single codeword has been checked
+    char *checkmask = malloc(1<<23);
+    memset(checkmask, 0, 1<<23);
+
+    //step through all possible messages
+    for (c = 0; c < (1<<12); c++) {
+        int g23 = golay23_encode(c);
+        int x,y,z;
+        checkmask[g23] = 1;
+        int c2 = golay23_decode(g23) >> 11;
+
+        printf("%03x -> %06x %03x\n", c, g23, c2);
+
+        if (c != c2) {
+            printf("Bad!\n");
+            exit(1);
+        }
+
+        //test the code by flipping every combination of one, two and three bits
+        for (x = 0; x < 23; x++) {
+            int flipped = g23 ^ (1<<x);
+            checkmask[flipped] = 1;
+            int c2 = golay23_decode(flipped) >> 11;
+            if (c != c2) {
+                printf("Bad!\n");
+                
+                exit(1);
+            }
+        }
+        
+        for (x = 0; x < 22; x++) {
+            for (y = x+1; y < 23; y++) {
+                int flipped = g23 ^ (1<<x) ^ (1<<y);
+                checkmask[flipped] = 1;
+                int c2 = golay23_decode(flipped) >> 11;
+                if (c != c2) {
+                    printf("Bad!\n");
+                    
+                    exit(1);
+                }
+            }
+        }
+        
+        for (x = 0; x < 21; x++) {
+            for (y = x+1; y < 22; y++) {
+                for (z = y+1; z < 23; z++) {
+                    int flipped = g23 ^ (1<<x) ^ (1<<y) ^ (1<<z);
+                    checkmask[flipped] = 1;
+                    int c2 = golay23_decode(flipped) >> 11;
+                    if (c != c2) {
+                        printf("Bad!\n");
+                        exit(1);
+                    }
+                }
+            }
+        }
+    }
+
+    //did we check every codeword?
+    for (c = 0; c < (1<<23); c++) {
+        if (checkmask[c] != 1) {
+            printf("%06x unchecked!\n", c);
+            exit(1);
+        }
+    }
+
+    printf("Everything checks out\n");
+    free(checkmask);
+    return 0;
+}
+#endif

src/golay23.h

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: golay23.h
+  AUTHOR......: David Rowe
+  DATE CREATED: 3 March 2013
+
+  Header file for Golay FEC.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2013 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __GOLAY23__
+#define __GOLAY23__
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void golay23_init(void);
+int  golay23_encode(int data);
+int  golay23_decode(int received_codeword);
+int  golay23_count_errors(int recd_codeword, int corrected_codeword);
+int  golay23_syndrome(int c);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

src/horus_api.c

@@ -0,0 +1,516 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: horus_api.c
+  AUTHOR......: David Rowe
+  DATE CREATED: March 2018
+
+  Library of API functions that implement High Altitude Balloon (HAB)
+  telemetry modems and protocols for Project Horus.  May also be useful for
+  other HAB projects.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2018 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "horus_api.h"
+#include "fsk.h"
+#include "horus_l2.h"
+
+#define MAX_UW_LENGTH                  100
+#define HORUS_API_VERSION                1    /* unique number that is bumped if API changes */
+#define HORUS_BINARY_NUM_BITS          360    /* fixed number of bytes in binary payload     */
+#define HORUS_BINARY_NUM_PAYLOAD_BYTES  22    /* fixed number of bytes in binary payload     */
+
+struct horus {
+    int         mode;
+    int         verbose;
+    struct FSK *fsk;                 /* states for FSK modem                */
+    int         Fs;                  /* sample rate in Hz                   */
+    int         mFSK;                /* number of FSK tones                 */
+    int         Rs;                  /* symbol rate in Hz                   */
+    int         uw[MAX_UW_LENGTH];   /* unique word bits mapped to +/-1     */
+    int         uw_thresh;           /* threshold for UW detection          */
+    int         uw_len;              /* length of unique word               */
+    int         max_packet_len;      /* max length of a telemetry packet    */
+    uint8_t    *rx_bits;             /* buffer of received bits             */
+    int         rx_bits_len;         /* length of rx_bits buffer            */
+    int         crc_ok;              /* most recent packet checksum results */
+    int         total_payload_bits;  /* num bits rx-ed in last RTTY packet  */
+};
+
+/* Unique word for Horus RTTY 7 bit '$' character, 3 sync bits,
+   repeated 5 times */
+
+int8_t uw_horus_rtty[] = {
+  0,0,1,0,0,1,0,1,1,0,
+  0,0,1,0,0,1,0,1,1,0,
+  0,0,1,0,0,1,0,1,1,0,
+  0,0,1,0,0,1,0,1,1,0,
+  0,0,1,0,0,1,0,1,1,0
+};
+
+/* Unique word for Horus Binary */
+
+int8_t uw_horus_binary[] = {
+    0,0,1,0,0,1,0,0,
+    0,0,1,0,0,1,0,0 
+};
+
+struct horus *horus_open (int mode, int Rs) {
+    int i;
+    assert((mode == HORUS_MODE_RTTY) || (mode == HORUS_MODE_BINARY));
+
+    struct horus *hstates = (struct horus *)malloc(sizeof(struct horus));
+    assert(hstates != NULL);
+
+    hstates->Fs = 48000; hstates->Rs = Rs; hstates->verbose = 0; hstates->mode = mode;
+
+    if (mode == HORUS_MODE_RTTY) {
+        hstates->mFSK = 2;
+        hstates->max_packet_len = 1000;
+
+        /* map UW to make it easier to search for */
+
+        for (i=0; i<sizeof(uw_horus_rtty); i++) {
+            hstates->uw[i] = 2*uw_horus_rtty[i] - 1;
+        }        
+        hstates->uw_len = sizeof(uw_horus_rtty);
+        hstates->uw_thresh = sizeof(uw_horus_rtty) - 2;  /* allow a few bit errors in UW detection */
+        hstates->rx_bits_len = hstates->max_packet_len;
+    }
+
+    if (mode == HORUS_MODE_BINARY) {
+        hstates->mFSK = 4;
+        hstates->max_packet_len = HORUS_BINARY_NUM_BITS;
+        for (i=0; i<sizeof(uw_horus_binary); i++) {
+            hstates->uw[i] = 2*uw_horus_binary[i] - 1;
+        }
+        hstates->uw_len = sizeof(uw_horus_binary);
+        hstates->uw_thresh = sizeof(uw_horus_binary) - 2; /* allow a few bit errors in UW detection */
+        horus_l2_init();
+        hstates->rx_bits_len = hstates->max_packet_len;
+    }
+   
+    hstates->fsk = fsk_create(hstates->Fs, hstates->Rs, hstates->mFSK, 1000, 2*hstates->Rs);
+
+    /* allocate enough room for two packets so we know there will be
+       one complete packet if we find a UW at start */
+    
+    hstates->rx_bits_len += hstates->fsk->Nbits;
+    hstates->rx_bits = (uint8_t*)malloc(hstates->rx_bits_len);
+    assert(hstates->rx_bits != NULL);
+    for(i=0; i<hstates->rx_bits_len; i++) {
+        hstates->rx_bits[i] = 0;
+    }
+
+    hstates->crc_ok = 0;
+    hstates->total_payload_bits = 0;
+    
+    return hstates;
+}
+
+void horus_close (struct horus *hstates) {
+    assert(hstates != NULL);
+    fsk_destroy(hstates->fsk);
+    free(hstates->rx_bits);
+    free(hstates);
+}
+
+uint32_t horus_nin(struct horus *hstates) {
+    assert(hstates != NULL);
+    int nin = fsk_nin(hstates->fsk);
+    assert(nin <= horus_get_max_demod_in(hstates));
+    return nin;
+}
+
+int horus_find_uw(struct horus *hstates, int n) {
+    int i, j, corr, mx, mx_ind;
+    int rx_bits_mapped[n+hstates->uw_len];
+    
+    /* map rx_bits to +/-1 for UW search */
+
+    for(i=0; i<n+hstates->uw_len; i++) {
+        rx_bits_mapped[i] = 2*hstates->rx_bits[i] - 1;
+    }
+    
+    /* look for UW  */
+
+    mx = 0; mx_ind = 0;
+    for(i=0; i<n; i++) {
+
+        /* calculate correlation between bit stream and UW */
+        
+        corr = 0;
+        for(j=0; j<hstates->uw_len; j++) {
+            corr += rx_bits_mapped[i+j]*hstates->uw[j];
+        }
+        
+        /* peak pick maximum */
+        
+        if (corr > mx) {
+            mx = corr;
+            mx_ind = i;
+        }
+    }
+
+    if (hstates->verbose) {
+        fprintf(stderr, "  horus_find_uw: mx_ind: %d mx: %d uw_thresh: %d n: %d\n",  mx_ind, mx, hstates->uw_thresh, n);
+    }
+    
+    if (mx >= hstates->uw_thresh) {
+        return mx_ind;
+    } else {
+        return -1;
+    }
+}
+
+int hex2int(char ch) {
+    if (ch >= '0' && ch <= '9')
+        return ch - '0';
+    if (ch >= 'A' && ch <= 'F')
+        return ch - 'A' + 10;
+    if (ch >= 'a' && ch <= 'f')
+        return ch - 'a' + 10;
+    return -1;
+}
+
+
+int extract_horus_rtty(struct horus *hstates, char ascii_out[], int uw_loc) {
+    const int nfield = 7;                               /* 7 bit ASCII                    */
+    const int npad   = 3;                               /* 3 sync bits between characters */
+    int st = uw_loc;                                    /* first bit of first char        */
+    int en = hstates->max_packet_len - nfield;          /* last bit of max length packet  */
+
+    int      i, j, endpacket, nout, crc_ok;
+    uint8_t  char_dec;
+    char    *pout, *ptx_crc;
+    uint16_t rx_crc, tx_crc;
+
+    pout = ascii_out; nout = 0; crc_ok = 0; endpacket = 0; rx_crc = tx_crc = 0;
+    
+    for (i=st; i<en; i+=nfield+npad) {
+
+        /* assemble char LSB to MSB */
+
+        char_dec = 0;
+        for(j=0; j<nfield; j++) {
+            assert(hstates->rx_bits[i+j] <= 1);
+            char_dec |= hstates->rx_bits[i+j] * (1<<j);
+        }
+        if (hstates->verbose) {
+            fprintf(stderr, "  extract_horus_rtty i: %4d 0x%02x %c ", i, char_dec, char_dec);
+            if ((nout % 6) == 0) {
+                fprintf(stderr, "\n");
+            }
+        }
+
+        /*  if we find a '*' that's the end of the packet for RX CRC calculations */
+
+        if (!endpacket && (char_dec == 42)) {
+            endpacket = 1;
+            rx_crc = horus_l2_gen_crc16((uint8_t*)&ascii_out[5], nout-5);
+            ptx_crc = pout + 1; /* start of tx CRC */
+        }
+
+        /* build up output array, really only need up to tx crc but
+           may end up going further */
+        
+        *pout++ = (char)char_dec;
+        nout++;
+        
+    }
+
+    /* if we found the end of packet flag and have enough chars to compute checksum ... */
+
+    //fprintf(stderr, "\n\ntx CRC...\n");
+    if (endpacket && (pout > (ptx_crc+3))) {
+        tx_crc = 0;
+        for(i=0; i<4; i++) {
+            tx_crc <<= 4;
+            tx_crc |= hex2int(ptx_crc[i]);
+            //fprintf(stderr, "ptx_crc[%d] %c 0x%02X tx_crc: 0x%04X\n", i, ptx_crc[i], hex2int(ptx_crc[i]), tx_crc);
+        }
+        crc_ok = (tx_crc == rx_crc);
+        *(ptx_crc+4) = 0;  /* terminate ASCII string */
+
+        if (crc_ok) {
+            hstates->total_payload_bits = strlen(ascii_out)*7;
+        }
+    }
+    else {
+        *ascii_out = 0;
+    }
+
+    if (hstates->verbose) {
+        fprintf(stderr, "\n  endpacket: %d nout: %d tx_crc: 0x%04x rx_crc: 0x%04x\n",
+                endpacket, nout, tx_crc, rx_crc);
+    }
+            
+    /* make sure we don't overrun storage */
+    
+    assert(nout <= horus_get_max_ascii_out_len(hstates));
+
+    hstates->crc_ok = crc_ok;
+    
+    return crc_ok;
+}
+
+
+int extract_horus_binary(struct horus *hstates, char hex_out[], int uw_loc) {
+    const int nfield = 8;                      /* 8 bit binary                   */
+    int st = uw_loc;                           /* first bit of first char        */
+    int en = uw_loc + hstates->max_packet_len; /* last bit of max length packet  */
+
+    int      j, b, nout;
+    uint8_t  rxpacket[hstates->max_packet_len];
+    uint8_t  rxbyte, *pout;
+ 
+    /* convert bits to a packet of bytes */
+    
+    pout = rxpacket; nout = 0;
+    
+    for (b=st; b<en; b+=nfield) {
+
+        /* assemble bytes MSB to LSB */
+
+        rxbyte = 0;
+        for(j=0; j<nfield; j++) {
+            assert(hstates->rx_bits[b+j] <= 1);
+            rxbyte <<= 1;
+            rxbyte |= hstates->rx_bits[b+j];
+        }
+        
+        /* build up output array */
+        
+        *pout++ = rxbyte;
+        nout++;
+    }
+
+    if (hstates->verbose) {
+        fprintf(stderr, "  extract_horus_binary nout: %d\n  Received Packet before decoding:\n  ", nout);
+        for (b=0; b<nout; b++) {
+            fprintf(stderr, "%02X", rxpacket[b]);
+        }
+        fprintf(stderr, "\n");
+    }
+    
+    uint8_t payload_bytes[HORUS_BINARY_NUM_PAYLOAD_BYTES];
+    horus_l2_decode_rx_packet(payload_bytes, rxpacket, HORUS_BINARY_NUM_PAYLOAD_BYTES);
+
+    uint16_t crc_tx, crc_rx;
+    crc_rx = horus_l2_gen_crc16(payload_bytes, HORUS_BINARY_NUM_PAYLOAD_BYTES-2);
+    crc_tx = (uint16_t)payload_bytes[HORUS_BINARY_NUM_PAYLOAD_BYTES-2] +
+        ((uint16_t)payload_bytes[HORUS_BINARY_NUM_PAYLOAD_BYTES-1]<<8);
+    
+    if (hstates->verbose) {
+        fprintf(stderr, "  extract_horus_binary crc_tx: %04X crc_rx: %04X\n", crc_tx, crc_rx);
+    }
+    
+    /* convert to ASCII string of hex characters */
+
+    hex_out[0] = 0;
+    char hex[3];
+    for (b=0; b<HORUS_BINARY_NUM_PAYLOAD_BYTES; b++) {
+        sprintf(hex, "%02X", payload_bytes[b]);
+        strcat(hex_out, hex);
+    }
+   
+    if (hstates->verbose) {
+        fprintf(stderr, "  nout: %d Decoded Payload bytes:\n  %s", nout, hex_out);
+    }
+
+    /* With noise input to FSK demod we can get occasinal UW matches,
+       so a good idea to only pass on any packets that pass CRC */
+    
+    hstates->crc_ok = (crc_tx == crc_rx);
+    if ( hstates->crc_ok) {
+        hstates->total_payload_bits += HORUS_BINARY_NUM_PAYLOAD_BYTES;
+    }
+    return hstates->crc_ok;
+}
+
+
+int horus_rx(struct horus *hstates, char ascii_out[], short demod_in[], int quadrature) {
+    int i, j, uw_loc, packet_detected;
+    
+    assert(hstates != NULL);
+    packet_detected = 0;
+
+    int Nbits = hstates->fsk->Nbits;
+    int rx_bits_len = hstates->rx_bits_len;
+    
+    if (hstates->verbose) {
+        fprintf(stderr, "  horus_rx max_packet_len: %d rx_bits_len: %d Nbits: %d nin: %d\n",
+                hstates->max_packet_len, rx_bits_len, Nbits, hstates->fsk->nin);
+    }
+    
+    /* shift buffer of bits to make room for new bits */
+
+    for(i=0,j=Nbits; j<rx_bits_len; i++,j++) {
+        hstates->rx_bits[i] = hstates->rx_bits[j];
+    }
+                   
+    /* demodulate latest bits */
+
+    /* Note: allocating this array as an automatic variable caused OSX to
+       "Bus Error 10" (segfault), so lets malloc() it. */
+    
+    COMP *demod_in_comp = (COMP*)malloc(sizeof(COMP)*hstates->fsk->nin);
+    
+    for (i=0; i<hstates->fsk->nin; i++) {
+	if (quadrature) {
+	    demod_in_comp[i].real = demod_in[i * 2];
+	    demod_in_comp[i].imag = demod_in[i * 2 + 1];
+	} else {
+	    demod_in_comp[i].real = demod_in[i];
+	    demod_in_comp[i].imag = 0;
+	}
+    }
+    fsk_demod(hstates->fsk, &hstates->rx_bits[rx_bits_len-Nbits], demod_in_comp);
+    free(demod_in_comp);
+    
+    /* UW search to see if we can find the start of a packet in the buffer */
+    
+    if ((uw_loc = horus_find_uw(hstates, Nbits)) != -1) {
+
+        if (hstates->verbose) {
+            fprintf(stderr, "  horus_rx uw_loc: %d mode: %d\n", uw_loc, hstates->mode);
+        }
+        
+        /* OK we have found a unique word, and therefore the start of
+           a packet, so lets try to extract valid packets */
+
+        if (hstates->mode == HORUS_MODE_RTTY) {
+            packet_detected = extract_horus_rtty(hstates, ascii_out, uw_loc);
+        }
+        if (hstates->mode == HORUS_MODE_BINARY) {
+            packet_detected = extract_horus_binary(hstates, ascii_out, uw_loc);
+            //#define DUMP_BINARY_PACKET
+            #ifdef DUMP_BINARY_PACKET
+            FILE *f = fopen("packetbits.txt", "wt"); assert(f != NULL);
+            for(i=0; i<hstates->max_packet_len; i++) {
+                fprintf(f,"%d ", hstates->rx_bits[uw_loc+i]);
+            }
+            fclose(f);
+            exit(0);
+            #endif
+        }
+    }
+     
+    return packet_detected;
+}
+
+int horus_get_version(void) {
+    return HORUS_API_VERSION;
+}
+
+int horus_get_mode(struct horus *hstates) {
+    assert(hstates != NULL);
+    return hstates->mode;
+}
+
+int horus_get_Fs(struct horus *hstates) {
+    assert(hstates != NULL);
+    return hstates->Fs;
+}
+
+int horus_get_mFSK(struct horus *hstates) {
+    assert(hstates != NULL);
+    return hstates->mFSK;
+}
+
+int horus_get_max_demod_in(struct horus *hstates) {
+    /* copied from fsk_demod.c, a nicer fsk_max_nin function would be useful */
+    return sizeof(short)*(hstates->fsk->N + hstates->fsk->Ts*2);
+}
+
+int horus_get_max_ascii_out_len(struct horus *hstates) {
+    assert(hstates != NULL);
+    if (hstates->mode == HORUS_MODE_RTTY) {
+        return hstates->max_packet_len/10;     /* 7 bit ASCII, plus 3 sync bits */
+    }
+    if (hstates->mode == HORUS_MODE_BINARY) {
+        return (HORUS_BINARY_NUM_PAYLOAD_BYTES*2+1);     /* Hexadecimal encoded */
+        //return HORUS_BINARY_NUM_PAYLOAD_BYTES;
+    }
+    assert(0); /* should never get here */
+    return 0;
+}
+
+void horus_get_modem_stats(struct horus *hstates, int *sync, float *snr_est) {
+    struct MODEM_STATS stats;
+    assert(hstates != NULL);
+
+    /* TODO set sync if UW found "recently", but WTF is recently? Maybe need a little state 
+       machine to "blink" sync when we get a packet */
+
+    *sync = 0;
+    
+    /* SNR scaled from Eb/No est returned by FSK to SNR in 3000 Hz */
+
+    fsk_get_demod_stats(hstates->fsk, &stats);
+    *snr_est = stats.snr_est + 10*log10((float)hstates->Rs*log2(hstates->mFSK)/3000);
+}
+
+void horus_get_modem_extended_stats (struct horus *hstates, struct MODEM_STATS *stats) {
+    int i;
+    
+    assert(hstates != NULL);
+
+    fsk_get_demod_stats(hstates->fsk, stats);
+    if (hstates->verbose) {
+        fprintf(stderr, "  horus_get_modem_extended_stats stats->snr_est: %f\n", stats->snr_est);
+    }
+    stats->snr_est = stats->snr_est + 10*log10((float)hstates->Rs*log2(hstates->mFSK)/3000);
+
+    assert(hstates->mFSK <= MODEM_STATS_MAX_F_EST);
+    for (i=0; i<hstates->mFSK; i++) {
+        stats->f_est[i] = hstates->fsk->f_est[i];
+    }
+}
+
+void horus_set_verbose(struct horus *hstates, int verbose) {
+    assert(hstates != NULL);
+    hstates->verbose = verbose;
+}
+
+int horus_crc_ok(struct horus *hstates) {
+    assert(hstates != NULL);
+    return hstates->crc_ok;
+}
+
+int horus_get_total_payload_bits(struct horus *hstates) {
+    assert(hstates != NULL);
+    return hstates->total_payload_bits;
+}
+
+void horus_set_total_payload_bits(struct horus *hstates, int val) {
+    assert(hstates != NULL);
+    hstates->total_payload_bits = val;
+}
+
+void horus_set_freq_est_limits(struct horus *hstates, float fsk_lower, float fsk_upper) {
+    assert(hstates != NULL);
+    assert(fsk_upper > fsk_lower);
+    hstates->fsk->est_min = fsk_lower;
+    hstates->fsk->est_max = fsk_upper;    
+}

src/horus_api.h

@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: horus_api.h
+  AUTHOR......: David Rowe
+  DATE CREATED: March 2018
+
+  Library of API functions that implement High Altitude Balloon (HAB)
+  telemetry modems and protocols.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2018 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifdef __cplusplus
+  extern "C" {
+#endif
+
+#ifndef __HORUS_API__
+
+#include <stdint.h>
+#include "modem_stats.h"
+      
+#define HORUS_MODE_BINARY            0
+#define HORUS_MODE_RTTY              1
+
+struct horus;
+struct MODEM_STATS;
+
+struct horus *horus_open  (int mode, int Rs);
+void          horus_close (struct horus *hstates);
+
+/* call before horus_rx() to determine how many shorts to pass in */
+
+uint32_t      horus_nin   (struct horus *hstates);
+
+/* returns 1 if ascii_out[] is valid */
+      
+int           horus_rx    (struct horus *hstates, char ascii_out[], short demod_in[], int quadrature);
+
+/* set verbose level */
+      
+void horus_set_verbose(struct horus *hstates, int verbose);
+      
+/* functions to get information from API  */
+      
+int           horus_get_version              (void);
+int           horus_get_mode                 (struct horus *hstates);
+int           horus_get_Fs                   (struct horus *hstates);      
+int           horus_get_mFSK                 (struct horus *hstates);      
+void          horus_get_modem_stats          (struct horus *hstates, int *sync, float *snr_est);
+void          horus_get_modem_extended_stats (struct horus *hstates, struct MODEM_STATS *stats);
+int           horus_crc_ok                   (struct horus *hstates);
+int           horus_get_total_payload_bits   (struct horus *hstates);
+void          horus_set_total_payload_bits   (struct horus *hstates, int val);
+void          horus_set_freq_est_limits      (struct horus *hstates, float fsk_lower, float fsk_upper);
+
+/* how much storage you need for demod_in[] and  ascii_out[] */
+      
+int           horus_get_max_demod_in         (struct horus *hstates);
+int           horus_get_max_ascii_out_len    (struct horus *hstates);
+
+#endif
+
+#ifdef __cplusplus
+}
+#endif

src/horus_demod.c

@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: horus_demod.c
+  AUTHOR......: David Rowe
+  DATE CREATED: April 2018
+
+  Command line demo program for the Horus API, that exercises
+  horus_api.c using file input/output (can be stdin/stdout for real
+  time operation).  Prints JSON stats, just like Brady's fsk_demod.c
+ 
+  Can operate in Horus RTTY or Binary mode.
+
+  Testing with a 8000Hz sample rate wave file:
+
+    $ sox ~/Desktop/horus.wav -r 48000 -t raw - | ./horus_demod -m RTTY -v - /dev/null
+
+    $ sox ~/Desktop/4FSK_binary_100Rb_8khzfs.wav -r 48000 -t raw - | ./horus_demod -m binary  - -
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2018 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <getopt.h>
+
+#include "horus_api.h"
+#include "fsk.h"
+#include "horus_l2.h"
+
+int main(int argc, char *argv[]) {
+    struct   horus *hstates;
+    struct   MODEM_STATS stats;
+    FILE    *fin,*fout;
+    int      i,j,Ndft,mode;
+    int      stats_ctr,stats_loop, stats_rate, verbose, crc_results;
+    float    loop_time;
+    int      enable_stats = 0;
+    int      quadrature = 0;
+    int fsk_lower = -1;
+    int fsk_upper = -1;
+    int Rs = 100;
+
+    stats_loop = 0;
+    stats_rate = 8;
+    mode = -1;
+    verbose = crc_results = 0;
+    
+    int o = 0;
+    int opt_idx = 0;
+    while ( o != -1 ) {
+        static struct option long_opts[] = {
+            {"help",      no_argument,        0, 'h'},
+            {"mode",      required_argument,  0, 'm'},
+            {"rate",      optional_argument,  0, 'r'},
+            {"stats",     optional_argument,  0, 't'},
+            {"fsk_lower", optional_argument,  0, 'b'},
+            {"fsk_upper", optional_argument,  0, 'u'},
+            {0, 0, 0, 0}
+        };
+        
+        o = getopt_long(argc,argv,"hvcqm:t::",long_opts,&opt_idx);
+        
+        switch(o) {
+            case 'm':
+                if ((strcmp(optarg, "RTTY") == 0) || (strcmp(optarg, "rtty") == 0)) {
+                    mode = HORUS_MODE_RTTY;
+                }
+                if ((strcmp(optarg, "BINARY") == 0) || (strcmp(optarg, "binary") == 0)) {
+                    mode = HORUS_MODE_BINARY;
+                }
+                if (mode == -1) {
+                    fprintf(stderr, "use --mode RTTY or --mode binary\n");
+                    exit(1);
+                }
+                break;
+            case 't':
+                enable_stats = 1;
+                if (optarg != NULL){
+                    stats_rate = atoi(optarg);
+                    if (stats_rate == 0) {
+                        stats_rate = 8;
+                    }
+                }
+                break;
+            case 'v':
+                verbose = 1;
+                break;    
+            case 'c':
+                crc_results = 1;
+                break;    
+            case 'h':
+            case '?':
+                goto helpmsg;
+                break;
+            case 'q':
+                quadrature = 1;
+                break;
+            case 'b':
+                if (optarg != NULL){
+                    fsk_lower = atoi(optarg);
+                }
+                break;
+            case 'u':
+                if (optarg != NULL){
+                    fsk_upper = atoi(optarg);
+                }
+                break;
+            case 'r':
+                if (optarg != NULL){
+                    Rs = atoi(optarg);
+                }
+                break;
+            break;
+        }
+    }
+    
+    int dx = optind;
+    
+    if( (argc - dx) < 2) {
+        fprintf(stderr, "Too few arguments\n");
+        goto helpmsg;
+    }
+    
+    if( (argc - dx) > 5) {
+        fprintf(stderr, "Too many arguments\n");
+    helpmsg:
+        fprintf(stderr,"usage: %s -m RTTY|binary [-q] [-v] [-c] [-t [r]] InputModemRawFile OutputAsciiFile\n",argv[0]);
+        fprintf(stderr,"\n");
+        fprintf(stderr,"InputModemRawFile      48 kHz 16 bit shorts real modem signal from radio\n");
+        fprintf(stderr," -m RTTY|binary\n"); 
+        fprintf(stderr,"--mode=RTTY|binary     RTTY or binary Horus protcols\n");
+        fprintf(stderr,"--rate=[Rs]            Modem baud rate. Default: 100\n");
+        fprintf(stderr," -t[r] --stats=[r]     Print out modem statistics to stderr in JSON.\n");
+        fprintf(stderr,"                       r, if provided, sets the number of modem frames\n"
+                       "                       between statistic printouts\n");
+        fprintf(stderr," -q                    use stereo (IQ) input\n");
+        fprintf(stderr," -v                    verbose debug info\n");
+        fprintf(stderr," -c                    display CRC results for each packet\n");
+        exit(1);
+    }
+        
+    /* Open files */
+
+    if (verbose) {
+         fprintf(stderr, "mode: %d verbose: %d stats_loop: %d stats_rate: %d  Rs: %d\n",mode, verbose, stats_loop, stats_rate, Rs);
+    }
+    if (strcmp(argv[dx],"-")==0) {
+        fin = stdin;
+    } else {
+        fin = fopen(argv[dx],"rb");
+    }
+    
+    if (strcmp(argv[dx + 1],"-")==0) {
+        fout = stdout;
+    } else {
+        fout = fopen(argv[dx + 1],"w");
+    }
+        
+    if ((fin==NULL) || (fout==NULL)) {
+        fprintf(stderr,"Couldn't open test vector files\n");
+        exit(1);
+    }
+
+    /* end command line processing */
+
+    hstates = horus_open(mode, Rs);
+    horus_set_verbose(hstates, verbose);
+    
+    if (hstates == NULL) {
+        fprintf(stderr, "Couldn't open Horus API\n");
+        exit(1);
+    }
+    
+    if (enable_stats) {
+        loop_time = (float)horus_nin(hstates)/horus_get_Fs(hstates);
+        stats_loop = (int)(1.0/(stats_rate*loop_time));
+        stats_ctr = 0;
+    }
+
+    if((fsk_lower> 0) && (fsk_upper > fsk_lower)){
+        horus_set_freq_est_limits(hstates, fsk_lower, fsk_upper);
+        fprintf(stderr,"Setting estimator limits to %d to %d Hz.\n",fsk_lower, fsk_upper);
+    }
+
+    
+    int   max_demod_in = horus_get_max_demod_in(hstates);
+    int   hsize = quadrature ? 2 : 1;
+    short demod_in[max_demod_in * hsize];
+    int   max_ascii_out = horus_get_max_ascii_out_len(hstates);
+    char  ascii_out[max_ascii_out];
+    
+    /* Main loop ----------------------------------------------------------------------- */
+
+    while(fread(demod_in, hsize * sizeof(short), horus_nin(hstates), fin) ==  horus_nin(hstates)) {
+
+        if (verbose) {
+            fprintf(stderr, "read nin %d\n", horus_nin(hstates));
+        }
+        if (horus_rx(hstates, ascii_out, demod_in, quadrature)) {
+            fprintf(stdout, "%s", ascii_out);
+            if (crc_results) {
+                if (horus_crc_ok(hstates)) {
+                    fprintf(stdout, "  CRC OK");
+                } else {
+                    fprintf(stdout, "  CRC BAD");
+                }
+            }
+            fprintf(stdout, "\n");
+        }
+        
+        if (enable_stats && stats_ctr <= 0) {
+
+            horus_get_modem_extended_stats(hstates, &stats);
+
+	    /* Print standard 2FSK stats */
+
+            fprintf(stderr,"{\"EbNodB\": %2.2f,\t\"ppm\": %d,",stats.snr_est, (int)stats.clock_offset);
+            fprintf(stderr,"\t\"f1_est\":%.1f,\t\"f2_est\":%.1f",stats.f_est[0], stats.f_est[1]);
+
+	    /* Print 4FSK stats if in 4FSK mode */
+
+            if (horus_get_mFSK(hstates) == 4) {
+                fprintf(stderr,",\t\"f3_est\":%.1f,\t\"f4_est\":%.1f", stats.f_est[2], stats.f_est[3]);
+            }
+	    
+	    /* Print the eye diagram */
+
+            fprintf(stderr,",\t\"eye_diagram\":[");
+            for(i=0;i<stats.neyetr;i++){
+                fprintf(stderr,"[");
+                for(j=0;j<stats.neyesamp;j++){
+                    fprintf(stderr,"%f ",stats.rx_eye[i][j]);
+                    if(j<stats.neyesamp-1) fprintf(stderr,",");
+                }
+                fprintf(stderr,"]");
+                if(i<stats.neyetr-1) fprintf(stderr,",");
+            }
+            fprintf(stderr,"],");
+	    
+	    fprintf(stderr,"\"samp_fft\":[");
+
+            #ifdef FIXME_LATER
+            /* TODO: need a horus_ function to dig into modem spectrum */
+            
+	    /* Print a sample of the FFT from the freq estimator */
+
+	    Ndft = hstates->fsk->Ndft/2;
+	    for(i=0; i<Ndft; i++){
+		fprintf(stderr,"%f ",(hstates->fsk->fft_est)[i]);
+		if(i<Ndft-1) fprintf(stderr,",");
+	    }
+            #else
+
+            /* All zero dummy data for now */
+
+ 	    Ndft = 128;
+	    for(i=0; i<Ndft; i++) {
+		fprintf(stderr,"%f ", 0.0);
+		if(i<Ndft-1) fprintf(stderr,",");
+	    }
+            
+            #endif
+
+	    fprintf(stderr,"]}\n");
+            stats_ctr = stats_loop;
+        }
+        stats_ctr--;
+
+        if (fin == stdin || fout == stdout){
+            fflush(fin);
+            fflush(fout);
+        }
+    }
+    
+    horus_close(hstates);
+
+    return 0;
+}
+

src/horus_gen_test_bits.c

@@ -0,0 +1,195 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: horus_gen_tx_bits.c
+  AUTHOR......: Mark Jessop
+  DATE CREATED: May 2020
+
+  Horus dummy packet generation, for use with fsk_demod.
+
+  Build:
+  gcc horus_gen_test_bits.c horus_l2.c golay23.c -o horus_get_test_bits -Wall -DSCRAMBLER -DINTERLEAVER
+
+  \*---------------------------------------------------------------------------*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <getopt.h>
+
+#include "horus_l2.h"
+
+// TODO: Move these packet format definitions to somehwere common.
+
+/* Horus Mode 0 (Legacy 22-byte) Binary Packet */
+struct TBinaryPacket0
+{
+    uint8_t     PayloadID;
+    uint16_t	Counter;
+    uint8_t	Hours;
+    uint8_t	Minutes;
+    uint8_t	Seconds;
+    float	Latitude;
+    float	Longitude;
+    uint16_t  	Altitude;
+    uint8_t     Speed;       // Speed in Knots (1-255 knots)
+    uint8_t     Sats;
+    int8_t      Temp;        // Twos Complement Temp value.
+    uint8_t     BattVoltage; // 0 = 0.5v, 255 = 2.0V, linear steps in-between.
+    uint16_t    Checksum;    // CRC16-CCITT Checksum.
+}  __attribute__ ((packed));
+
+/* Horus Mode 1 (32-byte) Binary Packet */
+struct TBinaryPacket1
+{
+    uint16_t     PayloadID;
+    uint16_t	Counter;
+    uint8_t	Hours;
+    uint8_t	Minutes;
+    uint8_t	Seconds;
+    float	Latitude;
+    float	Longitude;
+    uint16_t  	Altitude;
+    uint8_t     Speed;       // Speed in Knots (1-255 knots)
+    uint8_t     Sats;
+    int8_t      Temp;        // Twos Complement Temp value.
+    uint8_t     BattVoltage; // 0 = 0.5v, 255 = 2.0V, linear steps in-between.
+    uint8_t     dummy1;      // Dummy values for user-configurable section.
+    uint8_t     dummy2;
+    uint8_t     dummy3;
+    uint8_t     dummy4;
+    uint8_t     dummy5;
+    uint8_t     dummy6;
+    uint8_t     dummy7;
+    uint8_t     dummy8;
+    uint8_t     dummy9;
+    uint16_t    Checksum;    // CRC16-CCITT Checksum.
+}  __attribute__ ((packed));
+
+/* Horus Mode 2 (16-byte) Binary Packet */
+struct TBinaryPacket2
+{
+    uint8_t     PayloadID;
+    uint8_t	Counter;
+    uint16_t	BiSeconds;
+    uint8_t	  LatitudeMSB;
+    uint16_t	Latitude;
+    uint8_t	  LongitudeMSB;
+    uint16_t	Longitude;
+    uint16_t  	Altitude;
+    uint8_t     BattVoltage; // 0 = 0.5v, 255 = 2.0V, linear steps in-between.
+    uint8_t     flags;      // Dummy values for user-configurable section.
+    uint16_t    Checksum;    // CRC16-CCITT Checksum.
+}  __attribute__ ((packed));
+
+
+
+
+int main(int argc,char *argv[]) {
+    int i, framecnt;
+    int horus_mode = 0;
+
+    char usage[] = "usage: %s horus_mode numFrames\nMode 0 = Legacy 22-byte Golay FEC\nMode 1 = 32-byte LDPC FEC\nMode 2 = 16-byte LDPC FEC\n";
+
+    if (argc < 3) {
+        fprintf(stderr, usage, argv[0]);
+        exit(1);
+    }
+
+    horus_mode = atoi(argv[1]);
+    fprintf(stderr, "Using Horus Mode %d.\n", horus_mode);
+
+    framecnt = atoi(argv[2]);
+    fprintf(stderr, "Generating %d frames.\n", framecnt);
+
+    if(horus_mode == 0){
+      int nbytes = sizeof(struct TBinaryPacket0);
+      struct TBinaryPacket0 input_payload;
+      int num_tx_data_bytes = horus_l2_get_num_tx_data_bytes(nbytes);
+      unsigned char tx[num_tx_data_bytes];
+
+      /* all zeros is nastiest sequence for demod before scrambling */
+
+      memset(&input_payload, 0, nbytes);
+      input_payload.Checksum = horus_l2_gen_crc16((unsigned char*)&input_payload, nbytes-2);
+
+      horus_l2_encode_tx_packet(tx, (unsigned char*)&input_payload, nbytes);
+
+      int b;
+      uint8_t tx_bit;
+      while(framecnt > 0){
+          for(i=0; i<num_tx_data_bytes; i++) {
+              for(b=0; b<8; b++) {
+                  tx_bit = (tx[i] >> (7-b)) & 0x1; /* msb first */
+                  fwrite(&tx_bit,sizeof(uint8_t),1,stdout);
+                  fflush(stdout);
+              }
+          }
+          framecnt -= 1;
+      }
+
+    } else if(horus_mode == 1){
+      // 32-Byte LDPC Encoded mode.
+      int nbytes = sizeof(struct TBinaryPacket1);
+      struct TBinaryPacket1 input_payload;
+
+      // TODO: Add Calculation of expected number of TX bytes based on LDPC code.
+      int num_tx_data_bytes = nbytes;
+      unsigned char tx[num_tx_data_bytes];
+
+      /* all zeros is nastiest sequence for demod before scrambling */
+      memset(&input_payload, 0, nbytes);
+      input_payload.Checksum = horus_l2_gen_crc16((unsigned char*)&input_payload, nbytes-2);
+
+
+      // TODO: Replaced with LDPC Encoding
+      memcpy(tx, (unsigned char*)&input_payload, nbytes);
+
+      int b;
+      uint8_t tx_bit;
+      while(framecnt > 0){
+          for(i=0; i<num_tx_data_bytes; i++) {
+              for(b=0; b<8; b++) {
+                  tx_bit = (tx[i] >> (7-b)) & 0x1; /* msb first */
+                  fwrite(&tx_bit,sizeof(uint8_t),1,stdout);
+                  fflush(stdout);
+              }
+          }
+          framecnt -= 1;
+      }
+    } else if(horus_mode == 2){
+      // 16-Byte LDPC Encoded mode.
+      int nbytes = sizeof(struct TBinaryPacket2);
+      struct TBinaryPacket2 input_payload;
+
+      // TODO: Add Calculation of expected number of TX bytes based on LDPC code.
+      int num_tx_data_bytes = nbytes;
+      unsigned char tx[num_tx_data_bytes];
+
+      /* all zeros is nastiest sequence for demod before scrambling */
+      memset(&input_payload, 0, nbytes);
+      input_payload.Checksum = horus_l2_gen_crc16((unsigned char*)&input_payload, nbytes-2);
+
+
+      // TODO: Replaced with LDPC Encoding
+      memcpy(tx, (unsigned char*)&input_payload, nbytes);
+
+      int b;
+      uint8_t tx_bit;
+      while(framecnt > 0){
+          for(i=0; i<num_tx_data_bytes; i++) {
+              for(b=0; b<8; b++) {
+                  tx_bit = (tx[i] >> (7-b)) & 0x1; /* msb first */
+                  fwrite(&tx_bit,sizeof(uint8_t),1,stdout);
+                  fflush(stdout);
+              }
+          }
+          framecnt -= 1;
+      }
+    } else {
+      fprintf(stderr, "Unknown Mode!");
+    }
+
+    return 0;
+}

src/horus_l2.c

@@ -0,0 +1,946 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: horus_l2.c
+  AUTHOR......: David Rowe
+  DATE CREATED: Dec 2015
+
+  Horus telemetry layer 2 processing.  Takes an array of 8 bit payload
+  data, generates parity bits for a (23,12) Golay code, interleaves
+  data and parity bits, pre-pends a Unique Word for modem sync.
+  Caller is responsible for providing storage for output packet.
+
+  1/ Unit test on a PC:
+
+     $ gcc horus_l2.c -o horus_l2 -Wall -DHORUS_L2_UNITTEST
+     $ ./horus_l2
+
+     test 0: 22 bytes of payload data BER: 0.00 errors: 0
+     test 0: 22 bytes of payload data BER: 0.01 errors: 0
+     test 0: 22 bytes of payload data BER: 0.05 errors: 0
+     test 0: 22 bytes of payload data BER: 0.10 errors: 7
+     
+     This indicates it's correcting all channel errors for 22 bytes of
+     payload data, at bit error rate (BER) of 0, 0.01, 0.05.  It falls
+     over at a BER of 0.10 which is expected.
+
+  2/ To build with just the tx function, ie for linking with the payload
+  firmware:
+
+    $ gcc horus_l2.c golay23.c -c -Wall
+    
+  By default the RX side is #ifdef-ed out, leaving the minimal amount
+  of code for tx.
+
+  3/ Generate some tx_bits as input for testing with fsk_horus:
+ 
+    $ gcc horus_l2.c golay23.c -o horus_l2 -Wall -DGEN_TX_BITS -DSCRAMBLER -DINTERLEAVER
+    $ ./horus_l2
+    $ more ../octave/horus_tx_bits_binary.txt
+
+  4/ Streaming test bits to stdout, for 'live' testing with fsk_mod and horus_demod:
+
+    $ gcc horus_l2.c golay23.c -o horus_l2 -Wall -DGEN_TX_BITSTREAM -DSCRAMBLER -DINTERLEAVER
+    $ cp horus_l2 ../build/src/
+    $ cd ../build/src/
+    $ ./horus_l2 100 | ./fsk_mod 4 48000 100 750 250 - - | ./horus_demod -m binary - -
+   
+  5/ Unit testing interleaver:
+
+    $ gcc horus_l2.c golay23.c -o horus_l2 -Wall -DINTERLEAVER -DTEST_INTERLEAVER -DSCRAMBLER
+
+  6/ Compile for use as decoder called by fsk_horus.m and fsk_horus_stream.m:
+
+    $ gcc horus_l2.c golay23.c -o horus_l2 -Wall -DDEC_RX_BITS -DHORUS_L2_RX
+
+\*---------------------------------------------------------------------------*/
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdint.h>
+#include "horus_l2.h"
+#include "golay23.h"
+
+#ifdef HORUS_L2_UNITTEST
+#define HORUS_L2_RX
+#endif
+
+static char uw[] = {'$','$'};
+
+/* Function Prototypes ------------------------------------------------*/
+
+#ifdef INTERLEAVER
+static void interleave(unsigned char *inout, int nbytes, int dir);
+#endif
+#ifdef SCRAMBLER
+static void scramble(unsigned char *inout, int nbytes);
+#endif
+
+/* Functions ----------------------------------------------------------*/
+
+/*
+   We are using a Golay (23,12) code which has a codeword 23 bits
+   long.  The tx packet format is:
+
+      | Unique Word | payload data bits | parity bits |
+
+   This function works out how much storage the caller of
+   horus_l2_encode_tx_packet() will need to store the tx packet
+ */
+
+int horus_l2_get_num_tx_data_bytes(int num_payload_data_bytes) {
+    int num_payload_data_bits, num_golay_codewords;
+    int num_tx_data_bits, num_tx_data_bytes;
+    
+    num_payload_data_bits = num_payload_data_bytes*8;
+    num_golay_codewords = num_payload_data_bits/12;
+    if (num_payload_data_bits % 12) /* round up to 12 bits, may mean some unused bits */
+        num_golay_codewords++;
+
+    num_tx_data_bits = sizeof(uw)*8 + num_payload_data_bits + num_golay_codewords*11;
+    num_tx_data_bytes = num_tx_data_bits/8;
+    if (num_tx_data_bits % 8) /* round up to nearest byte, may mean some unused bits */
+        num_tx_data_bytes++;
+    
+    #ifdef DEBUG0
+    fprintf(stderr, "\nnum_payload_data_bytes: %d\n", num_payload_data_bytes);
+    fprintf(stderr, "num_golay_codewords...: %d\n", num_golay_codewords);
+    fprintf(stderr, "num_tx_data_bits......: %d\n", num_tx_data_bits);
+    fprintf(stderr, "num_tx_data_bytes.....: %d\n\n", num_tx_data_bytes);
+    #endif
+
+    return num_tx_data_bytes;
+}
+
+void horus_l2_init(void) {
+    golay23_init();
+}
+
+/*
+  Takes an array of payload data bytes, prepends a unique word and appends
+  parity bits.
+
+  The encoder will run on the payload on a small 8-bit uC.  As we are
+  memory constrained so we do a lot of burrowing for bits out of
+  packed arrays, and don't use a LUT for Golay encoding.  Hopefully it
+  will run fast enough.  This was quite difficult to get going,
+  suspect there is a better way to write this.  Oh well, have to start
+  somewhere.
+ */
+
+int horus_l2_encode_tx_packet(unsigned char *output_tx_data,
+                              unsigned char *input_payload_data,
+                              int            num_payload_data_bytes)
+{
+    int            num_tx_data_bytes, num_payload_data_bits;
+    unsigned char *pout = output_tx_data;
+    int            ninbit, ningolay, nparitybits;
+    int32_t        ingolay, paritybyte, inbit, golayparity;
+    int            ninbyte, shift, golayparitybit, i;
+
+    num_tx_data_bytes = horus_l2_get_num_tx_data_bytes(num_payload_data_bytes);
+    memcpy(pout, uw, sizeof(uw)); pout += sizeof(uw);
+    memcpy(pout, input_payload_data, num_payload_data_bytes); pout += num_payload_data_bytes;
+
+    /* Read input bits one at a time.  Fill input Golay codeword.  Find output Golay codeword.
+       Write this to parity bits.  Write parity bytes when we have 8 parity bits.  Bits are
+       written MSB first. */
+
+    num_payload_data_bits = num_payload_data_bytes*8;
+    ninbit = 0;
+    ingolay = 0;
+    ningolay = 0;
+    paritybyte = 0;
+    nparitybits = 0;
+
+    while (ninbit < num_payload_data_bits) {
+
+        /* extract input data bit */
+
+        ninbyte = ninbit/8;
+        shift = 7 - (ninbit % 8);
+        inbit = (input_payload_data[ninbyte] >> shift) & 0x1;
+        #ifdef DEBUG1
+        fprintf(stderr, "inbit %d ninbyte: %d inbyte: 0x%02x inbit: %d\n", 
+                ninbit, ninbyte, input_payload_data[ninbyte], inbit);
+        #endif
+        ninbit++;
+
+        /* build up input golay codeword */
+
+        ingolay = ingolay | inbit;
+        ningolay++;
+
+        /* when we get 12 bits do a Golay encode */
+
+        if (ningolay % 12) {
+            ingolay <<= 1;
+        }
+        else {
+            #ifdef DEBUG0
+            fprintf(stderr, "  ningolay: %d ingolay: 0x%04x\n", ningolay, ingolay);
+            #endif
+            golayparity = golay23_syndrome(ingolay<<11);
+            ingolay = 0;
+
+            #ifdef DEBUG0
+            fprintf(stderr, "  golayparity: 0x%04x\n", golayparity);
+            #endif
+
+            /* write parity bits to output data */
+
+            for (i=0; i<11; i++) {
+                golayparitybit = (golayparity >> (10-i)) & 0x1;
+                paritybyte = paritybyte | golayparitybit;
+                #ifdef DEBUG0
+                fprintf(stderr, "    i: %d golayparitybit: %d paritybyte: 0x%02x\n", 
+                        i, golayparitybit, paritybyte);
+                #endif
+                nparitybits++;
+                if (nparitybits % 8) {
+                   paritybyte <<= 1;
+                }
+                else {
+                    /* OK we have a full byte ready */
+                    *pout = paritybyte;
+                    #ifdef DEBUG0
+                    fprintf(stderr,"      Write paritybyte: 0x%02x\n", paritybyte);
+                    #endif
+                    pout++;
+                    paritybyte = 0;
+                }
+            }
+        }
+    } /* while(.... */
+
+
+    /* Complete final Golay encode, we may have partially finished ingolay, paritybyte */
+
+    #ifdef DEBUG0
+    fprintf(stderr, "finishing up .....\n");
+    #endif
+
+    if (ningolay % 12) {
+        ingolay >>= 1;
+        golayparity = golay23_syndrome(ingolay<<12);
+        #ifdef DEBUG0
+        fprintf(stderr, "  ningolay: %d ingolay: 0x%04x\n", ningolay, ingolay);
+        fprintf(stderr, "  golayparity: 0x%04x\n", golayparity);
+        #endif
+
+        /* write parity bits to output data */
+
+        for (i=0; i<11; i++) {
+            golayparitybit = (golayparity >> (10 - i)) & 0x1;
+            paritybyte = paritybyte | golayparitybit;
+            #ifdef DEBUG1
+            fprintf(stderr, "    i: %d golayparitybit: %d paritybyte: 0x%02x\n", 
+                    i, golayparitybit, paritybyte);
+            #endif
+            nparitybits++;
+            if (nparitybits % 8) {
+                paritybyte <<= 1;
+            }
+            else {
+                /* OK we have a full byte ready */
+                *pout++ = (unsigned char)paritybyte;
+                #ifdef DEBUG0
+                fprintf(stderr,"      Write paritybyte: 0x%02x\n", paritybyte);
+                #endif
+                paritybyte = 0;
+            }
+        }
+    }
+ 
+    /* and final, partially complete, parity byte */
+
+    if (nparitybits % 8) {
+        paritybyte <<= 7 - (nparitybits % 8);  // use MS bits first
+        *pout++ = (unsigned char)paritybyte;
+        #ifdef DEBUG0
+        fprintf(stderr,"      Write last paritybyte: 0x%02x nparitybits: %d \n", paritybyte, nparitybits);
+        #endif
+    }
+
+    #ifdef DEBUG0
+    fprintf(stderr, "\npout - output_tx_data: %ld num_tx_data_bytes: %d\n",
+            pout - output_tx_data, num_tx_data_bytes);
+    #endif
+    assert(pout == (output_tx_data + num_tx_data_bytes));
+
+    /* optional interleaver - we dont interleave UW */
+
+    #ifdef INTERLEAVER
+    interleave(&output_tx_data[sizeof(uw)], num_tx_data_bytes-2, 0);
+    #endif
+
+    /* optional scrambler to prevent long strings of the same symbol
+       which upsets the modem - we dont scramble UW */
+
+    #ifdef SCRAMBLER
+    scramble(&output_tx_data[sizeof(uw)], num_tx_data_bytes-2);
+    #endif
+
+    return num_tx_data_bytes;
+}
+
+
+#ifdef HORUS_L2_RX
+void horus_l2_decode_rx_packet(unsigned char *output_payload_data,
+                               unsigned char *input_rx_data,
+                               int            num_payload_data_bytes)
+{
+    int            num_payload_data_bits;
+    unsigned char *pout = output_payload_data;
+    unsigned char *pin  = input_rx_data;
+    int            ninbit, ingolay, ningolay, paritybyte, nparitybits;
+    int            ninbyte, shift, inbit, golayparitybit, i, outbit, outbyte, noutbits, outdata;
+    #if defined(SCRAMBLER) || defined(INTERLEAVER)
+    int num_tx_data_bytes = horus_l2_get_num_tx_data_bytes(num_payload_data_bytes);
+    #endif
+    
+    /* optional scrambler and interleaver - we dont interleave UW */
+
+    #ifdef SCRAMBLER
+    scramble(&input_rx_data[sizeof(uw)], num_tx_data_bytes-2);
+    #endif
+
+    #ifdef INTERLEAVER
+    interleave(&input_rx_data[sizeof(uw)], num_tx_data_bytes-2, 1);
+    #endif
+
+    pin = input_rx_data + sizeof(uw) + num_payload_data_bytes;
+
+    /* Read input data bits one at a time.  When we have 12 read 11 parity bits. Golay decode.
+       Write decoded (output data) bits every time we have 8 of them. */
+
+    num_payload_data_bits = num_payload_data_bytes*8;
+    ninbit = 0;
+    ingolay = 0;
+    ningolay = 0;
+    nparitybits = 0;
+    paritybyte = *pin++;
+    #ifdef DEBUG0
+    fprintf(stderr,"  Read paritybyte: 0x%02x\n", paritybyte);
+    #endif
+    pout = output_payload_data;
+    noutbits = 0;
+    outbyte = 0;
+
+    while (ninbit < num_payload_data_bits) {
+
+        /* extract input data bit */
+
+        ninbyte = ninbit/8 + sizeof(uw);
+        shift = 7 - (ninbit % 8);
+        inbit = (input_rx_data[ninbyte] >> shift) & 0x1;
+        #ifdef DEBUG1
+        fprintf(stderr, "inbit %d ninbyte: %d inbyte: 0x%02x inbit: %d\n", 
+                ninbit, ninbyte, input_rx_data[ninbyte], inbit);
+        #endif
+        ninbit++;
+
+        /* build up golay codeword */
+
+        ingolay = ingolay | inbit;
+        ningolay++;
+        ingolay <<= 1;
+
+        /* when we get 12 data bits start reading parity bits */
+
+        if ((ningolay % 12) == 0) {
+            #ifdef DEBUG0
+            fprintf(stderr, "  ningolay: %d ingolay: 0x%04x\n", ningolay, ingolay>>1);
+            #endif
+            for (i=0; i<11; i++) {
+                shift = 7 - (nparitybits % 8);
+                golayparitybit = (paritybyte >> shift) & 0x1;
+                ingolay |= golayparitybit;
+                if (i != 10)
+                    ingolay <<=1;
+                nparitybits++;
+                if ((nparitybits % 8) == 0) {
+                    /* OK grab a new byte */
+                    paritybyte = *pin++;
+                    #ifdef DEBUG0
+                    fprintf(stderr,"  Read paritybyte: 0x%02x\n", paritybyte);
+                    #endif
+                }
+            }
+
+            #ifdef DEBUG0
+            fprintf(stderr, "  golay code word: 0x%04x\n", ingolay);
+            fprintf(stderr, "  golay decode...: 0x%04x\n", golay23_decode(ingolay));
+            #endif
+           
+            /* write decoded/error corrected bits to output payload data */
+
+            outdata = golay23_decode(ingolay) >> 11;
+            #ifdef DEBUG0
+            fprintf(stderr, "  outdata...: 0x%04x\n", outdata);
+            #endif
+
+            for(i=0; i<12; i++) {   
+                shift = 11 - i;
+                outbit = (outdata >> shift) & 0x1;
+                outbyte |= outbit;
+                noutbits++;
+                if (noutbits % 8) {
+                    outbyte <<= 1;
+                }
+                else {
+                    #ifdef DEBUG0
+                    fprintf(stderr, "  output payload byte: 0x%02x\n", outbyte);
+                    #endif
+                    *pout++ = outbyte;
+                    outbyte = 0;
+                }
+            }
+
+            ingolay = 0;
+        }
+    } /* while(.... */
+
+
+    #ifdef DEBUG0
+    fprintf(stderr, "finishing up .....\n");
+    #endif
+
+    /* Complete final Golay decode  */
+
+    int golayparity = 0;
+    if (ningolay % 12) {
+        for (i=0; i<11; i++) {
+            shift = 7 - (nparitybits % 8);
+            golayparitybit = (paritybyte >> shift) & 0x1;
+            golayparity |= golayparitybit;
+            if (i != 10)
+                golayparity <<=1;
+            nparitybits++;
+            if ((nparitybits % 8) == 0) {
+                /* OK grab a new byte */
+                paritybyte = *pin++;
+                #ifdef DEBUG0
+                fprintf(stderr,"  Read paritybyte: 0x%02x\n", paritybyte);
+                #endif
+            }
+        }
+
+        ingolay >>= 1;
+        int codeword = (ingolay<<12) + golayparity;
+        #ifdef DEBUG0
+        fprintf(stderr, "  ningolay: %d ingolay: 0x%04x\n", ningolay, ingolay);
+        fprintf(stderr, "  golay code word: 0x%04x\n", codeword);
+        fprintf(stderr, "  golay decode...: 0x%04x\n", golay23_decode(codeword));
+        #endif
+
+        outdata = golay23_decode(codeword) >> 11;
+        #ifdef DEBUG0
+        fprintf(stderr, "  outdata...: 0x%04x\n", outdata);
+        fprintf(stderr, "  num_payload_data_bits: %d noutbits: %d\n", num_payload_data_bits, noutbits);
+        #endif
+
+        /* write final byte */
+
+        int ntogo = num_payload_data_bits - noutbits;
+        for(i=0; i<ntogo; i++) {   
+            shift = ntogo - i;
+            outbit = (outdata >> shift) & 0x1;
+            outbyte |= outbit;
+            noutbits++;
+            if (noutbits % 8) {
+                outbyte <<= 1;
+            }
+            else {
+                #ifdef DEBUG0
+                fprintf(stderr, "  output payload byte: 0x%02x\n", outbyte);
+                #endif
+                *pout++ = outbyte;
+                outbyte = 0;
+            }
+        }
+    }
+
+    #ifdef DEBUG0
+    fprintf(stderr, "\npin - output_payload_data: %ld num_payload_data_bytes: %d\n",
+            pout - output_payload_data, num_payload_data_bytes);
+    #endif
+
+    assert(pout == (output_payload_data + num_payload_data_bytes));
+
+}
+#endif
+
+#ifdef INTERLEAVER
+
+uint16_t primes[] = {
+    2,      3,      5,      7,      11,     13,     17,     19,     23,     29, 
+    31,     37,     41,     43,     47,     53,     59,     61,     67,     71, 
+    73,     79,     83,     89,     97,     101,    103,    107,    109,    113, 
+    127,    131,    137,    139,    149,    151,    157,    163,    167,    173, 
+    179,    181,    191,    193,    197,    199,    211,    223,    227,    229, 
+    233,    239,    241,    251,    257,    263,    269,    271,    277,    281, 
+    283,    293,    307,    311,    313,    317,    331,    337,    347
+};
+
+void interleave(unsigned char *inout, int nbytes, int dir)
+{
+    /* note: to work on small uCs (e.g. AVR) needed to declare specific words sizes */
+    uint16_t nbits = (uint16_t)nbytes*8;
+    uint32_t i, j, n, ibit, ibyte, ishift, jbyte, jshift;
+    uint32_t b;
+    unsigned char out[nbytes];
+
+    memset(out, 0, nbytes);
+           
+    /* b chosen to be co-prime with nbits, I'm cheating by just finding the 
+       nearest prime to nbits.  It also uses storage, is run on every call,
+       and has an upper limit.  Oh Well, still seems to interleave OK. */
+    i = 1;
+    uint16_t imax = sizeof(primes)/sizeof(uint16_t);
+    while ((primes[i] < nbits) && (i < imax))
+        i++;
+    b = primes[i-1];
+
+    for(n=0; n<nbits; n++) {
+
+        /*
+          "On the Analysis and Design of Good Algebraic Interleavers", Xie et al,eq (5)
+        */
+
+        i = n;
+        j = (b*i) % nbits; /* note these all need to be 32-bit ints to make multiply work without overflow */
+        
+        if (dir) {
+            uint16_t tmp = j;
+            j = i;
+            i = tmp;
+        }
+        
+        #ifdef DEBUG0
+        printf("i: %d j: %d\n",i, j);
+        #endif
+
+        /* read bit i and write to bit j postion */
+
+        ibyte = i/8;
+        ishift = i%8;
+        ibit = (inout[ibyte] >> ishift) & 0x1;
+
+        jbyte = j/8;
+        jshift = j%8;
+
+        /* write jbit to ibit position */
+
+        out[jbyte] |= ibit << jshift; // replace with i-th bit
+        //out[ibyte] |= ibit << ishift; // replace with i-th bit
+    }
+ 
+    memcpy(inout, out, nbytes);
+
+    #ifdef DEBUG0
+    printf("\nInterleaver Out:\n");
+    for (i=0; i<nbytes; i++)
+        printf("%02d 0x%02x\n", i, inout[i]);
+    #endif
+}
+#endif
+
+
+#ifdef TEST_INTERLEAVER
+int main(void) {
+    int nbytes = 43;
+    unsigned char inout[nbytes];
+    unsigned char inter[nbytes];
+    unsigned char incopy[nbytes];
+    int i;
+
+    /* copy of input for later comp   */
+
+    for(i=0; i<nbytes; i++)
+        inout[i] = incopy[i] = rand() & 0xff;    
+    
+    interleave(inout, nbytes, 0);    /* interleave                     */
+    memcpy(inter, inout, nbytes);    /* snap shot of interleaved bytes */
+    interleave(inout, nbytes, 1);    /* de-interleave                  */
+
+    /* all ones in last col means it worked! */
+
+    for(i=0; i<nbytes; i++) {
+        printf("%d 0x%02x 0x%02x 0x%02x %d\n", 
+               i, incopy[i], inter[i], inout[i],  incopy[i] == inout[i]);
+        assert(incopy[i] == inout[i]);
+    }
+    printf("Interleaver tested OK!\n");
+
+    return 0;
+}
+#endif
+
+
+#ifdef SCRAMBLER
+
+/* 16 bit DVB additive scrambler as per Wikpedia example */
+
+void scramble(unsigned char *inout, int nbytes)
+{
+    int nbits = nbytes*8;
+    int i, ibit, ibits, ibyte, ishift, mask;
+    uint16_t scrambler = 0x4a80;  /* init additive scrambler at start of every frame */
+    uint16_t scrambler_out;
+
+    /* in place modification of each bit */
+
+    for(i=0; i<nbits; i++) {
+
+        scrambler_out = ((scrambler & 0x2) >> 1) ^ (scrambler & 0x1);
+
+        /* modify i-th bit by xor-ing with scrambler output sequence */
+
+        ibyte = i/8;
+        ishift = i%8;
+        ibit = (inout[ibyte] >> ishift) & 0x1;
+        ibits = ibit ^ scrambler_out;                  // xor ibit with scrambler output
+
+        mask = 1 << ishift;
+        inout[ibyte] &= ~mask;                  // clear i-th bit
+        inout[ibyte] |= ibits << ishift;         // set to scrambled value
+
+        /* update scrambler */
+
+        scrambler >>= 1;
+        scrambler |= scrambler_out << 14;
+
+        #ifdef DEBUG0
+        printf("i: %02d ibyte: %d ishift: %d ibit: %d ibits: %d scrambler_out: %d\n", 
+               i, ibyte, ishift, ibit, ibits, scrambler_out);
+        #endif
+
+    }
+
+    #ifdef DEBUG0
+    printf("\nScrambler Out:\n");
+    for (i=0; i<nbytes; i++)
+        printf("%02d 0x%02x\n", i, inout[i]);
+    #endif
+}
+#endif
+
+#ifdef HORUS_L2_UNITTEST
+
+/*
+  Test function to construct a packet of payload data, encode, add
+  some bit errors, decode, count errors.
+*/
+
+int test_sending_bytes(int nbytes, float ber, int error_pattern) {
+    unsigned char input_payload[nbytes];
+    int num_tx_data_bytes = horus_l2_get_num_tx_data_bytes(sizeof(input_payload));
+    unsigned char tx[num_tx_data_bytes];
+    unsigned char output_payload[sizeof(input_payload)];
+    int b, nbiterrors = 0;
+    int i;
+
+    for(i=0; i<nbytes; i++)
+        input_payload[i] = i;
+
+    horus_l2_encode_tx_packet(tx, input_payload, sizeof(input_payload));
+
+    #ifdef DEBUG0
+    fprintf(stderr, "\nTx Data:\n");
+    for(i=0; i<num_tx_data_bytes; i++)
+        fprintf(stderr, "  %02d 0x%02x\n", i, tx[i]);
+    #endif
+
+    /* insert random bit errors */
+
+    if (error_pattern == 0) {
+        float r;
+        for(i=0; i<num_tx_data_bytes; i++) {
+            for (b=0; b<8; b++) {
+                r = (float)rand()/RAND_MAX;
+                if (r < ber) {
+                    unsigned char mask = (1<<b);
+                    #ifdef DEBUG1
+                    fprintf("mask: 0x%x tx[%d] = 0x%x ", mask, i, tx[i]);
+                    #endif
+                    tx[i] ^= mask;
+                    #ifdef DEBUG1
+                    fprintf("0x%x\n", tx[i]);
+                    #endif
+                    nbiterrors++;
+                }
+            }
+        }
+    }
+
+    /* insert and error burst */
+
+    if (error_pattern == 1) {
+        tx[2] ^= 0xff;
+        tx[3] ^= 0xff;
+    }
+
+    /* insert 1 error every 12 bits, this gives up to 3 errors per 23
+       bit codeword, which is the limit of the code */
+
+    if (error_pattern == 2) {
+        int bn = 0;
+        for(i=0; i<num_tx_data_bytes; i++) {
+            for (b=0; b<8; b++) {
+                bn++;
+                if ((bn % 12) == 0) {
+                    unsigned char mask = (1<<(7-b));
+                    #ifdef DEBUG1
+                    fprintf("mask: 0x%x tx[%d] = 0x%x ", mask, i, tx[i]);
+                    #endif
+                    tx[i] ^= mask;
+                    #ifdef DEBUG1
+                    fprintf("0x%x\n", tx[i]);
+                    #endif
+                    nbiterrors++;
+                }
+            }
+        }
+    }
+
+    #ifdef DEBUG0
+    fprintf(stderr, "\nTx Data after errors:\n");
+    for(i=0; i<num_tx_data_bytes; i++)
+        fprintf(stderr, "  %02d 0x%02x\n", i, tx[i]);
+    #endif
+
+    #ifdef DEBUG0
+    fprintf(stderr, "nbiterrors: %d BER: %3.2f\n", nbiterrors, (float)nbiterrors/(num_tx_data_bytes*8));
+    #endif
+
+    golay23_init();
+    horus_l2_decode_rx_packet(output_payload, tx, sizeof(input_payload));
+
+    #ifdef DEBUG0
+    fprintf(stderr, "\nOutput Payload:\n");
+    for(i=0; i<sizeof(input_payload); i++)
+        fprintf(stderr, "  %02d 0x%02x\n", i, output_payload[i]);
+    #endif
+
+    /* count bit errors */
+
+    int nerr = 0;
+    for(i=0; i<nbytes; i++) {
+        int error_pattern = input_payload[i] ^ output_payload[i];
+        for(b=0; b<8; b++)
+            nerr += (error_pattern>>b) & 0x1;
+    }
+    
+    return nerr;
+}
+
+/* unit test designed to run on a PC */
+
+int main(void) {
+    printf("test 0: BER: 0.00 ...........: %d\n", test_sending_bytes(22, 0.00, 0));
+    printf("test 1: BER: 0.01 ...........: %d\n", test_sending_bytes(22, 0.01, 0));
+    printf("test 2: BER: 0.05 ...........: %d\n", test_sending_bytes(22, 0.05, 0));
+
+    /* we expect this always to fail, as chance of > 3 errors/codeword is high */
+
+    printf("test 3: BER: 0.10 ...........: %d\n", test_sending_bytes(22, 0.10, 0));
+
+    /* -DINTERLEAVER will make this puppy pass */
+
+    printf("test 4: 8 bit burst error....: %d\n", test_sending_bytes(22, 0.00, 1));
+
+    /* Insert 2 errors in every codeword, the maximum correction
+       capability of a Golay (23,12) code. note this one will fail
+       with -DINTERLEAVER, as we can't guarantee <= 3 errors per
+       codeword after interleaving */
+
+    printf("test 5: 1 error every 12 bits: %d\n", test_sending_bytes(22, 0.00, 2));
+    return 0;
+}
+#endif
+
+/* Horus binary packet */
+
+struct TBinaryPacket
+{
+    uint8_t     PayloadID;
+    uint16_t	Counter;
+    uint8_t	Hours;
+    uint8_t	Minutes;
+    uint8_t	Seconds;
+    float	Latitude;
+    float	Longitude;
+    uint16_t  	Altitude;
+    uint8_t     Speed;       // Speed in Knots (1-255 knots)
+    uint8_t     Sats;
+    int8_t      Temp;        // Twos Complement Temp value.
+    uint8_t     BattVoltage; // 0 = 0.5v, 255 = 2.0V, linear steps in-between.
+    uint16_t    Checksum;    // CRC16-CCITT Checksum.
+}  __attribute__ ((packed));
+
+#ifdef GEN_TX_BITS
+/* generate a file of tx_bits to modulate using fsk_horus.m for modem simulations */
+
+int main(void) {
+    int nbytes = sizeof(struct TBinaryPacket);
+    struct TBinaryPacket input_payload;
+    int num_tx_data_bytes = horus_l2_get_num_tx_data_bytes(nbytes);
+    unsigned char tx[num_tx_data_bytes];
+    int i;
+
+    /* all zeros is nastiest sequence for demod before scrambling */
+
+    memset(&input_payload, 0, nbytes);
+    input_payload.Checksum = horus_l2_gen_crc16((unsigned char*)&input_payload, nbytes-2);
+
+    horus_l2_encode_tx_packet(tx, (unsigned char*)&input_payload, nbytes);
+    
+    FILE *f = fopen("../octave/horus_tx_bits_binary.txt","wt");
+    assert(f != NULL);
+    int b, tx_bit;
+    for(i=0; i<num_tx_data_bytes; i++) {
+        for(b=0; b<8; b++) {
+            tx_bit = (tx[i] >> (7-b)) & 0x1; /* msb first */
+            fprintf(f,"%d ", tx_bit);
+        }
+    }
+    fclose(f);
+    for(i=0; i<num_tx_data_bytes; i++) {
+        fprintf(stdout,"%02X", tx[i]);
+    }
+    fprintf(stdout, "\n");
+
+    return 0;
+}
+#endif
+
+
+#ifdef DEC_RX_BITS
+
+/* Decode a binary file rx_bytes, e.g. from fsk_horus.m */
+
+int main(void) {
+    int nbytes = 22;
+    unsigned char output_payload[nbytes];
+    int num_tx_data_bytes = horus_l2_get_num_tx_data_bytes(nbytes);
+
+    /* real world data horus payload generated when running tx above */
+    unsigned char rx[45] = {
+        0x24,0x24,0x01,0x0b,0x00,0x00,0x05,0x3b,0xf2,0xa7,0x0b,0xc2,0x1b,
+        0xaa,0x0a,0x43,0x7e,0x00,0x05,0x00,0x25,0xc0,0xce,0xbb,0x36,0x69,
+        0x50,0x00,0x41,0xb0,0xa6,0x5e,0x91,0xa2,0xa3,0xf8,0x1d,0x00,0x00,
+        0x0c,0x76,0xc6,0x05,0xb0,0xb8};
+    int i, ret;
+
+    assert(num_tx_data_bytes == 45);
+
+    #define READ_FILE /* overwrite tx[] above, that's OK */
+    #ifdef READ_FILE
+    FILE *f = fopen("../octave/horus_rx_bits_binary.bin","rb");
+    assert(f != NULL);
+    ret = fread(rx, sizeof(char), num_tx_data_bytes, f);
+    assert(ret == num_tx_data_bytes);
+    fclose(f);
+    #endif
+
+    golay23_init();
+    horus_l2_decode_rx_packet(output_payload, rx, nbytes);
+
+    #ifdef HEX_DUMP
+    fprintf(stderr, "\nOutput Payload:\n");
+    for(i=0; i<nbytes; i++)
+        fprintf(stderr, "  %02d 0x%02x 0x%02x\n", i, output_payload[i], rx[i+2]);
+    #endif
+
+    struct TBinaryPacket h;
+    assert(sizeof(h) == nbytes);
+    memcpy(&h, output_payload, nbytes);
+
+    uint16_t crc_rx = horus_l2_gen_crc16(output_payload, nbytes-2);
+    char crc_str[80];
+    
+    if (crc_rx == h.Checksum) {
+        sprintf(crc_str, "CRC OK");
+    } else {
+        sprintf(crc_str, "CRC BAD");
+    }
+
+    fprintf(stderr, "%d,%d,%02d:%02d:%02d,%f,%f,%d,%d,%d,%d,%d,%04x %s\n",
+        h.PayloadID, h.Counter, h.Hours, h.Minutes, h.Seconds,
+        h.Latitude, h.Longitude, h.Altitude, h.Speed, h.Sats, h.Temp, 
+            h.BattVoltage, h.Checksum, crc_str);
+    
+    /* Hex ASCII file output */
+
+    #define WRITE_HEX_FILE /* overwrite tx[] above, that's OK */
+    #ifdef WRITE_HEX_FILE
+    FILE *fh = fopen("../octave/horus_rx_bits_hex.txt","wt");
+    assert(fh != NULL);
+    for(i=0; i<nbytes; i++) {
+        fprintf(fh, "%02X", (unsigned int)output_payload[i]);
+    }
+    fclose(fh);
+    #endif
+
+    return 0;
+}
+#endif
+
+
+#ifdef GEN_TX_BITSTREAM
+/* Generate a stream of encoded Horus packets in a format suitable to feed into fsk_mod */
+
+int main(int argc,char *argv[]) {
+    int nbytes = sizeof(struct TBinaryPacket);
+    struct TBinaryPacket input_payload;
+    int num_tx_data_bytes = horus_l2_get_num_tx_data_bytes(nbytes);
+    unsigned char tx[num_tx_data_bytes];
+    int i, framecnt;
+
+    if(argc != 2){
+        fprintf(stderr,"usage: %s numFrames\n",argv[0]);
+        exit(1);
+    }
+
+    framecnt = atoi(argv[1]);
+
+    /* all zeros is nastiest sequence for demod before scrambling */
+
+    memset(&input_payload, 0, nbytes);
+    input_payload.Checksum = horus_l2_gen_crc16((unsigned char*)&input_payload, nbytes-2);
+
+    horus_l2_encode_tx_packet(tx, (unsigned char*)&input_payload, nbytes);
+
+    int b;
+    uint8_t tx_bit;
+    while(framecnt >= 0){
+        for(i=0; i<num_tx_data_bytes; i++) {
+            for(b=0; b<8; b++) {
+                tx_bit = (tx[i] >> (7-b)) & 0x1; /* msb first */
+                fwrite(&tx_bit,sizeof(uint8_t),1,stdout);
+                fflush(stdout);
+            }
+        }
+        framecnt -= 1;
+    }
+
+    return 0;
+}
+#endif
+
+// from http://stackoverflow.com/questions/10564491/function-to-calculate-a-crc16-checksum
+
+unsigned short horus_l2_gen_crc16(unsigned char* data_p, unsigned char length) {
+    unsigned char x;
+    unsigned short crc = 0xFFFF;
+
+    while (length--){
+        x = crc >> 8 ^ *data_p++;
+        x ^= x>>4;
+        crc = (crc << 8) ^ ((unsigned short)(x << 12)) ^ ((unsigned short)(x <<5)) ^ ((unsigned short)x);
+    }
+    return crc;
+}
+

src/horus_l2.h

@@ -0,0 +1,30 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: horus_l2.h
+  AUTHOR......: David Rowe
+  DATE CREATED: Dec 2015
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef __HORUS_L2__
+#define __HORUS_L2__
+
+int horus_l2_get_num_tx_data_bytes(int num_payload_data_bytes);
+
+/* call this first */
+
+void horus_l2_init(void);
+
+/* returns number of output bytes in output_tx_data */
+
+int horus_l2_encode_tx_packet(unsigned char *output_tx_data,
+                              unsigned char *input_payload_data,
+                              int            num_payload_data_bytes);
+
+void horus_l2_decode_rx_packet(unsigned char *output_payload_data,
+                               unsigned char *input_rx_data,
+                               int            num_payload_data_bytes);
+
+unsigned short horus_l2_gen_crc16(unsigned char* data_p, unsigned char length);
+
+#endif

src/kiss_fft.c

@@ -0,0 +1,408 @@
+/*
+Copyright (c) 2003-2010, Mark Borgerding
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+    * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#include "_kiss_fft_guts.h"
+/* The guts header contains all the multiplication and addition macros that are defined for
+ fixed or floating point complex numbers.  It also delares the kf_ internal functions.
+ */
+
+static void kf_bfly2(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        int m
+        )
+{
+    kiss_fft_cpx * Fout2;
+    kiss_fft_cpx * tw1 = st->twiddles;
+    kiss_fft_cpx t;
+    Fout2 = Fout + m;
+    do{
+        C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2);
+
+        C_MUL (t,  *Fout2 , *tw1);
+        tw1 += fstride;
+        C_SUB( *Fout2 ,  *Fout , t );
+        C_ADDTO( *Fout ,  t );
+        ++Fout2;
+        ++Fout;
+    }while (--m);
+}
+
+static void kf_bfly4(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        const size_t m
+        )
+{
+    kiss_fft_cpx *tw1,*tw2,*tw3;
+    kiss_fft_cpx scratch[6];
+    size_t k=m;
+    const size_t m2=2*m;
+    const size_t m3=3*m;
+
+
+    tw3 = tw2 = tw1 = st->twiddles;
+
+    do {
+        C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4);
+
+        C_MUL(scratch[0],Fout[m] , *tw1 );
+        C_MUL(scratch[1],Fout[m2] , *tw2 );
+        C_MUL(scratch[2],Fout[m3] , *tw3 );
+
+        C_SUB( scratch[5] , *Fout, scratch[1] );
+        C_ADDTO(*Fout, scratch[1]);
+        C_ADD( scratch[3] , scratch[0] , scratch[2] );
+        C_SUB( scratch[4] , scratch[0] , scratch[2] );
+        C_SUB( Fout[m2], *Fout, scratch[3] );
+        tw1 += fstride;
+        tw2 += fstride*2;
+        tw3 += fstride*3;
+        C_ADDTO( *Fout , scratch[3] );
+
+        if(st->inverse) {
+            Fout[m].r = scratch[5].r - scratch[4].i;
+            Fout[m].i = scratch[5].i + scratch[4].r;
+            Fout[m3].r = scratch[5].r + scratch[4].i;
+            Fout[m3].i = scratch[5].i - scratch[4].r;
+        }else{
+            Fout[m].r = scratch[5].r + scratch[4].i;
+            Fout[m].i = scratch[5].i - scratch[4].r;
+            Fout[m3].r = scratch[5].r - scratch[4].i;
+            Fout[m3].i = scratch[5].i + scratch[4].r;
+        }
+        ++Fout;
+    }while(--k);
+}
+
+static void kf_bfly3(
+         kiss_fft_cpx * Fout,
+         const size_t fstride,
+         const kiss_fft_cfg st,
+         size_t m
+         )
+{
+     size_t k=m;
+     const size_t m2 = 2*m;
+     kiss_fft_cpx *tw1,*tw2;
+     kiss_fft_cpx scratch[5];
+     kiss_fft_cpx epi3;
+     epi3 = st->twiddles[fstride*m];
+
+     tw1=tw2=st->twiddles;
+
+     do{
+         C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
+
+         C_MUL(scratch[1],Fout[m] , *tw1);
+         C_MUL(scratch[2],Fout[m2] , *tw2);
+
+         C_ADD(scratch[3],scratch[1],scratch[2]);
+         C_SUB(scratch[0],scratch[1],scratch[2]);
+         tw1 += fstride;
+         tw2 += fstride*2;
+
+         Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
+         Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
+
+         C_MULBYSCALAR( scratch[0] , epi3.i );
+
+         C_ADDTO(*Fout,scratch[3]);
+
+         Fout[m2].r = Fout[m].r + scratch[0].i;
+         Fout[m2].i = Fout[m].i - scratch[0].r;
+
+         Fout[m].r -= scratch[0].i;
+         Fout[m].i += scratch[0].r;
+
+         ++Fout;
+     }while(--k);
+}
+
+static void kf_bfly5(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        int m
+        )
+{
+    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
+    int u;
+    kiss_fft_cpx scratch[13];
+    kiss_fft_cpx * twiddles = st->twiddles;
+    kiss_fft_cpx *tw;
+    kiss_fft_cpx ya,yb;
+    ya = twiddles[fstride*m];
+    yb = twiddles[fstride*2*m];
+
+    Fout0=Fout;
+    Fout1=Fout0+m;
+    Fout2=Fout0+2*m;
+    Fout3=Fout0+3*m;
+    Fout4=Fout0+4*m;
+
+    tw=st->twiddles;
+    for ( u=0; u<m; ++u ) {
+        C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
+        scratch[0] = *Fout0;
+
+        C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
+        C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
+        C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
+        C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
+
+        C_ADD( scratch[7],scratch[1],scratch[4]);
+        C_SUB( scratch[10],scratch[1],scratch[4]);
+        C_ADD( scratch[8],scratch[2],scratch[3]);
+        C_SUB( scratch[9],scratch[2],scratch[3]);
+
+        Fout0->r += scratch[7].r + scratch[8].r;
+        Fout0->i += scratch[7].i + scratch[8].i;
+
+        scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
+        scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
+
+        scratch[6].r =  S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
+        scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
+
+        C_SUB(*Fout1,scratch[5],scratch[6]);
+        C_ADD(*Fout4,scratch[5],scratch[6]);
+
+        scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
+        scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
+        scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
+        scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
+
+        C_ADD(*Fout2,scratch[11],scratch[12]);
+        C_SUB(*Fout3,scratch[11],scratch[12]);
+
+        ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
+    }
+}
+
+/* perform the butterfly for one stage of a mixed radix FFT */
+static void kf_bfly_generic(
+        kiss_fft_cpx * Fout,
+        const size_t fstride,
+        const kiss_fft_cfg st,
+        int m,
+        int p
+        )
+{
+    int u,k,q1,q;
+    kiss_fft_cpx * twiddles = st->twiddles;
+    kiss_fft_cpx t;
+    int Norig = st->nfft;
+
+    kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p);
+
+    for ( u=0; u<m; ++u ) {
+        k=u;
+        for ( q1=0 ; q1<p ; ++q1 ) {
+            scratch[q1] = Fout[ k  ];
+            C_FIXDIV(scratch[q1],p);
+            k += m;
+        }
+
+        k=u;
+        for ( q1=0 ; q1<p ; ++q1 ) {
+            int twidx=0;
+            Fout[ k ] = scratch[0];
+            for (q=1;q<p;++q ) {
+                twidx += fstride * k;
+                if (twidx>=Norig) twidx-=Norig;
+                C_MUL(t,scratch[q] , twiddles[twidx] );
+                C_ADDTO( Fout[ k ] ,t);
+            }
+            k += m;
+        }
+    }
+    KISS_FFT_TMP_FREE(scratch);
+}
+
+static
+void kf_work(
+        kiss_fft_cpx * Fout,
+        const kiss_fft_cpx * f,
+        const size_t fstride,
+        int in_stride,
+        int * factors,
+        const kiss_fft_cfg st
+        )
+{
+    kiss_fft_cpx * Fout_beg=Fout;
+    const int p=*factors++; /* the radix  */
+    const int m=*factors++; /* stage's fft length/p */
+    const kiss_fft_cpx * Fout_end = Fout + p*m;
+
+#ifdef _OPENMP
+    // use openmp extensions at the
+    // top-level (not recursive)
+    if (fstride==1 && p<=5)
+    {
+        int k;
+
+        // execute the p different work units in different threads
+#       pragma omp parallel for
+        for (k=0;k<p;++k)
+            kf_work( Fout +k*m, f+ fstride*in_stride*k,fstride*p,in_stride,factors,st);
+        // all threads have joined by this point
+
+        switch (p) {
+            case 2: kf_bfly2(Fout,fstride,st,m); break;
+            case 3: kf_bfly3(Fout,fstride,st,m); break;
+            case 4: kf_bfly4(Fout,fstride,st,m); break;
+            case 5: kf_bfly5(Fout,fstride,st,m); break;
+            default: kf_bfly_generic(Fout,fstride,st,m,p); break;
+        }
+        return;
+    }
+#endif
+
+    if (m==1) {
+        do{
+            *Fout = *f;
+            f += fstride*in_stride;
+        }while(++Fout != Fout_end );
+    }else{
+        do{
+            // recursive call:
+            // DFT of size m*p performed by doing
+            // p instances of smaller DFTs of size m,
+            // each one takes a decimated version of the input
+            kf_work( Fout , f, fstride*p, in_stride, factors,st);
+            f += fstride*in_stride;
+        }while( (Fout += m) != Fout_end );
+    }
+
+    Fout=Fout_beg;
+
+    // recombine the p smaller DFTs
+    switch (p) {
+        case 2: kf_bfly2(Fout,fstride,st,m); break;
+        case 3: kf_bfly3(Fout,fstride,st,m); break;
+        case 4: kf_bfly4(Fout,fstride,st,m); break;
+        case 5: kf_bfly5(Fout,fstride,st,m); break;
+        default: kf_bfly_generic(Fout,fstride,st,m,p); break;
+    }
+}
+
+/*  facbuf is populated by p1,m1,p2,m2, ...
+    where
+    p[i] * m[i] = m[i-1]
+    m0 = n                  */
+static
+void kf_factor(int n,int * facbuf)
+{
+    int p=4;
+    double floor_sqrt;
+    floor_sqrt = floorf( sqrtf((double)n) );
+
+    /*factor out powers of 4, powers of 2, then any remaining primes */
+    do {
+        while (n % p) {
+            switch (p) {
+                case 4: p = 2; break;
+                case 2: p = 3; break;
+                default: p += 2; break;
+            }
+            if (p > floor_sqrt)
+                p = n;          /* no more factors, skip to end */
+        }
+        n /= p;
+        *facbuf++ = p;
+        *facbuf++ = n;
+    } while (n > 1);
+}
+
+/*
+ *
+ * User-callable function to allocate all necessary storage space for the fft.
+ *
+ * The return value is a contiguous block of memory, allocated with malloc.  As such,
+ * It can be freed with free(), rather than a kiss_fft-specific function.
+ * */
+kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
+{
+    kiss_fft_cfg st=NULL;
+    size_t memneeded = sizeof(struct kiss_fft_state)
+        + sizeof(kiss_fft_cpx)*(nfft-1); /* twiddle factors*/
+
+    if ( lenmem==NULL ) {
+        st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
+    }else{
+        if (mem != NULL && *lenmem >= memneeded)
+            st = (kiss_fft_cfg)mem;
+        *lenmem = memneeded;
+    }
+    if (st) {
+        int i;
+        st->nfft=nfft;
+        st->inverse = inverse_fft;
+
+        for (i=0;i<nfft;++i) {
+            const double pi=3.141592653589793238462643383279502884197169399375105820974944;
+            double phase = -2*pi*i / nfft;
+            if (st->inverse)
+                phase *= -1;
+            kf_cexp(st->twiddles+i, phase );
+        }
+
+        kf_factor(nfft,st->factors);
+    }
+    return st;
+}
+
+
+void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
+{
+    if (fin == fout) {
+        //NOTE: this is not really an in-place FFT algorithm.
+        //It just performs an out-of-place FFT into a temp buffer
+        kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft);
+        kf_work(tmpbuf,fin,1,in_stride, st->factors,st);
+        memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft);
+        KISS_FFT_TMP_FREE(tmpbuf);
+    }else{
+        kf_work( fout, fin, 1,in_stride, st->factors,st );
+    }
+}
+
+void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
+{
+    kiss_fft_stride(cfg,fin,fout,1);
+}
+
+
+void kiss_fft_cleanup(void)
+{
+    // nothing needed any more
+}
+
+int kiss_fft_next_fast_size(int n)
+{
+    while(1) {
+        int m=n;
+        while ( (m%2) == 0 ) m/=2;
+        while ( (m%3) == 0 ) m/=3;
+        while ( (m%5) == 0 ) m/=5;
+        if (m<=1)
+            break; /* n is completely factorable by twos, threes, and fives */
+        n++;
+    }
+    return n;
+}

src/kiss_fft.h

@@ -0,0 +1,124 @@
+#ifndef KISS_FFT_H
+#define KISS_FFT_H
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ ATTENTION!
+ If you would like a :
+ -- a utility that will handle the caching of fft objects
+ -- real-only (no imaginary time component ) FFT
+ -- a multi-dimensional FFT
+ -- a command-line utility to perform ffts
+ -- a command-line utility to perform fast-convolution filtering
+
+ Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c
+  in the tools/ directory.
+*/
+
+#ifdef USE_SIMD
+# include <xmmintrin.h>
+# define kiss_fft_scalar __m128
+#define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16)
+#define KISS_FFT_FREE _mm_free
+#else
+#define KISS_FFT_MALLOC malloc
+#define KISS_FFT_FREE free
+#endif
+
+
+#ifdef FIXED_POINT
+#include <sys/types.h>
+# if (FIXED_POINT == 32)
+#  define kiss_fft_scalar int32_t
+# else
+#  define kiss_fft_scalar int16_t
+# endif
+#else
+# ifndef kiss_fft_scalar
+/*  default is float */
+#   define kiss_fft_scalar float
+# endif
+#endif
+
+typedef struct {
+    kiss_fft_scalar r;
+    kiss_fft_scalar i;
+}kiss_fft_cpx;
+
+typedef struct kiss_fft_state* kiss_fft_cfg;
+
+/*
+ *  kiss_fft_alloc
+ *
+ *  Initialize a FFT (or IFFT) algorithm's cfg/state buffer.
+ *
+ *  typical usage:      kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL);
+ *
+ *  The return value from fft_alloc is a cfg buffer used internally
+ *  by the fft routine or NULL.
+ *
+ *  If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc.
+ *  The returned value should be free()d when done to avoid memory leaks.
+ *
+ *  The state can be placed in a user supplied buffer 'mem':
+ *  If lenmem is not NULL and mem is not NULL and *lenmem is large enough,
+ *      then the function places the cfg in mem and the size used in *lenmem
+ *      and returns mem.
+ *
+ *  If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough),
+ *      then the function returns NULL and places the minimum cfg
+ *      buffer size in *lenmem.
+ * */
+
+kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem);
+
+/*
+ * kiss_fft(cfg,in_out_buf)
+ *
+ * Perform an FFT on a complex input buffer.
+ * for a forward FFT,
+ * fin should be  f[0] , f[1] , ... ,f[nfft-1]
+ * fout will be   F[0] , F[1] , ... ,F[nfft-1]
+ * Note that each element is complex and can be accessed like
+    f[k].r and f[k].i
+ * */
+void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);
+
+/*
+ A more generic version of the above function. It reads its input from every Nth sample.
+ * */
+void kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride);
+
+/* If kiss_fft_alloc allocated a buffer, it is one contiguous
+   buffer and can be simply free()d when no longer needed*/
+#define kiss_fft_free free
+
+/*
+ Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up
+ your compiler output to call this before you exit.
+*/
+void kiss_fft_cleanup(void);
+
+
+/*
+ * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5)
+ */
+int kiss_fft_next_fast_size(int n);
+
+/* for real ffts, we need an even size */
+#define kiss_fftr_next_fast_size_real(n) \
+        (kiss_fft_next_fast_size( ((n)+1)>>1)<<1)
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

src/kiss_fftr.c

@@ -0,0 +1,154 @@
+/*
+Copyright (c) 2003-2004, Mark Borgerding
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+    * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include "kiss_fftr.h"
+#include "_kiss_fft_guts.h"
+#include "assert.h"
+
+struct kiss_fftr_state{
+    kiss_fft_cfg substate;
+    kiss_fft_cpx * tmpbuf;
+    kiss_fft_cpx * super_twiddles;
+#ifdef USE_SIMD
+    void * pad;
+#endif
+};
+
+kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)
+{
+    int i;
+    kiss_fftr_cfg st = NULL;
+    size_t subsize, memneeded;
+
+    if (nfft & 1) {
+        fprintf(stderr,"Real FFT optimization must be even.\n");
+        return NULL;
+    }
+    nfft >>= 1;
+
+    kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize);
+    memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2);
+
+    if (lenmem == NULL) {
+        st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded);
+    } else {
+        if (*lenmem >= memneeded)
+            st = (kiss_fftr_cfg) mem;
+        *lenmem = memneeded;
+    }
+    if (!st)
+        return NULL;
+
+    st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */
+    st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize);
+    st->super_twiddles = st->tmpbuf + nfft;
+    kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize);
+
+    for (i = 0; i < nfft/2; ++i) {
+        float phase =
+            -3.14159265358979323846264338327 * ((float) (i+1) / nfft + .5);
+        if (inverse_fft)
+            phase *= -1;
+        kf_cexp (st->super_twiddles+i,phase);
+    }
+    return st;
+}
+
+void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
+{
+    /* input buffer timedata is stored row-wise */
+    int k,ncfft;
+    kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc;
+
+    assert(st->substate->inverse==0);
+
+    ncfft = st->substate->nfft;
+
+    /*perform the parallel fft of two real signals packed in real,imag*/
+    kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf );
+    /* The real part of the DC element of the frequency spectrum in st->tmpbuf
+     * contains the sum of the even-numbered elements of the input time sequence
+     * The imag part is the sum of the odd-numbered elements
+     *
+     * The sum of tdc.r and tdc.i is the sum of the input time sequence. 
+     *      yielding DC of input time sequence
+     * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... 
+     *      yielding Nyquist bin of input time sequence
+     */
+ 
+    tdc.r = st->tmpbuf[0].r;
+    tdc.i = st->tmpbuf[0].i;
+    C_FIXDIV(tdc,2);
+    CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i);
+    CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i);
+    freqdata[0].r = tdc.r + tdc.i;
+    freqdata[ncfft].r = tdc.r - tdc.i;
+#ifdef USE_SIMD    
+    freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0);
+#else
+    freqdata[ncfft].i = freqdata[0].i = 0;
+#endif
+
+    for ( k=1;k <= ncfft/2 ; ++k ) {
+        fpk    = st->tmpbuf[k]; 
+        fpnk.r =   st->tmpbuf[ncfft-k].r;
+        fpnk.i = - st->tmpbuf[ncfft-k].i;
+        C_FIXDIV(fpk,2);
+        C_FIXDIV(fpnk,2);
+
+        C_ADD( f1k, fpk , fpnk );
+        C_SUB( f2k, fpk , fpnk );
+        C_MUL( tw , f2k , st->super_twiddles[k-1]);
+
+        freqdata[k].r = HALF_OF(f1k.r + tw.r);
+        freqdata[k].i = HALF_OF(f1k.i + tw.i);
+        freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r);
+        freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i);
+    }
+}
+
+void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
+{
+    /* input buffer timedata is stored row-wise */
+    int k, ncfft;
+
+    assert(st->substate->inverse == 1);
+
+    ncfft = st->substate->nfft;
+
+    st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r;
+    st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r;
+    C_FIXDIV(st->tmpbuf[0],2);
+
+    for (k = 1; k <= ncfft / 2; ++k) {
+        kiss_fft_cpx fk, fnkc, fek, fok, tmp;
+        fk = freqdata[k];
+        fnkc.r = freqdata[ncfft - k].r;
+        fnkc.i = -freqdata[ncfft - k].i;
+        C_FIXDIV( fk , 2 );
+        C_FIXDIV( fnkc , 2 );
+
+        C_ADD (fek, fk, fnkc);
+        C_SUB (tmp, fk, fnkc);
+        C_MUL (fok, tmp, st->super_twiddles[k-1]);
+        C_ADD (st->tmpbuf[k],     fek, fok);
+        C_SUB (st->tmpbuf[ncfft - k], fek, fok);
+#ifdef USE_SIMD        
+        st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0);
+#else
+        st->tmpbuf[ncfft - k].i *= -1;
+#endif
+    }
+    kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
+}

src/kiss_fftr.h

@@ -0,0 +1,46 @@
+#ifndef KISS_FTR_H
+#define KISS_FTR_H
+
+#include "kiss_fft.h"
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+    
+/* 
+ 
+ Real optimized version can save about 45% cpu time vs. complex fft of a real seq.
+
+ 
+ 
+ */
+
+typedef struct kiss_fftr_state *kiss_fftr_cfg;
+
+
+kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem);
+/*
+ nfft must be even
+
+ If you don't care to allocate space, use mem = lenmem = NULL 
+*/
+
+
+void kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata);
+/*
+ input timedata has nfft scalar points
+ output freqdata has nfft/2+1 complex points
+*/
+
+void kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata);
+/*
+ input freqdata has  nfft/2+1 complex points
+ output timedata has nfft scalar points
+*/
+
+#define kiss_fftr_free free
+
+#ifdef __cplusplus
+}
+#endif
+#endif

src/modem_probe.c

@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: modem_probe.c
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 9 January 2016
+
+  Library to easily extract debug traces from modems during development and
+  verification
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include "comp.h"
+#include "octave.h"
+
+#define TRACE_I 1
+#define TRACE_F 2
+#define TRACE_C 3
+
+
+typedef struct probe_trace_info_s probe_trace_info;
+typedef struct datlink_s datlink;
+
+struct datlink_s{
+	void * data;
+	size_t len;
+	datlink * next;
+};
+
+struct probe_trace_info_s{
+	int type;
+	char name[255];
+	datlink * data;
+	datlink * last;
+	probe_trace_info *next;
+};
+
+static char *run = NULL;
+static char *mod = NULL;
+static probe_trace_info *first_trace = NULL;
+
+/* Init the probing library */
+void modem_probe_init_int(char *modname, char *runname){
+	mod = malloc((strlen(modname)+1)*sizeof(char));
+	run = malloc((strlen(runname)+1)*sizeof(char));
+	strcpy(run,runname);
+	strcpy(mod,modname);
+}
+
+/* 
+ * Gather the data stored in the linked list into a single blob,
+ * freeing links and buffers as it goes
+ */
+void * gather_data(datlink * d,size_t * len){
+	size_t size = 0;
+	datlink * cur = d;
+	datlink * next;
+	while(cur!=NULL){
+		size += d->len;
+		cur = cur->next;
+	}
+	cur = d;
+	size_t i = 0;
+	void * newbuf = malloc(size);
+	
+	while(cur!=NULL){
+		memcpy(newbuf+i,cur->data,cur->len);
+		i += cur->len;
+		free(cur->data);
+		next = cur->next;
+		free(cur);
+		cur = next;
+	}
+	*len = size;
+	return newbuf;
+}
+
+/* Dump all of the traces into a nice octave-able dump file */
+void modem_probe_close_int(){
+	if(run==NULL)
+		return;
+	
+	probe_trace_info *cur,*next;
+	cur = first_trace;
+	FILE * dumpfile = fopen(run,"w");
+	void * dbuf;
+	size_t len;
+	
+	while(cur != NULL){
+		dbuf = gather_data(cur->data,&len);
+		switch(cur->type){
+			case TRACE_I:
+				octave_save_int(dumpfile,cur->name,(int32_t*)dbuf,1,len/sizeof(int32_t));
+	        	break;
+			case TRACE_F:
+				octave_save_float(dumpfile,cur->name,(float*)dbuf,1,len/sizeof(float),10);
+				break;
+			case TRACE_C:
+				octave_save_complex(dumpfile,cur->name,(COMP*)dbuf,1,len/sizeof(COMP),10);
+				break;
+		}
+		next = cur->next;
+		free(cur);
+		free(dbuf);
+		cur = next;
+	}
+	
+	fclose(dumpfile);
+	free(run);
+	free(mod);
+}
+
+/* Look up or create a trace by name */
+probe_trace_info * modem_probe_get_trace(char * tracename){
+	probe_trace_info *cur,*npti;
+	
+	/* Make sure probe session is open */
+	if(run==NULL)
+		return NULL;
+	
+	cur = first_trace;
+	/* Walk through list, find trace with matching name */
+	while(cur != NULL){
+		/* We got one! */
+		if(strcmp( cur->name, tracename) == 0){
+			return cur;
+		}
+		cur = cur->next;
+	}
+	/* None found, open a new trace */
+	
+	npti = (probe_trace_info *) malloc(sizeof(probe_trace_info));
+	npti->next = first_trace;
+	npti->data = NULL;
+	npti->last = NULL;
+	strcpy(npti->name,tracename);
+	first_trace = npti;
+	
+	return npti;
+	
+}
+
+
+void modem_probe_samp_i_int(char * tracename,int32_t samp[],size_t cnt){
+	probe_trace_info *pti;
+	datlink *ndat;
+	
+	pti = modem_probe_get_trace(tracename);
+	if(pti == NULL)
+		return;
+	
+	pti->type = TRACE_I;
+	
+	ndat = (datlink*) malloc(sizeof(datlink));
+	ndat->data = malloc(sizeof(int32_t)*cnt);
+	
+	ndat->len = cnt*sizeof(int32_t);
+	ndat->next = NULL;
+	memcpy(ndat->data,(void*)&(samp[0]),sizeof(int32_t)*cnt);
+	
+	if(pti->last!=NULL){
+		pti->last->next = ndat;
+		pti->last = ndat;
+	} else {
+		pti->data = ndat;
+		pti->last = ndat;
+	}
+	
+}
+
+void modem_probe_samp_f_int(char * tracename,float samp[],size_t cnt){
+	probe_trace_info *pti;
+	datlink *ndat;
+	
+	pti = modem_probe_get_trace(tracename);
+	if(pti == NULL)
+		return;
+	
+	pti->type = TRACE_F;
+	
+	ndat = (datlink*) malloc(sizeof(datlink));
+	ndat->data = malloc(sizeof(float)*cnt);
+	
+	ndat->len = cnt*sizeof(float);
+	ndat->next = NULL;
+	memcpy(ndat->data,(void*)&(samp[0]),sizeof(float)*cnt);
+	
+	if(pti->last!=NULL){
+		pti->last->next = ndat;
+		pti->last = ndat;
+	} else {
+		pti->data = ndat;
+		pti->last = ndat;
+	}
+}
+	
+void modem_probe_samp_c_int(char * tracename,COMP samp[],size_t cnt){
+	probe_trace_info *pti;
+	datlink *ndat;
+	
+	pti = modem_probe_get_trace(tracename);
+	if(pti == NULL)
+		return;
+	
+	pti->type = TRACE_C;
+	
+	ndat = (datlink*) malloc(sizeof(datlink));
+	ndat->data = malloc(sizeof(COMP)*cnt);
+	
+	ndat->len = cnt*sizeof(COMP);
+	ndat->next = NULL;
+	memcpy(ndat->data,(void*)&(samp[0]),sizeof(COMP)*cnt);
+	
+	if(pti->last!=NULL){
+		pti->last->next = ndat;
+		pti->last = ndat;
+	} else {
+		pti->data = ndat;
+		pti->last = ndat;
+	}
+}

src/modem_probe.h

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: modem_probe.h
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 9 January 2016
+
+  Library to easily extract debug traces from modems during development
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __MODEMPROBE_H
+#define __MODEMPROBE_H
+
+#include <stdint.h>
+#include <stdlib.h>
+#include <complex.h>
+#include "comp.h"
+
+#ifdef MODEMPROBE_ENABLE
+
+/* Internal functions */
+void modem_probe_init_int(char *modname, char *runname);
+void modem_probe_close_int();
+
+void modem_probe_samp_i_int(char * tracename,int samp[],size_t cnt);
+void modem_probe_samp_f_int(char * tracename,float samp[],size_t cnt);
+void modem_probe_samp_c_int(char * tracename,COMP samp[],size_t cnt);
+
+/* 
+ * Init the probe library.
+ * char *modname - Name of the modem under test
+ * char *runname - Name/path of the file data is dumped to
+ */
+static inline void modem_probe_init(char *modname,char *runname){
+        modem_probe_init_int(modname,runname);
+}
+
+/*
+ * Dump traces to a file and clean up
+ */
+static inline void modem_probe_close(){
+        modem_probe_close_int();
+}
+
+/*
+ * Save some number of int samples to a named trace
+ * char *tracename - name of trace being saved to
+ * int samp[] - int samples
+ * size_t cnt - how many samples to save
+ */
+static inline void modem_probe_samp_i(char *tracename,int samp[],size_t cnt){
+        modem_probe_samp_i_int(tracename,samp,cnt);
+}
+
+/*
+ * Save some number of float samples to a named trace
+ * char *tracename - name of trace being saved to
+ * float samp[] - int samples
+ * size_t cnt - how many samples to save
+ */
+static inline void modem_probe_samp_f(char *tracename,float samp[],size_t cnt){
+        modem_probe_samp_f_int(tracename,samp,cnt);
+}       
+
+/*
+ * Save some number of complex samples to a named trace
+ * char *tracename - name of trace being saved to
+ * COMP samp[] - int samples
+ * size_t cnt - how many samples to save
+ */
+static inline void modem_probe_samp_c(char *tracename,COMP samp[],size_t cnt){
+        modem_probe_samp_c_int(tracename,samp,cnt);
+}
+
+/*
+ * Save some number of complex samples to a named trace
+ * char *tracename - name of trace being saved to
+ * float complex samp[] - int samples
+ * size_t cnt - how many samples to save
+ */
+static inline void modem_probe_samp_cft(char *tracename,complex float samp[],size_t cnt){
+        modem_probe_samp_c_int(tracename,(COMP*)samp,cnt);
+}
+
+#else
+
+static inline void modem_probe_init(char *modname,char *runname){
+        return;
+}
+
+static inline void modem_probe_close(){
+        return;
+}
+
+static inline void modem_probe_samp_i(char *name,int samp[],size_t sampcnt){
+        return;
+}
+
+static inline void modem_probe_samp_f(char *name,float samp[],size_t cnt){
+        return;
+}
+
+static inline void modem_probe_samp_c(char *name,COMP samp[],size_t cnt){
+        return;
+}
+
+static inline void modem_probe_samp_cft(char *name,complex float samp[],size_t cnt){
+        return;
+}
+
+#endif
+
+#endif

src/modem_stats.c

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: modem_stats.c
+  AUTHOR......: David Rowe
+  DATE CREATED: June 2015
+
+  Common functions for returning demod stats from fdmdv and cohpsk modems.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2015 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <math.h>
+#include "modem_stats.h"
+#include "codec2_fdmdv.h"
+
+void modem_stats_open(struct MODEM_STATS *f)
+{
+    int i;
+
+    /* zero out all the stats */
+
+    memset(f, 0, sizeof(struct MODEM_STATS));
+
+    /* init the FFT */
+    
+#ifndef __EMBEDDED__
+    for(i=0; i<2*MODEM_STATS_NSPEC; i++)
+	f->fft_buf[i] = 0.0;
+    f->fft_cfg = kiss_fft_alloc (2*MODEM_STATS_NSPEC, 0, NULL, NULL);
+    assert(f->fft_cfg != NULL);
+#endif
+}
+
+void modem_stats_close(struct MODEM_STATS *f)
+{
+#ifndef __EMBEDDED__
+    KISS_FFT_FREE(f->fft_cfg);
+#endif
+}
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: modem_stats_get_rx_spectrum()
+  AUTHOR......: David Rowe
+  DATE CREATED: 9 June 2012
+
+  Returns the MODEM_STATS_NSPEC point magnitude spectrum of the rx signal in
+  dB. The spectral samples are scaled so that 0dB is the peak, a good
+  range for plotting is 0 to -40dB.
+
+  Note only the real part of the complex input signal is used at
+  present.  A complex variable is used for input for compatability
+  with the other rx signal procesing.
+
+  Successive calls can be used to build up a waterfall or spectrogram
+  plot, by mapping the received levels to colours.
+
+  The time-frequency resolution of the spectrum can be adjusted by varying
+  MODEM_STATS_NSPEC.  Note that a 2* MODEM_STATS_NSPEC size FFT is reqd to get
+  MODEM_STATS_NSPEC output points.  MODEM_STATS_NSPEC must be a power of 2.
+
+  See octave/tget_spec.m for a demo real time spectral display using
+  Octave. This demo averages the output over time to get a smoother
+  display:
+
+     av = 0.9*av + 0.1*mag_dB
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef __EMBEDDED__
+void modem_stats_get_rx_spectrum(struct MODEM_STATS *f, float mag_spec_dB[], COMP rx_fdm[], int nin)
+{
+    int   i,j;
+    COMP  fft_in[2*MODEM_STATS_NSPEC];
+    COMP  fft_out[2*MODEM_STATS_NSPEC];
+    float full_scale_dB;
+
+    /* update buffer of input samples */
+
+    for(i=0; i<2*MODEM_STATS_NSPEC-nin; i++)
+	f->fft_buf[i] = f->fft_buf[i+nin];
+    for(j=0; j<nin; j++,i++)
+	f->fft_buf[i] = rx_fdm[j].real;
+    assert(i == 2*MODEM_STATS_NSPEC);
+
+    /* window and FFT */
+
+    for(i=0; i<2*MODEM_STATS_NSPEC; i++) {
+	fft_in[i].real = f->fft_buf[i] * (0.5 - 0.5*cosf((float)i*2.0*M_PI/(2*MODEM_STATS_NSPEC)));
+	fft_in[i].imag = 0.0;
+    }
+
+    kiss_fft(f->fft_cfg, (kiss_fft_cpx *)fft_in, (kiss_fft_cpx *)fft_out);
+
+    /* FFT scales up a signal of level 1 FDMDV_NSPEC */
+
+    full_scale_dB = 20*log10(MODEM_STATS_NSPEC*FDMDV_SCALE);
+
+    /* scale and convert to dB */
+
+    for(i=0; i<MODEM_STATS_NSPEC; i++) {
+	mag_spec_dB[i]  = 10.0*log10f(fft_out[i].real*fft_out[i].real + fft_out[i].imag*fft_out[i].imag + 1E-12);
+	mag_spec_dB[i] -= full_scale_dB;
+    }
+}
+#endif

src/modem_stats.h

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: modem_stats.h
+  AUTHOR......: David Rowe
+  DATE CREATED: June 2015
+
+  Common structure for returning demod stats from fdmdv and cohpsk modems.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2015 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __MODEM_STATS__
+#define __MODEM_STATS__
+
+#include "comp.h"
+#include "kiss_fft.h"
+
+#ifdef __cplusplus
+  extern "C" {
+#endif
+
+#define MODEM_STATS_NC_MAX      50
+#define MODEM_STATS_NR_MAX      8
+#define MODEM_STATS_ET_MAX      8
+#define MODEM_STATS_EYE_IND_MAX 160     
+#define MODEM_STATS_NSPEC       512
+#define MODEM_STATS_MAX_F_HZ    4000
+#define MODEM_STATS_MAX_F_EST   4
+      
+struct MODEM_STATS {
+    int    Nc;
+    float  snr_est;                          /* estimated SNR of rx signal in dB (3 kHz noise BW)  */
+#ifndef __EMBEDDED__
+    COMP   rx_symbols[MODEM_STATS_NR_MAX][MODEM_STATS_NC_MAX+1];
+                                             /* latest received symbols, for scatter plot          */
+#endif
+  int    nr;                               /* number of rows in rx_symbols                       */
+    int    sync;                             /* demod sync state                                   */
+    float  foff;                             /* estimated freq offset in Hz                        */
+    float  rx_timing;                        /* estimated optimum timing offset in samples         */
+    float  clock_offset;                     /* Estimated tx/rx sample clock offset in ppm         */
+    float  sync_metric;                      /* number between 0 and 1 indicating quality of sync  */
+    
+    /* eye diagram traces */
+    /* Eye diagram plot -- first dim is trace number, second is the trace idx */
+#ifndef __EMBEDDED__
+    float  rx_eye[MODEM_STATS_ET_MAX][MODEM_STATS_EYE_IND_MAX];
+    int    neyetr;                           /* How many eye traces are plotted */
+    int    neyesamp;                         /* How many samples in the eye diagram */
+
+    /* optional for FSK modems - est tone freqs */
+
+    float f_est[MODEM_STATS_MAX_F_EST];
+#endif
+    
+    /* Buf for FFT/waterfall */
+
+#ifndef __EMBEDDED__
+   float        fft_buf[2*MODEM_STATS_NSPEC];
+   kiss_fft_cfg fft_cfg;
+#endif
+};
+
+void modem_stats_open(struct MODEM_STATS *f);
+void modem_stats_close(struct MODEM_STATS *f);
+void modem_stats_get_rx_spectrum(struct MODEM_STATS *f, float mag_spec_dB[], COMP rx_fdm[], int nin);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
+

src/mpdecode_core.c

@@ -0,0 +1,709 @@
+/*
+  FILE...: mpdecode_core.c
+  AUTHOR.: Matthew C. Valenti, Rohit Iyer Seshadri, David Rowe
+  CREATED: Sep 2016
+
+  C-callable core functions moved from MpDecode.c, so they can be used for
+  Octave and C programs.
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <assert.h>
+#include "mpdecode_core.h"
+#ifndef USE_ORIGINAL_PHI0
+#include "phi0.h"
+#endif
+
+#include "debug_alloc.h"
+
+#ifdef __EMBEDDED__
+#include "machdep.h"
+#endif
+
+#define QPSK_CONSTELLATION_SIZE 4
+#define QPSK_BITS_PER_SYMBOL    2
+
+
+/* QPSK constellation for symbol likelihood calculations */
+
+static COMP S_matrix[] = {
+    { 1.0f,  0.0f},
+    { 0.0f,  1.0f},
+    { 0.0f, -1.0f},
+    {-1.0f,  0.0f}
+};
+
+// c_nodes will be an array of NumberParityBits of struct c_node
+// Each c_node contains an array of <degree> c_sub_node elements
+// This structure reduces the indexing caluclations in SumProduct()
+
+struct c_sub_node { // Order is important here to keep total size small.
+  uint16_t index;   // Values from H_rows (except last 2 entries)
+  uint16_t socket;  // The socket number at the v_node
+  float    message; // modified during operation!
+};
+
+struct c_node {
+  int degree;       // A count of elements in the following arrays
+  struct c_sub_node *subs;
+};
+
+// v_nodes will be an array of CodeLength of struct v_node
+
+struct v_sub_node {
+  uint16_t index;  //    the index of a c_node it is connected to
+                   //    Filled with values from H_cols (except last 2 entries)
+  uint16_t socket; //    socket number at the c_node
+  float message;   //    Loaded with input data
+                   //    modified during operation!
+  uint8_t sign;    //    1 if input is negative
+                   //    modified during operation!
+};
+
+struct v_node {
+  int degree;       // A count of ???
+  float initial_value;
+  struct v_sub_node *subs;
+};
+
+void encode(struct LDPC *ldpc, unsigned char ibits[], unsigned char pbits[]) {
+    unsigned int p, i, tmp, par, prev=0;
+    int          ind;
+    uint16_t     *H_rows = ldpc->H_rows;
+
+    for (p=0; p<ldpc->NumberParityBits; p++) {
+        par = 0;
+
+        for (i=0; i<ldpc->max_row_weight; i++) {
+            ind = H_rows[p + i*ldpc->NumberParityBits];
+            par = par + ibits[ind-1];
+        }
+
+        tmp = par + prev;
+
+        tmp &= 1;    // only retain the lsb
+        prev = tmp;
+        pbits[p] = tmp;
+    }
+}
+
+#ifdef USE_ORIGINAL_PHI0
+/* Phi function */
+static float phi0(
+                  float x )
+{
+  float z;
+
+  if (x>10)
+    return( 0 );
+  else if (x< 9.08e-5 )
+    return( 10 );
+  else if (x > 9)
+    return( 1.6881e-4 );
+  /* return( 1.4970e-004 ); */
+  else if (x > 8)
+    return( 4.5887e-4 );
+  /* return( 4.0694e-004 ); */
+  else if (x > 7)
+    return( 1.2473e-3 );
+  /* return( 1.1062e-003 ); */
+  else if (x > 6)
+    return( 3.3906e-3 );
+  /* return( 3.0069e-003 ); */
+  else if (x > 5)
+    return( 9.2168e-3 );
+  /* return( 8.1736e-003 ); */
+  else {
+    z = (float) exp(x);
+    return( (float) log( (z+1)/(z-1) ) );
+  }
+}
+#endif
+
+
+/* Values for linear approximation (DecoderType=5) */
+
+#define AJIAN -0.24904163195436
+#define TJIAN 2.50681740420944
+
+/* The linear-log-MAP algorithm */
+
+static float max_star0(
+                       float delta1,
+                       float delta2 )
+{
+    register float diff;
+
+    diff = delta2 - delta1;
+
+    if ( diff > TJIAN )
+        return( delta2 );
+    else if ( diff < -TJIAN )
+        return( delta1 );
+    else if ( diff > 0 )
+        return( delta2 + AJIAN*(diff-TJIAN) );
+    else
+        return( delta1 - AJIAN*(diff+TJIAN) );
+}
+
+void init_c_v_nodes(struct c_node *c_nodes,
+                    int     shift,
+                    int     NumberParityBits,
+                    int     max_row_weight,
+                    uint16_t *H_rows,
+                    int     H1,
+                    int     CodeLength,
+                    struct v_node *v_nodes,
+                    int     NumberRowsHcols,
+                    uint16_t *H_cols,
+                    int     max_col_weight,
+                    int     dec_type,
+                    float  *input)
+{
+    int i, j, k, count, cnt, c_index, v_index;
+
+    /* first determine the degree of each c-node */
+
+    if (shift ==0){
+        for (i=0;i<NumberParityBits;i++) {
+            count = 0;
+            for (j=0;j<max_row_weight;j++) {
+                if ( H_rows[i+j*NumberParityBits] > 0 ) {
+                    count++;
+                }
+            }
+            c_nodes[i].degree = count;
+            if (H1){
+                if (i==0){
+                    c_nodes[i].degree=count+1;
+                }
+                else{
+                    c_nodes[i].degree=count+2;
+                }
+            }
+        }
+    }
+    else{
+        cnt=0;
+        for (i=0;i<(NumberParityBits/shift);i++) {
+            for (k=0;k<shift;k++){
+                count = 0;
+                for (j=0;j<max_row_weight;j++) {
+                    if ( H_rows[cnt+j*NumberParityBits] > 0 ) {
+                        count++;
+                    }
+                }
+                c_nodes[cnt].degree = count;
+                if ((i==0)||(i==(NumberParityBits/shift)-1)){
+                    c_nodes[cnt].degree=count+1;
+                }
+                else{
+                    c_nodes[cnt].degree=count+2;
+                }
+                cnt++;
+            }
+        }
+    }
+
+    if (H1){
+
+        if (shift ==0){
+            for (i=0;i<NumberParityBits;i++) {
+
+                // Allocate sub nodes
+                c_nodes[i].subs = CALLOC(c_nodes[i].degree, sizeof(struct c_sub_node));
+                assert(c_nodes[i].subs);
+
+	        // Populate sub nodes
+                for (j=0;j<c_nodes[i].degree-2;j++) {
+                    c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
+                }
+                j=c_nodes[i].degree-2;
+
+                if (i==0){
+                    c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
+                }
+                else {
+                    c_nodes[i].subs[j].index = (CodeLength-NumberParityBits)+i-1;
+                }
+
+                j=c_nodes[i].degree-1;
+                c_nodes[i].subs[j].index = (CodeLength-NumberParityBits)+i;
+
+            }
+        }
+        if (shift >0){
+            cnt=0;
+            for (i=0;i<(NumberParityBits/shift);i++){
+
+                for (k =0;k<shift;k++){
+
+                    // Allocate sub nodes
+                    c_nodes[cnt].subs = CALLOC(c_nodes[cnt].degree, sizeof(struct c_sub_node));
+                    assert(c_nodes[cnt].subs);
+
+	            // Populate sub nodes
+                    for (j=0;j<c_nodes[cnt].degree-2;j++) {
+                        c_nodes[cnt].subs[j].index = (H_rows[cnt+j*NumberParityBits] - 1);
+                    }
+                    j=c_nodes[cnt].degree-2;
+                    if ((i ==0)||(i==(NumberParityBits/shift-1))){
+                        c_nodes[cnt].subs[j].index = (H_rows[cnt+j*NumberParityBits] - 1);
+                    }
+                    else{
+                        c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i);
+                    }
+                    j=c_nodes[cnt].degree-1;
+                    c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i+1);
+                    if (i== (NumberParityBits/shift-1))
+                        {
+                            c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i);
+                        }
+                    cnt++;
+                }
+            }
+        }
+
+    } else {
+        for (i=0;i<NumberParityBits;i++) {
+            // Allocate sub nodes
+            c_nodes[i].subs = CALLOC(c_nodes[i].degree, sizeof(struct c_sub_node));
+            assert(c_nodes[i].subs);
+
+	    // Populate sub nodes
+            for (j=0;j<c_nodes[i].degree;j++){
+                c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
+            }
+        }
+    }
+
+
+    /* determine degree of each v-node */
+
+    for(i=0;i<(CodeLength-NumberParityBits+shift);i++){
+        count=0;
+        for (j=0;j<max_col_weight;j++) {
+            if ( H_cols[i+j*NumberRowsHcols] > 0 ) {
+                count++;
+            }
+        }
+        v_nodes[i].degree = count;
+    }
+
+    for(i=CodeLength-NumberParityBits+shift;i<CodeLength;i++){
+        count=0;
+        if (H1){
+            if(i!=CodeLength-1){
+                v_nodes[i].degree=2;
+            }  else{
+                v_nodes[i].degree=1;
+            }
+
+        } else{
+            for (j=0;j<max_col_weight;j++) {
+                if ( H_cols[i+j*NumberRowsHcols] > 0 ) {
+                    count++;
+                }
+            }
+            v_nodes[i].degree = count;
+        }
+    }
+
+    if (shift>0){
+        v_nodes[CodeLength-1].degree =v_nodes[CodeLength-1].degree+1;
+    }
+
+
+    /* set up v_nodes */
+
+    for (i=0;i<CodeLength;i++) {
+        // Allocate sub nodes
+        v_nodes[i].subs = CALLOC(v_nodes[i].degree, sizeof(struct v_sub_node));
+        assert(v_nodes[i].subs);
+
+	// Populate sub nodes
+
+        /* index tells which c-nodes this v-node is connected to */
+        v_nodes[i].initial_value = input[i];
+        count=0;
+
+        for (j=0;j<v_nodes[i].degree;j++) {
+            if ((H1)&& (i>=CodeLength-NumberParityBits+shift)){
+                v_nodes[i].subs[j].index=i-(CodeLength-NumberParityBits+shift)+count;
+                if (shift ==0){
+                    count=count+1;
+                }
+                else{
+                    count=count+shift;
+                }
+            } else  {
+                v_nodes[i].subs[j].index = (H_cols[i+j*NumberRowsHcols] - 1);
+            }
+
+            /* search the connected c-node for the proper message value */
+            for (c_index=0;c_index<c_nodes[ v_nodes[i].subs[j].index ].degree;c_index++)
+                if ( c_nodes[ v_nodes[i].subs[j].index ].subs[c_index].index == i ) {
+                    v_nodes[i].subs[j].socket = c_index;
+                    break;
+                }
+            /* initialize v-node with received LLR */
+            if ( dec_type == 1)
+                v_nodes[i].subs[j].message = fabs(input[i]);
+            else
+                v_nodes[i].subs[j].message = phi0( fabs(input[i]) );
+
+            if (input[i] < 0)
+                v_nodes[i].subs[j].sign = 1;
+        }
+
+    }
+
+
+
+    /* now finish setting up the c_nodes */
+    for (i=0;i<NumberParityBits;i++) {
+        /* index tells which v-nodes this c-node is connected to */
+        for (j=0;j<c_nodes[i].degree;j++) {
+            /* search the connected v-node for the proper message value */
+            for (v_index=0;v_index<v_nodes[ c_nodes[i].subs[j].index ].degree;v_index++)
+                if (v_nodes[ c_nodes[i].subs[j].index ].subs[v_index].index == i ) {
+                    c_nodes[i].subs[j].socket = v_index;
+                    break;
+                }
+        }
+    }
+
+}
+
+
+///////////////////////////////////////
+/* function for doing the MP decoding */
+// Returns the iteration count
+int SumProduct( int       *parityCheckCount,
+                char     DecodedBits[],
+                struct c_node c_nodes[],
+                struct v_node v_nodes[],
+                int       CodeLength,
+                int       NumberParityBits,
+                int       max_iter,
+                float    r_scale_factor,
+                float    q_scale_factor,
+                int      data[] )
+{
+  int result;
+  int bitErrors;
+  int i,j, iter;
+  float phi_sum;
+  int sign;
+  float temp_sum;
+  float Qi;
+  int   ssum;
+
+
+  result = max_iter;
+  for (iter=0;iter<max_iter;iter++) {
+
+    for(i=0; i<CodeLength; i++) DecodedBits[i] = 0; // Clear each pass!
+    bitErrors = 0;
+
+    /* update r */
+    ssum = 0;
+    for (j=0;j<NumberParityBits;j++) {
+      sign = v_nodes[ c_nodes[j].subs[0].index ].subs[ c_nodes[j].subs[0].socket ].sign;
+      phi_sum = v_nodes[ c_nodes[j].subs[0].index ].subs[ c_nodes[j].subs[0].socket ].message;
+
+      for (i=1;i<c_nodes[j].degree;i++) {
+        // Compiler should optomize this but write the best we can to start from.
+        struct c_sub_node *cp = &c_nodes[j].subs[i];
+        struct v_sub_node *vp = &v_nodes[ cp->index ].subs[ cp->socket ];
+	    phi_sum += vp->message;
+	    sign ^= vp->sign;
+      }
+
+      if (sign==0) ssum++;
+
+      for (i=0;i<c_nodes[j].degree;i++) {
+        struct c_sub_node *cp = &c_nodes[j].subs[i];
+        struct v_sub_node *vp = &v_nodes[ cp->index ].subs[ cp->socket ];
+	    if ( sign ^ vp->sign ) {
+	      cp->message = -phi0( phi_sum - vp->message ); // *r_scale_factor;
+        } else
+	      cp->message =  phi0( phi_sum - vp->message ); // *r_scale_factor;
+      }
+    }
+
+    /* update q */
+    for (i=0;i<CodeLength;i++) {
+
+      /* first compute the LLR */
+      Qi = v_nodes[i].initial_value;
+      for (j=0;j<v_nodes[i].degree;j++) {
+        struct v_sub_node *vp = &v_nodes[i].subs[j];
+	    Qi += c_nodes[ vp->index ].subs[ vp->socket ].message;
+      }
+
+      /* make hard decision */
+      if (Qi < 0) {
+            DecodedBits[i] = 1;
+      }
+
+      /* now subtract to get the extrinsic information */
+      for (j=0;j<v_nodes[i].degree;j++) {
+        struct v_sub_node *vp = &v_nodes[i].subs[j];
+	    temp_sum = Qi - c_nodes[ vp->index ].subs[ vp->socket ].message;
+
+	    vp->message = phi0( fabs( temp_sum ) ); // *q_scale_factor;
+        if (temp_sum > 0)
+	      vp->sign = 0;
+        else
+	      vp->sign = 1;
+      }
+    }
+
+    /* count data bit errors, assuming that it is systematic */
+    for (i=0;i<CodeLength-NumberParityBits;i++)
+      if ( DecodedBits[i] != data[i] )
+            bitErrors++;
+
+
+    /* Halt if zero errors */
+    if (bitErrors == 0) {
+      result = iter + 1;
+      break;
+    }
+
+    // count the number of PC satisfied and exit if all OK
+    *parityCheckCount = ssum;
+    if (ssum==NumberParityBits)  {
+      result = iter + 1;
+      break;
+    }
+
+
+  }
+
+return(result);
+}
+
+
+/* Convenience function to call LDPC decoder from C programs */
+
+int run_ldpc_decoder(struct LDPC *ldpc, uint8_t out_char[], float input[], int *parityCheckCount) {
+    int         max_iter, dec_type;
+    float       q_scale_factor, r_scale_factor;
+    int         max_row_weight, max_col_weight;
+    int         CodeLength, NumberParityBits, NumberRowsHcols, shift, H1;
+    int         i;
+    struct c_node *c_nodes;
+    struct v_node *v_nodes;
+
+    /* default values */
+
+    max_iter  = ldpc->max_iter;
+    dec_type  = ldpc->dec_type;
+    q_scale_factor = ldpc->q_scale_factor;
+    r_scale_factor = ldpc->r_scale_factor;
+
+    CodeLength = ldpc->CodeLength;                    /* length of entire codeword */
+    NumberParityBits = ldpc->NumberParityBits;
+    NumberRowsHcols = ldpc->NumberRowsHcols;
+
+    char *DecodedBits = CALLOC( CodeLength, sizeof( char ) );
+    assert(DecodedBits);
+
+    /* derive some parameters */
+
+    shift = (NumberParityBits + NumberRowsHcols) - CodeLength;
+    if (NumberRowsHcols == CodeLength) {
+        H1=0;
+        shift=0;
+    } else {
+        H1=1;
+    }
+
+    max_row_weight = ldpc->max_row_weight;
+    max_col_weight = ldpc->max_col_weight;
+
+    /* initialize c-node and v-node structures */
+
+    c_nodes = CALLOC( NumberParityBits, sizeof( struct c_node ) );
+    assert(c_nodes);
+    v_nodes = CALLOC( CodeLength, sizeof( struct v_node));
+    assert(v_nodes);
+
+    init_c_v_nodes(c_nodes, shift, NumberParityBits, max_row_weight, ldpc->H_rows, H1, CodeLength,
+                   v_nodes, NumberRowsHcols, ldpc->H_cols, max_col_weight, dec_type, input);
+
+    int DataLength = CodeLength - NumberParityBits;
+    int *data_int = CALLOC( DataLength, sizeof(int) );
+
+    /* need to clear these on each call */
+
+    for(i=0; i<CodeLength; i++) DecodedBits[i] = 0;
+
+    /* Call function to do the actual decoding */
+    int iter = SumProduct( parityCheckCount, DecodedBits, c_nodes, v_nodes, 
+                           CodeLength, NumberParityBits, max_iter, 
+                           r_scale_factor, q_scale_factor, data_int );
+
+    for (i=0; i<CodeLength; i++) out_char[i] = DecodedBits[i];
+
+    /* Clean up memory */
+
+    FREE(DecodedBits);
+    FREE( data_int );
+
+    for (i=0;i<NumberParityBits;i++) FREE( c_nodes[i].subs );
+    FREE( c_nodes );
+
+    for (i=0;i<CodeLength;i++) FREE( v_nodes[i].subs);
+    FREE( v_nodes );
+
+    return iter;
+}
+
+
+void sd_to_llr(float llr[], double sd[], int n) {
+    double sum, mean, sign, sumsq, estvar, estEsN0, x;
+    int i;
+
+    /* convert SD samples to LLRs -------------------------------*/
+
+    sum = 0.0;
+    for(i=0; i<n; i++)
+        sum += fabs(sd[i]);
+    mean = sum/n;
+
+    /* find variance from +/-1 symbol position */
+
+    sum = sumsq = 0.0;
+    for(i=0; i<n; i++) {
+        sign = (sd[i] > 0.0L) - (sd[i] < 0.0L);
+        x = (sd[i]/mean - sign);
+        sum += x;
+        sumsq += x*x;
+    }
+    estvar = (n * sumsq - sum * sum) / (n * (n - 1));
+    //fprintf(stderr, "mean: %f var: %f\n", mean, estvar);
+
+    estEsN0 = 1.0/(2.0L * estvar + 1E-3);
+    for(i=0; i<n; i++)
+        llr[i] = 4.0L * estEsN0 * sd[i];
+}
+
+
+/*
+   Determine symbol likelihood from received QPSK symbols.
+
+   Notes:
+
+   1) We assume fading[] is real, it is also possible to compute
+      with complex fading, see CML library Demod2D.c source code.
+   2) Using floats instead of doubles, for stm32.
+      Testing shows good BERs with floats.
+*/
+
+void Demod2D(float   symbol_likelihood[],       /* output, M*number_symbols              */
+             COMP    r[],                       /* received QPSK symbols, number_symbols */
+             COMP    S_matrix[],                /* constellation of size M               */
+             float   EsNo,
+             float   fading[],                  /* real fading values, number_symbols    */
+             float   mean_amp,
+             int     number_symbols)
+{
+    int     M=QPSK_CONSTELLATION_SIZE;
+    int     i,j;
+    float  tempsr, tempsi, Er, Ei;
+
+    /* determine output */
+
+    for (i=0;i<number_symbols;i++) {                /* go through each received symbol */
+        for (j=0;j<M;j++) {                         /* each postulated symbol          */
+            tempsr = fading[i]*S_matrix[j].real/mean_amp;
+            tempsi = fading[i]*S_matrix[j].imag/mean_amp;
+            Er = r[i].real/mean_amp - tempsr;
+            Ei = r[i].imag/mean_amp - tempsi;
+            symbol_likelihood[i*M+j] = -EsNo*(Er*Er+Ei*Ei);
+            //printf("symbol_likelihood[%d][%d] = %f\n", i,j,symbol_likelihood[i*M+j]);
+        }
+        //exit(0);
+    }
+
+}
+
+
+void Somap(float  bit_likelihood[],      /* number_bits, bps*number_symbols */
+           float  symbol_likelihood[],   /* M*number_symbols                */
+           int     number_symbols)
+{
+    int    M=QPSK_CONSTELLATION_SIZE, bps = QPSK_BITS_PER_SYMBOL;
+    int    n,i,j,k,mask;
+    float num[bps], den[bps];
+    float metric;
+
+    for (n=0; n<number_symbols; n++) { /* loop over symbols */
+        for (k=0;k<bps;k++) {
+            /* initialize */
+            num[k] = -1000000;
+            den[k] = -1000000;
+        }
+
+        for (i=0;i<M;i++) {
+            metric =  symbol_likelihood[n*M+i]; /* channel metric for this symbol */
+
+            mask = 1 << (bps - 1);
+            for (j=0;j<bps;j++) {
+                mask = mask >> 1;
+            }
+            mask = 1 << (bps - 1);
+
+            for (k=0;k<bps;k++) {       /* loop over bits */
+                if (mask&i) {
+                    /* this bit is a one */
+                    num[k] = max_star0( num[k], metric );
+                } else {
+                    /* this bit is a zero */
+                    den[k] = max_star0( den[k], metric );
+                }
+                mask = mask >> 1;
+            }
+        }
+        for (k=0;k<bps;k++) {
+            bit_likelihood[bps*n+k] = num[k] - den[k];
+        }
+    }
+}
+
+
+void symbols_to_llrs(float llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, float mean_amp, int nsyms) {
+    int i;
+
+    float symbol_likelihood[nsyms*QPSK_CONSTELLATION_SIZE];
+    float bit_likelihood[nsyms*QPSK_BITS_PER_SYMBOL];
+
+    Demod2D(symbol_likelihood, rx_qpsk_symbols, S_matrix, EsNo, rx_amps, mean_amp, nsyms);
+    Somap(bit_likelihood, symbol_likelihood, nsyms);
+    for(i=0; i<nsyms*QPSK_BITS_PER_SYMBOL; i++) {
+        llr[i] = -bit_likelihood[i];
+    }
+}
+
+void ldpc_print_info(struct LDPC *ldpc) {
+fprintf(stderr, "ldpc->max_iter = %d\n", ldpc->max_iter);
+fprintf(stderr, "ldpc->dec_type = %d\n", ldpc->dec_type);
+fprintf(stderr, "ldpc->q_scale_factor = %d\n", ldpc->q_scale_factor);
+fprintf(stderr, "ldpc->r_scale_factor = %d\n", ldpc->r_scale_factor);
+fprintf(stderr, "ldpc->CodeLength = %d\n", ldpc->CodeLength);
+fprintf(stderr, "ldpc->NumberParityBits = %d\n", ldpc->NumberParityBits);
+fprintf(stderr, "ldpc->NumberRowsHcols = %d\n", ldpc->NumberRowsHcols);
+fprintf(stderr, "ldpc->max_row_weight = %d\n", ldpc->max_row_weight);
+fprintf(stderr, "ldpc->max_col_weight = %d\n", ldpc->max_col_weight);
+fprintf(stderr, "ldpc->data_bits_per_frame = %d\n", ldpc->data_bits_per_frame);
+fprintf(stderr, "ldpc->coded_bits_per_frame = %d\n", ldpc->coded_bits_per_frame);
+fprintf(stderr, "ldpc->coded_syms_per_frame = %d\n", ldpc->coded_syms_per_frame);
+}
+
+/* vi:set ts=4 et sts=4: */

src/mpdecode_core.h

@@ -0,0 +1,54 @@
+/*
+  FILE...: mpdecode_core.h
+  AUTHOR.: David Rowe
+  CREATED: Sep 2016
+
+  C-callable core functions for MpDecode, so they can be used for
+  Octave and C programs.  Also some convenience functions to help use
+  the C-callable LDPC decoder in C programs.
+*/
+
+#ifndef __MPDECODE_CORE__
+#define __MPDECODE_CORE__
+
+#include <stdint.h>
+
+#include "comp.h"
+
+struct LDPC {
+    int max_iter;
+    int dec_type;
+    int q_scale_factor;
+    int r_scale_factor;
+    int CodeLength;
+    int NumberParityBits;
+    int NumberRowsHcols;
+    int max_row_weight;
+    int max_col_weight;
+
+    /* these two are fixed to code params */
+    int ldpc_data_bits_per_frame;
+    int ldpc_coded_bits_per_frame;
+
+    /* these three may vary if we don't use all data bits in code */
+    int data_bits_per_frame;
+    int coded_bits_per_frame;
+    int coded_syms_per_frame;
+
+    uint16_t *H_rows;
+    uint16_t *H_cols;
+};
+
+void encode(struct LDPC *ldpc, unsigned char ibits[], unsigned char pbits[]);
+
+int run_ldpc_decoder(struct LDPC *ldpc, uint8_t out_char[], float input[], int *parityCheckCount);
+
+void sd_to_llr(float llr[], double sd[], int n);
+void Demod2D(float symbol_likelihood[], COMP r[], COMP S_matrix[], float EsNo, float fading[], float mean_amp, int number_symbols);
+void Somap(float bit_likelihood[], float symbol_likelihood[], int number_symbols);
+void symbols_to_llrs(float llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, float mean_amp, int nsyms);
+
+void ldpc_print_info(struct LDPC *ldpc);
+
+
+#endif

src/mpdecode_core_test.c

@@ -0,0 +1,693 @@
+/*
+  FILE...: mpdecode_core.c
+  AUTHOR.: Matthew C. Valenti, Rohit Iyer Seshadri, David Rowe
+  CREATED: Sep 2016
+
+  C-callable core functions moved from MpDecode.c, so they can be used for
+  Octave and C programs.
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <assert.h>
+
+#include "mpdecode_core_test.h"
+#ifndef USE_ORIGINAL_PHI0
+#include "phi0.h"
+#endif
+
+#include "debug_alloc.h"
+
+#ifdef __EMBEDDED__
+#include "machdep.h"
+#endif
+
+#define QPSK_CONSTELLATION_SIZE 4
+#define QPSK_BITS_PER_SYMBOL    2
+
+/* QPSK constellation for symbol likelihood calculations */
+
+static COMP S_matrix[] = {
+    { 1.0f,  0.0f},
+    { 0.0f,  1.0f},
+    { 0.0f, -1.0f},
+    {-1.0f,  0.0f}
+};
+
+// c_nodes will be an array of NumberParityBits of struct c_node
+// Each c_node contains an array of <degree> c_sub_node elements
+// This structure reduces the indexing caluclations in SumProduct()
+
+struct c_sub_node { // Order is important here to keep total size small.
+  uint16_t index;   // Values from H_rows (except last 2 entries)
+  uint16_t socket;  // The socket number at the v_node
+  float    message; // modified during operation!
+};
+
+struct c_node {
+  int degree;       // A count of elements in the following arrays
+  struct c_sub_node *subs;
+};
+
+// v_nodes will be an array of CodeLength of struct v_node
+
+struct v_sub_node {
+  uint16_t index;  //    the index of a c_node it is connected to
+                   //    Filled with values from H_cols (except last 2 entries)
+  uint16_t socket; //    socket number at the c_node
+  float message;   //    Loaded with input data
+                   //    modified during operation!
+  uint8_t sign;    //    1 if input is negative
+                   //    modified during operation!
+};
+
+struct v_node {
+  int degree;       // A count of ???
+  float initial_value;
+  struct v_sub_node *subs;
+};
+
+void encode(struct LDPC *ldpc, unsigned char ibits[], unsigned char pbits[]) {
+    unsigned int p, i, tmp, par, prev=0;
+    int          ind;
+    uint16_t     *H_rows = ldpc->H_rows;
+
+    for (p=0; p<ldpc->NumberParityBits; p++) {
+        par = 0;
+
+        for (i=0; i<ldpc->max_row_weight; i++) {
+            ind = H_rows[p + i*ldpc->NumberParityBits];
+            par = par + ibits[ind-1];
+        }
+
+        tmp = par + prev;
+
+        tmp &= 1;    // only retain the lsb
+        prev = tmp;
+        pbits[p] = tmp;
+    }
+}
+
+#ifdef USE_ORIGINAL_PHI0
+/* Phi function */
+static float phi0(
+                  float x )
+{
+  float z;
+
+  if (x>10)
+    return( 0 );
+  else if (x< 9.08e-5 )
+    return( 10 );
+  else if (x > 9)
+    return( 1.6881e-4 );
+  /* return( 1.4970e-004 ); */
+  else if (x > 8)
+    return( 4.5887e-4 );
+  /* return( 4.0694e-004 ); */
+  else if (x > 7)
+    return( 1.2473e-3 );
+  /* return( 1.1062e-003 ); */
+  else if (x > 6)
+    return( 3.3906e-3 );
+  /* return( 3.0069e-003 ); */
+  else if (x > 5)
+    return( 9.2168e-3 );
+  /* return( 8.1736e-003 ); */
+  else {
+    z = (float) exp(x);
+    return( (float) log( (z+1)/(z-1) ) );
+  }
+}
+#endif
+
+
+/* Values for linear approximation (DecoderType=5) */
+
+
+#define AJIAN -0.24904163195436
+#define TJIAN 2.50681740420944
+
+/* The linear-log-MAP algorithm */
+
+static float max_star0(
+                       float delta1,
+                       float delta2 )
+{
+    register float diff;
+
+    diff = delta2 - delta1;
+
+    if ( diff > TJIAN )
+        return( delta2 );
+    else if ( diff < -TJIAN )
+        return( delta1 );
+    else if ( diff > 0 )
+        return( delta2 + AJIAN*(diff-TJIAN) );
+    else
+        return( delta1 - AJIAN*(diff+TJIAN) );
+}
+
+void init_c_v_nodes(struct c_node *c_nodes,
+                    int     shift,
+                    int     NumberParityBits,
+                    int     max_row_weight,
+                    uint16_t *H_rows,
+                    int     H1,
+                    int     CodeLength,
+                    struct v_node *v_nodes,
+                    int     NumberRowsHcols,
+                    uint16_t *H_cols,
+                    int     max_col_weight,
+                    int     dec_type,
+                    float  *input)
+{
+    int i, j, k, count, cnt, c_index, v_index;
+
+    /* first determine the degree of each c-node */
+
+    if (shift ==0){
+        for (i=0;i<NumberParityBits;i++) {
+            count = 0;
+            for (j=0;j<max_row_weight;j++) {
+                if ( H_rows[i+j*NumberParityBits] > 0 ) {
+                    count++;
+                }
+            }
+            c_nodes[i].degree = count;
+            if (H1){
+                if (i==0){
+                    c_nodes[i].degree=count+1;
+                }
+                else{
+                    c_nodes[i].degree=count+2;
+                }
+            }
+        }
+    }
+    else{
+        cnt=0;
+        for (i=0;i<(NumberParityBits/shift);i++) {
+            for (k=0;k<shift;k++){
+                count = 0;
+                for (j=0;j<max_row_weight;j++) {
+                    if ( H_rows[cnt+j*NumberParityBits] > 0 ) {
+                        count++;
+                    }
+                }
+                c_nodes[cnt].degree = count;
+                if ((i==0)||(i==(NumberParityBits/shift)-1)){
+                    c_nodes[cnt].degree=count+1;
+                }
+                else{
+                    c_nodes[cnt].degree=count+2;
+                }
+                cnt++;
+            }
+        }
+    }
+
+    if (H1) {
+
+        if (shift ==0){
+            for (i=0;i<NumberParityBits;i++) {
+
+                // Allocate sub nodes
+                c_nodes[i].subs = CALLOC(c_nodes[i].degree, sizeof(struct c_sub_node));
+                assert(c_nodes[i].subs);
+
+	        // Populate sub nodes
+                for (j=0;j<c_nodes[i].degree-2;j++) {
+                    c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
+                }
+                j=c_nodes[i].degree-2;
+
+                if (i==0){
+                    c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
+                }
+                else {
+                    c_nodes[i].subs[j].index = (CodeLength-NumberParityBits)+i-1;
+                }
+
+                j=c_nodes[i].degree-1;
+                c_nodes[i].subs[j].index = (CodeLength-NumberParityBits)+i;
+
+            }
+        }
+        if (shift >0){
+            cnt=0;
+            for (i=0;i<(NumberParityBits/shift);i++){
+
+                for (k =0;k<shift;k++){
+
+                    // Allocate sub nodes
+                    c_nodes[cnt].subs = CALLOC(c_nodes[cnt].degree, sizeof(struct c_sub_node));
+                    assert(c_nodes[cnt].subs);
+
+	            // Populate sub nodes
+                    for (j=0;j<c_nodes[cnt].degree-2;j++) {
+                        c_nodes[cnt].subs[j].index = (H_rows[cnt+j*NumberParityBits] - 1);
+                    }
+                    j=c_nodes[cnt].degree-2;
+                    if ((i ==0)||(i==(NumberParityBits/shift-1))){
+                        c_nodes[cnt].subs[j].index = (H_rows[cnt+j*NumberParityBits] - 1);
+                    }
+                    else{
+                        c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i);
+                    }
+                    j=c_nodes[cnt].degree-1;
+                    c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i+1);
+                    if (i== (NumberParityBits/shift-1))
+                        {
+                            c_nodes[cnt].subs[j].index = (CodeLength-NumberParityBits)+k+shift*(i);
+                        }
+                    cnt++;
+                }
+            }
+        }
+
+    } else {
+        for (i=0;i<NumberParityBits;i++) {
+
+            // Allocate sub nodes
+            c_nodes[i].subs = CALLOC(c_nodes[i].degree, sizeof(struct c_sub_node));
+            assert(c_nodes[i].subs);
+
+	    // Populate sub nodes
+            for (j=0;j<c_nodes[i].degree;j++){
+                c_nodes[i].subs[j].index = (H_rows[i+j*NumberParityBits] - 1);
+
+            }
+        }
+    }
+
+
+    /* determine degree of each v-node */
+
+    for(i=0;i<(CodeLength-NumberParityBits+shift);i++){
+        count=0;
+        for (j=0;j<max_col_weight;j++) {
+            if ( H_cols[i+j*NumberRowsHcols] > 0 ) {
+                count++;
+            }
+        }
+        v_nodes[i].degree = count;
+    }
+
+    for(i=CodeLength-NumberParityBits+shift;i<CodeLength;i++){
+        count=0;
+        if (H1){
+            if(i!=CodeLength-1){
+                v_nodes[i].degree=2;
+            }  else{
+                v_nodes[i].degree=1;
+            }
+
+        } else{
+            for (j=0;j<max_col_weight;j++) {
+                if ( H_cols[i+j*NumberRowsHcols] > 0 ) {
+                    count++;
+                }
+            }
+            v_nodes[i].degree = count;
+        }
+    }
+
+    if (shift>0){
+        v_nodes[CodeLength-1].degree =v_nodes[CodeLength-1].degree+1;
+    }
+
+
+    /* set up v_nodes */
+
+    for (i=0;i<CodeLength;i++) {
+
+        // Allocate sub nodes
+        v_nodes[i].subs = CALLOC(v_nodes[i].degree, sizeof(struct v_sub_node));
+        assert(v_nodes[i].subs);
+
+	// Populate sub nodes
+
+        /* index tells which c-nodes this v-node is connected to */
+        v_nodes[i].initial_value = input[i];
+        count=0;
+
+        for (j=0;j<v_nodes[i].degree;j++) {
+            if ((H1)&& (i>=CodeLength-NumberParityBits+shift)){
+                v_nodes[i].subs[j].index=i-(CodeLength-NumberParityBits+shift)+count;
+                if (shift ==0){
+                    count=count+1;
+                }
+                else{
+                    count=count+shift;
+                }
+            } else  {
+                v_nodes[i].subs[j].index = (H_cols[i+j*NumberRowsHcols] - 1);
+            }
+
+            /* search the connected c-node for the proper message value */
+            for (c_index=0;c_index<c_nodes[ v_nodes[i].subs[j].index ].degree;c_index++)
+                if ( c_nodes[ v_nodes[i].subs[j].index ].subs[c_index].index == i ) {
+                    v_nodes[i].subs[j].socket = c_index;
+                    break;
+                }
+            /* initialize v-node with received LLR */
+            if ( dec_type == 1)
+                v_nodes[i].subs[j].message = fabs(input[i]);
+            else
+                v_nodes[i].subs[j].message = phi0( fabs(input[i]) );
+
+            if (input[i] < 0)
+                v_nodes[i].subs[j].sign = 1;
+        }
+
+    }
+
+    /* now finish setting up the c_nodes */
+    for (i=0;i<NumberParityBits;i++) {
+        /* index tells which v-nodes this c-node is connected to */
+        for (j=0;j<c_nodes[i].degree;j++) {
+            /* search the connected v-node for the proper message value */
+            for (v_index=0;v_index<v_nodes[ c_nodes[i].subs[j].index ].degree;v_index++)
+                if (v_nodes[ c_nodes[i].subs[j].index ].subs[v_index].index == i ) {
+                    c_nodes[i].subs[j].socket = v_index;
+                    break;
+                }
+        }
+    }
+
+}
+
+
+///////////////////////////////////////
+/* function for doing the MP decoding */
+// Returns the iteration count
+int SumProduct(int	  *parityCheckCount,
+			 char     DecodedBits[],
+			 struct c_node c_nodes[],
+			 struct v_node v_nodes[],
+			 int	  CodeLength,
+			 int	  NumberParityBits,
+			 int	  max_iter,
+			 float    r_scale_factor,
+			 float    q_scale_factor,
+			 int      data[] )
+{
+  int result;
+  int bitErrors;
+  int i,j, iter;
+  float phi_sum;
+  int sign;
+  float temp_sum;
+  float Qi;
+  int   ssum;
+
+  result = max_iter;
+  for (iter=0;iter<max_iter;iter++) {
+    for(i=0; i<CodeLength; i++) DecodedBits[i] = 0; // Clear each pass!
+    bitErrors = 0;
+    /* update r */
+    ssum = 0;
+    for (j=0;j<NumberParityBits;j++) {
+      sign = v_nodes[ c_nodes[j].subs[0].index ].subs[ c_nodes[j].subs[0].socket ].sign;
+      phi_sum = v_nodes[ c_nodes[j].subs[0].index ].subs[ c_nodes[j].subs[0].socket ].message;
+
+      for (i=1;i<c_nodes[j].degree;i++) {
+        // Compiler should optomize this but write the best we can to start from.
+        struct c_sub_node *cp = &c_nodes[j].subs[i];
+        struct v_sub_node *vp = &v_nodes[ cp->index ].subs[ cp->socket ];
+	    phi_sum += vp->message;
+	    sign ^= vp->sign;
+      }
+
+      if (sign==0) ssum++;
+
+//fprintf(stderr, "  up-r: %d: sign=%d, phi_sum=%f\n", j, sign, (double)phi_sum);
+
+      for (i=0;i<c_nodes[j].degree;i++) {
+        struct c_sub_node *cp = &c_nodes[j].subs[i];
+        struct v_sub_node *vp = &v_nodes[ cp->index ].subs[ cp->socket ];
+	    if ( sign ^ vp->sign ) {
+	      cp->message = -phi0( phi_sum - vp->message ); // *r_scale_factor;
+	    } else
+	      cp->message =  phi0( phi_sum - vp->message ); // *r_scale_factor;
+      }
+    }
+
+    /* update q */
+//#ifdef __EMBEDDED__
+//PROFILE_SAMPLE_AND_LOG(ldpc_SP_upq, ldpc_SP_upr, "ldpc_SP_update_r");
+//#endif
+    for (i=0;i<CodeLength;i++) {
+
+      /* first compute the LLR */
+      Qi = v_nodes[i].initial_value;
+      for (j=0;j<v_nodes[i].degree;j++) {
+        struct v_sub_node *vp = &v_nodes[i].subs[j];
+	    Qi += c_nodes[ vp->index ].subs[ vp->socket ].message;
+      }
+
+      /* make hard decision */
+      if (Qi < 0) {
+	    DecodedBits[i] = 1;
+      }
+
+      /* now subtract to get the extrinsic information */
+      for (j=0;j<v_nodes[i].degree;j++) {
+        struct v_sub_node *vp = &v_nodes[i].subs[j];
+	    temp_sum = Qi - c_nodes[ vp->index ].subs[ vp->socket ].message;
+
+	    vp->message = phi0( fabs( temp_sum ) ); // *q_scale_factor;
+	    if (temp_sum > 0)
+	      vp->sign = 0;
+	    else
+	      vp->sign = 1;
+      }
+    }
+//#ifdef __EMBEDDED__
+//PROFILE_SAMPLE_AND_LOG(ldpc_SP_misc, ldpc_SP_upq, "ldpc_SP_update_q");
+//#endif
+
+    /* count data bit errors, assuming that it is systematic */
+    for (i=0;i<CodeLength-NumberParityBits;i++)
+      if ( DecodedBits[i] != data[i] )
+	    bitErrors++;
+
+    /* Halt if zero errors */
+    if (bitErrors == 0) {
+      result = iter + 1;
+      break;
+    }
+
+    // count the number of PC satisfied and exit if all OK
+    *parityCheckCount = ssum;
+    if (ssum==NumberParityBits)  {
+      result = iter + 1;
+      break;
+    }
+
+//#ifdef __EMBEDDED__
+//PROFILE_SAMPLE_AND_LOG2(ldpc_SP_misc, "ldpc_SP_misc");
+//PROFILE_SAMPLE_AND_LOG2(ldpc_SP_iter, "ldpc_SP_iter");
+//#endif
+  }
+
+return(result);
+}
+
+/* Convenience function to call LDPC decoder from C programs */
+
+int run_ldpc_decoder(struct LDPC *ldpc, uint8_t out_char[], float input[], int *parityCheckCount)
+    int		max_iter;
+    float       q_scale_factor, r_scale_factor;
+    int         CodeLength, NumberParityBits, NumberRowsHcols, shift, H1;
+    int         i;
+    struct c_node *c_nodes;
+    struct v_node *v_nodes;
+
+    #ifdef __EMBEDDED__
+    PROFILE_VAR(ldpc_init, ldpc_SP);
+    PROFILE_SAMPLE(ldpc_init);
+    #endif
+
+    /* default values */
+    max_iter  = ldpc->max_iter;
+    q_scale_factor = ldpc->q_scale_factor;
+    r_scale_factor = ldpc->r_scale_factor;
+
+    CodeLength = ldpc->CodeLength;         /* length of entire codeword */
+    NumberParityBits = ldpc->NumberParityBits;
+    NumberRowsHcols = ldpc->NumberRowsHcols;
+
+    char *DecodedBits = CALLOC( CodeLength, sizeof( char ) );
+    assert(DecodedBits);
+
+    /* derive some parameters */
+
+    shift = (NumberParityBits + NumberRowsHcols) - CodeLength;
+    if (NumberRowsHcols == CodeLength) {
+        H1=0;
+        shift=0;
+    } else {
+        H1=1;
+    }
+
+    /* initialize c-node and v-node structures */
+    c_nodes = CALLOC( NumberParityBits, sizeof( struct c_node ) );
+    assert(c_nodes);
+    v_nodes = CALLOC( CodeLength, sizeof( struct v_node));
+    assert(v_nodes);
+
+    init_c_v_nodes( c_nodes, shift,
+        NumberParityBits, ldpc->max_row_weight, ldpc->H_rows, H1, CodeLength,
+        v_nodes, NumberRowsHcols, ldpc->H_cols, ldpc->max_col_weight, ldpc->dec_type,
+        input);
+
+    int DataLength = CodeLength - NumberParityBits;
+    int *data_int = CALLOC( DataLength, sizeof(int) );
+
+    /* Call function to do the actual decoding */
+
+    int iter = SumProduct(parityCheckCount, DecodedBits, c_nodes, v_nodes, 
+                CodeLength, NumberParityBits, max_iter, 
+                r_scale_factor, q_scale_factor, data_int);
+
+    for (i=0; i<CodeLength; i++) out_char[i] = DecodedBits[i];
+
+    /* Clean up memory */
+    FREE(DecodedBits);
+    FREE( data_int );
+    for (i=0;i<NumberParityBits;i++) FREE( c_nodes[i].subs );
+    FREE( c_nodes );
+    for (i=0;i<CodeLength;i++) FREE( v_nodes[i].subs);
+    FREE( v_nodes );
+
+    return iter;
+}
+
+
+void sd_to_llr(float llr[], double sd[], int n) {
+    double sum, mean, sign, sumsq, estvar, estEsN0, x;
+    int i;
+
+    /* convert SD samples to LLRs -------------------------------*/
+
+    sum = 0.0;
+    for(i=0; i<n; i++)
+        sum += fabs(sd[i]);
+    mean = sum/n;
+
+    /* find variance from +/-1 symbol position */
+
+    sum = sumsq = 0.0;
+    for(i=0; i<n; i++) {
+        sign = (sd[i] > 0.0L) - (sd[i] < 0.0L);
+        x = (sd[i]/mean - sign);
+        sum += x;
+        sumsq += x*x;
+    }
+    estvar = (n * sumsq - sum * sum) / (n * (n - 1));
+    //fprintf(stderr, "mean: %f var: %f\n", mean, estvar);
+
+    estEsN0 = 1.0/(2.0L * estvar + 1E-3);
+    for(i=0; i<n; i++)
+        llr[i] = 4.0L * estEsN0 * sd[i];
+}
+
+
+/*
+   Determine symbol likelihood from received QPSK symbols.
+
+   Notes:
+
+   1) We assume fading[] is real, it is also possible to compute
+      with complex fading, see CML library Demod2D.c source code.
+   2) Using floats instead of doubles, for stm32.
+      Testing shows good BERs with floats.
+*/
+
+void Demod2D(float   symbol_likelihood[],       /* output, M*number_symbols              */
+             COMP    r[],                       /* received QPSK symbols, number_symbols */
+             COMP    S_matrix[],                /* constellation of size M               */
+             float   EsNo,
+             float   fading[],                  /* real fading values, number_symbols    */
+             float   mean_amp,
+             int     number_symbols)
+{
+    int     M=QPSK_CONSTELLATION_SIZE;
+    int     i,j;
+    float  tempsr, tempsi, Er, Ei;
+
+    /* determine output */
+
+    for (i=0;i<number_symbols;i++) {                /* go through each received symbol */
+        for (j=0;j<M;j++) {                         /* each postulated symbol          */
+            tempsr = fading[i]*S_matrix[j].real/mean_amp;
+            tempsi = fading[i]*S_matrix[j].imag/mean_amp;
+            Er = r[i].real/mean_amp - tempsr;
+            Ei = r[i].imag/mean_amp - tempsi;
+            symbol_likelihood[i*M+j] = -EsNo*(Er*Er+Ei*Ei);
+            //printf("symbol_likelihood[%d][%d] = %f\n", i,j,symbol_likelihood[i*M+j]);
+        }
+        //exit(0);
+    }
+
+}
+
+
+void Somap(float  bit_likelihood[],      /* number_bits, bps*number_symbols */
+           float  symbol_likelihood[],   /* M*number_symbols                */
+           int     number_symbols)
+{
+    int    M=QPSK_CONSTELLATION_SIZE, bps = QPSK_BITS_PER_SYMBOL;
+    int    n,i,j,k,mask;
+    float num[bps], den[bps];
+    float metric;
+
+    for (n=0; n<number_symbols; n++) { /* loop over symbols */
+        for (k=0;k<bps;k++) {
+            /* initialize */
+            num[k] = -1000000;
+            den[k] = -1000000;
+        }
+
+        for (i=0;i<M;i++) {
+            metric =  symbol_likelihood[n*M+i]; /* channel metric for this symbol */
+
+            mask = 1 << (bps - 1);
+            for (j=0;j<bps;j++) {
+                mask = mask >> 1;
+            }
+            mask = 1 << (bps - 1);
+
+            for (k=0;k<bps;k++) {	/* loop over bits */
+                if (mask&i) {
+                    /* this bit is a one */
+                    num[k] = max_star0( num[k], metric );
+                } else {
+                    /* this bit is a zero */
+                    den[k] = max_star0( den[k], metric );
+                }
+                mask = mask >> 1;
+            }
+        }
+        for (k=0;k<bps;k++) {
+            bit_likelihood[bps*n+k] = num[k] - den[k];
+        }
+    }
+}
+
+void symbols_to_llrs(float llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, float mean_amp, int nsyms) {
+    int i;
+
+    float symbol_likelihood[nsyms*QPSK_CONSTELLATION_SIZE];
+    float bit_likelihood[nsyms*QPSK_BITS_PER_SYMBOL];
+
+    Demod2D(symbol_likelihood, rx_qpsk_symbols, S_matrix, EsNo, rx_amps, mean_amp, nsyms);
+    Somap(bit_likelihood, symbol_likelihood, nsyms);
+    for(i=0; i<nsyms*QPSK_BITS_PER_SYMBOL; i++) {
+        llr[i] = -bit_likelihood[i];
+    }
+}

src/mpdecode_core_test.h

@@ -0,0 +1,49 @@
+/*
+  FILE...: mpdecode_core.h
+  AUTHOR.: David Rowe
+  CREATED: Sep 2016
+
+  C-callable core functions for MpDecode, so they can be used for
+  Octave and C programs.  Also some convenience functions to help use
+  the C-callable LDPC decoder in C programs.
+*/
+
+#ifndef __MPDECODE_CORE__
+#define __MPDECODE_CORE__
+
+#include <stdint.h>
+
+#include "comp.h"
+
+struct LDPC {
+    int max_iter;
+    int dec_type;
+    int q_scale_factor;
+    int r_scale_factor;
+    int CodeLength;
+    int NumberParityBits;
+    int NumberRowsHcols;
+    int max_row_weight;
+    int max_col_weight;
+    int data_bits_per_frame;
+    int coded_bits_per_frame;
+    int coded_syms_per_frame;
+    uint16_t *H_rows;
+    uint16_t *H_cols;
+};
+
+extern void ldpc_init(struct LDPC *ldpc, int *size_common);
+extern void ldpc_free_mem(struct LDPC *ldpc);
+
+extern void encode(struct LDPC *ldpc, unsigned char ibits[], unsigned char pbits[]);
+
+int run_ldpc_decoder(struct LDPC *ldpc, uint8_t out_char[], float input[], int *parityCheckCount);
+extern void ldpc_dump_nodes(struct LDPC *ldpc);
+
+extern void sd_to_llr(float llr[], double sd[], int n);
+
+extern void Demod2D(float symbol_likelihood[], COMP r[], COMP S_matrix[], float EsNo, float fading[], float mean_amp, int number_symbols);
+extern void Somap(float bit_likelihood[], float symbol_likelihood[], int number_symbols);
+extern void symbols_to_llrs(float llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, float mean_amp, int nsyms);
+
+#endif

src/octave.c

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: octave.c
+  AUTHOR......: David Rowe
+  DATE CREATED: April 28 2012
+
+  Functions to save C arrays in GNU Octave matrix format.  The output text
+  file can be directly read into Octave using "load filename".
+
+\*---------------------------------------------------------------------------*/
+
+
+/*
+  Copyright (C) 2012 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <stdio.h>
+#include <stdarg.h>
+
+#include "octave.h"
+
+#ifdef ARM_MATH_CM4
+#include "Trace.h"
+#endif
+
+#define OCTAVE_BUFSIZE 2048
+
+
+void flush_buffer(FILE* f, char* buffer,size_t* buf_idx_ptr)
+{
+#ifdef ARM_MATH_CM4
+    trace_write(buffer,*buf_idx_ptr);
+#else
+    fwrite(buffer,*buf_idx_ptr,1,f);
+#endif
+    *buf_idx_ptr = 0;
+}
+
+void handle_buffer(FILE* f, char* buffer,const size_t max_buf, size_t* buf_idx_ptr, size_t l)
+{
+    *buf_idx_ptr += l;
+    if (*buf_idx_ptr > max_buf - 64)
+    {
+        flush_buffer(f, buffer,buf_idx_ptr);
+    }
+}
+
+signed int printf_buffer(FILE* f, char* buffer,const size_t max_buf, size_t* buf_idx_ptr, const char *pFormat, ...)
+{
+    va_list    ap;
+    signed int rc;
+
+    va_start(ap, pFormat);
+    rc = vsnprintf(&buffer[*buf_idx_ptr], max_buf - *buf_idx_ptr, pFormat, ap);
+    va_end(ap);
+    if (rc>0)
+    {
+        handle_buffer(f, buffer,max_buf,buf_idx_ptr,rc);
+    }
+    return rc;
+}
+
+
+void printf_header(FILE* f, char* buffer,const size_t max_buf, size_t* buf_idx_ptr, const char *name, const char *dtype, int rows, int cols, int isFloat)
+{
+#ifdef ARM_MATH_CM4
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, buf_idx_ptr, "# hex: %s\n", isFloat?"true":"false");
+#endif    
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, buf_idx_ptr, "# name: %s\n", name);
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, buf_idx_ptr, "# type: %s\n",dtype);
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, buf_idx_ptr, "# rows: %d\n", rows);
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, buf_idx_ptr, "# columns: %d\n", cols);
+}
+void octave_save_int(FILE *f, char name[], int data[], int rows, int cols)
+{
+    int r,c;
+    char buffer[OCTAVE_BUFSIZE];
+    size_t buf_idx = 0;
+
+    printf_header(f, buffer, OCTAVE_BUFSIZE, &buf_idx, name, "matrix", rows, cols, 0);
+
+    for(r=0; r<rows; r++) {
+        for(c=0; c<cols; c++)
+            printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, " %d", data[r*cols+c]);
+        printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, "\n");
+    }
+
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, "\n\n");
+    flush_buffer(f, buffer, &buf_idx);
+}
+
+void octave_save_float(FILE *f, char name[], float data[], int rows, int cols, int col_len)
+{
+    int r,c;
+    char buffer[OCTAVE_BUFSIZE];
+    size_t buf_idx = 0;
+
+    printf_header(f, buffer, OCTAVE_BUFSIZE, &buf_idx, name, "matrix", rows, cols, 1);
+
+    for(r=0; r<rows; r++) {
+        for(c=0; c<cols; c++)
+            printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, " %f", data[r*col_len+c]);
+        printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, "\n");
+    }
+
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, "\n\n");
+    flush_buffer(f, buffer, &buf_idx);
+}
+
+
+void octave_save_complex(FILE *f, char name[], COMP data[], int rows, int cols, int col_len)
+{
+    int r,c;
+    char buffer[OCTAVE_BUFSIZE];
+    size_t buf_idx = 0;
+
+
+    printf_header(f, buffer, OCTAVE_BUFSIZE, &buf_idx, name, "complex matrix", rows, cols, 1);
+
+    for(r=0; r<rows; r++) {
+
+        for(c=0; c<cols; c++)
+        {
+            printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, " (%f,%f)", data[r*col_len+c].real, data[r*col_len+c].imag);
+        }
+        printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, "\n");
+    }
+    printf_buffer(f, buffer, OCTAVE_BUFSIZE, &buf_idx, "\n\n");
+    flush_buffer(f, buffer, &buf_idx);
+}

src/octave.h

@@ -0,0 +1,39 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: octave.h
+  AUTHOR......: David Rowe
+  DATE CREATED: April 28 2012
+
+  Functions to save C arrays in Octave matrix format.  the output text
+  file can be directly read into octave using "load filename".
+
+\*---------------------------------------------------------------------------*/
+
+
+/*
+  Copyright (C) 2012 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program 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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __OCTAVE__
+#define __OCTAVE__
+
+#include "comp.h"
+
+void octave_save_int(FILE *f, char name[], int data[], int rows, int cols);
+void octave_save_float(FILE *f, char name[], float data[], int rows, int cols, int col_len);
+void octave_save_complex(FILE *f, char name[], COMP data[], int rows, int cols, int col_len);
+
+#endif

src/phi0.c

@@ -0,0 +1,218 @@
+
+// phi0.c
+//
+// An approximation of the function
+//
+// This file is generated by the gen_phi0 scritps
+// Any changes should be made to that file, not this one
+
+#include <stdint.h>
+
+#define SI16(f) ((int32_t)(f * (1<<16)))
+
+float phi0( float xf ) {
+
+    int32_t x = SI16(xf);
+
+    if (x >= SI16(10.0f)) return(0.0f);
+    else {
+      if (x >= SI16(5.0f)) {
+        int i = 19 - (x >> 15);
+        switch (i) {
+          case 0: return(0.000116589f); // (9.5)
+          case 1: return(0.000192223f); // (9.0)
+          case 2: return(0.000316923f); // (8.5)
+          case 3: return(0.000522517f); // (8.0)
+          case 4: return(0.000861485f); // (7.5)
+          case 5: return(0.001420349f); // (7.0)
+          case 6: return(0.002341760f); // (6.5)
+          case 7: return(0.003860913f); // (6.0)
+          case 8: return(0.006365583f); // (5.5)
+          case 9: return(0.010495133f); // (5.0)
+        }
+      }
+      else {
+        if (x >= SI16(1.0f)) {
+          int i = 79 - (x >> 12);
+          switch (i) {
+            case 0: return(0.013903889f); // (4.9375)
+            case 1: return(0.014800644f); // (4.8750)
+            case 2: return(0.015755242f); // (4.8125)
+            case 3: return(0.016771414f); // (4.7500)
+            case 4: return(0.017853133f); // (4.6875)
+            case 5: return(0.019004629f); // (4.6250)
+            case 6: return(0.020230403f); // (4.5625)
+            case 7: return(0.021535250f); // (4.5000)
+            case 8: return(0.022924272f); // (4.4375)
+            case 9: return(0.024402903f); // (4.3750)
+            case 10: return(0.025976926f); // (4.3125)
+            case 11: return(0.027652501f); // (4.2500)
+            case 12: return(0.029436184f); // (4.1875)
+            case 13: return(0.031334956f); // (4.1250)
+            case 14: return(0.033356250f); // (4.0625)
+            case 15: return(0.035507982f); // (4.0000)
+            case 16: return(0.037798579f); // (3.9375)
+            case 17: return(0.040237016f); // (3.8750)
+            case 18: return(0.042832850f); // (3.8125)
+            case 19: return(0.045596260f); // (3.7500)
+            case 20: return(0.048538086f); // (3.6875)
+            case 21: return(0.051669874f); // (3.6250)
+            case 22: return(0.055003924f); // (3.5625)
+            case 23: return(0.058553339f); // (3.5000)
+            case 24: return(0.062332076f); // (3.4375)
+            case 25: return(0.066355011f); // (3.3750)
+            case 26: return(0.070637993f); // (3.3125)
+            case 27: return(0.075197917f); // (3.2500)
+            case 28: return(0.080052790f); // (3.1875)
+            case 29: return(0.085221814f); // (3.1250)
+            case 30: return(0.090725463f); // (3.0625)
+            case 31: return(0.096585578f); // (3.0000)
+            case 32: return(0.102825462f); // (2.9375)
+            case 33: return(0.109469985f); // (2.8750)
+            case 34: return(0.116545700f); // (2.8125)
+            case 35: return(0.124080967f); // (2.7500)
+            case 36: return(0.132106091f); // (2.6875)
+            case 37: return(0.140653466f); // (2.6250)
+            case 38: return(0.149757747f); // (2.5625)
+            case 39: return(0.159456024f); // (2.5000)
+            case 40: return(0.169788027f); // (2.4375)
+            case 41: return(0.180796343f); // (2.3750)
+            case 42: return(0.192526667f); // (2.3125)
+            case 43: return(0.205028078f); // (2.2500)
+            case 44: return(0.218353351f); // (2.1875)
+            case 45: return(0.232559308f); // (2.1250)
+            case 46: return(0.247707218f); // (2.0625)
+            case 47: return(0.263863255f); // (2.0000)
+            case 48: return(0.281099022f); // (1.9375)
+            case 49: return(0.299492155f); // (1.8750)
+            case 50: return(0.319127030f); // (1.8125)
+            case 51: return(0.340095582f); // (1.7500)
+            case 52: return(0.362498271f); // (1.6875)
+            case 53: return(0.386445235f); // (1.6250)
+            case 54: return(0.412057648f); // (1.5625)
+            case 55: return(0.439469363f); // (1.5000)
+            case 56: return(0.468828902f); // (1.4375)
+            case 57: return(0.500301872f); // (1.3750)
+            case 58: return(0.534073947f); // (1.3125)
+            case 59: return(0.570354566f); // (1.2500)
+            case 60: return(0.609381573f); // (1.1875)
+            case 61: return(0.651427083f); // (1.1250)
+            case 62: return(0.696805010f); // (1.0625)
+            case 63: return(0.745880827f); // (1.0000)
+          }
+        }
+        else {
+          if (x > SI16(0.007812f)) {
+            if (x > SI16(0.088388f)) {
+              if (x > SI16(0.250000f)) {
+                if (x > SI16(0.500000f)) {
+                  if (x > SI16(0.707107f)) {
+                    return(0.922449644f);
+                  } else {
+                    return(1.241248638f);
+                  }
+                } else {
+                  if (x > SI16(0.353553f)) {
+                    return(1.573515241f);
+                  } else {
+                    return(1.912825912f);
+                  }
+                }
+              } else {
+                if (x > SI16(0.125000f)) {
+                  if (x > SI16(0.176777f)) {
+                    return(2.255740095f);
+                  } else {
+                    return(2.600476919f);
+                  }
+                } else {
+                  return(2.946130351f);
+                }
+              }
+            } else {
+              if (x > SI16(0.022097f)) {
+                if (x > SI16(0.044194f)) {
+                  if (x > SI16(0.062500f)) {
+                    return(3.292243417f);
+                  } else {
+                    return(3.638586634f);
+                  }
+                } else {
+                  if (x > SI16(0.031250f)) {
+                    return(3.985045009f);
+                  } else {
+                    return(4.331560985f);
+                  }
+                }
+              } else {
+                if (x > SI16(0.011049f)) {
+                  if (x > SI16(0.015625f)) {
+                    return(4.678105767f);
+                  } else {
+                    return(5.024664952f);
+                  }
+                } else {
+                  return(5.371231340f);
+                }
+              }
+            }
+          } else {
+            if (x > SI16(0.000691f)) {
+              if (x > SI16(0.001953f)) {
+                if (x > SI16(0.003906f)) {
+                  if (x > SI16(0.005524f)) {
+                    return(5.717801329f);
+                  } else {
+                    return(6.064373119f);
+                  }
+                } else {
+                  if (x > SI16(0.002762f)) {
+                    return(6.410945809f);
+                  } else {
+                    return(6.757518949f);
+                  }
+                }
+              } else {
+                if (x > SI16(0.000977f)) {
+                  if (x > SI16(0.001381f)) {
+                    return(7.104092314f);
+                  } else {
+                    return(7.450665792f);
+                  }
+                } else {
+                  return(7.797239326f);
+                }
+              }
+            } else {
+              if (x > SI16(0.000173f)) {
+                if (x > SI16(0.000345f)) {
+                  if (x > SI16(0.000488f)) {
+                    return(8.143812888f);
+                  } else {
+                    return(8.490386464f);
+                  }
+                } else {
+                  if (x > SI16(0.000244f)) {
+                    return(8.836960047f);
+                  } else {
+                    return(9.183533634f);
+                  }
+                }
+              } else {
+                if (x > SI16(0.000086f)) {
+                  if (x > SI16(0.000122f)) {
+                    return(9.530107222f);
+                  } else {
+                    return(9.876680812f);
+                  }
+                } else {
+                  return(10.000000000f);
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+    return(10.0f);
+}

src/phi0.h

@@ -0,0 +1,7 @@
+// phi0.h
+#ifndef PHI0_H
+#define PHI0_H
+
+extern float phi0( float xf );
+
+#endif