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 */
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. Transmissions are made continuously on
36 * 2.5, 5, 10 and 15 MHz from WWV and WWVH, and 20 MHz from WWV. An
37 * ordinary AM shortwave receiver can be tuned manually to one of these
38 * frequencies or, in the case of ICOM receivers, the receiver can be
39 * tuned automatically using this program as propagation conditions
40 * change throughout the weasons, both day and night.
42 * The driver receives, demodulates and decodes the radio signals when
43 * connected to the audio codec of a workstation running Solaris, SunOS
44 * FreeBSD or Linux, and with a little help, other workstations with
45 * similar codecs or sound cards. In this implementation, only one audio
46 * driver 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.html. 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.
65 * Fudge flag4 causes the dubugging output described above to be
66 * recorded in the clockstats file. Fudge flag2 selects the audio input
67 * port, where 0 is the mike port (default) and 1 is the line-in port.
68 * It does not seem useful to select the compact disc player port. Fudge
69 * flag3 enables audio monitoring of the input signal. For this purpose,
70 * the monitor gain is set to a default value.
73 * General definitions. These ordinarily do not need to be changed.
75 #define DEVICE_AUDIO "/dev/audio" /* audio device name */
76 #define AUDIO_BUFSIZ 320 /* audio buffer size (50 ms) */
77 #define PRECISION (-10) /* precision assumed (about 1 ms) */
78 #define DESCRIPTION "WWV/H Audio Demodulator/Decoder" /* WRU */
79 #define SECOND 8000 /* second epoch (sample rate) (Hz) */
80 #define MINUTE (SECOND * 60) /* minute epoch */
81 #define OFFSET 128 /* companded sample offset */
82 #define SIZE 256 /* decompanding table size */
83 #define MAXAMP 6000. /* max signal level reference */
84 #define MAXCLP 100 /* max clips above reference per s */
85 #define MAXSNR 40. /* max SNR reference */
86 #define MAXFREQ 1.5 /* max frequency tolerance (187 PPM) */
87 #define DATCYC 170 /* data filter cycles */
88 #define DATSIZ (DATCYC * MS) /* data filter size */
89 #define SYNCYC 800 /* minute filter cycles */
90 #define SYNSIZ (SYNCYC * MS) /* minute filter size */
91 #define TCKCYC 5 /* tick filter cycles */
92 #define TCKSIZ (TCKCYC * MS) /* tick filter size */
93 #define NCHAN 5 /* number of radio channels */
94 #define AUDIO_PHI 5e-6 /* dispersion growth factor */
97 * Tunable parameters. The DGAIN parameter can be changed to fit the
98 * audio response of the radio at 100 Hz. The WWV/WWVH data subcarrier
99 * is transmitted at about 20 percent percent modulation; the matched
100 * filter boosts it by a factor of 17 and the receiver response does
101 * what it does. The compromise value works for ICOM radios. If the
102 * radio is not tunable, the DCHAN parameter can be changed to fit the
103 * expected best propagation frequency: higher if further from the
104 * transmitter, lower if nearer. The compromise value works for the US
105 * right coast. The FREQ_OFFSET parameter can be used as a frequency
106 * vernier to correct codec requency if greater than MAXFREQ.
108 #define DCHAN 3 /* default radio channel (15 Mhz) */
109 #define DGAIN 5. /* subcarrier gain */
110 #define FREQ_OFFSET 0. /* codec frequency correction (PPM) */
113 * General purpose status bits (status)
115 * SELV and/or SELH are set when WWV or WWVH have been heard and cleared
116 * on signal loss. SSYNC is set when the second sync pulse has been
117 * acquired and cleared by signal loss. MSYNC is set when the minute
118 * sync pulse has been acquired. DSYNC is set when the units digit has
119 * has reached the threshold and INSYNC is set when all nine digits have
120 * reached the threshold. The MSYNC, DSYNC and INSYNC bits are cleared
121 * only by timeout, upon which the driver starts over from scratch.
123 * DGATE is lit if the data bit amplitude or SNR is below thresholds and
124 * BGATE is lit if the pulse width amplitude or SNR is below thresolds.
125 * LEPSEC is set during the last minute of the leap day. At the end of
126 * this minute the driver inserts second 60 in the seconds state machine
127 * and the minute sync slips a second.
129 #define MSYNC 0x0001 /* minute epoch sync */
130 #define SSYNC 0x0002 /* second epoch sync */
131 #define DSYNC 0x0004 /* minute units sync */
132 #define INSYNC 0x0008 /* clock synchronized */
133 #define FGATE 0x0010 /* frequency gate */
134 #define DGATE 0x0020 /* data pulse amplitude error */
135 #define BGATE 0x0040 /* data pulse width error */
136 #define LEPSEC 0x1000 /* leap minute */
139 * Station scoreboard bits
141 * These are used to establish the signal quality for each of the five
142 * frequencies and two stations.
144 #define SELV 0x0100 /* WWV station select */
145 #define SELH 0x0200 /* WWVH station select */
148 * Alarm status bits (alarm)
150 * These bits indicate various alarm conditions, which are decoded to
151 * form the quality character included in the timecode.
153 #define CMPERR 1 /* digit or misc bit compare error */
154 #define LOWERR 2 /* low bit or digit amplitude or SNR */
155 #define NINERR 4 /* less than nine digits in minute */
156 #define SYNERR 8 /* not tracking second sync */
159 * Watchcat timeouts (watch)
161 * If these timeouts expire, the status bits are mashed to zero and the
162 * driver starts from scratch. Suitably more refined procedures may be
163 * developed in future. All these are in minutes.
165 #define ACQSN 6 /* station acquisition timeout */
166 #define DATA 15 /* unit minutes timeout */
167 #define SYNCH 40 /* station sync timeout */
168 #define PANIC (2 * 1440) /* panic timeout */
171 * Thresholds. These establish the minimum signal level, minimum SNR and
172 * maximum jitter thresholds which establish the error and false alarm
173 * rates of the driver. The values defined here may be on the
174 * adventurous side in the interest of the highest sensitivity.
176 #define MTHR 13. /* minute sync gate (percent) */
177 #define TTHR 50. /* minute sync threshold (percent) */
178 #define AWND 20 /* minute sync jitter threshold (ms) */
179 #define ATHR 2500. /* QRZ minute sync threshold */
180 #define ASNR 20. /* QRZ minute sync SNR threshold (dB) */
181 #define QTHR 2500. /* QSY minute sync threshold */
182 #define QSNR 20. /* QSY minute sync SNR threshold (dB) */
183 #define STHR 2500. /* second sync threshold */
184 #define SSNR 15. /* second sync SNR threshold (dB) */
185 #define SCMP 10 /* second sync compare threshold */
186 #define DTHR 1000. /* bit threshold */
187 #define DSNR 10. /* bit SNR threshold (dB) */
188 #define AMIN 3 /* min bit count */
189 #define AMAX 6 /* max bit count */
190 #define BTHR 1000. /* digit threshold */
191 #define BSNR 3. /* digit likelihood threshold (dB) */
192 #define BCMP 3 /* digit compare threshold */
193 #define MAXERR 40 /* maximum error alarm */
196 * Tone frequency definitions. The increments are for 4.5-deg sine
199 #define MS (SECOND / 1000) /* samples per millisecond */
200 #define IN100 ((100 * 80) / SECOND) /* 100 Hz increment */
201 #define IN1000 ((1000 * 80) / SECOND) /* 1000 Hz increment */
202 #define IN1200 ((1200 * 80) / SECOND) /* 1200 Hz increment */
205 * Acquisition and tracking time constants
207 #define MINAVG 8 /* min averaging time */
208 #define MAXAVG 1024 /* max averaging time */
209 #define FCONST 3 /* frequency time constant */
210 #define TCONST 16 /* data bit/digit time constant */
213 * Miscellaneous status bits (misc)
215 * These bits correspond to designated bits in the WWV/H timecode. The
216 * bit probabilities are exponentially averaged over several minutes and
217 * processed by a integrator and threshold.
219 #define DUT1 0x01 /* 56 DUT .1 */
220 #define DUT2 0x02 /* 57 DUT .2 */
221 #define DUT4 0x04 /* 58 DUT .4 */
222 #define DUTS 0x08 /* 50 DUT sign */
223 #define DST1 0x10 /* 55 DST1 leap warning */
224 #define DST2 0x20 /* 2 DST2 DST1 delayed one day */
225 #define SECWAR 0x40 /* 3 leap second warning */
228 * The on-time synchronization point for the driver is the second epoch
229 * sync pulse produced by the FIR matched filters. As the 5-ms delay of
230 * these filters is compensated, the program delay is 1.1 ms due to the
231 * 600-Hz IIR bandpass filter. The measured receiver delay is 4.7 ms and
232 * the codec delay less than 0.2 ms. The additional propagation delay
233 * specific to each receiver location can be programmed in the fudge
234 * time1 and time2 values for WWV and WWVH, respectively.
236 #define PDELAY (.0011 + .0047 + .0002) /* net system delay (s) */
239 * Table of sine values at 4.5-degree increments. This is used by the
240 * synchronous matched filter demodulators.
243 0.000000e+00, 7.845910e-02, 1.564345e-01, 2.334454e-01, /* 0-3 */
244 3.090170e-01, 3.826834e-01, 4.539905e-01, 5.224986e-01, /* 4-7 */
245 5.877853e-01, 6.494480e-01, 7.071068e-01, 7.604060e-01, /* 8-11 */
246 8.090170e-01, 8.526402e-01, 8.910065e-01, 9.238795e-01, /* 12-15 */
247 9.510565e-01, 9.723699e-01, 9.876883e-01, 9.969173e-01, /* 16-19 */
248 1.000000e+00, 9.969173e-01, 9.876883e-01, 9.723699e-01, /* 20-23 */
249 9.510565e-01, 9.238795e-01, 8.910065e-01, 8.526402e-01, /* 24-27 */
250 8.090170e-01, 7.604060e-01, 7.071068e-01, 6.494480e-01, /* 28-31 */
251 5.877853e-01, 5.224986e-01, 4.539905e-01, 3.826834e-01, /* 32-35 */
252 3.090170e-01, 2.334454e-01, 1.564345e-01, 7.845910e-02, /* 36-39 */
253 -0.000000e+00, -7.845910e-02, -1.564345e-01, -2.334454e-01, /* 40-43 */
254 -3.090170e-01, -3.826834e-01, -4.539905e-01, -5.224986e-01, /* 44-47 */
255 -5.877853e-01, -6.494480e-01, -7.071068e-01, -7.604060e-01, /* 48-51 */
256 -8.090170e-01, -8.526402e-01, -8.910065e-01, -9.238795e-01, /* 52-55 */
257 -9.510565e-01, -9.723699e-01, -9.876883e-01, -9.969173e-01, /* 56-59 */
258 -1.000000e+00, -9.969173e-01, -9.876883e-01, -9.723699e-01, /* 60-63 */
259 -9.510565e-01, -9.238795e-01, -8.910065e-01, -8.526402e-01, /* 64-67 */
260 -8.090170e-01, -7.604060e-01, -7.071068e-01, -6.494480e-01, /* 68-71 */
261 -5.877853e-01, -5.224986e-01, -4.539905e-01, -3.826834e-01, /* 72-75 */
262 -3.090170e-01, -2.334454e-01, -1.564345e-01, -7.845910e-02, /* 76-79 */
263 0.000000e+00}; /* 80 */
266 * Decoder operations at the end of each second are driven by a state
267 * machine. The transition matrix consists of a dispatch table indexed
268 * by second number. Each entry in the table contains a case switch
269 * number and argument.
272 int sw; /* case switch number */
273 int arg; /* argument */
277 * Case switch numbers
279 #define IDLE 0 /* no operation */
280 #define COEF 1 /* BCD bit */
281 #define COEF1 2 /* BCD bit for minute unit */
282 #define COEF2 3 /* BCD bit not used */
283 #define DECIM9 4 /* BCD digit 0-9 */
284 #define DECIM6 5 /* BCD digit 0-6 */
285 #define DECIM3 6 /* BCD digit 0-3 */
286 #define DECIM2 7 /* BCD digit 0-2 */
287 #define MSCBIT 8 /* miscellaneous bit */
288 #define MSC20 9 /* miscellaneous bit */
289 #define MSC21 10 /* QSY probe channel */
290 #define MIN1 11 /* latch time */
291 #define MIN2 12 /* leap second */
292 #define SYNC2 13 /* latch minute sync pulse */
293 #define SYNC3 14 /* latch data pulse */
296 * Offsets in decoding matrix
298 #define MN 0 /* minute digits (2) */
299 #define HR 2 /* hour digits (2) */
300 #define DA 4 /* day digits (3) */
301 #define YR 7 /* year digits (2) */
303 struct progx progx[] = {
304 {SYNC2, 0}, /* 0 latch minute sync pulse */
305 {SYNC3, 0}, /* 1 latch data pulse */
306 {MSCBIT, DST2}, /* 2 dst2 */
307 {MSCBIT, SECWAR}, /* 3 lw */
308 {COEF, 0}, /* 4 1 year units */
312 {DECIM9, YR}, /* 8 */
313 {IDLE, 0}, /* 9 p1 */
314 {COEF1, 0}, /* 10 1 minute units */
315 {COEF1, 1}, /* 11 2 */
316 {COEF1, 2}, /* 12 4 */
317 {COEF1, 3}, /* 13 8 */
318 {DECIM9, MN}, /* 14 */
319 {COEF, 0}, /* 15 10 minute tens */
320 {COEF, 1}, /* 16 20 */
321 {COEF, 2}, /* 17 40 */
322 {COEF2, 3}, /* 18 80 (not used) */
323 {DECIM6, MN + 1}, /* 19 p2 */
324 {COEF, 0}, /* 20 1 hour units */
325 {COEF, 1}, /* 21 2 */
326 {COEF, 2}, /* 22 4 */
327 {COEF, 3}, /* 23 8 */
328 {DECIM9, HR}, /* 24 */
329 {COEF, 0}, /* 25 10 hour tens */
330 {COEF, 1}, /* 26 20 */
331 {COEF2, 2}, /* 27 40 (not used) */
332 {COEF2, 3}, /* 28 80 (not used) */
333 {DECIM2, HR + 1}, /* 29 p3 */
334 {COEF, 0}, /* 30 1 day units */
335 {COEF, 1}, /* 31 2 */
336 {COEF, 2}, /* 32 4 */
337 {COEF, 3}, /* 33 8 */
338 {DECIM9, DA}, /* 34 */
339 {COEF, 0}, /* 35 10 day tens */
340 {COEF, 1}, /* 36 20 */
341 {COEF, 2}, /* 37 40 */
342 {COEF, 3}, /* 38 80 */
343 {DECIM9, DA + 1}, /* 39 p4 */
344 {COEF, 0}, /* 40 100 day hundreds */
345 {COEF, 1}, /* 41 200 */
346 {COEF2, 2}, /* 42 400 (not used) */
347 {COEF2, 3}, /* 43 800 (not used) */
348 {DECIM3, DA + 2}, /* 44 */
353 {IDLE, 0}, /* 49 p5 */
354 {MSCBIT, DUTS}, /* 50 dut+- */
355 {COEF, 0}, /* 51 10 year tens */
356 {COEF, 1}, /* 52 20 */
357 {COEF, 2}, /* 53 40 */
358 {COEF, 3}, /* 54 80 */
359 {MSC20, DST1}, /* 55 dst1 */
360 {MSCBIT, DUT1}, /* 56 0.1 dut */
361 {MSCBIT, DUT2}, /* 57 0.2 */
362 {MSC21, DUT4}, /* 58 0.4 QSY probe channel */
363 {MIN1, 0}, /* 59 p6 latch time */
364 {MIN2, 0} /* 60 leap second */
368 * BCD coefficients for maximum likelihood digit decode
370 #define P15 1. /* max positive number */
371 #define N15 -1. /* max negative number */
376 #define P9 (P15 / 4) /* mark (+1) */
377 #define N9 (N15 / 4) /* space (-1) */
380 {N9, N9, N9, N9}, /* 0 */
381 {P9, N9, N9, N9}, /* 1 */
382 {N9, P9, N9, N9}, /* 2 */
383 {P9, P9, N9, N9}, /* 3 */
384 {N9, N9, P9, N9}, /* 4 */
385 {P9, N9, P9, N9}, /* 5 */
386 {N9, P9, P9, N9}, /* 6 */
387 {P9, P9, P9, N9}, /* 7 */
388 {N9, N9, N9, P9}, /* 8 */
389 {P9, N9, N9, P9}, /* 9 */
390 {0, 0, 0, 0} /* backstop */
394 * Digits 0-6 (minute tens)
396 #define P6 (P15 / 3) /* mark (+1) */
397 #define N6 (N15 / 3) /* space (-1) */
400 {N6, N6, N6, 0}, /* 0 */
401 {P6, N6, N6, 0}, /* 1 */
402 {N6, P6, N6, 0}, /* 2 */
403 {P6, P6, N6, 0}, /* 3 */
404 {N6, N6, P6, 0}, /* 4 */
405 {P6, N6, P6, 0}, /* 5 */
406 {N6, P6, P6, 0}, /* 6 */
407 {0, 0, 0, 0} /* backstop */
411 * Digits 0-3 (day hundreds)
413 #define P3 (P15 / 2) /* mark (+1) */
414 #define N3 (N15 / 2) /* space (-1) */
417 {N3, N3, 0, 0}, /* 0 */
418 {P3, N3, 0, 0}, /* 1 */
419 {N3, P3, 0, 0}, /* 2 */
420 {P3, P3, 0, 0}, /* 3 */
421 {0, 0, 0, 0} /* backstop */
425 * Digits 0-2 (hour tens)
427 #define P2 (P15 / 2) /* mark (+1) */
428 #define N2 (N15 / 2) /* space (-1) */
431 {N2, N2, 0, 0}, /* 0 */
432 {P2, N2, 0, 0}, /* 1 */
433 {N2, P2, 0, 0}, /* 2 */
434 {0, 0, 0, 0} /* backstop */
438 * DST decode (DST2 DST1) for prettyprint
441 'S', /* 00 standard time */
442 'I', /* 01 set clock ahead at 0200 local */
443 'O', /* 10 set clock back at 0200 local */
444 'D' /* 11 daylight time */
448 * The decoding matrix consists of nine row vectors, one for each digit
449 * of the timecode. The digits are stored from least to most significant
450 * order. The maximum likelihood timecode is formed from the digits
451 * corresponding to the maximum likelihood values reading in the
452 * opposite order: yy ddd hh:mm.
455 int radix; /* radix (3, 4, 6, 10) */
456 int digit; /* current clock digit */
457 int mldigit; /* maximum likelihood digit */
458 int count; /* match count */
459 double digprb; /* max digit probability */
460 double digsnr; /* likelihood function (dB) */
461 double like[10]; /* likelihood integrator 0-9 */
465 * The station structure (sp) is used to acquire the minute pulse from
466 * WWV and/or WWVH. These stations are distinguished by the frequency
467 * used for the second and minute sync pulses, 1000 Hz for WWV and 1200
468 * Hz for WWVH. Other than frequency, the format is the same.
471 double epoch; /* accumulated epoch differences */
472 double maxeng; /* sync max energy */
473 double noieng; /* sync noise energy */
474 long pos; /* max amplitude position */
475 long lastpos; /* last max position */
476 long mepoch; /* minute synch epoch */
478 double amp; /* sync signal */
479 double syneng; /* sync signal max */
480 double synmax; /* sync signal max latched at 0 s */
481 double synsnr; /* sync signal SNR */
482 double metric; /* signal quality metric */
483 int reach; /* reachability register */
484 int count; /* bit counter */
485 int select; /* select bits */
486 char refid[5]; /* reference identifier */
490 * The channel structure (cp) is used to mitigate between channels.
493 int gain; /* audio gain */
494 struct sync wwv; /* wwv station */
495 struct sync wwvh; /* wwvh station */
499 * WWV unit control structure (up)
502 l_fp timestamp; /* audio sample timestamp */
503 l_fp tick; /* audio sample increment */
504 double phase, freq; /* logical clock phase and frequency */
505 double monitor; /* audio monitor point */
507 int fd_icom; /* ICOM file descriptor */
509 int errflg; /* error flags */
510 int watch; /* watchcat */
513 * Audio codec variables
515 double comp[SIZE]; /* decompanding table */
516 int port; /* codec port */
517 int gain; /* codec gain */
518 int mongain; /* codec monitor gain */
519 int clipcnt; /* sample clipped count */
522 * Variables used to establish basic system timing
524 int avgint; /* master time constant */
525 int yepoch; /* sync epoch */
526 int repoch; /* buffered sync epoch */
527 double epomax; /* second sync amplitude */
528 double eposnr; /* second sync SNR */
529 double irig; /* data I channel amplitude */
530 double qrig; /* data Q channel amplitude */
531 int datapt; /* 100 Hz ramp */
532 double datpha; /* 100 Hz VFO control */
533 int rphase; /* second sample counter */
534 long mphase; /* minute sample counter */
537 * Variables used to mitigate which channel to use
539 struct chan mitig[NCHAN]; /* channel data */
540 struct sync *sptr; /* station pointer */
541 int dchan; /* data channel */
542 int schan; /* probe channel */
543 int achan; /* active channel */
546 * Variables used by the clock state machine
548 struct decvec decvec[9]; /* decoding matrix */
549 int rsec; /* seconds counter */
550 int digcnt; /* count of digits synchronized */
553 * Variables used to estimate signal levels and bit/digit
556 double datsig; /* data signal max */
557 double datsnr; /* data signal SNR (dB) */
560 * Variables used to establish status and alarm conditions
562 int status; /* status bits */
563 int alarm; /* alarm flashers */
564 int misc; /* miscellaneous timecode bits */
565 int errcnt; /* data bit error counter */
569 * Function prototypes
571 static int wwv_start P((int, struct peer *));
572 static void wwv_shutdown P((int, struct peer *));
573 static void wwv_receive P((struct recvbuf *));
574 static void wwv_poll P((int, struct peer *));
577 * More function prototypes
579 static void wwv_epoch P((struct peer *));
580 static void wwv_rf P((struct peer *, double));
581 static void wwv_endpoc P((struct peer *, int));
582 static void wwv_rsec P((struct peer *, double));
583 static void wwv_qrz P((struct peer *, struct sync *, int));
584 static void wwv_corr4 P((struct peer *, struct decvec *,
585 double [], double [][4]));
586 static void wwv_gain P((struct peer *));
587 static void wwv_tsec P((struct peer *));
588 static int timecode P((struct wwvunit *, char *));
589 static double wwv_snr P((double, double));
590 static int carry P((struct decvec *));
591 static int wwv_newchan P((struct peer *));
592 static void wwv_newgame P((struct peer *));
593 static double wwv_metric P((struct sync *));
594 static void wwv_clock P((struct peer *));
596 static int wwv_qsy P((struct peer *, int));
599 static double qsy[NCHAN] = {2.5, 5, 10, 15, 20}; /* frequencies (MHz) */
604 struct refclock refclock_wwv = {
605 wwv_start, /* start up driver */
606 wwv_shutdown, /* shut down driver */
607 wwv_poll, /* transmit poll message */
608 noentry, /* not used (old wwv_control) */
609 noentry, /* initialize driver (not used) */
610 noentry, /* not used (old wwv_buginfo) */
611 NOFLAGS /* not used */
616 * wwv_start - open the devices and initialize data for processing
620 int unit, /* instance number (used by PCM) */
621 struct peer *peer /* peer structure pointer */
624 struct refclockproc *pp;
633 int fd; /* file descriptor */
635 double step; /* codec adjustment */
640 fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit);
649 * Allocate and initialize unit structure
651 if (!(up = (struct wwvunit *)emalloc(sizeof(struct wwvunit)))) {
655 memset(up, 0, sizeof(struct wwvunit));
657 pp->unitptr = (caddr_t)up;
658 pp->io.clock_recv = wwv_receive;
659 pp->io.srcclock = (caddr_t)peer;
662 if (!io_addclock(&pp->io)) {
669 * Initialize miscellaneous variables
671 peer->precision = PRECISION;
672 pp->clockdesc = DESCRIPTION;
675 * The companded samples are encoded sign-magnitude. The table
676 * contains all the 256 values in the interest of speed.
678 up->comp[0] = up->comp[OFFSET] = 0.;
679 up->comp[1] = 1.; up->comp[OFFSET + 1] = -1.;
680 up->comp[2] = 3.; up->comp[OFFSET + 2] = -3.;
682 for (i = 3; i < OFFSET; i++) {
683 up->comp[i] = up->comp[i - 1] + step;
684 up->comp[OFFSET + i] = -up->comp[i];
688 DTOLFP(1. / SECOND, &up->tick);
691 * Initialize the decoding matrix with the radix for each digit
694 up->decvec[MN].radix = 10; /* minutes */
695 up->decvec[MN + 1].radix = 6;
696 up->decvec[HR].radix = 10; /* hours */
697 up->decvec[HR + 1].radix = 3;
698 up->decvec[DA].radix = 10; /* days */
699 up->decvec[DA + 1].radix = 10;
700 up->decvec[DA + 2].radix = 4;
701 up->decvec[YR].radix = 10; /* years */
702 up->decvec[YR + 1].radix = 10;
706 * Initialize autotune if available. Note that the ICOM select
707 * code must be less than 128, so the high order bit can be used
708 * to select the line speed 0 (9600 bps) or 1 (1200 bps).
715 if (peer->ttl != 0) {
716 if (peer->ttl & 0x80)
717 up->fd_icom = icom_init("/dev/icom", B1200,
720 up->fd_icom = icom_init("/dev/icom", B9600,
722 if (up->fd_icom < 0) {
723 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
726 up->errflg = CEVNT_FAULT;
729 if (up->fd_icom > 0) {
730 if (wwv_qsy(peer, DCHAN) != 0) {
731 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
733 "icom: radio not found");
734 up->errflg = CEVNT_FAULT;
738 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
740 "icom: autotune enabled");
746 * Let the games begin.
754 * wwv_shutdown - shut down the clock
758 int unit, /* instance number (not used) */
759 struct peer *peer /* peer structure pointer */
762 struct refclockproc *pp;
766 up = (struct wwvunit *)pp->unitptr;
770 io_closeclock(&pp->io);
780 * wwv_receive - receive data from the audio device
782 * This routine reads input samples and adjusts the logical clock to
783 * track the A/D sample clock by dropping or duplicating codec samples.
784 * It also controls the A/D signal level with an AGC loop to mimimize
785 * quantization noise and avoid overload.
789 struct recvbuf *rbufp /* receive buffer structure pointer */
793 struct refclockproc *pp;
799 double sample; /* codec sample */
800 u_char *dpt; /* buffer pointer */
801 int bufcnt; /* buffer counter */
804 peer = (struct peer *)rbufp->recv_srcclock;
806 up = (struct wwvunit *)pp->unitptr;
809 * Main loop - read until there ain't no more. Note codec
810 * samples are bit-inverted.
812 DTOLFP((double)rbufp->recv_length / SECOND, <emp);
813 L_SUB(&rbufp->recv_time, <emp);
814 up->timestamp = rbufp->recv_time;
815 dpt = rbufp->recv_buffer;
816 for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
817 sample = up->comp[~*dpt++ & 0xff];
820 * Clip noise spikes greater than MAXAMP (6000) and
821 * record the number of clips to be used later by the
824 if (sample > MAXAMP) {
827 } else if (sample < -MAXAMP) {
833 * Variable frequency oscillator. The codec oscillator
834 * runs at the nominal rate of 8000 samples per second,
835 * or 125 us per sample. A frequency change of one unit
836 * results in either duplicating or deleting one sample
837 * per second, which results in a frequency change of
840 up->phase += up->freq / SECOND;
841 up->phase += FREQ_OFFSET / 1e6;
842 if (up->phase >= .5) {
844 } else if (up->phase < -.5) {
846 wwv_rf(peer, sample);
847 wwv_rf(peer, sample);
849 wwv_rf(peer, sample);
851 L_ADD(&up->timestamp, &up->tick);
855 * Set the input port and monitor gain for the next buffer.
857 if (pp->sloppyclockflag & CLK_FLAG2)
861 if (pp->sloppyclockflag & CLK_FLAG3)
862 up->mongain = MONGAIN;
869 * wwv_poll - called by the transmit procedure
871 * This routine keeps track of status. If no offset samples have been
872 * processed during a poll interval, a timeout event is declared. If
873 * errors have have occurred during the interval, they are reported as
878 int unit, /* instance number (not used) */
879 struct peer *peer /* peer structure pointer */
882 struct refclockproc *pp;
886 up = (struct wwvunit *)pp->unitptr;
887 if (pp->coderecv == pp->codeproc)
888 up->errflg = CEVNT_TIMEOUT;
890 refclock_report(peer, up->errflg);
897 * wwv_rf - process signals and demodulate to baseband
899 * This routine grooms and filters decompanded raw audio samples. The
900 * output signal is the 100-Hz filtered baseband data signal in
901 * quadrature phase. The routine also determines the minute synch epoch,
902 * as well as certain signal maxima, minima and related values.
904 * There are two 1-s ramps used by this program. Both count the 8000
905 * logical clock samples spanning exactly one second. The epoch ramp
906 * counts the samples starting at an arbitrary time. The rphase ramp
907 * counts the samples starting at the 5-ms second sync pulse found
908 * during the epoch ramp.
910 * There are two 1-m ramps used by this program. The mphase ramp counts
911 * the 480,000 logical clock samples spanning exactly one minute and
912 * starting at an arbitrary time. The rsec ramp counts the 60 seconds of
913 * the minute starting at the 800-ms minute sync pulse found during the
914 * mphase ramp. The rsec ramp drives the seconds state machine to
915 * determine the bits and digits of the timecode.
917 * Demodulation operations are based on three synthesized quadrature
918 * sinusoids: 100 Hz for the data signal, 1000 Hz for the WWV sync
919 * signal and 1200 Hz for the WWVH sync signal. These drive synchronous
920 * matched filters for the data signal (170 ms at 100 Hz), WWV minute
921 * sync signal (800 ms at 1000 Hz) and WWVH minute sync signal (800 ms
922 * at 1200 Hz). Two additional matched filters are switched in
923 * as required for the WWV second sync signal (5 cycles at 1000 Hz) and
924 * WWVH second sync signal (6 cycles at 1200 Hz).
928 struct peer *peer, /* peerstructure pointer */
929 double isig /* input signal */
932 struct refclockproc *pp;
934 struct sync *sp, *rp;
936 static double lpf[5]; /* 150-Hz lpf delay line */
937 double data; /* lpf output */
938 static double bpf[9]; /* 1000/1200-Hz bpf delay line */
939 double syncx; /* bpf output */
940 static double mf[41]; /* 1000/1200-Hz mf delay line */
941 double mfsync; /* mf output */
943 static int iptr; /* data channel pointer */
944 static double ibuf[DATSIZ]; /* data I channel delay line */
945 static double qbuf[DATSIZ]; /* data Q channel delay line */
947 static int jptr; /* sync channel pointer */
948 static int kptr; /* tick channel pointer */
950 static int csinptr; /* wwv channel phase */
951 static double cibuf[SYNSIZ]; /* wwv I channel delay line */
952 static double cqbuf[SYNSIZ]; /* wwv Q channel delay line */
953 static double ciamp; /* wwv I channel amplitude */
954 static double cqamp; /* wwv Q channel amplitude */
956 static double csibuf[TCKSIZ]; /* wwv I tick delay line */
957 static double csqbuf[TCKSIZ]; /* wwv Q tick delay line */
958 static double csiamp; /* wwv I tick amplitude */
959 static double csqamp; /* wwv Q tick amplitude */
961 static int hsinptr; /* wwvh channel phase */
962 static double hibuf[SYNSIZ]; /* wwvh I channel delay line */
963 static double hqbuf[SYNSIZ]; /* wwvh Q channel delay line */
964 static double hiamp; /* wwvh I channel amplitude */
965 static double hqamp; /* wwvh Q channel amplitude */
967 static double hsibuf[TCKSIZ]; /* wwvh I tick delay line */
968 static double hsqbuf[TCKSIZ]; /* wwvh Q tick delay line */
969 static double hsiamp; /* wwvh I tick amplitude */
970 static double hsqamp; /* wwvh Q tick amplitude */
972 static double epobuf[SECOND]; /* second sync comb filter */
973 static double epomax, nxtmax; /* second sync amplitude buffer */
974 static int epopos; /* epoch second sync position buffer */
976 static int iniflg; /* initialization flag */
977 int pdelay; /* propagation delay (samples) */
978 int epoch; /* comb filter index */
983 up = (struct wwvunit *)pp->unitptr;
987 memset((char *)lpf, 0, sizeof(lpf));
988 memset((char *)bpf, 0, sizeof(bpf));
989 memset((char *)mf, 0, sizeof(mf));
990 memset((char *)ibuf, 0, sizeof(ibuf));
991 memset((char *)qbuf, 0, sizeof(qbuf));
992 memset((char *)cibuf, 0, sizeof(cibuf));
993 memset((char *)cqbuf, 0, sizeof(cqbuf));
994 memset((char *)csibuf, 0, sizeof(csibuf));
995 memset((char *)csqbuf, 0, sizeof(csqbuf));
996 memset((char *)hibuf, 0, sizeof(hibuf));
997 memset((char *)hqbuf, 0, sizeof(hqbuf));
998 memset((char *)hsibuf, 0, sizeof(hsibuf));
999 memset((char *)hsqbuf, 0, sizeof(hsqbuf));
1000 memset((char *)epobuf, 0, sizeof(epobuf));
1004 * Baseband data demodulation. The 100-Hz subcarrier is
1005 * extracted using a 150-Hz IIR lowpass filter. This attenuates
1006 * the 1000/1200-Hz sync signals, as well as the 440-Hz and
1007 * 600-Hz tones and most of the noise and voice modulation
1010 * The subcarrier is transmitted 10 dB down from the carrier.
1011 * The DGAIN parameter can be adjusted for this and to
1012 * compensate for the radio audio response at 100 Hz.
1014 * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB
1015 * passband ripple, -50 dB stopband ripple.
1017 data = (lpf[4] = lpf[3]) * 8.360961e-01;
1018 data += (lpf[3] = lpf[2]) * -3.481740e+00;
1019 data += (lpf[2] = lpf[1]) * 5.452988e+00;
1020 data += (lpf[1] = lpf[0]) * -3.807229e+00;
1021 lpf[0] = isig * DGAIN - data;
1022 data = lpf[0] * 3.281435e-03
1023 + lpf[1] * -1.149947e-02
1024 + lpf[2] * 1.654858e-02
1025 + lpf[3] * -1.149947e-02
1026 + lpf[4] * 3.281435e-03;
1029 * The 100-Hz data signal is demodulated using a pair of
1030 * quadrature multipliers, matched filters and a phase lock
1031 * loop. The I and Q quadrature data signals are produced by
1032 * multiplying the filtered signal by 100-Hz sine and cosine
1033 * signals, respectively. The signals are processed by 170-ms
1034 * synchronous matched filters to produce the amplitude and
1035 * phase signals used by the demodulator. The signals are scaled
1036 * to produce unit energy at the maximum value.
1039 up->datapt = (up->datapt + IN100) % 80;
1040 dtemp = sintab[i] * data / (MS / 2. * DATCYC);
1041 up->irig -= ibuf[iptr];
1046 dtemp = sintab[i] * data / (MS / 2. * DATCYC);
1047 up->qrig -= qbuf[iptr];
1050 iptr = (iptr + 1) % DATSIZ;
1053 * Baseband sync demodulation. The 1000/1200 sync signals are
1054 * extracted using a 600-Hz IIR bandpass filter. This removes
1055 * the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz
1056 * tones and most of the noise and voice modulation components.
1058 * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB
1059 * passband ripple, -50 dB stopband ripple.
1061 syncx = (bpf[8] = bpf[7]) * 4.897278e-01;
1062 syncx += (bpf[7] = bpf[6]) * -2.765914e+00;
1063 syncx += (bpf[6] = bpf[5]) * 8.110921e+00;
1064 syncx += (bpf[5] = bpf[4]) * -1.517732e+01;
1065 syncx += (bpf[4] = bpf[3]) * 1.975197e+01;
1066 syncx += (bpf[3] = bpf[2]) * -1.814365e+01;
1067 syncx += (bpf[2] = bpf[1]) * 1.159783e+01;
1068 syncx += (bpf[1] = bpf[0]) * -4.735040e+00;
1069 bpf[0] = isig - syncx;
1070 syncx = bpf[0] * 8.203628e-03
1071 + bpf[1] * -2.375732e-02
1072 + bpf[2] * 3.353214e-02
1073 + bpf[3] * -4.080258e-02
1074 + bpf[4] * 4.605479e-02
1075 + bpf[5] * -4.080258e-02
1076 + bpf[6] * 3.353214e-02
1077 + bpf[7] * -2.375732e-02
1078 + bpf[8] * 8.203628e-03;
1081 * The 1000/1200 sync signals are demodulated using a pair of
1082 * quadrature multipliers and matched filters. However,
1083 * synchronous demodulation at these frequencies is impractical,
1084 * so only the signal amplitude is used. The I and Q quadrature
1085 * sync signals are produced by multiplying the filtered signal
1086 * by 1000-Hz (WWV) and 1200-Hz (WWVH) sine and cosine signals,
1087 * respectively. The WWV and WWVH signals are processed by 800-
1088 * ms synchronous matched filters and combined to produce the
1089 * minute sync signal and detect which one (or both) the WWV or
1090 * WWVH signal is present. The WWV and WWVH signals are also
1091 * processed by 5-ms synchronous matched filters and combined to
1092 * produce the second sync signal. The signals are scaled to
1093 * produce unit energy at the maximum value.
1095 * Note the master timing ramps, which run continuously. The
1096 * minute counter (mphase) counts the samples in the minute,
1097 * while the second counter (epoch) counts the samples in the
1100 up->mphase = (up->mphase + 1) % MINUTE;
1101 epoch = up->mphase % SECOND;
1107 csinptr = (csinptr + IN1000) % 80;
1109 dtemp = sintab[i] * syncx / (MS / 2.);
1110 ciamp -= cibuf[jptr];
1111 cibuf[jptr] = dtemp;
1113 csiamp -= csibuf[kptr];
1114 csibuf[kptr] = dtemp;
1118 dtemp = sintab[i] * syncx / (MS / 2.);
1119 cqamp -= cqbuf[jptr];
1120 cqbuf[jptr] = dtemp;
1122 csqamp -= csqbuf[kptr];
1123 csqbuf[kptr] = dtemp;
1126 sp = &up->mitig[up->achan].wwv;
1127 sp->amp = sqrt(ciamp * ciamp + cqamp * cqamp) / SYNCYC;
1128 if (!(up->status & MSYNC))
1129 wwv_qrz(peer, sp, (int)(pp->fudgetime1 * SECOND));
1135 hsinptr = (hsinptr + IN1200) % 80;
1137 dtemp = sintab[i] * syncx / (MS / 2.);
1138 hiamp -= hibuf[jptr];
1139 hibuf[jptr] = dtemp;
1141 hsiamp -= hsibuf[kptr];
1142 hsibuf[kptr] = dtemp;
1146 dtemp = sintab[i] * syncx / (MS / 2.);
1147 hqamp -= hqbuf[jptr];
1148 hqbuf[jptr] = dtemp;
1150 hsqamp -= hsqbuf[kptr];
1151 hsqbuf[kptr] = dtemp;
1154 rp = &up->mitig[up->achan].wwvh;
1155 rp->amp = sqrt(hiamp * hiamp + hqamp * hqamp) / SYNCYC;
1156 if (!(up->status & MSYNC))
1157 wwv_qrz(peer, rp, (int)(pp->fudgetime2 * SECOND));
1158 jptr = (jptr + 1) % SYNSIZ;
1159 kptr = (kptr + 1) % TCKSIZ;
1162 * The following section is called once per minute. It does
1163 * housekeeping and timeout functions and empties the dustbins.
1165 if (up->mphase == 0) {
1167 if (!(up->status & MSYNC)) {
1170 * If minute sync has not been acquired before
1171 * ACQSN timeout (6 min), or if no signal is
1172 * heard, the program cycles to the next
1173 * frequency and tries again.
1175 if (!wwv_newchan(peer))
1178 if (up->fd_icom > 0)
1179 wwv_qsy(peer, up->dchan);
1184 * If the leap bit is set, set the minute epoch
1185 * back one second so the station processes
1186 * don't miss a beat.
1188 if (up->status & LEPSEC) {
1189 up->mphase -= SECOND;
1191 up->mphase += MINUTE;
1197 * When the channel metric reaches threshold and the second
1198 * counter matches the minute epoch within the second, the
1199 * driver has synchronized to the station. The second number is
1200 * the remaining seconds until the next minute epoch, while the
1201 * sync epoch is zero. Watch out for the first second; if
1202 * already synchronized to the second, the buffered sync epoch
1205 * Note the guard interval is 200 ms; if for some reason the
1206 * clock drifts more than that, it might wind up in the wrong
1207 * second. If the maximum frequency error is not more than about
1208 * 1 PPM, the clock can go as much as two days while still in
1211 if (up->status & MSYNC) {
1213 } else if (up->sptr != NULL) {
1215 if (sp->metric >= TTHR && epoch == sp->mepoch % SECOND) {
1216 up->rsec = (60 - sp->mepoch / SECOND) % 60;
1218 up->status |= MSYNC;
1220 if (!(up->status & SSYNC))
1221 up->repoch = up->yepoch = epoch;
1223 up->repoch = up->yepoch;
1229 * The second sync pulse is extracted using 5-ms (40 sample) FIR
1230 * matched filters at 1000 Hz for WWV or 1200 Hz for WWVH. This
1231 * pulse is used for the most precise synchronization, since if
1232 * provides a resolution of one sample (125 us). The filters run
1233 * only if the station has been reliably determined.
1235 if (up->status & SELV) {
1236 pdelay = (int)(pp->fudgetime1 * SECOND);
1237 mfsync = sqrt(csiamp * csiamp + csqamp * csqamp) /
1239 } else if (up->status & SELH) {
1240 pdelay = (int)(pp->fudgetime2 * SECOND);
1241 mfsync = sqrt(hsiamp * hsiamp + hsqamp * hsqamp) /
1249 * Enhance the seconds sync pulse using a 1-s (8000-sample) comb
1250 * filter. Correct for the FIR matched filter delay, which is 5
1251 * ms for both the WWV and WWVH filters, and also for the
1252 * propagation delay. Once each second look for second sync. If
1253 * not in minute sync, fiddle the codec gain. Note the SNR is
1254 * computed from the maximum sample and the envelope of the
1255 * sample 6 ms before it, so if we slip more than a cycle the
1256 * SNR should plummet. The signal is scaled to produce unit
1257 * energy at the maximum value.
1259 dtemp = (epobuf[epoch] += (mfsync - epobuf[epoch]) /
1261 if (dtemp > epomax) {
1269 nxtmax = fabs(epobuf[j]);
1272 up->epomax = epomax;
1273 up->eposnr = wwv_snr(epomax, nxtmax);
1274 epopos -= pdelay + TCKCYC * MS;
1277 wwv_endpoc(peer, epopos);
1278 if (!(up->status & SSYNC))
1279 up->alarm |= SYNERR;
1281 if (!(up->status & MSYNC))
1288 * wwv_qrz - identify and acquire WWV/WWVH minute sync pulse
1290 * This routine implements a virtual station process used to acquire
1291 * minute sync and to mitigate among the ten frequency and station
1292 * combinations. During minute sync acquisition the process probes each
1293 * frequency and station in turn for the minute pulse, which
1294 * involves searching through the entire 480,000-sample minute. The
1295 * process finds the maximum signal and RMS noise plus signal. Then, the
1296 * actual noise is determined by subtracting the energy of the matched
1299 * Students of radar receiver technology will discover this algorithm
1300 * amounts to a range-gate discriminator. A valid pulse must have peak
1301 * amplitude at least QTHR (2500) and SNR at least QSNR (20) dB and the
1302 * difference between the current and previous epoch must be less than
1303 * AWND (20 ms). Note that the discriminator peak occurs about 800 ms
1304 * into the second, so the timing is retarded to the previous second
1309 struct peer *peer, /* peer structure pointer */
1310 struct sync *sp, /* sync channel structure */
1311 int pdelay /* propagation delay (samples) */
1314 struct refclockproc *pp;
1316 char tbuf[80]; /* monitor buffer */
1320 up = (struct wwvunit *)pp->unitptr;
1323 * Find the sample with peak amplitude, which defines the minute
1324 * epoch. Accumulate all samples to determine the total noise
1327 epoch = up->mphase - pdelay - SYNSIZ;
1330 if (sp->amp > sp->maxeng) {
1331 sp->maxeng = sp->amp;
1334 sp->noieng += sp->amp;
1337 * At the end of the minute, determine the epoch of the minute
1338 * sync pulse, as well as the difference between the current and
1339 * previous epoches due to the intrinsic frequency error plus
1340 * jitter. When calculating the SNR, subtract the pulse energy
1341 * from the total noise energy and then normalize.
1343 if (up->mphase == 0) {
1344 sp->synmax = sp->maxeng;
1345 sp->synsnr = wwv_snr(sp->synmax, (sp->noieng -
1346 sp->synmax) / MINUTE);
1348 sp->lastpos = sp->pos;
1349 epoch = (sp->pos - sp->lastpos) % MINUTE;
1351 if (sp->reach & (1 << AMAX))
1353 if (sp->synmax > ATHR && sp->synsnr > ASNR) {
1354 if (abs(epoch) < AWND * MS) {
1357 sp->mepoch = sp->lastpos = sp->pos;
1358 } else if (sp->count == 1) {
1359 sp->lastpos = sp->pos;
1362 if (up->watch > ACQSN)
1365 sp->metric = wwv_metric(sp);
1366 if (pp->sloppyclockflag & CLK_FLAG4) {
1368 "wwv8 %04x %3d %s %04x %.0f %.0f/%.1f %4ld %4ld",
1369 up->status, up->gain, sp->refid,
1370 sp->reach & 0xffff, sp->metric, sp->synmax,
1371 sp->synsnr, sp->pos % SECOND, epoch);
1372 record_clock_stats(&peer->srcadr, tbuf);
1375 printf("%s\n", tbuf);
1378 sp->maxeng = sp->noieng = 0;
1384 * wwv_endpoc - identify and acquire second sync pulse
1386 * This routine is called at the end of the second sync interval. It
1387 * determines the second sync epoch position within the second and
1388 * disciplines the sample clock using a frequency-lock loop (FLL).
1390 * Second sync is determined in the RF input routine as the maximum
1391 * over all 8000 samples in the second comb filter. To assure accurate
1392 * and reliable time and frequency discipline, this routine performs a
1393 * great deal of heavy-handed heuristic data filtering and grooming.
1397 struct peer *peer, /* peer structure pointer */
1398 int epopos /* epoch max position */
1401 struct refclockproc *pp;
1403 static int epoch_mf[3]; /* epoch median filter */
1404 static int tepoch; /* current second epoch */
1405 static int xepoch; /* last second epoch */
1406 static int zepoch; /* last run epoch */
1407 static int zcount; /* last run end time */
1408 static int scount; /* seconds counter */
1409 static int syncnt; /* run length counter */
1410 static int maxrun; /* longest run length */
1411 static int mepoch; /* longest run end epoch */
1412 static int mcount; /* longest run end time */
1413 static int avgcnt; /* averaging interval counter */
1414 static int avginc; /* averaging ratchet */
1415 static int iniflg; /* initialization flag */
1416 char tbuf[80]; /* monitor buffer */
1421 up = (struct wwvunit *)pp->unitptr;
1424 memset((char *)epoch_mf, 0, sizeof(epoch_mf));
1428 * If the signal amplitude or SNR fall below thresholds, dim the
1429 * second sync lamp and wait for hotter ions. If no stations are
1430 * heard, we are either in a probe cycle or the ions are really
1434 if (up->epomax < STHR || up->eposnr < SSNR) {
1435 up->status &= ~(SSYNC | FGATE);
1436 avgcnt = syncnt = maxrun = 0;
1439 if (!(up->status & (SELV | SELH)))
1443 * A three-stage median filter is used to help denoise the
1444 * second sync pulse. The median sample becomes the candidate
1447 epoch_mf[2] = epoch_mf[1];
1448 epoch_mf[1] = epoch_mf[0];
1449 epoch_mf[0] = epopos;
1450 if (epoch_mf[0] > epoch_mf[1]) {
1451 if (epoch_mf[1] > epoch_mf[2])
1452 tepoch = epoch_mf[1]; /* 0 1 2 */
1453 else if (epoch_mf[2] > epoch_mf[0])
1454 tepoch = epoch_mf[0]; /* 2 0 1 */
1456 tepoch = epoch_mf[2]; /* 0 2 1 */
1458 if (epoch_mf[1] < epoch_mf[2])
1459 tepoch = epoch_mf[1]; /* 2 1 0 */
1460 else if (epoch_mf[2] < epoch_mf[0])
1461 tepoch = epoch_mf[0]; /* 1 0 2 */
1463 tepoch = epoch_mf[2]; /* 1 2 0 */
1468 * If the epoch candidate is the same as the last one, increment
1469 * the run counter. If not, save the length, epoch and end
1470 * time of the current run for use later and reset the counter.
1471 * The epoch is considered valid if the run is at least SCMP
1472 * (10) s, the minute is synchronized and the interval since the
1473 * last epoch is not greater than the averaging interval. Thus,
1474 * after a long absence, the program will wait a full averaging
1475 * interval while the comb filter charges up and noise
1478 tmp2 = (tepoch - xepoch) % SECOND;
1481 if (syncnt > SCMP && up->status & MSYNC && (up->status &
1482 FGATE || scount - zcount <= up->avgint)) {
1483 up->status |= SSYNC;
1484 up->yepoch = tepoch;
1486 } else if (syncnt >= maxrun) {
1492 if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status & MSYNC))
1495 "wwv1 %04x %3d %4d %5.0f %5.1f %5d %4d %4d %4d",
1496 up->status, up->gain, tepoch, up->epomax,
1497 up->eposnr, tmp2, avgcnt, syncnt,
1499 record_clock_stats(&peer->srcadr, tbuf);
1502 printf("%s\n", tbuf);
1506 if (avgcnt < up->avgint) {
1512 * The sample clock frequency is disciplined using a first-order
1513 * feedback loop with time constant consistent with the Allan
1514 * intercept of typical computer clocks. During each averaging
1515 * interval the candidate epoch at the end of the longest run is
1516 * determined. If the longest run is zero, all epoches in the
1517 * interval are different, so the candidate epoch is the current
1518 * epoch. The frequency update is computed from the candidate
1519 * epoch difference (125-us units) and time difference (seconds)
1522 if (syncnt >= maxrun) {
1534 * The master clock runs at the codec sample frequency of 8000
1535 * Hz, so the intrinsic time resolution is 125 us. The frequency
1536 * resolution ranges from 18 PPM at the minimum averaging
1537 * interval of 8 s to 0.12 PPM at the maximum interval of 1024
1538 * s. An offset update is determined at the end of the longest
1539 * run in each averaging interval. The frequency adjustment is
1540 * computed from the difference between offset updates and the
1541 * interval between them.
1543 * The maximum frequency adjustment ranges from 187 PPM at the
1544 * minimum interval to 1.5 PPM at the maximum. If the adjustment
1545 * exceeds the maximum, the update is discarded and the
1546 * hysteresis counter is decremented. Otherwise, the frequency
1547 * is incremented by the adjustment, but clamped to the maximum
1548 * 187.5 PPM. If the update is less than half the maximum, the
1549 * hysteresis counter is incremented. If the counter increments
1550 * to +3, the averaging interval is doubled and the counter set
1551 * to zero; if it decrements to -3, the interval is halved and
1552 * the counter set to zero.
1554 dtemp = (mepoch - zepoch) % SECOND;
1555 if (up->status & FGATE) {
1556 if (abs(dtemp) < MAXFREQ * MINAVG) {
1557 up->freq += (dtemp / 2.) / ((mcount - zcount) *
1559 if (up->freq > MAXFREQ)
1561 else if (up->freq < -MAXFREQ)
1562 up->freq = -MAXFREQ;
1563 if (abs(dtemp) < MAXFREQ * MINAVG / 2.) {
1567 if (up->avgint < MAXAVG) {
1577 if (up->avgint > MINAVG) {
1584 if (pp->sloppyclockflag & CLK_FLAG4) {
1586 "wwv2 %04x %5.0f %5.1f %5d %4d %4d %4d %4.0f %7.2f",
1587 up->status, up->epomax, up->eposnr, mepoch,
1588 up->avgint, maxrun, mcount - zcount, dtemp,
1589 up->freq * 1e6 / SECOND);
1590 record_clock_stats(&peer->srcadr, tbuf);
1593 printf("%s\n", tbuf);
1598 * This is a valid update; set up for the next interval.
1600 up->status |= FGATE;
1603 avgcnt = syncnt = maxrun = 0;
1608 * wwv_epoch - epoch scanner
1610 * This routine extracts data signals from the 100-Hz subcarrier. It
1611 * scans the receiver second epoch to determine the signal amplitudes
1612 * and pulse timings. Receiver synchronization is determined by the
1613 * minute sync pulse detected in the wwv_rf() routine and the second
1614 * sync pulse detected in the wwv_epoch() routine. The transmitted
1615 * signals are delayed by the propagation delay, receiver delay and
1616 * filter delay of this program. Delay corrections are introduced
1617 * separately for WWV and WWVH.
1619 * Most communications radios use a highpass filter in the audio stages,
1620 * which can do nasty things to the subcarrier phase relative to the
1621 * sync pulses. Therefore, the data subcarrier reference phase is
1622 * disciplined using the hardlimited quadrature-phase signal sampled at
1623 * the same time as the in-phase signal. The phase tracking loop uses
1624 * phase adjustments of plus-minus one sample (125 us).
1628 struct peer *peer /* peer structure pointer */
1631 struct refclockproc *pp;
1634 static double sigmin, sigzer, sigone, engmax, engmin;
1637 up = (struct wwvunit *)pp->unitptr;
1640 * Find the maximum minute sync pulse energy for both the
1641 * WWV and WWVH stations. This will be used later for channel
1642 * and station mitigation. Also set the seconds epoch at 800 ms
1643 * well before the end of the second to make sure we never set
1644 * the epoch backwards.
1646 cp = &up->mitig[up->achan];
1647 if (cp->wwv.amp > cp->wwv.syneng)
1648 cp->wwv.syneng = cp->wwv.amp;
1649 if (cp->wwvh.amp > cp->wwvh.syneng)
1650 cp->wwvh.syneng = cp->wwvh.amp;
1651 if (up->rphase == 800 * MS)
1652 up->repoch = up->yepoch;
1655 * Use the signal amplitude at epoch 15 ms as the noise floor.
1656 * This gives a guard time of +-15 ms from the beginning of the
1657 * second until the second pulse rises at 30 ms. There is a
1658 * compromise here; we want to delay the sample as long as
1659 * possible to give the radio time to change frequency and the
1660 * AGC to stabilize, but as early as possible if the second
1661 * epoch is not exact.
1663 if (up->rphase == 15 * MS)
1664 sigmin = sigzer = sigone = up->irig;
1667 * Latch the data signal at 200 ms. Keep this around until the
1668 * end of the second. Use the signal energy as the peak to
1669 * compute the SNR. Use the Q sample to adjust the 100-Hz
1670 * reference oscillator phase.
1672 if (up->rphase == 200 * MS) {
1674 engmax = sqrt(up->irig * up->irig + up->qrig *
1676 up->datpha = up->qrig / up->avgint;
1677 if (up->datpha >= 0) {
1679 if (up->datapt >= 80)
1690 * Latch the data signal at 500 ms. Keep this around until the
1691 * end of the second.
1693 else if (up->rphase == 500 * MS)
1697 * At the end of the second crank the clock state machine and
1698 * adjust the codec gain. Note the epoch is buffered from the
1699 * center of the second in order to avoid jitter while the
1700 * seconds synch is diddling the epoch. Then, determine the true
1701 * offset and update the median filter in the driver interface.
1703 * Use the energy at the end of the second as the noise to
1704 * compute the SNR for the data pulse. This gives a better
1705 * measurement than the beginning of the second, especially when
1706 * returning from the probe channel. This gives a guard time of
1707 * 30 ms from the decay of the longest pulse to the rise of the
1711 if (up->mphase % SECOND == up->repoch) {
1712 up->status &= ~(DGATE | BGATE);
1713 engmin = sqrt(up->irig * up->irig + up->qrig *
1715 up->datsig = engmax;
1716 up->datsnr = wwv_snr(engmax, engmin);
1719 * If the amplitude or SNR is below threshold, average a
1720 * 0 in the the integrators; otherwise, average the
1721 * bipolar signal. This is done to avoid noise polution.
1723 if (engmax < DTHR || up->datsnr < DSNR) {
1724 up->status |= DGATE;
1729 wwv_rsec(peer, sigone - sigzer);
1731 if (up->status & (DGATE | BGATE))
1733 if (up->errcnt > MAXERR)
1734 up->alarm |= LOWERR;
1736 cp = &up->mitig[up->achan];
1738 cp->wwvh.syneng = 0;
1745 * wwv_rsec - process receiver second
1747 * This routine is called at the end of each receiver second to
1748 * implement the per-second state machine. The machine assembles BCD
1749 * digit bits, decodes miscellaneous bits and dances the leap seconds.
1751 * Normally, the minute has 60 seconds numbered 0-59. If the leap
1752 * warning bit is set, the last minute (1439) of 30 June (day 181 or 182
1753 * for leap years) or 31 December (day 365 or 366 for leap years) is
1754 * augmented by one second numbered 60. This is accomplished by
1755 * extending the minute interval by one second and teaching the state
1756 * machine to ignore it.
1760 struct peer *peer, /* peer structure pointer */
1764 static int iniflg; /* initialization flag */
1765 static double bcddld[4]; /* BCD data bits */
1766 static double bitvec[61]; /* bit integrator for misc bits */
1767 struct refclockproc *pp;
1770 struct sync *sp, *rp;
1771 char tbuf[80]; /* monitor buffer */
1775 up = (struct wwvunit *)pp->unitptr;
1778 memset((char *)bitvec, 0, sizeof(bitvec));
1782 * The bit represents the probability of a hit on zero (negative
1783 * values), a hit on one (positive values) or a miss (zero
1784 * value). The likelihood vector is the exponential average of
1785 * these probabilities. Only the bits of this vector
1786 * corresponding to the miscellaneous bits of the timecode are
1787 * used, but it's easier to do them all. After that, crank the
1788 * seconds state machine.
1792 bitvec[nsec] += (bit - bitvec[nsec]) / TCONST;
1793 sw = progx[nsec].sw;
1794 arg = progx[nsec].arg;
1797 * The minute state machine. Fly off to a particular section as
1798 * directed by the transition matrix and second number.
1803 * Ignore this second.
1805 case IDLE: /* 9, 45-49 */
1809 * Probe channel stuff
1811 * The WWV/H format contains data pulses in second 59 (position
1812 * identifier) and second 1, but not in second 0. The minute
1813 * sync pulse is contained in second 0. At the end of second 58
1814 * QSY to the probe channel, which rotates in turn over all
1815 * WWV/H frequencies. At the end of second 0 measure the minute
1816 * sync pulse. At the end of second 1 measure the data pulse and
1817 * QSY back to the data channel. Note that the actions commented
1818 * here happen at the end of the second numbered as shown.
1820 * At the end of second 0 save the minute sync amplitude latched
1821 * at 800 ms as the signal later used to calculate the SNR.
1824 cp = &up->mitig[up->achan];
1825 cp->wwv.synmax = cp->wwv.syneng;
1826 cp->wwvh.synmax = cp->wwvh.syneng;
1830 * At the end of second 1 use the minute sync amplitude latched
1831 * at 800 ms as the noise to calculate the SNR. If the minute
1832 * sync pulse and SNR are above thresholds and the data pulse
1833 * amplitude and SNR are above thresolds, shift a 1 into the
1834 * station reachability register; otherwise, shift a 0. The
1835 * number of 1 bits in the last six intervals is a component of
1836 * the channel metric computed by the wwv_metric() routine.
1837 * Finally, QSY back to the data channel.
1840 cp = &up->mitig[up->achan];
1846 sp->synsnr = wwv_snr(sp->synmax, sp->amp);
1848 if (sp->reach & (1 << AMAX))
1850 if (sp->synmax >= QTHR && sp->synsnr >= QSNR &&
1851 !(up->status & (DGATE | BGATE))) {
1855 sp->metric = wwv_metric(sp);
1861 rp->synsnr = wwv_snr(rp->synmax, rp->amp);
1863 if (rp->reach & (1 << AMAX))
1865 if (rp->synmax >= QTHR && rp->synsnr >= QSNR &&
1866 !(up->status & (DGATE | BGATE))) {
1870 rp->metric = wwv_metric(rp);
1871 if (pp->sloppyclockflag & CLK_FLAG4) {
1873 "wwv5 %04x %3d %4d %.0f/%.1f %.0f/%.1f %s %04x %.0f %.0f/%.1f %s %04x %.0f %.0f/%.1f",
1874 up->status, up->gain, up->yepoch,
1875 up->epomax, up->eposnr, up->datsig,
1877 sp->refid, sp->reach & 0xffff,
1878 sp->metric, sp->synmax, sp->synsnr,
1879 rp->refid, rp->reach & 0xffff,
1880 rp->metric, rp->synmax, rp->synsnr);
1881 record_clock_stats(&peer->srcadr, tbuf);
1884 printf("%s\n", tbuf);
1887 up->errcnt = up->digcnt = up->alarm = 0;
1890 * We now begin the minute scan. If not yet synchronized
1891 * to a station, restart if the units digit has not been
1892 * found within the DATA timeout (15 m) or if not
1893 * synchronized within the SYNCH timeout (40 m). After
1894 * synchronizing to a station, restart if no stations
1895 * are found within the PANIC timeout (2 days).
1897 if (up->status & INSYNC) {
1898 if (up->watch > PANIC) {
1903 if (!(up->status & DSYNC)) {
1904 if (up->watch > DATA) {
1909 if (up->watch > SYNCH) {
1916 if (up->fd_icom > 0)
1917 wwv_qsy(peer, up->dchan);
1922 * Save the bit probability in the BCD data vector at the index
1923 * given by the argument. Bits not used in the digit are forced
1926 case COEF1: /* 4-7 */
1930 case COEF: /* 10-13, 15-17, 20-23, 25-26,
1931 30-33, 35-38, 40-41, 51-54 */
1932 if (up->status & DSYNC)
1938 case COEF2: /* 18, 27-28, 42-43 */
1943 * Correlate coefficient vector with each valid digit vector and
1944 * save in decoding matrix. We step through the decoding matrix
1945 * digits correlating each with the coefficients and saving the
1946 * greatest and the next lower for later SNR calculation.
1948 case DECIM2: /* 29 */
1949 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd2);
1952 case DECIM3: /* 44 */
1953 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd3);
1956 case DECIM6: /* 19 */
1957 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd6);
1960 case DECIM9: /* 8, 14, 24, 34, 39 */
1961 wwv_corr4(peer, &up->decvec[arg], bcddld, bcd9);
1965 * Miscellaneous bits. If above the positive threshold, declare
1966 * 1; if below the negative threshold, declare 0; otherwise
1967 * raise the BGATE bit. The design is intended to avoid
1968 * integrating noise under low SNR conditions.
1970 case MSC20: /* 55 */
1971 wwv_corr4(peer, &up->decvec[YR + 1], bcddld, bcd9);
1974 case MSCBIT: /* 2-3, 50, 56-57 */
1975 if (bitvec[nsec] > BTHR) {
1976 if (!(up->misc & arg))
1977 up->alarm |= CMPERR;
1979 } else if (bitvec[nsec] < -BTHR) {
1981 up->alarm |= CMPERR;
1984 up->status |= BGATE;
1989 * Save the data channel gain, then QSY to the probe channel and
1990 * dim the seconds comb filters. The newchan() routine will
1991 * light them back up.
1993 case MSC21: /* 58 */
1994 if (bitvec[nsec] > BTHR) {
1995 if (!(up->misc & arg))
1996 up->alarm |= CMPERR;
1998 } else if (bitvec[nsec] < -BTHR) {
2000 up->alarm |= CMPERR;
2003 up->status |= BGATE;
2005 up->status &= ~(SELV | SELH);
2007 if (up->fd_icom > 0) {
2008 up->schan = (up->schan + 1) % NCHAN;
2009 wwv_qsy(peer, up->schan);
2011 up->mitig[up->achan].gain = up->gain;
2014 up->mitig[up->achan].gain = up->gain;
2021 * During second 59 the receiver and codec AGC are settling
2022 * down, so the data pulse is unusable as quality metric. If
2023 * LEPSEC is set on the last minute of 30 June or 31 December,
2024 * the transmitter and receiver insert an extra second (60) in
2025 * the timescale and the minute sync repeats the second. Once
2026 * leaps occurred at intervals of about 18 months, but the last
2027 * leap before the most recent leap in 1995 was in 1998.
2030 if (up->status & LEPSEC)
2036 up->status &= ~LEPSEC;
2042 if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status &
2045 "wwv3 %2d %04x %3d %4d %5.0f %5.1f %5.0f %5.1f %5.0f",
2046 nsec, up->status, up->gain, up->yepoch, up->epomax,
2047 up->eposnr, up->datsig, up->datsnr, bit);
2048 record_clock_stats(&peer->srcadr, tbuf);
2051 printf("%s\n", tbuf);
2054 pp->disp += AUDIO_PHI;
2058 * The radio clock is set if the alarm bits are all zero. After that,
2059 * the time is considered valid if the second sync bit is lit. It should
2060 * not be a surprise, especially if the radio is not tunable, that
2061 * sometimes no stations are above the noise and the integrators
2062 * discharge below the thresholds. We assume that, after a day of signal
2063 * loss, the minute sync epoch will be in the same second. This requires
2064 * the codec frequency be accurate within 6 PPM. Practical experience
2065 * shows the frequency typically within 0.1 PPM, so after a day of
2066 * signal loss, the time should be within 8.6 ms..
2070 struct peer *peer /* peer unit pointer */
2073 struct refclockproc *pp;
2075 l_fp offset; /* offset in NTP seconds */
2078 up = (struct wwvunit *)pp->unitptr;
2079 if (!(up->status & SSYNC))
2080 up->alarm |= SYNERR;
2082 up->alarm |= NINERR;
2084 up->status |= INSYNC;
2085 if (up->status & INSYNC && up->status & SSYNC) {
2086 if (up->misc & SECWAR)
2087 pp->leap = LEAP_ADDSECOND;
2089 pp->leap = LEAP_NOWARNING;
2090 pp->second = up->rsec;
2091 pp->minute = up->decvec[MN].digit + up->decvec[MN +
2093 pp->hour = up->decvec[HR].digit + up->decvec[HR +
2095 pp->day = up->decvec[DA].digit + up->decvec[DA +
2096 1].digit * 10 + up->decvec[DA + 2].digit * 100;
2097 pp->year = up->decvec[YR].digit + up->decvec[YR +
2101 if (!clocktime(pp->day, pp->hour, pp->minute,
2102 pp->second, GMT, up->timestamp.l_ui,
2103 &pp->yearstart, &offset.l_ui)) {
2104 up->errflg = CEVNT_BADTIME;
2108 pp->lastref = up->timestamp;
2109 refclock_process_offset(pp, offset,
2110 up->timestamp, PDELAY);
2111 refclock_receive(peer);
2114 pp->lencode = timecode(up, pp->a_lastcode);
2115 record_clock_stats(&peer->srcadr, pp->a_lastcode);
2118 printf("wwv: timecode %d %s\n", pp->lencode,
2125 * wwv_corr4 - determine maximum likelihood digit
2127 * This routine correlates the received digit vector with the BCD
2128 * coefficient vectors corresponding to all valid digits at the given
2129 * position in the decoding matrix. The maximum value corresponds to the
2130 * maximum likelihood digit, while the ratio of this value to the next
2131 * lower value determines the likelihood function. Note that, if the
2132 * digit is invalid, the likelihood vector is averaged toward a miss.
2136 struct peer *peer, /* peer unit pointer */
2137 struct decvec *vp, /* decoding table pointer */
2138 double data[], /* received data vector */
2139 double tab[][4] /* correlation vector array */
2142 struct refclockproc *pp;
2144 double topmax, nxtmax; /* metrics */
2145 double acc; /* accumulator */
2146 char tbuf[80]; /* monitor buffer */
2147 int mldigit; /* max likelihood digit */
2151 up = (struct wwvunit *)pp->unitptr;
2154 * Correlate digit vector with each BCD coefficient vector. If
2155 * any BCD digit bit is bad, consider all bits a miss. Until the
2156 * minute units digit has been resolved, don't to anything else.
2157 * Note the SNR is calculated as the ratio of the largest
2158 * likelihood value to the next largest likelihood value.
2161 topmax = nxtmax = -MAXAMP;
2162 for (i = 0; tab[i][0] != 0; i++) {
2164 for (j = 0; j < 4; j++)
2165 acc += data[j] * tab[i][j];
2166 acc = (vp->like[i] += (acc - vp->like[i]) / TCONST);
2171 } else if (acc > nxtmax) {
2175 vp->digprb = topmax;
2176 vp->digsnr = wwv_snr(topmax, nxtmax);
2179 * The current maximum likelihood digit is compared to the last
2180 * maximum likelihood digit. If different, the compare counter
2181 * and maximum likelihood digit are reset. When the compare
2182 * counter reaches the BCMP threshold (3), the digit is assumed
2183 * correct. When the compare counter of all nine digits have
2184 * reached threshold, the clock is assumed correct.
2186 * Note that the clock display digit is set before the compare
2187 * counter has reached threshold; however, the clock display is
2188 * not considered correct until all nine clock digits have
2189 * reached threshold. This is intended as eye candy, but avoids
2190 * mistakes when the signal is low and the SNR is very marginal.
2191 * once correctly set, the maximum likelihood digit is ignored
2192 * on the assumption the clock will always be correct unless for
2193 * some reason it drifts to a different second.
2195 vp->mldigit = mldigit;
2196 if (vp->digprb < BTHR || vp->digsnr < BSNR) {
2198 up->status |= BGATE;
2200 up->status |= DSYNC;
2201 if (vp->digit != mldigit) {
2203 up->alarm |= CMPERR;
2204 if (!(up->status & INSYNC))
2205 vp->digit = mldigit;
2207 if (vp->count < BCMP)
2213 if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status &
2216 "wwv4 %2d %04x %3d %4d %5.0f %2d %d %d %d %5.0f %5.1f",
2217 up->rsec - 1, up->status, up->gain, up->yepoch,
2218 up->epomax, vp->radix, vp->digit, vp->mldigit,
2219 vp->count, vp->digprb, vp->digsnr);
2220 record_clock_stats(&peer->srcadr, tbuf);
2223 printf("%s\n", tbuf);
2230 * wwv_tsec - transmitter minute processing
2232 * This routine is called at the end of the transmitter minute. It
2233 * implements a state machine that advances the logical clock subject to
2234 * the funny rules that govern the conventional clock and calendar.
2238 struct peer *peer /* driver structure pointer */
2241 struct refclockproc *pp;
2243 int minute, day, isleap;
2247 up = (struct wwvunit *)pp->unitptr;
2250 * Advance minute unit of the day. Don't propagate carries until
2251 * the unit minute digit has been found.
2253 temp = carry(&up->decvec[MN]); /* minute units */
2254 if (!(up->status & DSYNC))
2258 * Propagate carries through the day.
2260 if (temp == 0) /* carry minutes */
2261 temp = carry(&up->decvec[MN + 1]);
2262 if (temp == 0) /* carry hours */
2263 temp = carry(&up->decvec[HR]);
2265 temp = carry(&up->decvec[HR + 1]);
2268 * Decode the current minute and day. Set leap day if the
2269 * timecode leap bit is set on 30 June or 31 December. Set leap
2270 * minute if the last minute on leap day, but only if the clock
2271 * is syncrhronized. This code fails in 2400 AD.
2273 minute = up->decvec[MN].digit + up->decvec[MN + 1].digit *
2274 10 + up->decvec[HR].digit * 60 + up->decvec[HR +
2276 day = up->decvec[DA].digit + up->decvec[DA + 1].digit * 10 +
2277 up->decvec[DA + 2].digit * 100;
2280 * Set the leap bit on the last minute of the leap day.
2282 isleap = up->decvec[YR].digit & 0x3;
2283 if (up->misc & SECWAR && up->status & INSYNC) {
2284 if ((day == (isleap ? 182 : 183) || day == (isleap ?
2285 365 : 366)) && minute == 1439)
2286 up->status |= LEPSEC;
2290 * Roll the day if this the first minute and propagate carries
2297 while (carry(&up->decvec[HR]) != 0); /* advance to minute 0 */
2298 while (carry(&up->decvec[HR + 1]) != 0);
2300 temp = carry(&up->decvec[DA]); /* carry days */
2302 temp = carry(&up->decvec[DA + 1]);
2304 temp = carry(&up->decvec[DA + 2]);
2307 * Roll the year if this the first day and propagate carries
2308 * through the century.
2310 if (day != (isleap ? 365 : 366))
2314 while (carry(&up->decvec[DA]) != 1); /* advance to day 1 */
2315 while (carry(&up->decvec[DA + 1]) != 0);
2316 while (carry(&up->decvec[DA + 2]) != 0);
2317 temp = carry(&up->decvec[YR]); /* carry years */
2319 carry(&up->decvec[YR + 1]);
2324 * carry - process digit
2326 * This routine rotates a likelihood vector one position and increments
2327 * the clock digit modulo the radix. It returns the new clock digit or
2328 * zero if a carry occurred. Once synchronized, the clock digit will
2329 * match the maximum likelihood digit corresponding to that position.
2333 struct decvec *dp /* decoding table pointer */
2340 if (dp->digit == dp->radix)
2342 temp = dp->like[dp->radix - 1];
2343 for (j = dp->radix - 1; j > 0; j--)
2344 dp->like[j] = dp->like[j - 1];
2351 * wwv_snr - compute SNR or likelihood function
2355 double signal, /* signal */
2356 double noise /* noise */
2362 * This is a little tricky. Due to the way things are measured,
2363 * either or both the signal or noise amplitude can be negative
2364 * or zero. The intent is that, if the signal is negative or
2365 * zero, the SNR must always be zero. This can happen with the
2366 * subcarrier SNR before the phase has been aligned. On the
2367 * other hand, in the likelihood function the "noise" is the
2368 * next maximum down from the peak and this could be negative.
2369 * However, in this case the SNR is truly stupendous, so we
2370 * simply cap at MAXSNR dB (40).
2374 } else if (noise <= 0) {
2377 rval = 20. * log10(signal / noise);
2386 * wwv_newchan - change to new data channel
2388 * The radio actually appears to have ten channels, one channel for each
2389 * of five frequencies and each of two stations (WWV and WWVH), although
2390 * if not tunable only the DCHAN channel appears live. While the radio
2391 * is tuned to the working data channel frequency and station for most
2392 * of the minute, during seconds 59, 0 and 1 the radio is tuned to a
2393 * probe frequency in order to search for minute sync pulse and data
2394 * subcarrier from other transmitters.
2396 * The search for WWV and WWVH operates simultaneously, with WWV minute
2397 * sync pulse at 1000 Hz and WWVH at 1200 Hz. The probe frequency
2398 * rotates each minute over 2.5, 5, 10, 15 and 20 MHz in order and yes,
2399 * we all know WWVH is dark on 20 MHz, but few remember when WWV was lit
2402 * This routine selects the best channel using a metric computed from
2403 * the reachability register and minute pulse amplitude. Normally, the
2404 * award goes to the the channel with the highest metric; but, in case
2405 * of ties, the award goes to the channel with the highest minute sync
2406 * pulse amplitude and then to the highest frequency.
2408 * The routine performs an important squelch function to keep dirty data
2409 * from polluting the integrators. In order to consider a station valid,
2410 * the metric must be at least MTHR (13); otherwise, the station select
2411 * bits are cleared so the second sync is disabled and the data bit
2412 * integrators averaged to a miss.
2416 struct peer *peer /* peer structure pointer */
2419 struct refclockproc *pp;
2421 struct sync *sp, *rp;
2426 up = (struct wwvunit *)pp->unitptr;
2429 * Search all five station pairs looking for the channel with
2430 * maximum metric. If no station is found above thresholds, tune
2431 * to WWV on 15 MHz, set the reference ID to NONE and wait for
2437 for (i = 0; i < NCHAN; i++) {
2438 rp = &up->mitig[i].wwvh;
2440 if (dtemp >= rank) {
2445 rp = &up->mitig[i].wwv;
2447 if (dtemp >= rank) {
2455 * If the strongest signal is less than the MTHR threshold (13),
2456 * we are beneath the waves, so squelch the second sync. If the
2457 * strongest signal is greater than the threshold, tune to that
2458 * frequency and transmitter QTH.
2461 up->dchan = (up->dchan + 1) % NCHAN;
2462 up->status &= ~(SELV | SELH);
2466 up->status |= SELV | SELH;
2468 memcpy(&pp->refid, sp->refid, 4);
2469 peer->refid = pp->refid;
2475 * wwv_newgame - reset and start over
2477 * There are four conditions resulting in a new game:
2479 * 1 During initial acquisition (MSYNC dark) going 6 minutes (ACQSN)
2480 * without reliably finding the minute pulse (MSYNC lit).
2482 * 2 After finding the minute pulse (MSYNC lit), going 15 minutes
2483 * (DATA) without finding the unit seconds digit.
2485 * 3 After finding good data (DATA lit), going more than 40 minutes
2486 * (SYNCH) without finding station sync (INSYNC lit).
2488 * 4 After finding station sync (INSYNC lit), going more than 2 days
2489 * (PANIC) without finding any station.
2493 struct peer *peer /* peer structure pointer */
2496 struct refclockproc *pp;
2502 up = (struct wwvunit *)pp->unitptr;
2505 * Initialize strategic values. Note we set the leap bits
2506 * NOTINSYNC and the refid "NONE".
2508 peer->leap = LEAP_NOTINSYNC;
2509 up->watch = up->status = up->alarm = 0;
2510 up->avgint = MINAVG;
2512 up->gain = MAXGAIN / 2;
2515 * Initialize the station processes for audio gain, select bit,
2516 * station/frequency identifier and reference identifier. Start
2517 * probing at the next channel after the data channel.
2519 memset(up->mitig, 0, sizeof(up->mitig));
2520 for (i = 0; i < NCHAN; i++) {
2522 cp->gain = up->gain;
2523 cp->wwv.select = SELV;
2524 sprintf(cp->wwv.refid, "WV%.0f", floor(qsy[i]));
2525 cp->wwvh.select = SELH;
2526 sprintf(cp->wwvh.refid, "WH%.0f", floor(qsy[i]));
2528 up->dchan = (DCHAN + NCHAN - 1) % NCHAN;;
2530 up->achan = up->schan = up->dchan;
2532 if (up->fd_icom > 0)
2533 wwv_qsy(peer, up->dchan);
2538 * wwv_metric - compute station metric
2540 * The most significant bits represent the number of ones in the
2541 * station reachability register. The least significant bits represent
2542 * the minute sync pulse amplitude. The combined value is scaled 0-100.
2546 struct sync *sp /* station pointer */
2551 dtemp = sp->count * MAXAMP;
2552 if (sp->synmax < MAXAMP)
2553 dtemp += sp->synmax;
2555 dtemp += MAXAMP - 1;
2556 dtemp /= (AMAX + 1) * MAXAMP;
2557 return (dtemp * 100.);
2563 * wwv_qsy - Tune ICOM receiver
2565 * This routine saves the AGC for the current channel, switches to a new
2566 * channel and restores the AGC for that channel. If a tunable receiver
2567 * is not available, just fake it.
2571 struct peer *peer, /* peer structure pointer */
2572 int chan /* channel */
2576 struct refclockproc *pp;
2580 up = (struct wwvunit *)pp->unitptr;
2581 if (up->fd_icom > 0) {
2582 up->mitig[up->achan].gain = up->gain;
2583 rval = icom_freq(up->fd_icom, peer->ttl & 0x7f,
2586 up->gain = up->mitig[up->achan].gain;
2594 * timecode - assemble timecode string and length
2596 * Prettytime format - similar to Spectracom
2598 * sq yy ddd hh:mm:ss ld dut lset agc iden sig errs freq avgt
2600 * s sync indicator ('?' or ' ')
2601 * q error bits (hex 0-F)
2602 * yyyy year of century
2606 * ss second of minute)
2607 * l leap second warning (' ' or 'L')
2608 * d DST state ('S', 'D', 'I', or 'O')
2609 * dut DUT sign and magnitude (0.1 s)
2610 * lset minutes since last clock update
2611 * agc audio gain (0-255)
2612 * iden reference identifier (station and frequency)
2613 * sig signal quality (0-100)
2614 * errs bit errors in last minute
2615 * freq frequency offset (PPM)
2616 * avgt averaging time (s)
2620 struct wwvunit *up, /* driver structure pointer */
2621 char *ptr /* target string */
2625 int year, day, hour, minute, second, dut;
2626 char synchar, leapchar, dst;
2631 * Common fixed-format fields
2633 synchar = (up->status & INSYNC) ? ' ' : '?';
2634 year = up->decvec[YR].digit + up->decvec[YR + 1].digit * 10 +
2636 day = up->decvec[DA].digit + up->decvec[DA + 1].digit * 10 +
2637 up->decvec[DA + 2].digit * 100;
2638 hour = up->decvec[HR].digit + up->decvec[HR + 1].digit * 10;
2639 minute = up->decvec[MN].digit + up->decvec[MN + 1].digit * 10;
2641 leapchar = (up->misc & SECWAR) ? 'L' : ' ';
2642 dst = dstcod[(up->misc >> 4) & 0x3];
2643 dut = up->misc & 0x7;
2644 if (!(up->misc & DUTS))
2646 sprintf(ptr, "%c%1X", synchar, up->alarm);
2647 sprintf(cptr, " %4d %03d %02d:%02d:%02d %c%c %+d",
2648 year, day, hour, minute, second, leapchar, dst, dut);
2652 * Specific variable-format fields
2655 sprintf(cptr, " %d %d %s %.0f %d %.1f %d", up->watch,
2656 up->mitig[up->dchan].gain, sp->refid, sp->metric,
2657 up->errcnt, up->freq / SECOND * 1e6, up->avgint);
2659 return (strlen(ptr));
2664 * wwv_gain - adjust codec gain
2666 * This routine is called at the end of each second. During the second
2667 * the number of signal clips above the MAXAMP threshold (6000). If
2668 * there are no clips, the gain is bumped up; if there are more than
2669 * MAXCLP clips (100), it is bumped down. The decoder is relatively
2670 * insensitive to amplitude, so this crudity works just peachy. The
2671 * input port is set and the error flag is cleared, mostly to be ornery.
2675 struct peer *peer /* peer structure pointer */
2678 struct refclockproc *pp;
2682 up = (struct wwvunit *)pp->unitptr;
2685 * Apparently, the codec uses only the high order bits of the
2686 * gain control field. Thus, it may take awhile for changes to
2687 * wiggle the hardware bits.
2689 if (up->clipcnt == 0) {
2691 if (up->gain > MAXGAIN)
2693 } else if (up->clipcnt > MAXCLP) {
2698 audio_gain(up->gain, up->mongain, up->port);
2708 int refclock_wwv_bs;
2709 #endif /* REFCLOCK */