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