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