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