Annotation of embedaddon/ntp/ntpd/refclock_wwv.c, revision 1.1
1.1 ! misho 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 dubugging 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 SECOND 8000 /* second epoch (sample rate) (Hz) */
! 85: #define MINUTE (SECOND * 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 (SECOND / 1000) /* samples per millisecond */
! 205: #define IN100 ((100 * 80) / SECOND) /* 100 Hz increment */
! 206: #define IN1000 ((1000 * 80) / SECOND) /* 1000 Hz increment */
! 207: #define IN1200 ((1200 * 80) / SECOND) /* 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 *);
! 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: if (!(up = (struct wwvunit *)emalloc(sizeof(struct wwvunit)))) {
! 663: close(fd);
! 664: return (0);
! 665: }
! 666: memset(up, 0, sizeof(struct wwvunit));
! 667: pp = peer->procptr;
! 668: pp->unitptr = (caddr_t)up;
! 669: pp->io.clock_recv = wwv_receive;
! 670: pp->io.srcclock = (caddr_t)peer;
! 671: pp->io.datalen = 0;
! 672: pp->io.fd = fd;
! 673: if (!io_addclock(&pp->io)) {
! 674: close(fd);
! 675: free(up);
! 676: return (0);
! 677: }
! 678:
! 679: /*
! 680: * Initialize miscellaneous variables
! 681: */
! 682: peer->precision = PRECISION;
! 683: pp->clockdesc = DESCRIPTION;
! 684:
! 685: /*
! 686: * The companded samples are encoded sign-magnitude. The table
! 687: * contains all the 256 values in the interest of speed.
! 688: */
! 689: up->comp[0] = up->comp[OFFSET] = 0.;
! 690: up->comp[1] = 1.; up->comp[OFFSET + 1] = -1.;
! 691: up->comp[2] = 3.; up->comp[OFFSET + 2] = -3.;
! 692: step = 2.;
! 693: for (i = 3; i < OFFSET; i++) {
! 694: up->comp[i] = up->comp[i - 1] + step;
! 695: up->comp[OFFSET + i] = -up->comp[i];
! 696: if (i % 16 == 0)
! 697: step *= 2.;
! 698: }
! 699: DTOLFP(1. / SECOND, &up->tick);
! 700:
! 701: /*
! 702: * Initialize the decoding matrix with the radix for each digit
! 703: * position.
! 704: */
! 705: up->decvec[MN].radix = 10; /* minutes */
! 706: up->decvec[MN + 1].radix = 6;
! 707: up->decvec[HR].radix = 10; /* hours */
! 708: up->decvec[HR + 1].radix = 3;
! 709: up->decvec[DA].radix = 10; /* days */
! 710: up->decvec[DA + 1].radix = 10;
! 711: up->decvec[DA + 2].radix = 4;
! 712: up->decvec[YR].radix = 10; /* years */
! 713: up->decvec[YR + 1].radix = 10;
! 714:
! 715: #ifdef ICOM
! 716: /*
! 717: * Initialize autotune if available. Note that the ICOM select
! 718: * code must be less than 128, so the high order bit can be used
! 719: * to select the line speed 0 (9600 bps) or 1 (1200 bps). Note
! 720: * we don't complain if the ICOM device is not there; but, if it
! 721: * is, the radio better be working.
! 722: */
! 723: temp = 0;
! 724: #ifdef DEBUG
! 725: if (debug > 1)
! 726: temp = P_TRACE;
! 727: #endif /* DEBUG */
! 728: if (peer->ttl != 0) {
! 729: if (peer->ttl & 0x80)
! 730: up->fd_icom = icom_init("/dev/icom", B1200,
! 731: temp);
! 732: else
! 733: up->fd_icom = icom_init("/dev/icom", B9600,
! 734: temp);
! 735: }
! 736: if (up->fd_icom > 0) {
! 737: if (wwv_qsy(peer, DCHAN) != 0) {
! 738: msyslog(LOG_NOTICE, "icom: radio not found");
! 739: close(up->fd_icom);
! 740: up->fd_icom = 0;
! 741: } else {
! 742: msyslog(LOG_NOTICE, "icom: autotune enabled");
! 743: }
! 744: }
! 745: #endif /* ICOM */
! 746:
! 747: /*
! 748: * Let the games begin.
! 749: */
! 750: wwv_newgame(peer);
! 751: return (1);
! 752: }
! 753:
! 754:
! 755: /*
! 756: * wwv_shutdown - shut down the clock
! 757: */
! 758: static void
! 759: wwv_shutdown(
! 760: int unit, /* instance number (not used) */
! 761: struct peer *peer /* peer structure pointer */
! 762: )
! 763: {
! 764: struct refclockproc *pp;
! 765: struct wwvunit *up;
! 766:
! 767: pp = peer->procptr;
! 768: up = (struct wwvunit *)pp->unitptr;
! 769: if (up == NULL)
! 770: return;
! 771:
! 772: io_closeclock(&pp->io);
! 773: #ifdef ICOM
! 774: if (up->fd_icom > 0)
! 775: close(up->fd_icom);
! 776: #endif /* ICOM */
! 777: free(up);
! 778: }
! 779:
! 780:
! 781: /*
! 782: * wwv_receive - receive data from the audio device
! 783: *
! 784: * This routine reads input samples and adjusts the logical clock to
! 785: * track the A/D sample clock by dropping or duplicating codec samples.
! 786: * It also controls the A/D signal level with an AGC loop to mimimize
! 787: * quantization noise and avoid overload.
! 788: */
! 789: static void
! 790: wwv_receive(
! 791: struct recvbuf *rbufp /* receive buffer structure pointer */
! 792: )
! 793: {
! 794: struct peer *peer;
! 795: struct refclockproc *pp;
! 796: struct wwvunit *up;
! 797:
! 798: /*
! 799: * Local variables
! 800: */
! 801: double sample; /* codec sample */
! 802: u_char *dpt; /* buffer pointer */
! 803: int bufcnt; /* buffer counter */
! 804: l_fp ltemp;
! 805:
! 806: peer = (struct peer *)rbufp->recv_srcclock;
! 807: pp = peer->procptr;
! 808: up = (struct wwvunit *)pp->unitptr;
! 809:
! 810: /*
! 811: * Main loop - read until there ain't no more. Note codec
! 812: * samples are bit-inverted.
! 813: */
! 814: DTOLFP((double)rbufp->recv_length / SECOND, <emp);
! 815: L_SUB(&rbufp->recv_time, <emp);
! 816: up->timestamp = rbufp->recv_time;
! 817: dpt = rbufp->recv_buffer;
! 818: for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
! 819: sample = up->comp[~*dpt++ & 0xff];
! 820:
! 821: /*
! 822: * Clip noise spikes greater than MAXAMP (6000) and
! 823: * record the number of clips to be used later by the
! 824: * AGC.
! 825: */
! 826: if (sample > MAXAMP) {
! 827: sample = MAXAMP;
! 828: up->clipcnt++;
! 829: } else if (sample < -MAXAMP) {
! 830: sample = -MAXAMP;
! 831: up->clipcnt++;
! 832: }
! 833:
! 834: /*
! 835: * Variable frequency oscillator. The codec oscillator
! 836: * runs at the nominal rate of 8000 samples per second,
! 837: * or 125 us per sample. A frequency change of one unit
! 838: * results in either duplicating or deleting one sample
! 839: * per second, which results in a frequency change of
! 840: * 125 PPM.
! 841: */
! 842: up->phase += (up->freq + clock_codec) / SECOND;
! 843: if (up->phase >= .5) {
! 844: up->phase -= 1.;
! 845: } else if (up->phase < -.5) {
! 846: up->phase += 1.;
! 847: wwv_rf(peer, sample);
! 848: wwv_rf(peer, sample);
! 849: } else {
! 850: wwv_rf(peer, sample);
! 851: }
! 852: L_ADD(&up->timestamp, &up->tick);
! 853: }
! 854:
! 855: /*
! 856: * Set the input port and monitor gain for the next buffer.
! 857: */
! 858: if (pp->sloppyclockflag & CLK_FLAG2)
! 859: up->port = 2;
! 860: else
! 861: up->port = 1;
! 862: if (pp->sloppyclockflag & CLK_FLAG3)
! 863: up->mongain = MONGAIN;
! 864: else
! 865: up->mongain = 0;
! 866: }
! 867:
! 868:
! 869: /*
! 870: * wwv_poll - called by the transmit procedure
! 871: *
! 872: * This routine keeps track of status. If no offset samples have been
! 873: * processed during a poll interval, a timeout event is declared. If
! 874: * errors have have occurred during the interval, they are reported as
! 875: * well.
! 876: */
! 877: static void
! 878: wwv_poll(
! 879: int unit, /* instance number (not used) */
! 880: struct peer *peer /* peer structure pointer */
! 881: )
! 882: {
! 883: struct refclockproc *pp;
! 884: struct wwvunit *up;
! 885:
! 886: pp = peer->procptr;
! 887: up = (struct wwvunit *)pp->unitptr;
! 888: if (up->errflg)
! 889: refclock_report(peer, up->errflg);
! 890: up->errflg = 0;
! 891: pp->polls++;
! 892: }
! 893:
! 894:
! 895: /*
! 896: * wwv_rf - process signals and demodulate to baseband
! 897: *
! 898: * This routine grooms and filters decompanded raw audio samples. The
! 899: * output signal is the 100-Hz filtered baseband data signal in
! 900: * quadrature phase. The routine also determines the minute synch epoch,
! 901: * as well as certain signal maxima, minima and related values.
! 902: *
! 903: * There are two 1-s ramps used by this program. Both count the 8000
! 904: * logical clock samples spanning exactly one second. The epoch ramp
! 905: * counts the samples starting at an arbitrary time. The rphase ramp
! 906: * counts the samples starting at the 5-ms second sync pulse found
! 907: * during the epoch ramp.
! 908: *
! 909: * There are two 1-m ramps used by this program. The mphase ramp counts
! 910: * the 480,000 logical clock samples spanning exactly one minute and
! 911: * starting at an arbitrary time. The rsec ramp counts the 60 seconds of
! 912: * the minute starting at the 800-ms minute sync pulse found during the
! 913: * mphase ramp. The rsec ramp drives the seconds state machine to
! 914: * determine the bits and digits of the timecode.
! 915: *
! 916: * Demodulation operations are based on three synthesized quadrature
! 917: * sinusoids: 100 Hz for the data signal, 1000 Hz for the WWV sync
! 918: * signal and 1200 Hz for the WWVH sync signal. These drive synchronous
! 919: * matched filters for the data signal (170 ms at 100 Hz), WWV minute
! 920: * sync signal (800 ms at 1000 Hz) and WWVH minute sync signal (800 ms
! 921: * at 1200 Hz). Two additional matched filters are switched in
! 922: * as required for the WWV second sync signal (5 cycles at 1000 Hz) and
! 923: * WWVH second sync signal (6 cycles at 1200 Hz).
! 924: */
! 925: static void
! 926: wwv_rf(
! 927: struct peer *peer, /* peerstructure pointer */
! 928: double isig /* input signal */
! 929: )
! 930: {
! 931: struct refclockproc *pp;
! 932: struct wwvunit *up;
! 933: struct sync *sp, *rp;
! 934:
! 935: static double lpf[5]; /* 150-Hz lpf delay line */
! 936: double data; /* lpf output */
! 937: static double bpf[9]; /* 1000/1200-Hz bpf delay line */
! 938: double syncx; /* bpf output */
! 939: static double mf[41]; /* 1000/1200-Hz mf delay line */
! 940: double mfsync; /* mf output */
! 941:
! 942: static int iptr; /* data channel pointer */
! 943: static double ibuf[DATSIZ]; /* data I channel delay line */
! 944: static double qbuf[DATSIZ]; /* data Q channel delay line */
! 945:
! 946: static int jptr; /* sync channel pointer */
! 947: static int kptr; /* tick channel pointer */
! 948:
! 949: static int csinptr; /* wwv channel phase */
! 950: static double cibuf[SYNSIZ]; /* wwv I channel delay line */
! 951: static double cqbuf[SYNSIZ]; /* wwv Q channel delay line */
! 952: static double ciamp; /* wwv I channel amplitude */
! 953: static double cqamp; /* wwv Q channel amplitude */
! 954:
! 955: static double csibuf[TCKSIZ]; /* wwv I tick delay line */
! 956: static double csqbuf[TCKSIZ]; /* wwv Q tick delay line */
! 957: static double csiamp; /* wwv I tick amplitude */
! 958: static double csqamp; /* wwv Q tick amplitude */
! 959:
! 960: static int hsinptr; /* wwvh channel phase */
! 961: static double hibuf[SYNSIZ]; /* wwvh I channel delay line */
! 962: static double hqbuf[SYNSIZ]; /* wwvh Q channel delay line */
! 963: static double hiamp; /* wwvh I channel amplitude */
! 964: static double hqamp; /* wwvh Q channel amplitude */
! 965:
! 966: static double hsibuf[TCKSIZ]; /* wwvh I tick delay line */
! 967: static double hsqbuf[TCKSIZ]; /* wwvh Q tick delay line */
! 968: static double hsiamp; /* wwvh I tick amplitude */
! 969: static double hsqamp; /* wwvh Q tick amplitude */
! 970:
! 971: static double epobuf[SECOND]; /* second sync comb filter */
! 972: static double epomax, nxtmax; /* second sync amplitude buffer */
! 973: static int epopos; /* epoch second sync position buffer */
! 974:
! 975: static int iniflg; /* initialization flag */
! 976: int epoch; /* comb filter index */
! 977: double dtemp;
! 978: int i;
! 979:
! 980: pp = peer->procptr;
! 981: up = (struct wwvunit *)pp->unitptr;
! 982:
! 983: if (!iniflg) {
! 984: iniflg = 1;
! 985: memset((char *)lpf, 0, sizeof(lpf));
! 986: memset((char *)bpf, 0, sizeof(bpf));
! 987: memset((char *)mf, 0, sizeof(mf));
! 988: memset((char *)ibuf, 0, sizeof(ibuf));
! 989: memset((char *)qbuf, 0, sizeof(qbuf));
! 990: memset((char *)cibuf, 0, sizeof(cibuf));
! 991: memset((char *)cqbuf, 0, sizeof(cqbuf));
! 992: memset((char *)csibuf, 0, sizeof(csibuf));
! 993: memset((char *)csqbuf, 0, sizeof(csqbuf));
! 994: memset((char *)hibuf, 0, sizeof(hibuf));
! 995: memset((char *)hqbuf, 0, sizeof(hqbuf));
! 996: memset((char *)hsibuf, 0, sizeof(hsibuf));
! 997: memset((char *)hsqbuf, 0, sizeof(hsqbuf));
! 998: memset((char *)epobuf, 0, sizeof(epobuf));
! 999: }
! 1000:
! 1001: /*
! 1002: * Baseband data demodulation. The 100-Hz subcarrier is
! 1003: * extracted using a 150-Hz IIR lowpass filter. This attenuates
! 1004: * the 1000/1200-Hz sync signals, as well as the 440-Hz and
! 1005: * 600-Hz tones and most of the noise and voice modulation
! 1006: * components.
! 1007: *
! 1008: * The subcarrier is transmitted 10 dB down from the carrier.
! 1009: * The DGAIN parameter can be adjusted for this and to
! 1010: * compensate for the radio audio response at 100 Hz.
! 1011: *
! 1012: * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB
! 1013: * passband ripple, -50 dB stopband ripple, phase delay 0.97 ms.
! 1014: */
! 1015: data = (lpf[4] = lpf[3]) * 8.360961e-01;
! 1016: data += (lpf[3] = lpf[2]) * -3.481740e+00;
! 1017: data += (lpf[2] = lpf[1]) * 5.452988e+00;
! 1018: data += (lpf[1] = lpf[0]) * -3.807229e+00;
! 1019: lpf[0] = isig * DGAIN - data;
! 1020: data = lpf[0] * 3.281435e-03
! 1021: + lpf[1] * -1.149947e-02
! 1022: + lpf[2] * 1.654858e-02
! 1023: + lpf[3] * -1.149947e-02
! 1024: + lpf[4] * 3.281435e-03;
! 1025:
! 1026: /*
! 1027: * The 100-Hz data signal is demodulated using a pair of
! 1028: * quadrature multipliers, matched filters and a phase lock
! 1029: * loop. The I and Q quadrature data signals are produced by
! 1030: * multiplying the filtered signal by 100-Hz sine and cosine
! 1031: * signals, respectively. The signals are processed by 170-ms
! 1032: * synchronous matched filters to produce the amplitude and
! 1033: * phase signals used by the demodulator. The signals are scaled
! 1034: * to produce unit energy at the maximum value.
! 1035: */
! 1036: i = up->datapt;
! 1037: up->datapt = (up->datapt + IN100) % 80;
! 1038: dtemp = sintab[i] * data / (MS / 2. * DATCYC);
! 1039: up->irig -= ibuf[iptr];
! 1040: ibuf[iptr] = dtemp;
! 1041: up->irig += dtemp;
! 1042:
! 1043: i = (i + 20) % 80;
! 1044: dtemp = sintab[i] * data / (MS / 2. * DATCYC);
! 1045: up->qrig -= qbuf[iptr];
! 1046: qbuf[iptr] = dtemp;
! 1047: up->qrig += dtemp;
! 1048: iptr = (iptr + 1) % DATSIZ;
! 1049:
! 1050: /*
! 1051: * Baseband sync demodulation. The 1000/1200 sync signals are
! 1052: * extracted using a 600-Hz IIR bandpass filter. This removes
! 1053: * the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz
! 1054: * tones and most of the noise and voice modulation components.
! 1055: *
! 1056: * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB
! 1057: * passband ripple, -50 dB stopband ripple, phase delay 0.91 ms.
! 1058: */
! 1059: syncx = (bpf[8] = bpf[7]) * 4.897278e-01;
! 1060: syncx += (bpf[7] = bpf[6]) * -2.765914e+00;
! 1061: syncx += (bpf[6] = bpf[5]) * 8.110921e+00;
! 1062: syncx += (bpf[5] = bpf[4]) * -1.517732e+01;
! 1063: syncx += (bpf[4] = bpf[3]) * 1.975197e+01;
! 1064: syncx += (bpf[3] = bpf[2]) * -1.814365e+01;
! 1065: syncx += (bpf[2] = bpf[1]) * 1.159783e+01;
! 1066: syncx += (bpf[1] = bpf[0]) * -4.735040e+00;
! 1067: bpf[0] = isig - syncx;
! 1068: syncx = bpf[0] * 8.203628e-03
! 1069: + bpf[1] * -2.375732e-02
! 1070: + bpf[2] * 3.353214e-02
! 1071: + bpf[3] * -4.080258e-02
! 1072: + bpf[4] * 4.605479e-02
! 1073: + bpf[5] * -4.080258e-02
! 1074: + bpf[6] * 3.353214e-02
! 1075: + bpf[7] * -2.375732e-02
! 1076: + bpf[8] * 8.203628e-03;
! 1077:
! 1078: /*
! 1079: * The 1000/1200 sync signals are demodulated using a pair of
! 1080: * quadrature multipliers and matched filters. However,
! 1081: * synchronous demodulation at these frequencies is impractical,
! 1082: * so only the signal amplitude is used. The I and Q quadrature
! 1083: * sync signals are produced by multiplying the filtered signal
! 1084: * by 1000-Hz (WWV) and 1200-Hz (WWVH) sine and cosine signals,
! 1085: * respectively. The WWV and WWVH signals are processed by 800-
! 1086: * ms synchronous matched filters and combined to produce the
! 1087: * minute sync signal and detect which one (or both) the WWV or
! 1088: * WWVH signal is present. The WWV and WWVH signals are also
! 1089: * processed by 5-ms synchronous matched filters and combined to
! 1090: * produce the second sync signal. The signals are scaled to
! 1091: * produce unit energy at the maximum value.
! 1092: *
! 1093: * Note the master timing ramps, which run continuously. The
! 1094: * minute counter (mphase) counts the samples in the minute,
! 1095: * while the second counter (epoch) counts the samples in the
! 1096: * second.
! 1097: */
! 1098: up->mphase = (up->mphase + 1) % MINUTE;
! 1099: epoch = up->mphase % SECOND;
! 1100:
! 1101: /*
! 1102: * WWV
! 1103: */
! 1104: i = csinptr;
! 1105: csinptr = (csinptr + IN1000) % 80;
! 1106:
! 1107: dtemp = sintab[i] * syncx / (MS / 2.);
! 1108: ciamp -= cibuf[jptr];
! 1109: cibuf[jptr] = dtemp;
! 1110: ciamp += dtemp;
! 1111: csiamp -= csibuf[kptr];
! 1112: csibuf[kptr] = dtemp;
! 1113: csiamp += dtemp;
! 1114:
! 1115: i = (i + 20) % 80;
! 1116: dtemp = sintab[i] * syncx / (MS / 2.);
! 1117: cqamp -= cqbuf[jptr];
! 1118: cqbuf[jptr] = dtemp;
! 1119: cqamp += dtemp;
! 1120: csqamp -= csqbuf[kptr];
! 1121: csqbuf[kptr] = dtemp;
! 1122: csqamp += dtemp;
! 1123:
! 1124: sp = &up->mitig[up->achan].wwv;
! 1125: sp->amp = sqrt(ciamp * ciamp + cqamp * cqamp) / SYNCYC;
! 1126: if (!(up->status & MSYNC))
! 1127: wwv_qrz(peer, sp, (int)(pp->fudgetime1 * SECOND));
! 1128:
! 1129: /*
! 1130: * WWVH
! 1131: */
! 1132: i = hsinptr;
! 1133: hsinptr = (hsinptr + IN1200) % 80;
! 1134:
! 1135: dtemp = sintab[i] * syncx / (MS / 2.);
! 1136: hiamp -= hibuf[jptr];
! 1137: hibuf[jptr] = dtemp;
! 1138: hiamp += dtemp;
! 1139: hsiamp -= hsibuf[kptr];
! 1140: hsibuf[kptr] = dtemp;
! 1141: hsiamp += dtemp;
! 1142:
! 1143: i = (i + 20) % 80;
! 1144: dtemp = sintab[i] * syncx / (MS / 2.);
! 1145: hqamp -= hqbuf[jptr];
! 1146: hqbuf[jptr] = dtemp;
! 1147: hqamp += dtemp;
! 1148: hsqamp -= hsqbuf[kptr];
! 1149: hsqbuf[kptr] = dtemp;
! 1150: hsqamp += dtemp;
! 1151:
! 1152: rp = &up->mitig[up->achan].wwvh;
! 1153: rp->amp = sqrt(hiamp * hiamp + hqamp * hqamp) / SYNCYC;
! 1154: if (!(up->status & MSYNC))
! 1155: wwv_qrz(peer, rp, (int)(pp->fudgetime2 * SECOND));
! 1156: jptr = (jptr + 1) % SYNSIZ;
! 1157: kptr = (kptr + 1) % TCKSIZ;
! 1158:
! 1159: /*
! 1160: * The following section is called once per minute. It does
! 1161: * housekeeping and timeout functions and empties the dustbins.
! 1162: */
! 1163: if (up->mphase == 0) {
! 1164: up->watch++;
! 1165: if (!(up->status & MSYNC)) {
! 1166:
! 1167: /*
! 1168: * If minute sync has not been acquired before
! 1169: * ACQSN timeout (6 min), or if no signal is
! 1170: * heard, the program cycles to the next
! 1171: * frequency and tries again.
! 1172: */
! 1173: if (!wwv_newchan(peer))
! 1174: up->watch = 0;
! 1175: } else {
! 1176:
! 1177: /*
! 1178: * If the leap bit is set, set the minute epoch
! 1179: * back one second so the station processes
! 1180: * don't miss a beat.
! 1181: */
! 1182: if (up->status & LEPSEC) {
! 1183: up->mphase -= SECOND;
! 1184: if (up->mphase < 0)
! 1185: up->mphase += MINUTE;
! 1186: }
! 1187: }
! 1188: }
! 1189:
! 1190: /*
! 1191: * When the channel metric reaches threshold and the second
! 1192: * counter matches the minute epoch within the second, the
! 1193: * driver has synchronized to the station. The second number is
! 1194: * the remaining seconds until the next minute epoch, while the
! 1195: * sync epoch is zero. Watch out for the first second; if
! 1196: * already synchronized to the second, the buffered sync epoch
! 1197: * must be set.
! 1198: *
! 1199: * Note the guard interval is 200 ms; if for some reason the
! 1200: * clock drifts more than that, it might wind up in the wrong
! 1201: * second. If the maximum frequency error is not more than about
! 1202: * 1 PPM, the clock can go as much as two days while still in
! 1203: * the same second.
! 1204: */
! 1205: if (up->status & MSYNC) {
! 1206: wwv_epoch(peer);
! 1207: } else if (up->sptr != NULL) {
! 1208: sp = up->sptr;
! 1209: if (sp->metric >= TTHR && epoch == sp->mepoch % SECOND)
! 1210: {
! 1211: up->rsec = (60 - sp->mepoch / SECOND) % 60;
! 1212: up->rphase = 0;
! 1213: up->status |= MSYNC;
! 1214: up->watch = 0;
! 1215: if (!(up->status & SSYNC))
! 1216: up->repoch = up->yepoch = epoch;
! 1217: else
! 1218: up->repoch = up->yepoch;
! 1219:
! 1220: }
! 1221: }
! 1222:
! 1223: /*
! 1224: * The second sync pulse is extracted using 5-ms (40 sample) FIR
! 1225: * matched filters at 1000 Hz for WWV or 1200 Hz for WWVH. This
! 1226: * pulse is used for the most precise synchronization, since if
! 1227: * provides a resolution of one sample (125 us). The filters run
! 1228: * only if the station has been reliably determined.
! 1229: */
! 1230: if (up->status & SELV)
! 1231: mfsync = sqrt(csiamp * csiamp + csqamp * csqamp) /
! 1232: TCKCYC;
! 1233: else if (up->status & SELH)
! 1234: mfsync = sqrt(hsiamp * hsiamp + hsqamp * hsqamp) /
! 1235: TCKCYC;
! 1236: else
! 1237: mfsync = 0;
! 1238:
! 1239: /*
! 1240: * Enhance the seconds sync pulse using a 1-s (8000-sample) comb
! 1241: * filter. Correct for the FIR matched filter delay, which is 5
! 1242: * ms for both the WWV and WWVH filters, and also for the
! 1243: * propagation delay. Once each second look for second sync. If
! 1244: * not in minute sync, fiddle the codec gain. Note the SNR is
! 1245: * computed from the maximum sample and the envelope of the
! 1246: * sample 6 ms before it, so if we slip more than a cycle the
! 1247: * SNR should plummet. The signal is scaled to produce unit
! 1248: * energy at the maximum value.
! 1249: */
! 1250: dtemp = (epobuf[epoch] += (mfsync - epobuf[epoch]) /
! 1251: up->avgint);
! 1252: if (dtemp > epomax) {
! 1253: int j;
! 1254:
! 1255: epomax = dtemp;
! 1256: epopos = epoch;
! 1257: j = epoch - 6 * MS;
! 1258: if (j < 0)
! 1259: j += SECOND;
! 1260: nxtmax = fabs(epobuf[j]);
! 1261: }
! 1262: if (epoch == 0) {
! 1263: up->epomax = epomax;
! 1264: up->eposnr = wwv_snr(epomax, nxtmax);
! 1265: epopos -= TCKCYC * MS;
! 1266: if (epopos < 0)
! 1267: epopos += SECOND;
! 1268: wwv_endpoc(peer, epopos);
! 1269: if (!(up->status & SSYNC))
! 1270: up->alarm |= SYNERR;
! 1271: epomax = 0;
! 1272: if (!(up->status & MSYNC))
! 1273: wwv_gain(peer);
! 1274: }
! 1275: }
! 1276:
! 1277:
! 1278: /*
! 1279: * wwv_qrz - identify and acquire WWV/WWVH minute sync pulse
! 1280: *
! 1281: * This routine implements a virtual station process used to acquire
! 1282: * minute sync and to mitigate among the ten frequency and station
! 1283: * combinations. During minute sync acquisition the process probes each
! 1284: * frequency and station in turn for the minute pulse, which
! 1285: * involves searching through the entire 480,000-sample minute. The
! 1286: * process finds the maximum signal and RMS noise plus signal. Then, the
! 1287: * actual noise is determined by subtracting the energy of the matched
! 1288: * filter.
! 1289: *
! 1290: * Students of radar receiver technology will discover this algorithm
! 1291: * amounts to a range-gate discriminator. A valid pulse must have peak
! 1292: * amplitude at least QTHR (2500) and SNR at least QSNR (20) dB and the
! 1293: * difference between the current and previous epoch must be less than
! 1294: * AWND (20 ms). Note that the discriminator peak occurs about 800 ms
! 1295: * into the second, so the timing is retarded to the previous second
! 1296: * epoch.
! 1297: */
! 1298: static void
! 1299: wwv_qrz(
! 1300: struct peer *peer, /* peer structure pointer */
! 1301: struct sync *sp, /* sync channel structure */
! 1302: int pdelay /* propagation delay (samples) */
! 1303: )
! 1304: {
! 1305: struct refclockproc *pp;
! 1306: struct wwvunit *up;
! 1307: char tbuf[TBUF]; /* monitor buffer */
! 1308: long epoch;
! 1309:
! 1310: pp = peer->procptr;
! 1311: up = (struct wwvunit *)pp->unitptr;
! 1312:
! 1313: /*
! 1314: * Find the sample with peak amplitude, which defines the minute
! 1315: * epoch. Accumulate all samples to determine the total noise
! 1316: * energy.
! 1317: */
! 1318: epoch = up->mphase - pdelay - SYNSIZ;
! 1319: if (epoch < 0)
! 1320: epoch += MINUTE;
! 1321: if (sp->amp > sp->maxeng) {
! 1322: sp->maxeng = sp->amp;
! 1323: sp->pos = epoch;
! 1324: }
! 1325: sp->noieng += sp->amp;
! 1326:
! 1327: /*
! 1328: * At the end of the minute, determine the epoch of the minute
! 1329: * sync pulse, as well as the difference between the current and
! 1330: * previous epoches due to the intrinsic frequency error plus
! 1331: * jitter. When calculating the SNR, subtract the pulse energy
! 1332: * from the total noise energy and then normalize.
! 1333: */
! 1334: if (up->mphase == 0) {
! 1335: sp->synmax = sp->maxeng;
! 1336: sp->synsnr = wwv_snr(sp->synmax, (sp->noieng -
! 1337: sp->synmax) / MINUTE);
! 1338: if (sp->count == 0)
! 1339: sp->lastpos = sp->pos;
! 1340: epoch = (sp->pos - sp->lastpos) % MINUTE;
! 1341: sp->reach <<= 1;
! 1342: if (sp->reach & (1 << AMAX))
! 1343: sp->count--;
! 1344: if (sp->synmax > ATHR && sp->synsnr > ASNR) {
! 1345: if (abs(epoch) < AWND * MS) {
! 1346: sp->reach |= 1;
! 1347: sp->count++;
! 1348: sp->mepoch = sp->lastpos = sp->pos;
! 1349: } else if (sp->count == 1) {
! 1350: sp->lastpos = sp->pos;
! 1351: }
! 1352: }
! 1353: if (up->watch > ACQSN)
! 1354: sp->metric = 0;
! 1355: else
! 1356: sp->metric = wwv_metric(sp);
! 1357: if (pp->sloppyclockflag & CLK_FLAG4) {
! 1358: sprintf(tbuf,
! 1359: "wwv8 %04x %3d %s %04x %.0f %.0f/%.1f %ld %ld",
! 1360: up->status, up->gain, sp->refid,
! 1361: sp->reach & 0xffff, sp->metric, sp->synmax,
! 1362: sp->synsnr, sp->pos % SECOND, epoch);
! 1363: record_clock_stats(&peer->srcadr, tbuf);
! 1364: #ifdef DEBUG
! 1365: if (debug)
! 1366: printf("%s\n", tbuf);
! 1367: #endif /* DEBUG */
! 1368: }
! 1369: sp->maxeng = sp->noieng = 0;
! 1370: }
! 1371: }
! 1372:
! 1373:
! 1374: /*
! 1375: * wwv_endpoc - identify and acquire second sync pulse
! 1376: *
! 1377: * This routine is called at the end of the second sync interval. It
! 1378: * determines the second sync epoch position within the second and
! 1379: * disciplines the sample clock using a frequency-lock loop (FLL).
! 1380: *
! 1381: * Second sync is determined in the RF input routine as the maximum
! 1382: * over all 8000 samples in the second comb filter. To assure accurate
! 1383: * and reliable time and frequency discipline, this routine performs a
! 1384: * great deal of heavy-handed heuristic data filtering and grooming.
! 1385: */
! 1386: static void
! 1387: wwv_endpoc(
! 1388: struct peer *peer, /* peer structure pointer */
! 1389: int epopos /* epoch max position */
! 1390: )
! 1391: {
! 1392: struct refclockproc *pp;
! 1393: struct wwvunit *up;
! 1394: static int epoch_mf[3]; /* epoch median filter */
! 1395: static int tepoch; /* current second epoch */
! 1396: static int xepoch; /* last second epoch */
! 1397: static int zepoch; /* last run epoch */
! 1398: static int zcount; /* last run end time */
! 1399: static int scount; /* seconds counter */
! 1400: static int syncnt; /* run length counter */
! 1401: static int maxrun; /* longest run length */
! 1402: static int mepoch; /* longest run end epoch */
! 1403: static int mcount; /* longest run end time */
! 1404: static int avgcnt; /* averaging interval counter */
! 1405: static int avginc; /* averaging ratchet */
! 1406: static int iniflg; /* initialization flag */
! 1407: char tbuf[TBUF]; /* monitor buffer */
! 1408: double dtemp;
! 1409: int tmp2;
! 1410:
! 1411: pp = peer->procptr;
! 1412: up = (struct wwvunit *)pp->unitptr;
! 1413: if (!iniflg) {
! 1414: iniflg = 1;
! 1415: memset((char *)epoch_mf, 0, sizeof(epoch_mf));
! 1416: }
! 1417:
! 1418: /*
! 1419: * If the signal amplitude or SNR fall below thresholds, dim the
! 1420: * second sync lamp and wait for hotter ions. If no stations are
! 1421: * heard, we are either in a probe cycle or the ions are really
! 1422: * cold.
! 1423: */
! 1424: scount++;
! 1425: if (up->epomax < STHR || up->eposnr < SSNR) {
! 1426: up->status &= ~(SSYNC | FGATE);
! 1427: avgcnt = syncnt = maxrun = 0;
! 1428: return;
! 1429: }
! 1430: if (!(up->status & (SELV | SELH)))
! 1431: return;
! 1432:
! 1433: /*
! 1434: * A three-stage median filter is used to help denoise the
! 1435: * second sync pulse. The median sample becomes the candidate
! 1436: * epoch.
! 1437: */
! 1438: epoch_mf[2] = epoch_mf[1];
! 1439: epoch_mf[1] = epoch_mf[0];
! 1440: epoch_mf[0] = epopos;
! 1441: if (epoch_mf[0] > epoch_mf[1]) {
! 1442: if (epoch_mf[1] > epoch_mf[2])
! 1443: tepoch = epoch_mf[1]; /* 0 1 2 */
! 1444: else if (epoch_mf[2] > epoch_mf[0])
! 1445: tepoch = epoch_mf[0]; /* 2 0 1 */
! 1446: else
! 1447: tepoch = epoch_mf[2]; /* 0 2 1 */
! 1448: } else {
! 1449: if (epoch_mf[1] < epoch_mf[2])
! 1450: tepoch = epoch_mf[1]; /* 2 1 0 */
! 1451: else if (epoch_mf[2] < epoch_mf[0])
! 1452: tepoch = epoch_mf[0]; /* 1 0 2 */
! 1453: else
! 1454: tepoch = epoch_mf[2]; /* 1 2 0 */
! 1455: }
! 1456:
! 1457:
! 1458: /*
! 1459: * If the epoch candidate is the same as the last one, increment
! 1460: * the run counter. If not, save the length, epoch and end
! 1461: * time of the current run for use later and reset the counter.
! 1462: * The epoch is considered valid if the run is at least SCMP
! 1463: * (10) s, the minute is synchronized and the interval since the
! 1464: * last epoch is not greater than the averaging interval. Thus,
! 1465: * after a long absence, the program will wait a full averaging
! 1466: * interval while the comb filter charges up and noise
! 1467: * dissapates..
! 1468: */
! 1469: tmp2 = (tepoch - xepoch) % SECOND;
! 1470: if (tmp2 == 0) {
! 1471: syncnt++;
! 1472: if (syncnt > SCMP && up->status & MSYNC && (up->status &
! 1473: FGATE || scount - zcount <= up->avgint)) {
! 1474: up->status |= SSYNC;
! 1475: up->yepoch = tepoch;
! 1476: }
! 1477: } else if (syncnt >= maxrun) {
! 1478: maxrun = syncnt;
! 1479: mcount = scount;
! 1480: mepoch = xepoch;
! 1481: syncnt = 0;
! 1482: }
! 1483: if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status &
! 1484: MSYNC)) {
! 1485: sprintf(tbuf,
! 1486: "wwv1 %04x %3d %4d %5.0f %5.1f %5d %4d %4d %4d",
! 1487: up->status, up->gain, tepoch, up->epomax,
! 1488: up->eposnr, tmp2, avgcnt, syncnt,
! 1489: maxrun);
! 1490: record_clock_stats(&peer->srcadr, tbuf);
! 1491: #ifdef DEBUG
! 1492: if (debug)
! 1493: printf("%s\n", tbuf);
! 1494: #endif /* DEBUG */
! 1495: }
! 1496: avgcnt++;
! 1497: if (avgcnt < up->avgint) {
! 1498: xepoch = tepoch;
! 1499: return;
! 1500: }
! 1501:
! 1502: /*
! 1503: * The sample clock frequency is disciplined using a first-order
! 1504: * feedback loop with time constant consistent with the Allan
! 1505: * intercept of typical computer clocks. During each averaging
! 1506: * interval the candidate epoch at the end of the longest run is
! 1507: * determined. If the longest run is zero, all epoches in the
! 1508: * interval are different, so the candidate epoch is the current
! 1509: * epoch. The frequency update is computed from the candidate
! 1510: * epoch difference (125-us units) and time difference (seconds)
! 1511: * between updates.
! 1512: */
! 1513: if (syncnt >= maxrun) {
! 1514: maxrun = syncnt;
! 1515: mcount = scount;
! 1516: mepoch = xepoch;
! 1517: }
! 1518: xepoch = tepoch;
! 1519: if (maxrun == 0) {
! 1520: mepoch = tepoch;
! 1521: mcount = scount;
! 1522: }
! 1523:
! 1524: /*
! 1525: * The master clock runs at the codec sample frequency of 8000
! 1526: * Hz, so the intrinsic time resolution is 125 us. The frequency
! 1527: * resolution ranges from 18 PPM at the minimum averaging
! 1528: * interval of 8 s to 0.12 PPM at the maximum interval of 1024
! 1529: * s. An offset update is determined at the end of the longest
! 1530: * run in each averaging interval. The frequency adjustment is
! 1531: * computed from the difference between offset updates and the
! 1532: * interval between them.
! 1533: *
! 1534: * The maximum frequency adjustment ranges from 187 PPM at the
! 1535: * minimum interval to 1.5 PPM at the maximum. If the adjustment
! 1536: * exceeds the maximum, the update is discarded and the
! 1537: * hysteresis counter is decremented. Otherwise, the frequency
! 1538: * is incremented by the adjustment, but clamped to the maximum
! 1539: * 187.5 PPM. If the update is less than half the maximum, the
! 1540: * hysteresis counter is incremented. If the counter increments
! 1541: * to +3, the averaging interval is doubled and the counter set
! 1542: * to zero; if it decrements to -3, the interval is halved and
! 1543: * the counter set to zero.
! 1544: */
! 1545: dtemp = (mepoch - zepoch) % SECOND;
! 1546: if (up->status & FGATE) {
! 1547: if (abs(dtemp) < MAXFREQ * MINAVG) {
! 1548: up->freq += (dtemp / 2.) / ((mcount - zcount) *
! 1549: FCONST);
! 1550: if (up->freq > MAXFREQ)
! 1551: up->freq = MAXFREQ;
! 1552: else if (up->freq < -MAXFREQ)
! 1553: up->freq = -MAXFREQ;
! 1554: if (abs(dtemp) < MAXFREQ * MINAVG / 2.) {
! 1555: if (avginc < 3) {
! 1556: avginc++;
! 1557: } else {
! 1558: if (up->avgint < MAXAVG) {
! 1559: up->avgint <<= 1;
! 1560: avginc = 0;
! 1561: }
! 1562: }
! 1563: }
! 1564: } else {
! 1565: if (avginc > -3) {
! 1566: avginc--;
! 1567: } else {
! 1568: if (up->avgint > MINAVG) {
! 1569: up->avgint >>= 1;
! 1570: avginc = 0;
! 1571: }
! 1572: }
! 1573: }
! 1574: }
! 1575: if (pp->sloppyclockflag & CLK_FLAG4) {
! 1576: sprintf(tbuf,
! 1577: "wwv2 %04x %5.0f %5.1f %5d %4d %4d %4d %4.0f %7.2f",
! 1578: up->status, up->epomax, up->eposnr, mepoch,
! 1579: up->avgint, maxrun, mcount - zcount, dtemp,
! 1580: up->freq * 1e6 / SECOND);
! 1581: record_clock_stats(&peer->srcadr, tbuf);
! 1582: #ifdef DEBUG
! 1583: if (debug)
! 1584: printf("%s\n", tbuf);
! 1585: #endif /* DEBUG */
! 1586: }
! 1587:
! 1588: /*
! 1589: * This is a valid update; set up for the next interval.
! 1590: */
! 1591: up->status |= FGATE;
! 1592: zepoch = mepoch;
! 1593: zcount = mcount;
! 1594: avgcnt = syncnt = maxrun = 0;
! 1595: }
! 1596:
! 1597:
! 1598: /*
! 1599: * wwv_epoch - epoch scanner
! 1600: *
! 1601: * This routine extracts data signals from the 100-Hz subcarrier. It
! 1602: * scans the receiver second epoch to determine the signal amplitudes
! 1603: * and pulse timings. Receiver synchronization is determined by the
! 1604: * minute sync pulse detected in the wwv_rf() routine and the second
! 1605: * sync pulse detected in the wwv_epoch() routine. The transmitted
! 1606: * signals are delayed by the propagation delay, receiver delay and
! 1607: * filter delay of this program. Delay corrections are introduced
! 1608: * separately for WWV and WWVH.
! 1609: *
! 1610: * Most communications radios use a highpass filter in the audio stages,
! 1611: * which can do nasty things to the subcarrier phase relative to the
! 1612: * sync pulses. Therefore, the data subcarrier reference phase is
! 1613: * disciplined using the hardlimited quadrature-phase signal sampled at
! 1614: * the same time as the in-phase signal. The phase tracking loop uses
! 1615: * phase adjustments of plus-minus one sample (125 us).
! 1616: */
! 1617: static void
! 1618: wwv_epoch(
! 1619: struct peer *peer /* peer structure pointer */
! 1620: )
! 1621: {
! 1622: struct refclockproc *pp;
! 1623: struct wwvunit *up;
! 1624: struct chan *cp;
! 1625: static double sigmin, sigzer, sigone, engmax, engmin;
! 1626:
! 1627: pp = peer->procptr;
! 1628: up = (struct wwvunit *)pp->unitptr;
! 1629:
! 1630: /*
! 1631: * Find the maximum minute sync pulse energy for both the
! 1632: * WWV and WWVH stations. This will be used later for channel
! 1633: * and station mitigation. Also set the seconds epoch at 800 ms
! 1634: * well before the end of the second to make sure we never set
! 1635: * the epoch backwards.
! 1636: */
! 1637: cp = &up->mitig[up->achan];
! 1638: if (cp->wwv.amp > cp->wwv.syneng)
! 1639: cp->wwv.syneng = cp->wwv.amp;
! 1640: if (cp->wwvh.amp > cp->wwvh.syneng)
! 1641: cp->wwvh.syneng = cp->wwvh.amp;
! 1642: if (up->rphase == 800 * MS)
! 1643: up->repoch = up->yepoch;
! 1644:
! 1645: /*
! 1646: * Use the signal amplitude at epoch 15 ms as the noise floor.
! 1647: * This gives a guard time of +-15 ms from the beginning of the
! 1648: * second until the second pulse rises at 30 ms. There is a
! 1649: * compromise here; we want to delay the sample as long as
! 1650: * possible to give the radio time to change frequency and the
! 1651: * AGC to stabilize, but as early as possible if the second
! 1652: * epoch is not exact.
! 1653: */
! 1654: if (up->rphase == 15 * MS)
! 1655: sigmin = sigzer = sigone = up->irig;
! 1656:
! 1657: /*
! 1658: * Latch the data signal at 200 ms. Keep this around until the
! 1659: * end of the second. Use the signal energy as the peak to
! 1660: * compute the SNR. Use the Q sample to adjust the 100-Hz
! 1661: * reference oscillator phase.
! 1662: */
! 1663: if (up->rphase == 200 * MS) {
! 1664: sigzer = up->irig;
! 1665: engmax = sqrt(up->irig * up->irig + up->qrig *
! 1666: up->qrig);
! 1667: up->datpha = up->qrig / up->avgint;
! 1668: if (up->datpha >= 0) {
! 1669: up->datapt++;
! 1670: if (up->datapt >= 80)
! 1671: up->datapt -= 80;
! 1672: } else {
! 1673: up->datapt--;
! 1674: if (up->datapt < 0)
! 1675: up->datapt += 80;
! 1676: }
! 1677: }
! 1678:
! 1679:
! 1680: /*
! 1681: * Latch the data signal at 500 ms. Keep this around until the
! 1682: * end of the second.
! 1683: */
! 1684: else if (up->rphase == 500 * MS)
! 1685: sigone = up->irig;
! 1686:
! 1687: /*
! 1688: * At the end of the second crank the clock state machine and
! 1689: * adjust the codec gain. Note the epoch is buffered from the
! 1690: * center of the second in order to avoid jitter while the
! 1691: * seconds synch is diddling the epoch. Then, determine the true
! 1692: * offset and update the median filter in the driver interface.
! 1693: *
! 1694: * Use the energy at the end of the second as the noise to
! 1695: * compute the SNR for the data pulse. This gives a better
! 1696: * measurement than the beginning of the second, especially when
! 1697: * returning from the probe channel. This gives a guard time of
! 1698: * 30 ms from the decay of the longest pulse to the rise of the
! 1699: * next pulse.
! 1700: */
! 1701: up->rphase++;
! 1702: if (up->mphase % SECOND == up->repoch) {
! 1703: up->status &= ~(DGATE | BGATE);
! 1704: engmin = sqrt(up->irig * up->irig + up->qrig *
! 1705: up->qrig);
! 1706: up->datsig = engmax;
! 1707: up->datsnr = wwv_snr(engmax, engmin);
! 1708:
! 1709: /*
! 1710: * If the amplitude or SNR is below threshold, average a
! 1711: * 0 in the the integrators; otherwise, average the
! 1712: * bipolar signal. This is done to avoid noise polution.
! 1713: */
! 1714: if (engmax < DTHR || up->datsnr < DSNR) {
! 1715: up->status |= DGATE;
! 1716: wwv_rsec(peer, 0);
! 1717: } else {
! 1718: sigzer -= sigone;
! 1719: sigone -= sigmin;
! 1720: wwv_rsec(peer, sigone - sigzer);
! 1721: }
! 1722: if (up->status & (DGATE | BGATE))
! 1723: up->errcnt++;
! 1724: if (up->errcnt > MAXERR)
! 1725: up->alarm |= LOWERR;
! 1726: wwv_gain(peer);
! 1727: cp = &up->mitig[up->achan];
! 1728: cp->wwv.syneng = 0;
! 1729: cp->wwvh.syneng = 0;
! 1730: up->rphase = 0;
! 1731: }
! 1732: }
! 1733:
! 1734:
! 1735: /*
! 1736: * wwv_rsec - process receiver second
! 1737: *
! 1738: * This routine is called at the end of each receiver second to
! 1739: * implement the per-second state machine. The machine assembles BCD
! 1740: * digit bits, decodes miscellaneous bits and dances the leap seconds.
! 1741: *
! 1742: * Normally, the minute has 60 seconds numbered 0-59. If the leap
! 1743: * warning bit is set, the last minute (1439) of 30 June (day 181 or 182
! 1744: * for leap years) or 31 December (day 365 or 366 for leap years) is
! 1745: * augmented by one second numbered 60. This is accomplished by
! 1746: * extending the minute interval by one second and teaching the state
! 1747: * machine to ignore it.
! 1748: */
! 1749: static void
! 1750: wwv_rsec(
! 1751: struct peer *peer, /* peer structure pointer */
! 1752: double bit
! 1753: )
! 1754: {
! 1755: static int iniflg; /* initialization flag */
! 1756: static double bcddld[4]; /* BCD data bits */
! 1757: static double bitvec[61]; /* bit integrator for misc bits */
! 1758: struct refclockproc *pp;
! 1759: struct wwvunit *up;
! 1760: struct chan *cp;
! 1761: struct sync *sp, *rp;
! 1762: char tbuf[TBUF]; /* monitor buffer */
! 1763: int sw, arg, nsec;
! 1764:
! 1765: pp = peer->procptr;
! 1766: up = (struct wwvunit *)pp->unitptr;
! 1767: if (!iniflg) {
! 1768: iniflg = 1;
! 1769: memset((char *)bitvec, 0, sizeof(bitvec));
! 1770: }
! 1771:
! 1772: /*
! 1773: * The bit represents the probability of a hit on zero (negative
! 1774: * values), a hit on one (positive values) or a miss (zero
! 1775: * value). The likelihood vector is the exponential average of
! 1776: * these probabilities. Only the bits of this vector
! 1777: * corresponding to the miscellaneous bits of the timecode are
! 1778: * used, but it's easier to do them all. After that, crank the
! 1779: * seconds state machine.
! 1780: */
! 1781: nsec = up->rsec;
! 1782: up->rsec++;
! 1783: bitvec[nsec] += (bit - bitvec[nsec]) / TCONST;
! 1784: sw = progx[nsec].sw;
! 1785: arg = progx[nsec].arg;
! 1786:
! 1787: /*
! 1788: * The minute state machine. Fly off to a particular section as
! 1789: * directed by the transition matrix and second number.
! 1790: */
! 1791: switch (sw) {
! 1792:
! 1793: /*
! 1794: * Ignore this second.
! 1795: */
! 1796: case IDLE: /* 9, 45-49 */
! 1797: break;
! 1798:
! 1799: /*
! 1800: * Probe channel stuff
! 1801: *
! 1802: * The WWV/H format contains data pulses in second 59 (position
! 1803: * identifier) and second 1, but not in second 0. The minute
! 1804: * sync pulse is contained in second 0. At the end of second 58
! 1805: * QSY to the probe channel, which rotates in turn over all
! 1806: * WWV/H frequencies. At the end of second 0 measure the minute
! 1807: * sync pulse. At the end of second 1 measure the data pulse and
! 1808: * QSY back to the data channel. Note that the actions commented
! 1809: * here happen at the end of the second numbered as shown.
! 1810: *
! 1811: * At the end of second 0 save the minute sync amplitude latched
! 1812: * at 800 ms as the signal later used to calculate the SNR.
! 1813: */
! 1814: case SYNC2: /* 0 */
! 1815: cp = &up->mitig[up->achan];
! 1816: cp->wwv.synmax = cp->wwv.syneng;
! 1817: cp->wwvh.synmax = cp->wwvh.syneng;
! 1818: break;
! 1819:
! 1820: /*
! 1821: * At the end of second 1 use the minute sync amplitude latched
! 1822: * at 800 ms as the noise to calculate the SNR. If the minute
! 1823: * sync pulse and SNR are above thresholds and the data pulse
! 1824: * amplitude and SNR are above thresolds, shift a 1 into the
! 1825: * station reachability register; otherwise, shift a 0. The
! 1826: * number of 1 bits in the last six intervals is a component of
! 1827: * the channel metric computed by the wwv_metric() routine.
! 1828: * Finally, QSY back to the data channel.
! 1829: */
! 1830: case SYNC3: /* 1 */
! 1831: cp = &up->mitig[up->achan];
! 1832:
! 1833: /*
! 1834: * WWV station
! 1835: */
! 1836: sp = &cp->wwv;
! 1837: sp->synsnr = wwv_snr(sp->synmax, sp->amp);
! 1838: sp->reach <<= 1;
! 1839: if (sp->reach & (1 << AMAX))
! 1840: sp->count--;
! 1841: if (sp->synmax >= QTHR && sp->synsnr >= QSNR &&
! 1842: !(up->status & (DGATE | BGATE))) {
! 1843: sp->reach |= 1;
! 1844: sp->count++;
! 1845: }
! 1846: sp->metric = wwv_metric(sp);
! 1847:
! 1848: /*
! 1849: * WWVH station
! 1850: */
! 1851: rp = &cp->wwvh;
! 1852: rp->synsnr = wwv_snr(rp->synmax, rp->amp);
! 1853: rp->reach <<= 1;
! 1854: if (rp->reach & (1 << AMAX))
! 1855: rp->count--;
! 1856: if (rp->synmax >= QTHR && rp->synsnr >= QSNR &&
! 1857: !(up->status & (DGATE | BGATE))) {
! 1858: rp->reach |= 1;
! 1859: rp->count++;
! 1860: }
! 1861: rp->metric = wwv_metric(rp);
! 1862: if (pp->sloppyclockflag & CLK_FLAG4) {
! 1863: sprintf(tbuf,
! 1864: "wwv5 %04x %3d %4d %.0f/%.1f %.0f/%.1f %s %04x %.0f %.0f/%.1f %s %04x %.0f %.0f/%.1f",
! 1865: up->status, up->gain, up->yepoch,
! 1866: up->epomax, up->eposnr, up->datsig,
! 1867: up->datsnr,
! 1868: sp->refid, sp->reach & 0xffff,
! 1869: sp->metric, sp->synmax, sp->synsnr,
! 1870: rp->refid, rp->reach & 0xffff,
! 1871: rp->metric, rp->synmax, rp->synsnr);
! 1872: record_clock_stats(&peer->srcadr, tbuf);
! 1873: #ifdef DEBUG
! 1874: if (debug)
! 1875: printf("%s\n", tbuf);
! 1876: #endif /* DEBUG */
! 1877: }
! 1878: up->errcnt = up->digcnt = up->alarm = 0;
! 1879:
! 1880: /*
! 1881: * If synchronized to a station, restart if no stations
! 1882: * have been heard within the PANIC timeout (2 days). If
! 1883: * not and the minute digit has been found, restart if
! 1884: * not synchronized withing the SYNCH timeout (40 m). If
! 1885: * not, restart if the unit digit has not been found
! 1886: * within the DATA timeout (15 m).
! 1887: */
! 1888: if (up->status & INSYNC) {
! 1889: if (up->watch > PANIC) {
! 1890: wwv_newgame(peer);
! 1891: return;
! 1892: }
! 1893: } else if (up->status & DSYNC) {
! 1894: if (up->watch > SYNCH) {
! 1895: wwv_newgame(peer);
! 1896: return;
! 1897: }
! 1898: } else if (up->watch > DATA) {
! 1899: wwv_newgame(peer);
! 1900: return;
! 1901: }
! 1902: wwv_newchan(peer);
! 1903: break;
! 1904:
! 1905: /*
! 1906: * Save the bit probability in the BCD data vector at the index
! 1907: * given by the argument. Bits not used in the digit are forced
! 1908: * to zero.
! 1909: */
! 1910: case COEF1: /* 4-7 */
! 1911: bcddld[arg] = bit;
! 1912: break;
! 1913:
! 1914: case COEF: /* 10-13, 15-17, 20-23, 25-26,
! 1915: 30-33, 35-38, 40-41, 51-54 */
! 1916: if (up->status & DSYNC)
! 1917: bcddld[arg] = bit;
! 1918: else
! 1919: bcddld[arg] = 0;
! 1920: break;
! 1921:
! 1922: case COEF2: /* 18, 27-28, 42-43 */
! 1923: bcddld[arg] = 0;
! 1924: break;
! 1925:
! 1926: /*
! 1927: * Correlate coefficient vector with each valid digit vector and
! 1928: * save in decoding matrix. We step through the decoding matrix
! 1929: * digits correlating each with the coefficients and saving the
! 1930: * greatest and the next lower for later SNR calculation.
! 1931: */
! 1932: case DECIM2: /* 29 */
! 1933: wwv_corr4(peer, &up->decvec[arg], bcddld, bcd2);
! 1934: break;
! 1935:
! 1936: case DECIM3: /* 44 */
! 1937: wwv_corr4(peer, &up->decvec[arg], bcddld, bcd3);
! 1938: break;
! 1939:
! 1940: case DECIM6: /* 19 */
! 1941: wwv_corr4(peer, &up->decvec[arg], bcddld, bcd6);
! 1942: break;
! 1943:
! 1944: case DECIM9: /* 8, 14, 24, 34, 39 */
! 1945: wwv_corr4(peer, &up->decvec[arg], bcddld, bcd9);
! 1946: break;
! 1947:
! 1948: /*
! 1949: * Miscellaneous bits. If above the positive threshold, declare
! 1950: * 1; if below the negative threshold, declare 0; otherwise
! 1951: * raise the BGATE bit. The design is intended to avoid
! 1952: * integrating noise under low SNR conditions.
! 1953: */
! 1954: case MSC20: /* 55 */
! 1955: wwv_corr4(peer, &up->decvec[YR + 1], bcddld, bcd9);
! 1956: /* fall through */
! 1957:
! 1958: case MSCBIT: /* 2-3, 50, 56-57 */
! 1959: if (bitvec[nsec] > BTHR) {
! 1960: if (!(up->misc & arg))
! 1961: up->alarm |= CMPERR;
! 1962: up->misc |= arg;
! 1963: } else if (bitvec[nsec] < -BTHR) {
! 1964: if (up->misc & arg)
! 1965: up->alarm |= CMPERR;
! 1966: up->misc &= ~arg;
! 1967: } else {
! 1968: up->status |= BGATE;
! 1969: }
! 1970: break;
! 1971:
! 1972: /*
! 1973: * Save the data channel gain, then QSY to the probe channel and
! 1974: * dim the seconds comb filters. The www_newchan() routine will
! 1975: * light them back up.
! 1976: */
! 1977: case MSC21: /* 58 */
! 1978: if (bitvec[nsec] > BTHR) {
! 1979: if (!(up->misc & arg))
! 1980: up->alarm |= CMPERR;
! 1981: up->misc |= arg;
! 1982: } else if (bitvec[nsec] < -BTHR) {
! 1983: if (up->misc & arg)
! 1984: up->alarm |= CMPERR;
! 1985: up->misc &= ~arg;
! 1986: } else {
! 1987: up->status |= BGATE;
! 1988: }
! 1989: up->status &= ~(SELV | SELH);
! 1990: #ifdef ICOM
! 1991: if (up->fd_icom > 0) {
! 1992: up->schan = (up->schan + 1) % NCHAN;
! 1993: wwv_qsy(peer, up->schan);
! 1994: } else {
! 1995: up->mitig[up->achan].gain = up->gain;
! 1996: }
! 1997: #else
! 1998: up->mitig[up->achan].gain = up->gain;
! 1999: #endif /* ICOM */
! 2000: break;
! 2001:
! 2002: /*
! 2003: * The endgames
! 2004: *
! 2005: * During second 59 the receiver and codec AGC are settling
! 2006: * down, so the data pulse is unusable as quality metric. If
! 2007: * LEPSEC is set on the last minute of 30 June or 31 December,
! 2008: * the transmitter and receiver insert an extra second (60) in
! 2009: * the timescale and the minute sync repeats the second. Once
! 2010: * leaps occurred at intervals of about 18 months, but the last
! 2011: * leap before the most recent leap in 1995 was in 1998.
! 2012: */
! 2013: case MIN1: /* 59 */
! 2014: if (up->status & LEPSEC)
! 2015: break;
! 2016:
! 2017: /* fall through */
! 2018:
! 2019: case MIN2: /* 60 */
! 2020: up->status &= ~LEPSEC;
! 2021: wwv_tsec(peer);
! 2022: up->rsec = 0;
! 2023: wwv_clock(peer);
! 2024: break;
! 2025: }
! 2026: if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status &
! 2027: DSYNC)) {
! 2028: sprintf(tbuf,
! 2029: "wwv3 %2d %04x %3d %4d %5.0f %5.1f %5.0f %5.1f %5.0f",
! 2030: nsec, up->status, up->gain, up->yepoch, up->epomax,
! 2031: up->eposnr, up->datsig, up->datsnr, bit);
! 2032: record_clock_stats(&peer->srcadr, tbuf);
! 2033: #ifdef DEBUG
! 2034: if (debug)
! 2035: printf("%s\n", tbuf);
! 2036: #endif /* DEBUG */
! 2037: }
! 2038: pp->disp += AUDIO_PHI;
! 2039: }
! 2040:
! 2041: /*
! 2042: * The radio clock is set if the alarm bits are all zero. After that,
! 2043: * the time is considered valid if the second sync bit is lit. It should
! 2044: * not be a surprise, especially if the radio is not tunable, that
! 2045: * sometimes no stations are above the noise and the integrators
! 2046: * discharge below the thresholds. We assume that, after a day of signal
! 2047: * loss, the minute sync epoch will be in the same second. This requires
! 2048: * the codec frequency be accurate within 6 PPM. Practical experience
! 2049: * shows the frequency typically within 0.1 PPM, so after a day of
! 2050: * signal loss, the time should be within 8.6 ms..
! 2051: */
! 2052: static void
! 2053: wwv_clock(
! 2054: struct peer *peer /* peer unit pointer */
! 2055: )
! 2056: {
! 2057: struct refclockproc *pp;
! 2058: struct wwvunit *up;
! 2059: l_fp offset; /* offset in NTP seconds */
! 2060:
! 2061: pp = peer->procptr;
! 2062: up = (struct wwvunit *)pp->unitptr;
! 2063: if (!(up->status & SSYNC))
! 2064: up->alarm |= SYNERR;
! 2065: if (up->digcnt < 9)
! 2066: up->alarm |= NINERR;
! 2067: if (!(up->alarm))
! 2068: up->status |= INSYNC;
! 2069: if (up->status & INSYNC && up->status & SSYNC) {
! 2070: if (up->misc & SECWAR)
! 2071: pp->leap = LEAP_ADDSECOND;
! 2072: else
! 2073: pp->leap = LEAP_NOWARNING;
! 2074: pp->second = up->rsec;
! 2075: pp->minute = up->decvec[MN].digit + up->decvec[MN +
! 2076: 1].digit * 10;
! 2077: pp->hour = up->decvec[HR].digit + up->decvec[HR +
! 2078: 1].digit * 10;
! 2079: pp->day = up->decvec[DA].digit + up->decvec[DA +
! 2080: 1].digit * 10 + up->decvec[DA + 2].digit * 100;
! 2081: pp->year = up->decvec[YR].digit + up->decvec[YR +
! 2082: 1].digit * 10;
! 2083: pp->year += 2000;
! 2084: L_CLR(&offset);
! 2085: if (!clocktime(pp->day, pp->hour, pp->minute,
! 2086: pp->second, GMT, up->timestamp.l_ui,
! 2087: &pp->yearstart, &offset.l_ui)) {
! 2088: up->errflg = CEVNT_BADTIME;
! 2089: } else {
! 2090: up->watch = 0;
! 2091: pp->disp = 0;
! 2092: pp->lastref = up->timestamp;
! 2093: refclock_process_offset(pp, offset,
! 2094: up->timestamp, PDELAY + up->pdelay);
! 2095: refclock_receive(peer);
! 2096: }
! 2097: }
! 2098: pp->lencode = timecode(up, pp->a_lastcode);
! 2099: record_clock_stats(&peer->srcadr, pp->a_lastcode);
! 2100: #ifdef DEBUG
! 2101: if (debug)
! 2102: printf("wwv: timecode %d %s\n", pp->lencode,
! 2103: pp->a_lastcode);
! 2104: #endif /* DEBUG */
! 2105: }
! 2106:
! 2107:
! 2108: /*
! 2109: * wwv_corr4 - determine maximum-likelihood digit
! 2110: *
! 2111: * This routine correlates the received digit vector with the BCD
! 2112: * coefficient vectors corresponding to all valid digits at the given
! 2113: * position in the decoding matrix. The maximum value corresponds to the
! 2114: * maximum-likelihood digit, while the ratio of this value to the next
! 2115: * lower value determines the likelihood function. Note that, if the
! 2116: * digit is invalid, the likelihood vector is averaged toward a miss.
! 2117: */
! 2118: static void
! 2119: wwv_corr4(
! 2120: struct peer *peer, /* peer unit pointer */
! 2121: struct decvec *vp, /* decoding table pointer */
! 2122: double data[], /* received data vector */
! 2123: double tab[][4] /* correlation vector array */
! 2124: )
! 2125: {
! 2126: struct refclockproc *pp;
! 2127: struct wwvunit *up;
! 2128: double topmax, nxtmax; /* metrics */
! 2129: double acc; /* accumulator */
! 2130: char tbuf[TBUF]; /* monitor buffer */
! 2131: int mldigit; /* max likelihood digit */
! 2132: int i, j;
! 2133:
! 2134: pp = peer->procptr;
! 2135: up = (struct wwvunit *)pp->unitptr;
! 2136:
! 2137: /*
! 2138: * Correlate digit vector with each BCD coefficient vector. If
! 2139: * any BCD digit bit is bad, consider all bits a miss. Until the
! 2140: * minute units digit has been resolved, don't to anything else.
! 2141: * Note the SNR is calculated as the ratio of the largest
! 2142: * likelihood value to the next largest likelihood value.
! 2143: */
! 2144: mldigit = 0;
! 2145: topmax = nxtmax = -MAXAMP;
! 2146: for (i = 0; tab[i][0] != 0; i++) {
! 2147: acc = 0;
! 2148: for (j = 0; j < 4; j++)
! 2149: acc += data[j] * tab[i][j];
! 2150: acc = (vp->like[i] += (acc - vp->like[i]) / TCONST);
! 2151: if (acc > topmax) {
! 2152: nxtmax = topmax;
! 2153: topmax = acc;
! 2154: mldigit = i;
! 2155: } else if (acc > nxtmax) {
! 2156: nxtmax = acc;
! 2157: }
! 2158: }
! 2159: vp->digprb = topmax;
! 2160: vp->digsnr = wwv_snr(topmax, nxtmax);
! 2161:
! 2162: /*
! 2163: * The current maximum-likelihood digit is compared to the last
! 2164: * maximum-likelihood digit. If different, the compare counter
! 2165: * and maximum-likelihood digit are reset. When the compare
! 2166: * counter reaches the BCMP threshold (3), the digit is assumed
! 2167: * correct. When the compare counter of all nine digits have
! 2168: * reached threshold, the clock is assumed correct.
! 2169: *
! 2170: * Note that the clock display digit is set before the compare
! 2171: * counter has reached threshold; however, the clock display is
! 2172: * not considered correct until all nine clock digits have
! 2173: * reached threshold. This is intended as eye candy, but avoids
! 2174: * mistakes when the signal is low and the SNR is very marginal.
! 2175: */
! 2176: if (vp->digprb < BTHR || vp->digsnr < BSNR) {
! 2177: up->status |= BGATE;
! 2178: } else {
! 2179: if (vp->digit != mldigit) {
! 2180: up->alarm |= CMPERR;
! 2181: if (vp->count > 0)
! 2182: vp->count--;
! 2183: if (vp->count == 0)
! 2184: vp->digit = mldigit;
! 2185: } else {
! 2186: if (vp->count < BCMP)
! 2187: vp->count++;
! 2188: if (vp->count == BCMP) {
! 2189: up->status |= DSYNC;
! 2190: up->digcnt++;
! 2191: }
! 2192: }
! 2193: }
! 2194: if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status &
! 2195: INSYNC)) {
! 2196: sprintf(tbuf,
! 2197: "wwv4 %2d %04x %3d %4d %5.0f %2d %d %d %d %5.0f %5.1f",
! 2198: up->rsec - 1, up->status, up->gain, up->yepoch,
! 2199: up->epomax, vp->radix, vp->digit, mldigit,
! 2200: vp->count, vp->digprb, vp->digsnr);
! 2201: record_clock_stats(&peer->srcadr, tbuf);
! 2202: #ifdef DEBUG
! 2203: if (debug)
! 2204: printf("%s\n", tbuf);
! 2205: #endif /* DEBUG */
! 2206: }
! 2207: }
! 2208:
! 2209:
! 2210: /*
! 2211: * wwv_tsec - transmitter minute processing
! 2212: *
! 2213: * This routine is called at the end of the transmitter minute. It
! 2214: * implements a state machine that advances the logical clock subject to
! 2215: * the funny rules that govern the conventional clock and calendar.
! 2216: */
! 2217: static void
! 2218: wwv_tsec(
! 2219: struct peer *peer /* driver structure pointer */
! 2220: )
! 2221: {
! 2222: struct refclockproc *pp;
! 2223: struct wwvunit *up;
! 2224: int minute, day, isleap;
! 2225: int temp;
! 2226:
! 2227: pp = peer->procptr;
! 2228: up = (struct wwvunit *)pp->unitptr;
! 2229:
! 2230: /*
! 2231: * Advance minute unit of the day. Don't propagate carries until
! 2232: * the unit minute digit has been found.
! 2233: */
! 2234: temp = carry(&up->decvec[MN]); /* minute units */
! 2235: if (!(up->status & DSYNC))
! 2236: return;
! 2237:
! 2238: /*
! 2239: * Propagate carries through the day.
! 2240: */
! 2241: if (temp == 0) /* carry minutes */
! 2242: temp = carry(&up->decvec[MN + 1]);
! 2243: if (temp == 0) /* carry hours */
! 2244: temp = carry(&up->decvec[HR]);
! 2245: if (temp == 0)
! 2246: temp = carry(&up->decvec[HR + 1]);
! 2247:
! 2248: /*
! 2249: * Decode the current minute and day. Set leap day if the
! 2250: * timecode leap bit is set on 30 June or 31 December. Set leap
! 2251: * minute if the last minute on leap day, but only if the clock
! 2252: * is syncrhronized. This code fails in 2400 AD.
! 2253: */
! 2254: minute = up->decvec[MN].digit + up->decvec[MN + 1].digit *
! 2255: 10 + up->decvec[HR].digit * 60 + up->decvec[HR +
! 2256: 1].digit * 600;
! 2257: day = up->decvec[DA].digit + up->decvec[DA + 1].digit * 10 +
! 2258: up->decvec[DA + 2].digit * 100;
! 2259:
! 2260: /*
! 2261: * Set the leap bit on the last minute of the leap day.
! 2262: */
! 2263: isleap = up->decvec[YR].digit & 0x3;
! 2264: if (up->misc & SECWAR && up->status & INSYNC) {
! 2265: if ((day == (isleap ? 182 : 183) || day == (isleap ?
! 2266: 365 : 366)) && minute == 1439)
! 2267: up->status |= LEPSEC;
! 2268: }
! 2269:
! 2270: /*
! 2271: * Roll the day if this the first minute and propagate carries
! 2272: * through the year.
! 2273: */
! 2274: if (minute != 1440)
! 2275: return;
! 2276:
! 2277: minute = 0;
! 2278: while (carry(&up->decvec[HR]) != 0); /* advance to minute 0 */
! 2279: while (carry(&up->decvec[HR + 1]) != 0);
! 2280: day++;
! 2281: temp = carry(&up->decvec[DA]); /* carry days */
! 2282: if (temp == 0)
! 2283: temp = carry(&up->decvec[DA + 1]);
! 2284: if (temp == 0)
! 2285: temp = carry(&up->decvec[DA + 2]);
! 2286:
! 2287: /*
! 2288: * Roll the year if this the first day and propagate carries
! 2289: * through the century.
! 2290: */
! 2291: if (day != (isleap ? 365 : 366))
! 2292: return;
! 2293:
! 2294: day = 1;
! 2295: while (carry(&up->decvec[DA]) != 1); /* advance to day 1 */
! 2296: while (carry(&up->decvec[DA + 1]) != 0);
! 2297: while (carry(&up->decvec[DA + 2]) != 0);
! 2298: temp = carry(&up->decvec[YR]); /* carry years */
! 2299: if (temp == 0)
! 2300: carry(&up->decvec[YR + 1]);
! 2301: }
! 2302:
! 2303:
! 2304: /*
! 2305: * carry - process digit
! 2306: *
! 2307: * This routine rotates a likelihood vector one position and increments
! 2308: * the clock digit modulo the radix. It returns the new clock digit or
! 2309: * zero if a carry occurred. Once synchronized, the clock digit will
! 2310: * match the maximum-likelihood digit corresponding to that position.
! 2311: */
! 2312: static int
! 2313: carry(
! 2314: struct decvec *dp /* decoding table pointer */
! 2315: )
! 2316: {
! 2317: int temp;
! 2318: int j;
! 2319:
! 2320: dp->digit++;
! 2321: if (dp->digit == dp->radix)
! 2322: dp->digit = 0;
! 2323: temp = dp->like[dp->radix - 1];
! 2324: for (j = dp->radix - 1; j > 0; j--)
! 2325: dp->like[j] = dp->like[j - 1];
! 2326: dp->like[0] = temp;
! 2327: return (dp->digit);
! 2328: }
! 2329:
! 2330:
! 2331: /*
! 2332: * wwv_snr - compute SNR or likelihood function
! 2333: */
! 2334: static double
! 2335: wwv_snr(
! 2336: double signal, /* signal */
! 2337: double noise /* noise */
! 2338: )
! 2339: {
! 2340: double rval;
! 2341:
! 2342: /*
! 2343: * This is a little tricky. Due to the way things are measured,
! 2344: * either or both the signal or noise amplitude can be negative
! 2345: * or zero. The intent is that, if the signal is negative or
! 2346: * zero, the SNR must always be zero. This can happen with the
! 2347: * subcarrier SNR before the phase has been aligned. On the
! 2348: * other hand, in the likelihood function the "noise" is the
! 2349: * next maximum down from the peak and this could be negative.
! 2350: * However, in this case the SNR is truly stupendous, so we
! 2351: * simply cap at MAXSNR dB (40).
! 2352: */
! 2353: if (signal <= 0) {
! 2354: rval = 0;
! 2355: } else if (noise <= 0) {
! 2356: rval = MAXSNR;
! 2357: } else {
! 2358: rval = 20. * log10(signal / noise);
! 2359: if (rval > MAXSNR)
! 2360: rval = MAXSNR;
! 2361: }
! 2362: return (rval);
! 2363: }
! 2364:
! 2365:
! 2366: /*
! 2367: * wwv_newchan - change to new data channel
! 2368: *
! 2369: * The radio actually appears to have ten channels, one channel for each
! 2370: * of five frequencies and each of two stations (WWV and WWVH), although
! 2371: * if not tunable only the DCHAN channel appears live. While the radio
! 2372: * is tuned to the working data channel frequency and station for most
! 2373: * of the minute, during seconds 59, 0 and 1 the radio is tuned to a
! 2374: * probe frequency in order to search for minute sync pulse and data
! 2375: * subcarrier from other transmitters.
! 2376: *
! 2377: * The search for WWV and WWVH operates simultaneously, with WWV minute
! 2378: * sync pulse at 1000 Hz and WWVH at 1200 Hz. The probe frequency
! 2379: * rotates each minute over 2.5, 5, 10, 15 and 20 MHz in order and yes,
! 2380: * we all know WWVH is dark on 20 MHz, but few remember when WWV was lit
! 2381: * on 25 MHz.
! 2382: *
! 2383: * This routine selects the best channel using a metric computed from
! 2384: * the reachability register and minute pulse amplitude. Normally, the
! 2385: * award goes to the the channel with the highest metric; but, in case
! 2386: * of ties, the award goes to the channel with the highest minute sync
! 2387: * pulse amplitude and then to the highest frequency.
! 2388: *
! 2389: * The routine performs an important squelch function to keep dirty data
! 2390: * from polluting the integrators. In order to consider a station valid,
! 2391: * the metric must be at least MTHR (13); otherwise, the station select
! 2392: * bits are cleared so the second sync is disabled and the data bit
! 2393: * integrators averaged to a miss.
! 2394: */
! 2395: static int
! 2396: wwv_newchan(
! 2397: struct peer *peer /* peer structure pointer */
! 2398: )
! 2399: {
! 2400: struct refclockproc *pp;
! 2401: struct wwvunit *up;
! 2402: struct sync *sp, *rp;
! 2403: double rank, dtemp;
! 2404: int i, j, rval;
! 2405:
! 2406: pp = peer->procptr;
! 2407: up = (struct wwvunit *)pp->unitptr;
! 2408:
! 2409: /*
! 2410: * Search all five station pairs looking for the channel with
! 2411: * maximum metric.
! 2412: */
! 2413: sp = NULL;
! 2414: j = 0;
! 2415: rank = 0;
! 2416: for (i = 0; i < NCHAN; i++) {
! 2417: rp = &up->mitig[i].wwvh;
! 2418: dtemp = rp->metric;
! 2419: if (dtemp >= rank) {
! 2420: rank = dtemp;
! 2421: sp = rp;
! 2422: j = i;
! 2423: }
! 2424: rp = &up->mitig[i].wwv;
! 2425: dtemp = rp->metric;
! 2426: if (dtemp >= rank) {
! 2427: rank = dtemp;
! 2428: sp = rp;
! 2429: j = i;
! 2430: }
! 2431: }
! 2432:
! 2433: /*
! 2434: * If the strongest signal is less than the MTHR threshold (13),
! 2435: * we are beneath the waves, so squelch the second sync and
! 2436: * advance to the next station. This makes sure all stations are
! 2437: * scanned when the ions grow dim. If the strongest signal is
! 2438: * greater than the threshold, tune to that frequency and
! 2439: * transmitter QTH.
! 2440: */
! 2441: up->status &= ~(SELV | SELH);
! 2442: if (rank < MTHR) {
! 2443: up->dchan = (up->dchan + 1) % NCHAN;
! 2444: if (up->status & METRIC) {
! 2445: up->status &= ~METRIC;
! 2446: refclock_report(peer, CEVNT_PROP);
! 2447: }
! 2448: rval = FALSE;
! 2449: } else {
! 2450: up->dchan = j;
! 2451: up->sptr = sp;
! 2452: memcpy(&pp->refid, sp->refid, 4);
! 2453: peer->refid = pp->refid;
! 2454: up->status |= METRIC;
! 2455: if (sp->select & SELV) {
! 2456: up->status |= SELV;
! 2457: up->pdelay = pp->fudgetime1;
! 2458: } else if (sp->select & SELH) {
! 2459: up->status |= SELH;
! 2460: up->pdelay = pp->fudgetime2;
! 2461: } else {
! 2462: up->pdelay = 0;
! 2463: }
! 2464: rval = TRUE;
! 2465: }
! 2466: #ifdef ICOM
! 2467: if (up->fd_icom > 0)
! 2468: wwv_qsy(peer, up->dchan);
! 2469: #endif /* ICOM */
! 2470: return (rval);
! 2471: }
! 2472:
! 2473:
! 2474: /*
! 2475: * wwv_newgame - reset and start over
! 2476: *
! 2477: * There are three conditions resulting in a new game:
! 2478: *
! 2479: * 1 After finding the minute pulse (MSYNC lit), going 15 minutes
! 2480: * (DATA) without finding the unit seconds digit.
! 2481: *
! 2482: * 2 After finding good data (DSYNC lit), going more than 40 minutes
! 2483: * (SYNCH) without finding station sync (INSYNC lit).
! 2484: *
! 2485: * 3 After finding station sync (INSYNC lit), going more than 2 days
! 2486: * (PANIC) without finding any station.
! 2487: */
! 2488: static void
! 2489: wwv_newgame(
! 2490: struct peer *peer /* peer structure pointer */
! 2491: )
! 2492: {
! 2493: struct refclockproc *pp;
! 2494: struct wwvunit *up;
! 2495: struct chan *cp;
! 2496: int i;
! 2497:
! 2498: pp = peer->procptr;
! 2499: up = (struct wwvunit *)pp->unitptr;
! 2500:
! 2501: /*
! 2502: * Initialize strategic values. Note we set the leap bits
! 2503: * NOTINSYNC and the refid "NONE".
! 2504: */
! 2505: if (up->status)
! 2506: up->errflg = CEVNT_TIMEOUT;
! 2507: peer->leap = LEAP_NOTINSYNC;
! 2508: up->watch = up->status = up->alarm = 0;
! 2509: up->avgint = MINAVG;
! 2510: up->freq = 0;
! 2511: up->gain = MAXGAIN / 2;
! 2512:
! 2513: /*
! 2514: * Initialize the station processes for audio gain, select bit,
! 2515: * station/frequency identifier and reference identifier. Start
! 2516: * probing at the strongest channel or the default channel if
! 2517: * nothing heard.
! 2518: */
! 2519: memset(up->mitig, 0, sizeof(up->mitig));
! 2520: for (i = 0; i < NCHAN; i++) {
! 2521: cp = &up->mitig[i];
! 2522: cp->gain = up->gain;
! 2523: cp->wwv.select = SELV;
! 2524: sprintf(cp->wwv.refid, "WV%.0f", floor(qsy[i]));
! 2525: cp->wwvh.select = SELH;
! 2526: sprintf(cp->wwvh.refid, "WH%.0f", floor(qsy[i]));
! 2527: }
! 2528: up->dchan = (DCHAN + NCHAN - 1) % NCHAN;
! 2529: wwv_newchan(peer);
! 2530: up->schan = up->dchan;
! 2531: }
! 2532:
! 2533: /*
! 2534: * wwv_metric - compute station metric
! 2535: *
! 2536: * The most significant bits represent the number of ones in the
! 2537: * station reachability register. The least significant bits represent
! 2538: * the minute sync pulse amplitude. The combined value is scaled 0-100.
! 2539: */
! 2540: double
! 2541: wwv_metric(
! 2542: struct sync *sp /* station pointer */
! 2543: )
! 2544: {
! 2545: double dtemp;
! 2546:
! 2547: dtemp = sp->count * MAXAMP;
! 2548: if (sp->synmax < MAXAMP)
! 2549: dtemp += sp->synmax;
! 2550: else
! 2551: dtemp += MAXAMP - 1;
! 2552: dtemp /= (AMAX + 1) * MAXAMP;
! 2553: return (dtemp * 100.);
! 2554: }
! 2555:
! 2556:
! 2557: #ifdef ICOM
! 2558: /*
! 2559: * wwv_qsy - Tune ICOM receiver
! 2560: *
! 2561: * This routine saves the AGC for the current channel, switches to a new
! 2562: * channel and restores the AGC for that channel. If a tunable receiver
! 2563: * is not available, just fake it.
! 2564: */
! 2565: static int
! 2566: wwv_qsy(
! 2567: struct peer *peer, /* peer structure pointer */
! 2568: int chan /* channel */
! 2569: )
! 2570: {
! 2571: int rval = 0;
! 2572: struct refclockproc *pp;
! 2573: struct wwvunit *up;
! 2574:
! 2575: pp = peer->procptr;
! 2576: up = (struct wwvunit *)pp->unitptr;
! 2577: if (up->fd_icom > 0) {
! 2578: up->mitig[up->achan].gain = up->gain;
! 2579: rval = icom_freq(up->fd_icom, peer->ttl & 0x7f,
! 2580: qsy[chan]);
! 2581: up->achan = chan;
! 2582: up->gain = up->mitig[up->achan].gain;
! 2583: }
! 2584: return (rval);
! 2585: }
! 2586: #endif /* ICOM */
! 2587:
! 2588:
! 2589: /*
! 2590: * timecode - assemble timecode string and length
! 2591: *
! 2592: * Prettytime format - similar to Spectracom
! 2593: *
! 2594: * sq yy ddd hh:mm:ss ld dut lset agc iden sig errs freq avgt
! 2595: *
! 2596: * s sync indicator ('?' or ' ')
! 2597: * q error bits (hex 0-F)
! 2598: * yyyy year of century
! 2599: * ddd day of year
! 2600: * hh hour of day
! 2601: * mm minute of hour
! 2602: * ss second of minute)
! 2603: * l leap second warning (' ' or 'L')
! 2604: * d DST state ('S', 'D', 'I', or 'O')
! 2605: * dut DUT sign and magnitude (0.1 s)
! 2606: * lset minutes since last clock update
! 2607: * agc audio gain (0-255)
! 2608: * iden reference identifier (station and frequency)
! 2609: * sig signal quality (0-100)
! 2610: * errs bit errors in last minute
! 2611: * freq frequency offset (PPM)
! 2612: * avgt averaging time (s)
! 2613: */
! 2614: static int
! 2615: timecode(
! 2616: struct wwvunit *up, /* driver structure pointer */
! 2617: char *ptr /* target string */
! 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: sprintf(ptr, "%c%1X", synchar, up->alarm);
! 2643: sprintf(cptr, " %4d %03d %02d:%02d:%02d %c%c %+d",
! 2644: year, day, hour, minute, second, leapchar, dst, dut);
! 2645: strcat(ptr, cptr);
! 2646:
! 2647: /*
! 2648: * Specific variable-format fields
! 2649: */
! 2650: sp = up->sptr;
! 2651: sprintf(cptr, " %d %d %s %.0f %d %.1f %d", up->watch,
! 2652: up->mitig[up->dchan].gain, sp->refid, sp->metric,
! 2653: up->errcnt, up->freq / SECOND * 1e6, up->avgint);
! 2654: strcat(ptr, cptr);
! 2655: return (strlen(ptr));
! 2656: }
! 2657:
! 2658:
! 2659: /*
! 2660: * wwv_gain - adjust codec gain
! 2661: *
! 2662: * This routine is called at the end of each second. During the second
! 2663: * the number of signal clips above the MAXAMP threshold (6000). If
! 2664: * there are no clips, the gain is bumped up; if there are more than
! 2665: * MAXCLP clips (100), it is bumped down. The decoder is relatively
! 2666: * insensitive to amplitude, so this crudity works just peachy. The
! 2667: * routine also jiggles the input port and selectively mutes the
! 2668: * monitor.
! 2669: */
! 2670: static void
! 2671: wwv_gain(
! 2672: struct peer *peer /* peer structure pointer */
! 2673: )
! 2674: {
! 2675: struct refclockproc *pp;
! 2676: struct wwvunit *up;
! 2677:
! 2678: pp = peer->procptr;
! 2679: up = (struct wwvunit *)pp->unitptr;
! 2680:
! 2681: /*
! 2682: * Apparently, the codec uses only the high order bits of the
! 2683: * gain control field. Thus, it may take awhile for changes to
! 2684: * wiggle the hardware bits.
! 2685: */
! 2686: if (up->clipcnt == 0) {
! 2687: up->gain += 4;
! 2688: if (up->gain > MAXGAIN)
! 2689: up->gain = MAXGAIN;
! 2690: } else if (up->clipcnt > MAXCLP) {
! 2691: up->gain -= 4;
! 2692: if (up->gain < 0)
! 2693: up->gain = 0;
! 2694: }
! 2695: audio_gain(up->gain, up->mongain, up->port);
! 2696: up->clipcnt = 0;
! 2697: #if DEBUG
! 2698: if (debug > 1)
! 2699: audio_show();
! 2700: #endif
! 2701: }
! 2702:
! 2703:
! 2704: #else
! 2705: int refclock_wwv_bs;
! 2706: #endif /* REFCLOCK */
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to refclock_wwv.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / ntp / ntpd