updated ui added new features

This commit is contained in:
zach
2025-12-27 15:32:32 -07:00
parent 02ca7801ea
commit a2cfae3a22
589 changed files with 181780 additions and 569 deletions
+123
View File
@@ -0,0 +1,123 @@
add_definitions(-DFLOATING_POINT -DVAR_ARRAYS)
include_directories(../src)
add_executable(tfdmdv tfdmdv.c ../src/fdmdv.c ../src/kiss_fft.c ../src/octave.c)
target_link_libraries(tfdmdv codec2)
add_executable(tcohpsk tcohpsk.c ../src/cohpsk.c ../src/octave.c)
target_link_libraries(tcohpsk codec2)
add_executable(tfsk tfsk.c ../src/kiss_fft.c ../src/kiss_fftr.c ../src/octave.c ../src/modem_probe.c)
target_link_libraries(tfsk m)
add_executable(tfreedv_data_channel tfreedv_data_channel.c ../src/freedv_data_channel.c)
add_executable(tfmfsk tfmfsk.c ../src/octave.c ../src/modem_probe.c)
target_link_libraries(tfmfsk m)
add_definitions(-DMODEMPROBE_ENABLE -DXXXXX)
add_executable(tofdm tofdm.c ../src/octave.c)
target_link_libraries(tofdm m codec2)
add_executable(tofdm_acq tofdm_acq.c ../src/octave.c)
target_link_libraries(tofdm_acq m codec2)
if(UNIX) # Uses pthreads
add_executable(tfifo tfifo.c ../src/codec2_fifo.c)
target_link_libraries(tfifo codec2 ${CMAKE_THREAD_LIBS_INIT})
endif()
add_definitions(-D__UNITTEST__)
add_executable(tnewamp1 tnewamp1.c ../src/quantise.c ../src/newamp1.c ../src/mbest.c ../src/kiss_fft.c ../src/sine.c ../src/nlp.c ../src/dump.c ../src/octave.c ${CODEBOOKS})
target_link_libraries(tnewamp1 codec2)
add_executable(compare_ints compare_ints.c)
add_executable(compare_floats compare_floats.c)
add_executable(tvq_mbest tvq_mbest.c)
add_executable(tfreedv_800XA_rawdata tfreedv_800XA_rawdata.c)
target_link_libraries(tfreedv_800XA_rawdata codec2)
add_executable(tfreedv_2400A_rawdata tfreedv_2400A_rawdata.c)
target_link_libraries(tfreedv_2400A_rawdata codec2)
add_executable(tfreedv_2400B_rawdata tfreedv_2400B_rawdata.c)
target_link_libraries(tfreedv_2400B_rawdata codec2)
add_executable(tfsk_llr tfsk_llr.c)
target_link_libraries(tfsk_llr codec2 m)
add_executable(thash thash.c)
target_link_libraries(thash codec2 m)
add_executable(tqam16 tqam16.c)
target_link_libraries(tqam16 codec2 m)
add_executable(t16_8 t16_8.c ../src/fdmdv.c ../src/kiss_fft.c)
target_link_libraries(t16_8 codec2)
add_executable(t16_8_short t16_8_short.c ../src/fdmdv.c ../src/kiss_fft.c)
target_link_libraries(t16_8_short codec2)
add_executable(t48_8 t48_8.c ../src/fdmdv.c ../src/kiss_fft.c)
target_link_libraries(t48_8 codec2)
add_executable(t48_8_short t48_8_short.c ../src/fdmdv.c ../src/kiss_fft.c)
target_link_libraries(t48_8_short codec2)
add_executable(tquisk_filter tquisk_filter.c)
target_link_libraries(tquisk_filter codec2)
# Build CML as part of unit test setup
find_program(OCTAVE_CMD octave-cli REQUIRED)
message("Octave command: ${OCTAVE_CMD}")
include(ExternalProject)
set(CML_PATH ${CMAKE_CURRENT_BINARY_DIR}/../cml)
ExternalProject_Add(cml
GIT_REPOSITORY https://github.com/drowe67/cml.git
SOURCE_DIR ${CML_PATH}
BUILD_IN_SOURCE 1
CONFIGURE_COMMAND true # No configuration required
BUILD_COMMAND cd ${CMAKE_CURRENT_BINARY_DIR}/../cml && make
INSTALL_COMMAND true # No installation required
)
# Create fading files (used for channel simulation) as part of unit test setup
add_custom_target(fading_files ALL
DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/fast_fading_samples.float
DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/faster_fading_samples.float
)
add_custom_command(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/fast_fading_samples.float
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/faster_fading_samples.float
COMMAND cd ${CMAKE_CURRENT_SOURCE_DIR} && ./fading_files.sh ${CMAKE_CURRENT_BINARY_DIR}
)
add_executable(freedv_700d_comptx freedv_700d_comptx.c)
add_executable(freedv_700d_comprx freedv_700d_comprx.c)
if(LPCNET AND lpcnetfreedv_FOUND)
target_link_libraries(freedv_700d_comptx m codec2 lpcnetfreedv)
target_link_libraries(freedv_700d_comprx m codec2 lpcnetfreedv)
else()
target_link_libraries(freedv_700d_comptx m codec2)
target_link_libraries(freedv_700d_comprx m codec2)
endif()
add_executable(golay23 ../src/golay23.c)
target_compile_options(golay23 PUBLIC -DGOLAY23_UNITTEST)
add_executable(golay23_runtime_tables ../src/golay23.c)
target_compile_options(golay23_runtime_tables PUBLIC -DGOLAY23_UNITTEST -DRUN_TIME_TABLES)
add_executable(mksine mksine.c)
target_link_libraries(mksine m)
add_executable(vq_mbest vq_mbest.c)
target_link_libraries(vq_mbest codec2)
add_executable(tesno_est tesno_est.c)
target_link_libraries(tesno_est m codec2)
+29
View File
@@ -0,0 +1,29 @@
#!/usr/bin/env bash
#
# Check Octave and C raw data mode waveforms have about the same
# compression - sanity check for C port of raw data modes.
#
# For manual run outside of ctest:
# cd codec/build_linux
# ../unittest/check_comp.sh ${CODEC2} ${PATH}:${CODEC2}/build_linux/src
CODEC2=$1
PATH=$2:$PATH
set -x
octave_log=$(mktemp)
ch_log=$(mktemp)
echo "warning ('off', 'Octave:data-file-in-path');
ofdm_ldpc_tx('test_datac0.raw','datac0',1,100,'awgn','bursts',10,'txclip');
quit" | DISPLAY="" octave-cli -p ${CODEC2}/octave 1>${octave_log}
oct_rms=$(cat ${octave_log} | tr -s ' ' | grep 'RMS:' | cut -d' ' -f4)
oct_cpapr=$(cat ${octave_log} | grep 'RMS:' | tr -s ' ' | cut -d' ' -f6)
freedv_data_raw_tx datac0 /dev/zero - --delay 1000 --testframes 10 --bursts 10 --clip 1 --txbpf 1 | \
ch - /dev/null 2>${ch_log}
ch_rms=$(cat ${ch_log} | grep RMS | tr -s ' ' | cut -d' ' -f5)
ch_cpapr=$(cat ${ch_log} | grep RMS | tr -s ' ' | cut -d' ' -f7)
# Allow 5% difference
python3 -c "import sys; sys.exit(0) if abs((${oct_rms} - ${ch_rms})/${oct_rms}) < 0.05 else sys.exit(1)"
python3 -c "import sys; sys.exit(0) if abs((${oct_cpapr} - ${ch_cpapr})/${oct_cpapr}) < 0.05 else sys.exit(1)"
+58
View File
@@ -0,0 +1,58 @@
#!/usr/bin/env bash
#
# Check peak level of each FreeDV waveform is about the same to present
# consistent drive to transmitters.
#
# For manual run outside of ctest:
# cd codec/build_linux
# PATH=${PATH}:${HOME}/codec2/build_linux/src
# ./unittest/check_peak.sh
# OR:
#
voice_test() {
mode=$1
echo -n "$mode "
f=$(mktemp)
freedv_tx $mode ../raw/ve9qrp_10s.raw $f --clip 1
octave_cmd="cd ../octave;
t=load_raw('${f}');
mx=max(t); printf('%d ',max(t));
if (mx > 16000) && (mx < 17000) printf('PASS\n') else printf('FAIL\n') end"
octave-cli -qf --eval "$octave_cmd"
}
data_test() {
mode=$1
echo -n "$mode "
f=$(mktemp)
freedv_data_raw_tx --framesperburst 2 --bursts 3 --testframes 6 $mode /dev/zero $f 2>/dev/null
octave_cmd="cd ../octave;
t=load_raw('${f}');
mx=max(t); printf('%d ',max(t));
if (mx > 16000) && (mx < 17000) printf('PASS\n') else printf('FAIL\n') end"
octave-cli -qf --eval "$octave_cmd"
}
if [ "$1" == "LPCNet" ]; then
# these don't get run unless we build with LPCNet
voice_test "2020"
voice_test "2020B"
else
voice_test "1600"
voice_test "700C"
voice_test "700D"
voice_test "700E"
voice_test "800XA"
voice_test "2400A"
voice_test "2400B"
data_test "datac0"
data_test "datac1"
data_test "datac3"
data_test "datac4"
data_test "datac13"
data_test "datac14"
fi
exit 0
+15
View File
@@ -0,0 +1,15 @@
#!/usr/bin/env bash
# check_real_comp.sh
# Check the output of freedv_tx() and the real part of freedv_comptx() match,
# as they use different code paths. Run from codec2/unittest, set path to
# include codec2/build/misc and codec2/build/unittest
set -x
cat ../raw/ve9qrp_10s.raw | freedv_700d_tx > tx_700d.int16
cat ../raw/ve9qrp_10s.raw | freedv_700d_comptx > tx_700d.iq16
echo "tx_real=load_raw('tx_700d.int16'); tx_comp=load_raw('tx_700d.iq16'); \
tx_comp=tx_comp(1:2:end)+j*tx_comp(2:2:end); \
diff = sum(real(tx_comp)-tx_real); printf('diff: %f\n', diff); \
if diff < 1, quit(0), end; \
quit(1)" | octave-cli -p ../octave -qf
+87
View File
@@ -0,0 +1,87 @@
/* compare floats - a test utility */
#include <errno.h>
#include <getopt.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
/* Declarations */
/* Globals */
/* Main */
int main(int argc, char *argv[]) {
char usage[] = "Usage: %s [-t tolerance] file1 file2\n";
float tol = .001;
int opt;
while ((opt = getopt(argc, argv, "t:")) != -1) {
switch (opt) {
case 't':
tol = atof(optarg);
break;
default:
fprintf(stderr, usage, argv[0]);
exit(1);
}
}
if ((optind + 2) > argc) {
fprintf(stderr, usage, argv[0]);
exit(1);
}
char *fname1 = argv[optind++];
char *fname2 = argv[optind++];
FILE *f1 = fopen(fname1, "rb");
if (f1 == NULL) {
fprintf(stderr, "Error opening file1 \"%s\": ", fname1);
perror(NULL);
exit(1);
}
FILE *f2 = fopen(fname2, "rb");
if (f2 == NULL) {
fprintf(stderr, "Error opening file2 \"%s\": ", fname2);
perror(NULL);
exit(1);
}
float data1, data2;
int count = 0;
int errors = 0;
double rms_sum = 0;
while (fread(&data1, sizeof(float), 1, f1)) {
if (!fread(&data2, sizeof(float), 1, f2)) {
fprintf(stderr, "Error: file2 is shorter!");
exit(1);
}
float err = fabsf((data1 - data2) / data1);
if (err > tol) {
errors++;
printf("%d %g %g %g\n", count, data1, data2, err);
}
rms_sum += (err * err);
count++;
}
if (fread(&data2, sizeof(float), 1, f2)) {
fprintf(stderr, "Error: file1 is shorter\n");
exit(1);
}
if (errors) {
printf("Fail: %d errors\n", errors);
printf(" rms error = %g\n", ((double)rms_sum / count));
exit(1);
} else
printf("Pass\n");
exit(0);
} // main
/* vi:set ts=4 et sts=4: */
+160
View File
@@ -0,0 +1,160 @@
/* compare ints - a test utility */
#include <errno.h>
#include <getopt.h>
#include <inttypes.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
/* Declarations */
/* Globals */
/* Functions */
int get_data(FILE *f, int64_t *dd, int signed_flag, int bytes) {
int res;
int8_t d_8;
int16_t d_16;
uint8_t d_u8;
uint16_t d_u16;
// TODO Loop on reads until, but catch EOF!!
if (signed_flag) {
switch (bytes) {
case 1:
res = fread(&d_8, bytes, 1, f);
*dd = d_8;
break;
case 2:
res = fread(&d_16, bytes, 1, f);
*dd = d_16;
break;
default:
fprintf(stderr, "Error: unsupported size %d bytes\n", bytes);
exit(1);
}
} else { // unsigned
switch (bytes) {
case 1:
res = fread(&d_u8, bytes, 1, f);
*dd = d_u8;
break;
case 2:
res = fread(&d_u16, bytes, 1, f);
*dd = d_u16;
break;
default:
fprintf(stderr, "Error: unsupported size %d bytes\n", bytes);
exit(1);
}
}
if (res != 1)
return (0);
else
return (1);
}
/* Main */
int main(int argc, char *argv[]) {
char usage[] =
"Usage: %s [-b size_in_bytes] [-c] [-s] [-t tolerance] [-n "
"numerrorstoexit] file1 file2\n";
int bytes = 1;
int count_errors = 0;
int signed_flag = 0;
int tol = 1;
int numerrorstoexit = -1;
int opt;
while ((opt = getopt(argc, argv, "b:cst:n:")) != -1) {
switch (opt) {
case 'b':
bytes = atoi(optarg);
break;
case 'c':
count_errors = 1;
break;
case 's':
signed_flag = 1;
break;
case 'n':
numerrorstoexit = atoi(optarg);
break;
case 't':
tol = atof(optarg);
break;
default:
fprintf(stderr, usage, argv[0]);
exit(1);
}
}
if ((optind + 2) > argc) {
fprintf(stderr, usage, argv[0]);
exit(1);
}
char *fname1 = argv[optind++];
char *fname2 = argv[optind++];
FILE *f1 = fopen(fname1, "rb");
if (f1 == NULL) {
fprintf(stderr, "Error opening file1 \"%s\": ", fname1);
perror(NULL);
exit(1);
}
FILE *f2 = fopen(fname2, "rb");
if (f2 == NULL) {
fprintf(stderr, "Error opening file2 \"%s\": ", fname2);
perror(NULL);
exit(1);
}
// Convert inputs to SIGNED long values
int64_t data1, data2;
int count = 0;
int errors = 0;
int rms_sum = 0;
while (get_data(f1, &data1, signed_flag, bytes)) {
if (!get_data(f2, &data2, signed_flag, bytes)) {
fprintf(stderr, "Error: file2 is shorter\n");
exit(1);
}
uint64_t err = llabs(data1 - data2);
if (err > tol) {
errors++;
printf("%d %" PRId64 " %" PRId64 "\n", count, data1, data2);
if (numerrorstoexit != -1)
if (errors > numerrorstoexit) {
printf("reached errors: %d, bailing!", numerrorstoexit);
exit(1);
}
}
rms_sum += (err * err);
count++;
}
if (get_data(f2, &data2, signed_flag, bytes)) {
fprintf(stderr, "Error: file1 is shorter\n");
exit(1);
}
if (count_errors)
exit(errors);
else {
if (errors) {
printf("Fail: %d errors\n", errors);
printf(" rms error = %f\n", ((double)rms_sum / count));
exit(1);
} else
printf("Pass\n");
exit(0);
}
} // main
/* vi:set ts=4 et sts=4: */
+14
View File
@@ -0,0 +1,14 @@
#!/usr/bin/env bash
#
# Generate fading files used for channel simulation
output_path=$1
echo "Generating fading files ......"
cmd='cd ../octave; pkg load signal; ch_fading("'${output_path}'/fast_fading_samples.float", 8000, 1.0, 8000*60)'
octave --no-gui -qf --eval "$cmd"
[ ! $? -eq 0 ] && { echo "octave failed to run correctly .... exiting"; exit 1; }
cmd='cd ../octave; pkg load signal; ch_fading("'${output_path}'/faster_fading_samples.float", 8000, 2.0, 8000*60)'
octave --no-gui -qf --eval "$cmd"
[ ! $? -eq 0 ] && { echo "octave failed to run correctly .... exiting"; exit 1; }
exit 0
+140
View File
@@ -0,0 +1,140 @@
/*---------------------------------------------------------------------------*\
FILE........: freedv_700d_comprx.c
AUTHOR......: David Rowe
DATE CREATED: July 2022
Complex valued rx to support ctests. Includes a few operations that will
only work if complex Tx and Rx signals are being handled correctly.
\*---------------------------------------------------------------------------*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "codec2_cohpsk.h"
#include "comp_prim.h"
#include "freedv_api.h"
#include "freedv_api_internal.h"
#include "ofdm_internal.h"
int main(int argc, char *argv[]) {
/* with no arguments then run with no test code */
int test_num = 0;
if (argc == 2) {
if (strcmp(argv[1], "tx") == 0) {
test_num = 1;
}
if (strcmp(argv[1], "rx") == 0) {
test_num = 2;
}
}
fprintf(stderr, "%d\n", test_num);
struct freedv *freedv;
freedv = freedv_open(FREEDV_MODE_700D);
assert(freedv != NULL);
/* note API functions to tell us how big our buffers need to be */
short speech_out[freedv_get_n_max_speech_samples(freedv)];
short demod_in[2 * freedv_get_n_max_modem_samples(freedv)];
COMP demod_in_comp[2 * freedv_get_n_max_modem_samples(freedv)];
/* set up small freq offset */
float foff_hz = 25;
COMP phase_ch;
phase_ch.real = 1.0;
phase_ch.imag = 0.0;
/* set complex sine wave interferer at -fc */
COMP interferer_phase = {1.0, 0.0};
COMP interferer_freq;
interferer_freq.real =
cos(2.0 * M_PI * freedv->ofdm->tx_centre / FREEDV_FS_8000);
interferer_freq.imag =
sin(2.0 * M_PI * freedv->ofdm->tx_centre / FREEDV_FS_8000);
interferer_freq = cconj(interferer_freq);
/* log a file of demod input samples for plotting in Octave */
FILE *fdemod = fopen("demod.f32", "wb");
assert(fdemod != NULL);
/* measure demod input power, interferer input power */
float power_d = 0.0;
float power_interferer = 0.0;
int frames = 0, sum_sync = 0, frames_snr = 0;
float sum_snr = 0.0;
size_t nin, nout;
nin = freedv_nin(freedv);
while (fread(demod_in, sizeof(short), 2 * nin, stdin) == 2 * nin) {
for (int i = 0; i < nin; i++) {
demod_in_comp[i].real = (float)demod_in[2 * i];
demod_in_comp[i].imag = (float)demod_in[2 * i + 1];
// demod_in_comp[i].imag = 0;
}
if (test_num == 1) {
/* So Tx is a complex OFDM signal centered at +fc. A small
shift fd followed by Re{} will only work if Tx is complex.
If Tx is real, neg freq components at -fc+fd will be
aliased on top of fc+fd wanted signal by Re{} operation.
This can be tested by setting demod_in_comp[i].imag = 0
above */
fdmdv_freq_shift_coh(demod_in_comp, demod_in_comp, foff_hz,
FREEDV_FS_8000, &phase_ch, nin);
for (int i = 0; i < nin; i++) demod_in_comp[i].imag = 0.0;
}
if (test_num == 2) {
/* a complex sinewave (carrier) at -fc will only be ignored if
Rx is treating signal as complex, otherwise if real a +fc
alias will appear in the middle of our wanted signal at
+fc, this can be tested by setting demod_in_comp[i].imag =
0 below */
for (int i = 0; i < nin; i++) {
COMP a = fcmult(2E4, interferer_phase);
interferer_phase = cmult(interferer_phase, interferer_freq);
power_interferer += a.real * a.real + a.imag * a.imag;
COMP d = demod_in_comp[i];
power_d += d.real * d.real + d.imag * d.imag;
demod_in_comp[i] = cadd(d, a);
// demod_in_comp[i].imag = 0;
}
}
/* useful to take a look at this with Octave */
fwrite(demod_in_comp, sizeof(COMP), nin, fdemod);
nout = freedv_comprx(freedv, speech_out, demod_in_comp);
nin = freedv_nin(freedv); /* call me on every loop! */
fwrite(speech_out, sizeof(short), nout, stdout);
int sync;
float snr_est;
freedv_get_modem_stats(freedv, &sync, &snr_est);
fprintf(stderr, "sync: %d snr_est: %f\n", sync, snr_est);
frames++;
sum_sync += sync;
if (sync) {
sum_snr += snr_est;
frames_snr++;
}
}
fclose(fdemod);
freedv_close(freedv);
if (test_num == 2)
fprintf(stderr, "Demod/Interferer power ratio: %3.2f dB\n",
10 * log10(power_d / power_interferer));
float snr_av = sum_snr / frames_snr;
fprintf(stderr, "frames: %d sum_sync: %d snr_av: %3.2f dB\n", frames,
sum_sync, snr_av);
if (snr_av > 8.0)
return 0;
else
return 1;
}
+44
View File
@@ -0,0 +1,44 @@
/*---------------------------------------------------------------------------*\
freedv_comptx.c
Complex valued Tx to support ctests.
\*---------------------------------------------------------------------------*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "freedv_api.h"
int main(int argc, char *argv[]) {
struct freedv *freedv;
freedv = freedv_open(FREEDV_MODE_700D);
assert(freedv != NULL);
/* handy functions to set buffer sizes */
int n_speech_samples = freedv_get_n_speech_samples(freedv);
short speech_in[n_speech_samples];
int n_nom_modem_samples = freedv_get_n_nom_modem_samples(freedv);
COMP mod_out[n_nom_modem_samples];
short mod_out_short[2 * n_nom_modem_samples];
/* OK main loop --------------------------------------- */
while (fread(speech_in, sizeof(short), n_speech_samples, stdin) ==
n_speech_samples) {
freedv_comptx(freedv, mod_out, speech_in);
for (int i = 0; i < n_nom_modem_samples; i++) {
mod_out_short[2 * i] = mod_out[i].real;
mod_out_short[2 * i + 1] = mod_out[i].imag;
}
fwrite(mod_out_short, sizeof(short), 2 * n_nom_modem_samples, stdout);
}
freedv_close(freedv);
return 0;
}
+890
View File
@@ -0,0 +1,890 @@
short hts1a_raw[] = {
-14, -14, -8, -7, -11, -15, -14, -16, -24,
-26, -25, -26, -22, -22, -24, -19, -19, -19,
-26, -28, -28, -21, -16, -14, -19, -19, -18,
-18, -16, -18, -26, -28, -35, -28, -19, -12,
-12, -14, -15, -21, -16, -12, -9, -11, -5,
-8, -7, -5, -8, -8, -7, 3, 3, -1,
-2, -5, -1, -7, -5, -4, -4, -7, -5,
-9, -8, -12, -21, -21, -28, -28, -24, -25,
-29, -29, -31, -32, -28, -31, -35, -26, -35,
-31, -28, -32, -26, -21, -22, -16, -15, -14,
-18, -12, -19, -12, -12, -16, -15, -16, -16,
-16, -12, -15, -12, -18, -18, -15, -19, -18,
-16, -14, -15, -16, -16, -16, -14, -16, -11,
-4, -8, -8, -9, -8, -15, -12, -11, -12,
-9, -11, -8, -11, -14, -11, -18, -16, -14,
-14, -9, -5, -8, -15, -16, -14, -16, -18,
-15, -31, -32, -19, -15, -12, -16, -15, -18,
-14, -12, -12, -16, -24, -25, -19, -18, -22,
-21, -19, -16, -16, -14, -16, -24, -24, -19,
-24, -24, -19, -21, -24, -25, -28, -25, -25,
-26, -18, -12, -22, -25, -26, -25, -24, -24,
-22, -21, -19, -16, -15, -12, -12, -14, -8,
-12, -11, -5, -1, 0, 0, -1, -4, -4,
-5, -7, -7, -11, -8, -11, -5, -4, -2,
-8, -5, -12, -14, -14, -14, -12, -12, -7,
-16, -21, -22, -22, -25, -28, -24, -31, -32,
-33, -33, -35, -31, -29, -32, -36, -36, -35,
-35, -32, -26, -35, -29, -22, -18, -11, -16,
-14, -14, -11, -8, -8, -12, -11, -24, -25,
-12, -8, -7, -9, -5, -9, -8, -4, -7,
-5, -11, -11, -4, -5, -5, -8, -4, -8,
-4, -8, -14, -8, -9, -12, -11, -15, -22,
-21, -19, -22, -15, -22, -18, -15, -16, -18,
-12, -14, -21, -19, -16, -19, -21, -22, -21,
-25, -19, -26, -21, -19, -21, -19, -16, -15,
-18, -9, -8, -9, -8, -14, -19, -15, -16,
-16, -16, -12, -11, -12, -11, -11, -7, -19,
-18, -14, -28, -28, -26, -28, -31, -31, -24,
-25, -26, -29, -28, -31, -28, -24, -21, -19,
-21, -16, -24, -24, -18, -16, -19, -12, -9,
-12, -14, -12, -7, -9, -9, -5, -5, -9,
-7, -5, -5, -4, -21, -25, -12, -5, 5,
10, 10, -1, -2, 3, -4, -7, -8, -5,
-11, -12, -25, -26, -24, -33, -32, -29, -35,
-36, -33, -38, -42, -38, -38, -36, -36, -33,
-32, -38, -32, -28, -22, -18, -14, -9, -15,
-9, -8, -11, -2, -2, -5, -2, 3, -2,
-1, -4, -7, -12, -12, -12, -16, -15, -11,
-11, -8, -7, -5, -9, -12, -12, -18, -19,
-26, -26, -22, -22, -24, -21, -15, -12, -14,
-18, -16, -15, -26, -26, -28, -26, -26, -25,
-28, -25, -14, -12, -14, -18, -24, -14, -9,
-5, -7, -9, -7, -8, -14, -11, -8, -15,
-8, -7, -5, -2, -5, -5, -8, -15, -15,
-16, -33, -35, -25, -15, -14, -18, -22, -18,
-22, -24, -29, -31, -32, -33, -31, -36, -31,
-25, -31, -33, -28, -26, -22, -25, -25, -16,
-18, -16, -15, -15, -14, -11, -7, -5, -4,
-8, -4, -5, -8, -4, 0, -7, -1, -1,
-9, -11, -12, -14, -15, -9, -8, -7, -5,
-11, -12, -12, -19, -15, -16, -16, -12, -18,
-19, -19, -18, -15, -19, -24, -19, -25, -28,
-26, -26, -26, -28, -26, -22, -21, -25, -19,
-16, -16, -9, -9, -12, -7, -11, -7, -12,
-15, -18, -24, -16, -16, -12, -19, -18, -16,
-25, -24, -22, -25, -21, -25, -24, -24, -19,
-15, -26, -19, -14, -19, -16, -18, -26, -28,
-16, -4, 10, 15, 12, 13, 10, 8, 8,
-7, -11, -19, -29, -31, -33, -29, -26, -26,
-26, -22, -15, -16, -15, -19, -18, -19, -12,
-15, -21, -21, -24, -22, -19, -19, -15, -16,
-16, -24, -18, -25, -24, -24, -26, -21, -22,
-22, -18, -16, -19, -22, -22, -19, -24, -19,
-16, -19, -15, -16, -12, -4, -9, -19, -12,
-15, -19, -16, -16, -15, -14, -12, -12, -11,
-8, -9, -9, -12, -11, -11, -8, -9, -5,
5, -1, -1, -4, -8, -8, -8, -9, -8,
-7, -11, -19, -19, -18, -26, -21, -24, -24,
-26, -35, -32, -36, -31, -26, -28, -26, -25,
-22, -26, -35, -35, -36, -45, -45, -33, -28,
-19, -16, -9, -14, -12, -5, -11, -8, -7,
-5, -5, 2, -2, -5, 0, -7, -5, -11,
-14, -14, -8, -8, -7, -9, -11, -11, -15,
-15, -14, -12, -14, -18, -16, -9, -12, -12,
-14, -16, -22, -25, -26, -31, -29, -26, -29,
-25, -22, -19, -18, -24, -21, -24, -19, -12,
-15, -15, -16, -15, -14, -16, -15, -16, -24,
-19, -22, -26, -24, -19, -18, -19, -15, -11,
-5, -2, -2, -4, -7, -4, -8, -8, -11,
-15, -11, -9, -7, -9, -4, 2, -8, -1,
-2, -12, -9, -15, -21, -31, -38, -32, -32,
-35, -31, -28, -33, -32, -35, -33, -33, -36,
-36, -42, -45, -43, -33, -29, -25, -21, -14,
-12, -11, -16, -15, -12, -14, -15, -15, -9,
-2, -2, 0, -5, -4, -2, 0, 8, 9,
10, 12, 3, 6, 5, 9, 3, -7, -9,
-25, -32, -25, -11, -1, 2, -1, -9, -15,
-12, -15, -21, -29, -35, -39, -39, -31, -33,
-33, -26, -28, -29, -31, -33, -26, -24, -22,
-24, -21, -18, -15, -18, -26, -25, -22, -18,
-21, -24, -26, -35, -28, -26, -26, -24, -22,
-16, -18, -22, -15, -22, -24, -16, -14, -11,
-4, 3, 5, 3, 8, 8, 6, 0, 6,
3, -5, 0, 6, -5, -5, -8, -11, -14,
-19, -21, -24, -25, -28, -28, -22, -28, -38,
-38, -26, -26, -22, -32, -31, -26, -18, -12,
-11, -9, -16, -21, -19, -16, -16, -18, -12,
-12, -8, -7, -15, -16, -16, -19, -21, -22,
-22, -22, -25, -32, -29, -35, -32, -33, -33,
-29, -25, -18, -15, -16, -12, -12, -8, -2,
0, 6, 2, 0, 5, 2, -5, -1, -5,
-8, -8, -8, -9, -15, -12, -11, -5, -9,
-18, -19, -22, -16, -14, -11, -2, -2, -8,
-16, -21, -22, -19, -24, -24, -16, -16, -12,
-7, -8, -5, -8, -16, -24, -29, -35, -36,
-36, -38, -41, -38, -43, -41, -39, -35, -32,
-31, -32, -32, -28, -24, -24, -19, -16, -18,
-12, -2, -1, 5, 5, -9, -25, -21, -15,
-14, -15, -14, -18, -16, -12, -11, -5, -2,
-7, -11, -7, -18, -11, -5, -2, 0, -2,
0, 0, -5, -11, -14, -12, -14, -14, -16,
-18, -18, -29, -38, -38, -42, -46, -38, -35,
-32, -31, -16, -5, 0, 12, 19, 20, 17,
20, 20, 16, 16, 12, 5, -7, -15, -14,
-22, -25, -26, -26, -32, -38, -43, -45, -49,
-55, -56, -52, -56, -48, -39, -33, -28, -31,
-24, -29, -32, -26, -16, -7, -14, -11, -11,
-16, -18, -29, -28, -24, -14, -11, -15, -4,
-11, -12, -5, -2, 3, 5, 2, 10, 5,
6, 2, -1, -7, -7, -12, -14, -16, -22,
-21, -18, -19, -22, -11, -1, 3, 2, 2,
3, 6, -2, -7, -15, -18, -24, -26, -31,
-38, -41, -39, -36, -39, -33, -26, -24, -18,
-19, -21, -19, -18, -16, -21, -21, -15, -14,
-18, -24, -25, -31, -38, -43, -45, -46, -43,
-39, -33, -28, -19, -11, -8, -4, 5, 12,
12, 17, 16, 9, 9, 10, 6, 8, -4,
3, 0, -5, -11, -14, -28, -26, -15, -24,
-32, -32, -28, -32, -28, -18, -22, -22, -15,
-15, -24, -25, -26, -25, -16, -16, -18, -22,
-21, -26, -29, -25, -22, -19, -16, -9, -4,
2, 6, 10, 3, 2, 0, -7, -7, -14,
-16, -15, -22, -26, -29, -25, -25, -22, -29,
-35, -25, -19, -14, -15, -12, -15, -26, -24,
-29, -28, -26, -26, -29, -32, -38, -42, -38,
-33, -29, -25, -25, -21, -14, -5, 5, 6,
8, 6, 2, 0, -1, -9, -16, -18, -19,
-22, -22, -21, -15, -18, -22, -11, -8, -11,
-7, 5, 2, -2, -4, -5, -7, -5, 6,
-7, -9, -8, -19, -22, -24, -26, -29, -33,
-29, -25, -24, -21, -21, -24, -29, -28, -26,
-25, -21, -26, -26, -25, -31, -31, -38, -39,
-38, -33, -21, -9, -5, -4, 5, 3, -2,
-7, -8, -9, -15, -16, -16, -22, -21, -24,
-24, -16, -22, -21, -24, -26, -22, -16, -9,
-11, -2, 6, 2, -4, -9, -16, -21, -21,
-21, -24, -22, -22, -19, -18, -16, -12, -14,
-12, -4, -2, -4, -8, -5, -5, -11, 0,
2, -7, -8, -12, -14, -19, -24, -25, -28,
-9, -5, -29, -33, -22, -22, -21, -15, -18,
-24, -21, -26, -29, -25, -33, -29, -29, -29,
-31, -28, -28, -25, -24, -21, -25, -14, -7,
-5, -16, -19, -21, -28, -33, -38, -36, -26,
-25, -22, -8, -5, 0, 5, 10, 16, 13,
10, 8, 5, -4, -1, -7, -11, -18, -28,
-31, -42, -43, -38, -38, -22, -11, -8, -7,
6, 6, 3, 13, 8, -7, -2, -9, -16,
-11, -15, -18, -24, -28, -24, -25, -22, -25,
-28, -25, -38, -39, -35, -36, -14, -25, -42,
-16, -9, -29, -28, -26, -31, -29, -39, -29,
-12, -7, -1, -2, -4, 0, 2, -5, -15,
-21, -35, -32, -22, -19, -18, -19, -5, 6,
0, -14, -26, -11, 0, -19, -24, -24, -25,
-25, -31, -35, -32, -19, -5, -7, -1, 0,
6, 8, 0, -8, -5, 9, 9, 6, 2,
9, 13, -7, -26, -36, -35, -42, -56, -49,
-42, -42, -36, -28, -12, -12, -21, -18, -24,
-19, -22, -25, -24, -21, -18, -15, -15, -8,
-7, 3, -4, -11, -22, -22, -16, -24, -21,
-7, -22, -31, -16, -21, -11, -12, -21, -26,
-28, -19, -28, -31, -25, -38, -38, -29, -33,
-38, -33, -9, 10, 19, 5, -4, -4, -1,
-12, -21, -18, -16, -16, -19, -8, -5, -7,
-2, 0, 12, 16, 15, 2, -14, -12, -22,
-29, -42, -36, -25, -16, -18, -19, -12, 6,
2, 2, 9, -4, -11, -19, -25, -24, -28,
-35, -43, -35, -25, -41, -45, -42, -39, -35,
-41, -33, -29, -18, -2, -19, -32, -12, 3,
-8, -11, -26, -35, -29, -29, -24, -15, -9,
-8, -1, 2, 0, -1, -2, 2, -1, -8,
-9, -18, -21, -26, -32, -35, -26, -15, -9,
-7, -2, -2, 2, 13, 12, -2, -11, -15,
-18, -28, -29, -24, -28, -32, -31, -28, -15,
-9, -21, -4, 10, -7, -5, -14, -18, -9,
-21, -25, -29, -32, -36, -31, -16, -24, -18,
-18, -22, -25, -32, -31, -12, -8, -8, 2,
-2, 15, 9, -14, -9, -9, 0, -1, -22,
-24, -18, -32, -29, -29, -39, -33, -24, -41,
-33, -18, -33, -28, -32, -28, -24, -35, -32,
-24, -21, -14, -21, -25, -18, -9, -8, -11,
0, 3, 5, -5, -14, -19, -15, -9, 2,
-2, 0, 12, 10, 19, 22, 8, -2, -5,
-9, -22, -16, -5, -16, -21, -19, -16, -16,
-26, -29, -28, -25, -1, -2, -15, -19, -29,
-35, -35, -38, -36, -38, -29, -46, -56, -38,
-45, -62, -55, -48, -33, -42, -79, -93, 33,
207, 203, 112, 30, -29, -28, -70, -148, -176,
-96, -69, -87, -32, -45, -24, 46, 64, 76,
46, 43, 60, 39, 20, 12, 6, -35, -26,
-28, -67, -8, -18, -55, 6, -38, -75, -24,
-5, -5, -25, -25, -35, 22, 93, 23, -28,
-55, -83, -42, -83, -103, -56, -22, 44, 3,
-33, 6, 8, 12, 17, -5, -131, -189, 100,
265, 40, 49, 135, -59, 12, 51, -123, -87,
-182, -298, -226, -192, -205, -198, -49, 30, 39,
183, 238, 183, 141, 187, 132, 83, 176, 16,
-116, -90, -118, -138, -189, -123, -137, -147, -9,
-28, 49, 125, 66, 29, 43, 46, -70, -100,
-75, -121, -117, -109, -58, -28, 29, 77, 74,
128, 118, 63, 12, 56, 158, 90, 0, -29,
-127, -103, -33, -137, -140, -18, -43, -66, -7,
-45, -73, 6, 42, -35, -111, -104, -92, -147,
-526, -1264, -1575, -765, 438, 997, 1207, 1339, 1320,
1748, 2310, 2055, 1176, 299, -434, -1009, -1257, -1632,
-2439, -2960, -2928, -2615, -2054, -1308, -736, -303, 445,
1377, 2112, 2587, 2728, 2464, 2061, 1838, 1525, 888,
182, -511, -1203, -1516, -1425, -1447, -1556, -1386, -1076,
-642, -16, 437, 547, 662, 915, 1092, 1055, 883,
519, 53, -259, -472, -771, -1037, -1070, -981, -814,
-511, -203, 2, 251, 519, 616, 631, 648, 582,
476, 391, 271, -18, -305, -348, -399, -478, -389,
-348, -351, -222, -121, -386, -1735, -3442, -2533, 573,
2049, 2601, 3434, 2708, 2736, 4175, 4209, 2551, 296,
-1119, -2477, -3494, -3361, -4179, -5485, -5636, -4580, -3316,
-2031, -22, 1181, 1763, 3271, 4874, 5510, 5224, 4612,
3186, 1425, 652, -147, -1598, -2834, -3484, -3818, -3613,
-2537, -1594, -1278, -574, 734, 1997, 2818, 3281, 3233,
2488, 1960, 1561, 520, -672, -1616, -2288, -2651, -2450,
-2115, -2137, -1772, -998, -341, 312, 891, 1208, 1322,
1574, 1687, 1285, 810, 427, -18, -404, -664, -1030,
-1306, -1064, -656, -397, -182, -16, 170, 514, 922,
1023, 956, 915, 626, -42, -1701, -4508, -5074, -1033,
2559, 2488, 3256, 3525, 2940, 5224, 6208, 3992, -50,
-2480, -3034, -4628, -4714, -5427, -7560, -7667, -5645, -2855,
-1414, 621, 2572, 3066, 5038, 7469, 7951, 6155, 4438,
2841, 544, -373, -1233, -3376, -5177, -5271, -4414, -3653,
-2121, -788, -392, 720, 2701, 3972, 4022, 3914, 3451,
2206, 1426, 796, -632, -2163, -3381, -4159, -4113, -3045,
-1776, -1312, -743, 275, 1098, 1851, 2335, 2131, 1527,
1259, 1200, 713, 73, -516, -1176, -1677, -1755, -1672,
-1670, -1336, -625, 37, 718, 1435, 1620, 1412, 1445,
1137, 438, 388, 645, 505, 194, -907, -3965, -7195,
-5524, 1111, 4509, 3145, 3535, 3294, 3846, 7424, 7315,
3111, -1601, -3095, -3366, -4972, -5281, -7332, -9647, -8219,
-4862, -1956, -808, 1125, 2677, 3750, 7120, 9202, 7949,
5515, 3979, 2538, 716, -4, -1881, -5029, -6072, -5264,
-4448, -3573, -2254, -1472, -802, 1513, 3981, 4536, 4289,
4083, 3336, 2484, 2141, 963, -1417, -3088, -3937, -4523,
-3991, -2694, -2023, -1694, -593, 580, 1431, 2310, 2417,
1777, 1433, 1598, 1241, 374, -69, -753, -1490, -1449,
-1357, -1625, -1628, -1176, -845, -356, 723, 1418, 1370,
1530, 1663, 1222, 1057, 1140, 345, -576, -628, -1435,
-4343, -6010, -2139, 2712, 2651, 2832, 3744, 2781, 5433,
7295, 4184, 117, -2377, -2990, -4553, -4853, -5022, -7143,
-6439, -4074, -2314, -1340, -260, 822, 720, 2570, 4741,
4967, 4932, 4619, 4026, 2856, 2222, 1221, -911, -1898,
-2510, -2962, -2593, -1870, -1447, -1553, -1085, -648, -342,
240, 418, 505, 623, 1054, 1416, 1337, 1193, 544,
-101, -348, -407, -426, -604, -573, -519, -297, 9,
-235, -589, -927, -952, -563, -310, 9, 168, 219,
369, 245, 40, -235, -451, -287, 46, 415, 683,
836, 842, 670, 489, 316, 166, 54, -36, -65,
36, -96, -883, -1693, -1570, -756, -117, 340, 975,
1918, 2596, 2086, 1068, 199, -509, -625, -662, -529,
71, -33, -597, -1067, -1512, -1670, -1870, -1901, -1547,
-974, -269, 142, 260, 272, 56, -164, -24, 509,
1019, 1292, 1496, 1616, 1660, 1613, 1357, 966, 660,
605, 597, 480, 275, -134, -671, -1134, -1376, -1473,
-1534, -1564, -1496, -1237, -812, -402, -128, 64, 267,
431, 522, 479, 332, 320, 468, 565, 548, 407,
180, -63, -264, -455, -638, -685, -574, -358, -32,
238, 328, 329, 371, 427, 438, 448, 424, 383,
441, 482, 29, -1189, -2337, -1898, -396, 261, 599,
1708, 2529, 2347, 1500, 312, -430, -148, 278, 227,
432, 471, -532, -1710, -2265, -2405, -2323, -2014, -1597,
-1182, -690, -406, -567, -617, -358, -89, 343, 1088,
1780, 1952, 1833, 1659, 1418, 1391, 1493, 1446, 1302,
1180, 911, 415, -45, -492, -962, -1186, -1142, -1095,
-1100, -1196, -1422, -1556, -1353, -995, -586, -59, 257,
287, 287, 197, 3, -36, 156, 418, 737, 932,
759, 427, 165, -84, -252, -225, -148, -84, 34,
39, -96, -169, -168, -72, 63, 179, 360, 539,
485, -67, -1439, -2707, -1793, 379, 1014, 1030, 2144,
2711, 2106, 1211, -4, -525, 364, 955, 628, 584,
170, -1432, -2769, -2892, -2525, -2075, -1551, -1390, -1363,
-1025, -944, -1087, -615, 145, 594, 1099, 1758, 1969,
1712, 1453, 1289, 1384, 1843, 2083, 1800, 1387, 975,
400, -86, -260, -446, -681, -781, -918, -1183, -1481,
-1768, -1877, -1558, -1020, -515, -28, 192, 9, -261,
-322, -178, 98, 444, 706, 805, 761, 492, 145,
-5, 6, 20, 53, 81, 27, -103, -249, -325,
-219, -9, 139, 291, 473, 547, 400, 42, -904,
-2367, -2330, -46, 1336, 924, 1518, 2371, 1895, 1372,
564, -271, 454, 1397, 1020, 578, 328, -986, -2542,
-2821, -2398, -1994, -1564, -1564, -1956, -1922, -1598, -1434,
-879, 175, 771, 983, 1429, 1726, 1615, 1534, 1616,
1721, 2049, 2369, 2073, 1426, 888, 304, -192, -219,
-94, -284, -630, -961, -1359, -1618, -1696, -1748, -1454,
-778, -324, -128, -18, -195, -423, -382, -226, 70,
483, 701, 628, 428, 210, 10, -22, 114, 227,
272, 257, 169, 25, -75, -138, -118, 102, 328,
287, 159, 135, -135, -1224, -2588, -2139, 243, 1339,
861, 1548, 2238, 1746, 1333, 575, 20, 917, 1538,
977, 662, 405, -964, -2463, -2640, -2228, -1890, -1544,
-1616, -1970, -1980, -1750, -1567, -867, 258, 771, 834,
1153, 1445, 1448, 1465, 1600, 1732, 2025, 2283, 1981,
1387, 883, 305, -131, -60, 112, -138, -472, -713,
-1204, -1632, -1686, -1616, -1274, -676, -513, -574, -314,
-246, -516, -545, -349, -28, 512, 853, 680, 427,
319, 192, 159, 289, 343, 264, 254, 193, -29,
-143, -148, -107, 152, 371, 271, -62, -913, -2303,
-2531, -572, 992, 788, 1179, 2054, 1888, 1608, 1079,
285, 625, 1391, 1211, 911, 815, -239, -1829, -2473,
-2303, -1946, -1557, -1548, -1949, -2143, -1976, -1775, -1202,
-178, 458, 614, 910, 1272, 1374, 1395, 1466, 1518,
1745, 2114, 2090, 1680, 1224, 657, 53, -175, -52,
-45, -203, -454, -892, -1329, -1492, -1447, -1204, -814,
-655, -631, -406, -247, -334, -468, -526, -358, 83,
466, 547, 455, 343, 228, 230, 309, 335, 371,
394, 352, 231, 64, -60, -48, 104, 202, -11,
-747, -2004, -2576, -1166, 550, 643, 737, 1650, 1868,
1642, 1375, 674, 546, 1071, 1077, 943, 1115, 595,
-841, -1963, -2249, -2154, -1810, -1444, -1523, -1822, -1890,
-1813, -1512, -726, 32, 333, 650, 1095, 1306, 1343,
1429, 1450, 1429, 1627, 1802, 1692, 1440, 1016, 360,
-106, -130, -84, -157, -290, -645, -1068, -1270, -1301,
-1210, -984, -715, -532, -436, -353, -370, -471, -464,
-363, -165, 206, 496, 517, 400, 306, 192, 166,
284, 373, 383, 386, 383, 319, 224, 135, 127,
115, -491, -1799, -2388, -1161, 363, 522, 420, 1024,
1442, 1559, 1517, 1026, 750, 863, 628, 452, 800,
684, -348, -1405, -1942, -2143, -2044, -1704, -1476, -1461,
-1526, -1621, -1422, -716, 5, 325, 585, 1013, 1288,
1380, 1476, 1487, 1409, 1382, 1354, 1255, 1174, 968,
519, 129, -52, -240, -438, -441, -399, -604, -941,
-1258, -1284, -315, 394, -628, -1085, -532, -675, -126,
-277, -1573, -522, 598, 237, 407, 805, 1075, 653,
318, 362, -59, -56, -240, -233, 312, -489, -2514,
-4669, -4202, -712, 2020, 3805, 4136, 2529, 3440, 5718,
5515, 3867, 1626, -647, -2193, -2176, -1803, -3276, -5465,
-7012, -7661, -6402, -3818, -2200, -1338, 105, 1535, 3292,
5837, 7231, 6633, 5443, 4310, 3120, 2167, 863, -1674,
-4124, -5139, -5305, -4852, -3473, -1936, -1080, -52, 1583,
2812, 3833, 4309, 3278, 2191, 1920, 1330, 168, -1132,
-2559, -3548, -3413, -3045, -3000, -2440, -1105, 424, 1562,
2003, 1718, 1409, 1794, 1922, 959, -191, -649, -329,
-22, -2232, -7587, -10133, -6123, 3815, 11010, 10048, 7655,
4465, 7798, 10768, 5766, -1615, -9642, -12287, -10166, -6666,
-6229, -9606, -9109, -5832, -1326, 4475, 7162, 6468, 5722,
6975, 8255, 7576, 5333, 919, -3576, -5125, -4751, -4305,
-4140, -4257, -4346, -2541, 1238, 4948, 7386, 6832, 4496,
3434, 2708, 1527, -206, -3454, -6091, -5990, -4237, -2710,
-1963, -1253, -324, 1191, 2740, 2787, 2061, 1796, 1911,
1408, -65, -1632, -2642, -2612, -2111, -2132, -2119, -703,
1824, 2113, -3052, -10217, -9285, 2474, 15080, 14696, 9168,
3569, 1871, 8840, 7733, -863, -10981, -16373, -12835, -7172,
-2606, -4832, -7242, -3728, 1232, 7185, 9862, 7265, 2730,
1101, 4013, 4560, 2412, -866, -5642, -6610, -3630, -702,
380, 20, -726, -825, 2400, 6645, 7202, 4332, 728,
-1582, -1799, -535, -1312, -4539, -5701, -4256, -1013, 2246,
2495, 1214, 859, 1630, 1894, 788, -715, -1510, -1008,
-797, -1376, -2139, -2347, -1166, 64, 874, 1826, 2927,
1167, -7529, -13865, -7119, 8693, 17405, 11376, 5806, -706,
3910, 12384, 4898, -5841, -15567, -15268, -9199, -2915, -2,
-6828, -7298, -645, 4931, 9593, 8510, 3404, -1029, 1545,
5927, 3179, 704, -3058, -7267, -5092, -1456, 523, 578,
803, 476, 128, 4503, 7132, 4843, 2484, -314, -2661,
-1629, -312, -2667, -4845, -4645, -3556, -106, 3033, 2307,
815, 1327, 2202, 1608, 238, -1347, -2299, -1676, -1349,
-2026, -2123, -993, 185, 601, 1363, 2536, 1356, -6385,
-14216, -8118, 9341, 18113, 10028, 5714, -662, 2201, 13868,
5555, -6998, -16122, -13932, -7484, -3936, -511, -7031, -7875,
1569, 7356, 9852, 7233, 2382, -812, 2168, 7258, 2096,
-2579, -3745, -6101, -3546, -123, 70, -1707, -134, 2587,
1895, 4233, 5028, 2090, 1861, 1389, -720, -1444, -1278,
-3112, -4277, -3164, -2898, -1520, 861, 1640, 1528, 1967,
2726, 1695, 329, -492, -1697, -1969, -1929, -2346, -2041,
-600, 480, 449, 1232, 2001, -2082, -11654, -11999, 3894,
17853, 11683, 5279, 1994, -1245, 12112, 9934, -5184, -14616,
-13918, -5837, -3437, -2238, -7345, -9704, 1457, 9286, 8837,
4980, 1010, 1222, 4271, 7312, 1819, -5336, -4002, -3308,
-1866, -985, -3300, -3466, -199, 5187, 4104, 1751, 2997,
2126, 3600, 3670, -678, -3106, -1891, -1261, -2368, -2264,
-3919, -3476, 176, 1251, 1354, 400, 691, 2138, 1896,
1695, -373, -2231, -1989, -1816, -1622, -1393, -557, -261,
839, 1234, -6099, -13830, -6704, 12207, 17650, 5219, 4873,
2652, 5438, 14840, 2038, -11781, -14722, -7010, -2482, -5104,
-8182, -11393, -3796, 9197, 8387, 2298, 752, 2771, 7659,
8564, 3986, -5080, -5073, 893, -749, -2984, -6554, -6884,
-950, 4891, 4912, -866, 561, 5333, 7261, 5833, -1037,
-4326, -1758, 1554, -393, -5394, -5454, -3454, 179, 2665,
-1061, -2457, 677, 3635, 4058, 1347, -451, -846, -346,
-403, -2221, -2535, -1778, -671, 454, -2159, -10071, -12945,
-179, 17116, 12061, 403, 7550, 7610, 11287, 11279, -5843,
-14361, -7034, -1442, -6719, -11613, -12495, -8241, 4100, 9619,
-682, -1891, 5146, 9786, 11057, 5758, -974, -3229, 3852,
2712, -6304, -8209, -8290, -4368, 1412, 1119, -2455, -1243,
5528, 8521, 7421, 3537, -1505, 772, 3590, -341, -4257,
-4940, -4165, -2414, -583, -3011, -4147, 39, 2313, 2528,
2055, 1207, 1947, 2715, 1007, -2408, -2377, -414, -1230,
-1867, -3633, -9724, -12863, -3041, 14045, 10335, -2909, 9310,
13964, 12254, 12191, -3180, -10105, -345, -192, -11763, -16662,
-13357, -7749, 648, 2417, -6327, -892, 9776, 11498, 8643,
4306, 4143, 5477, 7138, 935, -8335, -6130, -4015, -5890,
-6088, -6166, -3580, 1133, 4766, 3210, 4129, 8211, 5810,
3999, 3155, 502, -186, -1683, -4503, -5169, -3667, -3607,
-4424, -3364, -1626, 496, 2096, 1967, 2167, 3503, 3642,
1586, -283, -332, -32, -1544, -4229, -7899, -12648, -8414,
7338, 10181, -4368, 4177, 17458, 13517, 13236, 4243, -5162,
1964, 2892, -11114, -17844, -12323, -7832, -6496, -4604, -7565,
-3331, 7284, 7605, 3611, 5653, 10209, 9934, 7845, 4885,
-679, -230, -329, -7610, -9785, -6780, -5146, -4804, -3114,
-808, 1753, 6140, 6461, 4152, 5323, 6386, 4387, 857,
-710, -1095, -2162, -3546, -6282, -6035, -3066, -2787, -2956,
-1881, 510, 2940, 3258, 2654, 2028, 2429, 2474, 1033,
-233, -2762, -7034, -11056, -7771, 4646, 5433, -5715, 3990,
16374, 11856, 10648, 5494, -18, 5462, 1544, -10835, -12677,
-6704, -7160, -10246, -8169, -7359, -3156, 2089, -420, 260,
6647, 8892, 6911, 7202, 8552, 5729, 4090, 1688, -3204,
-3246, -3532, -7259, -7538, -4539, -3069, -2792, -884, 932,
3421, 5464, 3526, 3536, 5345, 4133, 2168, 584, -281,
-1044, -2500, -4101, -4951, -3804, -3109, -3549, -2503, -910,
410, 1643, 2037, 1937, 2188, 2788, 2315, -572, -6503,
-9657, -1302, 7699, -1541, -6680, 9781, 13307, 6279, 7430,
3430, 3365, 5749, -4342, -10037, -4328, -4039, -10940, -9452,
-4717, -5445, -3381, -2392, -2448, 2723, 4173, 2037, 4710,
9437, 7904, 4223, 5554, 4404, 1385, -43, -3303, -4253,
-3654, -5194, -5691, -3290, -1456, -2081, -498, 1909, 2226,
3394, 3325, 2560, 4033, 3472, 1280, 1159, 354, -1428,
-2003, -2912, -3889, -3756, -3522, -3245, -2163, -801, -58,
556, 1671, 2327, 2137, 1292, -3214, -7879, -365, 8792,
-2058, -6971, 10764, 11376, 2080, 7475, 5194, 2926, 4953,
-5039, -7315, -715, -5646, -12680, -6937, -1755, -6613, -5901,
-2174, -1343, 2238, 575, -682, 5661, 9086, 4503, 2478,
7910, 6936, 1933, 1266, -19, -835, -2541, -5805, -4897,
-2034, -3317, -5455, -1761, 1654, -167, -301, 2165, 3287,
3091, 2318, 1799, 3063, 3081, -334, -481, 970, -1720,
-3691, -3239, -2820, -2991, -3206, -2470, -811, 529, 173,
716, 2623, 1564, -3845, -6581, 2529, 8041, -5370, -4499,
13582, 7503, -109, 8820, 5684, 4438, 3673, -6724, -2925,
1272, -9126, -12206, -2956, -2350, -9067, -5194, -1653, -1288,
390, -2551, 530, 6921, 4657, 798, 5477, 10107, 4739,
2049, 4361, 2961, 561, -2843, -4325, -1612, -2273, -6381,
-5260, -109, -1162, -3773, -522, 1928, 1069, 328, 1702,
3025, 2946, 2046, 955, 2895, 2338, -1672, -1136, -243,
-2378, -3644, -2902, -1987, -2112, -1495, -920, 270, 1683,
595, -992, -4240, -4096, 5163, 3365, -8693, 2644, 13871,
721, 2637, 10991, 4279, 5176, 1215, -4254, 1188, -1917,
-10538, -7706, -284, -6409, -9875, -1970, -2217, -3228, -2424,
-2575, 1869, 3423, 1156, 1937, 7458, 7611, 2627, 5239,
6736, 2968, 897, -186, -48, -1270, -3511, -4529, -2943,
-1507, -4284, -3616, -532, -870, -1332, -877, 1520, 2161,
605, 2045, 3158, 2852, 1636, 554, 1527, 479, -1720,
-2091, -1287, -1420, -2853, -1994, -642, -702, -651, -305,
809, -1228, -6310, -3284, 6444, 57, -9694, 6692, 12322,
-2264, 4655, 10553, 3995, 5217, -66, -1294, 2509, -3061,
-9247, -4914, 871, -7968, -9271, -1035, -3443, -5114, -3791,
-1539, 1162, 829, 660, 2385, 7135, 5285, 1660, 6675,
7080, 2594, 1368, 2909, 2358, -1581, -2121, -1548, -1810,
-3079, -4944, -2309, -1422, -3854, -2656, -14, 115, -634,
735, 2188, 2383, 1821, 1231, 2298, 1998, 6, -378,
210, -244, -1924, -1735, -317, -1118, -1969, -886, 5,
-368, -410, -812, -3948, -5540, 1489, 4735, -7002, -3001,
13628, 3594, -3100, 10227, 7590, 2375, 2916, 1211, 1479,
-631, -3134, -6350, -1802, -978, -10207, -6177, -1046, -4499,
-6925, -3423, 1443, -1421, -584, 2444, 2882, 5009, 3294,
3499, 5559, 5411, 3077, 1096, 4451, 2750, -1703, -631,
383, -1183, -4035, -2959, -1350, -3069, -3804, -2694, -236,
-904, -1864, 513, 1545, 1463, 455, 1019, 2475, 970,
230, 551, 798, 211, -916, -468, -216, -352, -1036,
-1298, -138, -232, -750, -580, 40, -1001, -5521, -4688,
5193, 3492, -9220, 1094, 13917, 1137, -3168, 7172, 9020,
2400, -2739, 3367, 3104, -2930, -2350, -4632, -2215, -2432,
-6993, -6593, -3793, 228, -5689, -6176, 2457, 1840, -1003,
-1275, 4305, 5439, 333, 2841, 4258, 4558, 3512, 1238,
3254, 2559, 1884, -421, -1639, 1150, -1170, -3158, -3000,
-1500, -1025, -3848, -2200, -277, -804, -1328, -773, 1228,
162, -96, 884, 812, 1234, 386, 645, 708, 381,
633, -624, -126, 146, -872, -957, -586, 223, -698,
-781, 544, 261, 200, -328, -2742, -4169, 1605, 6106,
-3718, -5298, 7649, 6832, -1261, -1390, 4067, 6767, -470,
-1764, 1864, -38, 257, -2297, -3900, -2558, -1958, -1561,
-6583, -4334, 516, -2266, -3226, -2200, 2206, 1794, -1454,
1915, 2658, 3157, 2576, 1949, 3410, 1896, 3450, 2310,
117, 1588, 1002, 917, -1493, -1639, 129, -2106, -2237,
-2421, -1808, -1254, -2307, -727, -923, -651, 449, 46,
471, 268, 1098, 710, -117, 1031, 362, -192, -90,
398, 471, -775, -349, 271, 125, -277, -403, -31,
162, 682, -104, -477, 441, 534, 682, 136, -219,
114, -447, -2159, -2436, 1852, 1779, -2721, -567, 1960,
1799, 1184, -179, 1457, 1266, 629, 1852, -82, -866,
-345, 61, -696, -3194, -1676, -726, -1898, -1857, -1936,
-516, -494, -253, 808, -301, 554, 1218, 1193, 1286,
398, 1408, 1364, 843, 1215, 506, 400, 94, 496,
632, -447, -192, -387, -386, -505, -1087, -390, -591,
-458, -196, -703, -257, -225, 87, 90, -470, 210,
258, -121, -270, -150, 267, -203, -43, 302, -101,
47, 244, 209, 63, 142, 353, 64, -158, -79,
312, -121, -985, -253, 224, -389, -603, -351, -28,
56, 376, 220, -222, 353, 635, 478, 441, 85,
-219, -768, -570, 876, 771, -478, -189, 482, 529,
329, 211, 618, 210, -94, 575, 101, -511, -508,
-454, -437, -1064, -996, -845, -1071, -617, -593, -579,
-382, -58, 471, 253, 135, 289, 926, 1221, 374,
432, 626, 645, 834, 325, 203, 84, -69, -36,
-344, -247, -287, -308, -257, -498, -158, -109, -198,
51, -60, -96, -138, -56, 98, -11, 54, -73,
-188, -216, -209, -5, -219, -244, -28, -33, 68,
-84, -84, 80, 73, 119, -16, 66, 277, 234,
102, -55, -16, 47, 2, -93, -193, -113, -55,
50, -53, -307, -24, 213, 98, 37, -31, -7,
158, 335, 243, 22, -138, -114, 100, -46, -12,
59, -209, -89, -387, -869, -338, 349, 214, -232,
57, 417, 335, 312, 414, 604, 151, 63, 478,
22, -329, -402, -315, -355, -754, -553, -526, -538,
-332, -424, -199, -28, 5, -42, 94, 415, 271,
316, 136, 90, 485, 233, 202, 189, -7, 119,
61, 238, 200, -237, -254, -49, 136, 64, -79,
20, 129, 93, 9, 163, 95, -203, -49, -24,
-295, -223, -110, -223, -396, -499, -400, -240, -33,
98, -1, 54, 179, 134, 97, 78, 77, 83,
197, 296, 85, 70, 180, -124, -99, 192, -128,
-218, 57, -193, -181, 74, -113, -140, -164, -338,
-179, 243, 366, -42, -277, -104, 326, 618, -14,
-451, -82, 427, 529, -158, -205, 119, 80, 131,
-334, -213, 29, -440, -203, -304, -836, -890, -169,
646, -67, -447, 172, 636, 871, 527, 643, 608,
262, 642, 401, -199, -569, -523, 23, -393, -1005,
-848, -468, -303, -376, -237, -406, -287, 127, -24,
32, 107, 294, 410, 155, 369, 735, 1157, 894,
228, 386, 77, -355, -358, -547, -555, -441, -222,
-250, -179, 108, 136, 345, 287, 162, -33, -328,
-152, -376, -598, -557, -479, -116, -18, -29, -179,
-58, 306, 129, 206, 138, 67, 337, 284, 462,
255, 127, 401, 97, 33, -124, -90, 194, -120,
-174, -185, -107, -46, -546, -587, -11, 101, -488,
-414, -205, -223, 604, 237, -778, -128, 197, -56,
346, -42, -471, 543, 625, 42, 398, 270, 352,
612, 226, -97, -545, -986, -945, 180, 328, -916,
-487, 557, 1254, 725, -351, 751, 1211, 679, 478,
-325, -431, -253, -48, -569, -1721, -1248, -523, -457,
-1384, -1587, 255, 138, -553, -114, -66, 592, 657,
497, 609, 570, 1263, 1060, 490, 405, 734, 1057,
-32, -239, 210, 68, -104, -819, -495, -219, -666,
-557, -737, -389, -352, -536, -222, -397, 50, 74,
-90, 59, -155, 548, 417, -70, 325, 411, 649,
248, 204, 645, 186, 219, 112, -70, -109, -446,
-278, -542, -702, -586, -596, -352, -617, -462, -216,
-222, 25, -247, -18, 456, 500, 543, 379, 531,
935, 1053, 335, 25, 672, 318, -28, -816, -2456,
-1516, 306, -121, -1684, -1292, 1694, 2342, 236, 60,
2252, 3274, 1020, -188, 1232, 1227, -252, -1636, -1516,
-904, -1924, -2779, -2848, -2019, -1439, -2150, -1854, -683,
309, 118, -339, 1142, 2092, 1617, 1234, 1834, 2791,
2007, 1289, 1677, 1789, 1272, 289, 217, 270, -317,
-828, -1357, -1162, -1102, -1626, -1775, -1529, -903, -1022,
-1267, -600, -113, -67, -274, 162, 766, 571, 349,
359, 1031, 1385, 534, 396, 1299, 1418, 144, -242,
669, 265, -870, -965, -404, -662, -1479, -1228, -777,
-889, -1008, -862, -344, 42, 110, 54, 404, 1167,
1098, 571, 942, 1528, 1010, 363, 766, 584, -897,
-1931, -1206, -25, -637, -1718, -533, 1422, 900, -67,
1286, 2403, 1732, 659, 1009, 1558, 228, -802, -877,
-925, -1503, -2632, -2514, -2009, -2026, -2343, -2281, -906,
-525, -892, -334, 718, 1351, 805, 1249, 2375, 2283,
1903, 1872, 2457, 2290, 1439, 1279, 1129, 860, 43,
-607, -550, -942, -1455, -1851, -1708, -1513, -1901, -1768,
-1313, -951, -920, -816, -144, 112, 93, 227, 495,
816, 757, 554, 582, 1324, 1474, 458, 711, 1433,
798, -86, 145, 626, -513, -1102, -339, -651, -1347,
-1325, -907, -1033, -1323, -828, -642, -444, -113, -33,
411, 789, 969, 827, 980, 1337, 1112, 1075, 727,
374, 27, -933, -1200, -696, -79, -913, -1343, 905,
997, -110, 963, 1732, 1617, 723, 866, 1238, 71,
-341, -839, -1040, -1214, -2238, -2186, -1978, -1808, -2140,
-2121, -819, -845, -894, -189, 489, 830, 605, 1397,
1894, 1789, 1928, 1903, 2286, 2038, 1530, 1394, 1204,
891, -12, -249, -227, -934, -1386, -1439, -1292, -1607,
-1796, -1363, -1213, -1135, -1015, -678, -307, -260, -90,
165, 466, 520, 431, 665, 667, 502, 727, 1115,
806, 367, 1021, 1166, 173, 83, 558, -42, -816,
-562, -464, -1091, -1190, -917, -975, -1025, -924, -692,
-450, -291, -193, 142, 595, 626, 686, 970, 1082,
990, 805, 867, 843, 319, 17, -195, -859, -952,
-270, -365, -1005, -184, 815, 199, 328, 1211, 1085,
894, 789, 633, 384, -89, -576, -976, -924, -1469,
-2074, -1662, -1558, -1801, -1713, -1090, -713, -758, -223,
244, 524, 791, 986, 1443, 1626, 1612, 1636, 1802,
1817, 1365, 1225, 1156, 713, 267, 23, -169, -603,
-894, -988, -1146, -1236, -1325, -1250, -1073, -1022, -940,
-717, -448, -366, -295, -28, 160, 192, 305, 522,
591, 553, 554, 502, 492, 703, 751, 373, 354,
735, 388, -164, 102, 87, -564, -675, -444, -706,
-992, -816, -678, -750, -689, -505, -259, -150, -116,
176, 438, 465, 560, 815, 880, 638, 615, 851,
466, -58, 207, 19, -777, -959, -577, -66, -379,
-678, 444, 939, 363, 565, 1241, 1237, 565, 407,
598, 155, -564, -972, -849, -1078, -1833, -1883, -1422,
-1428, -1737, -1373, -620, -553, -545, 46, 621, 725,
734, 1186, 1600, 1530, 1380, 1579, 1777, 1428, 1053,
1111, 992, 452, 61, 49, -175, -716, -967, -907,
-1025, -1335, -1370, -1115, -1063, -1161, -979, -622, -477,
-475, -219, 156, 240, 243, 503, 785, 734, 650,
795, 826, 557, 315, 452, 605, 60, -329, 190,
153, -550, -494, -135, -389, -811, -593, -314, -597,
-669, -386, -167, -210, -266, 78, 299, 187, 206,
471, 554, 353, 455, 582, 393, 261, 217, 119,
-53, -186, -361, -647, -732, -604, -236, -86, -325,
388, 1004, 599, 871, 1273, 1200, 945, 645, 653,
224, -325, -740, -1025, -1149, -1776, -1976, -1733, -1752,
-1781, -1622, -985, -698, -570, 27, 500, 829, 990,
1336, 1726, 1749, 1748, 1770, 1867, 1667, 1258, 1130,
910, 447, -38, -286, -509, -991, -1326, -1383, -1461,
-1615, -1615, -1384, -1196, -1080, -838, -482, -175, -14,
216, 546, 721, 771, 874, 1017, 977, 800, 710,
674, 476, 139, 50, 37, -290, -491, -485, -547,
-501, -539, -647, -450, -303, -443, -355, 5, 6,
-121, 158, 337, 213, 237, 445, 488, 373, 371,
309, 216, 257, 166, -93, -67, -25, -487, -577,
-237, -757, -1373, -832, -334, -368, -174, 226, 1135,
1470, 1102, 1724, 2099, 1709, 1214, 900, 893, -62,
-961, -1193, -1544, -2098, -2837, -2650, -2327, -2490, -2272,
-1748, -890, -508, -179, 742, 1418, 1785, 1945, 2451,
2852, 2549, 2342, 2263, 2071, 1439, 723, 422, -70,
-780, -1397, -1652, -1795, -2228, -2373, -2126, -1836, -1672,
-1380, -736, -179, 146, 496, 999, 1399, 1449, 1470,
1602, 1565, 1249, 881, 663, 364, -144, -546, -764,
-1006, -1224, -1274, -1323, -1236, -962, -746, -407, -120,
80, 463, 721, 747, 843, 993, 900, 669, 618,
533, 236, -48, -107, -201, -458, -556, -542, -529,
-465, -437, -307, -188, -438, -662, -305, -39, -849,
-1609, -263, 1428, 228, -375, 2323, 3098, 1664, 1745,
2553, 2535, 1112, 77, -77, -756, -1720, -3222, -3587,
-2878, -3439, -4100, -3459, -1755, -1227, -1445, 83, 1728,
2314, 2379, 2965, 3941, 3821, 3266, 2732, 2490, 2154,
755, -389, -760, -1169, -2230, -3181, -2864, -2619, -2872,
-2717, -1946, -927, -441, 68, 922, 1777, 2357, 2318,
2405, 2664, 2408, 1643, 895, 553, -104, -1131, -1816,
-2166, -2399, -2609, -2615, -2416, -1697, -736, -598, 40,
1632, 2103, 1854, 2137, 2549, 2345, 1456, 902, 564,
-128, -821, -1510, -1762, -1629, -1827, -1897, -1407, -639,
-298, -256, 342, 1237, 1516, 948, 881, 1300, 285,
-601, -106, -1005, -2152, -765, 366, -225, 105, 1800,
2667, 2263, 2358, 2609, 1902, 1309, 383, -1178, -2028,
-2431, -3415, -4550, -4205, -3344, -3426, -2845, -1330, 27,
949, 1939, 3073, 3622, 4112, 4210, 3486, 2987, 2535,
1477, 2, -942, -1447, -2432, -3129, -3225, -3139, -2816,
-2288, -1557, -826, 190, 1229, 1697, 2171, 2713, 2848,
2494, 1998, 1574, 810, -38, -777, -1573, -2133, -2414,
-2554, -2640, -2453, -1750, -916, -403, 165, 1489, 2218,
1984, 2950, 3358, 1847, 1268, 1173, 3, -1163, -1755,
-2176, -2690, -2518, -1972, -1908, -1192, 22, 618, 1019,
1695, 2107, 1666, 1636, 2010, 1092, 94, -263, -1376,
-1768, -457, -1731, -5151, -3981, 955, 2337, -400, 750,
5303, 6160, 5173, 3503, 1711, 2310, 1724, -1812, -6147,
-6323, -3861, -5816, -7464, -5990, -3493, -1001, 97, 1540,
2770, 5367, 7786, 5898, 4386, 4606, 4302, 1937, -1267,
-2217, -3112, -3790, -4215, -5298, -4526, -2551, -773, -546,
-32, 2599, 3965, 3955, 3518, 3064, 2937, 2131, 893,
-1169, -2477, -2268, -2844, -3763, -3940, -2899, -1616, -1023,
-181, 607, 1586, 2617, 3023, 3070, 2914, 2991, 2172,
623, -310, -1374, -2183, -2752, -3338, -3214, -2667, -1629,
-816, -278, 799, 1765, 2420, 2415, 2095, 1811, 1020,
509, 175, -542, -1258, -1676, -1622, -1721, -1740, -1068,
-852, -2054, -2617, 929, 5290, 4343, 2568, 4316, 5511,
5146, 3002, -412, -2276, -2717, -3204, -6159, -8896, -7253,
-4846, -3602, -2969, -1694, 1159, 3995, 6267, 6021, 4739,
5447, 5599, 3839, 720, -1621, -2390, -3349, -4039, -4975,
-5170, -3395, -1315, 139, 701, 2083, 4056, 4738, 4561,
3559, 2429, 1602, 507, -1149, -3232, -4042, -3810, -3787,
-3580, -3000, -1677, -99, 1293, 2238, 2416, 2580, 2706,
3382, 4128, 2082, -743, -502, -239, -2169, -3426, -3827,
-3585, -2181, -999, -751, -658, 1255, 3304, 2794, 1816,
1619, 1680, 1241, 8, -1322, -2153, -1180, -321, -1526,
-2302, -1478, -72, 363, -230, -426, -82, -103, -77,
2541, 4581, 3326, 3631, 4197, 3161, 2512, 785, -1567,
-3641, -4417, -4491, -6115, -6640, -5357, -3878, -2181, -630,
1004, 2187, 3811, 5725, 5413, 4276, 3600, 2869, 1625,
-298, -1819, -3206, -3817, -3270, -3022, -2734, -1844, -175,
1408, 2225, 2944, 3158, 3145, 3144, 2440, 1057, -421,
-1236, -1857, -2701, -3165, -3225, -2780, -1898, -937, -298,
-206, 87, 1419, 3719, 4816, 3523, 2068, 1482, 881,
-351, -2295, -3896, -4219, -3195, -2016, -1842, -1410, 9,
1300, 1732, 1636, 1591, 1442, 1292, 1060, 905, 951,
40, -693, -286, -671, -1437, -1264, -383, 371, 148,
623, 1538, -1352, -5744, -4090, 1312, 2921, 1799, 1479,
2488, 5347, 7009, 4782, 210, -1976, -191, -607, -4050,
-6603, -6961, -5451, -3585, -2599, -3395, -3123, 320, 3569,
4271, 3628, 3813, 4892, 5331, 4548, 2148, -569, -1301,
-1040, -1894, -3494, -4162, -3511, -2160, -805, -235, -93,
737, 2313, 3162, 2711, 2027, 1466, 942, 786, 684,
-712, -2331, -1677, -1330, -2395, -1867, -539, -52, -62,
434, 820, 233, 735, 1092, -206, -628, -8, 124,
-546, -681, -390, -814, -535, 90, -5, -80, 381,
1054, 919, 810, 759, 221, 373, 445, -358, -1278,
-1213, 15, 517, 240, 469, 618, 117, -366, -1994,
-3771, -679, 3229, 1528, -90, 1057, 1545, 2429, 2774,
1183, -1064, -1624, 94, -409, -2670, -2848, -2292, -1846,
-1057, -812, -1316, -1070, 454, 1003, 468, 747, 1275,
1268, 1598, 1792, 975, 808, 942, 410, 441, -32,
-622, -261, -436, -603, -474, -709, -641, -368, -305,
-310, -242, -165, -29, 114, 107, 193, 187, 158,
287, 194, 100, 51, -28, -31, -70, -19, -33,
-143, -165, -225, -240, -244, -273, -283, -269, -118,
-7, 34, 68, 94, 196, 267, 258, 202, 136,
102, 29, -72, -179, -281, -298, -274, -227, -161,
-15, 162, 227, 236, 268, 332, 335, 237, 132,
8, -77, -130, -179, -220, -267, -270, -218, -159,
-65, 46, 111, 144, 187, 226, 197, 129, 56,
-14, -93, -179, -259, -280, -256, -236, -198, -168,
-134, -41, 51, 66, 42, 50, 76, 107, 156,
168, 144, 124, 110, 76, -5, -50, -72, -130,
-219, -284, -253, -181, -118, -75, -56, -9, 77,
144, 175, 176, 206, 247, 240, 187, 85, -15,
-79, -158, -240, -274, -283, -257, -178, -96, -36,
44, 131, 196, 206, 180, 175, 134, 61, -26,
-116, -179, -195, -188, -215, -250, -209, -137, -97,
-62, -22, -5, 32, 84, 114, 93, 100, 142,
135, 84, 33, 3, -8, -9, -22, -60, -79,
-50, -5, 17, 2, -1, 6, 23, 40, 43,
25, 26, 42, 29, -15, -55, -104, -134, -147,
-203, -244, -220, -175, -118, -58, -5, 44, 95,
136, 135, 134, 132, 115, 60, -4, -46, -72,
-110, -124, -135, -159, -145, -93, -46, -21, 44,
107, 141, 151, 132, 111, 110, 100, 44, -15,
-43, -67, -92, -107, -110, -96, -80, -58, -15,
13, 44, 67, 68, 42, 23, 19, -15, -49,
-86, -101, -121, -131, -141, -159, -165, -152, -120,
-92, -52, 13, 67, 107, 139, 151, 141, 117,
76, 36, 10, -16, -35, -73, -75, -56, -46,
-49, -58, -46, -15, 17, 40, 43, 53, 59,
49, 34, 12, -9, -19, -43, -79, -116, -130,
-134, -123, -100, -96, -77, -48, -12, 30, 44,
59, 60, 47, 40, 25, 20, 22, 15, 6,
-14, -24, -29, -35, -36, -43, -45, -41, -33,
-18, -9, 3, 6, -2, -4, 3, 9, 17,
20, 17, 10, 3, 5, 3, -12, -24, -31,
-36, -45, -52, -48, -55, -59, -60, -55, -53,
-53, -63, -65, -60, -50, -43, -48, -35, -24,
-25, -35, -41, -39, -33, -16, -2, 12, 27,
42, 57, 68, 66, 61, 51, 44, 25, 0,
-9, -31, -49, -53, -52, -36, -39, -41, -7,
27, 25, 15, 27, 32, 23, 20, 5, -18,
-26, -36, -39, -50, -62, -63, -62, -55, -50,
-41, -33, -25, -14, -1, 5, 9, 15, 5,
-7, -18, -26, -33, -46, -56, -73, -82, -66,
-67, -50, -38, -22, -7, 17, 30, 37, 36,
29, 19, 13, 12, -7, -25, -41, -42, -38,
-49, -49, -35, -26, -8, 6, 0, 0, -2,
-5, -1, -4, -2, -12, -12, -15, -16, -19,
-24, -35, -39, -48, -46, -35, -25, -4, 13,
25, 27, 3, 2, 30, 27, -1, -1, -16,
-25, -18, -36, -42, -48, -62, -56, -48, -53,
-46, -36, -26, -15, -8, 12, 23, 29, 25,
13, 5, 5, 0, -16, -12, -19, -25, -24,
-29, -29, -39, -33, -25, -33, -19, -5, -11,
-22, -24, -18, -25, -32, -39, -39, -32, -41,
-36, -43, -39, -18, -11, -7, -12, -8, -4,
-7, 0, 5, -4, 2, 6, -2, -5, -19,
-25, -29, -33, -39, -36, -25, -24, -12, -7,
5, 13, 16, 34, 29, 3, -1, -1, -14,
-16, -18, -19, -19, -26, -22, -32, -26, -26,
-28, -14, -11, -5, -2, -4, 0, -9, -15,
-16, -28, -28, -32, -48, -56, -52, -55, -56,
-56, -48, -41, -31, -22, -8, -2, -1, 9,
10, -1, 3, -7, -1, 0, 2, -1, -11,
-16, -28, -38, -45, -41, -42, -24, -15, 5,
9, -9, -19, -62, -192, -328, 125, 1057, 514,
-579, 277, -138, -1499, -18, 595, -276, 67, -332,
-853, 761, 1310, -301, -237, 513, 762, 540, -164,
57, -1481, -1500, 396, -1350, -365, 1445, -933, 136,
1200, -686, 876, 1901, 166, -349, -244, -198, -1073,
-1085, 187, -481, 294, 893, -877, -4, 659, -7,
-7, -178, 73, -62, 463, 223, -523, 500, 209,
318, 199, -818, 398, -28, -508, 320, -635, -127,
226, -688, 83, -65, -426, 214, -318, -297, 425,
134, 29, 311, -21, 115, 316, -417, -84, 127,
-294, 407, 219, -470, -110, 221, -148, -176, 231,
-189, 180, 194, -562, 172, 71, 110, 90, -645,
107, -58, 267, 646, -467, 349, 628, -366, -205,
-16, -235, -386, -322, -31, 22, 54, 316, -226,
219, 689, -186, 272, -288, -679, 37, -351, 17,
-284, -130, 95, -312, 347, -15, 563, 439, -696,
54, -113, -233, -212, -246, 190, -33, 91, -76,
-106, 192, 88, 117, 196, 432, 56, -100, 279,
0, -65, -243, -412, -49, -193, -361, 121, 104,
-280, 102, -89, -120, 488, -7, -145, 194, 261,
-270, -314, 396, -185, -93, 271, 81, 134, -317,
13, -158, -549, 445, 388, -53, -2, -550, -332,
88, -237, -131, -45, 36, 386, -128, -243, 90,
145, 686, -259, -511, 217, -213, 422, -318, -277,
809, -339, 131, -73, -577, 354, 296, 527, -158,
-488, -89, -196, 203, -252, -86, 146, -242, 264,
-39, 56, 287, -470, -624, -666, 71, 632, 221,
292, 189, 390, 39, -600, -80, -182, -79, 182,
44, -33, -700, 101, 319, -417, 895, 459, -155,
550, -325, -645, -325, -94, -491, -460, 284, -193,
182, 669, 325, -58, 122, 520, -504, -423, 169,
169, 360, -165, -195, -404, -161, 42, -1156, -581,
449, 434, 301, -41, 405, 544, 298, -87, -716,
46, 322, -460, -590, -164, -28, -530, -140, 141,
138, 776, 514, -121, -461, 380, 1170, 257, -32,
10, -174, -114, -369, -771, -925, 148, 571, -413,
-89, 539, 17, -87, 6, -716, -392, 683, 124,
-436, 233, 473, 66, -198, -329, -478, -26, 588,
0, -475, 362, 483, -189, -145, -310, -516, -4,
63, -461, -269, 413, 360, 507, 805, 291, 447,
580, -383, -723, -312, -141, -270, -638, -824, -416,
-89, -277, -126, 476, 1006, 1010, 898, 1040, 320,
-185, -148, -794, -989, -1019, -972, -450, -210, 17,
-31, 83, 90, -138, 415, 219, -162, 597, 1156,
1003, 108, -376, -202, -291, -216, -644, -1095, -404,
282, 248, -150, 145, 371, -104, 305, 261, -512,
-267, -101, -549, -553, 59, 260, 350, 557, -114,
-223, 495, 388, 278, 414, 177, 22, 23, 29,
-12, -65, 90, -250, -756, -60, 533, 142, -76,
-157, -181, 193, -77, -1185, -1054, 54, 333, -124,
-750, -431, 700, 1449, 1017, -877, -1236, 415, 277,
-586, -126, 373, 326, 85, 410, 132, -361, 185,
-293, -440, 915, 1344, 155, -698, 262, 519, -270,
-431, -678, -693, -556, -294, 78, -182, -325, -478,
-675, -55, 248, -404, -395, 415, 907, 614, -206,
-14, 699, 387, -135, -276, 192, 759, 456, -263,
-420, 340, 585, -128, -208, -114, 111, 781, 93,
-675, -99, -7, -560, -835, -647, -638, -199, 853,
291, -198, 666, 439, -127, -266, -457, -229, -134,
-103, -678, -1277, -120, 782, 461, -101, 124, 985,
421, -76, -128, -329, 493, 806, 653, 363, 349,
687, -593, -945, -67, -495, -1020, -1158, -869, -436,
359, 757, -753, -1234, 42, 680, 217, -504, -137,
609, 934, 333, -750, -252, 867, 632, -327, -104,
951, 1184, 548, -529, -618, 565, 969, -179, -967,
-184, 680, 735, -87, -1227, -985, 125, -7, -1102,
-794, 676, 669, -28, -332, -1333, -1237, 513, 1204,
335, 318, 1758, 420, -2489, -1776, 76, -359, -1987,
-1489, 925, 2204, 2995, 2332, 306, 1254, 3049, 1171,
-1367, -1189, -468, -1397, -2357, -2116, -1536, -1060, -750,
-790, -788, -75, 638, 683, 699, 383, 91, 642,
911, 258, -376, -362, -271, -167, 241, 296, 308,
650, 945, 999, 883, 953, 946, 418, -453, -808,
-461, -706, -1258, -988, -434, -52, 66, -123, 151,
679, 428, -481, -509, 771, 714, -512, -1035, -1486,
-726, -31, -1281, -1212, 1309, 2854, -717, -5541, -2528,
4411, 4364, -965, -2135, 1605, 6128, 6542, 901, -2670,
955, 3056, -1970, -6739, -4304, -1117, -2714, -3398, -2501,
-995, 1508, 2242, 360, -608, 1734, 2558, 452, -246,
119, -35, -26, 448, 260, 60, 1385, 1602, 292,
541, 1344, 673, -386, -482, -556, -744, -250, -352,
-906, -236, 1023, 238, -1847, -1223, 598, 466, 107,
156, 44, 587, 1174, 431, -1296, -1827, -647, -169,
-780, -1298, -1319, -43, 2335, 3246, 1074, -2402, -6047,
-7079, -1059, 4745, 1480, -734, 4822, 9405, 9224, 5118,
-67, -1939, -822, -1629, -6357, -8638, -5218, -2163, -2139,
-2847, -1910, 727, 2801, 3025, 1442, 682, 1915, 2640,
1630, -328, -1203, -831, -587, -253, 102, 145, 388,
1416, 1997, 1157, 750, 902, 285, -634, -968, -935,
-1489, -1441, -559, -307, 98, 255, -79, 393, 757,
383, -92, 224, 977, 612, -250, -768, -1037, -1196,
-1434, -1257, -528, 53, -83, 128, 1382, 1759, 1323,
-341, -8125, -13054, -157, 15976, 6087, -6341, 5951, 14103,
11032, 4159, -6743, -10446, -3027, 1775, -9133, -16626, -6476,
3326, 1674, -3781, -1209, 4156, 6635, 7114, 2682, -1791,
1685, 5147, 898, -4427, -3907, -1369, -1268, -375, -530,
-1009, 1666, 3818, 2716, 1439, 2971, 2807, -334, -836,
548, -1022, -3429, -2739, -1867, -1887, -494, 346, 44,
972, 2167, 1438, -239, -396, -202, -971, -1043, -826,
-696, 13, 401, -335, -805, 129, 226, -515, 410,
1208, 717, 88, -5350, -13627, -8527, 10641, 17424, -276,
-2734, 12987, 14423, 8109, -1717, -11484, -8209, -128, -4410,
-16396, -14255, -1867, 5154, 1136, -2316, 3394, 8584, 10376,
6841, -603, -2497, 2138, 1684, -4646, -7354, -4663, -2188,
-675, 945, 795, 1892, 5103, 6024, 3291, 1198, 1634,
500, -1851, -2889, -3691, -3739, -1530, 228, -833, -215,
2780, 3205, 1946, 1246, 642, -103, -849, -1823, -3171,
-2957, -1587, -1554, -1490, -322, 710, 1122, 1530, 1998,
1508, 689, 476, 669, 754, -720, -4025, -6837, -7325,
-6367, -2698, 6174, 12964, 9699, 7442, 11835, 10703, 3935,
-2704, -7836, -10023, -10369, -10922, -12144, -10017, -3153, 3001,
3870, 4798, 9021, 10308, 7988, 4810, 1285, -1649};
+54
View File
@@ -0,0 +1,54 @@
/*
mksine.c
David Rowe
10 Nov 2010
Creates a file of sine wave samples.
*/
#include <assert.h>
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define TWO_PI 6.283185307
#define FS 8000.0
int main(int argc, char *argv[]) {
FILE *f;
int i, n;
float freq, length;
short *buf;
float amp = 1E4;
if (argc < 4) {
printf("usage: %s outputFile frequencyHz lengthSecs [PeakAmp]\n", argv[0]);
exit(1);
}
if (strcmp(argv[1], "-") == 0) {
f = stdout;
} else if ((f = fopen(argv[1], "wb")) == NULL) {
fprintf(stderr, "Error opening output file: %s: %s.\n", argv[3],
strerror(errno));
exit(1);
}
freq = atof(argv[2]);
length = atof(argv[3]);
if (argc == 5) amp = atof(argv[4]);
n = length * FS;
buf = (short *)malloc(sizeof(short) * n);
assert(buf != NULL);
for (i = 0; i < n; i++) buf[i] = amp * cos(freq * i * (TWO_PI / FS));
fwrite(buf, sizeof(short), n, f);
fclose(f);
free(buf);
return 0;
}
+68
View File
@@ -0,0 +1,68 @@
#! /bin/bash
# ofdm_check
#
# A series of checks of the ofdm functions, mostly decode.
#
# This uses ofdm_mod to supply test data to ofdm_demod and mostly
# assumes that the encode function is correct.
# Define macros to (later) allow testing alternate versions.
alias OFDM_MOD=ofdm_mod
alias OFDM_DEMOD=ofdm_demod
shopt -s expand_aliases
# PATH
PATH=$PATH:../src
PASS=1
###############################
echo
echo "Simple test, plain, ideal"
OFDM_MOD --in /dev/zero --testframes 100 |
OFDM_DEMOD --out /dev/null --testframes --verbose 1 2> tmp
cat tmp
p1=$(grep '^BER\.*: 0.000' tmp | wc -l)
p2=$(grep '^BER2\.*: 0.000' tmp | wc -l)
if [[ $p1 -eq 1 && $p2 -eq 1 ]]; then echo "OK"; else echo "BAD"; PASS=0; fi
###############################
echo
echo "Simple test, plain, AWGN"
OFDM_MOD --in /dev/zero --testframes 100 |
cohpsk_ch - - -20 -Fs 8000 -f -5 |
OFDM_DEMOD --out /dev/null --testframes --verbose 1 2>tmp
cat tmp
n=$(grep '^BER\.*:' tmp | cut -d ' ' -f 2)
p1=$(echo $n '<=' 0.10 | bc)
n=$(grep '^BER2\.*:' tmp | cut -d ' ' -f 2)
p2=$(echo $n '<=' 0.10 | bc)
if [[ $p1 -eq 1 && $p2 -eq 1 ]]; then echo "OK"; else echo "BAD"; PASS=0; fi
###############################
echo
echo "Simple test, LDPC, ideal"
OFDM_MOD --in /dev/zero --ldpc --testframes 100 |
OFDM_DEMOD --out /dev/null --ldpc --testframes --verbose 1 2>tmp
cat tmp
p1=$(grep '^BER\.*: 0.000' tmp | wc -l)
p2=$(grep '^Coded BER: 0.000' tmp | wc -l)
if [[ $p1 -eq 1 && $p2 -eq 1 ]]; then echo "OK"; else echo "BAD"; PASS=0; fi
###############################
echo
echo "Simple test, LDPC, AWGN"
OFDM_MOD --in /dev/zero --ldpc --testframes 100 |
cohpsk_ch - - -20 -Fs 8000 -f -5 --fading_dir ../../build_linux/unittest |
OFDM_DEMOD --out /dev/null --ldpc --testframes --verbose 1 2>tmp
cat tmp
n=$(grep '^BER\.*:' tmp | cut -d ' ' -f 2)
p1=$(echo $n '<=' 0.10 | bc)
n=$(grep '^Coded.*BER\.*:' tmp | cut -d ' ' -f 3)
p2=$(echo $n '<=' 0.01 | bc)
if [[ $p1 -eq 1 && $p2 -eq 1 ]]; then echo "OK"; else echo "BAD"; PASS=0; fi
echo
if [[ $PASS == 1 ]]; then echo "PASSED"; else echo "FAILED"; fi
+12
View File
@@ -0,0 +1,12 @@
#!/usr/bin/env bash
#
# David June 2019
# Tests 700D OFDM modem fading channel performance, using a simulated channel
results=$(mktemp)
fading_dir=$1
# BER should be around 4% for this test (it's better for larger interleavers but no one uses interleaving in practice)
ofdm_mod --in /dev/zero --ldpc 1 --testframes 60 --txbpf | ch - - --No -24 -f -10 --mpp --fading_dir $fading_dir | ofdm_demod --out /dev/null --testframes --verbose 2 --ldpc 1 2> $results
cat $results
cber=$(cat $results | sed -n "s/^Coded BER.* \([0-9..]*\) Tbits.*/\1/p")
python3 -c "import sys; sys.exit(0) if $cber<=0.05 else sys.exit(1)"
+24
View File
@@ -0,0 +1,24 @@
#!/usr/bin/env bash
#
# ofdm_phase_est_bw.sh
# David August 2019
# Tests 2020 OFDM modem phase est bandwidth mode option. Locking the
# phase est bandwidth to "high" is useful for High SNR channels with
# fast fading or high phase noise. In this test we show that with
# high bandwidth phase est mode, the BER is < 5% for the "--mpp" (1
# Hz fading) channel model on a fairly high SNR channel.
#
# To run manually outside of ctest:
# $ cd codec2/unittest
# $ PATH=$PATH:../build_linux/src ./ofdm_phase_est_bw.sh
results=$(mktemp)
fading_dir=$1
# BER should be < 5% for this test
ofdm_mod --in /dev/zero --testframes 300 --mode 2020 --ldpc --verbose 0 | \
ch - - --No -40 -f 10 --ssbfilt 1 --mpp --fading_dir $fading_dir | \
ofdm_demod --out /dev/null --testframes --mode 2020 --verbose 2 --ldpc --bandwidth 1 2> $results
cat $results
cber=$(cat $results | sed -n "s/^Coded BER.* \([0-9..]*\) Tbits.*/\1/p")
python3 -c "import sys; sys.exit(0) if $cber<=0.05 else sys.exit(1)"
+30
View File
@@ -0,0 +1,30 @@
#!/usr/bin/env bash
# Shell script version of ofdm_time_sync()
# David June 2019
# Tests ofdm modem sync time, using real, off air files
onerun=$(mktemp)
results=$(mktemp)
# bunch of runs at different time offsets into a sample off air file with low SNR and semi-staionary fading
for start_secs in `seq 0 29`;
do
if [ "$1" = "700D" ]; then
ofdm_demod --in ../wav/vk2tpm_004.wav --out /dev/null --verbose 2 --ldpc \
--start_secs $start_secs --len_secs 5 2>/dev/null > $onerun
fi;
if [ "$1" = "2020" ]; then
ofdm_demod --mode 2020 --in ../wav/david4.wav --out /dev/null --verbose 2 --ldpc \
--start_secs $start_secs --len_secs 5 2>/dev/null > $onerun
fi;
[ ! $? -eq 0 ] && { echo "error running ofdm_demod"; exit 1; }
cat $onerun | sed -n "s/time_to_sync: \([0-9..]*\)/\1/p" >> $results
done
# a pass is we never take longer than 5 secs to sync (mean is much smaller)
python3 -c "
import sys; import numpy as np
x=np.loadtxt(\"$results\")
fails=sum(x == -1)
print(\"fails: %d mean: %5.2f var: %5.2f \" % (fails, np.mean(x), np.var(x)))
sys.exit(0) if fails==0 else sys.exit(1)
"
+149
View File
@@ -0,0 +1,149 @@
# Makefile
# Dec 2022
#
# Automates PER/BER curve generation for raw data mode:
#
# 1. Compare "ch" noise injection/SNR measurement against reference Octave Tx
# 2. Compare C Tx against Octave Tx (with and without compression)
# 3. Plot curves for SNR estimated from C Rx against actual SNR
# 4. Plot AWGN/PER C Tx curves for end user documentation
#
# usage:
#
# 1. Create 20 minutes of MPP fading samples:
# $ cd codec2/octave/
# $ octave-cli
# octave:24> pkg load signal
# octave:24> time_secs=60*20
# octave:26> ch_fading("~/codec2/build_linux/unittest/fast_fading_samples.float", 8000, 1.0, 8000*time_secs)
#
# 2. Run scripts:
#
# $ make
SHELL := /bin/bash
CODEC2 := $(HOME)/codec2
all: test \
octave_ch_noise_awgn.png octave_c_tx_awgn.png octave_c_tx_comp_awgn.png \
octave_ch_noise_mpp.png octave_c_tx_mpp.png octave_c_tx_comp_mpp.png \
snrest_snr_ctx.png snrest_snr_ctxc.png \
c_tx_comp.png c_tx_comp_thruput.png
clean:
rm -f *.txt *.png *.raw
# run this first, traps common CML setup error
test:
source snr_curves.sh; test_ldpc
# subset of files generated, but enough to set up Makefile dependencies
snr_oct = snr_oct_datac0_awgn.txt snr_oct_datac1_awgn.txt snr_oct_datac3_awgn.txt
snr_ch = snr_ch_datac0_awgn.txt snr_ch_datac1_awgn.txt snr_ch_datac3_awgn.txt
snr_ctx = snr_ctx_datac0_awgn.txt snr_ctx_datac1_awgn.txt snr_ctx_datac3_awgn.txt
snr_ctxc = snr_ctxc_datac0_awgn.txt snr_ctxc_datac3_awgn.txt
snr_oct_mpp = snr_oct_datac0_mpp.txt snr_oct_datac1_mpp.txt snr_oct_datac3_mpp.txt
snr_ch_mpp = snr_ch_datac0_mpp.txt snr_ch_datac1_mpp.txt snr_ch_datac3_mpp.txt
snr_ctx_mpp = snr_ctx_datac0_mpp.txt snr_ctx_datac1_mpp.txt snr_ctx_datac3_mpp.txt
snr_ctxc_mpp = snr_ctxc_datac0_mpp.txt snr_ctxc_datac3_mpp.txt
$(snr_oct):
source snr_curves.sh; generate_octave_tx_data datac0 awgn
source snr_curves.sh; generate_octave_tx_data datac1 awgn
source snr_curves.sh; generate_octave_tx_data datac3 awgn
$(snr_oct_mpp):
source snr_curves.sh; generate_octave_tx_data datac0 mpp
source snr_curves.sh; generate_octave_tx_data datac1 mpp
source snr_curves.sh; generate_octave_tx_data datac3 mpp
$(snr_ch):
source snr_curves.sh; generate_ch_data datac0 awgn
source snr_curves.sh; generate_ch_data datac1 awgn
source snr_curves.sh; generate_ch_data datac3 awgn
$(snr_ch_mpp):
source snr_curves.sh; generate_ch_data datac0 mpp
source snr_curves.sh; generate_ch_data datac1 mpp
source snr_curves.sh; generate_ch_data datac3 mpp
# C without compression
$(snr_ctx):
source snr_curves.sh; generate_snrest_v_snr_data datac0 awgn
source snr_curves.sh; generate_snrest_v_snr_data datac1 awgn
source snr_curves.sh; generate_snrest_v_snr_data datac3 awgn
source snr_curves.sh; generate_snrest_v_snr_data datac4 awgn
source snr_curves.sh; generate_snrest_v_snr_data datac13 awgn
source snr_curves.sh; generate_snrest_v_snr_data datac14 awgn
$(snr_ctx_mpp):
source snr_curves.sh; generate_snrest_v_snr_data datac0 mpp
source snr_curves.sh; generate_snrest_v_snr_data datac1 mpp
source snr_curves.sh; generate_snrest_v_snr_data datac3 mpp
source snr_curves.sh; generate_snrest_v_snr_data datac4 mpp
source snr_curves.sh; generate_snrest_v_snr_data datac13 mpp
source snr_curves.sh; generate_snrest_v_snr_data datac14 mpp
# C with compression
$(snr_ctxc):
source snr_curves.sh; generate_snrest_v_snr_data datac0 awgn 1
source snr_curves.sh; generate_snrest_v_snr_data datac1 awgn 1
source snr_curves.sh; generate_snrest_v_snr_data datac3 awgn 1
source snr_curves.sh; generate_snrest_v_snr_data datac4 awgn 1
source snr_curves.sh; generate_snrest_v_snr_data datac13 awgn 1
source snr_curves.sh; generate_snrest_v_snr_data datac14 awgn 1
$(snr_ctxc_mpp):
source snr_curves.sh; generate_snrest_v_snr_data datac0 mpp 1
source snr_curves.sh; generate_snrest_v_snr_data datac1 mpp 1
source snr_curves.sh; generate_snrest_v_snr_data datac3 mpp 1
source snr_curves.sh; generate_snrest_v_snr_data datac4 mpp 1
source snr_curves.sh; generate_snrest_v_snr_data datac13 mpp 1
source snr_curves.sh; generate_snrest_v_snr_data datac14 mpp 1
# Octave and C curves should be on top of each other, indicating Octave
# and ch noise injection/SNR measurement are equivalent
octave_ch_noise_awgn.png: $(snr_oct) $(snr_ch)
echo "snr_curves_plot; octave_ch_noise_print('awgn'); quit" | \
octave-cli -p $(CODEC2)/octave
octave_ch_noise_mpp.png: $(snr_oct_mpp) $(snr_ch_mpp)
echo "snr_curves_plot; octave_ch_noise_print('mpp'); quit" | \
octave-cli -p $(CODEC2)/octave
# Octave Tx and C Tx curves should be on top of each other
octave_c_tx_awgn.png: $(snr_oct) $(snr_ctx)
echo "snr_curves_plot; octave_c_tx_print('awgn'); quit" | \
octave-cli -p $(CODEC2)/octave
octave_c_tx_mpp.png: $(snr_oct_mpp) $(snr_ctx_mpp)
echo "snr_curves_plot; octave_c_tx_print('mpp'); quit" | \
octave-cli -p $(CODEC2)/octave
# Octave Tx and C Tx (compressed) curves should be close, but C may be 1dB
# poorer
octave_c_tx_comp_awgn.png: $(snr_oc) $(snr_ctxc)
echo "snr_curves_plot; octave_c_tx_comp_print('awgn'); quit" | \
octave-cli -p $(CODEC2)/octave
octave_c_tx_comp_mpp.png: $(snr_oct_mpp) $(snr_ctxc_mpp)
echo "snr_curves_plot; octave_c_tx_comp_print('mpp'); quit" | \
octave-cli -p $(CODEC2)/octave
# combined AWGN and MPP from C Tx (compressed) - what end users would run
c_tx_comp.png: $(snr_ctxc) $(snr_ctxc_mpp)
echo "snr_curves_plot; c_tx_comp_print; quit" | \
octave-cli -p $(CODEC2)/octave
# Curves of SNR estimates from C Rx compared to actual SNR, useful for "gear shifting"
snrest_snr_ctx.png: $(snr_ctx)
echo "snr_curves_plot; snrest_snr_print('ctx', 'awgn'); quit" | \
octave-cli -p $(CODEC2)/octave
snrest_snr_ctxc.png: $(snr_ctxc)
echo "snr_curves_plot; snrest_snr_print('ctxc', 'awgn'); quit" | \
octave-cli -p $(CODEC2)/octave
# Throughput of payload data in bits/s of modes against SNR
c_tx_comp_thruput.png: $(snr_ctxc) $(snr_ctxc_mpp)
echo "snr_curves_plot; c_tx_comp_thruput_print; quit" | \
octave-cli -p $(CODEC2)/octave
+191
View File
@@ -0,0 +1,191 @@
# snr_curves.sh
#
# Library of bash functions to generate data for SNR curves.
#
# testing a function example:
# $ bash -c "source ./snr_curves.sh; generate_octave_tx_data datac0 awgn"
set -x
PATH=${PATH}:${HOME}/codec2/build_linux/src
CODEC2=${HOME}/codec2
FADING_DIR=${CODEC2}/build_linux/unittest
snr_list='-5 -4 -3 -2 0 1 2 4'
No_list='-13 -14 -15 -16 -18 -20 -22 -24 -26'
Nbursts_awgn=20
Nbursts_mpp=100
# Octave Tx injects noise and is source of truth for SNR, measure BER/PER v SNR
function generate_octave_tx_data {
mode=$1
channel=$2
Nbursts=$Nbursts_awgn
snr_nudge=0
if [ "$channel" == "mpp" ]; then
Nbursts=$Nbursts_mpp
snr_nudge=4
fi
rx_log=$(mktemp)
i=1
rm -f snr_oct_${mode}_${channel}*.txt
rm -f ber_oct_${mode}_${channel}*.txt
rm -f per_oct_${mode}_${channel}*.txt
for snr in $snr_list
do
snr_adj=$((${snr}+${snr_nudge}))
echo "warning ('off', 'Octave:data-file-in-path');
ofdm_ldpc_tx('test_${mode}.raw','${mode}',1,${snr_adj},'${channel}','bursts',${Nbursts},'crc');
quit" | DISPLAY="" octave-cli -p ${CODEC2}/octave
freedv_data_raw_rx --testframes $mode test_${mode}.raw /dev/null 2>${rx_log} -v
BERmeas=$(cat ${rx_log} | grep 'BER......:' | cut -d' ' -f2)
PERmeas=$(cat ${rx_log} | grep 'Coded FER' | cut -d' ' -f3)
echo ${snr_adj} >> snr_oct_${mode}_${channel}.txt
echo ${BERmeas} >> ber_oct_${mode}_${channel}.txt
echo ${PERmeas} >> per_oct_${mode}_${channel}.txt
i=$((i+1))
done
echo 0 > offset_oct_${mode}_${channel}.txt
}
# ch injects noise and is source of truth for SNR, measure BER/PER v SNR
# Octave Tx
function generate_ch_data {
mode=$1
channel=$2
ch_multipath=''
Nbursts=$Nbursts_awgn
snr_nudge=0
if [ "$channel" == "mpp" ]; then
ch_multipath='--mpp'
Nbursts=$Nbursts_mpp
snr_nudge=4
fi
octave_log=$(mktemp)
ch_log=$(mktemp)
rx_log=$(mktemp)
i=1
rm -f snr_ch_${mode}_${channel}*.txt
rm -f ber_ch_${mode}_${channel}*.txt
rm -f per_ch_${mode}_${channel}*.txt
for No in $No_list
do
No_adj=$((${No}-${snr_nudge}))
echo "warning ('off', 'Octave:data-file-in-path');
ofdm_ldpc_tx('test_${mode}.raw','${mode}',1,100,'awgn','bursts',${Nbursts},'crc');
quit" | DISPLAY="" octave-cli -p ${CODEC2}/octave 1>${octave_log}
SNRoffset=$(cat ${octave_log} | grep 'Burst offset:' | cut -d' ' -f5)
ch test_${mode}.raw - --No $No_adj ${ch_multipath} --fading_dir ${FADING_DIR} 2>>${ch_log} | \
freedv_data_raw_rx --testframes $mode - /dev/null -v 2>${rx_log}
BERmeas=$(cat ${rx_log} | grep 'BER......:' | cut -d' ' -f2)
PERmeas=$(cat ${rx_log} | grep 'Coded FER' | cut -d' ' -f3)
echo ${BERmeas} >> ber_ch_${mode}_${channel}.txt
echo ${PERmeas} >> per_ch_${mode}_${channel}.txt
i=$((i+1))
# trap not enough fading file samples (with mpp)
grep "Fading file finished" ${ch_log}
if [ $? -eq 0 ]; then
cat ${ch_log}
exit 1
fi
done
echo ${SNRoffset} > offset_ch_${mode}_${channel}.txt
SNRch=$(cat ${ch_log} | grep SNR3k | tr -s ' ' | cut -d' ' -f3)
echo ${SNRch} > snr_ch_${mode}_${channel}.txt
}
# ch injects noise and is source of truth for SNR, measure BER/PER v SNR and
# SNR estimates v SNR from rx, C Tx
function generate_snrest_v_snr_data {
mode=$1
channel=$2
snr_nudge=0
aNo_list=$No_list
# nudge SNR test range to get meaningful results for these tests
if [ "$mode" == "datac1" ]; then
snr_nudge=4
fi
if [[ "$mode" == "datac4" || "$mode" == "datac13" || "$mode" == "datac14" ]]; then
snr_nudge=-6
fi
ch_multipath=''
Nbursts=$Nbursts_awgn
if [ "$channel" == "mpp" ]; then
ch_multipath='--mpp'
Nbursts=$Nbursts_mpp
snr_nudge=$((${snr_nudge}+4))
fi
clip=0
id='ctx'
if [ "$#" -eq 3 ]; then
clip=$3
id='ctxc'
snr_nudge=$((${snr_nudge}-4))
fi
tx_log=$(mktemp)
ch_log=$(mktemp)
rx_log=$(mktemp)
i=1
rm -f snrest_${id}_${mode}_${channel}*.txt
rm -f ber_${id}_${mode}_${channel}*.txt
rm -f per_${id}_${mode}_${channel}*.txt
for No in $aNo_list
do
No_adj=$((${No}-${snr_nudge}))
freedv_data_raw_tx --clip ${clip} --delay 1000 --txbpf ${clip} --bursts $Nbursts --testframes $Nbursts $mode /dev/zero - 2>${tx_log} | \
ch - - --No $No_adj ${ch_multipath} --fading_dir ${FADING_DIR} 2>>${ch_log} | \
freedv_data_raw_rx --testframes $mode - /dev/null 2>${rx_log} -v
SNRoffset=$(cat ${tx_log} | grep "mark:space" | tr -s ' ' | cut -d' ' -f 5)
SNRest=$(cat ${rx_log} | grep '\-BS\-' | tr -s ' ' | cut -d' ' -f17)
if [ ! -z "$SNRest" ]; then
echo ${SNRest} > snrest_${id}_${mode}_${channel}_${i}.txt
fi
BERmeas=$(cat ${rx_log} | grep 'BER......:' | cut -d' ' -f2)
PERmeas=$(cat ${rx_log} | grep 'Coded FER' | cut -d' ' -f3)
echo ${BERmeas} >> ber_${id}_${mode}_${channel}.txt
echo ${PERmeas} >> per_${id}_${mode}_${channel}.txt
i=$((i+1))
done
echo ${SNRoffset} > offset_${id}_${mode}_${channel}.txt
# trap not enough fading file samples (with mpp)
grep "Fading file finished" ${ch_log}
if [ $? -eq 0 ]; then
cat ${ch_log}
exit 1
fi
SNRch=$(cat ${ch_log} | grep SNR3k | tr -s ' ' | cut -d' ' -f3)
echo ${SNRch} > snr_${id}_${mode}_${channel}.txt
}
# Sanity check to make sure Octave/CML is set up OK
function test_ldpc {
echo "ldpcut; quit" | DISPLAY="" octave-cli -p ${CODEC2}/octave
if [ "$?" -ne 0 ]; then
echo "basic octave test failed, you may need to"
echo "(a) run ctests to create build_xxx/cml"
echo "(b) set up ~/.octaverc as per octave/ldpc.m"
exit 1
else
echo "OK"
fi
}
+27
View File
@@ -0,0 +1,27 @@
#!/usr/bin/env bash
#
# Tests reliable_text fading channel performance, using a simulated channel
results=$(mktemp -d)
mode=$1
snr=$2
min_text_packets=$3
clip=$4
build_folder=$5
fading_dir=${build_folder}/../unittest
rx=$build_folder/freedv_rx
tx=$build_folder/freedv_tx
if [ $clip -eq 1 ]; then
clip_args="--txbpf 1 --clip 1"
else
clip_args=
fi
$tx $mode ../raw/ve9qrp.raw - --reliabletext AB1CDEF $clip_args | $build_folder/ch - - --No $snr --mpp -f -5 --fading_dir $fading_dir > $results/reliable_fade.raw
$rx $mode $results/reliable_fade.raw /dev/null --txtrx $results/reliable_fade.txt --reliabletext
if [ `cat $results/reliable_fade.txt | wc -l` -ge $min_text_packets ]; then
exit 0
else
exit -1
fi
+79
View File
@@ -0,0 +1,79 @@
#!/usr/bin/env python3
""" sum_debug_alloc
Sum lines reported from codec2-dev/src/debug_alloc.h and report
Lines like:
CALLOC: run_ldpc_decoder(112, 32)
MALLOC: other_func(238)
"""
import fileinput
by_func = {}
by_addr = {}
def new_func():
rec = {}
rec['cur'] = 0
rec['max'] = 0
return(rec)
def new_addr():
rec = {}
rec['func'] = 0 # Where allocated, may not be same as where free'd!
rec['size'] = 0
return(rec)
for line in fileinput.input():
if ((line.startswith("MALLOC:")) or (line.startswith("CALLOC:"))):
if (line.startswith("MALLOC:")):
words = line.translate( str.maketrans("()", " ") ).strip().split()
func = words[1]
addr = words[2]
size = int(words[3])
elif (line.startswith("CALLOC:")):
words = line.translate( str.maketrans("(,)", " ") ).strip().split()
func = words[1]
addr = words[2]
size = int(words[3]) * int(words[4])
#print("Alloc: {} to {} in {}".format(size, addr, func))
if (not (func in by_func)): by_func[func] = new_func()
data = by_func[func]
data['cur'] += size
if (data['cur'] > data['max']):
data['max'] = data['cur']
if (addr in by_addr):
print("Error: duplicate allocation to {} in {}".format(addr, func))
else:
by_addr[addr] = new_addr()
by_addr[addr]['func'] = func
by_addr[addr]['size'] = size
elif (line.startswith("FREE:")):
words = line.translate( str.maketrans("(,)", " ") ).strip().split()
func = words[1]
addr = words[2]
if (addr in by_addr):
func_a = by_addr[addr]['func']
by_func[func_a]['cur'] -= by_addr[addr]['size']
del(by_addr[addr])
else:
print("Error: free of unallocated location {} in {}".format(addr, func))
#print("Free: {}:{} in {} to {}".format(func_a, addr, func, by_func[func_a]['cur']))
#####
total = 0
for func, sum in by_func.items():
max = by_func[func]['max']
print("{} = {}".format(func, max))
total += max
print("Total = {}".format(total))
+99
View File
@@ -0,0 +1,99 @@
/*
t16_8.c
David Rowe
May 10 2012
Unit test for 16 <-> 8 kHz sample rate conversion functions. I
Evaluated output by plotting using Octave and looking for jaggies:
pl("../unittest/out16.raw",1,3000)
pl("../unittest/out8.raw",1,3000)
Listening to it also shows up anything nasty:
$ play -s -2 -r 16000 out16.raw
$ play -s -2 -r 8000 out8.raw
*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "codec2_fdmdv.h"
#define N8 159 /* processing buffer size at 8 kHz (odd number deliberate) */
#define N16 (N8 * FDMDV_OS)
#define FRAMES 50
#define TWO_PI 6.283185307
#define FS 16000
#define SINE
int main() {
float in8k[FDMDV_OS_TAPS_8K + N8];
short in8k_short[N8];
float out16k[N16];
short out16k_short[N16];
FILE *f16;
float in16k[FDMDV_OS_TAPS_16K + N16];
float out8k[N16];
short out8k_short[N8];
FILE *f8, *f8in;
int i, f, t, t1;
float freq = 800.0;
f16 = fopen("out16.raw", "wb");
assert(f16 != NULL);
f8 = fopen("out8.raw", "wb");
assert(f8 != NULL);
f8in = fopen("in8.raw", "wb");
assert(f8in != NULL);
/* clear filter memories */
for (i = 0; i < FDMDV_OS_TAPS_8K; i++) in8k[i] = 0.0;
for (i = 0; i < FDMDV_OS_TAPS_16K; i++) in16k[i] = 0.0;
t = t1 = 0;
for (f = 0; f < FRAMES; f++) {
#ifdef DC
for (i = 0; i < N8; i++) in8k[FDMDV_OS_TAPS_8K + i] = 16000.0;
#endif
#ifdef SINE
for (i = 0; i < N8; i++, t++)
in8k[FDMDV_OS_TAPS_8K + i] =
16000.0 * cos(TWO_PI * t * freq / (FS / FDMDV_OS));
#endif
for (i = 0; i < N8; i++) in8k_short[i] = (short)in8k[i];
fwrite(in8k_short, sizeof(short), N8, f8in);
/* upsample */
fdmdv_8_to_16(out16k, &in8k[FDMDV_OS_TAPS_8K], N8);
/* save 16k to disk for plotting and check out */
for (i = 0; i < N16; i++) out16k_short[i] = (short)out16k[i];
fwrite(out16k_short, sizeof(short), N16, f16);
/* add a 6 kHz spurious signal, down sampler should
knock this out */
for (i = 0; i < N16; i++, t1++)
in16k[i + FDMDV_OS_TAPS_16K] =
out16k[i] + 16000.0 * cos(TWO_PI * t1 * 6000.0 / FS);
/* downsample */
fdmdv_16_to_8(out8k, &in16k[FDMDV_OS_TAPS_16K], N8);
/* save 8k to disk for plotting and check out */
for (i = 0; i < N8; i++) out8k_short[i] = (short)out8k[i];
fwrite(out8k_short, sizeof(short), N8, f8);
}
fclose(f16);
fclose(f8);
fclose(f8in);
return 0;
}
+90
View File
@@ -0,0 +1,90 @@
/*
t16_8_short.c
David Rowe
19 August 2014
Unit test for 16 <-> 8 kHz sample rate conversion functions. I
evaluated output by plotting using Octave and looking for jaggies:
pl("../unittest/out16_short.raw",1,3000)
pl("../unittest/out8.raw",1,3000)
Listening to it also shows up anything nasty:
$ play -s -2 -r 16000 out16_short.raw
$ play -s -2 -r 8000 out8.raw
*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "codec2_fdmdv.h"
#define N8 159 /* procssing buffer size at 8 kHz */
#define N16 (N8 * FDMDV_OS)
#define FRAMES 100
#define TWO_PI 6.283185307
#define FS 16000
#define SINE
int main() {
short in8k_short[FDMDV_OS_TAPS_8K + N8];
short out16k_short[N16];
FILE *f16;
short in16k_short[FDMDV_OS_TAPS_16K + N16];
short out8k_short[N16];
FILE *f8, *f8in;
int i, f, t, t1;
float freq = 800.0;
f16 = fopen("out16_short.raw", "wb");
assert(f16 != NULL);
f8 = fopen("out8_short.raw", "wb");
assert(f8 != NULL);
f8in = fopen("in8_short.raw", "wb");
assert(f8in != NULL);
/* clear filter memories */
for (i = 0; i < FDMDV_OS_TAPS_8K; i++) in8k_short[i] = 0;
for (i = 0; i < FDMDV_OS_TAPS_16K; i++) in16k_short[i] = 0;
t = t1 = 0;
for (f = 0; f < FRAMES; f++) {
#ifdef DC
for (i = 0; i < N8; i++) in8k_short[FDMDV_OS_TAPS_8K + i] = 16000.0;
#endif
#ifdef SINE
for (i = 0; i < N8; i++, t++)
in8k_short[FDMDV_OS_TAPS_8K + i] = 8000.0 * cos(TWO_PI * t * freq / FS);
#endif
fwrite(in8k_short, sizeof(short), N8, f8in);
/* upsample */
fdmdv_8_to_16_short(out16k_short, &in8k_short[FDMDV_OS_TAPS_8K], N8);
fwrite(out16k_short, sizeof(short), N16, f16);
/* add a 6 kHz spurious signal for fun, we want down sampler to
knock this out */
for (i = 0; i < N16; i++, t1++)
in16k_short[i + FDMDV_OS_TAPS_16K] =
out16k_short[i] + 8000.0 * cos(TWO_PI * t1 * 6000.0 / FS);
/* downsample */
fdmdv_16_to_8_short(out8k_short, &in16k_short[FDMDV_OS_TAPS_16K], N8);
/* save 8k to disk for plotting and check out */
fwrite(out8k_short, sizeof(short), N8, f8);
}
fclose(f16);
fclose(f8);
fclose(f8in);
return 0;
}
+104
View File
@@ -0,0 +1,104 @@
/*
t48_8.c
David Rowe
May 10 2012
Unit test for 48 to 8 kHz sample rate conversion functions. I
evaluated output by plotting using Octave and looking for jaggies:
pl("../unittest/out48.raw",1,3000)
pl("../unittest/out8.raw",1,3000)
Listening to it also shows up anything nasty:
$ play -s -2 -r 48000 out48.raw
$ play -s -2 -r 8000 out8.raw
*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "codec2_fdmdv.h"
#define N8 180 /* processing buffer size at 8 kHz */
#define N48 (N8 * FDMDV_OS_48)
#define MEM8 FDMDV_OS_TAPS_48_8K
#define FRAMES 50
#define TWO_PI 6.283185307
#define FS 48000
#define SINE
int main() {
float in8k[MEM8 + N8];
short in8k_short[N8];
float out48k[N48];
short out48k_short[N48];
FILE *f48;
float in48k[FDMDV_OS_TAPS_48K + N48];
float out8k[N48];
short out8k_short[N8];
FILE *f8, *f8in;
int i, f, t, t1;
float freq = 800.0;
f48 = fopen("out48.raw", "wb");
assert(f48 != NULL);
f8 = fopen("out8.raw", "wb");
assert(f8 != NULL);
f8in = fopen("in8.raw", "wb");
assert(f8in != NULL);
/* clear filter memories */
for (i = 0; i < MEM8; i++) in8k[i] = 0.0;
for (i = 0; i < FDMDV_OS_TAPS_48K; i++) in48k[i] = 0.0;
t = t1 = 0;
for (f = 0; f < FRAMES; f++) {
#ifdef DC
for (i = 0; i < N8; i++) in8k[MEM8 + i] = 16000.0;
#endif
#ifdef SINE
for (i = 0; i < N8; i++, t++)
in8k[MEM8 + i] = 16000.0 * cos(TWO_PI * t * freq / (FS / FDMDV_OS_48));
#endif
for (i = 0; i < N8; i++) in8k_short[i] = (short)in8k[i];
fwrite(in8k_short, sizeof(short), N8, f8in);
/* upsample */
fdmdv_8_to_48(out48k, &in8k[MEM8], N8);
/* save 48k to disk for plotting and check out */
for (i = 0; i < N48; i++) out48k_short[i] = (short)out48k[i];
fwrite(out48k_short, sizeof(short), N48, f48);
/* add a 10 kHz spurious signal for fun, we want down sampler to
knock this out */
for (i = 0; i < N48; i++, t1++)
in48k[i + FDMDV_OS_TAPS_48K] =
out48k[i] + 16000.0 * cos(TWO_PI * t1 * 1E4 / FS);
/* downsample */
fdmdv_48_to_8(out8k, &in48k[FDMDV_OS_TAPS_48K], N8);
/* save 8k to disk for plotting and check out */
for (i = 0; i < N8; i++) out8k_short[i] = (short)out8k[i];
fwrite(out8k_short, sizeof(short), N8, f8);
}
fclose(f48);
fclose(f8);
fclose(f8in);
return 0;
}
+81
View File
@@ -0,0 +1,81 @@
/*
t48_8_short.c
David Rowe
Dec 2021
*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "codec2_fdmdv.h"
#define N8 180 /* processing buffer size at 8 kHz */
#define N48 (N8 * FDMDV_OS_48)
#define MEM8 FDMDV_OS_TAPS_48_8K
#define FRAMES 50
#define TWO_PI 6.283185307
#define FS 48000
#define SINE
int main() {
short in8k[MEM8 + N8];
short out48k[N48];
FILE *f48;
short in48k[FDMDV_OS_TAPS_48K + N48];
short out8k[N48];
FILE *f8, *f8in;
int i, f, t, t1;
float freq = 800.0;
f48 = fopen("out48.raw", "wb");
assert(f48 != NULL);
f8 = fopen("out8.raw", "wb");
assert(f8 != NULL);
f8in = fopen("in8.raw", "wb");
assert(f8in != NULL);
/* clear filter memories */
for (i = 0; i < MEM8; i++) in8k[i] = 0.0;
for (i = 0; i < FDMDV_OS_TAPS_48K; i++) in48k[i] = 0.0;
t = t1 = 0;
for (f = 0; f < FRAMES; f++) {
#ifdef DC
for (i = 0; i < N8; i++) in8k[MEM8 + i] = 16000.0;
#endif
#ifdef SINE
for (i = 0; i < N8; i++, t++)
in8k[MEM8 + i] = 16000.0 * cos(TWO_PI * t * freq / (FS / FDMDV_OS_48));
#endif
fwrite(&in8k[MEM8], sizeof(short), N8, f8in);
/* upsample */
fdmdv_8_to_48_short(out48k, &in8k[MEM8], N8);
/* save 48k to disk for plotting and check out */
fwrite(out48k, sizeof(short), N48, f48);
/* add a 10 kHz spurious signal for fun, we want down sampler to
knock this out */
for (i = 0; i < N48; i++, t1++)
in48k[i + FDMDV_OS_TAPS_48K] =
out48k[i] + 16000.0 * cos(TWO_PI * t1 * 1E4 / FS);
/* downsample */
fdmdv_48_to_8_short(out8k, &in48k[FDMDV_OS_TAPS_48K], N8);
/* save 8k to disk for plotting and check out */
fwrite(out8k, sizeof(short), N8, f8);
}
fclose(f48);
fclose(f8);
fclose(f8in);
return 0;
}
+324
View File
@@ -0,0 +1,324 @@
/*---------------------------------------------------------------------------*\
FILE........: tcohpsk.c
AUTHOR......: David Rowe
DATE CREATED: March 2015
Tests for the C version of the coherent PSK FDM modem. This program
outputs a file of Octave vectors that are loaded and automatically
tested against the Octave version of the modem by the Octave script
tcohpsk.m
\*---------------------------------------------------------------------------*/
/*
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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "codec2_cohpsk.h"
#include "codec2_fdmdv.h"
#include "cohpsk_defs.h"
#include "cohpsk_internal.h"
#include "comp_prim.h"
#include "fdmdv_internal.h"
#include "noise_samples.h"
#include "octave.h"
#define FRAMES \
30 /* LOG_FRAMES is #defined in cohpsk_internal.h */
#define SYNC_FRAMES \
12 /* sync state uses up extra log storage as we reprocess several times */
#define FRAMESL \
(SYNC_FRAMES * FRAMES) /* worst case is every frame is out of sync */
#define FOFF 58.7
#define DFOFF (-0.5 / (float)COHPSK_FS)
#define ESNODB 8
#define PPM -1500
extern float pilots_coh[][PILOTS_NC];
int main(int argc, char *argv[]) {
struct COHPSK *coh;
int tx_bits[COHPSK_BITS_PER_FRAME];
COMP tx_symb[NSYMROWPILOT][COHPSK_NC * COHPSK_ND];
COMP tx_fdm_frame[COHPSK_M * NSYMROWPILOT];
COMP ch_fdm_frame[COHPSK_M * NSYMROWPILOT];
// COMP rx_fdm_frame_bb[M*NSYMROWPILOT];
// COMP ch_symb[NSYMROWPILOT][COHPSK_NC*COHPSK_ND];
float rx_bits_sd[COHPSK_BITS_PER_FRAME];
int rx_bits[COHPSK_BITS_PER_FRAME];
int tx_bits_log[COHPSK_BITS_PER_FRAME * FRAMES];
COMP tx_symb_log[NSYMROWPILOT * FRAMES][COHPSK_NC * COHPSK_ND];
COMP tx_fdm_frame_log[COHPSK_M * NSYMROWPILOT * FRAMES];
COMP ch_fdm_frame_log[COHPSK_M * NSYMROWPILOT * FRAMES];
COMP ch_fdm_frame_log_out[(COHPSK_M * NSYMROWPILOT + 1) * FRAMES];
// COMP rx_fdm_frame_bb_log[M*NSYMROWPILOT*FRAMES];
// COMP ch_symb_log[NSYMROWPILOT*FRAMES][COHPSK_NC*COHPSK_ND];
COMP ct_symb_ff_log[NSYMROWPILOT * FRAMES][COHPSK_NC * COHPSK_ND];
float rx_amp_log[NSYMROW * FRAMES][COHPSK_NC * COHPSK_ND];
float rx_phi_log[NSYMROW * FRAMES][COHPSK_NC * COHPSK_ND];
COMP rx_symb_log[NSYMROW * FRAMES][COHPSK_NC * COHPSK_ND];
int rx_bits_log[COHPSK_BITS_PER_FRAME * FRAMES];
FILE *fout;
int f, r, c, log_r, log_data_r, noise_r, ff_log_r, i;
double foff;
COMP foff_rect, phase_ch;
struct FDMDV *fdmdv;
// COMP rx_filt[COHPSK_NC*COHPSK_ND][P+1];
// int rx_filt_log_col_index = 0;
// float env[NT*P];
// float __attribute__((unused)) rx_timing;
COMP tx_onesym[COHPSK_NC * COHPSK_ND];
// COMP rx_onesym[COHPSK_NC*COHPSK_ND];
// int rx_baseband_log_col_index = 0;
// COMP rx_baseband_log[COHPSK_NC*COHPSK_ND][(M+M/P)*NSYMROWPILOT*FRAMES];
float f_est_log[FRAMES], sig_rms_log[FRAMES], noise_rms_log[FRAMES];
int f_est_samples;
int log_bits;
float EsNo, variance;
COMP scaled_noise;
int reliable_sync_bit;
int ch_fdm_frame_log_index, nin_frame, tmp, nout;
coh = cohpsk_create();
fdmdv = coh->fdmdv;
assert(coh != NULL);
cohpsk_set_verbose(coh, 1);
/* these puppies are used for logging data in the bowels on the modem */
coh->rx_baseband_log_col_sz =
(COHPSK_M + COHPSK_M / P) * NSYMROWPILOT * FRAMESL;
coh->rx_baseband_log = (COMP *)malloc(sizeof(COMP) * COHPSK_NC * COHPSK_ND *
coh->rx_baseband_log_col_sz);
coh->rx_filt_log_col_sz = (P + 1) * NSYMROWPILOT * FRAMESL;
coh->rx_filt_log = (COMP *)malloc(sizeof(COMP) * COHPSK_NC * COHPSK_ND *
coh->rx_filt_log_col_sz);
coh->ch_symb_log_col_sz = COHPSK_NC * COHPSK_ND;
coh->ch_symb_log = (COMP *)malloc(sizeof(COMP) * NSYMROWPILOT * FRAMESL *
coh->ch_symb_log_col_sz);
coh->rx_timing_log = (float *)malloc(sizeof(float) * NSYMROWPILOT * FRAMESL);
/* init stuff */
log_r = log_data_r = noise_r = log_bits = ff_log_r = f_est_samples = 0;
phase_ch.real = 1.0;
phase_ch.imag = 0.0;
foff = FOFF;
/* each carrier has power = 2, total power 2Nc, total symbol rate
NcRs, noise BW B=Fs Es/No = (C/Rs)/(N/B), N = var =
2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No) */
EsNo = pow(10.0, ESNODB / 10.0);
variance = 2.0 * COHPSK_FS / (COHPSK_RS * EsNo);
// fprintf(stderr, "doff: %e\n", DFOFF);
/* Main Loop
* ---------------------------------------------------------------------*/
for (f = 0; f < FRAMES; f++) {
/* --------------------------------------------------------*\
Mod
\*---------------------------------------------------------*/
cohpsk_get_test_bits(coh, tx_bits);
bits_to_qpsk_symbols(tx_symb, (int *)tx_bits, COHPSK_BITS_PER_FRAME);
for (r = 0; r < NSYMROWPILOT; r++) {
for (c = 0; c < COHPSK_NC * COHPSK_ND; c++) tx_onesym[c] = tx_symb[r][c];
tx_filter_and_upconvert_coh(
&tx_fdm_frame[r * COHPSK_M], COHPSK_NC * COHPSK_ND, tx_onesym,
fdmdv->tx_filter_memory, fdmdv->phase_tx, fdmdv->freq,
&fdmdv->fbb_phase_tx, fdmdv->fbb_rect);
}
cohpsk_clip(tx_fdm_frame, COHPSK_CLIP, NSYMROWPILOT * COHPSK_M);
/* --------------------------------------------------------*\
Channel
\*---------------------------------------------------------*/
for (r = 0; r < NSYMROWPILOT * COHPSK_M; r++) {
foff_rect.real = cos(2.0 * M_PI * foff / COHPSK_FS);
foff_rect.imag = sin(2.0 * M_PI * foff / COHPSK_FS);
foff += DFOFF;
phase_ch = cmult(phase_ch, foff_rect);
ch_fdm_frame[r] = cmult(tx_fdm_frame[r], phase_ch);
}
phase_ch.real /= cabsolute(phase_ch);
phase_ch.imag /= cabsolute(phase_ch);
// fprintf(stderr, "%f\n", foff);
for (r = 0; r < NSYMROWPILOT * COHPSK_M; r++, noise_r++) {
scaled_noise = fcmult(sqrt(variance), noise[noise_r]);
ch_fdm_frame[r] = cadd(ch_fdm_frame[r], scaled_noise);
}
/* optional gain to test level sensitivity */
for (r = 0; r < NSYMROWPILOT * COHPSK_M; r++) {
ch_fdm_frame[r] = fcmult(1.0, ch_fdm_frame[r]);
}
/* tx vector logging */
memcpy(&tx_bits_log[COHPSK_BITS_PER_FRAME * f], tx_bits,
sizeof(int) * COHPSK_BITS_PER_FRAME);
memcpy(&tx_fdm_frame_log[COHPSK_M * NSYMROWPILOT * f], tx_fdm_frame,
sizeof(COMP) * COHPSK_M * NSYMROWPILOT);
memcpy(&ch_fdm_frame_log[COHPSK_M * NSYMROWPILOT * f], ch_fdm_frame,
sizeof(COMP) * COHPSK_M * NSYMROWPILOT);
for (r = 0; r < NSYMROWPILOT; r++, log_r++) {
for (c = 0; c < COHPSK_NC * COHPSK_ND; c++) {
tx_symb_log[log_r][c] = tx_symb[r][c];
}
}
}
/* Fs offset simulation */
nout = cohpsk_fs_offset(ch_fdm_frame_log_out, ch_fdm_frame_log,
COHPSK_M * NSYMROWPILOT * FRAMES, PPM);
assert(nout < (COHPSK_M * NSYMROWPILOT + 1) * FRAMES);
nin_frame = COHPSK_NOM_SAMPLES_PER_FRAME;
ch_fdm_frame_log_index = 0;
/* --------------------------------------------------------*\
Demod
\*---------------------------------------------------------*/
for (f = 0; f < FRAMES; f++) {
coh->frame = f;
// printf("nin_frame: %d\n", nin_frame);
assert(ch_fdm_frame_log_index < COHPSK_M * NSYMROWPILOT * FRAMES);
tmp = nin_frame;
cohpsk_demod(coh, rx_bits_sd, &reliable_sync_bit,
&ch_fdm_frame_log_out[ch_fdm_frame_log_index], &nin_frame);
for (i = 0; i < COHPSK_BITS_PER_FRAME; i++)
rx_bits[i] = rx_bits_sd[i] < 0.0;
ch_fdm_frame_log_index += tmp;
/* --------------------------------------------------------*\
Log each vector
\*---------------------------------------------------------*/
if (coh->sync == 1) {
for (r = 0; r < NSYMROWPILOT; r++, ff_log_r++) {
for (c = 0; c < COHPSK_NC * COHPSK_ND; c++) {
ct_symb_ff_log[ff_log_r][c] = coh->ct_symb_ff_buf[r][c];
}
}
for (r = 0; r < NSYMROW; r++, log_data_r++) {
for (c = 0; c < COHPSK_NC * COHPSK_ND; c++) {
rx_amp_log[log_data_r][c] = coh->amp_[r][c];
rx_phi_log[log_data_r][c] = coh->phi_[r][c];
rx_symb_log[log_data_r][c] = coh->rx_symb[r][c];
}
}
memcpy(&rx_bits_log[COHPSK_BITS_PER_FRAME * log_bits], rx_bits,
sizeof(int) * COHPSK_BITS_PER_FRAME);
log_bits++;
f_est_log[f_est_samples] = coh->f_est;
sig_rms_log[f_est_samples] = coh->sig_rms;
noise_rms_log[f_est_samples] = coh->noise_rms;
f_est_samples++;
;
}
assert(log_r <= NSYMROWPILOT * FRAMES);
assert(noise_r <= NSYMROWPILOT * COHPSK_M * FRAMES);
assert(log_data_r <= NSYMROW * FRAMES);
printf("\r [%d]", f + 1);
}
printf("\n");
/*---------------------------------------------------------*\
Dump logs to Octave file for evaluation
by tcohpsk.m Octave script
\*---------------------------------------------------------*/
fout = fopen("tcohpsk_out.txt", "wt");
assert(fout != NULL);
fprintf(fout, "# Created by tcohpsk.c\n");
octave_save_int(fout, "tx_bits_log_c", tx_bits_log, 1,
COHPSK_BITS_PER_FRAME * FRAMES);
octave_save_complex(fout, "tx_symb_log_c", (COMP *)tx_symb_log,
NSYMROWPILOT * FRAMES, COHPSK_NC * COHPSK_ND,
COHPSK_NC * COHPSK_ND);
octave_save_complex(fout, "tx_fdm_frame_log_c", (COMP *)tx_fdm_frame_log, 1,
COHPSK_M * NSYMROWPILOT * FRAMES,
COHPSK_M * NSYMROWPILOT * FRAMES);
octave_save_complex(fout, "ch_fdm_frame_log_c", (COMP *)ch_fdm_frame_log_out,
1, nout - 1, nout - 1);
// octave_save_complex(fout, "rx_fdm_frame_bb_log_c",
// (COMP*)rx_fdm_frame_bb_log, 1, M*NSYMROWPILOT*FRAMES,
// M*NSYMROWPILOT*FRAMES);
octave_save_complex(fout, "rx_baseband_log_c", (COMP *)coh->rx_baseband_log,
COHPSK_NC * COHPSK_ND, coh->rx_baseband_log_col_index,
coh->rx_baseband_log_col_sz);
octave_save_complex(fout, "rx_filt_log_c", (COMP *)coh->rx_filt_log,
COHPSK_NC * COHPSK_ND, coh->rx_filt_log_col_index,
coh->rx_filt_log_col_sz);
octave_save_complex(fout, "ch_symb_log_c", (COMP *)coh->ch_symb_log,
coh->ch_symb_log_r, COHPSK_NC * COHPSK_ND,
COHPSK_NC * COHPSK_ND);
octave_save_float(fout, "rx_timing_log_c", (float *)coh->rx_timing_log, 1,
coh->rx_timing_log_index, coh->rx_timing_log_index);
octave_save_complex(fout, "ct_symb_ff_log_c", (COMP *)ct_symb_ff_log,
NSYMROWPILOT * FRAMES, COHPSK_NC * COHPSK_ND,
COHPSK_NC * COHPSK_ND);
octave_save_float(fout, "rx_amp_log_c", (float *)rx_amp_log, log_data_r,
COHPSK_NC * COHPSK_ND, COHPSK_NC * COHPSK_ND);
octave_save_float(fout, "rx_phi_log_c", (float *)rx_phi_log, log_data_r,
COHPSK_NC * COHPSK_ND, COHPSK_NC * COHPSK_ND);
octave_save_complex(fout, "rx_symb_log_c", (COMP *)rx_symb_log, log_data_r,
COHPSK_NC * COHPSK_ND, COHPSK_NC * COHPSK_ND);
octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1,
COHPSK_BITS_PER_FRAME * log_bits);
octave_save_float(fout, "f_est_log_c", &f_est_log[1], 1, f_est_samples - 1,
f_est_samples - 1);
octave_save_float(fout, "sig_rms_log_c", sig_rms_log, 1, f_est_samples,
f_est_samples - 1);
octave_save_float(fout, "noise_rms_log_c", noise_rms_log, 1, f_est_samples,
f_est_samples);
#ifdef XX
#endif
fclose(fout);
cohpsk_destroy(coh);
return 0;
}
+32
View File
@@ -0,0 +1,32 @@
/*---------------------------------------------------------------------------*\
FILE........: tesno_est.c
AUTHORS.....: David Rowe
DATE CREATED: Mar 2021
Test for C port of Es/No estimator.
\*---------------------------------------------------------------------------*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "ofdm_internal.h"
int main(int argc, char *argv[]) {
FILE *fin = fopen(argv[1], "rb");
assert(fin != NULL);
size_t nsym = atoi(argv[2]);
assert(nsym >= 0);
complex float rx_sym[nsym];
size_t nread = fread(rx_sym, sizeof(complex float), nsym, fin);
assert(nread == nsym);
fclose(fin);
float EsNodB = ofdm_esno_est_calc(rx_sym, nsym);
printf("%f\n", EsNodB);
return 0;
}
+12
View File
@@ -0,0 +1,12 @@
#!/usr/bin/env bash
# test_700c_eq.sh
# make sure 700C EQ is reducing VQ distortion
results=$(mktemp)
c2enc 700C ../raw/kristoff.raw /dev/null --var 2> $results
var=$(cat $results | sed -n "s/.*var: \([0-9..]*\) .*/\1/p")
c2enc 700C ../raw/kristoff.raw /dev/null --var --eq 2> $results
var_eq=$(cat $results | sed -n "s/.*var: \([0-9..]*\) .*/\1/p")
printf "var: %5.2f var_eq: %5.2f\n" $var $var_eq
python3 -c "import sys; sys.exit(0) if $var_eq<=$var else sys.exit(1)"
+319
View File
@@ -0,0 +1,319 @@
/*---------------------------------------------------------------------------*\
FILE........: tfdmdv.c
AUTHOR......: David Rowe
DATE CREATED: April 16 2012
Tests for the C version of the FDMDV modem. This program outputs a
file of Octave vectors that are loaded and automatically tested
against the Octave version of the modem by the Octave script
tfmddv.m
\*---------------------------------------------------------------------------*/
/*
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 <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "codec2_fdmdv.h"
#include "fdmdv_internal.h"
#include "octave.h"
#define FRAMES 35
#define CHANNEL_BUF_SIZE (10 * M_FAC)
extern float pilot_coeff[];
int main(int argc, char *argv[]) {
struct FDMDV *fdmdv;
int tx_bits[FDMDV_BITS_PER_FRAME];
COMP tx_symbols[FDMDV_NC + 1];
COMP tx_fdm[M_FAC];
float channel[CHANNEL_BUF_SIZE];
int channel_count;
COMP rx_fdm[M_FAC + M_FAC / P];
float foff_coarse;
int nin, next_nin;
COMP rx_fdm_fcorr[M_FAC + M_FAC / P];
COMP rx_fdm_filter[M_FAC + M_FAC / P];
COMP rx_filt[NC + 1][P + 1];
float rx_timing;
float env[NT * P];
COMP rx_symbols[FDMDV_NC + 1];
int rx_bits[FDMDV_BITS_PER_FRAME];
float foff_fine;
int sync_bit, reliable_sync_bit;
int tx_bits_log[FDMDV_BITS_PER_FRAME * FRAMES];
COMP tx_symbols_log[(FDMDV_NC + 1) * FRAMES];
COMP tx_fdm_log[M_FAC * FRAMES];
COMP pilot_baseband1_log[NPILOTBASEBAND * FRAMES];
COMP pilot_baseband2_log[NPILOTBASEBAND * FRAMES];
COMP pilot_lpf1_log[NPILOTLPF * FRAMES];
COMP pilot_lpf2_log[NPILOTLPF * FRAMES];
COMP S1_log[MPILOTFFT * FRAMES];
COMP S2_log[MPILOTFFT * FRAMES];
float foff_coarse_log[FRAMES];
float foff_log[FRAMES];
COMP rx_fdm_filter_log[(M_FAC + M_FAC / P) * FRAMES];
int rx_fdm_filter_log_index;
COMP rx_filt_log[NC + 1][(P + 1) * FRAMES];
int rx_filt_log_col_index;
float env_log[NT * P * FRAMES];
float rx_timing_log[FRAMES];
COMP rx_symbols_log[FDMDV_NC + 1][FRAMES];
COMP phase_difference_log[FDMDV_NC + 1][FRAMES];
float sig_est_log[FDMDV_NC + 1][FRAMES];
float noise_est_log[FDMDV_NC + 1][FRAMES];
int rx_bits_log[FDMDV_BITS_PER_FRAME * FRAMES];
float foff_fine_log[FRAMES];
int sync_bit_log[FRAMES];
int sync_log[FRAMES];
int nin_log[FRAMES];
FILE *fout;
int f, c, i, j;
fdmdv = fdmdv_create(FDMDV_NC);
next_nin = M_FAC;
channel_count = 0;
rx_fdm_filter_log_index = 0;
rx_filt_log_col_index = 0;
printf("sizeof FDMDV states: %zd bytes\n", sizeof(struct FDMDV));
for (f = 0; f < FRAMES; f++) {
/* --------------------------------------------------------*\
Modulator
\*---------------------------------------------------------*/
fdmdv_get_test_bits(fdmdv, tx_bits);
bits_to_dqpsk_symbols(tx_symbols, FDMDV_NC, fdmdv->prev_tx_symbols, tx_bits,
&fdmdv->tx_pilot_bit, 0);
memcpy(fdmdv->prev_tx_symbols, tx_symbols, sizeof(COMP) * (FDMDV_NC + 1));
tx_filter_and_upconvert(tx_fdm, FDMDV_NC, tx_symbols,
fdmdv->tx_filter_memory, fdmdv->phase_tx,
fdmdv->freq, &fdmdv->fbb_phase_tx, fdmdv->fbb_rect);
/* --------------------------------------------------------*\
Channel
\*---------------------------------------------------------*/
nin = next_nin;
// nin = M_FAC; // when debugging good idea to uncomment this to "open
// loop"
/* add M_FAC tx samples to end of buffer */
assert((channel_count + M_FAC) < CHANNEL_BUF_SIZE);
for (i = 0; i < M_FAC; i++) channel[channel_count + i] = tx_fdm[i].real;
channel_count += M_FAC;
/* take nin samples from start of buffer */
for (i = 0; i < nin; i++) {
rx_fdm[i].real = channel[i];
rx_fdm[i].imag = 0;
}
/* shift buffer back */
for (i = 0, j = nin; j < channel_count; i++, j++) channel[i] = channel[j];
channel_count -= nin;
/* --------------------------------------------------------*\
Demodulator
\*---------------------------------------------------------*/
/* shift down to complex baseband */
fdmdv_freq_shift(rx_fdm, rx_fdm, -FDMDV_FCENTRE, &fdmdv->fbb_phase_rx, nin);
/* freq offset estimation and correction */
// fdmdv->sync = 0; // when debugging good idea to uncomment this to "open
// loop"
foff_coarse = rx_est_freq_offset(fdmdv, rx_fdm, nin, !fdmdv->sync);
if (fdmdv->sync == 0) fdmdv->foff = foff_coarse;
fdmdv_freq_shift(rx_fdm_fcorr, rx_fdm, -fdmdv->foff,
&fdmdv->foff_phase_rect, nin);
/* baseband processing */
rxdec_filter(rx_fdm_filter, rx_fdm_fcorr, fdmdv->rxdec_lpf_mem, nin);
down_convert_and_rx_filter(rx_filt, fdmdv->Nc, rx_fdm_filter,
fdmdv->rx_fdm_mem, fdmdv->phase_rx, fdmdv->freq,
fdmdv->freq_pol, nin, M_FAC / Q);
rx_timing = rx_est_timing(rx_symbols, FDMDV_NC, rx_filt,
fdmdv->rx_filter_mem_timing, env, nin, M_FAC);
foff_fine =
qpsk_to_bits(rx_bits, &sync_bit, FDMDV_NC, fdmdv->phase_difference,
fdmdv->prev_rx_symbols, rx_symbols, 0);
// for(i=0; i<FDMDV_NC;i++)
// printf("rx_symbols: %f %f prev_rx_symbols: %f %f phase_difference: %f
// %f\n", rx_symbols[i].real, rx_symbols[i].imag,
// fdmdv->prev_rx_symbols[i].real, fdmdv->prev_rx_symbols[i].imag,
// fdmdv->phase_difference[i].real,
// fdmdv->phase_difference[i].imag);
// if (f==1)
// exit(0);
snr_update(fdmdv->sig_est, fdmdv->noise_est, FDMDV_NC,
fdmdv->phase_difference);
memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP) * (FDMDV_NC + 1));
next_nin = M_FAC;
if (rx_timing > 2 * M_FAC / P) next_nin += M_FAC / P;
if (rx_timing < 0) next_nin -= M_FAC / P;
fdmdv->sync = freq_state(&reliable_sync_bit, sync_bit, &fdmdv->fest_state,
&fdmdv->timer, fdmdv->sync_mem);
fdmdv->foff -= TRACK_COEFF * foff_fine;
/* --------------------------------------------------------*\
Log each vector
\*---------------------------------------------------------*/
memcpy(&tx_bits_log[FDMDV_BITS_PER_FRAME * f], tx_bits,
sizeof(int) * FDMDV_BITS_PER_FRAME);
memcpy(&tx_symbols_log[(FDMDV_NC + 1) * f], tx_symbols,
sizeof(COMP) * (FDMDV_NC + 1));
memcpy(&tx_fdm_log[M_FAC * f], tx_fdm, sizeof(COMP) * M_FAC);
memcpy(&pilot_baseband1_log[f * NPILOTBASEBAND], fdmdv->pilot_baseband1,
sizeof(COMP) * NPILOTBASEBAND);
memcpy(&pilot_baseband2_log[f * NPILOTBASEBAND], fdmdv->pilot_baseband2,
sizeof(COMP) * NPILOTBASEBAND);
memcpy(&pilot_lpf1_log[f * NPILOTLPF], fdmdv->pilot_lpf1,
sizeof(COMP) * NPILOTLPF);
memcpy(&pilot_lpf2_log[f * NPILOTLPF], fdmdv->pilot_lpf2,
sizeof(COMP) * NPILOTLPF);
memcpy(&S1_log[f * MPILOTFFT], fdmdv->S1, sizeof(COMP) * MPILOTFFT);
memcpy(&S2_log[f * MPILOTFFT], fdmdv->S2, sizeof(COMP) * MPILOTFFT);
foff_coarse_log[f] = foff_coarse;
foff_log[f] = fdmdv->foff;
/* rx filtering */
for (i = 0; i < nin; i++)
rx_fdm_filter_log[rx_fdm_filter_log_index + i] = rx_fdm_filter[i];
rx_fdm_filter_log_index += nin;
for (c = 0; c < NC + 1; c++) {
for (i = 0; i < (P * nin) / M_FAC; i++)
rx_filt_log[c][rx_filt_log_col_index + i] = rx_filt[c][i];
}
rx_filt_log_col_index += (P * nin) / M_FAC;
/* timing estimation */
memcpy(&env_log[NT * P * f], env, sizeof(float) * NT * P);
rx_timing_log[f] = rx_timing;
nin_log[f] = nin;
for (c = 0; c < FDMDV_NC + 1; c++) {
rx_symbols_log[c][f] = rx_symbols[c];
phase_difference_log[c][f] = fdmdv->phase_difference[c];
}
/* qpsk_to_bits() */
memcpy(&rx_bits_log[FDMDV_BITS_PER_FRAME * f], rx_bits,
sizeof(int) * FDMDV_BITS_PER_FRAME);
for (c = 0; c < FDMDV_NC + 1; c++) {
sig_est_log[c][f] = fdmdv->sig_est[c];
noise_est_log[c][f] = fdmdv->noise_est[c];
}
foff_fine_log[f] = foff_fine;
sync_bit_log[f] = sync_bit;
sync_log[f] = fdmdv->sync;
}
/*---------------------------------------------------------*\
Dump logs to Octave file for evaluation
by tfdmdv.m Octave script
\*---------------------------------------------------------*/
fout = fopen("tfdmdv_out.txt", "wt");
assert(fout != NULL);
fprintf(fout, "# Created by tfdmdv.c\n");
octave_save_int(fout, "tx_bits_log_c", tx_bits_log, 1,
FDMDV_BITS_PER_FRAME * FRAMES);
octave_save_complex(fout, "tx_symbols_log_c", tx_symbols_log, 1,
(FDMDV_NC + 1) * FRAMES, (FDMDV_NC + 1) * FRAMES);
octave_save_complex(fout, "tx_fdm_log_c", (COMP *)tx_fdm_log, 1,
M_FAC * FRAMES, M_FAC * FRAMES);
octave_save_complex(fout, "pilot_lut_c", (COMP *)fdmdv->pilot_lut, 1,
NPILOT_LUT, NPILOT_LUT);
octave_save_complex(fout, "pilot_baseband1_log_c", pilot_baseband1_log, 1,
NPILOTBASEBAND * FRAMES, NPILOTBASEBAND * FRAMES);
octave_save_complex(fout, "pilot_baseband2_log_c", pilot_baseband2_log, 1,
NPILOTBASEBAND * FRAMES, NPILOTBASEBAND * FRAMES);
octave_save_float(fout, "pilot_coeff_c", pilot_coeff, 1, NPILOTCOEFF,
NPILOTCOEFF);
octave_save_complex(fout, "pilot_lpf1_log_c", pilot_lpf1_log, 1,
NPILOTLPF * FRAMES, NPILOTLPF * FRAMES);
octave_save_complex(fout, "pilot_lpf2_log_c", pilot_lpf2_log, 1,
NPILOTLPF * FRAMES, NPILOTLPF * FRAMES);
octave_save_complex(fout, "S1_log_c", S1_log, 1, MPILOTFFT * FRAMES,
MPILOTFFT * FRAMES);
octave_save_complex(fout, "S2_log_c", S2_log, 1, MPILOTFFT * FRAMES,
MPILOTFFT * FRAMES);
octave_save_float(fout, "foff_log_c", foff_log, 1, FRAMES, FRAMES);
octave_save_float(fout, "foff_coarse_log_c", foff_coarse_log, 1, FRAMES,
FRAMES);
octave_save_complex(fout, "rx_fdm_filter_log_c", (COMP *)rx_fdm_filter_log, 1,
rx_fdm_filter_log_index, rx_fdm_filter_log_index);
octave_save_complex(fout, "rx_filt_log_c", (COMP *)rx_filt_log,
(FDMDV_NC + 1), rx_filt_log_col_index, (P + 1) * FRAMES);
octave_save_float(fout, "env_log_c", env_log, 1, NT * P * FRAMES,
NT * P * FRAMES);
octave_save_float(fout, "rx_timing_log_c", rx_timing_log, 1, FRAMES, FRAMES);
octave_save_complex(fout, "rx_symbols_log_c", (COMP *)rx_symbols_log,
(FDMDV_NC + 1), FRAMES, FRAMES);
octave_save_complex(fout, "phase_difference_log_c",
(COMP *)phase_difference_log, (FDMDV_NC + 1), FRAMES,
FRAMES);
octave_save_float(fout, "sig_est_log_c", (float *)sig_est_log, (FDMDV_NC + 1),
FRAMES, FRAMES);
octave_save_float(fout, "noise_est_log_c", (float *)noise_est_log,
(FDMDV_NC + 1), FRAMES, FRAMES);
octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1,
FDMDV_BITS_PER_FRAME * FRAMES);
octave_save_float(fout, "foff_fine_log_c", foff_fine_log, 1, FRAMES, FRAMES);
octave_save_int(fout, "sync_bit_log_c", sync_bit_log, 1, FRAMES);
octave_save_int(fout, "sync_log_c", sync_log, 1, FRAMES);
octave_save_int(fout, "nin_log_c", nin_log, 1, FRAMES);
fclose(fout);
fdmdv_destroy(fdmdv);
return 0;
}
+103
View File
@@ -0,0 +1,103 @@
/*
tfifo.c
David Rowe
Nov 19 2012
Tests FIFOs, in particular thread safety.
*/
#include <assert.h>
#include <pthread.h>
#include <stdio.h>
#include "codec2_fifo.h"
#define FIFO_SZ 1024
#define WRITE_SZ 10
#define READ_SZ 8
#define N_MAX 100
#define LOOPS 1000000
int run_thread = 1;
struct FIFO *f;
void writer(void);
void *writer_thread(void *data);
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
#define USE_THREADS
//#define USE_MUTEX
int main() {
pthread_t awriter_thread;
int i, j;
short read_buf[READ_SZ];
int n_out = 0;
int success;
f = codec2_fifo_create(FIFO_SZ);
#ifdef USE_THREADS
pthread_create(&awriter_thread, NULL, writer_thread, NULL);
#endif
for (i = 0; i < LOOPS;) {
#ifndef USE_THREADS
writer();
#endif
#ifdef USE_MUTEX
pthread_mutex_lock(&mutex);
#endif
success = (codec2_fifo_read(f, read_buf, READ_SZ) == 0);
#ifdef USE_MUTEX
pthread_mutex_unlock(&mutex);
#endif
if (success) {
for (j = 0; j < READ_SZ; j++) {
if (read_buf[j] != n_out) {
printf("error: %d %d\n", read_buf[j], n_out);
return (1);
}
n_out++;
if (n_out == N_MAX) n_out = 0;
}
i++;
}
}
#ifdef USE_THREADS
run_thread = 0;
pthread_join(awriter_thread, NULL);
#endif
printf("%d loops tested OK\n", LOOPS);
return 0;
}
int n_in = 0;
void writer(void) {
short write_buf[WRITE_SZ];
int i;
if ((FIFO_SZ - codec2_fifo_used(f)) > WRITE_SZ) {
for (i = 0; i < WRITE_SZ; i++) {
write_buf[i] = n_in++;
if (n_in == N_MAX) n_in = 0;
}
#ifdef USE_MUTEX
pthread_mutex_lock(&mutex);
#endif
codec2_fifo_write(f, write_buf, WRITE_SZ);
pthread_mutex_unlock(&mutex);
}
}
void *writer_thread(void *data) {
while (run_thread) {
writer();
}
return NULL;
}
+200
View File
@@ -0,0 +1,200 @@
/*---------------------------------------------------------------------------*\
FILE........: tfmfsk.c
AUTHOR......: Brady O'Brien
DATE CREATED: 8 February 2016
C test driver for fmfsk_mod and fmfsk_demod in fmfsk.c. Reads a file with
input bits/rf and spits out modulated/demoduladed samples and a dump of internal
state. To run unit test, see octave/tfmfsk.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/>.
*/
//#define MODEMPROBE_ENABLE
#include <stdio.h>
#include "modem_probe.h"
/* Note: This is a dirty hack to force fsk.c to compile with modem probing
* enabled */
#include "fmfsk.c"
#define ST_BITS 10000
#define ST_FS 48000
#define ST_RS 2400
#define ST_EBNO 8
#define TEST_SELF_FULL 1 /* No-arg self test */
#define TEST_MOD 2 /* Test modulator using in and out file */
#define TEST_DEMOD 3 /* Test demodulator using in and out file */
int main(int argc, char *argv[]) {
struct FMFSK *fmfsk;
int Fs, Rs;
FILE *fin, *fout;
uint8_t *bitbuf = NULL;
float *modbuf = NULL;
uint8_t *bitbufp;
float *modbufp;
size_t bitbufsize = 0;
size_t modbufsize = 0;
int test_type;
int i;
fin = NULL;
fout = NULL;
/* Set up full self-test */
if (argc == 1) {
test_type = TEST_SELF_FULL;
modem_probe_init("fmfsk", "fmfsk_tfmfsk_log.txt");
Fs = ST_FS;
Rs = ST_RS;
} else if (argc < 7) {
/* Not running any test */
printf(
"Usage: %s [(M|D) SampleRate BitRate InputFile OutputFile "
"OctaveLogFile]\n",
argv[0]);
exit(1);
} else {
/* Running stim-drivin test */
/* Mod test */
if (strcmp(argv[1], "M") == 0 || strcmp(argv[1], "m") == 0) {
test_type = TEST_MOD;
/* Demod test */
} else if (strcmp(argv[1], "D") == 0 || strcmp(argv[1], "d") == 0) {
test_type = TEST_DEMOD;
} else {
printf("Must specify mod or demod test with M or D\n");
exit(1);
}
/* Extract parameters */
Fs = atoi(argv[2]);
Rs = atoi(argv[3]);
/* Open files */
fin = fopen(argv[4], "r");
fout = fopen(argv[5], "w");
if (fin == NULL || fout == NULL) {
printf("Couldn't open test vector files\n");
exit(1);
}
/* Init modem probing */
modem_probe_init("fmfsk", argv[6]);
}
srand(1);
/* set up FSK */
fmfsk = fmfsk_create(Fs, Rs);
/* Modulate! */
if (test_type == TEST_MOD || test_type == TEST_SELF_FULL) {
/* Generate random bits for self test */
if (test_type == TEST_SELF_FULL) {
bitbufsize = ST_BITS;
bitbuf = (uint8_t *)malloc(sizeof(uint8_t) * ST_BITS);
for (i = 0; i < ST_BITS; i++) {
/* Generate a randomish bit */
bitbuf[i] = (uint8_t)(rand() & 0x01);
}
} else { /* Load bits from a file */
/* Figure out how many bits are in the input file */
fseek(fin, 0L, SEEK_END);
bitbufsize = ftell(fin);
fseek(fin, 0L, SEEK_SET);
bitbuf = malloc(sizeof(uint8_t) * bitbufsize);
i = 0;
/* Read in some bits */
bitbufp = bitbuf;
while (fread(bitbufp, sizeof(uint8_t), fmfsk->nbit, fin) == fmfsk->nbit) {
i++;
bitbufp += fmfsk->nbit;
/* Make sure we don't break the buffer */
if (i * fmfsk->nbit > bitbufsize) {
bitbuf =
realloc(bitbuf, sizeof(uint8_t) * (bitbufsize + fmfsk->nbit));
bitbufsize += fmfsk->nbit;
}
}
}
/* Allocate modulation buffer */
modbuf = (float *)malloc(sizeof(float) * (bitbufsize / fmfsk->nbit) *
fmfsk->N * 4);
modbufsize = (bitbufsize / fmfsk->nbit) * fmfsk->N;
/* Do the modulation */
modbufp = modbuf;
bitbufp = bitbuf;
while (bitbufp < bitbuf + bitbufsize) {
fmfsk_mod(fmfsk, modbufp, bitbufp);
modbufp += fmfsk->N;
bitbufp += fmfsk->nbit;
}
/* For a mod-only test, write out the result */
if (test_type == TEST_MOD) {
fwrite(modbuf, sizeof(float), modbufsize, fout);
free(modbuf);
}
/* Free bit buffer */
free(bitbuf);
}
/* Add channel imp here */
/* Now test the demod */
if (test_type == TEST_DEMOD || test_type == TEST_SELF_FULL) {
free(modbuf);
modbuf = malloc(sizeof(float) * (fmfsk->N + fmfsk->Ts * 2));
bitbuf = malloc(sizeof(uint8_t) * fmfsk->nbit);
/* Demod-only test */
if (test_type == TEST_DEMOD) {
// fprintf(stderr,"%d\n",(fmfsk->N+fmfsk->Ts*2));
while (fread(modbuf, sizeof(float), fmfsk_nin(fmfsk), fin) ==
fmfsk_nin(fmfsk)) {
fmfsk_demod(fmfsk, bitbuf, modbuf);
fwrite(bitbuf, sizeof(uint8_t), fmfsk->nbit, fout);
}
}
/* Demod after channel imp. and mod */
else {
bitbufp = bitbuf;
modbufp = modbuf;
while (modbufp < modbuf + modbufsize) {
fmfsk_demod(fmfsk, bitbuf, modbuf);
modbufp += fmfsk_nin(fmfsk);
}
}
free(bitbuf);
}
modem_probe_close();
if (test_type == TEST_DEMOD || test_type == TEST_MOD) {
fclose(fin);
fclose(fout);
}
fmfsk_destroy(fmfsk);
exit(0);
}
+113
View File
@@ -0,0 +1,113 @@
/*---------------------------------------------------------------------------*\
FILE........: tfreedv_2400A_rawdata.c
AUTHOR......: Jeroen Vreeken
DATE CREATED: 24 May 2020
FreeDV 2400A rawdata test.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2020 Jeroen Vreeken <jeroen@vreeken.net>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser 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 "assert.h"
#include "freedv_api.h"
int main(int argc, char **argv) {
struct freedv *f;
int i;
printf("freedv_api tests for mode 2400A\n");
printf("freedv_open(FREEDV_MODE_2400A) ");
f = freedv_open(FREEDV_MODE_2400A);
assert(f != NULL);
printf("Passed\n");
printf("freedv_get_mode() ");
int mode = freedv_get_mode(f);
assert(mode == FREEDV_MODE_2400A);
printf("Passed\n");
printf("freedv_get_n_max_modem_samples() ");
int max_samples = freedv_get_n_max_modem_samples(f);
assert(max_samples == 2040);
printf("%d Passed\n", max_samples);
printf("freedv_get_n_nom_modem_samples() ");
int nom_samples = freedv_get_n_nom_modem_samples(f);
assert(nom_samples == 2000);
printf("%d Passed\n", nom_samples);
printf("freedv_get_n_speech_samples() ");
int speech_samples = freedv_get_n_speech_samples(f);
assert(speech_samples == 320);
printf("%d Passed\n", speech_samples);
printf("freedv_get_n_bits_per_codec_frame() ");
int codec_bits = freedv_get_bits_per_codec_frame(f);
assert(codec_bits == 52);
printf("%d Passed\n", codec_bits);
printf("freedv_get_n_bits_per_modem_frame() ");
int frame_bits = freedv_get_bits_per_modem_frame(f);
assert(frame_bits == 52);
printf("%d Passed\n", frame_bits);
printf("freedv_rawdatatx()/freedv_rawdatarx() ");
int frames = 0;
{
short mod[nom_samples * 10];
/* Note: A codec frame is only 6.5 bytes!
so the seventh byte will be half empty!
*/
unsigned char payload[7] = {0x11, 0x22, 0x33, 0x44, 0x55, 0x66, 0x70};
for (i = 0; i < 10; i++) {
freedv_rawdatatx(f, mod + i * nom_samples, payload);
}
int nin = 0;
for (i = 0; i < nom_samples * 9; i += nin) {
nin = freedv_nin(f);
unsigned char payload_rx[7] = {0};
int r = freedv_rawdatarx(f, payload_rx, mod + i);
if (r) {
int b;
for (b = 0; b < 7; b++) {
if (payload[b] != payload_rx[b]) {
printf("Received codec bits 0x%02x do not match expected 0x%02x\n",
payload_rx[b], payload[b]);
}
}
frames++;
}
}
}
if (!frames) {
printf("Did not decode any frames successfully\n");
goto fail;
}
printf("Tests passed\n");
return 0;
fail:
printf("Test failed\n");
return 1;
}
+113
View File
@@ -0,0 +1,113 @@
/*---------------------------------------------------------------------------*\
FILE........: tfreedv_2400B_rawdata.c
AUTHOR......: Jeroen Vreeken
DATE CREATED: 24 May 2020
FreeDV 2400B rawdata test.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2020 Jeroen Vreeken <jeroen@vreeken.net>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser 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 "assert.h"
#include "freedv_api.h"
int main(int argc, char **argv) {
struct freedv *f;
int i;
printf("freedv_api tests for mode 2400B\n");
printf("freedv_open(FREEDV_MODE_2400B) ");
f = freedv_open(FREEDV_MODE_2400B);
assert(f != NULL);
printf("Passed\n");
printf("freedv_get_mode() ");
int mode = freedv_get_mode(f);
assert(mode == FREEDV_MODE_2400B);
printf("Passed\n");
printf("freedv_get_n_max_modem_samples() ");
int max_samples = freedv_get_n_max_modem_samples(f);
assert(max_samples == 1930);
printf("%d Passed\n", max_samples);
printf("freedv_get_n_nom_modem_samples() ");
int nom_samples = freedv_get_n_nom_modem_samples(f);
assert(nom_samples == 1920);
printf("%d Passed\n", nom_samples);
printf("freedv_get_n_speech_samples() ");
int speech_samples = freedv_get_n_speech_samples(f);
assert(speech_samples == 320);
printf("%d Passed\n", speech_samples);
printf("freedv_get_n_bits_per_codec_frame() ");
int codec_bits = freedv_get_bits_per_codec_frame(f);
assert(codec_bits == 52);
printf("%d Passed\n", codec_bits);
printf("freedv_get_n_bits_per_modem_frame() ");
int frame_bits = freedv_get_bits_per_modem_frame(f);
assert(frame_bits == 52);
printf("%d Passed\n", frame_bits);
printf("freedv_rawdatatx()/freedv_rawdatarx() ");
int frames = 0;
{
short mod[nom_samples * 10];
/* Note: A codec frame is only 6.5 bytes!
so the seventh byte will be half empty!
*/
unsigned char payload[7] = {0x11, 0x22, 0x33, 0x44, 0x55, 0x66, 0x70};
for (i = 0; i < 10; i++) {
freedv_rawdatatx(f, mod + i * nom_samples, payload);
}
int nin = 0;
for (i = 0; i < nom_samples * 9; i += nin) {
nin = freedv_nin(f);
unsigned char payload_rx[7] = {0};
int r = freedv_rawdatarx(f, payload_rx, mod + i);
if (r) {
int b;
for (b = 0; b < 7; b++) {
if (payload[b] != payload_rx[b]) {
printf("Received codec bits 0x%02x do not match expected 0x%02x\n",
payload_rx[b], payload[b]);
}
}
frames++;
}
}
}
if (!frames) {
printf("Did not decode any frames successfully\n");
goto fail;
}
printf("Tests passed\n");
return 0;
fail:
printf("Test failed\n");
return 1;
}
+147
View File
@@ -0,0 +1,147 @@
/*---------------------------------------------------------------------------*\
FILE........: tfreedv_800XA_rawdata.c
AUTHOR......: Jeroen Vreeken
DATE CREATED: 24 May 2020
FreeDV 800XA rawdata test.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2020 Jeroen Vreeken <jeroen@vreeken.net>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser 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 "assert.h"
#include "freedv_api.h"
int main(int argc, char **argv) {
struct freedv *f;
int i;
printf("freedv_api tests for mode 800XA\n");
printf("freedv_open(FREEDV_MODE_800XA) ");
f = freedv_open(FREEDV_MODE_800XA);
assert(f != NULL);
printf("Passed\n");
printf("freedv_get_mode() ");
int mode = freedv_get_mode(f);
assert(mode == FREEDV_MODE_800XA);
printf("Passed\n");
printf("freedv_get_n_max_modem_samples() ");
int max_samples = freedv_get_n_max_modem_samples(f);
assert(max_samples == 660);
printf("%d Passed\n", max_samples);
printf("freedv_get_n_nom_modem_samples() ");
int nom_samples = freedv_get_n_nom_modem_samples(f);
assert(nom_samples == 640);
printf("%d Passed\n", nom_samples);
printf("freedv_get_n_speech_samples() ");
int speech_samples = freedv_get_n_speech_samples(f);
assert(speech_samples == 640);
printf("%d Passed\n", speech_samples);
printf("freedv_get_n_bits_per_codec_frame() ");
int codec_bits = freedv_get_bits_per_codec_frame(f);
assert(codec_bits == 28);
printf("%d Passed\n", codec_bits);
printf("freedv_get_n_bits_per_modem_frame() ");
int frame_bits = freedv_get_bits_per_modem_frame(f);
assert(frame_bits == 56);
printf("%d Passed\n", frame_bits);
/* Note: A codec frame is only 3.5 bytes!
so the fourth and eight bytes will be half empty!
*/
unsigned char payload[8] = {0x12, 0x34, 0x56, 0x70, 0x89, 0xab, 0xcd, 0xe0};
unsigned char payload_tx[7] = {0x12, 0x34, 0x56, 0x78, 0x9a, 0xbc, 0xde};
printf("freedv_codec_frames_from_rawdata() ");
unsigned char codec_frames[8] = {0};
freedv_codec_frames_from_rawdata(f, codec_frames, payload_tx);
int fails = 0;
for (i = 0; i < 8; i++) {
if (codec_frames[i] != payload[i]) {
printf("byte %d: 0x%02x does not match expected 0x%02x\n", i,
codec_frames[i], payload[i]);
fails++;
}
}
if (fails) goto fail;
printf("Passed\n");
printf("freedv_rawdata_from_codec_frames() ");
unsigned char rawdata[7] = {0};
freedv_rawdata_from_codec_frames(f, rawdata, payload);
fails = 0;
for (i = 0; i < 7; i++) {
if (rawdata[i] != payload_tx[i]) {
printf("byte %d: 0x%02x does not match expected 0x%02x\n", i, rawdata[i],
payload_tx[i]);
fails++;
}
}
if (fails) goto fail;
printf("Passed\n");
printf("freedv_rawdatatx()/freedv_rawdatarx() ");
int frames = 0;
fails = 0;
{
short mod[nom_samples * 10];
for (i = 0; i < 10; i++) {
freedv_rawdatatx(f, mod + i * nom_samples, payload_tx);
}
int nin = 0;
for (i = 0; i < nom_samples * 9; i += nin) {
nin = freedv_nin(f);
unsigned char payload_rx[8] = {0};
int r = freedv_rawdatarx(f, payload_rx, mod + i);
if (r == 7) {
int b;
for (b = 0; b < 7; b++) {
if (payload_tx[b] != payload_rx[b]) {
printf("Received codec bits 0x%02x do not match expected 0x%02x\n",
payload_rx[b], payload_tx[b]);
fails++;
}
}
frames++;
}
}
}
if (!frames) {
printf("Did not decode any frames successfully\n");
goto fail;
}
if (fails) goto fail;
printf("Passed\n");
printf("Tests passed\n");
return 0;
fail:
printf("Test failed\n");
return 1;
}
+282
View File
@@ -0,0 +1,282 @@
/*---------------------------------------------------------------------------*\
FILE........: tfreedv_data_channel
AUTHOR......: Jeroen Vreeken
DATE CREATED: May 3 2016
Tests for the data channel code.
Data channel frame behaviour is tested with test vectors.
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2016 Jeroen Vreeken
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 "freedv_data_channel.h"
unsigned char test_header[] = {0x11, 0x22, 0x33, 0x44, 0x55, 0x66};
unsigned char bcast_header[] = {0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
struct testvec {
char *testname;
unsigned char *data;
size_t data_size;
size_t frame_size;
unsigned char *frame_data;
size_t frame_data_size;
unsigned char *flags;
} testvec[] = {
{
"Regular packet, does not match header and no broadcast",
(unsigned char[]){0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09,
0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10, 0x11, 0x12},
0x12,
8,
(unsigned char[]){0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x01,
0x02, 0x03, 0x04, 0x05, 0x06, 0x0d, 0x0e,
0x0f, 0x10, 0x11, 0x12, 0x47, 0x6e},
0x14,
(unsigned char[]){0x00, 0x00, 0x04},
},
{
"Header",
NULL,
0,
8,
(unsigned char[]){0x11, 0x22, 0x33, 0x44, 0x55, 0x66, 0x5a, 0x60},
0x08,
(unsigned char[]){0x08},
},
{
"Broadcast packet",
(unsigned char[]){0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x11, 0x22, 0x33,
0x44, 0x55, 0x66, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a,
0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10, 0x11},
0x19,
8,
(unsigned char[]){0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d,
0x0e, 0x0f, 0x10, 0x11, 0x3c, 0xbe},
0x0f,
(unsigned char[]){0xc0, 0x07},
},
{
"Broadcast packet, header does not match",
(unsigned char[]){0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xaa, 0x22, 0xbb,
0xcc, 0xdd, 0xee, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a,
0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10, 0x11},
0x19,
8,
(unsigned char[]){0xaa, 0x22, 0xbb, 0xcc, 0xdd, 0xee, 0x05,
0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c,
0x0d, 0x0e, 0x0f, 0x10, 0x11, 0x1a, 0x68},
0x15,
(unsigned char[]){0x40, 0x00, 0x05},
},
{
"Header 6 bytes",
NULL,
0,
6,
(unsigned char[]){0x11, 0x22, 0x33, 0x44, 0x55, 0x66},
0x06,
(unsigned char[]){0x2f},
},
{
"Broadcast packet (6 byte frames)",
(unsigned char[]){0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x11, 0x22, 0x33,
0x44, 0x55, 0x66, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a,
0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10, 0x11},
0x19,
6,
(unsigned char[]){0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d,
0x0e, 0x0f, 0x10, 0x11, 0x3c, 0xbe},
0x0f,
(unsigned char[]){0xc0, 0x00, 0x03},
},
{
"Broadcast packet, header does not match (6 byte frames)",
(unsigned char[]){0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xaa, 0x22, 0xbb,
0xcc, 0xdd, 0xee, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a,
0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10, 0x11},
0x19,
6,
(unsigned char[]){0xaa, 0x22, 0xbb, 0xcc, 0xdd, 0xee, 0x05,
0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c,
0x0d, 0x0e, 0x0f, 0x10, 0x11, 0x1a, 0x68},
0x15,
(unsigned char[]){0x40, 0x00, 0x00, 0x03},
},
};
static int ret = 0;
static int vector = 0;
static size_t frame_data_pos = 0;
static int rx_done = 0;
void *tx_cb_arg = (void *)0xaa55;
void *rx_cb_arg = (void *)0xbb44;
void tfreedv_data_callback_tx(void *arg, unsigned char *packet, size_t *size) {
if (tx_cb_arg != arg) {
ret++;
printf("FAIL: %s called with wrong argument value\n", __func__);
}
printf("--------------------------------------\n");
printf("TX callback called for test %zd bytes data for test %d:\n'%s'\n",
testvec[vector].data_size, vector, testvec[vector].testname);
memcpy(packet, testvec[vector].data, testvec[vector].data_size);
*size = testvec[vector].data_size;
return;
}
void tfreedv_data_callback_rx(void *arg, unsigned char *packet, size_t size) {
if (rx_cb_arg != arg) {
ret++;
printf("FAIL: %s called with wrong argument value\n", __func__);
}
printf("RX callback called with %zd bytes\n", size);
if (testvec[vector].data_size) {
size_t data_size = testvec[vector].data_size;
if (data_size != size) {
printf("FAIL: Received size does not match test vector: %zd != %zd\n",
size, data_size);
ret++;
} else {
size_t i;
for (i = 0; i < data_size; i++) {
if (packet[i] != testvec[vector].data[i]) {
printf("FAIL: byte %zd does not match 0x%02x != 0x%02x\n", i,
packet[i], testvec[vector].data[i]);
ret++;
}
}
}
} else {
if (size != 12) {
printf("FAIL: Received header is not 12 bytes: %zd\n", size);
ret++;
} else {
if (memcmp(packet, bcast_header, 6)) {
printf("FAIL: Header is not a broadcast!\n");
ret++;
}
if (memcmp(packet + 6, test_header, 6)) {
printf("FAIL: Header does not match!\n");
ret++;
}
}
}
rx_done = 1;
}
int main(int argc, char **argv) {
struct freedv_data_channel *fdc;
fdc = freedv_data_channel_create();
freedv_data_set_header(fdc, test_header);
freedv_data_set_cb_tx(fdc, tfreedv_data_callback_tx, tx_cb_arg);
freedv_data_set_cb_rx(fdc, tfreedv_data_callback_rx, rx_cb_arg);
while (vector < sizeof(testvec) / sizeof(struct testvec)) {
size_t frame_size = testvec[vector].frame_size;
unsigned char frame[frame_size];
int from, bcast, crc, end;
int i;
size_t check_size;
unsigned char flags;
int nr_frames;
freedv_data_channel_tx_frame(fdc, frame, frame_size, &from, &bcast, &crc,
&end);
check_size = frame_size;
if (frame_data_pos + check_size > testvec[vector].frame_data_size)
check_size = testvec[vector].frame_data_size - frame_data_pos;
flags = from * 0x80 + bcast * 0x40 + crc * 0x20 + end;
printf("0x%02x:", flags);
for (i = 0; i < check_size; i++) {
if (frame[i] != testvec[vector].frame_data[frame_data_pos + i]) {
printf(" [0x%02x!=0x%02x]", frame[i],
testvec[vector].frame_data[frame_data_pos + i]);
ret++;
} else {
printf(" 0x%02x", frame[i]);
}
}
printf("\n");
if (flags != testvec[vector].flags[frame_data_pos / frame_size]) {
printf("FAIL: Flags byte does not match 0x%02x != 0x%02x\n", flags,
testvec[vector].flags[frame_data_pos / frame_size]);
ret++;
}
freedv_data_channel_rx_frame(fdc, frame, frame_size, from, bcast, crc, end);
frame_data_pos += frame_size;
nr_frames = freedv_data_get_n_tx_frames(fdc, frame_size);
if (frame_data_pos >= testvec[vector].frame_data_size) {
if (nr_frames) {
printf("FAIL: nr_frames is not zero: %d\n", nr_frames);
ret++;
}
vector++;
frame_data_pos = 0;
if (!rx_done) {
printf("FAIL: RX callback not executed\n");
ret++;
}
rx_done = 0;
} else {
int vec_frames = (testvec[vector].frame_data_size - frame_data_pos);
vec_frames /= frame_size;
vec_frames++;
if (nr_frames != vec_frames) {
printf("FAIL: nr_frames != vec_frames: %d != %d\n", nr_frames,
vec_frames);
ret++;
}
}
}
freedv_data_channel_destroy(fdc);
printf("--------------------------------------\n");
printf("tfreedv_data_channel test result: ");
if (ret) {
printf("Failed %d\n", ret);
} else {
printf("Passed\n");
}
return ret;
}
+234
View File
@@ -0,0 +1,234 @@
/*---------------------------------------------------------------------------*\
FILE........: tfsk.c
AUTHOR......: Brady O'Brien
DATE CREATED: 20 January 2016
C test driver for fsk_mod and fsk_demod in fsk.c. Reads a file with input
bits/rf and spits out modulated/demoduladed samples and a dump of internal
state. To run unit test, see octave/tfsk.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/>.
*/
//#define MODEMPROBE_ENABLE
#include <stdio.h>
#include "modem_probe.h"
/* Note: This is a dirty hack to force fsk.c to compile with modem probing
* enabled */
#include "fsk.c"
#define ST_BITS 10000
#define ST_FS 8000
#define ST_RS 100
#define ST_F1 1200
#define ST_Fs 400
#define ST_EBNO 8
#define ST_M 2
#define TEST_SELF_FULL 1 /* No-arg self test */
#define TEST_MOD 2 /* Test modulator using in and out file */
#define TEST_MOD_H 3 /* Test modulator using in and out file */
#define TEST_DEMOD 4 /* Test demodulator using in and out file */
#define TEST_DEMOD_H 5 /* Test demodulator using in and out file */
int main(int argc, char *argv[]) {
struct FSK *fsk;
int Fs, Rs, f1, fs, M, lock_nin;
FILE *fin, *fout;
uint8_t *bitbuf = NULL;
float *modbuf = NULL;
uint8_t *bitbufp;
float *modbufp;
size_t bitbufsize = 0;
size_t modbufsize = 0;
int test_type;
int i;
fin = NULL;
fout = NULL;
/* Set up full self-test */
if (argc == 1) {
test_type = TEST_SELF_FULL;
modem_probe_init("fsk2", "fsk2_tfsk_log.txt");
Fs = ST_FS;
Rs = ST_RS;
f1 = ST_F1;
fs = ST_Fs;
M = ST_M;
lock_nin = 0;
} else if (argc < 10) {
/* Not running any test */
printf(
"Usage: %s [(M|D|DX) Mode TXFreq1 TXFreqSpace SampleRate SymbolRate "
"lock_nin InputFile OutputFile OctaveLogFile]\n",
argv[0]);
exit(1);
} else {
/* Running stim-drivin test */
/* Mod test */
if (strcmp(argv[1], "MX") == 0) {
test_type = TEST_MOD_H;
} else if (strcmp(argv[1], "M") == 0 || strcmp(argv[1], "m") == 0) {
test_type = TEST_MOD;
/* Demod test */
} else if (strcmp(argv[1], "DX") == 0) {
test_type = TEST_DEMOD_H;
} else if (strcmp(argv[1], "D") == 0 || strcmp(argv[1], "d") == 0) {
test_type = TEST_DEMOD;
} else {
printf("Must specify mod or demod test with M or D\n");
exit(1);
}
/* Extract parameters */
M = atoi(argv[2]);
f1 = atoi(argv[3]);
fs = atoi(argv[4]);
Fs = atoi(argv[5]);
Rs = atoi(argv[6]);
lock_nin = atoi(argv[7]);
/* Open files */
fin = fopen(argv[8], "r");
fout = fopen(argv[9], "w");
if (fin == NULL || fout == NULL) {
printf("Couldn't open test vector files\n");
exit(1);
}
/* Init modem probing */
modem_probe_init("fsk", argv[10]);
}
srand(1);
/* set up FSK */
if (test_type == TEST_DEMOD_H || test_type == TEST_MOD_H) {
fsk = fsk_create_hbr(Fs, Rs, M, 10, FSK_DEFAULT_NSYM, f1, fs);
if (test_type == TEST_DEMOD_H)
test_type = TEST_DEMOD;
else
test_type = TEST_MOD;
} else {
fsk = fsk_create(Fs, Rs, M, f1, fs);
}
fsk_set_freq_est_limits(fsk, 300, 2800);
fsk->lock_nin = lock_nin;
/* Modulate! */
if (test_type == TEST_MOD || test_type == TEST_SELF_FULL) {
/* Generate random bits for self test */
if (test_type == TEST_SELF_FULL) {
bitbufsize = ST_BITS;
bitbuf = (uint8_t *)malloc(sizeof(uint8_t) * ST_BITS);
for (i = 0; i < ST_BITS; i++) {
/* Generate a randomish bit */
bitbuf[i] = (uint8_t)(rand() & 0x01);
}
} else { /* Load bits from a file */
/* Figure out how many bits are in the input file */
fseek(fin, 0L, SEEK_END);
bitbufsize = ftell(fin);
fseek(fin, 0L, SEEK_SET);
bitbuf = malloc(sizeof(uint8_t) * bitbufsize);
i = 0;
/* Read in some bits */
bitbufp = bitbuf;
while (fread(bitbufp, sizeof(uint8_t), fsk->Nbits, fin) == fsk->Nbits) {
i++;
bitbufp += fsk->Nbits;
/* Make sure we don't break the buffer */
if (i * fsk->Nbits > bitbufsize) {
bitbuf = realloc(bitbuf, sizeof(uint8_t) * (bitbufsize + fsk->Nbits));
bitbufsize += fsk->Nbits;
}
}
}
/* Allocate modulation buffer */
modbuf =
(float *)malloc(sizeof(float) * (bitbufsize / fsk->Nbits) * fsk->N * 4);
modbufsize = (bitbufsize / fsk->Nbits) * fsk->N;
/* Do the modulation */
modbufp = modbuf;
bitbufp = bitbuf;
while (bitbufp < bitbuf + bitbufsize) {
fsk_mod(fsk, modbufp, bitbufp, fsk->Nbits);
modbufp += fsk->N;
bitbufp += fsk->Nbits;
}
/* For a mod-only test, write out the result */
if (test_type == TEST_MOD) {
fwrite(modbuf, sizeof(float), modbufsize, fout);
free(modbuf);
}
/* Free bit buffer */
free(bitbuf);
}
/* Now test the demod */
if (test_type == TEST_DEMOD || test_type == TEST_SELF_FULL) {
free(modbuf);
modbuf = malloc(sizeof(float) * (fsk->N + fsk->Ts * 2));
bitbuf = malloc(sizeof(uint8_t) * fsk->Nbits);
/* Demod-only test */
if (test_type == TEST_DEMOD) {
while (fread(modbuf, sizeof(float), fsk_nin(fsk), fin) == fsk_nin(fsk)) {
int n = fsk_nin(fsk);
COMP modbuf_comp[n];
for (i = 0; i < n; i++) {
modbuf_comp[i].real = modbuf[i];
modbuf_comp[i].imag = 0.0;
}
fsk_demod(fsk, bitbuf, modbuf_comp);
fwrite(bitbuf, sizeof(uint8_t), fsk->Nbits, fout);
}
}
/* Demod after channel imp. and mod */
else {
bitbufp = bitbuf;
modbufp = modbuf;
while (modbufp < modbuf + modbufsize) {
int n = fsk_nin(fsk);
COMP modbuf_comp[n];
for (i = 0; i < n; i++) {
modbuf_comp[i].real = modbuf[i];
modbuf_comp[i].imag = 0.0;
}
fsk_demod(fsk, bitbuf, modbuf_comp);
modbufp += fsk_nin(fsk);
}
}
free(bitbuf);
}
modem_probe_close();
if (test_type == TEST_DEMOD || test_type == TEST_MOD) {
fclose(fin);
fclose(fout);
}
fsk_destroy(fsk);
exit(0);
}
+60
View File
@@ -0,0 +1,60 @@
/*---------------------------------------------------------------------------*\
FILE........: tfsk_llr.c
AUTHOR......: David Rowe
DATE CREATED: July 2020
Simple test program for 4FSK LLR routines.
\*---------------------------------------------------------------------------*/
#include <math.h>
#include <stdio.h>
#include "mpdecode_core.h"
#define M 4
#define BPS 2
#define NSYM 5
#define V_EST 2
#define SNR_EST 10
/* Generated test vectors with:
octave:100> init_cml('~/cml/');
octave:101> rx_filt=[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
octave:102> symL = DemodFSK(rx_filt,10,1); -Somap(symL)
*/
/* one col per symbol:
0 1 2 3 4 */
float rx_filt[] = {
1.0, 0.0, 0.0, 0.0, 1.0, /* filter 0 */
0.0, 1.0, 0.0, 0.0, 0.0, /* filter 1 */
0.0, 0.0, 1.0, 0.0, 0.0, /* filter 2 */
0.0, 0.0, 0.0, 1.0, 0.0 /* filter 3 */
};
float llrs_target[] = {7.3252, 7.3252, /* bit 0, 1 */
7.3252, -7.3252, /* 2, 3, ... */
-7.3252, 7.3252, -7.3252, -7.3252, 7.3252, 7.3252};
int main(void) {
float llrs[BPS * NSYM] = {0};
fsk_rx_filt_to_llrs(llrs, rx_filt, V_EST, SNR_EST, M, NSYM);
float error = 0.0;
for (int i = 0; i < NSYM * BPS; i++) {
fprintf(stderr, "% f\n", llrs[i]);
error += pow(llrs[i] - llrs_target[i], 2.0);
}
if (error < 1E-3) {
fprintf(stderr, "PASS\n");
return 0;
} else {
fprintf(stderr, "FAIL\n");
return 1;
}
}
+18
View File
@@ -0,0 +1,18 @@
/*---------------------------------------------------------------------------*\
FILE........: thash.c
AUTHOR......: David Rowe
DATE CREATED: July 2020
Simple test program for freeDV API get hash function
\*---------------------------------------------------------------------------*/
#include <stdio.h>
#include "freedv_api.h"
int main(void) {
printf("%s\n", freedv_get_hash());
return 0;
}
+295
View File
@@ -0,0 +1,295 @@
/*---------------------------------------------------------------------------*\
FILE........: tnewamp1.c
AUTHOR......: David Rowe
DATE CREATED: Jan 2017
Tests for the C version of the newamp1 amplitude modelling used for
700c. This program outputs a file of Octave vectors that are loaded
and automatically tested against the Octave version of the modem by
the Octave script tnewamp1.m
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2017 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 "codec2_fft.h"
#include "defines.h"
#include "dump.h"
#include "newamp1.h"
#include "nlp.h"
#include "octave.h"
#include "quantise.h"
#include "sine.h"
#define FRAMES 300
int main(int argc, char *argv[]) {
int Fs = 8000;
C2CONST c2const = c2const_create(Fs, N_S);
int n_samp = c2const.n_samp;
int m_pitch = c2const.m_pitch;
short buf[n_samp]; /* input/output buffer */
float Sn[m_pitch]; /* float input speech samples */
COMP Sw[FFT_ENC]; /* DFT of Sn[] */
codec2_fft_cfg fft_fwd_cfg; /* fwd FFT states */
float w[m_pitch]; /* time domain hamming window */
float W[FFT_ENC]; /* DFT of w[] */
MODEL model;
void *nlp_states;
codec2_fft_cfg phase_fft_fwd_cfg, phase_fft_inv_cfg;
float pitch, prev_f0;
int i, m, f, k;
if (argc != 2) {
printf("usage: ./tnewamp1 RawFile\n");
exit(1);
}
nlp_states = nlp_create(&c2const);
prev_f0 = 1.0 / P_MAX_S;
fft_fwd_cfg = codec2_fft_alloc(FFT_ENC, 0, NULL, NULL);
make_analysis_window(&c2const, fft_fwd_cfg, w, W);
phase_fft_fwd_cfg = codec2_fft_alloc(NEWAMP1_PHASE_NFFT, 0, NULL, NULL);
phase_fft_inv_cfg = codec2_fft_alloc(NEWAMP1_PHASE_NFFT, 1, NULL, NULL);
for (i = 0; i < m_pitch; i++) {
Sn[i] = 1.0;
}
int K = 20;
float rate_K_sample_freqs_kHz[K];
float model_octave[FRAMES][MAX_AMP + 2]; // model params in matrix format,
// useful for C <-> Octave
float rate_K_surface[FRAMES][K]; // rate K vecs for each frame, form a
// surface that makes pretty graphs
float rate_K_surface_no_mean[FRAMES][K]; // mean removed surface
float rate_K_surface_no_mean_[FRAMES][K]; // quantised mean removed surface
float mean[FRAMES];
float mean_[FRAMES];
float rate_K_surface_[FRAMES][K]; // quantised rate K vecs for each frame
float interpolated_surface_[FRAMES][K]; // dec/interpolated surface
// int voicing[FRAMES];
int voicing_[FRAMES];
float model_octave_[FRAMES][MAX_AMP + 2];
COMP H[FRAMES][MAX_AMP];
int indexes[FRAMES][NEWAMP1_N_INDEXES];
float se = 0.0;
float eq[K];
for (k = 0; k < K; k++) eq[k] = 0.0;
for (f = 0; f < FRAMES; f++) {
for (m = 0; m < MAX_AMP + 2; m++) {
model_octave[f][m] = 0.0;
model_octave_[f][m] = 0.0;
}
for (m = 0; m < MAX_AMP; m++) {
H[f][m].real = 0.0;
H[f][m].imag = 0.0;
}
for (k = 0; k < K; k++) interpolated_surface_[f][k] = 0.0;
voicing_[f] = 0;
}
mel_sample_freqs_kHz(rate_K_sample_freqs_kHz, K, ftomel(200.0),
ftomel(3700.0));
// for(int k=0; k<K; k++)
// printf("k: %d sf: %f\n", k, rate_K_sample_freqs_kHz[k]);
FILE *fin = fopen(argv[1], "rb");
if (fin == NULL) {
fprintf(stderr, "Problem opening hts1.raw\n");
exit(1);
}
int M = 4;
for (f = 0; f < FRAMES; f++) {
assert(fread(buf, sizeof(short), n_samp, fin) == n_samp);
/* shift buffer of input samples, and insert new samples */
for (i = 0; i < m_pitch - n_samp; i++) {
Sn[i] = Sn[i + n_samp];
}
for (i = 0; i < n_samp; i++) {
Sn[i + m_pitch - n_samp] = buf[i];
}
/* Estimate Sinusoidal Model Parameters ----------------------*/
nlp(nlp_states, Sn, n_samp, &pitch, Sw, W, &prev_f0);
model.Wo = TWO_PI / pitch;
dft_speech(&c2const, fft_fwd_cfg, Sw, Sn, w);
two_stage_pitch_refinement(&c2const, &model, Sw);
estimate_amplitudes(&model, Sw, W, 1);
est_voicing_mbe(&c2const, &model, Sw, W);
// voicing[f] = model.voiced;
/* newamp1 processing ----------------------------------------*/
newamp1_model_to_indexes(&c2const, &indexes[f][0], &model,
&rate_K_surface[f][0], rate_K_sample_freqs_kHz, K,
&mean[f], &rate_K_surface_no_mean[f][0],
&rate_K_surface_no_mean_[f][0], &se, eq, 0);
newamp1_indexes_to_rate_K_vec(
&rate_K_surface_[f][0], &rate_K_surface_no_mean_[f][0],
rate_K_sample_freqs_kHz, K, &mean_[f], &indexes[f][0], NULL, 1);
#ifdef VERBOSE
fprintf(stderr, "f: %d Wo: %4.3f L: %d v: %d\n", f, model.Wo, model.L,
model.voiced);
if ((f % M) == 0) {
for (i = 0; i < 5; i++) {
fprintf(stderr, " %5.3f", rate_K_surface_[f][i]);
}
fprintf(stderr, "\n");
fprintf(stderr, " %d %d %d %d\n", indexes[f][0], indexes[f][1],
indexes[f][2], indexes[f][3]);
}
#endif
/* log vectors */
model_octave[f][0] = model.Wo;
model_octave[f][1] = model.L;
for (m = 1; m <= model.L; m++) {
model_octave[f][m + 1] = model.A[m];
}
}
/* Decoder */
MODEL model__[M];
float prev_rate_K_vec_[K];
COMP HH[M][MAX_AMP + 1];
float Wo_left;
int voicing_left;
/* initial conditions */
for (k = 0; k < K; k++) prev_rate_K_vec_[k] = 0.0;
voicing_left = 0;
Wo_left = 2.0 * M_PI / 100.0;
/* decoder runs on every M-th frame, 25Hz frame rate, offset at
start is to minimise processing delay (thanks Jeroen!) */
fprintf(stderr, "\n");
for (f = M - 1; f < FRAMES; f += M) {
float a_interpolated_surface_[M][K];
newamp1_indexes_to_model(
&c2const, model__, (COMP *)HH, (float *)a_interpolated_surface_,
prev_rate_K_vec_, &Wo_left, &voicing_left, rate_K_sample_freqs_kHz, K,
phase_fft_fwd_cfg, phase_fft_inv_cfg, &indexes[f][0], NULL, 1);
#ifdef VERBOSE
fprintf(stderr, "f: %d\n", f);
fprintf(stderr, " %d %d %d %d\n", indexes[f][0], indexes[f][1],
indexes[f][2], indexes[f][3]);
for (i = 0; i < M; i++) {
fprintf(stderr, " Wo: %4.3f L: %d v: %d\n", model__[i].Wo, model__[i].L,
model__[i].voiced);
}
fprintf(stderr, " rate_K_vec: ");
for (i = 0; i < 5; i++) {
fprintf(stderr, "%5.3f ", prev_rate_K_vec_[i]);
}
fprintf(stderr, "\n");
fprintf(stderr, " Am H:\n");
for (m = 0; m < M; m++) {
fprintf(stderr, " ");
for (i = 1; i <= 5; i++) {
fprintf(stderr, "%5.1f (%5.3f %5.3f) ", model__[m].A[i], HH[m][i].real,
HH[m][i].imag);
}
fprintf(stderr, "\n");
}
fprintf(stderr, "\n\n");
#endif
// if (f == 7)
// exit(0);
/* with f == 0, we don't store output, but memories are updated, helps to
match what happens in Codec 2 mode */
if (f >= M) {
for (i = 0; i < M; i++) {
for (k = 0; k < K; k++) {
interpolated_surface_[f - M + i][k] = a_interpolated_surface_[i][k];
}
}
/* store test vectors */
for (i = f - M, m = 0; i < f; i++, m++) {
model_octave_[i][0] = model__[m].Wo;
model_octave_[i][1] = model__[m].L;
voicing_[i] = model__[m].voiced;
}
int j;
for (i = f - M, j = 0; i < f; i++, j++) {
for (m = 1; m <= model__[j].L; m++) {
model_octave_[i][m + 1] = model__[j].A[m];
H[i][m - 1] = HH[j][m]; // aH[m];
}
}
}
}
fclose(fin);
/* save vectors in Octave format */
FILE *fout = fopen("tnewamp1_out.txt", "wt");
assert(fout != NULL);
fprintf(fout, "# Created by tnewamp1.c\n");
octave_save_float(fout, "rate_K_surface_c", (float *)rate_K_surface, FRAMES,
K, K);
octave_save_float(fout, "mean_c", (float *)mean, 1, FRAMES, 1);
octave_save_float(fout, "eq_c", eq, 1, K, K);
octave_save_float(fout, "rate_K_surface_no_mean_c",
(float *)rate_K_surface_no_mean, FRAMES, K, K);
octave_save_float(fout, "rate_K_surface_no_mean__c",
(float *)rate_K_surface_no_mean_, FRAMES, K, K);
octave_save_float(fout, "mean__c", (float *)mean_, FRAMES, 1, 1);
octave_save_float(fout, "rate_K_surface__c", (float *)rate_K_surface_, FRAMES,
K, K);
octave_save_float(fout, "interpolated_surface__c",
(float *)interpolated_surface_, FRAMES, K, K);
octave_save_float(fout, "model_c", (float *)model_octave, FRAMES, MAX_AMP + 2,
MAX_AMP + 2);
octave_save_float(fout, "model__c", (float *)model_octave_, FRAMES,
MAX_AMP + 2, MAX_AMP + 2);
octave_save_int(fout, "voicing__c", (int *)voicing_, 1, FRAMES);
octave_save_complex(fout, "H_c", (COMP *)H, FRAMES, MAX_AMP, MAX_AMP);
fclose(fout);
printf(
"Done! Now run\n octave:1> "
"tnewamp1(\"../path/to/build_linux/src/hts1a\", "
"\"../path/to/build_linux/unittest\")\n");
return 0;
}
+626
View File
@@ -0,0 +1,626 @@
/*---------------------------------------------------------------------------*\
FILE........: tofdm.c
AUTHORS.....: David Rowe & Steve Sampson
DATE CREATED: June 2017
Tests for the C version of the OFDM modem. This program outputs a
file of Octave vectors that are loaded and automatically tested
against the Octave version of the modem by the Octave script tofdm.m
\*---------------------------------------------------------------------------*/
/*
Copyright (C) 2017 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, 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 <complex.h>
#include <getopt.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "HRA_112_112.h" /* generated by ldpc_fsk_lib.m:ldpc_decode() */
#include "codec2_ofdm.h"
#include "comp_prim.h"
#include "mpdecode_core.h"
#include "octave.h"
#include "ofdm_internal.h"
#include "test_bits_ofdm.h"
#define NFRAMES 10
#define SAMPLE_CLOCK_OFFSET_PPM 100
#define FOFF_HZ 0.5f
#define ASCALE (2E5f * 1.1491f / 2.0f) /* scale from shorts back to floats */
#define CODED_BITSPERFRAME 224 /* number of LDPC codeword bits/frame */
/* 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}};
/* constants we use a lot and don't want to have to deference all the time */
static float ofdm_tx_centre; /* TX Center frequency */
static float ofdm_rx_centre; /* RX Center frequency */
static float ofdm_fs; /* Sample rate */
static float ofdm_ts; /* Symbol cycle time */
static float ofdm_rs; /* Symbol rate */
static float ofdm_tcp; /* Cyclic prefix duration */
static float
ofdm_timing_mx_thresh; /* See 700D Part 4 Acquisition blog post and
ofdm_dev.m routines for how this was set */
static int ofdm_nc; /* NS-1 data symbols between pilots */
static int ofdm_ns;
static int ofdm_bps; /* Bits per symbol */
static int ofdm_m; /* duration of each symbol in samples */
static int ofdm_ncp; /* duration of CP in samples */
static int ofdm_ftwindowwidth;
static int ofdm_bitsperframe;
static int ofdm_rowsperframe;
static int ofdm_samplesperframe;
static int ofdm_samplespersymbol;
static int ofdm_max_samplesperframe;
static int ofdm_nrxbuf;
static int
ofdm_ntxtbits; /* reserve bits/frame for auxiliary text information */
static int ofdm_nuwbits; /* Unique word, used for positive indication of lock */
/*---------------------------------------------------------------------------*\
FUNCTION....: fs_offset()
AUTHOR......: David Rowe
DATE CREATED: May 2015
Simulates small Fs offset between mod and demod.
(Note: Won't work with float, works OK with double)
\*---------------------------------------------------------------------------*/
static int fs_offset(COMP out[], COMP in[], int n, float sample_rate_ppm) {
double f;
double tin = 0.0;
double step = 1.0 + sample_rate_ppm / 1E6;
int t1, t2;
int tout = 0;
while (tin < (double)(n - 1)) {
t1 = (int)floor(tin);
t2 = (int)ceil(tin);
assert(t2 < n);
f = tin - (double)t1;
out[tout].real =
((double)1.0 - f) * (double)in[t1].real + f * (double)in[t2].real;
out[tout].imag =
((double)1.0 - f) * (double)in[t1].imag + f * (double)in[t2].imag;
tin += step;
tout++;
// printf("tin: %f tout: %d f: %f\n", tin, tout, f);
}
return tout;
}
/*---------------------------------------------------------------------------*\
FUNCTION....: freq_shift()
AUTHOR......: David Rowe
DATE CREATED: 26/4/2012
Frequency shift modem signal. The use of complex input and output allows
single sided frequency shifting (no images).
\*---------------------------------------------------------------------------*/
static void freq_shift(COMP rx_fdm_fcorr[], COMP rx_fdm[], float foff,
COMP *foff_phase_rect, int nin) {
float temp = (TAU * foff / ofdm_fs);
COMP foff_rect = {cosf(temp), sinf(temp)};
int i;
for (i = 0; i < nin; i++) {
*foff_phase_rect = cmult(*foff_phase_rect, foff_rect);
rx_fdm_fcorr[i] = cmult(rx_fdm[i], *foff_phase_rect);
}
/* normalise digital oscillator as the magnitude can drift over time */
float mag = cabsolute(*foff_phase_rect);
foff_phase_rect->real /= mag;
foff_phase_rect->imag /= mag;
}
int main(int argc, char *argv[]) {
int opt_Nc = 0;
int ldpc_enable = 1;
struct OFDM *ofdm;
struct OFDM_CONFIG *ofdm_config;
static struct option long_options[] = {{"nc", required_argument, 0, 'n'},
{"noldpc", no_argument, 0, 'l'},
{0, 0, 0, 0}};
int opt_index = 0;
char c;
while ((c = getopt_long(argc, argv, "n:l", long_options, &opt_index)) != -1) {
switch (c) {
case 'n':
opt_Nc = atoi(optarg);
fprintf(stderr, "Nc = %d\n", opt_Nc);
break;
case 'l':
ldpc_enable = 0;
fprintf(stderr, "LDPC disabled\n");
break;
default:
fprintf(stderr,
"usage: %s [Options]:\n [-l --noldpc]\n [-n --nc "
"NumberoFCarriers]\n",
argv[0]);
exit(1);
}
}
// init once to get a copy of default config params
ofdm = ofdm_create(NULL);
assert(ofdm != NULL);
struct OFDM_CONFIG ofdm_config_default;
memcpy(&ofdm_config_default, ofdm_get_config_param(ofdm),
sizeof(struct OFDM_CONFIG));
ofdm_destroy(ofdm);
// now do a little customisation on default config, and re-create modem
// instance
if (opt_Nc) ofdm_config_default.nc = opt_Nc;
// printf("ofdm_create() 2\n");
ofdm = ofdm_create(&ofdm_config_default);
assert(ofdm != NULL);
ofdm_config = ofdm_get_config_param(ofdm);
ofdm_set_tx_bpf(ofdm, false);
// same levels as Octave sim
ofdm->amp_scale = 1.0;
// make local copies for convenience
ofdm_tx_centre = ofdm_config->tx_centre;
ofdm_rx_centre = ofdm_config->rx_centre;
ofdm_fs = ofdm_config->fs;
ofdm_ts = ofdm_config->ts;
ofdm_rs = ofdm_config->rs;
ofdm_tcp = ofdm_config->tcp;
ofdm_timing_mx_thresh = ofdm_config->timing_mx_thresh;
ofdm_nc = ofdm_config->nc;
ofdm_ns = ofdm_config->ns;
ofdm_bps = ofdm_config->bps;
ofdm_m = (int)(ofdm_config->fs / ofdm_config->rs);
ofdm_ncp = (int)(ofdm_config->tcp * ofdm_config->fs);
ofdm_ftwindowwidth = ofdm_config->ftwindowwidth;
ofdm_bitsperframe = ofdm_get_bits_per_frame(ofdm);
ofdm_rowsperframe = ofdm_bitsperframe / (ofdm_config->nc * ofdm_config->bps);
ofdm_samplesperframe = ofdm_get_samples_per_frame(ofdm);
ofdm_samplespersymbol = (ofdm->m + ofdm->ncp);
ofdm_max_samplesperframe = ofdm_get_max_samples_per_frame(ofdm);
ofdm_nrxbuf = ofdm->nrxbuf;
ofdm_ntxtbits = ofdm_config->txtbits;
ofdm_nuwbits = ofdm_config->nuwbits;
int tx_bits[ofdm_samplesperframe];
COMP tx[ofdm_samplesperframe]; /* one frame of tx samples */
int rx_bits[ofdm_bitsperframe]; /* one frame of rx bits */
printf("Nc = %d ofdm_bitsperframe: %d\n", ofdm_nc, ofdm_bitsperframe);
/* log arrays */
int tx_bits_log[ofdm_bitsperframe * NFRAMES];
COMP tx_log[ofdm_samplesperframe * NFRAMES];
COMP rx_log[ofdm_samplesperframe * NFRAMES];
COMP rxbuf_in_log[ofdm_max_samplesperframe * NFRAMES];
COMP rxbuf_log[ofdm_nrxbuf * NFRAMES];
COMP rx_sym_log[(ofdm_ns + 3) * NFRAMES][ofdm_nc + 2];
float phase_est_pilot_log[ofdm_rowsperframe * NFRAMES][ofdm_nc];
COMP rx_np_log[ofdm_rowsperframe * ofdm_nc * NFRAMES];
float rx_amp_log[ofdm_rowsperframe * ofdm_nc * NFRAMES];
float foff_hz_log[NFRAMES];
int rx_bits_log[ofdm_bitsperframe * NFRAMES];
int timing_est_log[NFRAMES];
int timing_valid_log[NFRAMES];
float timing_mx_log[NFRAMES];
float coarse_foff_est_hz_log[NFRAMES];
int sample_point_log[NFRAMES];
float symbol_likelihood_log[(CODED_BITSPERFRAME / ofdm_bps) *
(1 << ofdm_bps) * NFRAMES];
float bit_likelihood_log[CODED_BITSPERFRAME * NFRAMES];
int detected_data_log[CODED_BITSPERFRAME * NFRAMES];
float mean_amp_log[NFRAMES];
float snr_log[NFRAMES];
FILE *fout;
int f, i, j;
/* set up LDPC code */
struct LDPC ldpc;
ldpc.max_iter = HRA_112_112_MAX_ITER;
ldpc.dec_type = 0;
ldpc.q_scale_factor = 1;
ldpc.r_scale_factor = 1;
ldpc.CodeLength = HRA_112_112_CODELENGTH;
ldpc.NumberParityBits = HRA_112_112_NUMBERPARITYBITS;
ldpc.NumberRowsHcols = HRA_112_112_NUMBERROWSHCOLS;
ldpc.max_row_weight = HRA_112_112_MAX_ROW_WEIGHT;
ldpc.max_col_weight = HRA_112_112_MAX_COL_WEIGHT;
ldpc.H_rows = (uint16_t *)HRA_112_112_H_rows;
ldpc.H_cols = (uint16_t *)HRA_112_112_H_cols;
/* Main Loop
* ---------------------------------------------------------------------*/
for (f = 0; f < NFRAMES; f++) {
/* --------------------------------------------------------*\
Mod
\*---------------------------------------------------------*/
/* See CML startup code in tofdm.m */
for (i = 0; i < ofdm_nuwbits; i++) {
tx_bits[i] = ofdm->tx_uw[i];
}
for (i = ofdm_nuwbits; i < ofdm_nuwbits + ofdm_ntxtbits; i++) {
tx_bits[i] = 0;
}
if (ldpc_enable) {
unsigned char ibits[HRA_112_112_NUMBERROWSHCOLS];
unsigned char pbits[HRA_112_112_NUMBERPARITYBITS];
assert(HRA_112_112_NUMBERROWSHCOLS == ldpc.CodeLength / 2);
for (i = 0; i < ldpc.CodeLength / 2; i++) {
ibits[i] = payload_data_bits[i];
}
encode(&ldpc, ibits, pbits);
for (j = 0, i = ofdm_nuwbits + ofdm_ntxtbits; j < ldpc.CodeLength / 2;
i++, j++) {
tx_bits[i] = ibits[j];
}
for (j = 0; j < ldpc.CodeLength / 2; i++, j++) {
tx_bits[i] = pbits[j];
}
assert(i == ofdm_bitsperframe);
} else {
int Npayloadbits = ofdm_bitsperframe - (ofdm_nuwbits + ofdm_ntxtbits);
uint16_t r[Npayloadbits];
uint8_t payload_bits[Npayloadbits];
ofdm_rand(r, Npayloadbits);
for (i = 0; i < Npayloadbits; i++) payload_bits[i] = r[i] > 16384;
uint8_t txt_bits[ofdm_ntxtbits];
for (i = 0; i < ofdm_ntxtbits; i++) txt_bits[i] = 0;
uint8_t tx_bits_char[ofdm_bitsperframe];
ofdm_assemble_qpsk_modem_packet(ofdm, tx_bits_char, payload_bits,
txt_bits);
for (i = 0; i < ofdm_bitsperframe; i++) tx_bits[i] = tx_bits_char[i];
}
ofdm_mod(ofdm, (COMP *)tx, tx_bits);
/* tx vector logging */
memcpy(&tx_bits_log[ofdm_bitsperframe * f], tx_bits,
sizeof(int) * ofdm_bitsperframe);
memcpy(&tx_log[ofdm_samplesperframe * f], tx,
sizeof(COMP) * ofdm_samplesperframe);
}
/* --------------------------------------------------------*\
Channel
\*---------------------------------------------------------*/
fs_offset(rx_log, tx_log, ofdm_samplesperframe * NFRAMES,
SAMPLE_CLOCK_OFFSET_PPM);
COMP foff_phase_rect = {1.0f, 0.0f};
freq_shift(rx_log, rx_log, FOFF_HZ, &foff_phase_rect,
ofdm_samplesperframe * NFRAMES);
/* --------------------------------------------------------*\
Demod
\*---------------------------------------------------------*/
/* Init/pre-load rx with ideal timing so we can test with timing estimation
* disabled */
int Nsam = ofdm_samplesperframe * NFRAMES;
int prx = 0;
int nin = ofdm_samplesperframe + 2 * ofdm_samplespersymbol;
int lnew;
COMP rxbuf_in[ofdm_max_samplesperframe];
#define FRONT_LOAD
#ifdef FRONT_LOAD
for (i = 0; i < nin; i++, prx++) {
ofdm->rxbuf[ofdm_nrxbuf - nin + i] =
rx_log[prx].real + rx_log[prx].imag * I;
}
#endif
int nin_tot = 0;
/* disable estimators for initial testing */
ofdm_set_verbose(ofdm, false);
ofdm_set_timing_enable(ofdm, true);
ofdm_set_foff_est_enable(ofdm, true);
ofdm_set_phase_est_enable(ofdm, true);
//#define TESTING_FILE
#ifdef TESTING_FILE
FILE *fin = fopen("~/codec2-dev/octave/ofdm_test.raw", "rb");
assert(fin != NULL);
int Nbitsperframe = ofdm_bitsperframe;
int Nmaxsamperframe = ofdm_max_samplesperframe;
short rx_scaled[Nmaxsamperframe];
#endif
/* start this with something sensible otherwise LDPC decode fails in tofdm.m
*/
ofdm->mean_amp = 1.0;
for (f = 0; f < NFRAMES; f++) {
/* For initial testing, timing est is off, so nin is always
fixed. TODO: we need a constant for rxbuf_in[] size that
is the maximum possible nin */
nin = ofdm_get_nin(ofdm);
assert(nin <= ofdm_max_samplesperframe);
/* Insert samples at end of buffer, set to zero if no samples
available to disable phase estimation on future pilots on
last frame of simulation. */
if ((Nsam - prx) < nin) {
lnew = Nsam - prx;
} else {
lnew = nin;
}
// printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew);
for (i = 0; i < nin; i++) {
rxbuf_in[i].real = 0.0;
rxbuf_in[i].imag = 0.0;
}
if (lnew) {
for (i = 0; i < lnew; i++, prx++) {
rxbuf_in[i] = rx_log[prx];
}
}
assert(prx <= ofdm_max_samplesperframe * NFRAMES);
#ifdef TESTING_FILE
fread(rx_scaled, sizeof(short), nin, fin);
for (i = 0; i < nin; i++) {
rxbuf_in[i].real = (float)rx_scaled[i] / ASCALE;
rxbuf_in[i].imag = 0.0;
}
#endif
/* uncoded OFDM modem ---------------------------------------*/
ofdm_demod(ofdm, rx_bits, rxbuf_in);
#ifdef TESTING_FILE
int Nerrs = 0;
for (i = 0; i < Nbitsperframe; i++) {
if (test_bits_ofdm[i] != rx_bits[i]) {
Nerrs++;
}
}
printf("f: %d Nerr: %d\n", f, Nerrs);
#endif
float symbol_likelihood[(CODED_BITSPERFRAME / ofdm_bps) * (1 << ofdm_bps)];
float bit_likelihood[CODED_BITSPERFRAME];
uint8_t out_char[CODED_BITSPERFRAME];
if (ldpc_enable) {
/* LDPC functions --------------------------------------*/
float EsNo = 10;
/* first few symbols are used for UW and txt bits, find start of (224,112)
* LDPC codeword */
assert((ofdm_nuwbits + ofdm_ntxtbits + CODED_BITSPERFRAME) ==
ofdm_bitsperframe);
COMP ldpc_codeword_symbols[(CODED_BITSPERFRAME / ofdm_bps)];
for (i = 0, j = (ofdm_nuwbits + ofdm_ntxtbits) / ofdm_bps;
i < (CODED_BITSPERFRAME / ofdm_bps); i++, j++) {
ldpc_codeword_symbols[i].real = crealf(ofdm->rx_np[j]);
ldpc_codeword_symbols[i].imag = cimagf(ofdm->rx_np[j]);
}
float *ldpc_codeword_symbol_amps =
&ofdm->rx_amp[(ofdm_nuwbits + ofdm_ntxtbits) / ofdm_bps];
Demod2D(symbol_likelihood, ldpc_codeword_symbols, S_matrix, EsNo,
ldpc_codeword_symbol_amps, ofdm->mean_amp,
CODED_BITSPERFRAME / ofdm_bps);
Somap(bit_likelihood, symbol_likelihood, 1 << ofdm_bps, ofdm_bps,
CODED_BITSPERFRAME / ofdm_bps);
float llr[CODED_BITSPERFRAME];
int parityCheckCount;
// fprintf(stderr, "\n");
for (i = 0; i < CODED_BITSPERFRAME; i++) {
llr[i] = -bit_likelihood[i];
// fprintf(stderr, "%f ", llr[i]);
}
// fprintf(stderr, "\n");
run_ldpc_decoder(&ldpc, out_char, llr, &parityCheckCount);
/*
fprintf(stderr, "iter: %d parityCheckCount: %d\n", iter,
parityCheckCount); for(i=0; i<CODED_BITSPERFRAME; i++) { fprintf(stderr,
"%d ", out_char[i]);
}
*/
}
/* rx vector logging -----------------------------------*/
assert(nin_tot < ofdm_samplesperframe * NFRAMES);
memcpy(&rxbuf_in_log[nin_tot], rxbuf_in, sizeof(COMP) * nin);
nin_tot += nin;
for (i = 0; i < ofdm_nrxbuf; i++) {
rxbuf_log[ofdm_nrxbuf * f + i].real = crealf(ofdm->rxbuf[i]);
rxbuf_log[ofdm_nrxbuf * f + i].imag = cimagf(ofdm->rxbuf[i]);
}
for (i = 0; i < (ofdm_ns + 3); i++) {
for (j = 0; j < (ofdm_nc + 2); j++) {
rx_sym_log[(ofdm_ns + 3) * f + i][j].real = crealf(ofdm->rx_sym[i][j]);
rx_sym_log[(ofdm_ns + 3) * f + i][j].imag = cimagf(ofdm->rx_sym[i][j]);
}
}
/* note corrected phase (rx no phase) is one big linear array for frame */
for (i = 0; i < ofdm_rowsperframe * ofdm_nc; i++) {
rx_np_log[ofdm_rowsperframe * ofdm_nc * f + i].real =
crealf(ofdm->rx_np[i]);
rx_np_log[ofdm_rowsperframe * ofdm_nc * f + i].imag =
cimagf(ofdm->rx_np[i]);
}
/* note phase/amp ests the same for each col, but check them all anyway */
for (i = 0; i < ofdm_rowsperframe; i++) {
for (j = 0; j < ofdm_nc; j++) {
phase_est_pilot_log[ofdm_rowsperframe * f + i][j] =
ofdm->aphase_est_pilot_log[ofdm_nc * i + j];
rx_amp_log[ofdm_rowsperframe * ofdm_nc * f + ofdm_nc * i + j] =
ofdm->rx_amp[ofdm_nc * i + j];
}
}
foff_hz_log[f] = ofdm->foff_est_hz;
timing_est_log[f] = ofdm->timing_est + 1; /* offset by 1 to match Octave */
timing_valid_log[f] = ofdm->timing_valid;
timing_mx_log[f] = ofdm->timing_mx;
coarse_foff_est_hz_log[f] = ofdm->coarse_foff_est_hz;
sample_point_log[f] =
ofdm->sample_point + 1; /* offset by 1 to match Octave */
float EsNodB = ofdm_esno_est_calc(ofdm->rx_np, ofdm_rowsperframe * ofdm_nc);
snr_log[f] = ofdm_snr_from_esno(ofdm, EsNodB);
mean_amp_log[f] = ofdm->mean_amp;
memcpy(&rx_bits_log[ofdm_bitsperframe * f], rx_bits, sizeof(rx_bits));
if (ldpc_enable) {
for (i = 0; i < (CODED_BITSPERFRAME / ofdm_bps) * (1 << ofdm_bps); i++) {
symbol_likelihood_log[(CODED_BITSPERFRAME / ofdm_bps) *
(1 << ofdm_bps) * f +
i] = symbol_likelihood[i];
}
for (i = 0; i < CODED_BITSPERFRAME; i++) {
bit_likelihood_log[CODED_BITSPERFRAME * f + i] = bit_likelihood[i];
detected_data_log[CODED_BITSPERFRAME * f + i] = out_char[i];
}
}
}
/*---------------------------------------------------------*\
Dump logs to Octave file for evaluation
by tofdm.m Octave script
\*---------------------------------------------------------*/
fout = fopen("tofdm_out.txt", "wt");
assert(fout != NULL);
fprintf(fout, "# Created by tofdm.c\n");
octave_save_complex(fout, "pilot_samples_c", (COMP *)ofdm->pilot_samples, 1,
ofdm_samplespersymbol, ofdm_samplespersymbol);
octave_save_int(fout, "tx_bits_log_c", tx_bits_log, 1,
ofdm_bitsperframe * NFRAMES);
octave_save_complex(fout, "tx_log_c", (COMP *)tx_log, 1,
ofdm_samplesperframe * NFRAMES,
ofdm_samplesperframe * NFRAMES);
octave_save_complex(fout, "rx_log_c", (COMP *)rx_log, 1,
ofdm_samplesperframe * NFRAMES,
ofdm_samplesperframe * NFRAMES);
octave_save_complex(fout, "rxbuf_in_log_c", (COMP *)rxbuf_in_log, 1, nin_tot,
nin_tot);
octave_save_complex(fout, "rxbuf_log_c", (COMP *)rxbuf_log, 1,
ofdm_nrxbuf * NFRAMES, ofdm_nrxbuf * NFRAMES);
octave_save_complex(fout, "rx_sym_log_c", (COMP *)rx_sym_log,
(ofdm_ns + 3) * NFRAMES, ofdm_nc + 2, ofdm_nc + 2);
octave_save_float(fout, "phase_est_pilot_log_c", (float *)phase_est_pilot_log,
ofdm_rowsperframe * NFRAMES, ofdm_nc, ofdm_nc);
octave_save_float(fout, "rx_amp_log_c", (float *)rx_amp_log, 1,
ofdm_rowsperframe * ofdm_nc * NFRAMES,
ofdm_rowsperframe * ofdm_nc * NFRAMES);
octave_save_float(fout, "foff_hz_log_c", foff_hz_log, NFRAMES, 1, 1);
octave_save_int(fout, "timing_est_log_c", timing_est_log, NFRAMES, 1);
octave_save_int(fout, "timing_valid_log_c", timing_valid_log, NFRAMES, 1);
octave_save_float(fout, "timing_mx_log_c", timing_mx_log, NFRAMES, 1, 1);
octave_save_float(fout, "coarse_foff_est_hz_log_c", coarse_foff_est_hz_log,
NFRAMES, 1, 1);
octave_save_int(fout, "sample_point_log_c", sample_point_log, NFRAMES, 1);
octave_save_complex(fout, "rx_np_log_c", (COMP *)rx_np_log, 1,
ofdm_rowsperframe * ofdm_nc * NFRAMES,
ofdm_rowsperframe * ofdm_nc * NFRAMES);
octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1,
ofdm_bitsperframe * NFRAMES);
octave_save_float(fout, "symbol_likelihood_log_c", symbol_likelihood_log,
(CODED_BITSPERFRAME / ofdm_bps) * (1 << ofdm_bps) * NFRAMES,
1, 1);
octave_save_float(fout, "bit_likelihood_log_c", bit_likelihood_log,
CODED_BITSPERFRAME * NFRAMES, 1, 1);
octave_save_int(fout, "detected_data_log_c", detected_data_log, 1,
CODED_BITSPERFRAME * NFRAMES);
octave_save_float(fout, "snr_log_c", snr_log, NFRAMES, 1, 1);
octave_save_float(fout, "mean_amp_log_c", mean_amp_log, NFRAMES, 1, 1);
fclose(fout);
#ifdef TESTING_FILE
fclose(fin);
#endif
ofdm_destroy(ofdm);
return 0;
}
+94
View File
@@ -0,0 +1,94 @@
/*---------------------------------------------------------------------------*\
FILE........: tofdm_acq.c
AUTHORS.....: David Rowe
DATE CREATED: Mar 2021
Tests for the acquistion (sync) parts of the C version of the OFDM modem.
This program outputs a file of Octave vectors that are loaded and
automatically tested against the Octave version of the modem by the Octave
script tofdm_acq.m
\*---------------------------------------------------------------------------*/
#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "codec2_ofdm.h"
#include "octave.h"
#include "ofdm_internal.h"
#define MAX_FRAMES 500
int main(int argc, char *argv[]) {
struct OFDM *ofdm;
struct OFDM_CONFIG ofdm_config;
ofdm_init_mode("datac0", &ofdm_config);
ofdm = ofdm_create(&ofdm_config);
ofdm->data_mode = "burst";
ofdm->verbose = 2;
ofdm->timing_mx_thresh = 0.15;
ofdm->postambledetectoren = 1;
assert(ofdm != NULL);
int nin = ofdm_get_nin(ofdm);
int rxbufst = ofdm->rxbufst;
FILE *fin = fopen(argv[1], "rb");
assert(fin != NULL);
short rx_scaled[ofdm_get_max_samples_per_frame(ofdm)];
int f = 0;
float timing_mx_log[MAX_FRAMES];
int ct_est_log[MAX_FRAMES];
float foff_est_log[MAX_FRAMES];
int timing_valid_log[MAX_FRAMES];
int nin_log[MAX_FRAMES];
while (fread(rx_scaled, sizeof(short), nin, fin) == nin) {
fprintf(stderr, "%3d ", f);
ofdm_sync_search_shorts(ofdm, rx_scaled, ofdm->amp_scale / 2.0f);
if (f < MAX_FRAMES) {
timing_mx_log[f] = ofdm->timing_mx;
ct_est_log[f] = ofdm->ct_est;
foff_est_log[f] = ofdm->foff_est_hz;
timing_valid_log[f] = ofdm->timing_valid;
nin_log[f] = ofdm->nin;
}
f++;
// reset these to defaults, as they get modified when timing_valid asserted
ofdm->nin = nin;
ofdm->rxbufst = rxbufst;
}
fclose(fin);
/*---------------------------------------------------------*\
Dump logs to Octave file for evaluation
by tofdm_acq.m Octave script
\*---------------------------------------------------------*/
FILE *fout = fopen("tofdm_acq_out.txt", "wt");
assert(fout != NULL);
fprintf(fout, "# Created by tofdm_acq.c\n");
octave_save_complex(fout, "tx_preamble_c", (COMP *)ofdm->tx_preamble, 1,
ofdm->samplesperframe, ofdm->samplesperframe);
octave_save_complex(fout, "tx_postamble_c", (COMP *)ofdm->tx_postamble, 1,
ofdm->samplesperframe, ofdm->samplesperframe);
octave_save_float(fout, "timing_mx_log_c", timing_mx_log, 1, f, f);
octave_save_float(fout, "foff_est_log_c", foff_est_log, 1, f, f);
octave_save_int(fout, "ct_est_log_c", ct_est_log, 1, f);
octave_save_int(fout, "timing_valid_log_c", timing_valid_log, 1, f);
octave_save_int(fout, "nin_log_c", nin_log, 1, f);
fclose(fout);
ofdm_destroy(ofdm);
return 0;
}
+35
View File
@@ -0,0 +1,35 @@
/*---------------------------------------------------------------------------*\
FILE........: tqam16.c
AUTHOR......: David Rowe
DATE CREATED: August 2020
Simple sanity check for QAM16 symbol mapping.
\*---------------------------------------------------------------------------*/
#include <stdio.h>
#include <string.h>
#include "ofdm_internal.h"
int main(void) {
int c;
for (c = 0; c < 16; c++) {
int tx_bits[4], rx_bits[4];
for (int i = 0; i < 4; i++) tx_bits[i] = (c >> (3 - i)) & 0x1;
complex float symbol = qam16_mod(tx_bits);
qam16_demod(symbol, rx_bits);
if (memcmp(tx_bits, rx_bits, 4)) {
fprintf(stderr, "FAIL on %d!\ntx_bits: ", c);
for (int i = 0; i < 4; i++) fprintf(stderr, "%d ", tx_bits[i]);
fprintf(stderr, "%f %f\nrx_bits: ", creal(symbol), cimag(symbol));
for (int i = 0; i < 4; i++) fprintf(stderr, "%d ", rx_bits[i]);
fprintf(stderr, "%f %f\n", creal(symbol), cimag(symbol));
return 1;
}
}
fprintf(stderr, "%d tested OK...\nPASS!\n", c);
return 0;
}
+46
View File
@@ -0,0 +1,46 @@
/*
tquisk_filter.c
Unit test for complex band pass filters in src/filter.c
cd codec2/build_linux
./misc/mksine - 1500 2 | unittest/tquisk_filter | aplay
By adjusting the frequency you can audibly test filter response.
*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "filter.h"
#include "filter_coef.h"
#define N 159 /* processing buffer size (odd number deliberate) */
#define CENTRE 1500.0
#define FS 8000.0
int main() {
short buf_short[N];
complex float buf[N];
struct quisk_cfFilter *bpf;
int i;
bpf = malloc(sizeof(struct quisk_cfFilter));
assert(bpf != NULL);
quisk_filt_cfInit(bpf, filtP200S400, sizeof(filtP200S400) / sizeof(float));
quisk_cfTune(bpf, CENTRE / FS);
while (fread(buf_short, sizeof(short), N, stdin) == N) {
for (i = 0; i < N; i++) buf[i] = buf_short[i];
quisk_ccfFilter(buf, buf, N, bpf);
/* we only output the real part in this test */
for (i = 0; i < N; i++) buf_short[i] = creal(buf[i]);
fwrite(buf_short, sizeof(short), N, stdout);
}
quisk_filt_destroy(bpf);
free(bpf);
return 0;
}
+32
View File
@@ -0,0 +1,32 @@
/*
tvq_mbest.c
David Rowe Dec 2019
Generate some test vectors to exercise misc/vq_mbest.c
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
void write_float_file(char fn[], float *values, int n) {
FILE *f = fopen(fn, "wb");
assert(f != NULL);
assert(fwrite(values, sizeof(float), n, f) == n);
fclose(f);
}
int main(void) {
/* we're only interested in searching the inner 2 values, outer elements
should be ignored */
float target[] = {0.0, 1.0, 1.0, 0.0};
write_float_file("target.f32", target, 4);
float vq1[] = {
1.0, 0.9, 0.9, 1.0, /* this will be a better match on first stage */
2.0, 0.8, 0.8, 2.0}; /* but after second stage should choose this */
write_float_file("vq1.f32", vq1, 8);
float vq2[] = {10.0, 0.3, 0.3, 10.0, 20.0,
0.2, 0.2, 20.0}; /* 0.8+0.2 == 1.0 so best 2nd stage entry */
write_float_file("vq2.f32", vq2, 8);
return 0;
}
+300
View File
@@ -0,0 +1,300 @@
/*
vq_mbest.c
David Rowe Dec 2019
Utility to perform a mbest VQ search on vectors from stdin, sending
quantised vectors to stdout.
*/
#include <assert.h>
#include <getopt.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mbest.h"
#define MAX_K 20
#define MAX_ENTRIES 4096
#define MAX_STAGES 5
void quant_mbest(float vec_out[], int indexes[], float vec_in[], int num_stages,
float vqw[], float vq[], int m[], int k, int mbest_survivors);
int verbose = 0;
int main(int argc, char *argv[]) {
float vq[MAX_STAGES * MAX_K * MAX_ENTRIES];
float vqw[MAX_STAGES * MAX_K * MAX_ENTRIES];
int m[MAX_STAGES];
int k = 0, mbest_survivors = 1, num_stages = 0;
char fnames[256], fn[256], *comma, *p;
FILE *fq;
float lower = -1E32;
int st = -1;
int en = -1;
int num = INT_MAX;
int output_vec_usage = 0;
int o = 0;
int opt_idx = 0;
while (o != -1) {
static struct option long_opts[] = {{"k", required_argument, 0, 'k'},
{"quant", required_argument, 0, 'q'},
{"mbest", required_argument, 0, 'm'},
{"lower", required_argument, 0, 'l'},
{"verbose", required_argument, 0, 'v'},
{"st", required_argument, 0, 't'},
{"en", required_argument, 0, 'e'},
{"num", required_argument, 0, 'n'},
{"vec_usage", no_argument, 0, 'u'},
{0, 0, 0, 0}};
o = getopt_long(argc, argv, "hk:q:m:vt:e:n:u", long_opts, &opt_idx);
switch (o) {
case 'k':
k = atoi(optarg);
assert(k <= MAX_K);
break;
case 'q':
/* load up list of comma delimited file names */
strcpy(fnames, optarg);
p = fnames;
num_stages = 0;
do {
assert(num_stages < MAX_STAGES);
strcpy(fn, p);
comma = strchr(fn, ',');
if (comma) {
*comma = 0;
p = comma + 1;
}
/* load quantiser file */
fprintf(stderr, "stage: %d loading %s ... ", num_stages, fn);
fq = fopen(fn, "rb");
if (fq == NULL) {
fprintf(stderr, "Couldn't open: %s\n", fn);
exit(1);
}
/* count how many entries m of dimension k are in this VQ file */
m[num_stages] = 0;
float dummy[k];
while (fread(dummy, sizeof(float), k, fq) == (size_t)k)
m[num_stages]++;
assert(m[num_stages] <= MAX_ENTRIES);
fprintf(stderr, "%d entries of vectors width %d\n", m[num_stages], k);
/* now load VQ into memory */
rewind(fq);
int rd = fread(&vq[num_stages * k * MAX_ENTRIES], sizeof(float),
m[num_stages] * k, fq);
assert(rd == m[num_stages] * k);
num_stages++;
fclose(fq);
} while (comma);
break;
case 'm':
mbest_survivors = atoi(optarg);
fprintf(stderr, "mbest_survivors = %d\n", mbest_survivors);
break;
case 'n':
num = atoi(optarg);
break;
case 'l':
lower = atof(optarg);
break;
case 't':
st = atoi(optarg);
break;
case 'e':
en = atoi(optarg);
break;
case 'u':
output_vec_usage = 1;
break;
case 'v':
verbose = 1;
break;
help:
fprintf(stderr, "\n");
fprintf(stderr,
"usage: %s -k dimension -q vq1.f32,vq2.f32,.... [Options]\n",
argv[0]);
fprintf(stderr, "\n");
fprintf(stderr,
"input vectors on stdin, output quantised vectors on stdout\n");
fprintf(stderr, "\n");
fprintf(stderr,
"--lower lowermeanLimit Only count vectors with average "
"above this level in distortion calculations\n");
fprintf(stderr,
"--mbest N number of survivors at each stage, "
"set to 0 for standard VQ search\n");
fprintf(stderr,
"--st Kst start vector element for error "
"calculation (default 0)\n");
fprintf(stderr,
"--en Ken end vector element for error "
"calculation (default K-1)\n");
fprintf(stderr,
"--num numToProcess number of vectors to quantise "
"(default to EOF)\n");
fprintf(stderr,
"--vec_usage Output a record of how many times "
"each vector is used\n");
fprintf(stderr, "-v Verbose\n");
exit(1);
}
}
if ((num_stages == 0) || (k == 0)) goto help;
/* default to measuring error on entire vector */
if (st == -1) st = 0;
if (en == -1) en = k - 1;
float w[k];
for (int i = 0; i < st; i++) w[i] = 0.0;
for (int i = st; i <= en; i++) w[i] = 1.0;
for (int i = en + 1; i < k; i++) w[i] = 0.0;
/* apply weighting to codebook (rather than in search) */
memcpy(vqw, vq, sizeof(vq));
for (int s = 0; s < num_stages; s++) {
mbest_precompute_weight(&vqw[s * k * MAX_ENTRIES], w, k, m[s]);
}
int indexes[num_stages], nvecs = 0;
int vec_usage[m[0]];
for (int i = 0; i < m[0]; i++) vec_usage[i] = 0;
float target[k], quantised[k];
float sqe = 0.0;
while (fread(&target, sizeof(float), k, stdin) && (nvecs < num)) {
for (int i = 0; i < k; i++) target[i] *= w[i];
int dont_count = 0;
/* optional clamping to lower limit or mean */
float mean = 0.0;
for (int i = 0; i < k; i++) mean += target[i];
mean /= k;
float difference = mean - lower;
if (difference < 0.0) {
/* bring target up to lower clamping limit */
for (int i = 0; i < k; i++) target[i] += -difference;
dont_count = 1;
}
quant_mbest(quantised, indexes, target, num_stages, vqw, vq, m, k,
mbest_survivors);
if (dont_count == 0) {
for (int i = st; i <= en; i++) sqe += pow(target[i] - quantised[i], 2.0);
}
fwrite(&quantised, sizeof(float), k, stdout);
nvecs++;
// count number f time each vector is used (just for first stage)
vec_usage[indexes[0]]++;
}
fprintf(stderr, "MSE: %4.2f\n", sqe / (nvecs * (en - st + 1)));
if (output_vec_usage) {
for (int i = 0; i < m[0]; i++) fprintf(stderr, "%d\n", vec_usage[i]);
}
return 0;
}
// print vector debug function
void pv(char s[], float v[], int k) {
int i;
if (verbose) {
fprintf(stderr, "%s", s);
for (i = 0; i < k; i++) fprintf(stderr, "%4.2f ", v[i]);
fprintf(stderr, "\n");
}
}
// mbest algorithm version, backported from LPCNet/src
void quant_mbest(float vec_out[], int indexes[], float vec_in[], int num_stages,
float vqw[], float vq[], int m[], int k, int mbest_survivors) {
float err[k], se1;
int i, j, s, s1, ind;
struct MBEST *mbest_stage[num_stages];
int index[num_stages];
float target[k];
for (i = 0; i < num_stages; i++) {
mbest_stage[i] = mbest_create(mbest_survivors);
index[i] = 0;
}
se1 = 0.0;
for (i = 0; i < k; i++) {
err[i] = vec_in[i];
se1 += err[i] * err[i];
}
se1 /= k;
/* now quantise err[] using multi-stage mbest search, preserving
mbest_survivors at each stage */
mbest_search(vqw, err, k, m[0], mbest_stage[0], index);
if (verbose) mbest_print("Stage 1:", mbest_stage[0]);
for (s = 1; s < num_stages; s++) {
/* for each candidate in previous stage, try to find best vector in next
* stage */
for (j = 0; j < mbest_survivors; j++) {
/* indexes that lead us this far */
for (s1 = 0; s1 < s; s1++) {
index[s1 + 1] = mbest_stage[s - 1]->list[j].index[s1];
}
/* target is residual err[] vector given path to this candidate */
for (i = 0; i < k; i++) target[i] = err[i];
for (s1 = 0; s1 < s; s1++) {
ind = index[s - s1];
if (verbose)
fprintf(stderr, " s: %d s1: %d s-s1: %d ind: %d\n", s, s1, s - s1,
ind);
for (i = 0; i < k; i++) {
target[i] -= vqw[s1 * k * MAX_ENTRIES + ind * k + i];
}
}
pv(" target: ", target, k);
mbest_search(&vqw[s * k * MAX_ENTRIES], target, k, m[s], mbest_stage[s],
index);
}
char str[80];
sprintf(str, "Stage %d:", s + 1);
if (verbose) mbest_print(str, mbest_stage[s]);
}
for (s = 0; s < num_stages; s++) {
indexes[s] = mbest_stage[num_stages - 1]->list[0].index[num_stages - 1 - s];
}
/* OK put it all back together using best survivor */
for (i = 0; i < k; i++) vec_out[i] = 0.0;
for (s = 0; s < num_stages; s++) {
int ind = indexes[s];
float se2 = 0.0;
for (i = 0; i < k; i++) {
err[i] -= vqw[s * k * MAX_ENTRIES + ind * k + i];
vec_out[i] += vq[s * k * MAX_ENTRIES + ind * k + i];
se2 += err[i] * err[i];
}
se2 /= k;
pv(" err: ", err, k);
if (verbose) fprintf(stderr, " se2: %f\n", se2);
}
pv(" vec_out: ", vec_out, k);
pv("\n vec_in: ", vec_in, k);
pv(" vec_out: ", vec_out, k);
pv(" err: ", err, k);
if (verbose) fprintf(stderr, " se1: %f\n", se1);
for (i = 0; i < num_stages; i++) mbest_destroy(mbest_stage[i]);
}