contrib/ntp/ntpd/refclock_wwv.c

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