/ decoder.h
decoder.h
1 // SPDX-FileCopyrightText: 2021 Jeff Epler 2 // 3 // SPDX-License-Identifier: GPL-3.0-only 4 5 #pragma once 6 7 #include <array> 8 #include <cassert> 9 #include <cstddef> 10 #include <cstdint> 11 #include <cstdio> 12 #include <ctime> 13 14 template <int N> int mod_diff(int a, int b) { 15 int c = a - b; 16 if (c > N / 2) 17 c -= N; 18 if (c < -N / 2) 19 c += N; 20 return c; 21 } 22 23 template <int N> int mod_between(int lo, int hi, int val) { 24 return mod_diff<N>(lo, val) < 0 && mod_diff<N>(val, hi) < 0; 25 } 26 27 // An efficient circular buffer of N bits. This is used to store the sampled 28 // WWVB signal for statistical purposes. It's also used as the basis of the 29 // circular_symbol_array, which handles multi-bit values 30 template <int N> struct circular_bit_array { 31 bool at(int i) const { 32 assert(i >= 0 && i < N); 33 i += shift; 34 if (i > N) 35 i -= N; 36 int j = i % 32; 37 i /= 32; 38 return data[i] & (1 << j); 39 } 40 41 bool put(bool b) { 42 int i = shift / 32; 43 int j = shift % 32; 44 45 bool result = data[i] & (1 << j); 46 47 if (b) { 48 data[i] |= (1 << j); 49 } else { 50 data[i] &= ~(1 << j); 51 } 52 53 shift += 1; 54 if (shift == N) 55 shift = 0; 56 57 return result; 58 } 59 60 std::array<uint32_t, (N + 31) / 32> data{}; 61 uint16_t shift{}; 62 }; 63 64 // An efficient circular buffer of N M-bit values. This is used to store the 65 // decoded WWVB signals, which can be 0, 1, or 2 (mark). 66 template <int N, int M> struct circular_symbol_array { 67 int at(int i) const { 68 int result = 0; 69 for (int j = 0; j < M; j++) { 70 result = (result << 1) | data.at(i * M + j); 71 } 72 return result; 73 } 74 75 int put(int v) { 76 int result = 0; 77 for (int j = 0; j < M; j++) { 78 result = (result << 1) | data.put(v & (1 << (M - 1))); 79 v <<= 1; 80 } 81 return result; 82 } 83 84 circular_bit_array<M * N> data{}; 85 }; 86 87 // I guess we need to represent time... 88 struct wwvb_time { 89 int16_t yday; 90 int8_t year, hour, minute, second; 91 int8_t ls, ly, dst, dut1; 92 93 time_t to_utc() const; 94 struct tm apply_zone_and_dst(int zone_offset, bool observe_dst) const; 95 96 int seconds_in_minute() const; 97 void advance_seconds(int n = 1); 98 void advance_minutes(int n = 1); 99 100 bool operator==(const wwvb_time &other) const; 101 }; 102 103 template <size_t SUBSEC_ = 50, size_t SYMBOLS_ = 60, size_t HISTORY_ = 40> 104 struct WWVBDecoder { 105 // The second is divided into units of SUBSEC 106 static constexpr size_t SUBSEC = SUBSEC_; 107 108 // This many WWVB symbols are accumulated 109 static constexpr size_t SYMBOLS = SYMBOLS_; 110 111 // This many whole seconds of symbols are accumulated for statistics. 112 // 5 seconds is too little history, 60 is plenty. 40 seems okay. 113 static constexpr size_t HISTORY = HISTORY_; 114 static constexpr size_t BUFFER = SUBSEC * HISTORY_; 115 116 typedef circular_symbol_array<SYMBOLS, 2> symbol_buffer_type; 117 typedef circular_bit_array<BUFFER> signal_buffer_type; 118 119 // Total number of samples ever received 120 size_t sample_count{}; 121 122 // Total number of symbols ever decoded 123 size_t symbol_count{}; 124 125 // Raw samples from the receiver 126 signal_buffer_type signal{}; 127 128 // Statistical information about the raw samples 129 std::array<int16_t, SUBSEC> counts{}; 130 std::array<int16_t, SUBSEC> edges{}; 131 132 // Statistical information about the symbols 133 int health{}; 134 std::array<uint8_t, SYMBOLS> health_history{}; 135 136 // subsec counts the position modulo SUBSEC; sos is the start-of-second 137 // modulo SUBSEC. tss is the time in ticks since the last second. 138 uint16_t subsec{}, sos{}, tss{}; 139 140 // Decoded symbols 141 symbol_buffer_type symbols{}; 142 143 // Receive a sample `b` from the receiver and process: 144 // * update statistics (counts and edges) incrementally 145 // * check all edges values to update the start-of-second value 146 // * updates the symbols buffer at the start of a new second 147 // Returns true if it is the START of a new WWVB second 148 149 bool update(bool b) { 150 // Put the new bit & extract the old bit 151 sample_count++; 152 auto ob = signal.put(b); 153 154 // Update the counts array 155 if (b && !ob) { 156 counts[subsec]++; 157 } else if (!b && ob) { 158 counts[subsec]--; 159 } 160 161 // Update the edges array 162 auto subsec1 = subsec == SUBSEC - 1 ? 0 : subsec + 1; 163 edges[subsec] = counts[subsec1] - counts[subsec]; 164 165 // Check for sharpest edge 166 // can this be done without a whole array scan? 167 int bi = 0, best = 0; 168 for (size_t i = 0; i < SUBSEC; i++) { 169 if (edges[i] > best) { 170 bi = i; 171 best = edges[i]; 172 } 173 } 174 int osos = sos; 175 sos = bi == SUBSEC - 1 ? 0 : bi + 1; 176 177 subsec = subsec1; 178 179 bool result = false; 180 // If it's been a long time since the last second, fake one. 181 if (tss > SUBSEC) { 182 result = true; 183 } else if (tss > SUBSEC / 2) { 184 // Otherwise, sos may be wandering, so don't repeat a second too 185 // soon 186 result = subsec == sos || subsec == osos; 187 } 188 189 // either reset or increment time-since-second 190 if (result) { 191 tss = 0; 192 decode_symbol(); 193 } else { 194 tss++; 195 } 196 197 return result; 198 } 199 200 // Return how many items from i..j in the raw data array are true 201 // (true represents the reduced-carrier state) 202 int count(int i, int j) { 203 int result = 0; 204 for (; i < j; i++) { 205 result += signal.at(i); 206 } 207 return result; 208 } 209 210 static constexpr int ms_in_subsec(int ms) { 211 return (ms * SUBSEC + SUBSEC / 2) / 1000; 212 } 213 214 static constexpr auto p0 = ms_in_subsec(0); 215 static constexpr auto p1 = ms_in_subsec(200); 216 static constexpr auto p2 = ms_in_subsec(500); 217 static constexpr auto p3 = ms_in_subsec(800); 218 static constexpr auto p4 = ms_in_subsec(1000); 219 220 static constexpr auto la = p1 - p0; 221 static constexpr auto lb = p2 - p1; 222 static constexpr auto lc = p3 - p2; 223 static constexpr auto ld = p4 - p3; 224 225 static constexpr auto MAX_HEALTH = SYMBOLS * SUBSEC; 226 // In around 300 hours of logs from the WWVB observatory, the current 227 // algorithm decoded 16004 minutes (at all, not back-checked for 228 // correctness). Of those, minutes about 86% had health above 97%. That 229 // makes 97% a plausible threshold for a healthy signal. 230 static constexpr auto HEALTH_97PCT = MAX_HEALTH * 97 / 100; 231 232 int check_health(int count, int length, int expect) { 233 return expect ? count : length - count; 234 } 235 236 // A second just concluded, so signal.at(BUFFER-1) is the last sample of 237 // the second, and signal.at(BUFFER-SUBSEC) is the first sample of the 238 // second 239 void decode_symbol() { 240 constexpr auto OFFSET = BUFFER - SUBSEC; 241 #if 0 242 for(size_t i=0; i<SUBSEC; i++) { 243 printf("%c", signal.at(OFFSET + i) ? '_' : '#'); 244 } 245 #endif 246 int count_a = count(OFFSET + p0, OFFSET + p1); 247 int count_b = count(OFFSET + p1, OFFSET + p2); 248 int count_c = count(OFFSET + p2, OFFSET + p3); 249 int count_d = count(OFFSET + p3, OFFSET + p4); 250 251 int result = 0; 252 253 if (count_c > lc / 2) { 254 if (count_b > lb / 2) { 255 result = 2; 256 } else { 257 result = 3; // a nonsense symbol 258 } 259 } else if (count_b > lb / 2) { 260 result = 1; 261 } 262 263 int h = 0; 264 if (result != 3) { 265 h += check_health(count_a, la, 1); 266 h += check_health(count_b, lb, result != 0); 267 h += check_health(count_c, lc, result == 2); 268 h += check_health(count_d, ld, 0); 269 } 270 271 int sc = symbol_count++; 272 int si = sc % SYMBOLS; 273 int oh = health_history[si]; 274 health_history[si] = h; 275 health += (h - oh); 276 277 symbols.put(result); 278 279 #if 0 280 printf(" %d\n", result); 281 for(size_t i=0; i<SYMBOLS; i++) { 282 printf("%d", symbols.at(i)); 283 } 284 printf("\n", result); 285 #endif 286 } 287 288 mutable bool bcderr; 289 290 // Simple BCD-decoder 291 int decode_bcd(int d, int c = -1, int b = -1, int a = -1) const { 292 int r = (a >= 0 ? (symbols.at(SYMBOLS - 60 + a) * 8) : 0) + 293 (b >= 0 ? (symbols.at(SYMBOLS - 60 + b) * 4) : 0) + 294 (c >= 0 ? (symbols.at(SYMBOLS - 60 + c) * 2) : 0) + 295 symbols.at(SYMBOLS - 60 + d) * 1; 296 if (r > 9) 297 bcderr = true; 298 return r; 299 } 300 301 template <class... Ints> 302 int decode_bcd(int d, int c, int b, int a, Ints... rest) const { 303 return decode_bcd(d, c, b, a) + 10 * decode_bcd(rest...); 304 } 305 306 // barebones decoding of some minute-fields 307 inline bool decode_minute(wwvb_time &m) const { 308 for (int i = 0; i < 60; i++) { 309 int sym = symbols.at(SYMBOLS - 60 + i); 310 bool is_mark = (i == 0) || (i % 10 == 9); 311 if (is_mark != (sym == 2)) 312 return false; 313 bool is_zero = (i % 10 == 4) || i == 10 || i == 11 || i == 20 || 314 i == 21 || i == 35; 315 if (is_zero && sym != 0) 316 return false; 317 } 318 319 bcderr = false; 320 m.year = decode_bcd(53, 52, 51, 50, 48, 47, 46, 45); 321 m.yday = decode_bcd(33, 32, 31, 30, 28, 27, 26, 25, 23, 22); 322 m.hour = decode_bcd(18, 17, 16, 15, 13, 12); 323 m.minute = decode_bcd(8, 7, 6, 5, 3, 2, 1); 324 m.ly = decode_bcd(55); 325 m.ls = decode_bcd(56); 326 m.dst = decode_bcd(58, 57); 327 m.second = 0; 328 int abs_dut1 = decode_bcd(43, 42, 41, 40); 329 int dut1_sign = decode_bcd(38, 37, 36); 330 switch (dut1_sign) { 331 case 2: 332 m.dut1 = -abs_dut1; 333 break; 334 case 5: 335 m.dut1 = abs_dut1; 336 break; 337 default: 338 bcderr = true; 339 } 340 return !bcderr; 341 } 342 };