PageRenderTime 121ms CodeModel.GetById 20ms app.highlight 87ms RepoModel.GetById 1ms app.codeStats 1ms

/contrib/ntp/ntpd/refclock_wwv.c

https://bitbucket.org/freebsd/freebsd-head/
C | 2709 lines | 1594 code | 204 blank | 911 comment | 264 complexity | d5531729a944f68437d2d9bd4047b8da MD5 | raw file

Large files files are truncated, but you can click here to view the full file

   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 receives, demodulates and decodes the radio signals when
  43 * connected to the audio codec of a workstation running Solaris, SunOS
  44 * FreeBSD or Linux, and with a little help, other workstations with
  45 * similar codecs or sound cards. In this implementation, only one audio
  46 * driver and codec can be supported on a single machine.
  47 *
  48 * The demodulation and decoding algorithms used in this driver are
  49 * based on those developed for the TAPR DSP93 development board and the
  50 * TI 320C25 digital signal processor described in: Mills, D.L. A
  51 * precision radio clock for WWV transmissions. Electrical Engineering
  52 * Report 97-8-1, University of Delaware, August 1997, 25 pp., available
  53 * from www.eecis.udel.edu/~mills/reports.html. The algorithms described
  54 * in this report have been modified somewhat to improve performance
  55 * under weak signal conditions and to provide an automatic station
  56 * identification feature.
  57 *
  58 * The ICOM code is normally compiled in the driver. It isn't used,
  59 * unless the mode keyword on the server configuration command specifies
  60 * a nonzero ICOM ID select code. The C-IV trace is turned on if the
  61 * debug level is greater than one.
  62 *
  63 * Fudge factors
  64 *
  65 * Fudge flag4 causes the dubugging output described above to be
  66 * recorded in the clockstats file. Fudge flag2 selects the audio input
  67 * port, where 0 is the mike port (default) and 1 is the line-in port.
  68 * It does not seem useful to select the compact disc player port. Fudge
  69 * flag3 enables audio monitoring of the input signal. For this purpose,
  70 * the monitor gain is set to a default value.
  71 */
  72/*
  73 * General definitions. These ordinarily do not need to be changed.
  74 */
  75#define	DEVICE_AUDIO	"/dev/audio" /* audio device name */
  76#define	AUDIO_BUFSIZ	320	/* audio buffer size (50 ms) */
  77#define	PRECISION	(-10)	/* precision assumed (about 1 ms) */
  78#define	DESCRIPTION	"WWV/H Audio Demodulator/Decoder" /* WRU */
  79#define SECOND		8000	/* second epoch (sample rate) (Hz) */
  80#define MINUTE		(SECOND * 60) /* minute epoch */
  81#define OFFSET		128	/* companded sample offset */
  82#define SIZE		256	/* decompanding table size */
  83#define	MAXAMP		6000.	/* max signal level reference */
  84#define	MAXCLP		100	/* max clips above reference per s */
  85#define MAXSNR		40.	/* max SNR reference */
  86#define MAXFREQ		1.5	/* max frequency tolerance (187 PPM) */
  87#define DATCYC		170	/* data filter cycles */
  88#define DATSIZ		(DATCYC * MS) /* data filter size */
  89#define SYNCYC		800	/* minute filter cycles */
  90#define SYNSIZ		(SYNCYC * MS) /* minute filter size */
  91#define TCKCYC		5	/* tick filter cycles */
  92#define TCKSIZ		(TCKCYC * MS) /* tick filter size */
  93#define NCHAN		5	/* number of radio channels */
  94#define	AUDIO_PHI	5e-6	/* dispersion growth factor */
  95
  96/*
  97 * Tunable parameters. The DGAIN parameter can be changed to fit the
  98 * audio response of the radio at 100 Hz. The WWV/WWVH data subcarrier
  99 * is transmitted at about 20 percent percent modulation; the matched
 100 * filter boosts it by a factor of 17 and the receiver response does
 101 * what it does. The compromise value works for ICOM radios. If the
 102 * radio is not tunable, the DCHAN parameter can be changed to fit the
 103 * expected best propagation frequency: higher if further from the
 104 * transmitter, lower if nearer. The compromise value works for the US
 105 * right coast. The FREQ_OFFSET parameter can be used as a frequency
 106 * vernier to correct codec requency if greater than MAXFREQ.
 107 */
 108#define DCHAN		3	/* default radio channel (15 Mhz) */
 109#define DGAIN		5.	/* subcarrier gain */
 110#define	FREQ_OFFSET	0.	/* codec frequency correction (PPM) */
 111
 112/*
 113 * General purpose status bits (status)
 114 *
 115 * SELV and/or SELH are set when WWV or WWVH have been heard and cleared
 116 * on signal loss. SSYNC is set when the second sync pulse has been
 117 * acquired and cleared by signal loss. MSYNC is set when the minute
 118 * sync pulse has been acquired. DSYNC is set when the units digit has
 119 * has reached the threshold and INSYNC is set when all nine digits have
 120 * reached the threshold. The MSYNC, DSYNC and INSYNC bits are cleared
 121 * only by timeout, upon which the driver starts over from scratch.
 122 *
 123 * DGATE is lit if the data bit amplitude or SNR is below thresholds and
 124 * BGATE is lit if the pulse width amplitude or SNR is below thresolds.
 125 * LEPSEC is set during the last minute of the leap day. At the end of
 126 * this minute the driver inserts second 60 in the seconds state machine
 127 * and the minute sync slips a second.
 128 */
 129#define MSYNC		0x0001	/* minute epoch sync */
 130#define SSYNC		0x0002	/* second epoch sync */
 131#define DSYNC		0x0004	/* minute units sync */
 132#define INSYNC		0x0008	/* clock synchronized */
 133#define FGATE		0x0010	/* frequency gate */
 134#define DGATE		0x0020	/* data pulse amplitude error */
 135#define BGATE		0x0040	/* data pulse width error */
 136#define LEPSEC		0x1000	/* leap minute */
 137
 138/*
 139 * Station scoreboard bits
 140 *
 141 * These are used to establish the signal quality for each of the five
 142 * frequencies and two stations.
 143 */
 144#define SELV		0x0100	/* WWV station select */
 145#define SELH		0x0200	/* WWVH station select */
 146
 147/*
 148 * Alarm status bits (alarm)
 149 *
 150 * These bits indicate various alarm conditions, which are decoded to
 151 * form the quality character included in the timecode.
 152 */
 153#define CMPERR		1	/* digit or misc bit compare error */
 154#define LOWERR		2	/* low bit or digit amplitude or SNR */
 155#define NINERR		4	/* less than nine digits in minute */
 156#define SYNERR		8	/* not tracking second sync */
 157
 158/*
 159 * Watchcat timeouts (watch)
 160 *
 161 * If these timeouts expire, the status bits are mashed to zero and the
 162 * driver starts from scratch. Suitably more refined procedures may be
 163 * developed in future. All these are in minutes.
 164 */
 165#define ACQSN		6	/* station acquisition timeout */
 166#define DATA		15	/* unit minutes timeout */
 167#define SYNCH		40	/* station sync timeout */
 168#define PANIC		(2 * 1440) /* panic timeout */
 169
 170/*
 171 * Thresholds. These establish the minimum signal level, minimum SNR and
 172 * maximum jitter thresholds which establish the error and false alarm
 173 * rates of the driver. The values defined here may be on the
 174 * adventurous side in the interest of the highest sensitivity.
 175 */
 176#define MTHR		13.	/* minute sync gate (percent) */
 177#define TTHR		50.	/* minute sync threshold (percent) */
 178#define AWND		20	/* minute sync jitter threshold (ms) */
 179#define ATHR		2500.	/* QRZ minute sync threshold */
 180#define ASNR		20.	/* QRZ minute sync SNR threshold (dB) */
 181#define QTHR		2500.	/* QSY minute sync threshold */
 182#define QSNR		20.	/* QSY minute sync SNR threshold (dB) */
 183#define STHR		2500.	/* second sync threshold */
 184#define	SSNR		15.	/* second sync SNR threshold (dB) */
 185#define SCMP		10 	/* second sync compare threshold */
 186#define DTHR		1000.	/* bit threshold */
 187#define DSNR		10.	/* bit SNR threshold (dB) */
 188#define AMIN		3	/* min bit count */
 189#define AMAX		6	/* max bit count */
 190#define BTHR		1000.	/* digit threshold */
 191#define BSNR		3.	/* digit likelihood threshold (dB) */
 192#define BCMP		3	/* digit compare threshold */
 193#define	MAXERR		40	/* maximum error alarm */
 194
 195/*
 196 * Tone frequency definitions. The increments are for 4.5-deg sine
 197 * table.
 198 */
 199#define MS		(SECOND / 1000) /* samples per millisecond */
 200#define IN100		((100 * 80) / SECOND) /* 100 Hz increment */
 201#define IN1000		((1000 * 80) / SECOND) /* 1000 Hz increment */
 202#define IN1200		((1200 * 80) / SECOND) /* 1200 Hz increment */
 203
 204/*
 205 * Acquisition and tracking time constants
 206 */
 207#define MINAVG		8	/* min averaging time */
 208#define MAXAVG		1024	/* max averaging time */
 209#define FCONST		3	/* frequency time constant */
 210#define TCONST		16	/* data bit/digit time constant */
 211
 212/*
 213 * Miscellaneous status bits (misc)
 214 *
 215 * These bits correspond to designated bits in the WWV/H timecode. The
 216 * bit probabilities are exponentially averaged over several minutes and
 217 * processed by a integrator and threshold.
 218 */
 219#define DUT1		0x01	/* 56 DUT .1 */
 220#define DUT2		0x02	/* 57 DUT .2 */
 221#define DUT4		0x04	/* 58 DUT .4 */
 222#define DUTS		0x08	/* 50 DUT sign */
 223#define DST1		0x10	/* 55 DST1 leap warning */
 224#define DST2		0x20	/* 2 DST2 DST1 delayed one day */
 225#define SECWAR		0x40	/* 3 leap second warning */
 226
 227/*
 228 * The on-time synchronization point for the driver is the second epoch
 229 * sync pulse produced by the FIR matched filters. As the 5-ms delay of
 230 * these filters is compensated, the program delay is 1.1 ms due to the
 231 * 600-Hz IIR bandpass filter. The measured receiver delay is 4.7 ms and
 232 * the codec delay less than 0.2 ms. The additional propagation delay
 233 * specific to each receiver location can be programmed in the fudge
 234 * time1 and time2 values for WWV and WWVH, respectively.
 235 */
 236#define PDELAY	(.0011 + .0047 + .0002)	/* net system delay (s) */
 237
 238/*
 239 * Table of sine values at 4.5-degree increments. This is used by the
 240 * synchronous matched filter demodulators.
 241 */
 242double sintab[] = {
 243 0.000000e+00,  7.845910e-02,  1.564345e-01,  2.334454e-01, /* 0-3 */
 244 3.090170e-01,  3.826834e-01,  4.539905e-01,  5.224986e-01, /* 4-7 */
 245 5.877853e-01,  6.494480e-01,  7.071068e-01,  7.604060e-01, /* 8-11 */
 246 8.090170e-01,  8.526402e-01,  8.910065e-01,  9.238795e-01, /* 12-15 */
 247 9.510565e-01,  9.723699e-01,  9.876883e-01,  9.969173e-01, /* 16-19 */
 248 1.000000e+00,  9.969173e-01,  9.876883e-01,  9.723699e-01, /* 20-23 */
 249 9.510565e-01,  9.238795e-01,  8.910065e-01,  8.526402e-01, /* 24-27 */
 250 8.090170e-01,  7.604060e-01,  7.071068e-01,  6.494480e-01, /* 28-31 */
 251 5.877853e-01,  5.224986e-01,  4.539905e-01,  3.826834e-01, /* 32-35 */
 252 3.090170e-01,  2.334454e-01,  1.564345e-01,  7.845910e-02, /* 36-39 */
 253-0.000000e+00, -7.845910e-02, -1.564345e-01, -2.334454e-01, /* 40-43 */
 254-3.090170e-01, -3.826834e-01, -4.539905e-01, -5.224986e-01, /* 44-47 */
 255-5.877853e-01, -6.494480e-01, -7.071068e-01, -7.604060e-01, /* 48-51 */
 256-8.090170e-01, -8.526402e-01, -8.910065e-01, -9.238795e-01, /* 52-55 */
 257-9.510565e-01, -9.723699e-01, -9.876883e-01, -9.969173e-01, /* 56-59 */
 258-1.000000e+00, -9.969173e-01, -9.876883e-01, -9.723699e-01, /* 60-63 */
 259-9.510565e-01, -9.238795e-01, -8.910065e-01, -8.526402e-01, /* 64-67 */
 260-8.090170e-01, -7.604060e-01, -7.071068e-01, -6.494480e-01, /* 68-71 */
 261-5.877853e-01, -5.224986e-01, -4.539905e-01, -3.826834e-01, /* 72-75 */
 262-3.090170e-01, -2.334454e-01, -1.564345e-01, -7.845910e-02, /* 76-79 */
 263 0.000000e+00};						    /* 80 */
 264
 265/*
 266 * Decoder operations at the end of each second are driven by a state
 267 * machine. The transition matrix consists of a dispatch table indexed
 268 * by second number. Each entry in the table contains a case switch
 269 * number and argument.
 270 */
 271struct progx {
 272	int sw;			/* case switch number */
 273	int arg;		/* argument */
 274};
 275
 276/*
 277 * Case switch numbers
 278 */
 279#define IDLE		0	/* no operation */
 280#define COEF		1	/* BCD bit */
 281#define COEF1		2	/* BCD bit for minute unit */
 282#define COEF2		3	/* BCD bit not used */
 283#define DECIM9		4	/* BCD digit 0-9 */
 284#define DECIM6		5	/* BCD digit 0-6 */
 285#define DECIM3		6	/* BCD digit 0-3 */
 286#define DECIM2		7	/* BCD digit 0-2 */
 287#define MSCBIT		8	/* miscellaneous bit */
 288#define MSC20		9	/* miscellaneous bit */		
 289#define MSC21		10	/* QSY probe channel */		
 290#define MIN1		11	/* latch time */		
 291#define MIN2		12	/* leap second */
 292#define SYNC2		13	/* latch minute sync pulse */		
 293#define SYNC3		14	/* latch data pulse */		
 294
 295/*
 296 * Offsets in decoding matrix
 297 */
 298#define MN		0	/* minute digits (2) */
 299#define HR		2	/* hour digits (2) */
 300#define DA		4	/* day digits (3) */
 301#define YR		7	/* year digits (2) */
 302
 303struct progx progx[] = {
 304	{SYNC2,	0},		/* 0 latch minute sync pulse */
 305	{SYNC3,	0},		/* 1 latch data pulse */
 306	{MSCBIT, DST2},		/* 2 dst2 */
 307	{MSCBIT, SECWAR},	/* 3 lw */
 308	{COEF,	0},		/* 4 1 year units */
 309	{COEF,	1},		/* 5 2 */
 310	{COEF,	2},		/* 6 4 */
 311	{COEF,	3},		/* 7 8 */
 312	{DECIM9, YR},		/* 8 */
 313	{IDLE,	0},		/* 9 p1 */
 314	{COEF1,	0},		/* 10 1 minute units */
 315	{COEF1,	1},		/* 11 2 */
 316	{COEF1,	2},		/* 12 4 */
 317	{COEF1,	3},		/* 13 8 */
 318	{DECIM9, MN},		/* 14 */
 319	{COEF,	0},		/* 15 10 minute tens */
 320	{COEF,	1},		/* 16 20 */
 321	{COEF,	2},		/* 17 40 */
 322	{COEF2,	3},		/* 18 80 (not used) */
 323	{DECIM6, MN + 1},	/* 19 p2 */
 324	{COEF,	0},		/* 20 1 hour units */
 325	{COEF,	1},		/* 21 2 */
 326	{COEF,	2},		/* 22 4 */
 327	{COEF,	3},		/* 23 8 */
 328	{DECIM9, HR},		/* 24 */
 329	{COEF,	0},		/* 25 10 hour tens */
 330	{COEF,	1},		/* 26 20 */
 331	{COEF2,	2},		/* 27 40 (not used) */
 332	{COEF2,	3},		/* 28 80 (not used) */
 333	{DECIM2, HR + 1},	/* 29 p3 */
 334	{COEF,	0},		/* 30 1 day units */
 335	{COEF,	1},		/* 31 2 */
 336	{COEF,	2},		/* 32 4 */
 337	{COEF,	3},		/* 33 8 */
 338	{DECIM9, DA},		/* 34 */
 339	{COEF,	0},		/* 35 10 day tens */
 340	{COEF,	1},		/* 36 20 */
 341	{COEF,	2},		/* 37 40 */
 342	{COEF,	3},		/* 38 80 */
 343	{DECIM9, DA + 1},	/* 39 p4 */
 344	{COEF,	0},		/* 40 100 day hundreds */
 345	{COEF,	1},		/* 41 200 */
 346	{COEF2,	2},		/* 42 400 (not used) */
 347	{COEF2,	3},		/* 43 800 (not used) */
 348	{DECIM3, DA + 2},	/* 44 */
 349	{IDLE,	0},		/* 45 */
 350	{IDLE,	0},		/* 46 */
 351	{IDLE,	0},		/* 47 */
 352	{IDLE,	0},		/* 48 */
 353	{IDLE,	0},		/* 49 p5 */
 354	{MSCBIT, DUTS},		/* 50 dut+- */
 355	{COEF,	0},		/* 51 10 year tens */
 356	{COEF,	1},		/* 52 20 */
 357	{COEF,	2},		/* 53 40 */
 358	{COEF,	3},		/* 54 80 */
 359	{MSC20, DST1},		/* 55 dst1 */
 360	{MSCBIT, DUT1},		/* 56 0.1 dut */
 361	{MSCBIT, DUT2},		/* 57 0.2 */
 362	{MSC21, DUT4},		/* 58 0.4 QSY probe channel */
 363	{MIN1,	0},		/* 59 p6 latch time */
 364	{MIN2,	0}		/* 60 leap second */
 365};
 366
 367/*
 368 * BCD coefficients for maximum likelihood digit decode
 369 */
 370#define P15	1.		/* max positive number */
 371#define N15	-1.		/* max negative number */
 372
 373/*
 374 * Digits 0-9
 375 */
 376#define P9	(P15 / 4)	/* mark (+1) */
 377#define N9	(N15 / 4)	/* space (-1) */
 378
 379double bcd9[][4] = {
 380	{N9, N9, N9, N9}, 	/* 0 */
 381	{P9, N9, N9, N9}, 	/* 1 */
 382	{N9, P9, N9, N9}, 	/* 2 */
 383	{P9, P9, N9, N9}, 	/* 3 */
 384	{N9, N9, P9, N9}, 	/* 4 */
 385	{P9, N9, P9, N9}, 	/* 5 */
 386	{N9, P9, P9, N9}, 	/* 6 */
 387	{P9, P9, P9, N9}, 	/* 7 */
 388	{N9, N9, N9, P9}, 	/* 8 */
 389	{P9, N9, N9, P9}, 	/* 9 */
 390	{0, 0, 0, 0}		/* backstop */
 391};
 392
 393/*
 394 * Digits 0-6 (minute tens)
 395 */
 396#define P6	(P15 / 3)	/* mark (+1) */
 397#define N6	(N15 / 3)	/* space (-1) */
 398
 399double bcd6[][4] = {
 400	{N6, N6, N6, 0}, 	/* 0 */
 401	{P6, N6, N6, 0}, 	/* 1 */
 402	{N6, P6, N6, 0}, 	/* 2 */
 403	{P6, P6, N6, 0}, 	/* 3 */
 404	{N6, N6, P6, 0}, 	/* 4 */
 405	{P6, N6, P6, 0}, 	/* 5 */
 406	{N6, P6, P6, 0}, 	/* 6 */
 407	{0, 0, 0, 0}		/* backstop */
 408};
 409
 410/*
 411 * Digits 0-3 (day hundreds)
 412 */
 413#define P3	(P15 / 2)	/* mark (+1) */
 414#define N3	(N15 / 2)	/* space (-1) */
 415
 416double bcd3[][4] = {
 417	{N3, N3, 0, 0}, 	/* 0 */
 418	{P3, N3, 0, 0}, 	/* 1 */
 419	{N3, P3, 0, 0}, 	/* 2 */
 420	{P3, P3, 0, 0}, 	/* 3 */
 421	{0, 0, 0, 0}		/* backstop */
 422};
 423
 424/*
 425 * Digits 0-2 (hour tens)
 426 */
 427#define P2	(P15 / 2)	/* mark (+1) */
 428#define N2	(N15 / 2)	/* space (-1) */
 429
 430double bcd2[][4] = {
 431	{N2, N2, 0, 0}, 	/* 0 */
 432	{P2, N2, 0, 0}, 	/* 1 */
 433	{N2, P2, 0, 0}, 	/* 2 */
 434	{0, 0, 0, 0}		/* backstop */
 435};
 436
 437/*
 438 * DST decode (DST2 DST1) for prettyprint
 439 */
 440char dstcod[] = {
 441	'S',			/* 00 standard time */
 442	'I',			/* 01 set clock ahead at 0200 local */
 443	'O',			/* 10 set clock back at 0200 local */
 444	'D'			/* 11 daylight time */
 445};
 446
 447/*
 448 * The decoding matrix consists of nine row vectors, one for each digit
 449 * of the timecode. The digits are stored from least to most significant
 450 * order. The maximum likelihood timecode is formed from the digits
 451 * corresponding to the maximum likelihood values reading in the
 452 * opposite order: yy ddd hh:mm.
 453 */
 454struct decvec {
 455	int radix;		/* radix (3, 4, 6, 10) */
 456	int digit;		/* current clock digit */
 457	int mldigit;		/* maximum likelihood digit */
 458	int count;		/* match count */
 459	double digprb;		/* max digit probability */
 460	double digsnr;		/* likelihood function (dB) */
 461	double like[10];	/* likelihood integrator 0-9 */
 462};
 463
 464/*
 465 * The station structure (sp) is used to acquire the minute pulse from
 466 * WWV and/or WWVH. These stations are distinguished by the frequency
 467 * used for the second and minute sync pulses, 1000 Hz for WWV and 1200
 468 * Hz for WWVH. Other than frequency, the format is the same.
 469 */
 470struct sync {
 471	double	epoch;		/* accumulated epoch differences */
 472	double	maxeng;		/* sync max energy */
 473	double	noieng;		/* sync noise energy */
 474	long	pos;		/* max amplitude position */
 475	long	lastpos;	/* last max position */
 476	long	mepoch;		/* minute synch epoch */
 477
 478	double	amp;		/* sync signal */
 479	double	syneng;		/* sync signal max */
 480	double	synmax;		/* sync signal max latched at 0 s */
 481	double	synsnr;		/* sync signal SNR */
 482	double	metric;		/* signal quality metric */
 483	int	reach;		/* reachability register */
 484	int	count;		/* bit counter */
 485	int	select;		/* select bits */
 486	char	refid[5];	/* reference identifier */
 487};
 488
 489/*
 490 * The channel structure (cp) is used to mitigate between channels.
 491 */
 492struct chan {
 493	int	gain;		/* audio gain */
 494	struct sync wwv;	/* wwv station */
 495	struct sync wwvh;	/* wwvh station */
 496};
 497
 498/*
 499 * WWV unit control structure (up)
 500 */
 501struct wwvunit {
 502	l_fp	timestamp;	/* audio sample timestamp */
 503	l_fp	tick;		/* audio sample increment */
 504	double	phase, freq;	/* logical clock phase and frequency */
 505	double	monitor;	/* audio monitor point */
 506#ifdef ICOM
 507	int	fd_icom;	/* ICOM file descriptor */
 508#endif /* ICOM */
 509	int	errflg;		/* error flags */
 510	int	watch;		/* watchcat */
 511
 512	/*
 513	 * Audio codec variables
 514	 */
 515	double	comp[SIZE];	/* decompanding table */
 516	int	port;		/* codec port */
 517	int	gain;		/* codec gain */
 518	int	mongain;	/* codec monitor gain */
 519	int	clipcnt;	/* sample clipped count */
 520
 521	/*
 522	 * Variables used to establish basic system timing
 523	 */
 524	int	avgint;		/* master time constant */
 525	int	yepoch;		/* sync epoch */
 526	int	repoch;		/* buffered sync epoch */
 527	double	epomax;		/* second sync amplitude */
 528	double	eposnr;		/* second sync SNR */
 529	double	irig;		/* data I channel amplitude */
 530	double	qrig;		/* data Q channel amplitude */
 531	int	datapt;		/* 100 Hz ramp */
 532	double	datpha;		/* 100 Hz VFO control */
 533	int	rphase;		/* second sample counter */
 534	long	mphase;		/* minute sample counter */
 535
 536	/*
 537	 * Variables used to mitigate which channel to use
 538	 */
 539	struct chan mitig[NCHAN]; /* channel data */
 540	struct sync *sptr;	/* station pointer */
 541	int	dchan;		/* data channel */
 542	int	schan;		/* probe channel */
 543	int	achan;		/* active channel */
 544
 545	/*
 546	 * Variables used by the clock state machine
 547	 */
 548	struct decvec decvec[9]; /* decoding matrix */
 549	int	rsec;		/* seconds counter */
 550	int	digcnt;		/* count of digits synchronized */
 551
 552	/*
 553	 * Variables used to estimate signal levels and bit/digit
 554	 * probabilities
 555	 */
 556	double	datsig;		/* data signal max */
 557	double	datsnr;		/* data signal SNR (dB) */
 558
 559	/*
 560	 * Variables used to establish status and alarm conditions
 561	 */
 562	int	status;		/* status bits */
 563	int	alarm;		/* alarm flashers */
 564	int	misc;		/* miscellaneous timecode bits */
 565	int	errcnt;		/* data bit error counter */
 566};
 567
 568/*
 569 * Function prototypes
 570 */
 571static	int	wwv_start	P((int, struct peer *));
 572static	void	wwv_shutdown	P((int, struct peer *));
 573static	void	wwv_receive	P((struct recvbuf *));
 574static	void	wwv_poll	P((int, struct peer *));
 575
 576/*
 577 * More function prototypes
 578 */
 579static	void	wwv_epoch	P((struct peer *));
 580static	void	wwv_rf		P((struct peer *, double));
 581static	void	wwv_endpoc	P((struct peer *, int));
 582static	void	wwv_rsec	P((struct peer *, double));
 583static	void	wwv_qrz		P((struct peer *, struct sync *, int));
 584static	void	wwv_corr4	P((struct peer *, struct decvec *,
 585				    double [], double [][4]));
 586static	void	wwv_gain	P((struct peer *));
 587static	void	wwv_tsec	P((struct peer *));
 588static	int	timecode	P((struct wwvunit *, char *));
 589static	double	wwv_snr		P((double, double));
 590static	int	carry		P((struct decvec *));
 591static	int	wwv_newchan	P((struct peer *));
 592static	void	wwv_newgame	P((struct peer *));
 593static	double	wwv_metric	P((struct sync *));
 594static	void	wwv_clock	P((struct peer *));
 595#ifdef ICOM
 596static	int	wwv_qsy		P((struct peer *, int));
 597#endif /* ICOM */
 598
 599static double qsy[NCHAN] = {2.5, 5, 10, 15, 20}; /* frequencies (MHz) */
 600
 601/*
 602 * Transfer vector
 603 */
 604struct	refclock refclock_wwv = {
 605	wwv_start,		/* start up driver */
 606	wwv_shutdown,		/* shut down driver */
 607	wwv_poll,		/* transmit poll message */
 608	noentry,		/* not used (old wwv_control) */
 609	noentry,		/* initialize driver (not used) */
 610	noentry,		/* not used (old wwv_buginfo) */
 611	NOFLAGS			/* not used */
 612};
 613
 614
 615/*
 616 * wwv_start - open the devices and initialize data for processing
 617 */
 618static int
 619wwv_start(
 620	int	unit,		/* instance number (used by PCM) */
 621	struct peer *peer	/* peer structure pointer */
 622	)
 623{
 624	struct refclockproc *pp;
 625	struct wwvunit *up;
 626#ifdef ICOM
 627	int	temp;
 628#endif /* ICOM */
 629
 630	/*
 631	 * Local variables
 632	 */
 633	int	fd;		/* file descriptor */
 634	int	i;		/* index */
 635	double	step;		/* codec adjustment */
 636
 637	/*
 638	 * Open audio device
 639	 */
 640	fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit);
 641	if (fd < 0)
 642		return (0);
 643#ifdef DEBUG
 644	if (debug)
 645		audio_show();
 646#endif /* DEBUG */
 647
 648	/*
 649	 * Allocate and initialize unit structure
 650	 */
 651	if (!(up = (struct wwvunit *)emalloc(sizeof(struct wwvunit)))) {
 652		close(fd);
 653		return (0);
 654	}
 655	memset(up, 0, sizeof(struct wwvunit));
 656	pp = peer->procptr;
 657	pp->unitptr = (caddr_t)up;
 658	pp->io.clock_recv = wwv_receive;
 659	pp->io.srcclock = (caddr_t)peer;
 660	pp->io.datalen = 0;
 661	pp->io.fd = fd;
 662	if (!io_addclock(&pp->io)) {
 663		close(fd);
 664		free(up);
 665		return (0);
 666	}
 667
 668	/*
 669	 * Initialize miscellaneous variables
 670	 */
 671	peer->precision = PRECISION;
 672	pp->clockdesc = DESCRIPTION;
 673
 674	/*
 675	 * The companded samples are encoded sign-magnitude. The table
 676	 * contains all the 256 values in the interest of speed.
 677	 */
 678	up->comp[0] = up->comp[OFFSET] = 0.;
 679	up->comp[1] = 1.; up->comp[OFFSET + 1] = -1.;
 680	up->comp[2] = 3.; up->comp[OFFSET + 2] = -3.;
 681	step = 2.;
 682	for (i = 3; i < OFFSET; i++) {
 683		up->comp[i] = up->comp[i - 1] + step;
 684		up->comp[OFFSET + i] = -up->comp[i];
 685                if (i % 16 == 0)
 686		    step *= 2.;
 687	}
 688	DTOLFP(1. / SECOND, &up->tick);
 689
 690	/*
 691	 * Initialize the decoding matrix with the radix for each digit
 692	 * position.
 693	 */
 694	up->decvec[MN].radix = 10;	/* minutes */
 695	up->decvec[MN + 1].radix = 6;
 696	up->decvec[HR].radix = 10;	/* hours */
 697	up->decvec[HR + 1].radix = 3;
 698	up->decvec[DA].radix = 10;	/* days */
 699	up->decvec[DA + 1].radix = 10;
 700	up->decvec[DA + 2].radix = 4;
 701	up->decvec[YR].radix = 10;	/* years */
 702	up->decvec[YR + 1].radix = 10;
 703
 704#ifdef ICOM
 705	/*
 706	 * Initialize autotune if available. Note that the ICOM select
 707	 * code must be less than 128, so the high order bit can be used
 708	 * to select the line speed 0 (9600 bps) or 1 (1200 bps).
 709	 */
 710	temp = 0;
 711#ifdef DEBUG
 712	if (debug > 1)
 713		temp = P_TRACE;
 714#endif /* DEBUG */
 715	if (peer->ttl != 0) {
 716		if (peer->ttl & 0x80)
 717			up->fd_icom = icom_init("/dev/icom", B1200,
 718			    temp);
 719		else
 720			up->fd_icom = icom_init("/dev/icom", B9600,
 721			    temp);
 722		if (up->fd_icom < 0) {
 723			NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
 724			    msyslog(LOG_NOTICE,
 725			    "icom: %m");
 726			up->errflg = CEVNT_FAULT;
 727		}
 728	}
 729	if (up->fd_icom > 0) {
 730		if (wwv_qsy(peer, DCHAN) != 0) {
 731			NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
 732			    msyslog(LOG_NOTICE,
 733			    "icom: radio not found");
 734			up->errflg = CEVNT_FAULT;
 735			close(up->fd_icom);
 736			up->fd_icom = 0;
 737		} else {
 738			NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
 739			    msyslog(LOG_NOTICE,
 740			    "icom: autotune enabled");
 741		}
 742	}
 743#endif /* ICOM */
 744
 745	/*
 746	 * Let the games begin.
 747	 */
 748	wwv_newgame(peer);
 749	return (1);
 750}
 751
 752
 753/*
 754 * wwv_shutdown - shut down the clock
 755 */
 756static void
 757wwv_shutdown(
 758	int	unit,		/* instance number (not used) */
 759	struct peer *peer	/* peer structure pointer */
 760	)
 761{
 762	struct refclockproc *pp;
 763	struct wwvunit *up;
 764
 765	pp = peer->procptr;
 766	up = (struct wwvunit *)pp->unitptr;
 767	if (up == NULL)
 768		return;
 769
 770	io_closeclock(&pp->io);
 771#ifdef ICOM
 772	if (up->fd_icom > 0)
 773		close(up->fd_icom);
 774#endif /* ICOM */
 775	free(up);
 776}
 777
 778
 779/*
 780 * wwv_receive - receive data from the audio device
 781 *
 782 * This routine reads input samples and adjusts the logical clock to
 783 * track the A/D sample clock by dropping or duplicating codec samples.
 784 * It also controls the A/D signal level with an AGC loop to mimimize
 785 * quantization noise and avoid overload.
 786 */
 787static void
 788wwv_receive(
 789	struct recvbuf *rbufp	/* receive buffer structure pointer */
 790	)
 791{
 792	struct peer *peer;
 793	struct refclockproc *pp;
 794	struct wwvunit *up;
 795
 796	/*
 797	 * Local variables
 798	 */
 799	double	sample;		/* codec sample */
 800	u_char	*dpt;		/* buffer pointer */
 801	int	bufcnt;		/* buffer counter */
 802	l_fp	ltemp;
 803
 804	peer = (struct peer *)rbufp->recv_srcclock;
 805	pp = peer->procptr;
 806	up = (struct wwvunit *)pp->unitptr;
 807
 808	/*
 809	 * Main loop - read until there ain't no more. Note codec
 810	 * samples are bit-inverted.
 811	 */
 812	DTOLFP((double)rbufp->recv_length / SECOND, &ltemp);
 813	L_SUB(&rbufp->recv_time, &ltemp);
 814	up->timestamp = rbufp->recv_time;
 815	dpt = rbufp->recv_buffer;
 816	for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
 817		sample = up->comp[~*dpt++ & 0xff];
 818
 819		/*
 820		 * Clip noise spikes greater than MAXAMP (6000) and
 821		 * record the number of clips to be used later by the
 822		 * AGC.
 823		 */
 824		if (sample > MAXAMP) {
 825			sample = MAXAMP;
 826			up->clipcnt++;
 827		} else if (sample < -MAXAMP) {
 828			sample = -MAXAMP;
 829			up->clipcnt++;
 830		}
 831
 832		/*
 833		 * Variable frequency oscillator. The codec oscillator
 834		 * runs at the nominal rate of 8000 samples per second,
 835		 * or 125 us per sample. A frequency change of one unit
 836		 * results in either duplicating or deleting one sample
 837		 * per second, which results in a frequency change of
 838		 * 125 PPM.
 839		 */
 840		up->phase += up->freq / SECOND;
 841		up->phase += FREQ_OFFSET / 1e6;
 842		if (up->phase >= .5) {
 843			up->phase -= 1.;
 844		} else if (up->phase < -.5) {
 845			up->phase += 1.;
 846			wwv_rf(peer, sample);
 847			wwv_rf(peer, sample);
 848		} else {
 849			wwv_rf(peer, sample);
 850		}
 851		L_ADD(&up->timestamp, &up->tick);
 852	}
 853
 854	/*
 855	 * Set the input port and monitor gain for the next buffer.
 856	 */
 857	if (pp->sloppyclockflag & CLK_FLAG2)
 858		up->port = 2;
 859	else
 860		up->port = 1;
 861	if (pp->sloppyclockflag & CLK_FLAG3)
 862		up->mongain = MONGAIN;
 863	else
 864		up->mongain = 0;
 865}
 866
 867
 868/*
 869 * wwv_poll - called by the transmit procedure
 870 *
 871 * This routine keeps track of status. If no offset samples have been
 872 * processed during a poll interval, a timeout event is declared. If
 873 * errors have have occurred during the interval, they are reported as
 874 * well.
 875 */
 876static void
 877wwv_poll(
 878	int	unit,		/* instance number (not used) */
 879	struct peer *peer	/* peer structure pointer */
 880	)
 881{
 882	struct refclockproc *pp;
 883	struct wwvunit *up;
 884
 885	pp = peer->procptr;
 886	up = (struct wwvunit *)pp->unitptr;
 887	if (pp->coderecv == pp->codeproc)
 888		up->errflg = CEVNT_TIMEOUT;
 889	if (up->errflg)
 890		refclock_report(peer, up->errflg);
 891	up->errflg = 0;
 892	pp->polls++;
 893}
 894
 895
 896/*
 897 * wwv_rf - process signals and demodulate to baseband
 898 *
 899 * This routine grooms and filters decompanded raw audio samples. The
 900 * output signal is the 100-Hz filtered baseband data signal in
 901 * quadrature phase. The routine also determines the minute synch epoch,
 902 * as well as certain signal maxima, minima and related values.
 903 *
 904 * There are two 1-s ramps used by this program. Both count the 8000
 905 * logical clock samples spanning exactly one second. The epoch ramp
 906 * counts the samples starting at an arbitrary time. The rphase ramp
 907 * counts the samples starting at the 5-ms second sync pulse found
 908 * during the epoch ramp.
 909 *
 910 * There are two 1-m ramps used by this program. The mphase ramp counts
 911 * the 480,000 logical clock samples spanning exactly one minute and
 912 * starting at an arbitrary time. The rsec ramp counts the 60 seconds of
 913 * the minute starting at the 800-ms minute sync pulse found during the
 914 * mphase ramp. The rsec ramp drives the seconds state machine to
 915 * determine the bits and digits of the timecode. 
 916 *
 917 * Demodulation operations are based on three synthesized quadrature
 918 * sinusoids: 100 Hz for the data signal, 1000 Hz for the WWV sync
 919 * signal and 1200 Hz for the WWVH sync signal. These drive synchronous
 920 * matched filters for the data signal (170 ms at 100 Hz), WWV minute
 921 * sync signal (800 ms at 1000 Hz) and WWVH minute sync signal (800 ms
 922 * at 1200 Hz). Two additional matched filters are switched in
 923 * as required for the WWV second sync signal (5 cycles at 1000 Hz) and
 924 * WWVH second sync signal (6 cycles at 1200 Hz).
 925 */
 926static void
 927wwv_rf(
 928	struct peer *peer,	/* peerstructure pointer */
 929	double isig		/* input signal */
 930	)
 931{
 932	struct refclockproc *pp;
 933	struct wwvunit *up;
 934	struct sync *sp, *rp;
 935
 936	static double lpf[5];	/* 150-Hz lpf delay line */
 937	double data;		/* lpf output */
 938	static double bpf[9];	/* 1000/1200-Hz bpf delay line */
 939	double syncx;		/* bpf output */
 940	static double mf[41];	/* 1000/1200-Hz mf delay line */
 941	double mfsync;		/* mf output */
 942
 943	static int iptr;	/* data channel pointer */
 944	static double ibuf[DATSIZ]; /* data I channel delay line */
 945	static double qbuf[DATSIZ]; /* data Q channel delay line */
 946
 947	static int jptr;	/* sync channel pointer */
 948	static int kptr;	/* tick channel pointer */
 949
 950	static int csinptr;	/* wwv channel phase */
 951	static double cibuf[SYNSIZ]; /* wwv I channel delay line */
 952	static double cqbuf[SYNSIZ]; /* wwv Q channel delay line */
 953	static double ciamp;	/* wwv I channel amplitude */
 954	static double cqamp;	/* wwv Q channel amplitude */
 955
 956	static double csibuf[TCKSIZ]; /* wwv I tick delay line */
 957	static double csqbuf[TCKSIZ]; /* wwv Q tick delay line */
 958	static double csiamp;	/* wwv I tick amplitude */
 959	static double csqamp;	/* wwv Q tick amplitude */
 960
 961	static int hsinptr;	/* wwvh channel phase */
 962	static double hibuf[SYNSIZ]; /* wwvh I channel delay line */
 963	static double hqbuf[SYNSIZ]; /* wwvh Q channel delay line */
 964	static double hiamp;	/* wwvh I channel amplitude */
 965	static double hqamp;	/* wwvh Q channel amplitude */
 966
 967	static double hsibuf[TCKSIZ]; /* wwvh I tick delay line */
 968	static double hsqbuf[TCKSIZ]; /* wwvh Q tick delay line */
 969	static double hsiamp;	/* wwvh I tick amplitude */
 970	static double hsqamp;	/* wwvh Q tick amplitude */
 971
 972	static double epobuf[SECOND]; /* second sync comb filter */
 973	static double epomax, nxtmax; /* second sync amplitude buffer */
 974	static int epopos;	/* epoch second sync position buffer */
 975
 976	static int iniflg;	/* initialization flag */
 977	int	pdelay;		/* propagation delay (samples) */
 978	int	epoch;		/* comb filter index */
 979	double	dtemp;
 980	int	i;
 981
 982	pp = peer->procptr;
 983	up = (struct wwvunit *)pp->unitptr;
 984
 985	if (!iniflg) {
 986		iniflg = 1;
 987		memset((char *)lpf, 0, sizeof(lpf));
 988		memset((char *)bpf, 0, sizeof(bpf));
 989		memset((char *)mf, 0, sizeof(mf));
 990		memset((char *)ibuf, 0, sizeof(ibuf));
 991		memset((char *)qbuf, 0, sizeof(qbuf));
 992		memset((char *)cibuf, 0, sizeof(cibuf));
 993		memset((char *)cqbuf, 0, sizeof(cqbuf));
 994		memset((char *)csibuf, 0, sizeof(csibuf));
 995		memset((char *)csqbuf, 0, sizeof(csqbuf));
 996		memset((char *)hibuf, 0, sizeof(hibuf));
 997		memset((char *)hqbuf, 0, sizeof(hqbuf));
 998		memset((char *)hsibuf, 0, sizeof(hsibuf));
 999		memset((char *)hsqbuf, 0, sizeof(hsqbuf));
1000		memset((char *)epobuf, 0, sizeof(epobuf));
1001	}
1002
1003	/*
1004	 * Baseband data demodulation. The 100-Hz subcarrier is
1005	 * extracted using a 150-Hz IIR lowpass filter. This attenuates
1006	 * the 1000/1200-Hz sync signals, as well as the 440-Hz and
1007	 * 600-Hz tones and most of the noise and voice modulation
1008	 * components.
1009	 *
1010	 * The subcarrier is transmitted 10 dB down from the carrier.
1011	 * The DGAIN parameter can be adjusted for this and to
1012	 * compensate for the radio audio response at 100 Hz.
1013	 *
1014	 * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB
1015	 * passband ripple, -50 dB stopband ripple.
1016	 */
1017	data = (lpf[4] = lpf[3]) * 8.360961e-01;
1018	data += (lpf[3] = lpf[2]) * -3.481740e+00;
1019	data += (lpf[2] = lpf[1]) * 5.452988e+00;
1020	data += (lpf[1] = lpf[0]) * -3.807229e+00;
1021	lpf[0] = isig * DGAIN - data;
1022	data = lpf[0] * 3.281435e-03
1023	    + lpf[1] * -1.149947e-02
1024	    + lpf[2] * 1.654858e-02
1025	    + lpf[3] * -1.149947e-02
1026	    + lpf[4] * 3.281435e-03;
1027
1028	/*
1029	 * The 100-Hz data signal is demodulated using a pair of
1030	 * quadrature multipliers, matched filters and a phase lock
1031	 * loop. The I and Q quadrature data signals are produced by
1032	 * multiplying the filtered signal by 100-Hz sine and cosine
1033	 * signals, respectively. The signals are processed by 170-ms
1034	 * synchronous matched filters to produce the amplitude and
1035	 * phase signals used by the demodulator. The signals are scaled
1036	 * to produce unit energy at the maximum value.
1037	 */
1038	i = up->datapt;
1039	up->datapt = (up->datapt + IN100) % 80;
1040	dtemp = sintab[i] * data / (MS / 2. * DATCYC);
1041	up->irig -= ibuf[iptr];
1042	ibuf[iptr] = dtemp;
1043	up->irig += dtemp;
1044
1045	i = (i + 20) % 80;
1046	dtemp = sintab[i] * data / (MS / 2. * DATCYC);
1047	up->qrig -= qbuf[iptr];
1048	qbuf[iptr] = dtemp;
1049	up->qrig += dtemp;
1050	iptr = (iptr + 1) % DATSIZ;
1051
1052	/*
1053	 * Baseband sync demodulation. The 1000/1200 sync signals are
1054	 * extracted using a 600-Hz IIR bandpass filter. This removes
1055	 * the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz
1056	 * tones and most of the noise and voice modulation components.
1057	 *
1058	 * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB
1059	 * passband ripple, -50 dB stopband ripple.
1060	 */
1061	syncx = (bpf[8] = bpf[7]) * 4.897278e-01;
1062	syncx += (bpf[7] = bpf[6]) * -2.765914e+00;
1063	syncx += (bpf[6] = bpf[5]) * 8.110921e+00;
1064	syncx += (bpf[5] = bpf[4]) * -1.517732e+01;
1065	syncx += (bpf[4] = bpf[3]) * 1.975197e+01;
1066	syncx += (bpf[3] = bpf[2]) * -1.814365e+01;
1067	syncx += (bpf[2] = bpf[1]) * 1.159783e+01;
1068	syncx += (bpf[1] = bpf[0]) * -4.735040e+00;
1069	bpf[0] = isig - syncx;
1070	syncx = bpf[0] * 8.203628e-03
1071	    + bpf[1] * -2.375732e-02
1072	    + bpf[2] * 3.353214e-02
1073	    + bpf[3] * -4.080258e-02
1074	    + bpf[4] * 4.605479e-02
1075	    + bpf[5] * -4.080258e-02
1076	    + bpf[6] * 3.353214e-02
1077	    + bpf[7] * -2.375732e-02
1078	    + bpf[8] * 8.203628e-03;
1079
1080	/*
1081	 * The 1000/1200 sync signals are demodulated using a pair of
1082	 * quadrature multipliers and matched filters. However,
1083	 * synchronous demodulation at these frequencies is impractical,
1084	 * so only the signal amplitude is used. The I and Q quadrature
1085	 * sync signals are produced by multiplying the filtered signal
1086	 * by 1000-Hz (WWV) and 1200-Hz (WWVH) sine and cosine signals,
1087	 * respectively. The WWV and WWVH signals are processed by 800-
1088	 * ms synchronous matched filters and combined to produce the
1089	 * minute sync signal and detect which one (or both) the WWV or
1090	 * WWVH signal is present. The WWV and WWVH signals are also
1091	 * processed by 5-ms synchronous matched filters and combined to
1092	 * produce the second sync signal. The signals are scaled to
1093	 * produce unit energy at the maximum value.
1094	 *
1095	 * Note the master timing ramps, which run continuously. The
1096	 * minute counter (mphase) counts the samples in the minute,
1097	 * while the second counter (epoch) counts the samples in the
1098	 * second.
1099	 */
1100	up->mphase = (up->mphase + 1) % MINUTE;
1101	epoch = up->mphase % SECOND;
1102
1103	/*
1104	 * WWV
1105	 */
1106	i = csinptr;
1107	csinptr = (csinptr + IN1000) % 80;
1108
1109	dtemp = sintab[i] * syncx / (MS / 2.);
1110	ciamp -= cibuf[jptr];
1111	cibuf[jptr] = dtemp;
1112	ciamp += dtemp;
1113	csiamp -= csibuf[kptr];
1114	csibuf[kptr] = dtemp;
1115	csiamp += dtemp;
1116
1117	i = (i + 20) % 80;
1118	dtemp = sintab[i] * syncx / (MS / 2.);
1119	cqamp -= cqbuf[jptr];
1120	cqbuf[jptr] = dtemp;
1121	cqamp += dtemp;
1122	csqamp -= csqbuf[kptr];
1123	csqbuf[kptr] = dtemp;
1124	csqamp += dtemp;
1125
1126	sp = &up->mitig[up->achan].wwv;
1127	sp->amp = sqrt(ciamp * ciamp + cqamp * cqamp) / SYNCYC;
1128	if (!(up->status & MSYNC))
1129		wwv_qrz(peer, sp, (int)(pp->fudgetime1 * SECOND));
1130
1131	/*
1132	 * WWVH
1133	 */
1134	i = hsinptr;
1135	hsinptr = (hsinptr + IN1200) % 80;
1136
1137	dtemp = sintab[i] * syncx / (MS / 2.);
1138	hiamp -= hibuf[jptr];
1139	hibuf[jptr] = dtemp;
1140	hiamp += dtemp;
1141	hsiamp -= hsibuf[kptr];
1142	hsibuf[kptr] = dtemp;
1143	hsiamp += dtemp;
1144
1145	i = (i + 20) % 80;
1146	dtemp = sintab[i] * syncx / (MS / 2.);
1147	hqamp -= hqbuf[jptr];
1148	hqbuf[jptr] = dtemp;
1149	hqamp += dtemp;
1150	hsqamp -= hsqbuf[kptr];
1151	hsqbuf[kptr] = dtemp;
1152	hsqamp += dtemp;
1153
1154	rp = &up->mitig[up->achan].wwvh;
1155	rp->amp = sqrt(hiamp * hiamp + hqamp * hqamp) / SYNCYC;
1156	if (!(up->status & MSYNC))
1157		wwv_qrz(peer, rp, (int)(pp->fudgetime2 * SECOND));
1158	jptr = (jptr + 1) % SYNSIZ;
1159	kptr = (kptr + 1) % TCKSIZ;
1160
1161	/*
1162	 * The following section is called once per minute. It does
1163	 * housekeeping and timeout functions and empties the dustbins.
1164	 */
1165	if (up->mphase == 0) {
1166		up->watch++;
1167		if (!(up->status & MSYNC)) {
1168
1169			/*
1170			 * If minute sync has not been acquired before
1171			 * ACQSN timeout (6 min), or if no signal is
1172			 * heard, the program cycles to the next
1173			 * frequency and tries again.
1174			 */
1175			if (!wwv_newchan(peer))
1176				up->watch = 0;
1177#ifdef ICOM
1178			if (up->fd_icom > 0)
1179				wwv_qsy(peer, up->dchan);
1180#endif /* ICOM */
1181		} else {
1182
1183			/*
1184			 * If the leap bit is set, set the minute epoch
1185			 * back one second so the station processes
1186			 * don't miss a beat.
1187			 */
1188			if (up->status & LEPSEC) {
1189				up->mphase -= SECOND;
1190				if (up->mphase < 0)
1191					up->mphase += MINUTE;
1192			}
1193		}
1194	}
1195
1196	/*
1197	 * When the channel metric reaches threshold and the second
1198	 * counter matches the minute epoch within the second, the
1199	 * driver has synchronized to the station. The second number is
1200	 * the remaining seconds until the next minute epoch, while the
1201	 * sync epoch is zero. Watch out for the first second; if
1202	 * already synchronized to the second, the buffered sync epoch
1203	 * must be set.
1204	 *
1205	 * Note the guard interval is 200 ms; if for some reason the
1206	 * clock drifts more than that, it might wind up in the wrong
1207	 * second. If the maximum frequency error is not more than about
1208	 * 1 PPM, the clock can go as much as two days while still in
1209	 * the same second.
1210	 */
1211	if (up->status & MSYNC) {
1212		wwv_epoch(peer);
1213	} else if (up->sptr != NULL) {
1214		sp = up->sptr;
1215		if (sp->metric >= TTHR && epoch == sp->mepoch % SECOND) 		    {
1216			up->rsec = (60 - sp->mepoch / SECOND) % 60;
1217			up->rphase = 0;
1218			up->status |= MSYNC;
1219			up->watch = 0;
1220			if (!(up->status & SSYNC))
1221				up->repoch = up->yepoch = epoch;
1222			else
1223				up->repoch = up->yepoch;
1224			
1225		}
1226	}
1227
1228	/*
1229	 * The second sync pulse is extracted using 5-ms (40 sample) FIR
1230	 * matched filters at 1000 Hz for WWV or 1200 Hz for WWVH. This
1231	 * pulse is used for the most precise synchronization, since if
1232	 * provides a resolution of one sample (125 us). The filters run
1233	 * only if the station has been reliably determined.
1234	 */
1235	if (up->status & SELV) {
1236		pdelay = (int)(pp->fudgetime1 * SECOND);
1237		mfsync = sqrt(csiamp * csiamp + csqamp * csqamp) /
1238		    TCKCYC;
1239	} else if (up->status & SELH) {
1240		pdelay = (int)(pp->fudgetime2 * SECOND);
1241		mfsync = sqrt(hsiamp * hsiamp + hsqamp * hsqamp) /
1242		    TCKCYC;
1243	} else {
1244		pdelay = 0;
1245		mfsync = 0;
1246	}
1247
1248	/*
1249	 * Enhance the seconds sync pulse using a 1-s (8000-sample) comb
1250	 * filter. Correct for the FIR matched filter delay, which is 5
1251	 * ms for both the WWV and WWVH filters, and also for the
1252	 * propagation delay. Once each second look for second sync. If
1253	 * not in minute sync, fiddle the codec gain. Note the SNR is
1254	 * computed from the maximum sample and the envelope of the
1255	 * sample 6 ms before it, so if we slip more than a cycle the
1256	 * SNR should plummet. The signal is scaled to produce unit
1257	 * energy at the maximum value.
1258	 */
1259	dtemp = (epobuf[epoch] += (mfsync - epobuf[epoch]) /
1260	    up->avgint);
1261	if (dtemp > epomax) {
1262		int	j;
1263
1264		epomax = dtemp;
1265		epopos = epoch;
1266		j = epoch - 6 * MS;
1267		if (j < 0)
1268			j += SECOND;
1269		nxtmax = fabs(epobuf[j]);
1270	}
1271	if (epoch == 0) {
1272		up->epomax = epomax;
1273		up->eposnr = wwv_snr(epomax, nxtmax);
1274		epopos -= pdelay + TCKCYC * MS;
1275		if (epopos < 0)
1276			epopos += SECOND;
1277		wwv_endpoc(peer, epopos);
1278		if (!(up->status & SSYNC))
1279			up->alarm |= SYNERR;
1280		epomax = 0;
1281		if (!(up->status & MSYNC))
1282			wwv_gain(peer);
1283	}
1284}
1285
1286
1287/*
1288 * wwv_qrz - identify and acquire WWV/WWVH minute sync pulse
1289 *
1290 * This routine implements a virtual station process used to acquire
1291 * minute sync and to mitigate among the ten frequency and station
1292 * combinations. During minute sync acquisition the process probes each
1293 * frequency and station in turn for the minute pulse, which
1294 * involves searching through the entire 480,000-sample minute. The
1295 * process finds the maximum signal and RMS noise plus signal. Then, the
1296 * actual noise is determined by subtracting the energy of the matched
1297 * filter.
1298 *
1299 * Students of radar receiver technology will discover this algorithm
1300 * amounts to a range-gate discriminator. A valid pulse must have peak
1301 * amplitude at least QTHR (2500) and SNR at least QSNR (20) dB and the
1302 * difference between the current and previous epoch must be less than
1303 * AWND (20 ms). Note that the discriminator peak occurs about 800 ms
1304 * into the second, so the timing is retarded to the previous second
1305 * epoch.
1306 */
1307static void
1308wwv_qrz(
1309	struct peer *peer,	/* peer structure pointer */
1310	struct sync *sp,	/* sync channel structure */
1311	int	pdelay		/* propagation delay (samples) */
1312	)
1313{
1314	struct refclockproc *pp;
1315	struct wwvunit *up;
1316	char	tbuf[80];	/* monitor buffer */
1317	long	epoch;
1318
1319	pp = peer->procptr;
1320	up = (struct wwvunit *)pp->unitptr;
1321
1322	/*
1323	 * Find the sample with peak amplitude, which defines the minute
1324	 * epoch. Accumulate all samples to determine the total noise
1325	 * energy.
1326	 */
1327	epoch = up->mphase - pdelay - SYNSIZ;
1328	if (epoch < 0)
1329		epoch += MINUTE;
1330	if (sp->amp > sp->maxeng) {
1331		sp->maxeng = sp->amp;
1332		sp->pos = epoch;
1333	}
1334	sp->noieng += sp->amp;
1335
1336	/*
1337	 * At the end of the minute, determine the epoch of the minute
1338	 * sync pulse, as well as the difference between the current and
1339	 * previous epoches due to the intrinsic frequency error plus
1340	 * jitter. When calculating the SNR, subtract the pulse energy
1341	 * from the total noise energy and then normalize.
1342	 */
1343	if (up->mphase == 0) {
1344		sp->synmax = sp->maxeng;
1345		sp->synsnr = wwv_snr(sp->synmax, (sp->noieng -
1346		    sp->synmax) / MINUTE);
1347		if (sp->count == 0)
1348			sp->lastpos = sp->pos;
1349		epoch = (sp->pos - sp->lastpos) % MINUTE;
1350		sp->reach <<= 1;
1351		if (sp->reach & (1 << AMAX))
1352			sp->count--;
1353		if (sp->synmax > ATHR && sp->synsnr > ASNR) {
1354			if (abs(epoch) < AWND * MS) {
1355				sp->reach |= 1;
1356				sp->count++;
1357				sp->mepoch = sp->lastpos = sp->pos;
1358			} else if (sp->count == 1) {
1359				sp->lastpos = sp->pos;
1360			}
1361		}
1362		if (up->watch > ACQSN)
1363			sp->metric = 0;
1364		else
1365			sp->metric = wwv_metric(sp);
1366		if (pp->sloppyclockflag & CLK_FLAG4) {
1367			sprintf(tbuf,
1368			    "wwv8 %04x %3d %s %04x %.0f %.0f/%.1f %4ld %4ld",
1369			    up->status, up->gain, sp->refid,
1370			    sp->reach & 0xffff, sp->metric, sp->synmax,
1371			    sp->synsnr, sp->pos % SECOND, epoch);
1372			record_clock_stats(&peer->srcadr, tbuf);
1373#ifdef DEBUG
1374			if (debug)
1375				printf("%s\n", tbuf);
1376#endif /* DEBUG */
1377		}
1378		sp->maxeng = sp->noieng = 0;
1379	}
1380}
1381
1382
1383/*
1384 * wwv_endpoc - identify and acquire second sync pulse
1385 *
1386 * This routine is called at the end of the second sync interval. It
1387 * determines the second sync epoch position within the second and
1388 * disciplines the sample clock using a frequency-lock loop (FLL).
1389 *
1390 * Second sync is determined in the RF input routine as the maximum
1391 * over all 8000 samples in the second comb filter. To assure accurate
1392 * and reliable time and frequency discipline, this routine performs a
1393 * great deal of heavy-handed heuristic data filtering and grooming.
1394 */
1395static void
1396wwv_endpoc(
1397	struct peer *peer,	/* peer structure pointer */
1398	int epopos		/* epoch max position */
1399	)
1400{
1401	struct refclockproc *pp;
1402	struct wwvunit *up;
1403	static int epoch_mf[3]; /* epoch median filter */
1404	static int tepoch;	/* current second epoch */
1405 	static int xepoch;	/* last second epoch */
1406 	static int zepoch;	/* last run epoch */
1407	static int zcount;	/* last run end time */
1408	static int scount;	/* seconds counter */
1409	static int syncnt;	/* run length counter */
1410	static int maxrun;	/* longest run length */
1411	static int mepoch;	/* longest run end epoch */
1412	static int mcount;	/* longest run end time */
1413	static int avgcnt;	/* averaging interval counter */
1414	static int avginc;	/* averaging ratchet */
1415	static int iniflg;	/* initialization flag */
1416	char tbuf[80];		/* monitor buffer */
1417	double dtemp;
1418	int tmp2;
1419
1420	pp = peer->procptr;
1421	up = (struct wwvunit *)pp->unitptr;
1422	if (!iniflg) {
1423		iniflg = 1;
1424		memset((char *)epoch_mf, 0, sizeof(epoch_mf));
1425	}
1426
1427	/*
1428	 * If the signal amplitude or SNR fall below thresholds, dim the
1429	 * second sync lamp and wait for hotter ions. If no stations are
1430	 * heard, we are either in a probe cycle or the ions are really
1431	 * cold. 
1432	 */
1433	scount++;
1434	if (up->epomax < STHR || up->eposnr < SSNR) {
1435		up->status &= ~(SSYNC | FGATE);
1436		avgcnt = syncnt = maxrun = 0;
1437		return;
1438	}
1439	if (!(up->status & (SELV | SELH)))
1440		return;
1441
1442	/*
1443	 * A three-stage median filter is used to help denoise the
1444	 * second sync pulse. The median sample becomes the candidate
1445	 * epoch.
1446	 */
1447	epoch_mf[2] = epoch_mf[1];
1448	epoch_mf[1] = epoch_mf[0];
1449	epoch_mf[0] = epopos;
1450	if (epoch_mf[0] > epoch_mf[1]) {
1451		if (epoch_mf[1] > epoch_mf[2])
1452			tepoch = epoch_mf[1];	/* 0 1 2 */
1453		else if (epoch_mf[2] > epoch_mf[0])
1454			tepoch = epoch_mf[0];	/* 2 0 1 */
1455		else
1456			tepoch = epoch_mf[2];	/* 0 2 1 */
1457	} else {
1458		if (epoch_mf[1] < epoch_mf[2])
1459			tepoch = epoch_mf[1];	/* 2 1 0 */
1460		else if (epoch_mf[2] < epoch_mf[0])
1461			tepoch = epoch_mf[0];	/* 1 0 2 */
1462		else
1463			tepoch = epoch_mf[2];	/* 1 2 0 */
1464	}
1465
1466
1467	/*
1468	 * If the epoch candidate is the same as the last one, increment
1469	 * the run counter. If not, save the length, epoch and end
1470	 * time of the current run for use later and reset the counter.
1471	 * The epoch is considered valid if the run is at least SCMP
1472	 * (10) s, the minute is synchronized and the interval since the
1473	 * last epoch  is not greater than the averaging interval. Thus,
1474	 * after a long absence, the program will wait a full averaging
1475	 * interval while the comb filter charges up and noise
1476	 * dissapates..
1477	 */
1478	tmp2 = (tepoch - xepoch) % SECOND;
1479	if (tmp2 == 0) {
1480		syncnt++;
1481		if (syncnt > SCMP && up->status & MSYNC && (up->status &
1482		    FGATE || scount - zcount <= up->avgint)) {
1483			up->status |= SSYNC;
1484			up->yepoch = tepoch;
1485		}
1486	} else if (syncnt >= maxrun) {
1487		maxrun = syncnt;
1488		mcount = scount;
1489		mepoch = xepoch;
1490		syncnt = 0;
1491	}
1492	if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status & MSYNC))
1493	    {
1494		sprintf(tbuf,
1495		    "wwv1 %04x %3d %4d %5.0f %5.1f %5d %4d %4d %4d",
1496		    up->status, up->gain, tepoch, up->epomax,
1497		    up->eposnr, tmp2, avgcnt, syncnt,
1498		    maxrun);
1499		record_clock_stats(&peer->srcadr, tbuf);
1500#ifdef DEBUG
1501		if (debug)
1502			printf("%s\n", tbuf);
1503#endif /* DEBUG */
1504	}
1505	avgcnt++;
1506	if (avgcnt < up->avgint) {
1507		xepoch = tepoch;
1508		return;
1509	}
1510
1511	/*
1512	 * The sample clock frequency is disciplined using a first-order
1513	 * feedback loop with time constant consistent with the Allan
1514	 * intercept of typical computer clocks. During each averaging
1515	 * interval the candidate epoch at the end of the longest run is
1516	 * determined. If the longest run is zero, all epoches in the
1517	 * interval are different, so the candidate epoch is the current
1518	 * epoch. The frequency update is computed from the candidate
1519	 * epoch difference (125-us units) and time difference (seconds)
1520	 * between updates.
1521	 */
1522	if (syncnt >= maxrun) {
1523		maxrun = syncnt;
1524		mcount = scount;
1525		mepoch = xepoch;
1526	}
1527	xepoch = tepoch;
1528	if (maxrun == 0) {
1529		mepoch = tepoch;
1530		mcount = scount;
1531	}
1532
1533	/*
1534	 * The master clock runs at the codec sample frequency of 8000
1535	 * Hz, so the intrinsic time resolution is 125 us. The frequency
1536	 * resolution ranges from 18 PPM at the minimum averaging
1537	 * interval of 8 s to 0.12 PPM at the maximum interval of 1024
1538	 * s. An offset update is determined at the end of the longest
1539	 * run in each averaging interval. The frequency adjustment is
1540	 * computed from the difference between offset updates and the
1541	 * interval between them.
1542	 *
1543	 * The maximum frequency adjustment ranges from 187 PPM at the
1544	 * minimum interval to 1.5 PPM at the maximum. If the adjustment
1545	 * exceeds the maximum, the update is discarded and the
1546	 * hysteresis counter is decremented. Otherwise, the frequency
1547	 * is incremented by the adjustment, but clamped to the maximum
1548	 * 187.5 PPM. If the update is less than half the maximum, the
1549	 * hysteresis counter is incremented. If the counter increments
1550	 * to +3, the averaging interval is doubled and the counter set
1551	 * to zero; if it decrements to -3, the interval is halved and
1552	 * the counter set to zero.
1553	 */
1554	dtemp = (mepoch - zepoch) % SECOND;
1555	if (up->status & FGATE) {
1556		if (abs(dtemp) < MAXFREQ * MINAVG) {
1557			up->freq += (dtemp / 2.) / ((mcount - zcount) *
1558			    FCONST);
1559			if (up->freq > MAXFREQ)
1560				up->freq = MAXFREQ;
1561			else if (up->freq < -MAXFREQ)
1562				up->freq = -MAXFREQ;
1563			if (abs(dtemp) < MAXFREQ * MINAVG / 2.) {
1564				if (avginc < 3) {
1565					avginc++;
1566				} else {
1567					if (up->avgint < MAXAVG) {
1568						up->avgint <<= 1;
1569						avginc = 0;
1570					}
1571				}
1572			}
1573		} else {
1574			if (avginc > -3) {
1575				avginc--;
1576			} else {
1577				if (up->avgint > MINAVG) {
1578					up->avgint >>= 1;
1579					avginc = 0;
1580				}
1581			}
1582		}
1583	}
1584	if (pp->sloppyclockflag & CLK_FLAG4) {
1585		sprintf(tbuf,
1586		    "wwv2 %04x %5.0f %5.1f %5d %4d %4d %4d %4.0f %7.2f",
1587		    up->status, up->epomax, up->eposnr, mepoch,
1588		    up->avgint, maxrun, mcount - zcount, dtemp,
1589		    up->freq * 1e6 / SECOND);
1590		r

Large files files are truncated, but you can click here to view the full file