1 /* 2 * refclock_wwv - clock driver for NIST WWV/H time/frequency station 3 */ 4 #ifdef HAVE_CONFIG_H 5 #include <config.h> 6 #endif 7 8 #if defined(REFCLOCK) && defined(CLOCK_WWV) 9 10 #include "ntpd.h" 11 #include "ntp_io.h" 12 #include "ntp_refclock.h" 13 #include "ntp_calendar.h" 14 #include "ntp_stdlib.h" 15 #include "audio.h" 16 17 #include <stdio.h> 18 #include <ctype.h> 19 #include <math.h> 20 #ifdef HAVE_SYS_IOCTL_H 21 # include <sys/ioctl.h> 22 #endif /* HAVE_SYS_IOCTL_H */ 23 24 #define ICOM 1 25 26 #ifdef ICOM 27 #include "icom.h" 28 #endif /* ICOM */ 29 30 /* 31 * Audio WWV/H demodulator/decoder 32 * 33 * This driver synchronizes the computer time using data encoded in 34 * radio transmissions from NIST time/frequency stations WWV in Boulder, 35 * CO, and WWVH in Kauai, HI. Transmissions are made continuously on 36 * 2.5, 5, 10, 15 and 20 MHz in AM mode. An ordinary shortwave receiver 37 * can be tuned manually to one of these frequencies or, in the case of 38 * ICOM receivers, the receiver can be tuned automatically using this 39 * program as propagation conditions change throughout the day and 40 * night. 41 * 42 * The driver receives, demodulates and decodes the radio signals when 43 * connected to the audio codec of a workstation running Solaris, SunOS 44 * FreeBSD or Linux, and with a little help, other workstations with 45 * similar codecs or sound cards. In this implementation, only one audio 46 * driver and codec can be supported on a single machine. 47 * 48 * The demodulation and decoding algorithms used in this driver are 49 * based on those developed for the TAPR DSP93 development board and the 50 * TI 320C25 digital signal processor described in: Mills, D.L. A 51 * precision radio clock for WWV transmissions. Electrical Engineering 52 * Report 97-8-1, University of Delaware, August 1997, 25 pp., available 53 * from www.eecis.udel.edu/~mills/reports.htm. The algorithms described 54 * in this report have been modified somewhat to improve performance 55 * under weak signal conditions and to provide an automatic station 56 * identification feature. 57 * 58 * The ICOM code is normally compiled in the driver. It isn't used, 59 * unless the mode keyword on the server configuration command specifies 60 * a nonzero ICOM ID select code. The C-IV trace is turned on if the 61 * debug level is greater than one. 62 */ 63 /* 64 * Interface definitions 65 */ 66 #define DEVICE_AUDIO "/dev/audio" /* audio device name */ 67 #define AUDIO_BUFSIZ 320 /* audio buffer size (50 ms) */ 68 #define PRECISION (-10) /* precision assumed (about 1 ms) */ 69 #define DESCRIPTION "WWV/H Audio Demodulator/Decoder" /* WRU */ 70 #define SECOND 8000 /* second epoch (sample rate) (Hz) */ 71 #define MINUTE (SECOND * 60) /* minute epoch */ 72 #define OFFSET 128 /* companded sample offset */ 73 #define SIZE 256 /* decompanding table size */ 74 #define MAXSIG 6000. /* max signal level reference */ 75 #define MAXCLP 100 /* max clips above reference per s */ 76 #define MAXSNR 30. /* max SNR reference */ 77 #define DGAIN 20. /* data channel gain reference */ 78 #define SGAIN 10. /* sync channel gain reference */ 79 #define MAXFREQ 1. /* max frequency tolerance (125 PPM) */ 80 #define PI 3.1415926535 /* the real thing */ 81 #define DATSIZ (170 * MS) /* data matched filter size */ 82 #define SYNSIZ (800 * MS) /* minute sync matched filter size */ 83 #define MAXERR 30 /* max data bit errors in minute */ 84 #define NCHAN 5 /* number of radio channels */ 85 #define AUDIO_PHI 5e-6 /* dispersion growth factor */ 86 #ifdef IRIG_SUCKS 87 #define WIGGLE 11 /* wiggle filter length */ 88 #endif /* IRIG_SUCKS */ 89 90 /* 91 * General purpose status bits (status) 92 * 93 * SELV and/or SELH are set when WWV or WWVH has been heard and cleared 94 * on signal loss. SSYNC is set when the second sync pulse has been 95 * acquired and cleared by signal loss. MSYNC is set when the minute 96 * sync pulse has been acquired. DSYNC is set when a digit reaches the 97 * threshold and INSYNC is set when all nine digits have reached the 98 * threshold. The MSYNC, DSYNC and INSYNC bits are cleared only by 99 * timeout, upon which the driver starts over from scratch. 100 * 101 * DGATE is set if a data bit is invalid and BGATE is set if a BCD digit 102 * bit is invalid. SFLAG is set when during seconds 59, 0 and 1 while 103 * probing alternate frequencies. LEPDAY is set when SECWAR of the 104 * timecode is set on 30 June or 31 December. LEPSEC is set during the 105 * last minute of the day when LEPDAY is set. At the end of this minute 106 * the driver inserts second 60 in the seconds state machine and the 107 * minute sync slips a second. The SLOSS and SJITR bits are for monitor 108 * only. 109 */ 110 #define MSYNC 0x0001 /* minute epoch sync */ 111 #define SSYNC 0x0002 /* second epoch sync */ 112 #define DSYNC 0x0004 /* minute units sync */ 113 #define INSYNC 0x0008 /* clock synchronized */ 114 #define FGATE 0x0010 /* frequency gate */ 115 #define DGATE 0x0020 /* data bit error */ 116 #define BGATE 0x0040 /* BCD digit bit error */ 117 #define SFLAG 0x1000 /* probe flag */ 118 #define LEPDAY 0x2000 /* leap second day */ 119 #define LEPSEC 0x4000 /* leap second minute */ 120 121 /* 122 * Station scoreboard bits 123 * 124 * These are used to establish the signal quality for each of the five 125 * frequencies and two stations. 126 */ 127 #define SYNCNG 0x0001 /* sync or SNR below threshold */ 128 #define DATANG 0x0002 /* data or SNR below threshold */ 129 #define ERRRNG 0x0004 /* data error */ 130 #define SELV 0x0100 /* WWV station select */ 131 #define SELH 0x0200 /* WWVH station select */ 132 133 /* 134 * Alarm status bits (alarm) 135 * 136 * These bits indicate various alarm conditions, which are decoded to 137 * form the quality character included in the timecode. If not tracking 138 * second sync, the SYNERR alarm is raised. The data error counter is 139 * incremented for each invalid data bit. If too many data bit errors 140 * are encountered in one minute, the MODERR alarm is raised. The DECERR 141 * alarm is raised if a maximum likelihood digit fails to compare with 142 * the current clock digit. If the probability of any miscellaneous bit 143 * or any digit falls below the threshold, the SYMERR alarm is raised. 144 */ 145 #define DECERR 1 /* BCD digit compare error */ 146 #define SYMERR 2 /* low bit or digit probability */ 147 #define MODERR 4 /* too many data bit errors */ 148 #define SYNERR 8 /* not synchronized to station */ 149 150 /* 151 * Watchcat timeouts (watch) 152 * 153 * If these timeouts expire, the status bits are mashed to zero and the 154 * driver starts from scratch. Suitably more refined procedures may be 155 * developed in future. All these are in minutes. 156 */ 157 #define ACQSN 5 /* station acquisition timeout */ 158 #define DIGIT 30 /* minute unit digit timeout */ 159 #define HOLD 30 /* reachable timeout */ 160 #define PANIC (2 * 1440) /* panic timeout */ 161 162 /* 163 * Thresholds. These establish the minimum signal level, minimum SNR and 164 * maximum jitter thresholds which establish the error and false alarm 165 * rates of the driver. The values defined here may be on the 166 * adventurous side in the interest of the highest sensitivity. 167 */ 168 #define MTHR 13. /* acquisition signal gate (percent) */ 169 #define TTHR 50. /* tracking signal gate (percent) */ 170 #define ATHR 2000. /* acquisition amplitude threshold */ 171 #define ASNR 6. /* acquisition SNR threshold (dB) */ 172 #define AWND 20. /* acquisition jitter threshold (ms) */ 173 #define AMIN 3 /* min bit count */ 174 #define AMAX 6 /* max bit count */ 175 #define QTHR 2000 /* QSY sync threshold */ 176 #define QSNR 20. /* QSY sync SNR threshold (dB) */ 177 #define XTHR 1000. /* QSY data threshold */ 178 #define XSNR 10. /* QSY data SNR threshold (dB) */ 179 #define STHR 500 /* second sync amplitude threshold */ 180 #define SSNR 10. /* second sync SNR threshold */ 181 #define SCMP 10 /* second sync compare threshold */ 182 #define DTHR 1000 /* bit amplitude threshold */ 183 #define DSNR 10. /* bit SNR threshold (dB) */ 184 #define BTHR 1000 /* digit amplitude threshold */ 185 #define BSNR 3. /* digit likelihood threshold (dB) */ 186 #define BCMP 5 /* digit compare threshold */ 187 188 /* 189 * Tone frequency definitions. The increments are for 4.5-deg sine 190 * table. 191 */ 192 #define MS (SECOND / 1000) /* samples per millisecond */ 193 #define IN100 ((100 * 80) / SECOND) /* 100 Hz increment */ 194 #define IN1000 ((1000 * 80) / SECOND) /* 1000 Hz increment */ 195 #define IN1200 ((1200 * 80) / SECOND) /* 1200 Hz increment */ 196 197 /* 198 * Acquisition and tracking time constants. Usually powers of 2. 199 */ 200 #define MINAVG 8 /* min time constant */ 201 #define MAXAVG 1024 /* max time constant */ 202 #define TCONST 16 /* data bit/digit time constant */ 203 204 /* 205 * Miscellaneous status bits (misc) 206 * 207 * These bits correspond to designated bits in the WWV/H timecode. The 208 * bit probabilities are exponentially averaged over several minutes and 209 * processed by a integrator and threshold. 210 */ 211 #define DUT1 0x01 /* 56 DUT .1 */ 212 #define DUT2 0x02 /* 57 DUT .2 */ 213 #define DUT4 0x04 /* 58 DUT .4 */ 214 #define DUTS 0x08 /* 50 DUT sign */ 215 #define DST1 0x10 /* 55 DST1 leap warning */ 216 #define DST2 0x20 /* 2 DST2 DST1 delayed one day */ 217 #define SECWAR 0x40 /* 3 leap second warning */ 218 219 /* 220 * The on-time synchronization point for the driver is the second epoch 221 * sync pulse produced by the FIR matched filters. As the 5-ms delay of 222 * these filters is compensated, the program delay is 1.1 ms due to the 223 * 600-Hz IIR bandpass filter. The measured receiver delay is 4.7 ms and 224 * the codec delay less than 0.2 ms. The additional propagation delay 225 * specific to each receiver location can be programmed in the fudge 226 * time1 and time2 values for WWV and WWVH, respectively. 227 */ 228 #define PDELAY (.0011 + .0047 + .0002) /* net system delay (s) */ 229 230 /* 231 * Table of sine values at 4.5-degree increments. This is used by the 232 * synchronous matched filter demodulators. The integral of sine-squared 233 * over one complete cycle is PI, so the table is normallized by 1 / PI. 234 */ 235 double sintab[] = { 236 0.000000e+00, 2.497431e-02, 4.979464e-02, 7.430797e-02, /* 0-3 */ 237 9.836316e-02, 1.218119e-01, 1.445097e-01, 1.663165e-01, /* 4-7 */ 238 1.870979e-01, 2.067257e-01, 2.250791e-01, 2.420447e-01, /* 8-11 */ 239 2.575181e-01, 2.714038e-01, 2.836162e-01, 2.940800e-01, /* 12-15 */ 240 3.027307e-01, 3.095150e-01, 3.143910e-01, 3.173286e-01, /* 16-19 */ 241 3.183099e-01, 3.173286e-01, 3.143910e-01, 3.095150e-01, /* 20-23 */ 242 3.027307e-01, 2.940800e-01, 2.836162e-01, 2.714038e-01, /* 24-27 */ 243 2.575181e-01, 2.420447e-01, 2.250791e-01, 2.067257e-01, /* 28-31 */ 244 1.870979e-01, 1.663165e-01, 1.445097e-01, 1.218119e-01, /* 32-35 */ 245 9.836316e-02, 7.430797e-02, 4.979464e-02, 2.497431e-02, /* 36-39 */ 246 -0.000000e+00, -2.497431e-02, -4.979464e-02, -7.430797e-02, /* 40-43 */ 247 -9.836316e-02, -1.218119e-01, -1.445097e-01, -1.663165e-01, /* 44-47 */ 248 -1.870979e-01, -2.067257e-01, -2.250791e-01, -2.420447e-01, /* 48-51 */ 249 -2.575181e-01, -2.714038e-01, -2.836162e-01, -2.940800e-01, /* 52-55 */ 250 -3.027307e-01, -3.095150e-01, -3.143910e-01, -3.173286e-01, /* 56-59 */ 251 -3.183099e-01, -3.173286e-01, -3.143910e-01, -3.095150e-01, /* 60-63 */ 252 -3.027307e-01, -2.940800e-01, -2.836162e-01, -2.714038e-01, /* 64-67 */ 253 -2.575181e-01, -2.420447e-01, -2.250791e-01, -2.067257e-01, /* 68-71 */ 254 -1.870979e-01, -1.663165e-01, -1.445097e-01, -1.218119e-01, /* 72-75 */ 255 -9.836316e-02, -7.430797e-02, -4.979464e-02, -2.497431e-02, /* 76-79 */ 256 0.000000e+00}; /* 80 */ 257 258 /* 259 * Decoder operations at the end of each second are driven by a state 260 * machine. The transition matrix consists of a dispatch table indexed 261 * by second number. Each entry in the table contains a case switch 262 * number and argument. 263 */ 264 struct progx { 265 int sw; /* case switch number */ 266 int arg; /* argument */ 267 }; 268 269 /* 270 * Case switch numbers 271 */ 272 #define IDLE 0 /* no operation */ 273 #define COEF 1 /* BCD bit */ 274 #define COEF2 2 /* BCD bit ignored */ 275 #define DECIM9 3 /* BCD digit 0-9 */ 276 #define DECIM6 4 /* BCD digit 0-6 */ 277 #define DECIM3 5 /* BCD digit 0-3 */ 278 #define DECIM2 6 /* BCD digit 0-2 */ 279 #define MSCBIT 7 /* miscellaneous bit */ 280 #define MSC20 8 /* miscellaneous bit */ 281 #define MSC21 9 /* QSY probe channel */ 282 #define MIN1 10 /* minute */ 283 #define MIN2 11 /* leap second */ 284 #define SYNC2 12 /* QSY data channel */ 285 #define SYNC3 13 /* QSY data channel */ 286 287 /* 288 * Offsets in decoding matrix 289 */ 290 #define MN 0 /* minute digits (2) */ 291 #define HR 2 /* hour digits (2) */ 292 #define DA 4 /* day digits (3) */ 293 #define YR 7 /* year digits (2) */ 294 295 struct progx progx[] = { 296 {SYNC2, 0}, /* 0 latch sync max */ 297 {SYNC3, 0}, /* 1 QSY data channel */ 298 {MSCBIT, DST2}, /* 2 dst2 */ 299 {MSCBIT, SECWAR}, /* 3 lw */ 300 {COEF, 0}, /* 4 1 year units */ 301 {COEF, 1}, /* 5 2 */ 302 {COEF, 2}, /* 6 4 */ 303 {COEF, 3}, /* 7 8 */ 304 {DECIM9, YR}, /* 8 */ 305 {IDLE, 0}, /* 9 p1 */ 306 {COEF, 0}, /* 10 1 minute units */ 307 {COEF, 1}, /* 11 2 */ 308 {COEF, 2}, /* 12 4 */ 309 {COEF, 3}, /* 13 8 */ 310 {DECIM9, MN}, /* 14 */ 311 {COEF, 0}, /* 15 10 minute tens */ 312 {COEF, 1}, /* 16 20 */ 313 {COEF, 2}, /* 17 40 */ 314 {COEF2, 3}, /* 18 80 (not used) */ 315 {DECIM6, MN + 1}, /* 19 p2 */ 316 {COEF, 0}, /* 20 1 hour units */ 317 {COEF, 1}, /* 21 2 */ 318 {COEF, 2}, /* 22 4 */ 319 {COEF, 3}, /* 23 8 */ 320 {DECIM9, HR}, /* 24 */ 321 {COEF, 0}, /* 25 10 hour tens */ 322 {COEF, 1}, /* 26 20 */ 323 {COEF2, 2}, /* 27 40 (not used) */ 324 {COEF2, 3}, /* 28 80 (not used) */ 325 {DECIM2, HR + 1}, /* 29 p3 */ 326 {COEF, 0}, /* 30 1 day units */ 327 {COEF, 1}, /* 31 2 */ 328 {COEF, 2}, /* 32 4 */ 329 {COEF, 3}, /* 33 8 */ 330 {DECIM9, DA}, /* 34 */ 331 {COEF, 0}, /* 35 10 day tens */ 332 {COEF, 1}, /* 36 20 */ 333 {COEF, 2}, /* 37 40 */ 334 {COEF, 3}, /* 38 80 */ 335 {DECIM9, DA + 1}, /* 39 p4 */ 336 {COEF, 0}, /* 40 100 day hundreds */ 337 {COEF, 1}, /* 41 200 */ 338 {COEF2, 2}, /* 42 400 (not used) */ 339 {COEF2, 3}, /* 43 800 (not used) */ 340 {DECIM3, DA + 2}, /* 44 */ 341 {IDLE, 0}, /* 45 */ 342 {IDLE, 0}, /* 46 */ 343 {IDLE, 0}, /* 47 */ 344 {IDLE, 0}, /* 48 */ 345 {IDLE, 0}, /* 49 p5 */ 346 {MSCBIT, DUTS}, /* 50 dut+- */ 347 {COEF, 0}, /* 51 10 year tens */ 348 {COEF, 1}, /* 52 20 */ 349 {COEF, 2}, /* 53 40 */ 350 {COEF, 3}, /* 54 80 */ 351 {MSC20, DST1}, /* 55 dst1 */ 352 {MSCBIT, DUT1}, /* 56 0.1 dut */ 353 {MSCBIT, DUT2}, /* 57 0.2 */ 354 {MSC21, DUT4}, /* 58 0.4 QSY probe channel */ 355 {MIN1, 0}, /* 59 p6 latch sync min */ 356 {MIN2, 0} /* 60 leap second */ 357 }; 358 359 /* 360 * BCD coefficients for maximum likelihood digit decode 361 */ 362 #define P15 1. /* max positive number */ 363 #define N15 -1. /* max negative number */ 364 365 /* 366 * Digits 0-9 367 */ 368 #define P9 (P15 / 4) /* mark (+1) */ 369 #define N9 (N15 / 4) /* space (-1) */ 370 371 double bcd9[][4] = { 372 {N9, N9, N9, N9}, /* 0 */ 373 {P9, N9, N9, N9}, /* 1 */ 374 {N9, P9, N9, N9}, /* 2 */ 375 {P9, P9, N9, N9}, /* 3 */ 376 {N9, N9, P9, N9}, /* 4 */ 377 {P9, N9, P9, N9}, /* 5 */ 378 {N9, P9, P9, N9}, /* 6 */ 379 {P9, P9, P9, N9}, /* 7 */ 380 {N9, N9, N9, P9}, /* 8 */ 381 {P9, N9, N9, P9}, /* 9 */ 382 {0, 0, 0, 0} /* backstop */ 383 }; 384 385 /* 386 * Digits 0-6 (minute tens) 387 */ 388 #define P6 (P15 / 3) /* mark (+1) */ 389 #define N6 (N15 / 3) /* space (-1) */ 390 391 double bcd6[][4] = { 392 {N6, N6, N6, 0}, /* 0 */ 393 {P6, N6, N6, 0}, /* 1 */ 394 {N6, P6, N6, 0}, /* 2 */ 395 {P6, P6, N6, 0}, /* 3 */ 396 {N6, N6, P6, 0}, /* 4 */ 397 {P6, N6, P6, 0}, /* 5 */ 398 {N6, P6, P6, 0}, /* 6 */ 399 {0, 0, 0, 0} /* backstop */ 400 }; 401 402 /* 403 * Digits 0-3 (day hundreds) 404 */ 405 #define P3 (P15 / 2) /* mark (+1) */ 406 #define N3 (N15 / 2) /* space (-1) */ 407 408 double bcd3[][4] = { 409 {N3, N3, 0, 0}, /* 0 */ 410 {P3, N3, 0, 0}, /* 1 */ 411 {N3, P3, 0, 0}, /* 2 */ 412 {P3, P3, 0, 0}, /* 3 */ 413 {0, 0, 0, 0} /* backstop */ 414 }; 415 416 /* 417 * Digits 0-2 (hour tens) 418 */ 419 #define P2 (P15 / 2) /* mark (+1) */ 420 #define N2 (N15 / 2) /* space (-1) */ 421 422 double bcd2[][4] = { 423 {N2, N2, 0, 0}, /* 0 */ 424 {P2, N2, 0, 0}, /* 1 */ 425 {N2, P2, 0, 0}, /* 2 */ 426 {0, 0, 0, 0} /* backstop */ 427 }; 428 429 /* 430 * DST decode (DST2 DST1) for prettyprint 431 */ 432 char dstcod[] = { 433 'S', /* 00 standard time */ 434 'I', /* 01 set clock ahead at 0200 local */ 435 'O', /* 10 set clock back at 0200 local */ 436 'D' /* 11 daylight time */ 437 }; 438 439 /* 440 * The decoding matrix consists of nine row vectors, one for each digit 441 * of the timecode. The digits are stored from least to most significant 442 * order. The maximum likelihood timecode is formed from the digits 443 * corresponding to the maximum likelihood values reading in the 444 * opposite order: yy ddd hh:mm. 445 */ 446 struct decvec { 447 int radix; /* radix (3, 4, 6, 10) */ 448 int digit; /* current clock digit */ 449 int mldigit; /* maximum likelihood digit */ 450 int phase; /* maximum likelihood digit phase */ 451 int count; /* match count */ 452 double digprb; /* max digit probability */ 453 double digsnr; /* likelihood function (dB) */ 454 double like[10]; /* likelihood integrator 0-9 */ 455 }; 456 457 /* 458 * The station structure is used to acquire the minute pulse from WWV 459 * and/or WWVH. These stations are distinguished by the frequency used 460 * for the second and minute sync pulses, 1000 Hz for WWV and 1200 Hz 461 * for WWVH. Other than frequency, the format is the same. 462 */ 463 struct sync { 464 double epoch; /* accumulated epoch differences */ 465 double maxamp; /* sync max envelope (square) */ 466 double noiamp; /* sync noise envelope (square) */ 467 long pos; /* max amplitude position */ 468 long lastpos; /* last max position */ 469 long mepoch; /* minute synch epoch */ 470 471 double amp; /* sync amplitude (I, Q squares) */ 472 double synamp; /* sync max envelope at 800 ms */ 473 double synmax; /* sync envelope at 0 s */ 474 double synmin; /* sync envelope at 59, 1 s */ 475 double synsnr; /* sync signal SNR */ 476 int count; /* bit counter */ 477 char refid[5]; /* reference identifier */ 478 int select; /* select bits */ 479 int reach; /* reachability register */ 480 }; 481 482 /* 483 * The channel structure is used to mitigate between channels. 484 */ 485 struct chan { 486 int gain; /* audio gain */ 487 double sigamp; /* data max envelope (square) */ 488 double noiamp; /* data noise envelope (square) */ 489 double datsnr; /* data signal SNR */ 490 struct sync wwv; /* wwv station */ 491 struct sync wwvh; /* wwvh station */ 492 }; 493 494 /* 495 * WWV unit control structure 496 */ 497 struct wwvunit { 498 l_fp timestamp; /* audio sample timestamp */ 499 l_fp tick; /* audio sample increment */ 500 double phase, freq; /* logical clock phase and frequency */ 501 double monitor; /* audio monitor point */ 502 int fd_icom; /* ICOM file descriptor */ 503 int errflg; /* error flags */ 504 int watch; /* watchcat */ 505 506 /* 507 * Audio codec variables 508 */ 509 double comp[SIZE]; /* decompanding table */ 510 int port; /* codec port */ 511 int gain; /* codec gain */ 512 int mongain; /* codec monitor gain */ 513 int clipcnt; /* sample clipped count */ 514 #ifdef IRIG_SUCKS 515 l_fp wigwag; /* wiggle accumulator */ 516 int wp; /* wiggle filter pointer */ 517 l_fp wiggle[WIGGLE]; /* wiggle filter */ 518 l_fp wigbot[WIGGLE]; /* wiggle bottom fisher*/ 519 #endif /* IRIG_SUCKS */ 520 521 /* 522 * Variables used to establish basic system timing 523 */ 524 int avgint; /* master time constant */ 525 int tepoch; /* sync epoch median */ 526 int yepoch; /* sync epoch */ 527 int repoch; /* buffered sync epoch */ 528 double epomax; /* second sync amplitude */ 529 double eposnr; /* second sync SNR */ 530 double irig; /* data I channel amplitude */ 531 double qrig; /* data Q channel amplitude */ 532 int datapt; /* 100 Hz ramp */ 533 double datpha; /* 100 Hz VFO control */ 534 int rphase; /* second sample counter */ 535 long mphase; /* minute sample counter */ 536 537 /* 538 * Variables used to mitigate which channel to use 539 */ 540 struct chan mitig[NCHAN]; /* channel data */ 541 struct sync *sptr; /* station pointer */ 542 int dchan; /* data channel */ 543 int schan; /* probe channel */ 544 int achan; /* active channel */ 545 546 /* 547 * Variables used by the clock state machine 548 */ 549 struct decvec decvec[9]; /* decoding matrix */ 550 int rsec; /* seconds counter */ 551 int digcnt; /* count of digits synchronized */ 552 553 /* 554 * Variables used to estimate signal levels and bit/digit 555 * probabilities 556 */ 557 double sigsig; /* data max signal */ 558 double sigamp; /* data max envelope (square) */ 559 double noiamp; /* data noise envelope (square) */ 560 double datsnr; /* data SNR (dB) */ 561 562 /* 563 * Variables used to establish status and alarm conditions 564 */ 565 int status; /* status bits */ 566 int alarm; /* alarm flashers */ 567 int misc; /* miscellaneous timecode bits */ 568 int errcnt; /* data bit error counter */ 569 int errbit; /* data bit errors in minute */ 570 }; 571 572 /* 573 * Function prototypes 574 */ 575 static int wwv_start P((int, struct peer *)); 576 static void wwv_shutdown P((int, struct peer *)); 577 static void wwv_receive P((struct recvbuf *)); 578 static void wwv_poll P((int, struct peer *)); 579 580 /* 581 * More function prototypes 582 */ 583 static void wwv_epoch P((struct peer *)); 584 static void wwv_rf P((struct peer *, double)); 585 static void wwv_endpoc P((struct peer *, int)); 586 static void wwv_rsec P((struct peer *, double)); 587 static void wwv_qrz P((struct peer *, struct sync *, 588 double, int)); 589 static void wwv_corr4 P((struct peer *, struct decvec *, 590 double [], double [][4])); 591 static void wwv_gain P((struct peer *)); 592 static void wwv_tsec P((struct wwvunit *)); 593 static double wwv_data P((struct wwvunit *, double)); 594 static int timecode P((struct wwvunit *, char *)); 595 static double wwv_snr P((double, double)); 596 static int carry P((struct decvec *)); 597 static void wwv_newchan P((struct peer *)); 598 static void wwv_newgame P((struct peer *)); 599 static double wwv_metric P((struct sync *)); 600 #ifdef ICOM 601 static int wwv_qsy P((struct peer *, int)); 602 #endif /* ICOM */ 603 604 static double qsy[NCHAN] = {2.5, 5, 10, 15, 20}; /* frequencies (MHz) */ 605 606 /* 607 * Transfer vector 608 */ 609 struct refclock refclock_wwv = { 610 wwv_start, /* start up driver */ 611 wwv_shutdown, /* shut down driver */ 612 wwv_poll, /* transmit poll message */ 613 noentry, /* not used (old wwv_control) */ 614 noentry, /* initialize driver (not used) */ 615 noentry, /* not used (old wwv_buginfo) */ 616 NOFLAGS /* not used */ 617 }; 618 619 620 /* 621 * wwv_start - open the devices and initialize data for processing 622 */ 623 static int 624 wwv_start( 625 int unit, /* instance number (used by PCM) */ 626 struct peer *peer /* peer structure pointer */ 627 ) 628 { 629 struct refclockproc *pp; 630 struct wwvunit *up; 631 #ifdef ICOM 632 int temp; 633 #endif /* ICOM */ 634 635 /* 636 * Local variables 637 */ 638 int fd; /* file descriptor */ 639 int i; /* index */ 640 double step; /* codec adjustment */ 641 642 /* 643 * Open audio device 644 */ 645 fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit); 646 if (fd < 0) 647 return (0); 648 #ifdef DEBUG 649 if (debug) 650 audio_show(); 651 #endif 652 653 /* 654 * Allocate and initialize unit structure 655 */ 656 if (!(up = (struct wwvunit *)emalloc(sizeof(struct wwvunit)))) { 657 close(fd); 658 return (0); 659 } 660 memset(up, 0, sizeof(struct wwvunit)); 661 pp = peer->procptr; 662 pp->unitptr = (caddr_t)up; 663 pp->io.clock_recv = wwv_receive; 664 pp->io.srcclock = (caddr_t)peer; 665 pp->io.datalen = 0; 666 pp->io.fd = fd; 667 if (!io_addclock(&pp->io)) { 668 close(fd); 669 free(up); 670 return (0); 671 } 672 673 /* 674 * Initialize miscellaneous variables 675 */ 676 peer->precision = PRECISION; 677 pp->clockdesc = DESCRIPTION; 678 679 /* 680 * The companded samples are encoded sign-magnitude. The table 681 * contains all the 256 values in the interest of speed. 682 */ 683 up->comp[0] = up->comp[OFFSET] = 0.; 684 up->comp[1] = 1; up->comp[OFFSET + 1] = -1.; 685 up->comp[2] = 3; up->comp[OFFSET + 2] = -3.; 686 step = 2.; 687 for (i = 3; i < OFFSET; i++) { 688 up->comp[i] = up->comp[i - 1] + step; 689 up->comp[OFFSET + i] = -up->comp[i]; 690 if (i % 16 == 0) 691 step *= 2.; 692 } 693 DTOLFP(1. / SECOND, &up->tick); 694 695 /* 696 * Initialize the decoding matrix with the radix for each digit 697 * position. 698 */ 699 up->decvec[MN].radix = 10; /* minutes */ 700 up->decvec[MN + 1].radix = 6; 701 up->decvec[HR].radix = 10; /* hours */ 702 up->decvec[HR + 1].radix = 3; 703 up->decvec[DA].radix = 10; /* days */ 704 up->decvec[DA + 1].radix = 10; 705 up->decvec[DA + 2].radix = 4; 706 up->decvec[YR].radix = 10; /* years */ 707 up->decvec[YR + 1].radix = 10; 708 wwv_newgame(peer); 709 up->schan = up->achan = 3; 710 711 /* 712 * Initialize autotune if available. Start out at 15 MHz. Note 713 * that the ICOM select code must be less than 128, so the high 714 * order bit can be used to select the line speed. 715 */ 716 #ifdef ICOM 717 temp = 0; 718 #ifdef DEBUG 719 if (debug > 1) 720 temp = P_TRACE; 721 #endif 722 if (peer->ttl != 0) { 723 if (peer->ttl & 0x80) 724 up->fd_icom = icom_init("/dev/icom", B1200, 725 temp); 726 else 727 up->fd_icom = icom_init("/dev/icom", B9600, 728 temp); 729 } 730 if (up->fd_icom > 0) { 731 if ((temp = wwv_qsy(peer, up->schan)) != 0) { 732 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT) 733 msyslog(LOG_NOTICE, 734 "icom: radio not found"); 735 up->errflg = CEVNT_FAULT; 736 close(up->fd_icom); 737 up->fd_icom = 0; 738 } else { 739 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT) 740 msyslog(LOG_NOTICE, 741 "icom: autotune enabled"); 742 } 743 } 744 #endif /* ICOM */ 745 return (1); 746 } 747 748 749 /* 750 * wwv_shutdown - shut down the clock 751 */ 752 static void 753 wwv_shutdown( 754 int unit, /* instance number (not used) */ 755 struct peer *peer /* peer structure pointer */ 756 ) 757 { 758 struct refclockproc *pp; 759 struct wwvunit *up; 760 761 pp = peer->procptr; 762 up = (struct wwvunit *)pp->unitptr; 763 io_closeclock(&pp->io); 764 if (up->fd_icom > 0) 765 close(up->fd_icom); 766 free(up); 767 } 768 769 770 /* 771 * wwv_receive - receive data from the audio device 772 * 773 * This routine reads input samples and adjusts the logical clock to 774 * track the A/D sample clock by dropping or duplicating codec samples. 775 * It also controls the A/D signal level with an AGC loop to mimimize 776 * quantization noise and avoid overload. 777 */ 778 static void 779 wwv_receive( 780 struct recvbuf *rbufp /* receive buffer structure pointer */ 781 ) 782 { 783 struct peer *peer; 784 struct refclockproc *pp; 785 struct wwvunit *up; 786 787 /* 788 * Local variables 789 */ 790 double sample; /* codec sample */ 791 u_char *dpt; /* buffer pointer */ 792 int bufcnt; /* buffer counter */ 793 l_fp ltemp; 794 795 peer = (struct peer *)rbufp->recv_srcclock; 796 pp = peer->procptr; 797 up = (struct wwvunit *)pp->unitptr; 798 799 /* 800 * Main loop - read until there ain't no more. Note codec 801 * samples are bit-inverted. 802 */ 803 DTOLFP((double)rbufp->recv_length / SECOND, <emp); 804 L_SUB(&rbufp->recv_time, <emp); 805 up->timestamp = rbufp->recv_time; 806 dpt = rbufp->recv_buffer; 807 for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) { 808 sample = up->comp[~*dpt++ & 0xff]; 809 810 /* 811 * Clip noise spikes greater than MAXSIG. If no clips, 812 * increase the gain a tad; if the clips are too high, 813 * decrease a tad. 814 */ 815 if (sample > MAXSIG) { 816 sample = MAXSIG; 817 up->clipcnt++; 818 } else if (sample < -MAXSIG) { 819 sample = -MAXSIG; 820 up->clipcnt++; 821 } 822 823 /* 824 * Variable frequency oscillator. The codec oscillator 825 * runs at the nominal rate of 8000 samples per second, 826 * or 125 us per sample. A frequency change of one unit 827 * results in either duplicating or deleting one sample 828 * per second, which results in a frequency change of 829 * 125 PPM. 830 */ 831 up->phase += up->freq / SECOND; 832 if (up->phase >= .5) { 833 up->phase -= 1.; 834 } else if (up->phase < -.5) { 835 up->phase += 1.; 836 wwv_rf(peer, sample); 837 wwv_rf(peer, sample); 838 } else { 839 wwv_rf(peer, sample); 840 } 841 L_ADD(&up->timestamp, &up->tick); 842 } 843 844 /* 845 * Set the input port and monitor gain for the next buffer. 846 */ 847 if (pp->sloppyclockflag & CLK_FLAG2) 848 up->port = 2; 849 else 850 up->port = 1; 851 if (pp->sloppyclockflag & CLK_FLAG3) 852 up->mongain = MONGAIN; 853 else 854 up->mongain = 0; 855 } 856 857 858 /* 859 * wwv_poll - called by the transmit procedure 860 * 861 * This routine keeps track of status. If no offset samples have been 862 * processed during a poll interval, a timeout event is declared. If 863 * errors have have occurred during the interval, they are reported as 864 * well. Once the clock is set, it always appears reachable, unless 865 * reset by watchdog timeout. 866 */ 867 static void 868 wwv_poll( 869 int unit, /* instance number (not used) */ 870 struct peer *peer /* peer structure pointer */ 871 ) 872 { 873 struct refclockproc *pp; 874 struct wwvunit *up; 875 876 pp = peer->procptr; 877 up = (struct wwvunit *)pp->unitptr; 878 if (pp->coderecv == pp->codeproc) 879 up->errflg = CEVNT_TIMEOUT; 880 if (up->errflg) 881 refclock_report(peer, up->errflg); 882 up->errflg = 0; 883 pp->polls++; 884 } 885 886 887 /* 888 * wwv_rf - process signals and demodulate to baseband 889 * 890 * This routine grooms and filters decompanded raw audio samples. The 891 * output signals include the 100-Hz baseband data signal in quadrature 892 * form, plus the epoch index of the second sync signal and the second 893 * index of the minute sync signal. 894 * 895 * There are two 1-s ramps used by this program. Both count the 8000 896 * logical clock samples spanning exactly one second. The epoch ramp 897 * counts the samples starting at an arbitrary time. The rphase ramp 898 * counts the samples starting at the 5-ms second sync pulse found 899 * during the epoch ramp. 900 * 901 * There are two 1-m ramps used by this program. The mphase ramp counts 902 * the 480,000 logical clock samples spanning exactly one minute and 903 * starting at an arbitrary time. The rsec ramp counts the 60 seconds of 904 * the minute starting at the 800-ms minute sync pulse found during the 905 * mphase ramp. The rsec ramp drives the seconds state machine to 906 * determine the bits and digits of the timecode. 907 * 908 * Demodulation operations are based on three synthesized quadrature 909 * sinusoids: 100 Hz for the data signal, 1000 Hz for the WWV sync 910 * signal and 1200 Hz for the WWVH sync signal. These drive synchronous 911 * matched filters for the data signal (170 ms at 100 Hz), WWV minute 912 * sync signal (800 ms at 1000 Hz) and WWVH minute sync signal (800 ms 913 * at 1200 Hz). Two additional matched filters are switched in 914 * as required for the WWV second sync signal (5 ms at 1000 Hz) and 915 * WWVH second sync signal (5 ms at 1200 Hz). 916 */ 917 static void 918 wwv_rf( 919 struct peer *peer, /* peerstructure pointer */ 920 double isig /* input signal */ 921 ) 922 { 923 struct refclockproc *pp; 924 struct wwvunit *up; 925 struct sync *sp; 926 927 static double lpf[5]; /* 150-Hz lpf delay line */ 928 double data; /* lpf output */ 929 static double bpf[9]; /* 1000/1200-Hz bpf delay line */ 930 double syncx; /* bpf output */ 931 static double mf[41]; /* 1000/1200-Hz mf delay line */ 932 double mfsync; /* mf output */ 933 934 static int iptr; /* data channel pointer */ 935 static double ibuf[DATSIZ]; /* data I channel delay line */ 936 static double qbuf[DATSIZ]; /* data Q channel delay line */ 937 938 static int jptr; /* sync channel pointer */ 939 static double cibuf[SYNSIZ]; /* wwv I channel delay line */ 940 static double cqbuf[SYNSIZ]; /* wwv Q channel delay line */ 941 static double ciamp; /* wwv I channel amplitude */ 942 static double cqamp; /* wwv Q channel amplitude */ 943 static int csinptr; /* wwv channel phase */ 944 static double hibuf[SYNSIZ]; /* wwvh I channel delay line */ 945 static double hqbuf[SYNSIZ]; /* wwvh Q channel delay line */ 946 static double hiamp; /* wwvh I channel amplitude */ 947 static double hqamp; /* wwvh Q channel amplitude */ 948 static int hsinptr; /* wwvh channels phase */ 949 950 static double epobuf[SECOND]; /* epoch sync comb filter */ 951 static double epomax; /* epoch sync amplitude buffer */ 952 static int epopos; /* epoch sync position buffer */ 953 954 static int iniflg; /* initialization flag */ 955 int epoch; /* comb filter index */ 956 int pdelay; /* propagation delay (samples) */ 957 double dtemp; 958 int i; 959 960 pp = peer->procptr; 961 up = (struct wwvunit *)pp->unitptr; 962 963 if (!iniflg) { 964 iniflg = 1; 965 memset((char *)lpf, 0, sizeof(lpf)); 966 memset((char *)bpf, 0, sizeof(bpf)); 967 memset((char *)mf, 0, sizeof(mf)); 968 memset((char *)ibuf, 0, sizeof(ibuf)); 969 memset((char *)qbuf, 0, sizeof(qbuf)); 970 memset((char *)cibuf, 0, sizeof(cibuf)); 971 memset((char *)cqbuf, 0, sizeof(cqbuf)); 972 memset((char *)hibuf, 0, sizeof(hibuf)); 973 memset((char *)hqbuf, 0, sizeof(hqbuf)); 974 memset((char *)epobuf, 0, sizeof(epobuf)); 975 } 976 977 /* 978 * Baseband data demodulation. The 100-Hz subcarrier is 979 * extracted using a 150-Hz IIR lowpass filter. This attenuates 980 * the 1000/1200-Hz sync signals, as well as the 440-Hz and 981 * 600-Hz tones and most of the noise and voice modulation 982 * components. 983 * 984 * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB 985 * passband ripple, -50 dB stopband ripple. 986 */ 987 data = (lpf[4] = lpf[3]) * 8.360961e-01; 988 data += (lpf[3] = lpf[2]) * -3.481740e+00; 989 data += (lpf[2] = lpf[1]) * 5.452988e+00; 990 data += (lpf[1] = lpf[0]) * -3.807229e+00; 991 lpf[0] = isig - data; 992 data = lpf[0] * 3.281435e-03 993 + lpf[1] * -1.149947e-02 994 + lpf[2] * 1.654858e-02 995 + lpf[3] * -1.149947e-02 996 + lpf[4] * 3.281435e-03; 997 998 /* 999 * The I and Q quadrature data signals are produced by 1000 * multiplying the filtered signal by 100-Hz sine and cosine 1001 * signals, respectively. The data signals are demodulated by 1002 * 170-ms synchronous matched filters to produce the amplitude 1003 * and phase signals used by the decoder. 1004 */ 1005 i = up->datapt; 1006 up->datapt = (up->datapt + IN100) % 80; 1007 dtemp = sintab[i] * data / DATSIZ * DGAIN; 1008 up->irig -= ibuf[iptr]; 1009 ibuf[iptr] = dtemp; 1010 up->irig += dtemp; 1011 i = (i + 20) % 80; 1012 dtemp = sintab[i] * data / DATSIZ * DGAIN; 1013 up->qrig -= qbuf[iptr]; 1014 qbuf[iptr] = dtemp; 1015 up->qrig += dtemp; 1016 iptr = (iptr + 1) % DATSIZ; 1017 1018 /* 1019 * Baseband sync demodulation. The 1000/1200 sync signals are 1020 * extracted using a 600-Hz IIR bandpass filter. This removes 1021 * the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz 1022 * tones and most of the noise and voice modulation components. 1023 * 1024 * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB 1025 * passband ripple, -50 dB stopband ripple. 1026 */ 1027 syncx = (bpf[8] = bpf[7]) * 4.897278e-01; 1028 syncx += (bpf[7] = bpf[6]) * -2.765914e+00; 1029 syncx += (bpf[6] = bpf[5]) * 8.110921e+00; 1030 syncx += (bpf[5] = bpf[4]) * -1.517732e+01; 1031 syncx += (bpf[4] = bpf[3]) * 1.975197e+01; 1032 syncx += (bpf[3] = bpf[2]) * -1.814365e+01; 1033 syncx += (bpf[2] = bpf[1]) * 1.159783e+01; 1034 syncx += (bpf[1] = bpf[0]) * -4.735040e+00; 1035 bpf[0] = isig - syncx; 1036 syncx = bpf[0] * 8.203628e-03 1037 + bpf[1] * -2.375732e-02 1038 + bpf[2] * 3.353214e-02 1039 + bpf[3] * -4.080258e-02 1040 + bpf[4] * 4.605479e-02 1041 + bpf[5] * -4.080258e-02 1042 + bpf[6] * 3.353214e-02 1043 + bpf[7] * -2.375732e-02 1044 + bpf[8] * 8.203628e-03; 1045 1046 /* 1047 * The I and Q quadrature minute sync signals are produced by 1048 * multiplying the filtered signal by 1000-Hz (WWV) and 1200-Hz 1049 * (WWVH) sine and cosine signals, respectively. The resulting 1050 * signals are demodulated by 800-ms synchronous matched filters 1051 * to synchronize the second and minute and to detect which one 1052 * (or both) the WWV or WWVH signal is present. 1053 * 1054 * Note the master timing ramps, which run continuously. The 1055 * minute counter (mphase) counts the samples in the minute, 1056 * while the second counter (epoch) counts the samples in the 1057 * second. 1058 */ 1059 up->mphase = (up->mphase + 1) % MINUTE; 1060 epoch = up->mphase % SECOND; 1061 i = csinptr; 1062 csinptr = (csinptr + IN1000) % 80; 1063 dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; 1064 ciamp = ciamp - cibuf[jptr] + dtemp; 1065 cibuf[jptr] = dtemp; 1066 i = (i + 20) % 80; 1067 dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; 1068 cqamp = cqamp - cqbuf[jptr] + dtemp; 1069 cqbuf[jptr] = dtemp; 1070 sp = &up->mitig[up->schan].wwv; 1071 dtemp = ciamp * ciamp + cqamp * cqamp; 1072 sp->amp = dtemp; 1073 if (!(up->status & MSYNC)) 1074 wwv_qrz(peer, sp, dtemp, (int)(pp->fudgetime1 * 1075 SECOND)); 1076 i = hsinptr; 1077 hsinptr = (hsinptr + IN1200) % 80; 1078 dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; 1079 hiamp = hiamp - hibuf[jptr] + dtemp; 1080 hibuf[jptr] = dtemp; 1081 i = (i + 20) % 80; 1082 dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; 1083 hqamp = hqamp - hqbuf[jptr] + dtemp; 1084 hqbuf[jptr] = dtemp; 1085 sp = &up->mitig[up->schan].wwvh; 1086 dtemp = hiamp * hiamp + hqamp * hqamp; 1087 sp->amp = dtemp; 1088 if (!(up->status & MSYNC)) 1089 wwv_qrz(peer, sp, dtemp, (int)(pp->fudgetime2 * 1090 SECOND)); 1091 jptr = (jptr + 1) % SYNSIZ; 1092 1093 /* 1094 * The following section is called once per minute. It does 1095 * housekeeping and timeout functions and empties the dustbins. 1096 */ 1097 if (up->mphase == 0) { 1098 up->watch++; 1099 if (!(up->status & MSYNC)) { 1100 1101 /* 1102 * If minute sync has not been acquired before 1103 * timeout, or if no signal is heard, the 1104 * program cycles to the next frequency and 1105 * tries again. 1106 */ 1107 wwv_newchan(peer); 1108 if (!(up->status & (SELV | SELH)) || up->watch > 1109 ACQSN) { 1110 wwv_newgame(peer); 1111 #ifdef ICOM 1112 if (up->fd_icom > 0) { 1113 up->schan = (up->schan + 1) % 1114 NCHAN; 1115 wwv_qsy(peer, up->schan); 1116 } 1117 #endif /* ICOM */ 1118 } 1119 } else { 1120 1121 /* 1122 * If the leap bit is set, set the minute epoch 1123 * back one second so the station processes 1124 * don't miss a beat. 1125 */ 1126 if (up->status & LEPSEC) { 1127 up->mphase -= SECOND; 1128 if (up->mphase < 0) 1129 up->mphase += MINUTE; 1130 } 1131 } 1132 } 1133 1134 /* 1135 * When the channel metric reaches threshold and the second 1136 * counter matches the minute epoch within the second, the 1137 * driver has synchronized to the station. The second number is 1138 * the remaining seconds until the next minute epoch, while the 1139 * sync epoch is zero. Watch out for the first second; if 1140 * already synchronized to the second, the buffered sync epoch 1141 * must be set. 1142 */ 1143 if (up->status & MSYNC) { 1144 wwv_epoch(peer); 1145 } else if ((sp = up->sptr) != NULL) { 1146 struct chan *cp; 1147 1148 if (sp->count >= AMIN && epoch == sp->mepoch % SECOND) { 1149 up->rsec = 60 - sp->mepoch / SECOND; 1150 up->rphase = 0; 1151 up->status |= MSYNC; 1152 up->watch = 0; 1153 if (!(up->status & SSYNC)) 1154 up->repoch = up->yepoch = epoch; 1155 else 1156 up->repoch = up->yepoch; 1157 for (i = 0; i < NCHAN; i++) { 1158 cp = &up->mitig[i]; 1159 cp->wwv.count = cp->wwv.reach = 0; 1160 cp->wwvh.count = cp->wwvh.reach = 0; 1161 } 1162 } 1163 } 1164 1165 /* 1166 * The second sync pulse is extracted using 5-ms (40 sample) FIR 1167 * matched filters at 1000 Hz for WWV or 1200 Hz for WWVH. This 1168 * pulse is used for the most precise synchronization, since if 1169 * provides a resolution of one sample (125 us). The filters run 1170 * only if the station has been reliably determined. 1171 */ 1172 if (up->status & SELV) { 1173 pdelay = (int)(pp->fudgetime1 * SECOND); 1174 1175 /* 1176 * WWV FIR matched filter, five cycles of 1000-Hz 1177 * sinewave. 1178 */ 1179 mf[40] = mf[39]; 1180 mfsync = (mf[39] = mf[38]) * 4.224514e-02; 1181 mfsync += (mf[38] = mf[37]) * 5.974365e-02; 1182 mfsync += (mf[37] = mf[36]) * 4.224514e-02; 1183 mf[36] = mf[35]; 1184 mfsync += (mf[35] = mf[34]) * -4.224514e-02; 1185 mfsync += (mf[34] = mf[33]) * -5.974365e-02; 1186 mfsync += (mf[33] = mf[32]) * -4.224514e-02; 1187 mf[32] = mf[31]; 1188 mfsync += (mf[31] = mf[30]) * 4.224514e-02; 1189 mfsync += (mf[30] = mf[29]) * 5.974365e-02; 1190 mfsync += (mf[29] = mf[28]) * 4.224514e-02; 1191 mf[28] = mf[27]; 1192 mfsync += (mf[27] = mf[26]) * -4.224514e-02; 1193 mfsync += (mf[26] = mf[25]) * -5.974365e-02; 1194 mfsync += (mf[25] = mf[24]) * -4.224514e-02; 1195 mf[24] = mf[23]; 1196 mfsync += (mf[23] = mf[22]) * 4.224514e-02; 1197 mfsync += (mf[22] = mf[21]) * 5.974365e-02; 1198 mfsync += (mf[21] = mf[20]) * 4.224514e-02; 1199 mf[20] = mf[19]; 1200 mfsync += (mf[19] = mf[18]) * -4.224514e-02; 1201 mfsync += (mf[18] = mf[17]) * -5.974365e-02; 1202 mfsync += (mf[17] = mf[16]) * -4.224514e-02; 1203 mf[16] = mf[15]; 1204 mfsync += (mf[15] = mf[14]) * 4.224514e-02; 1205 mfsync += (mf[14] = mf[13]) * 5.974365e-02; 1206 mfsync += (mf[13] = mf[12]) * 4.224514e-02; 1207 mf[12] = mf[11]; 1208 mfsync += (mf[11] = mf[10]) * -4.224514e-02; 1209 mfsync += (mf[10] = mf[9]) * -5.974365e-02; 1210 mfsync += (mf[9] = mf[8]) * -4.224514e-02; 1211 mf[8] = mf[7]; 1212 mfsync += (mf[7] = mf[6]) * 4.224514e-02; 1213 mfsync += (mf[6] = mf[5]) * 5.974365e-02; 1214 mfsync += (mf[5] = mf[4]) * 4.224514e-02; 1215 mf[4] = mf[3]; 1216 mfsync += (mf[3] = mf[2]) * -4.224514e-02; 1217 mfsync += (mf[2] = mf[1]) * -5.974365e-02; 1218 mfsync += (mf[1] = mf[0]) * -4.224514e-02; 1219 mf[0] = syncx; 1220 } else if (up->status & SELH) { 1221 pdelay = (int)(pp->fudgetime2 * SECOND); 1222 1223 /* 1224 * WWVH FIR matched filter, six cycles of 1200-Hz 1225 * sinewave. 1226 */ 1227 mf[40] = mf[39]; 1228 mfsync = (mf[39] = mf[38]) * 4.833363e-02; 1229 mfsync += (mf[38] = mf[37]) * 5.681959e-02; 1230 mfsync += (mf[37] = mf[36]) * 1.846180e-02; 1231 mfsync += (mf[36] = mf[35]) * -3.511644e-02; 1232 mfsync += (mf[35] = mf[34]) * -5.974365e-02; 1233 mfsync += (mf[34] = mf[33]) * -3.511644e-02; 1234 mfsync += (mf[33] = mf[32]) * 1.846180e-02; 1235 mfsync += (mf[32] = mf[31]) * 5.681959e-02; 1236 mfsync += (mf[31] = mf[30]) * 4.833363e-02; 1237 mf[30] = mf[29]; 1238 mfsync += (mf[29] = mf[28]) * -4.833363e-02; 1239 mfsync += (mf[28] = mf[27]) * -5.681959e-02; 1240 mfsync += (mf[27] = mf[26]) * -1.846180e-02; 1241 mfsync += (mf[26] = mf[25]) * 3.511644e-02; 1242 mfsync += (mf[25] = mf[24]) * 5.974365e-02; 1243 mfsync += (mf[24] = mf[23]) * 3.511644e-02; 1244 mfsync += (mf[23] = mf[22]) * -1.846180e-02; 1245 mfsync += (mf[22] = mf[21]) * -5.681959e-02; 1246 mfsync += (mf[21] = mf[20]) * -4.833363e-02; 1247 mf[20] = mf[19]; 1248 mfsync += (mf[19] = mf[18]) * 4.833363e-02; 1249 mfsync += (mf[18] = mf[17]) * 5.681959e-02; 1250 mfsync += (mf[17] = mf[16]) * 1.846180e-02; 1251 mfsync += (mf[16] = mf[15]) * -3.511644e-02; 1252 mfsync += (mf[15] = mf[14]) * -5.974365e-02; 1253 mfsync += (mf[14] = mf[13]) * -3.511644e-02; 1254 mfsync += (mf[13] = mf[12]) * 1.846180e-02; 1255 mfsync += (mf[12] = mf[11]) * 5.681959e-02; 1256 mfsync += (mf[11] = mf[10]) * 4.833363e-02; 1257 mf[10] = mf[9]; 1258 mfsync += (mf[9] = mf[8]) * -4.833363e-02; 1259 mfsync += (mf[8] = mf[7]) * -5.681959e-02; 1260 mfsync += (mf[7] = mf[6]) * -1.846180e-02; 1261 mfsync += (mf[6] = mf[5]) * 3.511644e-02; 1262 mfsync += (mf[5] = mf[4]) * 5.974365e-02; 1263 mfsync += (mf[4] = mf[3]) * 3.511644e-02; 1264 mfsync += (mf[3] = mf[2]) * -1.846180e-02; 1265 mfsync += (mf[2] = mf[1]) * -5.681959e-02; 1266 mfsync += (mf[1] = mf[0]) * -4.833363e-02; 1267 mf[0] = syncx; 1268 } else { 1269 mfsync = 0; 1270 pdelay = 0; 1271 } 1272 1273 /* 1274 * Enhance the seconds sync pulse using a 1-s (8000-sample) comb 1275 * filter. Correct for the FIR matched filter delay, which is 5 1276 * ms for both the WWV and WWVH filters, and also for the 1277 * propagation delay. Once each second look for second sync. If 1278 * not in minute sync, fiddle the codec gain. Note the SNR is 1279 * computed from the maximum sample and the two samples 6 ms 1280 * before and 6 ms after it, so if we slip more than a cycle the 1281 * SNR should plummet. 1282 */ 1283 dtemp = (epobuf[epoch] += (mfsync - epobuf[epoch]) / 1284 up->avgint); 1285 if (dtemp > epomax) { 1286 epomax = dtemp; 1287 epopos = epoch; 1288 } 1289 if (epoch == 0) { 1290 int k, j; 1291 1292 up->epomax = epomax; 1293 k = epopos - 6 * MS; 1294 if (k < 0) 1295 k += SECOND; 1296 j = epopos + 6 * MS; 1297 if (j >= SECOND) 1298 i -= SECOND; 1299 up->eposnr = wwv_snr(epomax, max(abs(epobuf[k]), 1300 abs(epobuf[j]))); 1301 epopos -= pdelay + 5 * MS; 1302 if (epopos < 0) 1303 epopos += SECOND; 1304 wwv_endpoc(peer, epopos); 1305 if (!(up->status & SSYNC)) 1306 up->alarm |= SYNERR; 1307 epomax = 0; 1308 if (!(up->status & MSYNC)) 1309 wwv_gain(peer); 1310 } 1311 } 1312 1313 1314 /* 1315 * wwv_qrz - identify and acquire WWV/WWVH minute sync pulse 1316 * 1317 * This routine implements a virtual station process used to acquire 1318 * minute sync and to mitigate among the ten frequency and station 1319 * combinations. During minute sync acquisition the process probes each 1320 * frequency in turn for the minute pulse from either station, which 1321 * involves searching through the entire minute of samples. After 1322 * finding a candidate, the process searches only the seconds before and 1323 * after the candidate for the signal and all other seconds for the 1324 * noise. 1325 * 1326 * Students of radar receiver technology will discover this algorithm 1327 * amounts to a range gate discriminator. The discriminator requires 1328 * that the peak minute pulse amplitude be at least 2000 and the SNR be 1329 * at least 6 dB. In addition after finding a candidate, The peak second 1330 * pulse amplitude must be at least 2000, the SNR at least 6 dB and the 1331 * difference between the current and previous epoch must be less than 1332 * 7.5 ms, which corresponds to a frequency error of 125 PPM.. A compare 1333 * counter keeps track of the number of successive intervals which 1334 * satisfy these criteria. 1335 * 1336 * Note that, while the minute pulse is found by by the discriminator, 1337 * the actual value is determined from the second epoch. The assumption 1338 * is that the discriminator peak occurs about 800 ms into the second, 1339 * so the timing is retarted to the previous second epoch. 1340 */ 1341 static void 1342 wwv_qrz( 1343 struct peer *peer, /* peer structure pointer */ 1344 struct sync *sp, /* sync channel structure */ 1345 double syncx, /* bandpass filtered sync signal */ 1346 int pdelay /* propagation delay (samples) */ 1347 ) 1348 { 1349 struct refclockproc *pp; 1350 struct wwvunit *up; 1351 char tbuf[80]; /* monitor buffer */ 1352 double snr; /* on-pulse/off-pulse ratio (dB) */ 1353 long epoch, fpoch; 1354 int isgood; 1355 1356 pp = peer->procptr; 1357 up = (struct wwvunit *)pp->unitptr; 1358 1359 /* 1360 * Find the sample with peak energy, which defines the minute 1361 * epoch. If a sample has been found with good amplitude, 1362 * accumulate the noise squares for all except the second before 1363 * and after that position. 1364 */ 1365 isgood = up->epomax > STHR && up->eposnr > SSNR; 1366 if (isgood) { 1367 fpoch = up->mphase % SECOND - up->tepoch; 1368 if (fpoch < 0) 1369 fpoch += SECOND; 1370 } else { 1371 fpoch = pdelay + SYNSIZ; 1372 } 1373 epoch = up->mphase - fpoch; 1374 if (epoch < 0) 1375 epoch += MINUTE; 1376 if (syncx > sp->maxamp) { 1377 sp->maxamp = syncx; 1378 sp->pos = epoch; 1379 } 1380 if (abs((epoch - sp->lastpos) % MINUTE) > SECOND) 1381 sp->noiamp += syncx; 1382 1383 /* 1384 * At the end of the minute, determine the epoch of the 1385 * sync pulse, as well as the SNR and difference between 1386 * the current and previous epoch, which represents the 1387 * intrinsic frequency error plus jitter. 1388 */ 1389 if (up->mphase == 0) { 1390 sp->synmax = sqrt(sp->maxamp); 1391 sp->synmin = sqrt(sp->noiamp / (MINUTE - 2 * SECOND)); 1392 epoch = (sp->pos - sp->lastpos) % MINUTE; 1393 1394 /* 1395 * If not yet in minute sync, we have to do a little 1396 * dance to find a valid minute sync pulse, emphasis 1397 * valid. 1398 */ 1399 snr = wwv_snr(sp->synmax, sp->synmin); 1400 isgood = isgood && sp->synmax > ATHR && snr > ASNR; 1401 switch (sp->count) { 1402 1403 /* 1404 * In state 0 the station was not heard during the 1405 * previous probe. Look for the biggest blip greater 1406 * than the amplitude threshold in the minute and assume 1407 * that the minute sync pulse. We're fishing here, since 1408 * the range gate has not yet been determined. If found, 1409 * bump to state 1. 1410 */ 1411 case 0: 1412 if (sp->synmax >= ATHR) 1413 sp->count++; 1414 break; 1415 1416 /* 1417 * In state 1 a candidate blip has been found and the 1418 * next minute has been searched for another blip. If 1419 * none are found acceptable, drop back to state 0 and 1420 * hunt some more. Otherwise, a legitimate minute pulse 1421 * may have been found, so bump to state 2. 1422 */ 1423 case 1: 1424 if (!isgood) { 1425 sp->count = 0; 1426 break; 1427 } 1428 sp->count++; 1429 break; 1430 1431 /* 1432 * In states 2 and above, continue to groom samples as 1433 * before and drop back to state 0 if the groom fails. 1434 * If it succeeds, set the epoch and bump to the next 1435 * state until reaching the threshold, if ever. 1436 */ 1437 default: 1438 if (!isgood || abs(epoch) > AWND * MS) { 1439 sp->count = 0; 1440 break; 1441 } 1442 sp->mepoch = sp->pos; 1443 sp->count++; 1444 break; 1445 } 1446 if (pp->sloppyclockflag & CLK_FLAG4) { 1447 sprintf(tbuf, 1448 "wwv8 %d %3d %s %d %5.0f %5.1f %5ld %5d %ld", 1449 up->port, up->gain, sp->refid, sp->count, 1450 sp->synmax, snr, sp->pos, up->tepoch, 1451 epoch); 1452 record_clock_stats(&peer->srcadr, tbuf); 1453 #ifdef DEBUG 1454 if (debug) 1455 printf("%s\n", tbuf); 1456 #endif 1457 } 1458 sp->lastpos = sp->pos; 1459 sp->maxamp = sp->noiamp = 0; 1460 } 1461 } 1462 1463 1464 /* 1465 * wwv_endpoc - identify and acquire second sync pulse 1466 * 1467 * This routine is called at the end of the second sync interval. It 1468 * determines the second sync epoch position within the interval and 1469 * disciplines the sample clock using a frequency-lock loop (FLL). 1470 * 1471 * Second sync is determined in the RF input routine as the maximum 1472 * over all 8000 samples in the second comb filter. To assure accurate 1473 * and reliable time and frequency discipline, this routine performs a 1474 * great deal of heavy-handed heuristic data filtering and grooming. 1475 * 1476 * Note that, since the minute sync pulse is very wide (800 ms), precise 1477 * minute sync epoch acquisition requires at least a rough estimate of 1478 * the second sync pulse (5 ms). This becomes more important in choppy 1479 * conditions at the lower frequencies at night, since sferics and 1480 * cochannel crude can badly distort the minute pulse. 1481 */ 1482 static void 1483 wwv_endpoc( 1484 struct peer *peer, /* peer structure pointer */ 1485 int epopos /* epoch max position */ 1486 ) 1487 { 1488 struct refclockproc *pp; 1489 struct wwvunit *up; 1490 static int epoch_mf[3]; /* epoch median filter */ 1491 static int xepoch; /* last second epoch */ 1492 static int zepoch; /* last averaging interval epoch */ 1493 static int syncnt; /* run length counter */ 1494 static int maxrun; /* longest run length */ 1495 static int mepoch; /* longest run epoch */ 1496 static int avgcnt; /* averaging interval counter */ 1497 static int avginc; /* averaging ratchet */ 1498 static int iniflg; /* initialization flag */ 1499 char tbuf[80]; /* monitor buffer */ 1500 double dtemp; 1501 int tmp2; 1502 1503 pp = peer->procptr; 1504 up = (struct wwvunit *)pp->unitptr; 1505 if (!iniflg) { 1506 iniflg = 1; 1507 memset((char *)epoch_mf, 0, sizeof(epoch_mf)); 1508 } 1509 1510 /* 1511 * A three-stage median filter is used to help denoise the 1512 * second sync pulse. The median sample becomes the candidate 1513 * epoch. 1514 */ 1515 epoch_mf[2] = epoch_mf[1]; 1516 epoch_mf[1] = epoch_mf[0]; 1517 epoch_mf[0] = epopos; 1518 if (epoch_mf[0] > epoch_mf[1]) { 1519 if (epoch_mf[1] > epoch_mf[2]) 1520 up->tepoch = epoch_mf[1]; /* 0 1 2 */ 1521 else if (epoch_mf[2] > epoch_mf[0]) 1522 up->tepoch = epoch_mf[0]; /* 2 0 1 */ 1523 else 1524 up->tepoch = epoch_mf[2]; /* 0 2 1 */ 1525 } else { 1526 if (epoch_mf[1] < epoch_mf[2]) 1527 up->tepoch = epoch_mf[1]; /* 2 1 0 */ 1528 else if (epoch_mf[2] < epoch_mf[0]) 1529 up->tepoch = epoch_mf[0]; /* 1 0 2 */ 1530 else 1531 up->tepoch = epoch_mf[2]; /* 1 2 0 */ 1532 } 1533 1534 /* 1535 * If the signal amplitude or SNR fall below thresholds or if no 1536 * stations are heard, dim the second sync lamp and start over. 1537 */ 1538 if (!(up->status & (SELV | SELH)) || up->epomax < STHR || 1539 up->eposnr < SSNR) { 1540 up->status &= ~(SSYNC | FGATE); 1541 avgcnt = syncnt = maxrun = 0; 1542 return; 1543 } 1544 avgcnt++; 1545 1546 /* 1547 * If the epoch candidate is the same as the last one, increment 1548 * the compare counter. If not, save the length and epoch of the 1549 * current run for use later and reset the counter. 1550 */ 1551 tmp2 = (up->tepoch - xepoch) % SECOND; 1552 if (tmp2 == 0) { 1553 syncnt++; 1554 } else { 1555 if (maxrun > 0 && mepoch == xepoch) { 1556 maxrun += syncnt; 1557 } else if (syncnt > maxrun) { 1558 maxrun = syncnt; 1559 mepoch = xepoch; 1560 } 1561 syncnt = 0; 1562 } 1563 if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status & (SSYNC | 1564 MSYNC))) { 1565 sprintf(tbuf, 1566 "wwv1 %04x %5.0f %5.1f %5d %5d %4d %4d", 1567 up->status, up->epomax, up->eposnr, up->tepoch, 1568 tmp2, avgcnt, syncnt); 1569 record_clock_stats(&peer->srcadr, tbuf); 1570 #ifdef DEBUG 1571 if (debug) 1572 printf("%s\n", tbuf); 1573 #endif /* DEBUG */ 1574 } 1575 1576 /* 1577 * The sample clock frequency is disciplined using a first order 1578 * feedback loop with time constant consistent with the Allan 1579 * intercept of typical computer clocks. 1580 * 1581 * The frequency update is calculated from the epoch change in 1582 * 125-us units divided by the averaging interval in seconds. 1583 * The averaging interval affects other receiver functions, 1584 * including the the 1000/1200-Hz comb filter and codec clock 1585 * loop. It also affects the 100-Hz subcarrier loop and the bit 1586 * and digit comparison counter thresholds. 1587 */ 1588 if (avgcnt < up->avgint) { 1589 xepoch = up->tepoch; 1590 return; 1591 } 1592 1593 /* 1594 * During the averaging interval the longest run of identical 1595 * epoches is determined. If the longest run is at least 10 1596 * seconds, the SSYNC bit is lit and the value becomes the 1597 * reference epoch for the next interval. If not, the second 1598 * synd lamp is dark and flashers set. 1599 */ 1600 if (maxrun > 0 && mepoch == xepoch) { 1601 maxrun += syncnt; 1602 } else if (syncnt > maxrun) { 1603 maxrun = syncnt; 1604 mepoch = xepoch; 1605 } 1606 xepoch = up->tepoch; 1607 if (maxrun > SCMP) { 1608 up->status |= SSYNC; 1609 up->yepoch = mepoch; 1610 } else { 1611 up->status &= ~SSYNC; 1612 } 1613 1614 /* 1615 * If the epoch change over the averaging interval is less than 1616 * 1 ms, the frequency is adjusted, but clamped at +-125 PPM. If 1617 * greater than 1 ms, the counter is decremented. If the epoch 1618 * change is less than 0.5 ms, the counter is incremented. If 1619 * the counter increments to +3, the averaging interval is 1620 * doubled and the counter set to zero; if it increments to -3, 1621 * the interval is halved and the counter set to zero. 1622 * 1623 * Here be spooks. From careful observations, the epoch 1624 * sometimes makes a long run of identical samples, then takes a 1625 * lurch due apparently to lost interrupts or spooks. If this 1626 * happens, the epoch change times the maximum run length will 1627 * be greater than the averaging interval, so the lurch should 1628 * be believed but the frequency left alone. Really intricate 1629 * here. 1630 */ 1631 if (maxrun == 0) 1632 mepoch = up->tepoch; 1633 dtemp = (mepoch - zepoch) % SECOND; 1634 if (up->status & FGATE) { 1635 if (abs(dtemp) < MAXFREQ * MINAVG) { 1636 if (maxrun * abs(mepoch - zepoch) < 1637 avgcnt) { 1638 up->freq += dtemp / avgcnt; 1639 if (up->freq > MAXFREQ) 1640 up->freq = MAXFREQ; 1641 else if (up->freq < -MAXFREQ) 1642 up->freq = -MAXFREQ; 1643 } 1644 if (abs(dtemp) < MAXFREQ * MINAVG / 2) { 1645 if (avginc < 3) { 1646 avginc++; 1647 } else { 1648 if (up->avgint < MAXAVG) { 1649 up->avgint <<= 1; 1650 avginc = 0; 1651 } 1652 } 1653 } 1654 } else { 1655 if (avginc > -3) { 1656 avginc--; 1657 } else { 1658 if (up->avgint > MINAVG) { 1659 up->avgint >>= 1; 1660 avginc = 0; 1661 } 1662 } 1663 } 1664 } 1665 if (pp->sloppyclockflag & CLK_FLAG4) { 1666 sprintf(tbuf, 1667 "wwv2 %04x %4.0f %4d %4d %2d %4d %4.0f %6.1f", 1668 up->status, up->epomax, mepoch, maxrun, avginc, 1669 avgcnt, dtemp, up->freq * 1e6 / SECOND); 1670 record_clock_stats(&peer->srcadr, tbuf); 1671 #ifdef DEBUG 1672 if (debug) 1673 printf("%s\n", tbuf); 1674 #endif /* DEBUG */ 1675 } 1676 up->status |= FGATE; 1677 zepoch = mepoch; 1678 avgcnt = syncnt = maxrun = 0; 1679 } 1680 1681 1682 /* 1683 * wwv_epoch - epoch scanner 1684 * 1685 * This routine scans the receiver second epoch to determine the signal 1686 * amplitudes and pulse timings. Receiver synchronization is determined 1687 * by the minute sync pulse detected in the wwv_rf() routine and the 1688 * second sync pulse detected in the wwv_epoch() routine. A pulse width 1689 * discriminator extracts data signals from the 100-Hz subcarrier. The 1690 * transmitted signals are delayed by the propagation delay, receiver 1691 * delay and filter delay of this program. Delay corrections are 1692 * introduced separately for WWV and WWVH. 1693 * 1694 * Most communications radios use a highpass filter in the audio stages, 1695 * which can do nasty things to the subcarrier phase relative to the 1696 * sync pulses. Therefore, the data subcarrier reference phase is 1697 * disciplined using the hardlimited quadrature-phase signal sampled at 1698 * the same time as the in-phase signal. The phase tracking loop uses 1699 * phase adjustments of plus-minus one sample (125 us). 1700 */ 1701 static void 1702 wwv_epoch( 1703 struct peer *peer /* peer structure pointer */ 1704 ) 1705 { 1706 struct refclockproc *pp; 1707 struct wwvunit *up; 1708 struct chan *cp; 1709 static double dpulse; /* data pulse length */ 1710 double dtemp; 1711 1712 pp = peer->procptr; 1713 up = (struct wwvunit *)pp->unitptr; 1714 1715 /* 1716 * Sample the minute sync pulse envelopes at epoch 800 for both 1717 * the WWV and WWVH stations. This will be used later for 1718 * channel and station mitigation. Note that the seconds epoch 1719 * is set here well before the end of the second to make sure we 1720 * never seet the epoch backwards. 1721 */ 1722 if (up->rphase == 800 * MS) { 1723 up->repoch = up->yepoch; 1724 cp = &up->mitig[up->achan]; 1725 cp->wwv.synamp = cp->wwv.amp; 1726 cp->wwvh.synamp = cp->wwvh.amp; 1727 } 1728 1729 /* 1730 * Sample the data subcarrier at epoch 15 ms, giving a guard 1731 * time of +-15 ms from the beginning of the second until the 1732 * pulse rises at 30 ms. The I-channel amplitude is used to 1733 * calculate the slice level. The envelope amplitude is used 1734 * during the probe seconds to determine the SNR. There is a 1735 * compromise here; we want to delay the sample as long as 1736 * possible to give the radio time to change frequency and the 1737 * AGC to stabilize, but as early as possible if the second 1738 * epoch is not exact. 1739 */ 1740 if (up->rphase == 15 * MS) { 1741 up->noiamp = up->irig * up->irig + up->qrig * up->qrig; 1742 1743 /* 1744 * Sample the data subcarrier at epoch 215 ms, giving a guard 1745 * time of +-15 ms from the earliest the pulse peak can be 1746 * reached to the earliest it can begin to fall. For the data 1747 * channel latch the I-channel amplitude for all except the 1748 * probe seconds and adjust the 100-Hz reference oscillator 1749 * phase using the Q-channel amplitude at this epoch. For the 1750 * probe channel latch the envelope amplitude. 1751 */ 1752 } else if (up->rphase == 215 * MS) { 1753 up->sigsig = up->irig; 1754 if (up->sigsig < 0) 1755 up->sigsig = 0; 1756 up->datpha = up->qrig / up->avgint; 1757 if (up->datpha >= 0) { 1758 up->datapt++; 1759 if (up->datapt >= 80) 1760 up->datapt -= 80; 1761 } else { 1762 up->datapt--; 1763 if (up->datapt < 0) 1764 up->datapt += 80; 1765 } 1766 up->sigamp = up->irig * up->irig + up->qrig * up->qrig; 1767 1768 /* 1769 * The slice level is set half way between the peak signal and 1770 * noise levels. Sample the negative zero crossing after epoch 1771 * 200 ms and record the epoch at that time. This defines the 1772 * length of the data pulse, which will later be converted into 1773 * scaled bit probabilities. 1774 */ 1775 } else if (up->rphase > 200 * MS) { 1776 dtemp = (up->sigsig + sqrt(up->noiamp)) / 2; 1777 if (up->irig < dtemp && dpulse == 0) 1778 dpulse = up->rphase; 1779 } 1780 1781 /* 1782 * At the end of the second crank the clock state machine and 1783 * adjust the codec gain. Note the epoch is buffered from the 1784 * center of the second in order to avoid jitter while the 1785 * seconds synch is diddling the epoch. Then, determine the true 1786 * offset and update the median filter in the driver interface. 1787 * 1788 * Sample the data subcarrier envelope at the end of the second 1789 * to determine the SNR for the pulse. This gives a guard time 1790 * of +-30 ms from the decay of the longest pulse to the rise of 1791 * the next pulse. 1792 */ 1793 up->rphase++; 1794 if (up->mphase % SECOND == up->repoch) { 1795 up->datsnr = wwv_snr(up->sigsig, sqrt(up->noiamp)); 1796 wwv_rsec(peer, dpulse); 1797 wwv_gain(peer); 1798 up->rphase = dpulse = 0; 1799 } 1800 } 1801 1802 1803 /* 1804 * wwv_rsec - process receiver second 1805 * 1806 * This routine is called at the end of each receiver second to 1807 * implement the per-second state machine. The machine assembles BCD 1808 * digit bits, decodes miscellaneous bits and dances the leap seconds. 1809 * 1810 * Normally, the minute has 60 seconds numbered 0-59. If the leap 1811 * warning bit is set, the last minute (1439) of 30 June (day 181 or 182 1812 * for leap years) or 31 December (day 365 or 366 for leap years) is 1813 * augmented by one second numbered 60. This is accomplished by 1814 * extending the minute interval by one second and teaching the state 1815 * machine to ignore it. 1816 */ 1817 static void 1818 wwv_rsec( 1819 struct peer *peer, /* peer structure pointer */ 1820 double dpulse 1821 ) 1822 { 1823 static int iniflg; /* initialization flag */ 1824 static double bcddld[4]; /* BCD data bits */ 1825 static double bitvec[61]; /* bit integrator for misc bits */ 1826 struct refclockproc *pp; 1827 struct wwvunit *up; 1828 struct chan *cp; 1829 struct sync *sp, *rp; 1830 l_fp offset; /* offset in NTP seconds */ 1831 double bit; /* bit likelihood */ 1832 char tbuf[80]; /* monitor buffer */ 1833 int sw, arg, nsec; 1834 #ifdef IRIG_SUCKS 1835 int i; 1836 l_fp ltemp; 1837 #endif /* IRIG_SUCKS */ 1838 1839 pp = peer->procptr; 1840 up = (struct wwvunit *)pp->unitptr; 1841 if (!iniflg) { 1842 iniflg = 1; 1843 memset((char *)bitvec, 0, sizeof(bitvec)); 1844 } 1845 1846 /* 1847 * The bit represents the probability of a hit on zero (negative 1848 * values), a hit on one (positive values) or a miss (zero 1849 * value). The likelihood vector is the exponential average of 1850 * these probabilities. Only the bits of this vector 1851 * corresponding to the miscellaneous bits of the timecode are 1852 * used, but it's easier to do them all. After that, crank the 1853 * seconds state machine. 1854 */ 1855 nsec = up->rsec + 1; 1856 bit = wwv_data(up, dpulse); 1857 bitvec[up->rsec] += (bit - bitvec[up->rsec]) / TCONST; 1858 sw = progx[up->rsec].sw; 1859 arg = progx[up->rsec].arg; 1860 switch (sw) { 1861 1862 /* 1863 * Ignore this second. 1864 */ 1865 case IDLE: /* 9, 45-49 */ 1866 break; 1867 1868 /* 1869 * Probe channel stuff 1870 * 1871 * The WWV/H format contains data pulses in second 59 (position 1872 * identifier), second 1 (not used) and the minute sync pulse in 1873 * second 0. At the end of second 58, QSY to the probe channel, 1874 * which rotates over all WWV/H frequencies. At the end of 1875 * second 1 QSY back to the data channel. 1876 * 1877 * At the end of second 0 save the minute sync pulse peak value 1878 * previously latched at 800 ms. 1879 */ 1880 case SYNC2: /* 0 */ 1881 cp = &up->mitig[up->achan]; 1882 cp->wwv.synmax = sqrt(cp->wwv.synamp); 1883 cp->wwvh.synmax = sqrt(cp->wwvh.synamp); 1884 break; 1885 1886 /* 1887 * At the end of second 1 determine the minute sync pulse 1888 * amplitude and SNR and set SYNCNG if these values are below 1889 * thresholds. Determine the data pulse amplitude and SNR and 1890 * set DATANG if these values are below thresholds. Set ERRRNG 1891 * if data pulses in second 59 and second 1 are decoded in 1892 * error. Shift a 1 into the reachability register if SYNCNG and 1893 * DATANG are both lit; otherwise shift a 0. Ignore ERRRNG for 1894 * the present. The number of 1 bits in the last six intervals 1895 * represents the channel metric used by the mitigation routine. 1896 * Finally, QSY back to the data channel. 1897 */ 1898 case SYNC3: /* 1 */ 1899 cp = &up->mitig[up->achan]; 1900 cp->sigamp = sqrt(up->sigamp); 1901 cp->noiamp = sqrt(up->noiamp); 1902 cp->datsnr = wwv_snr(cp->sigamp, cp->noiamp); 1903 1904 /* 1905 * WWV station 1906 */ 1907 sp = &cp->wwv; 1908 sp->synmin = sqrt((sp->synmin + sp->synamp) / 2.); 1909 sp->synsnr = wwv_snr(sp->synmax, sp->synmin); 1910 sp->select &= ~(SYNCNG | DATANG | ERRRNG); 1911 if (sp->synmax < QTHR || sp->synsnr < QSNR) 1912 sp->select |= SYNCNG; 1913 if (cp->sigamp < XTHR || cp->datsnr < XSNR) 1914 sp->select |= DATANG; 1915 if (up->errcnt > 2) 1916 sp->select |= ERRRNG; 1917 sp->reach <<= 1; 1918 if (sp->reach & (1 << AMAX)) 1919 sp->count--; 1920 if (!(sp->select & (SYNCNG | DATANG))) { 1921 sp->reach |= 1; 1922 sp->count++; 1923 } 1924 1925 /* 1926 * WWVH station 1927 */ 1928 rp = &cp->wwvh; 1929 rp->synmin = sqrt((rp->synmin + rp->synamp) / 2.); 1930 rp->synsnr = wwv_snr(rp->synmax, rp->synmin); 1931 rp->select &= ~(SYNCNG | DATANG | ERRRNG); 1932 if (rp->synmax < QTHR || rp->synsnr < QSNR) 1933 rp->select |= SYNCNG; 1934 if (cp->sigamp < XTHR || cp->datsnr < XSNR) 1935 rp->select |= DATANG; 1936 if (up->errcnt > 2) 1937 rp->select |= ERRRNG; 1938 rp->reach <<= 1; 1939 if (rp->reach & (1 << AMAX)) 1940 rp->count--; 1941 if (!(rp->select & (SYNCNG | DATANG | ERRRNG))) { 1942 rp->reach |= 1; 1943 rp->count++; 1944 } 1945 1946 /* 1947 * Set up for next minute. 1948 */ 1949 if (pp->sloppyclockflag & CLK_FLAG4) { 1950 sprintf(tbuf, 1951 "wwv5 %2d %04x %3d %4d %d %.0f/%.1f %s %04x %.0f %.0f/%.1f %s %04x %.0f %.0f/%.1f", 1952 up->port, up->status, up->gain, up->yepoch, 1953 up->errcnt, cp->sigamp, cp->datsnr, 1954 sp->refid, sp->reach & 0xffff, 1955 wwv_metric(sp), sp->synmax, sp->synsnr, 1956 rp->refid, rp->reach & 0xffff, 1957 wwv_metric(rp), rp->synmax, rp->synsnr); 1958 record_clock_stats(&peer->srcadr, tbuf); 1959 #ifdef DEBUG 1960 if (debug) 1961 printf("%s\n", tbuf); 1962 #endif /* DEBUG */ 1963 } 1964 #ifdef ICOM 1965 if (up->fd_icom > 0) 1966 wwv_qsy(peer, up->dchan); 1967 #endif /* ICOM */ 1968 up->status &= ~SFLAG; 1969 up->errcnt = 0; 1970 up->alarm = 0; 1971 wwv_newchan(peer); 1972 break; 1973 1974 /* 1975 * Save the bit probability in the BCD data vector at the index 1976 * given by the argument. Note that all bits of the vector have 1977 * to be above the data gate threshold for the digit to be 1978 * considered valid. Bits not used in the digit are forced to 1979 * zero and not checked for errors. 1980 */ 1981 case COEF: /* 4-7, 10-13, 15-17, 20-23, 1982 25-26, 30-33, 35-38, 40-41, 1983 51-54 */ 1984 if (up->status & DGATE) 1985 up->status |= BGATE; 1986 bcddld[arg] = bit; 1987 break; 1988 1989 case COEF2: /* 18, 27-28, 42-43 */ 1990 bcddld[arg] = 0; 1991 break; 1992 1993 /* 1994 * Correlate coefficient vector with each valid digit vector and 1995 * save in decoding matrix. We step through the decoding matrix 1996 * digits correlating each with the coefficients and saving the 1997 * greatest and the next lower for later SNR calculation. 1998 */ 1999 case DECIM2: /* 29 */ 2000 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd2); 2001 break; 2002 2003 case DECIM3: /* 44 */ 2004 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd3); 2005 break; 2006 2007 case DECIM6: /* 19 */ 2008 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd6); 2009 break; 2010 2011 case DECIM9: /* 8, 14, 24, 34, 39 */ 2012 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd9); 2013 break; 2014 2015 /* 2016 * Miscellaneous bits. If above the positive threshold, declare 2017 * 1; if below the negative threshold, declare 0; otherwise 2018 * raise the SYMERR alarm. At the end of second 58, QSY to the 2019 * probe channel. The design is intended to preserve the bits 2020 * over periods of signal loss. 2021 */ 2022 case MSC20: /* 55 */ 2023 wwv_corr4(peer, &up->decvec[YR + 1], bcddld, bcd9); 2024 /* fall through */ 2025 2026 case MSCBIT: /* 2-3, 50, 56-57 */ 2027 if (bitvec[up->rsec] > BTHR) 2028 up->misc |= arg; 2029 else if (bitvec[up->rsec] < -BTHR) 2030 up->misc &= ~arg; 2031 else 2032 up->alarm |= SYMERR; 2033 break; 2034 2035 /* 2036 * Save the data channel gain, then QSY to the probe channel. 2037 */ 2038 case MSC21: /* 58 */ 2039 if (bitvec[up->rsec] > BTHR) 2040 up->misc |= arg; 2041 else if (bitvec[up->rsec] < -BTHR) 2042 up->misc &= ~arg; 2043 else 2044 up->alarm |= SYMERR; 2045 up->mitig[up->dchan].gain = up->gain; 2046 #ifdef ICOM 2047 if (up->fd_icom > 0) { 2048 up->schan = (up->schan + 1) % NCHAN; 2049 wwv_qsy(peer, up->schan); 2050 } 2051 #endif /* ICOM */ 2052 up->status |= SFLAG | SELV | SELH; 2053 up->errbit = up->errcnt; 2054 up->errcnt = 0; 2055 break; 2056 2057 /* 2058 * The endgames 2059 * 2060 * During second 59 the receiver and codec AGC are settling 2061 * down, so the data pulse is unusable. At the end of this 2062 * second, latch the minute sync pulse noise floor. Then do the 2063 * minute processing and update the system clock. If a leap 2064 * second sail on to the next second (60); otherwise, set up for 2065 * the next minute. 2066 */ 2067 case MIN1: /* 59 */ 2068 cp = &up->mitig[up->achan]; 2069 cp->wwv.synmin = cp->wwv.synamp; 2070 cp->wwvh.synmin = cp->wwvh.synamp; 2071 2072 /* 2073 * Dance the leap if necessary and the kernel has the 2074 * right stuff. Then, wind up the clock and initialize 2075 * for the following minute. If the leap dance, note the 2076 * kernel is armed one second before the actual leap is 2077 * scheduled. 2078 */ 2079 if (up->status & SSYNC && up->digcnt >= 9) 2080 up->status |= INSYNC; 2081 if (up->status & LEPDAY) { 2082 pp->leap = LEAP_ADDSECOND; 2083 } else { 2084 pp->leap = LEAP_NOWARNING; 2085 wwv_tsec(up); 2086 nsec = up->digcnt = 0; 2087 } 2088 pp->lencode = timecode(up, pp->a_lastcode); 2089 record_clock_stats(&peer->srcadr, pp->a_lastcode); 2090 #ifdef DEBUG 2091 if (debug) 2092 printf("wwv: timecode %d %s\n", pp->lencode, 2093 pp->a_lastcode); 2094 #endif /* DEBUG */ 2095 if (up->status & INSYNC && up->watch < HOLD) 2096 refclock_receive(peer); 2097 break; 2098 2099 /* 2100 * If LEPDAY is set on the last minute of 30 June or 31 2101 * December, the LEPSEC bit is set. At the end of the minute in 2102 * which LEPSEC is set the transmitter and receiver insert an 2103 * extra second (60) in the timescale and the minute sync skips 2104 * a second. We only get to test this wrinkle at intervals of 2105 * about 18 months; the actual mileage may vary. 2106 */ 2107 case MIN2: /* 60 */ 2108 wwv_tsec(up); 2109 nsec = up->digcnt = 0; 2110 break; 2111 } 2112 2113 /* 2114 * If digit sync has not been acquired before timeout or if no 2115 * station has been heard, game over and restart from scratch. 2116 */ 2117 if (!(up->status & DSYNC) && (!(up->status & (SELV | SELH)) || 2118 up->watch > DIGIT)) { 2119 wwv_newgame(peer); 2120 return; 2121 } 2122 2123 /* 2124 * If no timestamps have been struck before timeout, game over 2125 * and restart from scratch. 2126 */ 2127 if (up->watch > PANIC) { 2128 wwv_newgame(peer); 2129 return; 2130 } 2131 pp->disp += AUDIO_PHI; 2132 up->rsec = nsec; 2133 2134 #ifdef IRIG_SUCKS 2135 /* 2136 * You really don't wanna know what comes down here. Leave it to 2137 * say Solaris 2.8 broke the nice clean audio stream, apparently 2138 * affected by a 5-ms sawtooth jitter. Sundown on Solaris. This 2139 * leaves a little twilight. 2140 * 2141 * The scheme involves differentiation, forward learning and 2142 * integration. The sawtooth has a period of 11 seconds. The 2143 * timestamp differences are integrated and subtracted from the 2144 * signal. 2145 */ 2146 ltemp = pp->lastrec; 2147 L_SUB(<emp, &pp->lastref); 2148 if (ltemp.l_f < 0) 2149 ltemp.l_i = -1; 2150 else 2151 ltemp.l_i = 0; 2152 pp->lastref = pp->lastrec; 2153 if (!L_ISNEG(<emp)) 2154 L_CLR(&up->wigwag); 2155 else 2156 L_ADD(&up->wigwag, <emp); 2157 L_SUB(&pp->lastrec, &up->wigwag); 2158 up->wiggle[up->wp] = ltemp; 2159 2160 /* 2161 * Bottom fisher. To understand this, you have to know about 2162 * velocity microphones and AM transmitters. No further 2163 * explanation is offered, as this is truly a black art. 2164 */ 2165 up->wigbot[up->wp] = pp->lastrec; 2166 for (i = 0; i < WIGGLE; i++) { 2167 if (i != up->wp) 2168 up->wigbot[i].l_ui++; 2169 L_SUB(&pp->lastrec, &up->wigbot[i]); 2170 if (L_ISNEG(&pp->lastrec)) 2171 L_ADD(&pp->lastrec, &up->wigbot[i]); 2172 else 2173 pp->lastrec = up->wigbot[i]; 2174 } 2175 up->wp++; 2176 up->wp %= WIGGLE; 2177 #endif /* IRIG_SUCKS */ 2178 2179 /* 2180 * If victory has been declared and seconds sync is lit, strike 2181 * a timestamp. It should not be a surprise, especially if the 2182 * radio is not tunable, that sometimes no stations are above 2183 * the noise and the reference ID set to NONE. 2184 */ 2185 if (up->status & INSYNC && up->status & SSYNC) { 2186 pp->second = up->rsec; 2187 pp->minute = up->decvec[MN].digit + up->decvec[MN + 2188 1].digit * 10; 2189 pp->hour = up->decvec[HR].digit + up->decvec[HR + 2190 1].digit * 10; 2191 pp->day = up->decvec[DA].digit + up->decvec[DA + 2192 1].digit * 10 + up->decvec[DA + 2].digit * 100; 2193 pp->year = up->decvec[YR].digit + up->decvec[YR + 2194 1].digit * 10; 2195 pp->year += 2000; 2196 L_CLR(&offset); 2197 if (!clocktime(pp->day, pp->hour, pp->minute, 2198 pp->second, GMT, up->timestamp.l_ui, 2199 &pp->yearstart, &offset.l_ui)) { 2200 up->errflg = CEVNT_BADTIME; 2201 } else { 2202 up->watch = 0; 2203 pp->disp = 0; 2204 refclock_process_offset(pp, offset, 2205 up->timestamp, PDELAY); 2206 } 2207 } 2208 if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status & 2209 DSYNC)) { 2210 sprintf(tbuf, 2211 "wwv3 %2d %04x %5.0f %5.1f %5.0f %5.1f %5.0f", 2212 up->rsec, up->status, up->epomax, up->eposnr, 2213 up->sigsig, up->datsnr, bit); 2214 record_clock_stats(&peer->srcadr, tbuf); 2215 #ifdef DEBUG 2216 if (debug) 2217 printf("%s\n", tbuf); 2218 #endif /* DEBUG */ 2219 } 2220 } 2221 2222 2223 /* 2224 * wwv_data - calculate bit probability 2225 * 2226 * This routine is called at the end of the receiver second to calculate 2227 * the probabilities that the previous second contained a zero (P0), one 2228 * (P1) or position indicator (P2) pulse. If not in sync or if the data 2229 * bit is bad, a bit error is declared and the probabilities are forced 2230 * to zero. Otherwise, the data gate is enabled and the probabilities 2231 * are calculated. Note that the data matched filter contributes half 2232 * the pulse width, or 85 ms. 2233 * 2234 * It's important to observe that there are three conditions to 2235 * determine: average to +1 (hit), average to -1 (miss) or average to 2236 * zero (erasure). The erasure condition results from insufficient 2237 * signal (noise) and has no bias toward either a hit or miss. 2238 */ 2239 static double 2240 wwv_data( 2241 struct wwvunit *up, /* driver unit pointer */ 2242 double pulse /* pulse length (sample units) */ 2243 ) 2244 { 2245 double p0, p1, p2; /* probabilities */ 2246 double dpulse; /* pulse length in ms */ 2247 2248 p0 = p1 = p2 = 0; 2249 dpulse = pulse - DATSIZ / 2; 2250 2251 /* 2252 * If no station is being tracked, if either the data amplitude 2253 * or SNR are below threshold or if the pulse length is less 2254 * than 170 ms, declare an erasure. 2255 */ 2256 if (!(up->status & (SELV | SELH)) || up->sigsig < DTHR || 2257 up->datsnr < DSNR || dpulse < DATSIZ) { 2258 up->status |= DGATE; 2259 up->errcnt++; 2260 if (up->errcnt > MAXERR) 2261 up->alarm |= MODERR; 2262 return (0); 2263 } 2264 2265 /* 2266 * The probability of P0 is one below 200 ms falling to zero at 2267 * 500 ms. The probability of P1 is zero below 200 ms rising to 2268 * one at 500 ms and falling to zero at 800 ms. The probability 2269 * of P2 is zero below 500 ms, rising to one above 800 ms. 2270 */ 2271 up->status &= ~DGATE; 2272 if (dpulse < (200 * MS)) { 2273 p0 = 1; 2274 } else if (dpulse < 500 * MS) { 2275 dpulse -= 200 * MS; 2276 p1 = dpulse / (300 * MS); 2277 p0 = 1 - p1; 2278 } else if (dpulse < 800 * MS) { 2279 dpulse -= 500 * MS; 2280 p2 = dpulse / (300 * MS); 2281 p1 = 1 - p2; 2282 } else { 2283 p2 = 1; 2284 } 2285 2286 /* 2287 * The ouput is a metric that ranges from -1 (P0), to +1 (P1) 2288 * scaled for convenience. An output of zero represents an 2289 * erasure, either because of a data error or pulse length 2290 * greater than 500 ms. At the moment, we don't use P2. 2291 */ 2292 return ((p1 - p0) * MAXSIG); 2293 } 2294 2295 2296 /* 2297 * wwv_corr4 - determine maximum likelihood digit 2298 * 2299 * This routine correlates the received digit vector with the BCD 2300 * coefficient vectors corresponding to all valid digits at the given 2301 * position in the decoding matrix. The maximum value corresponds to the 2302 * maximum likelihood digit, while the ratio of this value to the next 2303 * lower value determines the likelihood function. Note that, if the 2304 * digit is invalid, the likelihood vector is averaged toward a miss. 2305 */ 2306 static void 2307 wwv_corr4( 2308 struct peer *peer, /* peer unit pointer */ 2309 struct decvec *vp, /* decoding table pointer */ 2310 double data[], /* received data vector */ 2311 double tab[][4] /* correlation vector array */ 2312 ) 2313 { 2314 struct refclockproc *pp; 2315 struct wwvunit *up; 2316 2317 double topmax, nxtmax; /* metrics */ 2318 double acc; /* accumulator */ 2319 char tbuf[80]; /* monitor buffer */ 2320 int mldigit; /* max likelihood digit */ 2321 int diff; /* decoding difference */ 2322 int i, j; 2323 2324 pp = peer->procptr; 2325 up = (struct wwvunit *)pp->unitptr; 2326 2327 /* 2328 * Correlate digit vector with each BCD coefficient vector. If 2329 * any BCD digit bit is bad, consider all bits a miss. 2330 */ 2331 mldigit = 0; 2332 topmax = nxtmax = -MAXSIG; 2333 for (i = 0; tab[i][0] != 0; i++) { 2334 acc = 0; 2335 for (j = 0; j < 4; j++) { 2336 if (!(up->status & BGATE)) 2337 acc += data[j] * tab[i][j]; 2338 } 2339 acc = (vp->like[i] += (acc - vp->like[i]) / TCONST); 2340 if (acc > topmax) { 2341 nxtmax = topmax; 2342 topmax = acc; 2343 mldigit = i; 2344 } else if (acc > nxtmax) { 2345 nxtmax = acc; 2346 } 2347 } 2348 vp->mldigit = mldigit; 2349 vp->digprb = topmax; 2350 vp->digsnr = wwv_snr(topmax, nxtmax); 2351 2352 /* 2353 * The maximum likelihood digit is compared with the current 2354 * clock digit. The difference represents the decoding phase 2355 * error. If the clock is not yet synchronized, the phase error 2356 * is corrected even of the digit probability and likelihood are 2357 * below thresholds. This avoids lengthy averaging times should 2358 * a carry mistake occur. However, the digit is not declared 2359 * synchronized until these values are above thresholds and the 2360 * last five decoded values are identical. If the clock is 2361 * synchronized, the phase error is not corrected unless the 2362 * last five digits are all above thresholds and identical. This 2363 * avoids mistakes when the signal is coming out of the noise 2364 * and the SNR is very marginal. 2365 */ 2366 diff = mldigit - vp->digit; 2367 if (diff < 0) 2368 diff += vp->radix; 2369 if (diff != vp->phase) { 2370 vp->count = 0; 2371 vp->phase = diff; 2372 } 2373 if (vp->digsnr < BSNR) { 2374 vp->count = 0; 2375 up->alarm |= SYMERR; 2376 } else if (vp->digprb < BTHR) { 2377 vp->count = 0; 2378 up->alarm |= SYMERR; 2379 if (!(up->status & INSYNC)) { 2380 vp->phase = 0; 2381 vp->digit = mldigit; 2382 } 2383 } else if (vp->count < BCMP) { 2384 vp->count++; 2385 up->status |= DSYNC; 2386 if (!(up->status & INSYNC)) { 2387 vp->phase = 0; 2388 vp->digit = mldigit; 2389 } 2390 } else { 2391 vp->phase = 0; 2392 vp->digit = mldigit; 2393 up->digcnt++; 2394 } 2395 if (vp->digit != mldigit) 2396 up->alarm |= DECERR; 2397 if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status & 2398 INSYNC)) { 2399 sprintf(tbuf, 2400 "wwv4 %2d %04x %5.0f %2d %d %d %d %d %5.0f %5.1f", 2401 up->rsec, up->status, up->epomax, vp->radix, 2402 vp->digit, vp->mldigit, vp->phase, vp->count, 2403 vp->digprb, vp->digsnr); 2404 record_clock_stats(&peer->srcadr, tbuf); 2405 #ifdef DEBUG 2406 if (debug) 2407 printf("%s\n", tbuf); 2408 #endif /* DEBUG */ 2409 } 2410 up->status &= ~BGATE; 2411 } 2412 2413 2414 /* 2415 * wwv_tsec - transmitter minute processing 2416 * 2417 * This routine is called at the end of the transmitter minute. It 2418 * implements a state machine that advances the logical clock subject to 2419 * the funny rules that govern the conventional clock and calendar. 2420 */ 2421 static void 2422 wwv_tsec( 2423 struct wwvunit *up /* driver structure pointer */ 2424 ) 2425 { 2426 int minute, day, isleap; 2427 int temp; 2428 2429 /* 2430 * Advance minute unit of the day. 2431 */ 2432 temp = carry(&up->decvec[MN]); /* minute units */ 2433 2434 /* 2435 * Propagate carries through the day. 2436 */ 2437 if (temp == 0) /* carry minutes */ 2438 temp = carry(&up->decvec[MN + 1]); 2439 if (temp == 0) /* carry hours */ 2440 temp = carry(&up->decvec[HR]); 2441 if (temp == 0) 2442 temp = carry(&up->decvec[HR + 1]); 2443 2444 /* 2445 * Decode the current minute and day. Set leap day if the 2446 * timecode leap bit is set on 30 June or 31 December. Set leap 2447 * minute if the last minute on leap day. This code fails in 2448 * 2400 AD. 2449 */ 2450 minute = up->decvec[MN].digit + up->decvec[MN + 1].digit * 2451 10 + up->decvec[HR].digit * 60 + up->decvec[HR + 2452 1].digit * 600; 2453 day = up->decvec[DA].digit + up->decvec[DA + 1].digit * 10 + 2454 up->decvec[DA + 2].digit * 100; 2455 isleap = (up->decvec[YR].digit & 0x3) == 0; 2456 if (up->misc & SECWAR && (day == (isleap ? 182 : 183) || day == 2457 (isleap ? 365 : 366)) && up->status & INSYNC && up->status & 2458 SSYNC) 2459 up->status |= LEPDAY; 2460 else 2461 up->status &= ~LEPDAY; 2462 if (up->status & LEPDAY && minute == 1439) 2463 up->status |= LEPSEC; 2464 else 2465 up->status &= ~LEPSEC; 2466 2467 /* 2468 * Roll the day if this the first minute and propagate carries 2469 * through the year. 2470 */ 2471 if (minute != 1440) 2472 return; 2473 minute = 0; 2474 while (carry(&up->decvec[HR]) != 0); /* advance to minute 0 */ 2475 while (carry(&up->decvec[HR + 1]) != 0); 2476 day++; 2477 temp = carry(&up->decvec[DA]); /* carry days */ 2478 if (temp == 0) 2479 temp = carry(&up->decvec[DA + 1]); 2480 if (temp == 0) 2481 temp = carry(&up->decvec[DA + 2]); 2482 2483 /* 2484 * Roll the year if this the first day and propagate carries 2485 * through the century. 2486 */ 2487 if (day != (isleap ? 365 : 366)) 2488 return; 2489 day = 1; 2490 while (carry(&up->decvec[DA]) != 1); /* advance to day 1 */ 2491 while (carry(&up->decvec[DA + 1]) != 0); 2492 while (carry(&up->decvec[DA + 2]) != 0); 2493 temp = carry(&up->decvec[YR]); /* carry years */ 2494 if (temp) 2495 carry(&up->decvec[YR + 1]); 2496 } 2497 2498 2499 /* 2500 * carry - process digit 2501 * 2502 * This routine rotates a likelihood vector one position and increments 2503 * the clock digit modulo the radix. It returns the new clock digit or 2504 * zero if a carry occurred. Once synchronized, the clock digit will 2505 * match the maximum likelihood digit corresponding to that position. 2506 */ 2507 static int 2508 carry( 2509 struct decvec *dp /* decoding table pointer */ 2510 ) 2511 { 2512 int temp; 2513 int j; 2514 2515 dp->digit++; /* advance clock digit */ 2516 if (dp->digit == dp->radix) { /* modulo radix */ 2517 dp->digit = 0; 2518 } 2519 temp = dp->like[dp->radix - 1]; /* rotate likelihood vector */ 2520 for (j = dp->radix - 1; j > 0; j--) 2521 dp->like[j] = dp->like[j - 1]; 2522 dp->like[0] = temp; 2523 return (dp->digit); 2524 } 2525 2526 2527 /* 2528 * wwv_snr - compute SNR or likelihood function 2529 */ 2530 static double 2531 wwv_snr( 2532 double signal, /* signal */ 2533 double noise /* noise */ 2534 ) 2535 { 2536 double rval; 2537 2538 /* 2539 * This is a little tricky. Due to the way things are measured, 2540 * either or both the signal or noise amplitude can be negative 2541 * or zero. The intent is that, if the signal is negative or 2542 * zero, the SNR must always be zero. This can happen with the 2543 * subcarrier SNR before the phase has been aligned. On the 2544 * other hand, in the likelihood function the "noise" is the 2545 * next maximum down from the peak and this could be negative. 2546 * However, in this case the SNR is truly stupendous, so we 2547 * simply cap at MAXSNR dB. 2548 */ 2549 if (signal <= 0) { 2550 rval = 0; 2551 } else if (noise <= 0) { 2552 rval = MAXSNR; 2553 } else { 2554 rval = 20 * log10(signal / noise); 2555 if (rval > MAXSNR) 2556 rval = MAXSNR; 2557 } 2558 return (rval); 2559 } 2560 2561 2562 /* 2563 * wwv_newchan - change to new data channel 2564 * 2565 * The radio actually appears to have ten channels, one channel for each 2566 * of five frequencies and each of two stations (WWV and WWVH), although 2567 * if not tunable only the 15 MHz channels appear live. While the radio 2568 * is tuned to the working data channel frequency and station for most 2569 * of the minute, during seconds 59, 0 and 1 the radio is tuned to a 2570 * probe frequency in order to search for minute sync pulse and data 2571 * subcarrier from other transmitters. 2572 * 2573 * The search for WWV and WWVH operates simultaneously, with WWV minute 2574 * sync pulse at 1000 Hz and WWVH at 1200 Hz. The probe frequency 2575 * rotates each minute over 2.5, 5, 10, 15 and 20 MHz in order and yes, 2576 * we all know WWVH is dark on 20 MHz, but few remember when WWV was lit 2577 * on 25 MHz. 2578 * 2579 * This routine selects the best channel using a metric computed from 2580 * the reachability register and minute pulse amplitude. Normally, the 2581 * award goes to the the channel with the highest metric; but, in case 2582 * of ties, the award goes to the channel with the highest minute sync 2583 * pulse amplitude and then to the highest frequency. 2584 * 2585 * The routine performs an important squelch function to keep dirty data 2586 * from polluting the integrators. During acquisition and until the 2587 * clock is synchronized, the signal metric must be at least MTR (13); 2588 * after that the metrict must be at least TTHR (50). If either of these 2589 * is not true, the station select bits are cleared so the second sync 2590 * is disabled and the data bit integrators averaged to a miss. 2591 */ 2592 static void 2593 wwv_newchan( 2594 struct peer *peer /* peer structure pointer */ 2595 ) 2596 { 2597 struct refclockproc *pp; 2598 struct wwvunit *up; 2599 struct sync *sp, *rp; 2600 double rank, dtemp; 2601 int i, j; 2602 2603 pp = peer->procptr; 2604 up = (struct wwvunit *)pp->unitptr; 2605 2606 /* 2607 * Search all five station pairs looking for the channel with 2608 * maximum metric. If no station is found above thresholds, the 2609 * reference ID is set to NONE and we wait for hotter ions. 2610 */ 2611 j = 0; 2612 sp = NULL; 2613 rank = 0; 2614 for (i = 0; i < NCHAN; i++) { 2615 rp = &up->mitig[i].wwvh; 2616 dtemp = wwv_metric(rp); 2617 if (dtemp >= rank) { 2618 rank = dtemp; 2619 sp = rp; 2620 j = i; 2621 } 2622 rp = &up->mitig[i].wwv; 2623 dtemp = wwv_metric(rp); 2624 if (dtemp >= rank) { 2625 rank = dtemp; 2626 sp = rp; 2627 j = i; 2628 } 2629 } 2630 up->dchan = j; 2631 up->sptr = sp; 2632 up->status &= ~(SELV | SELH); 2633 memcpy(&pp->refid, "NONE", 4); 2634 if ((!(up->status & INSYNC) && rank >= MTHR) || ((up->status & 2635 INSYNC) && rank >= TTHR)) { 2636 up->status |= sp->select & (SELV | SELH); 2637 memcpy(&pp->refid, sp->refid, 4); 2638 } 2639 if (peer->stratum <= 1) 2640 memcpy(&peer->refid, &pp->refid, 4); 2641 } 2642 2643 2644 /* 2645 * www_newgame - reset and start over 2646 */ 2647 static void 2648 wwv_newgame( 2649 struct peer *peer /* peer structure pointer */ 2650 ) 2651 { 2652 struct refclockproc *pp; 2653 struct wwvunit *up; 2654 struct chan *cp; 2655 int i; 2656 2657 pp = peer->procptr; 2658 up = (struct wwvunit *)pp->unitptr; 2659 2660 /* 2661 * Initialize strategic values. Note we set the leap bits 2662 * NOTINSYNC and the refid "NONE". 2663 */ 2664 peer->leap = LEAP_NOTINSYNC; 2665 up->watch = up->status = up->alarm = 0; 2666 up->avgint = MINAVG; 2667 up->freq = 0; 2668 up->sptr = NULL; 2669 up->gain = MAXGAIN / 2; 2670 2671 /* 2672 * Initialize the station processes for audio gain, select bit, 2673 * station/frequency identifier and reference identifier. 2674 */ 2675 memset(up->mitig, 0, sizeof(up->mitig)); 2676 for (i = 0; i < NCHAN; i++) { 2677 cp = &up->mitig[i]; 2678 cp->gain = up->gain; 2679 cp->wwv.select = SELV; 2680 sprintf(cp->wwv.refid, "WV%.0f", floor(qsy[i])); 2681 cp->wwvh.select = SELH; 2682 sprintf(cp->wwvh.refid, "WH%.0f", floor(qsy[i])); 2683 } 2684 wwv_newchan(peer); 2685 } 2686 2687 /* 2688 * wwv_metric - compute station metric 2689 * 2690 * The most significant bits represent the number of ones in the 2691 * reachability register. The least significant bits represent the 2692 * minute sync pulse amplitude. The combined value is scaled 0-100. 2693 */ 2694 double 2695 wwv_metric( 2696 struct sync *sp /* station pointer */ 2697 ) 2698 { 2699 double dtemp; 2700 2701 dtemp = sp->count * MAXSIG; 2702 if (sp->synmax < MAXSIG) 2703 dtemp += sp->synmax; 2704 else 2705 dtemp += MAXSIG - 1; 2706 dtemp /= (AMAX + 1) * MAXSIG; 2707 return (dtemp * 100.); 2708 } 2709 2710 2711 #ifdef ICOM 2712 /* 2713 * wwv_qsy - Tune ICOM receiver 2714 * 2715 * This routine saves the AGC for the current channel, switches to a new 2716 * channel and restores the AGC for that channel. If a tunable receiver 2717 * is not available, just fake it. 2718 */ 2719 static int 2720 wwv_qsy( 2721 struct peer *peer, /* peer structure pointer */ 2722 int chan /* channel */ 2723 ) 2724 { 2725 int rval = 0; 2726 struct refclockproc *pp; 2727 struct wwvunit *up; 2728 2729 pp = peer->procptr; 2730 up = (struct wwvunit *)pp->unitptr; 2731 if (up->fd_icom > 0) { 2732 up->mitig[up->achan].gain = up->gain; 2733 rval = icom_freq(up->fd_icom, peer->ttl & 0x7f, 2734 qsy[chan]); 2735 up->achan = chan; 2736 up->gain = up->mitig[up->achan].gain; 2737 } 2738 return (rval); 2739 } 2740 #endif /* ICOM */ 2741 2742 2743 /* 2744 * timecode - assemble timecode string and length 2745 * 2746 * Prettytime format - similar to Spectracom 2747 * 2748 * sq yy ddd hh:mm:ss ld dut lset agc iden sig errs freq avgt 2749 * 2750 * s sync indicator ('?' or ' ') 2751 * q error bits (hex 0-F) 2752 * yyyy year of century 2753 * ddd day of year 2754 * hh hour of day 2755 * mm minute of hour 2756 * ss second of minute) 2757 * l leap second warning (' ' or 'L') 2758 * d DST state ('S', 'D', 'I', or 'O') 2759 * dut DUT sign and magnitude (0.1 s) 2760 * lset minutes since last clock update 2761 * agc audio gain (0-255) 2762 * iden reference identifier (station and frequency) 2763 * sig signal quality (0-100) 2764 * errs bit errors in last minute 2765 * freq frequency offset (PPM) 2766 * avgt averaging time (s) 2767 */ 2768 static int 2769 timecode( 2770 struct wwvunit *up, /* driver structure pointer */ 2771 char *ptr /* target string */ 2772 ) 2773 { 2774 struct sync *sp; 2775 int year, day, hour, minute, second, dut; 2776 char synchar, leapchar, dst; 2777 char cptr[50]; 2778 2779 2780 /* 2781 * Common fixed-format fields 2782 */ 2783 synchar = (up->status & INSYNC) ? ' ' : '?'; 2784 year = up->decvec[YR].digit + up->decvec[YR + 1].digit * 10 + 2785 2000; 2786 day = up->decvec[DA].digit + up->decvec[DA + 1].digit * 10 + 2787 up->decvec[DA + 2].digit * 100; 2788 hour = up->decvec[HR].digit + up->decvec[HR + 1].digit * 10; 2789 minute = up->decvec[MN].digit + up->decvec[MN + 1].digit * 10; 2790 second = 0; 2791 leapchar = (up->misc & SECWAR) ? 'L' : ' '; 2792 dst = dstcod[(up->misc >> 4) & 0x3]; 2793 dut = up->misc & 0x7; 2794 if (!(up->misc & DUTS)) 2795 dut = -dut; 2796 sprintf(ptr, "%c%1X", synchar, up->alarm); 2797 sprintf(cptr, " %4d %03d %02d:%02d:%02d %c%c %+d", 2798 year, day, hour, minute, second, leapchar, dst, dut); 2799 strcat(ptr, cptr); 2800 2801 /* 2802 * Specific variable-format fields 2803 */ 2804 sp = up->sptr; 2805 sprintf(cptr, " %d %d %s %.0f %d %.1f %d", up->watch, 2806 up->mitig[up->dchan].gain, sp->refid, wwv_metric(sp), 2807 up->errbit, up->freq / SECOND * 1e6, up->avgint); 2808 strcat(ptr, cptr); 2809 return (strlen(ptr)); 2810 } 2811 2812 2813 /* 2814 * wwv_gain - adjust codec gain 2815 * 2816 * This routine is called at the end of each second. It counts the 2817 * number of signal clips above the MAXSIG threshold during the previous 2818 * second. If there are no clips, the gain is bumped up; if too many 2819 * clips, it is bumped down. The decoder is relatively insensitive to 2820 * amplitude, so this crudity works just fine. The input port is set and 2821 * the error flag is cleared, mostly to be ornery. 2822 */ 2823 static void 2824 wwv_gain( 2825 struct peer *peer /* peer structure pointer */ 2826 ) 2827 { 2828 struct refclockproc *pp; 2829 struct wwvunit *up; 2830 2831 pp = peer->procptr; 2832 up = (struct wwvunit *)pp->unitptr; 2833 2834 /* 2835 * Apparently, the codec uses only the high order bits of the 2836 * gain control field. Thus, it may take awhile for changes to 2837 * wiggle the hardware bits. 2838 */ 2839 if (up->clipcnt == 0) { 2840 up->gain += 4; 2841 if (up->gain > MAXGAIN) 2842 up->gain = MAXGAIN; 2843 } else if (up->clipcnt > MAXCLP) { 2844 up->gain -= 4; 2845 if (up->gain < 0) 2846 up->gain = 0; 2847 } 2848 audio_gain(up->gain, up->mongain, up->port); 2849 up->clipcnt = 0; 2850 #if DEBUG 2851 if (debug > 1) 2852 audio_show(); 2853 #endif 2854 } 2855 2856 2857 #else 2858 int refclock_wwv_bs; 2859 #endif /* REFCLOCK */ 2860