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