2 * refclock_wwv - clock driver for NIST WWV/H time/frequency station
8 #if defined(REFCLOCK) && defined(CLOCK_WWV)
12 #include "ntp_refclock.h"
13 #include "ntp_calendar.h"
14 #include "ntp_stdlib.h"
20 #ifdef HAVE_SYS_IOCTL_H
21 # include <sys/ioctl.h>
22 #endif /* HAVE_SYS_IOCTL_H */
24 #define ICOM 1 /* undefine to suppress ICOM code */
31 * Audio WWV/H demodulator/decoder
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
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.
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.
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.
64 * Interface definitions
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 */
89 #define MOD(x, y) ((x) < 0 ? -(-(x) % (y)) : (x) % (y))
92 * General purpose status bits (status)
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.
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.
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 */
119 * Station scoreboard bits (select)
121 * These are used to establish the signal quality for each of the five
122 * frequencies and two stations.
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 */
131 * Alarm status bits (alarm)
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.
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.
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 */
156 * Watchdog timeouts (watch)
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.
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 */
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.
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) */
189 * Tone frequency definitions.
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 */
197 * Acquisition and tracking time constants. Usually powers of 2.
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) */
205 * Miscellaneous status bits (misc)
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.
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 */
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.
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.
238 #define PDELAY (.0036 + .0011 + .0135) /* net system delay (s) */
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.
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 */
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.
275 int sw; /* case switch number */
276 int arg; /* argument */
280 * Case switch numbers
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 */
299 * Offsets in decoding matrix
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) */
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 */
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 */
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 */
371 * BCD coefficients for maximum likelihood digit decode
373 #define P15 1. /* max positive number */
374 #define N15 -1. /* max negative number */
379 #define P9 (P15 / 4) /* mark (+1) */
380 #define N9 (N15 / 4) /* space (-1) */
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 */
397 * Digits 0-6 (minute tens)
399 #define P6 (P15 / 3) /* mark (+1) */
400 #define N6 (N15 / 3) /* space (-1) */
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 */
414 * Digits 0-3 (day hundreds)
416 #define P3 (P15 / 2) /* mark (+1) */
417 #define N3 (N15 / 2) /* space (-1) */
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 */
428 * Digits 0-2 (hour tens)
430 #define P2 (P15 / 2) /* mark (+1) */
431 #define N2 (N15 / 2) /* space (-1) */
434 {N2, N2, 0, 0}, /* 0 */
435 {P2, N2, 0, 0}, /* 1 */
436 {N2, P2, 0, 0}, /* 2 */
437 {0, 0, 0, 0} /* backstop */
441 * DST decode (DST2 DST1) for prettyprint
444 'S', /* 00 standard time */
445 'I', /* 01 daylight warning */
446 'O', /* 10 standard warning */
447 'D' /* 11 daylight time */
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.
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 */
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.
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 */
494 * The channel structure is used to mitigate between channels. At this
495 * point we have already decided which station to use.
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 */
506 * WWV unit control structure
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 */
527 * Variables used to establish basic system timing
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 */
544 * Variables used to mitigate which channel to use
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 */
553 * Variables used by the clock state machine
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 */
564 * Variables used to estimate signal levels and bit/digit
567 double sigamp; /* I-channel peak signal amplitude */
568 double noiamp; /* I-channel average noise amplitude */
569 double datsnr; /* data SNR (dB) */
572 * Variables used to establish status and alarm conditions
574 int status; /* status bits */
575 int alarm; /* alarm flashers */
576 int misc; /* miscellaneous timecode bits */
577 int errcnt; /* data bit error counter */
581 * Function prototypes
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 *));
589 * More function prototypes
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 *,
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) */
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 */
624 * wwv_start - open the devices and initialize data for processing
628 int unit, /* instance number (not used) */
629 struct peer *peer /* peer structure pointer */
632 struct refclockproc *pp;
642 int fd; /* file descriptor */
644 double step; /* codec adjustment */
649 fd = audio_init(DEVICE_AUDIO);
658 * Allocate and initialize unit structure
660 if (!(up = (struct wwvunit *)
661 emalloc(sizeof(struct wwvunit)))) {
665 memset((char *)up, 0, sizeof(struct wwvunit));
667 pp->unitptr = (caddr_t)up;
668 pp->io.clock_recv = wwv_receive;
669 pp->io.srcclock = (caddr_t)peer;
672 if (!io_addclock(&pp->io)) {
679 * Initialize miscellaneous variables
681 peer->precision = PRECISION;
682 pp->clockdesc = DESCRIPTION;
683 memcpy((char *)&pp->refid, REFID, 4);
684 DTOLFP(1. / SECOND, &up->tick);
687 * The companded samples are encoded sign-magnitude. The table
688 * contains all the 256 values in the interest of speed.
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.;
694 for (i = 3; i < OFFSET; i++) {
695 up->comp[i] = up->comp[i - 1] + step;
696 up->comp[OFFSET + i] = -up->comp[i];
702 * Initialize the decoding matrix with the radix for each digit
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;
716 * Initialize the station processes for audio gain, select bit,
717 * station/frequency identifier and reference identifier.
720 for (i = 0; i < NCHAN; i++) {
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]));
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.
742 if (peer->ttlmax != 0) {
743 if (peer->ttlmax & 0x80)
744 up->fd_icom = icom_init("/dev/icom", B1200,
747 up->fd_icom = icom_init("/dev/icom", B9600,
750 if (up->fd_icom > 0) {
752 if ((temp = wwv_qsy(peer, up->schan)) < 0) {
753 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
755 "ICOM bus error; autotune disabled");
756 up->errflg = CEVNT_FAULT;
767 * wwv_shutdown - shut down the clock
771 int unit, /* instance number (not used) */
772 struct peer *peer /* peer structure pointer */
775 struct refclockproc *pp;
779 up = (struct wwvunit *)pp->unitptr;
780 io_closeclock(&pp->io);
788 * wwv_receive - receive data from the audio device
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.
797 struct recvbuf *rbufp /* receive buffer structure pointer */
801 struct refclockproc *pp;
807 double sample; /* codec sample */
808 u_char *dpt; /* buffer pointer */
814 peer = (struct peer *)rbufp->recv_srcclock;
816 up = (struct wwvunit *)pp->unitptr;
819 * Main loop - read until there ain't no more. Note codec
820 * samples are bit-inverted.
822 up->timestamp = rbufp->recv_time;
823 up->bufcnt = rbufp->recv_length;
824 DTOLFP((double)up->bufcnt / SECOND, <emp);
825 L_SUB(&up->timestamp, <emp);
826 dpt = rbufp->recv_buffer;
827 for (up->bufptr = 0; up->bufptr < up->bufcnt; up->bufptr++) {
828 sample = up->comp[~*dpt & 0xff];
831 * Clip noise spikes greater than MAXSIG. If no clips,
832 * increase the gain a tad; if the clips are too high,
835 if (sample > MAXSIG) {
838 } else if (sample < -MAXSIG) {
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.
848 up->phase += up->freq / SECOND;
849 if (up->phase >= .5) {
851 } else if (up->phase < - .5) {
853 wwv_rf(peer, sample);
854 wwv_rf(peer, sample);
856 wwv_rf(peer, sample);
858 L_ADD(&up->timestamp, &up->tick);
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.
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)
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
892 if (dtemp > up->comp[i])
894 else if (dtemp < up->comp[i])
901 *dpt = ~(i + OFFSET);
908 * Squawk to the monitor speaker if enabled.
910 if (pp->sloppyclockflag & CLK_FLAG3)
911 if (write(pp->io.fd, (u_char *)&rbufp->recv_space,
912 (u_int)up->bufcnt) < 0)
918 * wwv_poll - called by the transmit procedure
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.
927 int unit, /* instance number (not used) */
928 struct peer *peer /* peer structure pointer */
931 struct refclockproc *pp;
935 up = (struct wwvunit *)pp->unitptr;
936 if (pp->coderecv == pp->codeproc)
937 up->errflg = CEVNT_TIMEOUT;
940 if (up->status & INSYNC)
943 refclock_report(peer, up->errflg);
949 * wwv_rf - process signals and demodulate to baseband
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.
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
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.
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.
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).
989 struct peer *peer, /* peerstructure pointer */
990 double isig /* input signal */
993 struct refclockproc *pp;
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 */
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 */
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 */
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 */
1023 static int iniflg; /* initialization flag */
1030 up = (struct wwvunit *)pp->unitptr;
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));
1044 up->monitor = isig; /* change for debug */
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
1053 * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB
1054 * passband ripple, -50 dB stopband ripple.
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;
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.
1076 i = up->datapt - up->pdelay % 80;
1079 up->datapt = (up->datapt + IN100) % 80;
1080 dtemp = sintab[i] * data / DATSIZ * DGAIN;
1081 up->irig -= ibuf[iptr];
1085 dtemp = sintab[i] * data / DATSIZ * DGAIN;
1086 up->qrig -= qbuf[iptr];
1089 iptr = (iptr + 1) % DATSIZ;
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.
1097 * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB
1098 * passband ripple, -50 dB stopband ripple.
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;
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.
1127 up->mphase = (up->mphase + 1) % MINUTE;
1130 csinptr = (csinptr + IN1000) % 80;
1131 dtemp = sintab[i] * syncx / SYNSIZ * SGAIN;
1132 ciamp = ciamp - cibuf[jptr] + dtemp;
1133 cibuf[jptr] = dtemp;
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);
1142 hsinptr = (hsinptr + IN1200) % 80;
1143 dtemp = sintab[i] * syncx / SYNSIZ * SGAIN;
1144 hiamp = hiamp - hibuf[jptr] + dtemp;
1145 hibuf[jptr] = dtemp;
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);
1153 jptr = (jptr + 1) % SYNSIZ;
1155 if (up->mphase == 0) {
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.
1167 if (up->rsec == 60) {
1168 up->mphase -= SECOND;
1170 up->mphase += MINUTE;
1171 } else if (!(up->status & MSYNC)) {
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.
1183 if (up->mitig[up->achan].wwv.count >=
1184 up->mitig[up->achan].wwvh.count)
1185 sp = &up->mitig[up->achan].wwv;
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;
1194 up->rsec = (MINUTE - ltemp) / SECOND;
1195 if (!(up->status & SSYNC)) {
1196 up->repoch = ltemp % SECOND;
1197 up->yepoch = up->repoch -
1200 up->yepoch += SECOND;
1203 } else if (sp->count == 0 || up->watch >= ACQSN)
1205 up->watch = sp->count = 0;
1206 up->schan = (up->schan + 1) % NCHAN;
1207 wwv_qsy(peer, up->schan);
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
1219 if (up->watch > DIGIT && !(up->status & DSYNC))
1220 up->watch = up->status = 0;
1223 * If the second sync times out, dim the sync
1224 * lamp and raise an alarm.
1227 if (up->swatch > HSPEC)
1228 up->status &= ~SSYNC;
1229 if (!(up->status & SSYNC))
1230 up->alarm |= 1 << SYNERR;
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).
1240 if (up->status & SELV) {
1241 up->pdelay = up->cdelay;
1244 * WWV FIR matched filter, five cycles of 1000-Hz
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
1288 } else if (up->status & SELH) {
1289 up->pdelay = up->hdelay;
1292 * WWVH FIR matched filter, six cycles of 1200-Hz
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;
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;
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;
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;
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
1346 up->epoch = (up->epoch + 1) % SECOND;
1347 if (up->epoch == 0) {
1348 wwv_endpoc(peer, epomax, epopos);
1349 up->epomax = epomax;
1351 if (!(up->status & MSYNC))
1354 dtemp = (epobuf[up->epoch] += (mfsync - epobuf[up->epoch]) /
1355 (MINAVG << up->avgint));
1356 if (dtemp > epomax) {
1358 epopos = up->epoch - up->pdelay - 5 * MS;
1362 if (up->status & MSYNC)
1368 * wwv_qrz - identify and acquire WWV/WWVH minute sync pulse
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
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.
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.
1395 struct peer *peer, /* peerstructure pointer */
1396 struct sync *sp, /* sync channel structure */
1397 double syncx /* bandpass filtered sync signal */
1400 struct refclockproc *pp;
1402 char tbuf[80]; /* monitor buffer */
1403 double snr; /* on-pulse/off-pulse ratio (dB) */
1408 up = (struct wwvunit *)pp->unitptr;
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.
1418 if (!(up->status & MSYNC) || up->status & SFLAG) {
1420 if (up->status & MSYNC)
1422 else if (sp->count > 1)
1425 epoch = sp->lastpos;
1426 if (syncx > sp->sigmax) {
1428 sp->pos = up->mphase;
1430 if (abs(MOD(up->mphase - epoch, MINUTE)) > SYNSIZ &&
1431 syncx > sp->noise) {
1435 if (up->mphase == 0) {
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).
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) {
1451 * If in minute sync, just count the runs up and
1454 if (sp->select & (DATANG | SYNCNG | JITRNG)) {
1458 if (sp->count < AMAX)
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.
1468 snr = wwv_snr(sp->sigmax, sp->noise);
1469 isgood = sp->sigmax > ATHR && snr > ASNR &&
1470 !(sp->select & JITRNG);
1471 switch (sp->count) {
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.
1481 if (sp->sigmax >= ATHR)
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
1497 if (sp->sigmax < ATHR) {
1500 } else if (!isgood) {
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.
1517 sp->mepoch = sp->pos;
1518 if (sp->count < AMAX)
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);
1531 printf("%s\n", tbuf);
1534 sp->lastmax = sp->sigmax;
1535 sp->lastpos = sp->pos;
1536 sp->sigmax = sp->noise = 0;
1542 * wwv_endpoc - process receiver epoch
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).
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.
1555 struct peer *peer, /* peer structure pointer */
1556 double epomax, /* epoch max */
1557 int epopos /* epoch max position */
1560 struct refclockproc *pp;
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 */
1573 static int iniflg; /* initialization flag */
1574 char tbuf[80]; /* monitor buffer */
1579 up = (struct wwvunit *)pp->unitptr;
1582 memset((char *)epoch_mf, 0, sizeof(epoch_mf));
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.
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];
1602 tepoch = epoch_mf[2]; /* 0 2 1 */
1603 tspan = epoch_mf[0] - epoch_mf[1];
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];
1613 tepoch = epoch_mf[2]; /* 1 2 0 */
1614 tspan = epoch_mf[1] - epoch_mf[0];
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.
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.
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))
1644 if (syncnt < SCMP) {
1647 up->status |= SSYNC;
1649 up->repoch = tepoch;
1650 up->yepoch = up->repoch;
1652 up->yepoch += SECOND;
1658 syncnt = avgcnt = 0;
1660 if (!(up->status & SSYNC) && 0) {
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);
1669 printf("%s\n", tbuf);
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
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.
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)
1694 else if (up->freq < -MAXFREQ)
1695 up->freq = -MAXFREQ;
1696 if (abs(tmp3) <= 1 && up->avgint < MAXAVG) {
1704 if (up->avgint < MAXAVG) {
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)
1713 &peer->srcadr, tbuf);
1716 printf("%s\n", tbuf);
1727 * wwv_epoch - main loop
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
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).
1750 struct peer *peer /* peer structure pointer */
1753 static double dpulse; /* data pulse length */
1754 struct refclockproc *pp;
1758 l_fp offset; /* NTP format offset */
1762 up = (struct wwvunit *)pp->unitptr;
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.
1769 cp = &up->mitig[up->achan];
1770 if (up->rphase == 800 * MS) {
1772 sp->synamp = SQRT(sp->amp);
1774 sp->synamp = SQRT(sp->amp);
1777 if (up->rsec == 0) {
1778 up->sigamp = up->datsnr = 0;
1782 * Estimate the noise level by integrating the I-channel
1783 * energy at epoch 30 ms.
1785 if (up->rphase == 30 * MS) {
1786 if (!(up->status & SFLAG))
1787 up->noiamp += (up->irig - up->noiamp) /
1788 (MINAVG << up->avgint);
1790 cp->noiamp += (SQRT(up->irig *
1791 up->irig + up->qrig * up->qrig) -
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
1801 } else if (up->rphase == 200 * MS) {
1802 if (!(up->status & SFLAG)) {
1803 up->sigamp = up->irig;
1806 up->datsnr = wwv_snr(up->sigamp,
1808 up->datpha = up->qrig / (MINAVG <<
1810 if (up->datpha >= 0) {
1812 if (up->datapt >= 80)
1820 up->sigamp = SQRT(up->irig * up->irig +
1821 up->qrig * up->qrig);
1822 up->datsnr = wwv_snr(up->sigamp,
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
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;
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.
1849 if (up->epoch == up->yepoch) {
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
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 *
1869 pp->year = up->decvec[YR].digit +
1870 up->decvec[YR + 1].digit * 10;
1871 if (pp->year < UTCYEAR)
1877 * We have to simulate refclock_process() here,
1878 * since the fudgetime gets added much earlier
1881 pp->lastrec = up->timestamp;
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;
1888 refclock_process_offset(pp, offset,
1894 * At the end of the receiver second, process the data bit and
1895 * update the decoding matrix probabilities.
1898 if (up->epoch == up->repoch) {
1899 wwv_rsec(peer, dpulse);
1901 up->rphase = dpulse = 0;
1907 * wwv_rsec - process receiver second
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.
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.
1924 struct peer *peer, /* peer structure pointer */
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;
1934 struct sync *sp, *rp;
1935 double bit; /* bit likelihood */
1936 char tbuf[80]; /* monitor buffer */
1940 up = (struct wwvunit *)pp->unitptr;
1943 memset((char *)bitvec, 0, sizeof(bitvec));
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.
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;
1963 * Ignore this second.
1965 case IDLE: /* 9, 45-49 */
1969 * Probe channel stuff
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
1980 cp = &up->mitig[up->achan];
1982 sp->synmax = sp->synamp;
1984 sp->synmax = sp->synamp;
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.
1995 cp = &up->mitig[up->achan];
1996 if (up->sigamp < DTHR || up->datsnr < DSNR)
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;
2006 sp->select |= DATANG;
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;
2015 rp->select |= DATANG;
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);
2028 printf("%s\n", tbuf);
2030 up->status &= ~SFLAG;
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.
2041 case COEF1: /* 10-13 */
2042 if (up->status & DGATE)
2043 up->status |= BGATE;
2047 case COEF2: /* 18, 27-28, 42-43 */
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;
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.
2064 case DECIM2: /* 29 */
2065 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd2);
2068 case DECIM3: /* 44 */
2069 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd3);
2072 case DECIM6: /* 19 */
2073 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd6);
2076 case DECIM9: /* 8, 14, 24, 34, 39 */
2077 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd9);
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
2086 case MSC20: /* 55 */
2087 wwv_corr4(peer, &up->decvec[YR + 1], bcddld, bcd9);
2090 case MSCBIT: /* 2, 3, 50, 56-57 */
2091 if (bitvec[up->rsec] > BTHR)
2093 else if (bitvec[up->rsec] < -BTHR)
2096 up->alarm |= 1 << SYMERR;
2099 case MSC21: /* 58 */
2100 if (bitvec[up->rsec] > BTHR)
2102 else if (bitvec[up->rsec] < -BTHR)
2105 up->alarm |= 1 << SYMERR;
2106 up->schan = (up->schan + 1) % NCHAN;
2107 wwv_qsy(peer, up->schan);
2108 up->status |= SFLAG;
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.
2119 cp = &up->mitig[up->achan];
2120 if (up->sigamp < DTHR || up->datsnr < DSNR)
2123 sp->synmin = sp->synamp;
2125 sp->synmin = sp->synamp;
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
2136 if (up->tsec == 60) {
2137 up->status &= ~LEPSEC;
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.
2149 case MIN2: /* 59/60 */
2150 up->minset = ((current_time - peer->update) + 30) / 60;
2152 up->status |= DSYNC;
2153 if (up->digcnt >= 9 && (up->alarm & (3 << SYNERR)) == 0)
2155 up->status |= INSYNC;
2158 pp->lencode = timecode(up, pp->a_lastcode);
2159 if (up->misc & SECWAR)
2160 pp->leap = LEAP_ADDSECOND;
2162 pp->leap = LEAP_NOWARNING;
2163 refclock_receive(peer);
2164 record_clock_stats(&peer->srcadr, pp->a_lastcode);
2167 printf("wwv: timecode %d %s\n", pp->lencode,
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.
2179 if (up->minset > PANIC)
2181 up->alarm = (up->alarm & ~0x8888) << 1;
2182 up->nepoch = (up->mphase + SYNSIZ) % MINUTE;
2183 up->errcnt = up->digcnt = nsec = 0;
2186 if (!(up->status & DSYNC)) {
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);
2195 printf("%s\n", tbuf);
2198 up->rsec = up->tsec = nsec;
2204 * wwv_data - calculate bit probability
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..
2216 struct wwvunit *up, /* driver unit pointer */
2217 double pulse /* pulse length (sample units) */
2220 double p0, p1, p2; /* probabilities */
2221 double dpulse; /* pulse length in ms */
2224 dpulse = pulse - DATSIZ / 2;
2227 * If the data amplitude or SNR are below threshold or if the
2228 * pulse length is less than 170 ms, declare an erasure.
2230 if (up->sigamp < DTHR || up->datsnr < DSNR || dpulse < DATSIZ) {
2231 up->status |= DGATE;
2233 if (up->errcnt > MAXERR)
2234 up->alarm |= 1 << MODERR;
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.
2244 up->status &= ~DGATE;
2245 if (dpulse < (200 * MS)) {
2247 } else if (dpulse < 500 * MS) {
2249 p1 = dpulse / (300 * MS);
2251 } else if (dpulse < 800 * MS) {
2253 p2 = dpulse / (300 * MS);
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.
2265 return ((p1 - p0) * MAXSIG);
2270 * wwv_corr4 - determine maximum likelihood digit
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.
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 */
2287 struct refclockproc *pp;
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 */
2298 up = (struct wwvunit *)pp->unitptr;
2301 * Correlate digit vector with each BCD coefficient vector. If
2302 * any BCD digit bit is bad, consider all bits a miss.
2305 topmax = nxtmax = -MAXSIG;
2306 for (i = 0; tab[i][0] != 0; i++) {
2308 for (j = 0; j < 4; j++) {
2309 if (!(up->status & BGATE))
2310 acc += data[j] * tab[i][j];
2312 acc = (vp->like[i] += (acc - vp->like[i]) / TCONST);
2317 } else if (acc > nxtmax) {
2321 vp->mldigit = mldigit;
2322 vp->digprb = topmax;
2323 vp->digsnr = wwv_snr(topmax, nxtmax);
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.
2332 diff = mldigit - vp->digit;
2335 if (diff != vp->phase) {
2339 if (vp->digprb < BTHR || vp->digsnr < BSNR) {
2341 up->alarm |= 1 << SYMERR;
2342 } else if (vp->count < BCMP) {
2343 if (!(up->status & INSYNC)) {
2345 vp->digit = mldigit;
2350 vp->digit = mldigit;
2353 if (vp->digit != mldigit)
2354 up->alarm |= 1 << DECERR;
2355 if (!(up->status & INSYNC)) {
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);
2365 printf("%s\n", tbuf);
2368 up->status &= ~BGATE;
2373 * wwv_tsec - transmitter second processing
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.
2383 struct wwvunit *up /* driver structure pointer */
2386 int minute, day, isleap;
2390 if (up->tsec < 60 || up->status & LEPSEC)
2395 * Advance minute unit of the day. If the minute unit is not
2396 * synchronized, go no further.
2398 temp = carry(&up->decvec[MN]); /* minute units */
2399 if (!(up->status & DSYNC))
2403 * Propagate carries through the day.
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]);
2410 temp = carry(&up->decvec[HR + 1]);
2413 * Decode the current minute and day. Set the leap second enable
2414 * bit on the last minute of 30 June and 31 December.
2416 minute = up->decvec[MN].digit + up->decvec[MN + 1].digit *
2417 10 + up->decvec[HR].digit * 60 + up->decvec[HR +
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;
2427 * Roll the day if this the first minute and propagate carries
2433 while (carry(&up->decvec[HR]) != 0); /* advance to minute 0 */
2434 while (carry(&up->decvec[HR + 1]) != 0);
2436 temp = carry(&up->decvec[DA]); /* carry days */
2438 temp = carry(&up->decvec[DA + 1]);
2440 temp = carry(&up->decvec[DA + 2]);
2443 * Roll the year if this the first day and propagate carries
2444 * through the century.
2446 if (day != (isleap ? 365 : 366))
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 */
2454 carry(&up->decvec[YR + 1]);
2459 * carry - process digit
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.
2468 struct decvec *dp /* decoding table pointer */
2474 dp->digit++; /* advance clock digit */
2475 if (dp->digit == dp->radix) { /* modulo radix */
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];
2487 * wwv_snr - compute SNR or likelihood function
2491 double signal, /* signal */
2492 double noise /* noise */
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.
2510 } else if (noise <= 0) {
2513 rval = 20 * log10(signal / noise);
2521 * wwv_newchan - change to new data channel
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.
2536 struct peer *peer /* peer structure pointer */
2539 struct refclockproc *pp;
2542 struct sync *sp, *rp;
2547 up = (struct wwvunit *)pp->unitptr;
2550 * Reset the matched filter selector and station pointer to
2551 * avoid fooling around should we lose this game.
2554 up->status &= ~(SELV | SELH);
2557 * Search all five station pairs looking for the station with
2558 * the maximum compare counter. Ties go to the highest frequency
2562 sp = (struct sync *)0;
2564 for (i = 0; i < NCHAN; i++) {
2567 if (rp->count >= rank) {
2573 if (rp->count >= rank) {
2581 * If we find a station, continue to track it. If not, X marks
2582 * the spot and we wait for better ions.
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);
2597 * wwv_qsy - Tune ICOM receiver
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.
2605 struct peer *peer, /* peer structure pointer */
2606 int chan /* channel */
2609 struct refclockproc *pp;
2614 up = (struct wwvunit *)pp->unitptr;
2615 up->mitig[up->achan].gain = up->gain;
2617 if (up->fd_icom > 0)
2618 rval = icom_freq(up->fd_icom, peer->ttlmax & 0x7f,
2622 up->gain = up->mitig[up->achan].gain;
2628 * timecode - assemble timecode string and length
2630 * Prettytime format - similar to Spectracom
2632 * sq yy ddd hh:mm:ss.fff ld dut lset agc stn comp errs freq avgt
2634 * s sync indicator ('?' or ' ')
2635 * q quality character (hex 0-F)
2636 * yyyy year of century
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)
2655 struct wwvunit *up, /* driver structure pointer */
2656 char *ptr /* target string */
2660 int year, day, hour, minute, second, frac, dut;
2661 char synchar, qual, leapchar, dst;
2666 * Common fixed-format fields
2668 synchar = (up->status & INSYNC) ? ' ' : '?';
2670 if (up->alarm & (3 << DECERR))
2672 if (up->alarm & (3 << SYMERR))
2674 if (up->alarm & (3 << MODERR))
2676 if (up->alarm & (3 << SYNERR))
2678 year = up->decvec[7].digit + up->decvec[7].digit * 10;
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;
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))
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);
2700 * Specific variable-format fields
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 <<
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);
2713 return (strlen(ptr));
2718 * wwv_gain - adjust codec gain
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.
2728 struct peer *peer /* peer structure pointer */
2731 struct refclockproc *pp;
2735 up = (struct wwvunit *)pp->unitptr;
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.
2742 if (up->clipcnt == 0) {
2746 } else if (up->clipcnt > SECOND / 100) {
2751 audio_gain(up->gain, up->port);
2757 int refclock_wwv_bs;
2758 #endif /* REFCLOCK */