Updated decoding to v2

Makefile

@@ -1,8 +1,7 @@
 CXXFLAGS = -std=c++14 -I.
 LDFLAGS = -lm
 
-TARGETS = gen_ft8
+TARGETS = gen_ft8 decode_ft8 test
-#TARGETS = gen_ft8 decode_ft8 test
 
 .PHONY: run_tests all clean
 
@@ -14,11 +13,11 @@ run_tests: test
 gen_ft8: gen_ft8.o ft8/constants.o ft8/text.o ft8/pack_v2.o ft8/encode_v2.o common/wave.o
 	$(CXX) $(LDFLAGS) -o $@ $^
 
-#test:  test.o ft8/encode.o ft8/pack.o ft8/text.o ft8/pack_v2.o ft8/encode_v2.o ft8/unpack.o
+test:  test.o ft8/v1/encode.o ft8/v1/pack.o ft8/v1/unpack.o ft8/pack_v2.o ft8/encode_v2.o ft8/text.o ft8/constants.o
-#	$(CXX) $(LDFLAGS) -o $@ $^
+	$(CXX) $(LDFLAGS) -o $@ $^
 
-#decode_ft8: decode_ft8.o fft/kiss_fftr.o fft/kiss_fft.o ft8/ldpc.o ft8/unpack.o ft8/text.o common/wave.o
+decode_ft8: decode_ft8.o fft/kiss_fftr.o fft/kiss_fft.o ft8/ldpc.o ft8/unpack_v2.o ft8/text.o ft8/constants.o common/wave.o
-#	$(CXX) $(LDFLAGS) -o $@ $^
+	$(CXX) $(LDFLAGS) -o $@ $^
 
 clean:
 	rm -f *.o ft8/*.o common/*.o fft/*.o $(TARGETS)

decode_ft8.cpp

@@ -3,16 +3,12 @@
 #include <cstdio>
 #include <cmath>
 
-#include "common/wave.h"
+#include "ft8/unpack_v2.h"
-#include "ft8/pack.h"
-#include "ft8/encode.h"
-#include "ft8/pack_v2.h"
-#include "ft8/encode_v2.h"
-#include "ft8/unpack.h"
-
 #include "ft8/ldpc.h"
-#include "fft/kiss_fftr.h"
+#include "ft8/constants.h"
 
+#include "common/wave.h"
+#include "fft/kiss_fftr.h"
 
 void usage() {
     printf("Decode a 15-second WAV file.\n");
@@ -79,11 +75,11 @@ void heapify_up(Candidate * heap, int heap_size) {
 
 // Find 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).
-void find_sync(const uint8_t *power, int num_blocks, int num_bins, const uint8_t *sync_map, int num_candidates, Candidate *heap) {
+int find_sync(const uint8_t *power, int num_blocks, int num_bins, const uint8_t *sync_map, int num_candidates, Candidate *heap) {
     int heap_size = 0;
 
     for (int alt = 0; alt < 4; ++alt) {
-        for (int time_offset = 0; time_offset < num_blocks - NN; ++time_offset) {
+        for (int time_offset = 0; time_offset < num_blocks - FT8_NN; ++time_offset) {
             for (int freq_offset = 0; freq_offset < num_bins - 8; ++freq_offset) {
                 int score = 0;
 
@@ -122,6 +118,8 @@ void find_sync(const uint8_t *power, int num_blocks, int num_bins, const uint8_t
             }
         }
     }
+
+    return heap_size;
 }
 
 
@@ -134,6 +132,9 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t *
     for (int i = 0; i < nfft; ++i) {
         window[i] = hann_i(i, nfft);
     }
+    // for (int i = 0; i < nfft; ++i) {
+    //     window[i] = (i < block_size) ? 2 * hann_i(i, block_size) : 0.0f;
+    // }
 
     size_t  fft_work_size;
     kiss_fftr_alloc(nfft, 0, 0, &fft_work_size);
@@ -141,11 +142,12 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t *
     printf("N_FFT = %d\n", nfft);
     printf("FFT work area = %lu\n", fft_work_size);
 
-    void *  fft_work = malloc(fft_work_size);
+    void *       fft_work = malloc(fft_work_size);
     kiss_fftr_cfg fft_cfg = kiss_fftr_alloc(nfft, 0, fft_work, &fft_work_size);
 
     int offset = 0;
     float fft_norm = 1.0f / nfft;
+    float max_mag = -100.0f;
     for (int i = 0; i < num_blocks; ++i) {
         // Loop over two possible time offsets (0 and block_size/2)
         for (int time_sub = 0; time_sub <= block_size/2; time_sub += block_size/2) {
@@ -162,8 +164,8 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t *
 
             // Compute log magnitude in decibels
             for (int j = 0; j < nfft/2 + 1; ++j) {
-                float mag2 = fft_norm * (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r);
+                float mag2 = (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r);
-                mag_db[j] = 10.0f * log10f(1.0E-10f + mag2);
+                mag_db[j] = 10.0f * log10f(1.0E-10f + mag2 * fft_norm);
             }
 
             // Loop over two possible frequency bin offsets (for averaging)
@@ -173,15 +175,19 @@ void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t *
                     float db2 = mag_db[j * 2 + freq_sub + 1];
                     float db = (db1 + db2) / 2;
 
-                    // Scale decibels to unsigned 8-bit range
+                    // Scale decibels to unsigned 8-bit range and clamp the value
                     int scaled = (int)(2 * (db + 100));
                     power[offset] = (scaled < 0) ? 0 : ((scaled > 255) ? 255 : scaled);
                     ++offset;
+
+                    if (db > max_mag) max_mag = db;
                 }
             }
         }
     }
 
+    printf("Max magnitude: %.1f dB\n", max_mag);
+
     free(fft_work);
 }
 
@@ -203,7 +209,7 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca
 
     int k = 0;
     // Go over FSK tones and skip Costas sync symbols
-    for (int i = 7; i < NN - 7; ++i) {
+    for (int i = 7; i < FT8_NN - 7; ++i) {
         if (i == 36) i += 7;
 
         // Pointer to 8 bins of the current symbol
@@ -222,14 +228,15 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca
         // printf("%d %d %d %d %d %d %d %d : %.0f %.0f %.0f\n", 
         //         ps[0], ps[1], ps[2], ps[3], ps[4], ps[5], ps[6], ps[7], 
         //         log174[k + 0], log174[k + 1], log174[k + 2]);
+
         k += 3;
     }
 
     // Compute the variance of log174
     float sum   = 0;
     float sum2  = 0;
-    float inv_n = 1.0f / (3 * ND);
+    float inv_n = 1.0f / (3 * FT8_ND);
-    for (int i = 0; i < 3 * ND; ++i) {
+    for (int i = 0; i < 3 * FT8_ND; ++i) {
         sum  += log174[i];
         sum2 += log174[i] * log174[i];
     }
@@ -238,7 +245,7 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca
     // Normalize log174 such that sigma = 2.83 (Why? It's in WSJT-X)
     float norm_factor = 2.83f / sqrtf(variance);
 
-    for (int i = 0; i < 3 * ND; ++i) {
+    for (int i = 0; i < 3 * FT8_ND; ++i) {
         log174[i] *= norm_factor;
         //printf("%.1f ", log174[i]);
     }
@@ -246,6 +253,35 @@ void extract_likelihood(const uint8_t *power, int num_bins, const Candidate & ca
 }
 
 
+void test_tones(float *log174) {
+    for (int i = 0; i < FT8_ND; ++i) {
+        const uint8_t inv_map[8] = {0, 1, 3, 2, 6, 4, 5, 7};
+        uint8_t tone = ("0000000011721762454112705354533170166234757420515470163426"[i]) - '0';
+        uint8_t b3 = inv_map[tone];
+        log174[3 * i]     = (b3 & 4) ? +1.0 : -1.0;
+        log174[3 * i + 1] = (b3 & 2) ? +1.0 : -1.0;
+        log174[3 * i + 2] = (b3 & 1) ? +1.0 : -1.0;
+    }    
+    // 3140652 00000000117217624541127053545 3140652 33170166234757420515470163426 3140652
+    // 0000000011721762454112705354533170166234757420515470163426
+    // 0000000011721762454112705454544170166344757430515470073537
+    // 0000000011711761444111704343433170166233747320414370072427
+    // 0000000011711761454111705353533170166233757320515370072527
+}
+
+
+void print_tones(const uint8_t *code_map, const float *log174) {
+    for (int k = 0; k < 3 * FT8_ND; k += 3) {
+        uint8_t max = 0;
+        if (log174[k + 0] > 0) max |= 4;
+        if (log174[k + 1] > 0) max |= 2;
+        if (log174[k + 2] > 0) max |= 1;
+        printf("%d", code_map[max]);
+    }
+    printf("\n");
+}
+
+
 int main(int argc, char ** argv) {
     // Expect one command-line argument
     if (argc < 2) {
@@ -264,13 +300,6 @@ int main(int argc, char ** argv) {
         return -1;
     }
 
-    // Costas 7x7 tone pattern
-    const uint8_t kCostas_map_v1[] = { 2,5,6,0,4,1,3 };
-    const uint8_t kCostas_map_v2[] = { 3,1,4,0,6,5,2 };
-    // Gray maps (used only in v2)
-    const uint8_t kGray_map_v1[8] = { 0,1,2,3,4,5,6,7 };    // identity map
-    const uint8_t kGray_map_v2[8] = { 0,1,3,2,5,6,4,7 };
-
     const float fsk_dev = 6.25f;
 
     const int num_bins = (int)(sample_rate / (2 * fsk_dev));
@@ -282,41 +311,45 @@ int main(int argc, char ** argv) {
 
     extract_power(signal, num_blocks, num_bins, power);
 
-    int num_candidates = 250;
+    int num_candidates = 100;
     Candidate heap[num_candidates];
 
-    find_sync(power, num_blocks, num_bins, kCostas_map_v1, num_candidates, heap);
+    find_sync(power, num_blocks, num_bins, kCostas_map, num_candidates, heap);
 
     for (int idx = 0; idx < num_candidates; ++idx) {
         Candidate &cand = heap[idx];
 
-        float log174[3 * ND];
+        float log174[3 * FT8_ND];
-        extract_likelihood(power, num_bins, cand, kGray_map_v1, log174);
+        extract_likelihood(power, num_bins, cand, kGray_map, log174);
 
-        const int num_iters = 20;
+        const int num_iters = 25;
-        int plain[3 * ND];
+        uint8_t plain[3 * FT8_ND];
-        int ok = 0;
+        int n_errors = 0;
 
-        bp_decode(log174, num_iters, plain, &ok);
-        //ldpc_decode(log174, num_iters, plain, &ok);
-        //printf("ldpc_decode() = %d\n", ok);
-
-        if (ok == 87) {
             float freq_hz  = (cand.freq_offset + cand.freq_sub / 2.0f) * fsk_dev;
             float time_sec = (cand.time_offset + cand.time_sub / 2.0f) / fsk_dev;
-            //printf("%03d: score = %d freq = %.1f time = %.2f\n", idx, 
-            //        cand.score, freq_hz, time_sec);
 
-            uint8_t a87[11];
+            printf("%03d: score = %d freq = %.1f time = %.2f\n", idx, 
+                    cand.score, freq_hz, time_sec);
+
+        bp_decode(log174, num_iters, plain, &n_errors);
+        //ldpc_decode(log174, num_iters, plain, &n_errors);
+        printf("ldpc_decode() = %d\n", n_errors);
+
+        if (n_errors == 0) {
+
+            //print_tones(kGray_map, log174);
+            
+            // Extract payload + CRC
+            uint8_t a91[12];
             uint8_t mask = 0x80;
             uint8_t position = 0;
-            for (int i = 0; i < 11; ++i) {
+            for (int i = 0; i < 12; ++i) {
-                a87[i] = 0;
+                a91[i] = 0;
             }
-            // Extract payload + CRC (last 87 bits)
+            for (int i = 0; i < FT8_K; ++i) {
-            for (int i = 174 - 87; i < 174; ++i) {
                 if (plain[i]) {
-                    a87[position] |= mask;
+                    a91[position] |= mask;
                 }
                 mask >>= 1;
                 if (!mask) {
@@ -325,14 +358,17 @@ int main(int argc, char ** argv) {
                 }
             }
 
-            for (int i = 0; i < 11; ++i) {
+            // TODO: check CRC
-                //printf("%02x ", a87[i]);
+
-            }
+            // for (int i = 0; i < 12; ++i) {
-            //printf("\n");
+            //     printf("%02x ", a91[i]);
+            // }
+            // printf("\n");
 
             char message[20];
-            unpack(a87, message);
+            unpack77(a91, message);
 
+            // Fake WSJT-X-like output for now
             printf("000000   0 %4.1f %4d ~  %s\n", time_sec, (int)(freq_hz + 0.5f), message);
         }
     }

ft8/ldpc.cpp

@@ -14,9 +14,190 @@
 #include <stdlib.h>
 #include "constants.h"
 
-int ldpc_check(int codeword[]);
+int ldpc_check(uint8_t codeword[]);
+float fast_tanh(float x);
+float fast_atanh(float x);
 
 
+// codeword is 174 log-likelihoods.
+// plain is a return value, 174 ints, to be 0 or 1.
+// max_iters is how hard to try.
+// 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;
+
+    for (int j = 0; j < FT8_M; j++) {
+        for (int i = 0; i < FT8_N; i++) {
+            m[j][i] = codeword[i];
+            e[j][i] = 0.0f;
+        }
+    }
+
+    for (int iter = 0; iter < max_iters; iter++) {
+        for (int j = 0; j < FT8_M; j++) {
+            for (int ii1 = 0; ii1 < kNrw[j]; ii1++) {
+                int i1 = kNm[j][ii1] - 1;
+                float a = 1.0f;
+                for (int ii2 = 0; ii2 < kNrw[j]; ii2++) {
+                    int i2 = kNm[j][ii2] - 1;
+                    if (i2 != i1) {
+                        a *= fast_tanh(-m[j][i2] / 2.0f);
+                    }
+                }
+                e[j][i1] = logf((1 - a) / (1 + a));
+            }
+        }
+
+        uint8_t cw[FT8_N];
+        for (int i = 0; i < FT8_N; i++) {
+            float l = codeword[i];
+            for (int j = 0; j < 3; j++)
+                l += e[kMn[i][j] - 1][i];
+            cw[i] = (l > 0) ? 1 : 0;
+        }
+
+        int errors = ldpc_check(cw);
+
+        if (errors < min_errors) {
+            // Update the current best result
+            for (int i = 0; i < FT8_N; i++) {
+                plain[i] = cw[i];
+            }
+            min_errors = errors;
+
+            if (errors == 0) {
+                break;  // Found a perfect answer
+            }
+        }
+
+        for (int i = 0; i < FT8_N; i++) {
+            for (int ji1 = 0; ji1 < 3; ji1++) {
+                int j1 = kMn[i][ji1] - 1;
+                float l = codeword[i];
+                for (int ji2 = 0; ji2 < 3; ji2++) {
+                    if (ji1 != ji2) {
+                        int j2 = kMn[i][ji2] - 1;
+                        l += e[j2][i];
+                    }
+                }
+                m[j1][i] = l;
+            }
+        }
+    }
+
+    *ok = min_errors;
+}
+
+
+//
+// does a 174-bit codeword pass the FT8's LDPC parity checks?
+// returns the number of parity errors.
+// 0 means total success.
+//
+int ldpc_check(uint8_t codeword[]) {
+    int errors = 0;
+
+    // kNm[87][7]
+    for (int j = 0; j < FT8_M; ++j) {
+        uint8_t x = 0;
+        for (int ii1 = 0; ii1 < kNrw[j]; ++ii1) {
+            x ^= codeword[kNm[j][ii1] - 1];
+        }
+        if (x != 0) {
+            ++errors;
+        }
+    }
+    return errors;
+}
+
+
+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;
+        }
+    }
+
+    for (int iter = 0; iter < max_iters; ++iter) {
+        float   zn[FT8_N];
+        uint8_t cw[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];
+            cw[i] = (zn[i] > 0) ? 1 : 0;
+        }
+
+        // Check to see if we have a codeword (check before we do any iter)
+        int errors = ldpc_check(cw);
+
+        if (errors < min_errors) {
+            // we have a better guess - update the result
+            for (int i = 0; i < FT8_N; i++) {
+                plain[i] = cw[i];
+            }
+            min_errors = errors;
+
+            if (errors == FT8_M) {
+                break;  // Found a perfect answer
+            }
+        }
+
+        // Send messages from bits to check nodes 
+        for (int i = 0; i < FT8_M; ++i) {
+            for (int j = 0; j < kNrw[i]; ++j) {
+                int ibj = kNm[i][j] - 1;
+                toc[i][j] = zn[ibj];
+                for (int kk = 0; kk < 3; ++kk) { 
+                    // subtract off what the bit had received from the check
+                    if (kMn[ibj][kk] - 1 == i) {
+                        toc[i][j] -= tov[ibj][kk];
+                    }
+                }
+            }
+        }
+        
+        // send messages from check nodes to variable nodes
+        for (int i = 0; i < FT8_M; ++i) {
+            for (int j = 0; j < kNrw[i]; ++j) {
+                toc[i][j] = fast_tanh(-toc[i][j] / 2);
+            }
+        }
+
+        for (int i = 0; i < FT8_N; ++i) {
+            for (int j = 0; j < 3; ++j) {
+                int ichk = kMn[i][j] - 1; // kMn(:,j) are the checks that include bit j
+                float Tmn = 1.0f;
+                for (int k = 0; k < kNrw[ichk]; ++k) {
+                    if (kNm[ichk][k] - 1 != i) {
+                        Tmn *= toc[ichk][k];
+                    }
+                }
+                tov[i][j] = 2 * fast_atanh(-Tmn);
+            }
+        }
+    }
+
+    *ok = min_errors;
+}
+
 // https://varietyofsound.wordpress.com/2011/02/14/efficient-tanh-computation-using-lamberts-continued-fraction/
 // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/10/0001/
 // https://mathr.co.uk/blog/2017-09-06_approximating_hyperbolic_tangent.html
@@ -95,192 +276,4 @@ float platanh(float x) {
         return isign * (x - 0.9914f) / 0.0012f;
     }
     return isign * 7.0f;
-}
+}
-
-
-// codeword is 174 log-likelihoods.
-// plain is a return value, 174 ints, to be 0 or 1.
-// max_iters is how hard to try.
-// ok == 87 means success.
-void ldpc_decode(float codeword[], int max_iters, int plain[], int *ok) {
-    float m[FT8_M][FT8_N];       // ~60 kB
-    float e[FT8_M][FT8_N];       // ~60 kB
-    int best_score = -1;
-
-    for (int j = 0; j < FT8_M; j++) {
-        for (int i = 0; i < FT8_N; i++) {
-            m[j][i] = codeword[i];
-            e[j][i] = 0.0f;
-        }
-    }
-
-    for (int iter = 0; iter < max_iters; iter++) {
-        for (int j = 0; j < FT8_M; j++) {
-            for (int ii1 = 0; ii1 < kNrw[j]; ii1++) {
-                int i1 = kNm[j][ii1] - 1;
-                float a = 1.0f;
-                for (int ii2 = 0; ii2 < kNrw[j]; ii2++) {
-                    int i2 = kNm[j][ii2] - 1;
-                    if (i2 != i1) {
-                        a *= fast_tanh(-m[j][i2] / 2.0f);
-                    }
-                }
-                e[j][i1] = logf((1 - a) / (1 + a));
-            }
-        }
-
-        int cw[FT8_N];
-        for (int i = 0; i < FT8_N; i++) {
-            float l = codeword[i];
-            for (int j = 0; j < 3; j++)
-                l += e[kMn[i][j] - 1][i];
-            cw[i] = (l > 0) ? 1 : 0;
-        }
-
-        int score = ldpc_check(cw);
-
-        if (score > best_score) {
-            // Update the current best result
-            for (int i = 0; i < FT8_N; i++) {
-                //plain[i] = cw[colorder[i]];   // Reverse the column permutation
-                plain[i] = cw[i];
-            }
-            best_score = score;
-
-            if (score == FT8_M) {
-                // Found a perfect answer
-                break;
-            }
-        }
-
-
-        for (int i = 0; i < FT8_N; i++) {
-            for (int ji1 = 0; ji1 < 3; ji1++) {
-                int j1 = kMn[i][ji1] - 1;
-                float l = codeword[i];
-                for (int ji2 = 0; ji2 < 3; ji2++) {
-                    if (ji1 != ji2) {
-                        int j2 = kMn[i][ji2] - 1;
-                        l += e[j2][i];
-                    }
-                }
-                m[j1][i] = l;
-            }
-        }
-    }
-
-    *ok = best_score;
-}
-
-
-//
-// does a 174-bit codeword pass the FT8's LDPC parity checks?
-// returns the number of parity checks that passed.
-// 87 means total success.
-//
-int ldpc_check(int codeword[]) {
-    int score = 0;
-
-    // kNm[87][7]
-    for (int j = 0; j < FT8_M; ++j) {
-        int x = 0;
-        for (int ii1 = 0; ii1 < kNrw[j]; ++ii1) {
-            x ^= codeword[kNm[j][ii1] - 1];
-        }
-        if (x == 0) {
-            ++score;
-        }
-    }
-    return score;
-}
-
-
-void bp_decode(float codeword[], int max_iters, int plain[], int *ok) {
-    float tov[FT8_N][3];
-    float toc[FT8_M][7];
-
-    int best_score = -1;
-
-    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;
-        }
-    }
-
-    for (int iter = 0; iter < max_iters; ++iter) {
-        float   zn[FT8_N];
-        int     cw[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];
-            cw[i] = (zn[i] > 0) ? 1 : 0;
-        }
-
-        // Check to see if we have a codeword (check before we do any iter)
-        int score = ldpc_check(cw);
-
-        if (score > best_score) {
-            // we have a better guess - update the result
-            for (int i = 0; i < FT8_N; i++) {
-                //plain[i] = cw[colorder[i]];   // Reverse the column permutation
-                plain[i] = cw[i];
-            }
-            best_score = score;
-
-            if (score == FT8_M) {
-                break;
-            }
-        }
-
-        // Send messages from bits to check nodes 
-        for (int i = 0; i < FT8_M; ++i) {
-            for (int j = 0; j < kNrw[i]; ++j) {
-                int ibj = kNm[i][j] - 1;
-                toc[i][j] = zn[ibj];
-                for (int kk = 0; kk < 3; ++kk) { 
-                    // subtract off what the bit had received from the check
-                    if (kMn[ibj][kk] - 1 == i) {
-                        toc[i][j] -= tov[ibj][kk];
-                    }
-                }
-            }
-        }
-        
-        // send messages from check nodes to variable nodes
-        for (int i = 0; i < FT8_M; ++i) {
-            for (int j = 0; j < kNrw[i]; ++j) {
-                //toc[i][j] = pltanh(-toc[i][j] / 2);
-                toc[i][j] = fast_tanh(-toc[i][j] / 2);
-                //toc[i][j] = tanhf(-toc[i][j] / 2);
-            }
-        }
-
-        for (int i = 0; i < FT8_N; ++i) {
-            for (int j = 0; j < 3; ++j) {
-                int ichk = kMn[i][j] - 1; // kMn(:,j) are the checks that include bit j
-                float Tmn = 1.0f;
-                for (int k = 0; k < kNrw[ichk]; ++k) {
-                    if (kNm[ichk][k] - 1 != i) {
-                        Tmn *= toc[ichk][k];
-                    }
-                }
-                //tov[i][j] = 2 * platanh(-Tmn);
-                tov[i][j] = 2 * fast_atanh(-Tmn);
-                //tov[i][j] = 2 * atanhf(-Tmn);
-            }
-        }
-    }
-
-    *ok = best_score;
-}

ft8/ldpc.h

@@ -4,6 +4,6 @@
 // plain is a return value, 174 ints, to be 0 or 1.
 // iters is how hard to try.
 // ok == 87 means success.
-void ldpc_decode(float codeword[], int max_iters, int plain[], int *ok);
+void ldpc_decode(float codeword[], int max_iters, uint8_t plain[], int *ok);
 
-void bp_decode(float codeword[], int max_iters, int plain[], int *ok);
+void bp_decode(float codeword[], int max_iters, uint8_t plain[], int *ok);

ft8/pack_v2.cpp

@@ -9,11 +9,12 @@
 namespace ft8_v2 {
 
 // TODO: This is wasteful, should figure out something more elegant
-const char *A0 = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ+-./?";
+const char A0[] = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ+-./?";
-const char *A1 = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
+const char A1[] = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
-const char *A2 = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
+const char A2[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
-const char *A3 = "0123456789";
+const char A3[] = "0123456789";
-const char *A4 = " ABCDEFGHIJKLMNOPQRSTUVWXYZ";
+const char A4[] = " ABCDEFGHIJKLMNOPQRSTUVWXYZ";
+
 
 int index(const char *string, char c) {
     for (int i = 0; *string; ++i, ++string) {
@@ -24,6 +25,7 @@ int index(const char *string, char c) {
     return -1;  // Not found
 }
 
+
 // Pack a special token, a 22-bit hash code, or a valid base call 
 // into a 28-bit integer.
 int32_t pack28(const char *callsign) {
@@ -44,8 +46,8 @@ int32_t pack28(const char *callsign) {
     }
 
     // TODO: Check for <...> callsign
-    // if(c13(1:1).eq.'<')then
+    // if(text(1:1).eq.'<')then
-    //   call save_hash_call(c13,n10,n12,n22)   !Save callsign in hash table
+    //   call save_hash_call(text,n10,n12,n22)   !Save callsign in hash table
     //   n28=NTOKENS + n22
 
     char c6[6] = {' ', ' ', ' ', ' ', ' ', ' '};
@@ -55,21 +57,24 @@ int32_t pack28(const char *callsign) {
         length++;
     }
 
-    // Work-around for Swaziland prefix:
+    // Copy callsign to 6 character buffer
     if (starts_with(callsign, "3DA0") && length <= 7) {
+        // Work-around for Swaziland prefix: 3DA0XYZ -> 3D0XYZ
         memcpy(c6, "3D0", 3);
         memcpy(c6 + 3, callsign + 4, length - 4);
     }
-    // Work-around for Guinea prefixes:
     else if (starts_with(callsign, "3X") && is_letter(callsign[2]) && length <= 7) {
+        // Work-around for Guinea prefixes: 3XA0XYZ -> QA0XYZ
         memcpy(c6, "Q", 1);
         memcpy(c6 + 1, callsign + 2, length - 2);
     }
     else {
         if (is_digit(callsign[2]) && length <= 6) {
+            // AB0XYZ
             memcpy(c6, callsign, length);
         }
         else if (is_digit(callsign[1]) && length <= 5) {
+            // A0XYZ -> " A0XYZ"
             memcpy(c6 + 1, callsign, length);
         }
     }
@@ -80,7 +85,7 @@ int32_t pack28(const char *callsign) {
         (i2 = index(A3, c6[2])) >= 0 && (i3 = index(A4, c6[3])) >= 0 &&
         (i4 = index(A4, c6[4])) >= 0 && (i5 = index(A4, c6[5])) >= 0) 
     {
-        printf("Pack28: idx=[%d, %d, %d, %d, %d, %d]\n", i0, i1, i2, i3, i4, i5);
+        //printf("Pack28: idx=[%d, %d, %d, %d, %d, %d]\n", i0, i1, i2, i3, i4, i5);
         // This is a standard callsign
         int32_t n28 = i0;
         n28 = n28 * 36 + i1;
@@ -88,17 +93,17 @@ int32_t pack28(const char *callsign) {
         n28 = n28 * 27 + i3;
         n28 = n28 * 27 + i4;
         n28 = n28 * 27 + i5;
-        printf("Pack28: n28=%d (%04xh)\n", n28, n28);
+        //printf("Pack28: n28=%d (%04xh)\n", n28, n28);
         return NTOKENS + MAX22 + n28;
     }
 
-    //char c13[13];
+    //char text[13];
 
     //if (length > 13) return -1;
 
     // TODO:
     // Treat this as a nonstandard callsign: compute its 22-bit hash
-    // call save_hash_call(c13,n10,n12,n22)   !Save callsign in hash table
+    // call save_hash_call(text,n10,n12,n22)   !Save callsign in hash table
     // n28=NTOKENS + n22
 
     // n28=iand(n28,ishft(1,28)-1)
@@ -213,10 +218,6 @@ int pack77_1(const char *msg, uint8_t *b77) {
     // write(c77,1000) n28a,ipa,n28b,ipb,ir,igrid4,i3
     // 1000 format(2(b28.28,b1),b1,b15.15,b3.3)  
 
-    // 00 00 00 27 b3 00 01 0b 27 8f b9 e0
-    // (00 00 00 2) 0 (f0 85 ab e) (0) 
-    // 7842D5F
-
     b77[0] = (n28a >> 21);
     b77[1] = (n28a >> 13);
     b77[2] = (n28a >> 5);
@@ -232,30 +233,43 @@ int pack77_1(const char *msg, uint8_t *b77) {
 }
 
 
-void packtext77(const char *c13, uint8_t *b71) {
+void packtext77(const char *text, uint8_t *b71) {
-    // TODO: w=adjustr(c13)
+    int length = strlen(text);
+
+    // Skip leading and trailing spaces
+    while (*text == ' ' && *text != 0) {
+        ++text;
+        --length;
+    }
+    while (length > 0 && text[length - 1] == ' ') {
+        --length;
+    }
 
     for (int i = 0; i < 9; ++i) {
         b71[i] = 0;
     }
 
     for (int j = 0; j < 13; ++j) {
-        int q = index(A0, c13[j]);
+        // Multiply the long integer in b71 by 42
-
-        // Multiply b71 by 42
         uint16_t x = 0;
         for (int i = 8; i >= 0; --i) {
             x += b71[i] * (uint16_t)42;
-            b71[i] = x;
+            b71[i] = (x & 0xFF);
             x >>= 8;
         }
 
         // Add index of the current char
-        x = (q > 0) ? q : 0;
+        if (j < length) {
+            int q = index(A0, text[j]);
+            x = (q > 0) ? q : 0;
+        }
+        else {
+            x = 0;
+        }
         for (int i = 8; i >= 0; --i) {
             if (x == 0) break;
             x += b71[i];
-            b71[i] = x;
+            b71[i] = (x & 0xFF);
             x >>= 8;
         }
     }

ft8/text.cpp

@@ -78,13 +78,13 @@ int dd_to_int(const char *str, int length) {
 
 
 // Convert a 2 digit integer to string
-void int_to_dd(char *str, int value, int width) {
+void int_to_dd(char *str, int value, int width, bool full_sign) {
     if (value < 0) {
         *str = '-';
         ++str;
         value = -value;
     }
-    else {
+    else if (full_sign) {
         *str = '+';
         ++str;
     }

ft8/text.h

@@ -17,4 +17,4 @@ void fmtmsg(char *msg_out, const char *msg_in);
 int dd_to_int(const char *str, int length);
 
 // Convert a 2 digit integer to string
-void int_to_dd(char *str, int value, int width);
+void int_to_dd(char *str, int value, int width, bool full_sign = false);

ft8/unpack_v2.cpp

@@ -0,0 +1,322 @@
+#include "unpack_v2.h"
+#include "text.h"
+
+#include <string.h>
+
+//const uint32_t NBASE = 37L*36L*10L*27L*27L*27L;
+const uint32_t MAX22    = 4194304L;
+const uint32_t NTOKENS  = 2063592L;
+const uint16_t MAXGRID4 = 32400L;
+
+
+// convert integer index to ASCII character according to one of 5 tables:
+// table 0: " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ+-./?"
+// table 1: " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+// table 2: "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+// table 3: "0123456789"
+// table 4: " ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+char charn(int c, int table_idx) {
+    if (table_idx == 0 || table_idx == 1 || table_idx == 4) {
+        if (c == 0) return ' ';
+        c -= 1;
+    }
+    if (table_idx != 4) {
+        if (c < 10) return '0' + c;
+        c -= 10;
+    }
+    if (table_idx != 3) {
+        if (c < 26) return 'A' + c;
+        c -= 26;
+    }
+    if (table_idx == 0) {
+        if (c < 5) return "+-./?" [c];
+    }
+
+    return '_'; // unknown character, should never get here
+}
+
+
+// n28 is a 28-bit integer, e.g. n28a or n28b, containing all the
+// call sign bits from a packed message.
+int unpack28(uint32_t n28, char *result) {
+    // Check for special tokens DE, QRZ, CQ, CQ_nnn, CQ_aaaa
+    if (n28 < NTOKENS) {
+        if (n28 <= 2) {
+            if (n28 == 0) strcpy(result, "DE");
+            if (n28 == 1) strcpy(result, "QRZ");
+            if (n28 == 2) strcpy(result, "CQ");
+            return 0;   // Success
+        }
+        if (n28 <= 1002) {
+            // CQ_nnn with 3 digits
+            strcpy(result, "CQ ");
+            int_to_dd(result + 3, n28 - 3, 3);
+            return 0;   // Success
+        }
+        if (n28 <= 532443L) {
+            // CQ_aaaa with 4 alphanumeric symbols
+            uint32_t n = n28 - 1003;
+            char aaaa[5];
+
+            aaaa[4] = '\0';
+            aaaa[3] = charn(n % 27, 4);
+            n /= 27;
+            aaaa[2] = charn(n % 27, 4);
+            n /= 27;
+            aaaa[1] = charn(n % 27, 4);
+            n /= 27;
+            aaaa[0] = charn(n % 27, 4);
+
+            // Skip leading whitespace
+            int ws_len = 0;
+            while (aaaa[ws_len] == ' ') {
+                ws_len++;
+            }
+
+            strcpy(result, "CQ ");
+            strcat(result, aaaa + ws_len);
+            return 0;   // Success
+        }
+        // ? TODO: unspecified in the WSJT-X code
+        return -1;
+    }
+
+    n28 = n28 - NTOKENS;
+    if (n28 < MAX22) {
+        // This is a 22-bit hash of a result
+        //n22=n28
+        //call hash22(n22,c13)     !Retrieve result from hash table
+        // TODO: implement
+        return -2;
+    }
+
+    // Standard callsign
+    uint32_t n = n28 - MAX22;
+
+    char callsign[7];
+    callsign[6] = '\0';
+    callsign[5] = charn(n % 27, 4);
+    n /= 27;
+    callsign[4] = charn(n % 27, 4);
+    n /= 27;
+    callsign[3] = charn(n % 27, 4);
+    n /= 27;
+    callsign[2] = charn(n % 10, 3);
+    n /= 10;
+    callsign[1] = charn(n % 36, 2);
+    n /= 36;
+    callsign[0] = charn(n % 37, 1);
+
+    // Skip trailing and leading whitespace in case of a short callsign
+    int ws_len = 0;
+    while (ws_len <= 5 && callsign[5 - ws_len] == ' ') {
+        callsign[5 - ws_len] = '\0';
+        ws_len++;
+    }
+    ws_len = 0;
+    while (callsign[ws_len] == ' ') {
+        ws_len++;
+    }
+
+    strcpy(result, callsign + ws_len);
+    return 0;   // Success
+}
+
+
+int unpack_type1(const uint8_t *a77, uint8_t i3, char *message) {
+    uint32_t n28a, n28b;
+    uint16_t igrid4;
+    uint8_t  ir;
+    
+    // Extract packed fields
+    // read(c77,1000) n28a,ipa,n28b,ipb,ir,igrid4,i3
+    // 1000 format(2(b28,b1),b1,b15,b3)
+    n28a  = (a77[0] << 21);
+    n28a |= (a77[1] << 13);
+    n28a |= (a77[2] << 5);
+    n28a |= (a77[3] >> 3);
+    n28b  = ((a77[3] & 0x07) << 26);
+    n28b |= (a77[4] << 18);
+    n28b |= (a77[5] << 10);
+    n28b |= (a77[6] << 2);
+    n28b |= (a77[7] >> 6);
+    ir      = ((a77[7] & 0x20) >> 5);
+    igrid4  = ((a77[7] & 0x1F) << 10);
+    igrid4 |= (a77[8] << 2);
+    igrid4 |= (a77[9] >> 6);
+
+    // Unpack both callsigns
+    char field_1[14];
+    char field_2[14];
+    if (unpack28(n28a >> 1, field_1) < 0) {
+        return -1;
+    }
+    if (unpack28(n28b >> 1, field_2) < 0) {
+        return -2;
+    }
+    // Fix "CQ_" to "CQ " -> already done in unpack28()
+    // if (starts_with(field_1, "CQ_")) {
+    //     field_1[2] = ' ';
+    // }
+
+    // Append first two fields to the result
+    strcpy(message, field_1);
+    strcat(message, " ");
+    strcat(message, field_2);
+
+    // TODO: add /P and /R suffixes
+    //  if(index(field_1,'<').le.0) then
+    //     ws_len=index(field_1,' ')
+    //     if(ws_len.ge.4 .and. ipa.eq.1 .and. i3.eq.1) field_1(ws_len:ws_len+1)='/R'
+    //     if(ws_len.ge.4 .and. ipa.eq.1 .and. i3.eq.2) field_1(ws_len:ws_len+1)='/P'
+    //     if(ws_len.ge.4) call add_call_to_recent_calls(field_1)
+    //  endif
+    //  if(index(field_2,'<').le.0) then
+    //     ws_len=index(field_2,' ')
+    //     if(ws_len.ge.4 .and. ipb.eq.1 .and. i3.eq.1) field_2(ws_len:ws_len+1)='/R'
+    //     if(ws_len.ge.4 .and. ipb.eq.1 .and. i3.eq.2) field_2(ws_len:ws_len+1)='/P'
+    //     if(ws_len.ge.4) call add_call_to_recent_calls(field_2)
+    //  endif
+
+    char field_3[5];
+    if (igrid4 <= MAXGRID4) {
+        // Extract 4 symbol grid locator
+        field_3[4] = '\0';
+
+        uint16_t n = igrid4;
+        field_3[3] = '0' + (n % 10);
+        n /= 10;
+        field_3[2] = '0' + (n % 10);
+        n /= 10;
+        field_3[1] = 'A' + (n % 18);
+        n /= 18;
+        field_3[0] = 'A' + (n % 18);
+
+        if (ir != 0) {
+            // In case of ir=1 add an " R " before grid
+            strcat(message, " R ");
+        }
+    }
+    else {
+        // Extract report
+        int irpt = igrid4 - MAXGRID4;
+
+        // Check special cases first
+        if (irpt == 1) field_3[0] = '\0';
+        else if (irpt == 2) strcpy(field_3, "RRR");
+        else if (irpt == 3) strcpy(field_3, "RR73");
+        else if (irpt == 4) strcpy(field_3, "73");
+        else if (irpt >= 5) {
+            // Extract signal report as a two digit number with a + or - sign
+            if (ir == 0) {
+                int_to_dd(field_3, irpt - 35, 2, true);
+            }
+            else {
+                field_3[0] = 'R';
+                int_to_dd(field_3 + 1, irpt - 35, 2, true);
+            }
+        }
+    }
+
+    // Append the last field to the result
+    if (strlen(field_3) > 0) {
+        strcat(message, " ");
+        strcat(message, field_3);
+    }
+
+    return 0;       // Success
+}
+
+
+int unpack_text(const uint8_t *a71, char *text) {
+    // TODO: test
+    uint8_t b71[9];
+
+    for (int i = 0; i < 9; ++i) {
+        b71[i] = a71[i];
+    }
+
+    for (int idx = 0; idx < 13; ++idx) {
+        // Divide the long integer in b71 by 42
+        uint16_t rem = 0;
+        for (int i = 8; i >= 0; --i) {
+            rem = (rem << 8) | b71[i];
+            b71[i] = rem / 42;
+            rem    = rem % 42;
+        }
+        text[idx] = charn(rem, 0);
+    }
+
+    text[13] = '\0';
+    return 0;       // Success
+}
+
+
+int unpack_telemetry(const uint8_t *a71, char *telemetry) {
+    uint8_t b71[9];
+
+    // Shift bits in a71 right by 1
+    uint8_t carry = 0;
+    for (int i = 0; i < 9; ++i) {
+        b71[i] = (carry << 7) | (a71[i] >> 1);
+        carry = (a71[i] & 0x01);
+    }
+
+    // Convert b71 to hexadecimal string
+    for (int i = 0; i < 9; ++i) {
+        uint8_t nibble1 = (b71[i] >> 4);
+        uint8_t nibble2 = (b71[i] & 0x0F);
+        char c1 = (nibble1 > 9) ? (nibble1 - 10 + 'A') : nibble1 + '0';
+        char c2 = (nibble2 > 9) ? (nibble2 - 10 + 'A') : nibble2 + '0';
+        telemetry[i * 2] = c1;
+        telemetry[i * 2 + 1] = c2;
+    }
+
+    telemetry[18] = '\0';
+    return 0;
+}
+
+
+int unpack77(const uint8_t *a77, char *message) {
+    uint8_t  n3, i3;
+
+    n3 = (a77[9] >> 6) & 0x07;
+    i3 = (a77[9] >> 3) & 0x07;
+
+    if (i3 == 0 && n3 == 0) {
+        // 0.0  Free text
+        return unpack_text(a77, message);
+    }
+    else if (i3 == 0 && n3 == 1) {
+        // 0.1  K1ABC RR73; W9XYZ <KH1/KH7Z> -11   28 28 10 5       71   DXpedition Mode
+    }
+    else if (i3 == 0 && n3 == 2) {
+        // 0.2  PA3XYZ/P R 590003 IO91NP           28 1 1 3 12 25   70   EU VHF contest
+    }
+    else if (i3 == 0 && (n3 == 3 || n3 == 4)) {
+        // 0.3   WA9XYZ KA1ABC R 16A EMA            28 28 1 4 3 7    71   ARRL Field Day
+        // 0.4   WA9XYZ KA1ABC R 32A EMA            28 28 1 4 3 7    71   ARRL Field Day
+    }
+    else if (i3 == 0 && n3 == 5) {
+        // 0.5   0123456789abcdef01                 71               71   Telemetry (18 hex)
+        return unpack_telemetry(a77, message);
+    }
+    else if (i3 == 1 || i3 == 2) {
+        // Type 1 (standard message) or Type 2 ("/P" form for EU VHF contest)
+        return unpack_type1(a77, i3, message);
+    }
+    else if (i3 == 3) {
+        // Type 3: ARRL RTTY Contest
+    }
+    else if (i3 == 4) {
+        // Type 4: Nonstandard calls, e.g. <WA9XYZ> PJ4/KA1ABC RR73
+        // One hashed call or "CQ"; one compound or nonstandard call with up
+        // to 11 characters; and (if not "CQ") an optional RRR, RR73, or 73.
+
+        // TODO: implement
+        // read(c77,1050) n12,n58,iflip,nrpt,icq
+        // 1050 format(b12,b58,b1,b2,b1)
+    }
+
+    return 0;
+}

ft8/unpack_v2.h

@@ -0,0 +1,7 @@
+#pragma once
+
+#include <stdint.h>
+
+
+// message should have at least 19 bytes allocated (18 characters + zero terminator)
+int unpack77(const uint8_t *a77, char *message);

ft8/v1/pack.cpp

@@ -2,7 +2,7 @@
 #include <cstdio>
 
 #include "pack.h"
-#include "text.h"
+#include "../text.h"
 
 constexpr int32_t NBASE = 37*36*10*27*27*27L;
 constexpr int32_t NGBASE = 180*180L;

ft8/v1/unpack.cpp

@@ -1,5 +1,5 @@
 #include "unpack.h"
-#include "text.h"
+#include "../text.h"
 
 #include <string.h>
 

gen_ft8.cpp

@@ -84,15 +84,17 @@ int main(int argc, char **argv) {
     printf("\n");
 
     // Third, convert the FSK tones into an audio signal
-    const int num_samples = (int)(0.5 + FT8_NN / 6.25 * 12000);
+    const int sample_rate = 12000;
-    const int num_silence = (15 * 12000 - num_samples) / 2;
+    const float symbol_rate = 6.25f;
+    const int num_samples = (int)(0.5f + FT8_NN / symbol_rate * sample_rate);
+    const int num_silence = (15 * sample_rate - num_samples) / 2;
     float signal[num_silence + num_samples + num_silence];
     for (int i = 0; i < num_silence + num_samples + num_silence; i++) {
         signal[i] = 0;
     }
 
-    synth_fsk(tones, FT8_NN, 1000, 6.25, 6.25, 12000, signal + num_silence);
+    synth_fsk(tones, FT8_NN, 1000, symbol_rate, symbol_rate, sample_rate, signal + num_silence);
-    save_wav(signal, num_silence + num_samples + num_silence, 12000, wav_path);
+    save_wav(signal, num_silence + num_samples + num_silence, sample_rate, wav_path);
 
     return 0;
 }

test.cpp

@@ -4,11 +4,11 @@
 #include <cmath>
 
 #include "ft8/text.h"
-#include "ft8/pack.h"
+#include "ft8/v1/pack.h"
+#include "ft8/v1/unpack.h"
+#include "ft8/v1/encode.h"
 #include "ft8/pack_v2.h"
-#include "ft8/encode.h"
 #include "ft8/encode_v2.h"
-#include "ft8/unpack.h"
 
 #include "common/debug.h"