Compare commits

..

3 commits

Author SHA1 Message Date
Your Name
91f2e648c8 Add library build and install 2021-01-13 09:44:19 +00:00
Karlis Goba
c6d79b961c Merge branch 'master' of https://github.com/kgoba/ft8_lib 2019-11-15 10:23:06 +02:00
Karlis Goba
2d1337d47e Fixed FT8 waveform generation for high frequencies 2019-11-15 10:22:45 +02:00
11 changed files with 117 additions and 1509 deletions

View file

@ -21,3 +21,6 @@ decode_ft8: decode_ft8.o fft/kiss_fftr.o fft/kiss_fft.o ft8/decode.o ft8/encode.
clean:
rm -f *.o ft8/*.o common/*.o fft/*.o $(TARGETS)
install:
$(AR) rc libft8.a ft8/constants.o ft8/encode.o ft8/pack.o ft8/text.o common/wave.o
install libft8.a /usr/lib/libft8.a

View file

@ -15,9 +15,6 @@
#define LOG_LEVEL LOG_INFO
const int kFreq_osr = 2;
const int kTime_osr = 2;
const int kMin_score = 40; // Minimum sync score threshold for candidates
const int kMax_candidates = 120;
const int kLDPC_iterations = 25;
@ -25,6 +22,11 @@ const int kLDPC_iterations = 25;
const int kMax_decoded_messages = 50;
const int kMax_message_length = 25;
const int kFreq_osr = 2;
const int kTime_osr = 2;
const float kFSK_dev = 6.25f; // tone deviation in Hz and symbol rate
void usage() {
fprintf(stderr, "Decode a 15-second WAV file.\n");
}
@ -58,63 +60,10 @@ float blackman_i(int i, int N) {
return a0 - a1*x1 + a2*x2;
}
class Monitor : public ft8::Monitor1Base {
public:
Monitor(float sample_rate);
protected:
virtual void fft_forward(const float *in, std::complex<float> *out) override;
private:
kiss_fftr_cfg fft_cfg;
};
Monitor::Monitor(float sample_rate) : Monitor1Base(sample_rate)
{
window_fn = new float[nfft]; // [nfft]
fft_frame = new float[nfft]; // [nfft]
last_frame = new float[nfft * 3/4]; // [nfft * 3/4]
freqdata = new std::complex<float>[nfft/2 + 1]; // [nfft/2 + 1]
fft_cfg = kiss_fftr_alloc(nfft, 0, NULL, NULL);
static float max2(float a, float b) {
return (a >= b) ? a : b;
}
void Monitor::fft_forward(const float *in, std::complex<float> *out) {
kiss_fftr(fft_cfg,
reinterpret_cast<const kiss_fft_scalar *>(in),
reinterpret_cast<kiss_fft_cpx *>(out));
}
// #include <complex>
// class FFT_r2c {
// public:
// // [N] real --> [N/2 + 1] complex
// virtual void forward(const float *in, std::complex<float> *out) = 0;
// };
// class KissFFT_r2c : public FFT_r2c {
// public:
// KissFFT_r2c(int N) {
// LOG(LOG_INFO, "N_FFT = %d\n", N);
// // size_t fft_work_size;
// // kiss_fftr_alloc(N, 0, 0, &fft_work_size);
// // LOG(LOG_INFO, "FFT work area = %lu\n", fft_work_size);
// // fft_work = malloc(fft_work_size);
// // fft_cfg = kiss_fftr_alloc(N, 0, fft_work, &fft_work_size);
// fft_cfg = kiss_fftr_alloc(N, 0, NULL, NULL);
// }
// ~KissFFT_r2c() {
// // free(fft_work);
// kiss_fftr_free(fft_cfg);
// }
// virtual void forward(const float *in, std::complex<float> *out) override {
// kiss_fftr(fft_cfg, in, reinterpret_cast<kiss_fft_cpx *>(out));
// }
// private:
// void *fft_work;
// kiss_fftr_cfg fft_cfg;
// };
// Compute FFT magnitudes (log power) for each timeslot in the signal
void extract_power(const float signal[], ft8::MagArray * power) {
const int block_size = 2 * power->num_bins; // Average over 2 bins per FSK tone
@ -157,16 +106,17 @@ void extract_power(const float signal[], ft8::MagArray * power) {
// Compute log magnitude in decibels
for (int j = 0; j < nfft/2 + 1; ++j) {
float mag2 = (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r);
mag_db[j] = 10.0f * log10f(1E-12f + mag2 * fft_norm * fft_norm);
mag_db[j] = 10.0f * log10f(1E-10f + mag2 * fft_norm * fft_norm);
}
// Loop over two possible frequency bin offsets (for averaging)
for (int freq_sub = 0; freq_sub < power->freq_osr; ++freq_sub) {
for (int j = 0; j < power->num_bins; ++j) {
//float db1 = mag_db[j * power->freq_osr + freq_sub];
float db1 = mag_db[j * power->freq_osr + freq_sub];
//float db2 = mag_db[j * 2 + freq_sub + 1];
//float db = (db1 + db2) / 2;
float db = mag_db[j * power->freq_osr + freq_sub];
float db = db1;
//float db = sqrtf(db1 * db2);
// Scale decibels to unsigned 8-bit range and clamp the value
int scaled = (int)(2 * (db + 120));
@ -230,7 +180,7 @@ int main(int argc, char **argv) {
normalize_signal(signal, num_samples);
// Compute DSP parameters that depend on the sample rate
const int num_bins = (int)(sample_rate / (2 * ft8::FSK_dev));
const int num_bins = (int)(sample_rate / (2 * kFSK_dev));
const int block_size = 2 * num_bins;
const int subblock_size = block_size / kTime_osr;
const int nfft = block_size * kFreq_osr;
@ -262,8 +212,8 @@ int main(int argc, char **argv) {
ft8::Candidate &cand = candidate_list[idx];
if (cand.score < kMin_score) continue;
float freq_hz = (cand.freq_offset + (float)cand.freq_sub / kFreq_osr) * ft8::FSK_dev;
float time_sec = (cand.time_offset + (float)cand.time_sub / kTime_osr) / ft8::FSK_dev;
float freq_hz = (cand.freq_offset + (float)cand.freq_sub / kFreq_osr) * kFSK_dev;
float time_sec = (cand.time_offset + (float)cand.time_sub / kTime_osr) / kFSK_dev;
float log174[ft8::N];
ft8::extract_likelihood(&power, cand, ft8::kGray_map, log174);
@ -272,9 +222,7 @@ int main(int argc, char **argv) {
uint8_t plain[ft8::N];
int n_errors = 0;
ft8::bp_decode(log174, kLDPC_iterations, plain, &n_errors);
if (n_errors > 0) {
// ft8::ldpc_decode(log174, kLDPC_iterations, plain, &n_errors);
}
//ft8::ldpc_decode(log174, kLDPC_iterations, plain, &n_errors);
if (n_errors > 0) {
LOG(LOG_DEBUG, "ldpc_decode() = %d (%.0f Hz)\n", n_errors, freq_hz);

View file

@ -3,7 +3,6 @@
#include <stdint.h>
namespace ft8 {
constexpr float FSK_dev = 6.25f; // tone deviation in Hz and symbol rate
// Define FT8 symbol counts
constexpr int ND = 58; // Data symbols

View file

@ -17,57 +17,6 @@ static int get_index(const MagArray *power, int block, int time_sub, int freq_su
return ((((block * power->time_osr) + time_sub) * power->freq_osr + freq_sub) * power->num_bins) + bin;
}
Monitor1Base::Monitor1Base(float sample_rate, int time_osr, int freq_osr, float fmin, float fmax) {
int block_size = (int)(0.5f + sample_rate / ft8::FSK_dev); // Samples per FSK tone
nfft = block_size * freq_osr; // FFT over symbols with frequency oversampling
int bin1 = (block_size * fmin) / sample_rate;
int bin2 = (block_size * fmax) / sample_rate;
power.time_osr = time_osr;
power.freq_osr = freq_osr;
power.num_bins = bin2 - bin1;
power.num_blocks = 0;
}
void Monitor1Base::feed(const float *frame) {
// Fill the first 3/4 of analysis frame
for (int i = 0; i < 3 * nfft / 4; ++i) {
fft_frame[i] = window_fn[i] * last_frame[i];
}
// Shift the frame history
for (int i = 0; i < nfft / 2; ++i) {
last_frame[i] = last_frame[i + nfft / 4];
}
// Now fill the last_frame array
for (int i = 0; i < nfft / 4; ++i) {
last_frame[i + nfft / 2] = frame[i];
}
// Fill the last 1/4 of analysis frame
for (int i = 3 * nfft / 4, j = nfft / 2; i < nfft; ++i, ++j) {
fft_frame[i] = window_fn[i] * last_frame[j];
}
fft_forward(fft_frame, freqdata);
for (int freq_sub = 0; freq_sub < power.freq_osr; ++freq_sub) {
for (int i = 0; i < power.num_bins; i += power.freq_osr) {
float mag2 = std::norm(freqdata[i]); // re^2 + im^2
float mag_db = 10.0f * log10f(1E-12f + mag2);
int scaled = (int)(2 * (mag_db + 120));
power.mag[offset] = (scaled < 0) ? 0 : ((scaled > 255) ? 255 : scaled);
++offset;
}
}
if (++time_sub >= power.time_osr) {
time_sub = 0;
++power.num_blocks;
}
}
// Localize top N candidates in frequency and time according to their sync strength (looking at Costas symbols)
// We treat and organize the candidate list as a min-heap (empty initially).
int find_sync(const MagArray *power, const uint8_t *sync_map, int num_candidates, Candidate *heap, int min_score) {
@ -86,14 +35,13 @@ int find_sync(const MagArray *power, const uint8_t *sync_map, int num_candidates
// Compute average score over sync symbols (m+k = 0-7, 36-43, 72-79)
int num_symbols = 0;
for (int m = 0; m <= 72; m += 36) {
// Iterate over 7 Costas synchronisation symbols
for (int k = 0; k < 7; ++k) {
int n = time_offset + k + m;
// Check for time boundaries
if (n < 0) continue;
if (n >= power->num_blocks) break;
if (time_offset + k + m < 0) continue;
if (time_offset + k + m >= power->num_blocks) break;
int offset = get_index(power, n, time_sub, freq_sub, freq_offset);
// int offset = ((time_offset + k + m) * num_alt + alt) * power->num_bins + freq_offset;
int offset = get_index(power, time_offset + k + m, time_sub, freq_sub, freq_offset);
const uint8_t *p8 = power->mag + offset;
// Weighted difference between the expected and all other symbols
@ -112,18 +60,18 @@ int find_sync(const MagArray *power, const uint8_t *sync_map, int num_candidates
// look at one frequency bin higher
score += p8[sm] - p8[sm + 1];
}
if (k > 0 && n - 1 >= 0) {
if (k > 0) {
// look one symbol back in time
score += p8[sm] - p8[sm - num_alt * power->num_bins];
}
if (k < 6 && n + 1 < power->num_blocks) {
if (k < 6) {
// look one symbol forward in time
score += p8[sm] - p8[sm + num_alt * power->num_bins];
}
++num_symbols;
}
}
if (num_symbols > 0)
score /= num_symbols;
if (score < min_score) continue;
@ -173,11 +121,10 @@ void extract_likelihood(const MagArray *power, const Candidate & cand, const uin
int sym_idx = (k < ft8::ND / 2) ? (k + 7) : (k + 14);
int bit_idx = 3 * k;
// Index of the 8 bins of the current symbol
int sym_offset = offset + sym_idx * num_alt * power->num_bins;
// Pointer to 8 bins of the current symbol
const uint8_t *ps = power->mag + (offset + sym_idx * num_alt * power->num_bins);
decode_symbol(power->mag + sym_offset, code_map, bit_idx, log174);
// decode_multi_symbols(power->mag + sym_offset, power->num_bins, n_syms, code_map, bit_idx, log174);
decode_symbol(ps, code_map, bit_idx, log174);
}
// Compute the variance of log174

View file

@ -1,7 +1,6 @@
#pragma once
#include <stdint.h>
#include <complex>
namespace ft8 {
@ -21,28 +20,6 @@ struct Candidate {
uint8_t freq_sub;
};
class Monitor1Base {
public:
Monitor1Base(float sample_rate, int time_osr = 2, int freq_osr = 2, float fmin = 300, float fmax = 3000);
void feed(const float *frame);
void search();
void reset();
protected:
float *window_fn; // [nfft]
float *fft_frame; // [nfft]
float *last_frame; // [nfft * 3/4]
std::complex<float> *freqdata; // [nfft/2 + 1]
int nfft;
int offset;
int time_sub;
ft8::MagArray power;
// [N] real --> [N/2 + 1] log magnitudes (decibels)
// virtual void fft_forward_mag_db(const float *frame, uint8_t *mag_db) = 0;
virtual void fft_forward(const float *in, std::complex<float> *out) = 0;
};
// Localize top N candidates in frequency and time according to their sync strength (looking at Costas symbols)
// We treat and organize the candidate list as a min-heap (empty initially).

View file

@ -8,11 +8,11 @@
// from Sarah Johnson's Iterative Error Correction book.
// codeword[i] = log ( P(x=0) / P(x=1) )
//
#include "ldpc.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "constants.h"
namespace ft8 {
@ -47,13 +47,11 @@ void pack_bits(const uint8_t plain[], int num_bits, uint8_t packed[]) {
// codeword is 174 log-likelihoods.
// plain is a return value, 174 ints, to be 0 or 1.
// max_iters is how hard to try.
// n_errors == 87 means success.
void ldpc_decode(const float codeword[], int max_iters, uint8_t plain[], int *n_errors) {
// ok == 87 means success.
void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok) {
float m[ft8::M][ft8::N]; // ~60 kB
float e[ft8::M][ft8::N]; // ~60 kB
int min_errors = ft8::M;
int n_err_last = 0;
int n_cnt = 0;
for (int j = 0; j < ft8::M; j++) {
for (int i = 0; i < ft8::N; i++) {
@ -95,22 +93,6 @@ void ldpc_decode(const float codeword[], int max_iters, uint8_t plain[], int *n_
}
}
// Early stopping criterion
if (iter > 0) {
int nd = errors - n_err_last;
if (nd < 0) {
n_cnt = 0;
}
else {
++n_cnt;
}
if (n_cnt >= 5 && iter >= 10 && errors >= 15) {
*n_errors = errors;
return;
}
}
n_err_last = errors;
for (int i = 0; i < ft8::N; i++) {
for (int ji1 = 0; ji1 < 3; ji1++) {
int j1 = kMn[i][ji1] - 1;
@ -126,7 +108,7 @@ void ldpc_decode(const float codeword[], int max_iters, uint8_t plain[], int *n_
}
}
*n_errors = min_errors;
*ok = min_errors;
}
@ -151,19 +133,32 @@ static int ldpc_check(uint8_t codeword[]) {
}
void BPDecoderState::reset() {
void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok) {
float tov[ft8::N][3];
float toc[ft8::M][7];
int min_errors = ft8::M;
int nclast = 0;
int ncnt = 0;
// initialize messages to checks
for (int i = 0; i < ft8::M; ++i) {
for (int j = 0; j < kNrw[i]; ++j) {
toc[i][j] = codeword[kNm[i][j] - 1];
}
}
for (int i = 0; i < ft8::N; ++i) {
for (int j = 0; j < 3; ++j) {
tov[i][j] = 0;
}
}
}
int BPDecoderState::iterate(const float codeword[], uint8_t plain[]) {
// Update bit log likelihood ratios (tov=0 in iter 0)
for (int iter = 0; iter < max_iters; ++iter) {
float zn[ft8::N];
// Update bit log likelihood ratios (tov=0 in iter 0)
for (int i = 0; i < ft8::N; ++i) {
zn[i] = codeword[i] + tov[i][0] + tov[i][1] + tov[i][2];
plain[i] = (zn[i] > 0) ? 1 : 0;
@ -172,8 +167,13 @@ int BPDecoderState::iterate(const float codeword[], uint8_t plain[]) {
// Check to see if we have a codeword (check before we do any iter)
int errors = ldpc_check(plain);
if (errors < min_errors) {
// we have a better guess - update the result
min_errors = errors;
if (errors == 0) {
return 0; // Found a perfect answer
break; // Found a perfect answer
}
}
// Send messages from bits to check nodes
@ -209,44 +209,9 @@ int BPDecoderState::iterate(const float codeword[], uint8_t plain[]) {
tov[i][j] = 2 * fast_atanh(-Tmn);
}
}
return errors;
}
void bp_decode(const float codeword[], int max_iters, uint8_t plain[], int *n_errors) {
BPDecoderState state;
int n_err_last = 0;
int n_cnt = 0;
state.reset();
*n_errors = ft8::M;
for (int iter = 0; iter < max_iters; ++iter) {
int errors = state.iterate(codeword, plain);
if (errors < *n_errors) {
*n_errors = errors;
}
if (errors == 0) return;
// Early stopping criterion
if (iter > 0) {
int nd = errors - n_err_last;
if (nd < 0) {
n_cnt = 0;
}
else {
++n_cnt;
}
if (n_cnt >= 5 && iter >= 10 && errors >= 15) {
*n_errors = errors;
return;
}
}
n_err_last = errors;
}
*ok = min_errors;
}
// https://varietyofsound.wordpress.com/2011/02/14/efficient-tanh-computation-using-lamberts-continued-fraction/

View file

@ -1,26 +1,14 @@
#pragma once
#include "constants.h"
namespace ft8 {
class BPDecoderState {
public:
void reset();
int iterate(const float codeword[], uint8_t plain[]);
private:
float tov[ft8::N][3];
float toc[ft8::M][7];
};
// codeword is 174 log-likelihoods.
// plain is a return value, 174 ints, to be 0 or 1.
// iters is how hard to try.
// n_errors == 0 means success.
void ldpc_decode(const float codeword[], int max_iters, uint8_t plain[], int *n_errors);
// ok == 87 means success.
void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok);
void bp_decode(const float codeword[], int max_iters, uint8_t plain[], int *n_errors);
void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok);
// Packs a string of bits each represented as a zero/non-zero byte in plain[],
// as a string of packed bits starting from the MSB of the first byte of packed[]

View file

@ -1,49 +0,0 @@
#include "message77.h"
#include "unpack.h"
#include <string.h>
namespace ft8 {
Message77::Message77() {
i3 = n3 = 0;
field1[0] = field2[0] = field3[0] = '\0';
}
int Message77::str(char *buf, int buf_sz) const {
// Calculate the available space sans the '\0' terminator
int rem_sz = buf_sz - 1;
int field1_sz = strlen(field1);
int field2_sz = strlen(field2);
int field3_sz = strlen(field3);
int msg_sz = field1_sz + (field2_sz > 0 ? 1 : 0) +
field2_sz + (field3_sz > 0 ? 1 : 0) +
field3_sz;
if (rem_sz < msg_sz) return -1;
char *dst = buf;
dst = stpcpy(dst, field1);
*dst++ = ' ';
dst = stpcpy(dst, field2);
*dst++ = ' ';
dst = stpcpy(dst, field3);
*dst = '\0';
return msg_sz;
}
int Message77::unpack(const uint8_t *a77) {
// Extract n3 (bits 71..73) and i3 (bits 74..76)
n3 = ((a77[8] << 2) & 0x04) | ((a77[9] >> 6) & 0x03);
i3 = (a77[9] >> 3) & 0x07;
int rc = unpack77_fields(a77, field1, field2, field3);
if (rc < 0) {
field1[0] = field2[0] = field3[0] = '\0';
}
return rc;
}
} // namespace ft8

View file

@ -1,48 +0,0 @@
#pragma once
#include <stdint.h>
namespace ft8 {
class CallsignHasher {
virtual void save_callsign(void *obj, const char *callsign) = 0;
virtual bool hash10(void *obj, uint16_t hash, char *result) = 0;
virtual bool hash12(void *obj, uint16_t hash, char *result) = 0;
virtual bool hash22(void *obj, uint32_t hash, char *result) = 0;
};
class EmptyHasher : public CallsignHasher {
virtual void save_callsign(void *obj, const char *callsign) override {
}
virtual bool hash10(void *obj, uint16_t hash, char *result) override {
strcpy(result, "...");
return true;
}
virtual bool hash12(void *obj, uint16_t hash, char *result) override {
strcpy(result, "...");
return true;
}
virtual bool hash22(void *obj, uint32_t hash, char *result) override {
strcpy(result, "...");
return true;
}
};
struct Message77 {
uint8_t i3, n3;
// 11 chars nonstd call + 2 chars <...>
// 6 chars for grid/report/courtesy
char field1[13 + 1];
char field2[13 + 1];
char field3[6 + 1];
Message77();
int unpack(const uint8_t *packed77);
void pack(uint8_t *packed77);
int str(char *buf, int buf_sz) const;
private:
};
} // namespace ft8

View file

@ -53,11 +53,11 @@ int unpack28(uint32_t n28, uint8_t ip, uint8_t i3, char *result) {
// This is a 22-bit hash of a result
//call hash22(n22,c13) !Retrieve result from hash table
// TODO: implement
strcpy(result, "<...>");
// result[0] = '<';
// int_to_dd(result + 1, n28, 7);
// result[8] = '>';
// result[9] = '\0';
// strcpy(result, "<...>");
result[0] = '<';
int_to_dd(result + 1, n28, 7);
result[8] = '>';
result[9] = '\0';
return 0;
}
@ -181,6 +181,7 @@ int unpack_type1(const uint8_t *a77, uint8_t i3, char *field1, char *field2, cha
int unpack_text(const uint8_t *a71, char *text) {
// TODO: test
uint8_t b71[9];
uint8_t carry = 0;
@ -271,11 +272,11 @@ int unpack_nonstandard(const uint8_t *a77, char *field1, char *field2, char *fie
char call_3[15];
// should replace with hash12(n12, call_3);
strcpy(call_3, "<...>");
// call_3[0] = '<';
// int_to_dd(call_3 + 1, n12, 4);
// call_3[5] = '>';
// call_3[6] = '\0';
// strcpy(call_3, "<...>");
call_3[0] = '<';
int_to_dd(call_3 + 1, n12, 4);
call_3[5] = '>';
call_3[6] = '\0';
char * call_1 = (iflip) ? c11 : call_3;
char * call_2 = (iflip) ? call_3 : c11;
@ -369,4 +370,4 @@ int unpack77(const uint8_t *a77, char *message) {
return 0;
}
} // namespace ft8
} // namespace

File diff suppressed because it is too large Load diff