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