xref: /netbsd-src/external/bsd/ntp/dist/ntpd/refclock_irig.c (revision b1c86f5f087524e68db12794ee9c3e3da1ab17a0)
1 /*	$NetBSD: refclock_irig.c,v 1.1.1.1 2009/12/13 16:55:50 kardel Exp $	*/
2 
3 /*
4  * refclock_irig - audio IRIG-B/E demodulator/decoder
5  */
6 #ifdef HAVE_CONFIG_H
7 #include <config.h>
8 #endif
9 
10 #if defined(REFCLOCK) && defined(CLOCK_IRIG)
11 
12 #include "ntpd.h"
13 #include "ntp_io.h"
14 #include "ntp_refclock.h"
15 #include "ntp_calendar.h"
16 #include "ntp_stdlib.h"
17 
18 #include <stdio.h>
19 #include <ctype.h>
20 #include <math.h>
21 #ifdef HAVE_SYS_IOCTL_H
22 #include <sys/ioctl.h>
23 #endif /* HAVE_SYS_IOCTL_H */
24 
25 #include "audio.h"
26 
27 /*
28  * Audio IRIG-B/E demodulator/decoder
29  *
30  * This driver synchronizes the computer time using data encoded in
31  * IRIG-B/E signals commonly produced by GPS receivers and other timing
32  * devices. The IRIG signal is an amplitude-modulated carrier with
33  * pulse-width modulated data bits. For IRIG-B, the carrier frequency is
34  * 1000 Hz and bit rate 100 b/s; for IRIG-E, the carrier frequenchy is
35  * 100 Hz and bit rate 10 b/s. The driver automatically recognizes which
36  & format is in use.
37  *
38  * The driver requires an audio codec or sound card with sampling rate 8
39  * kHz and mu-law companding. This is the same standard as used by the
40  * telephone industry and is supported by most hardware and operating
41  * systems, including Solaris, SunOS, FreeBSD, NetBSD and Linux. In this
42  * implementation, only one audio driver and codec can be supported on a
43  * single machine.
44  *
45  * The program processes 8000-Hz mu-law companded samples using separate
46  * signal filters for IRIG-B and IRIG-E, a comb filter, envelope
47  * detector and automatic threshold corrector. Cycle crossings relative
48  * to the corrected slice level determine the width of each pulse and
49  * its value - zero, one or position identifier.
50  *
51  * The data encode 20 BCD digits which determine the second, minute,
52  * hour and day of the year and sometimes the year and synchronization
53  * condition. The comb filter exponentially averages the corresponding
54  * samples of successive baud intervals in order to reliably identify
55  * the reference carrier cycle. A type-II phase-lock loop (PLL) performs
56  * additional integration and interpolation to accurately determine the
57  * zero crossing of that cycle, which determines the reference
58  * timestamp. A pulse-width discriminator demodulates the data pulses,
59  * which are then encoded as the BCD digits of the timecode.
60  *
61  * The timecode and reference timestamp are updated once each second
62  * with IRIG-B (ten seconds with IRIG-E) and local clock offset samples
63  * saved for later processing. At poll intervals of 64 s, the saved
64  * samples are processed by a trimmed-mean filter and used to update the
65  * system clock.
66  *
67  * An automatic gain control feature provides protection against
68  * overdriven or underdriven input signal amplitudes. It is designed to
69  * maintain adequate demodulator signal amplitude while avoiding
70  * occasional noise spikes. In order to assure reliable capture, the
71  * decompanded input signal amplitude must be greater than 100 units and
72  * the codec sample frequency error less than 250 PPM (.025 percent).
73  *
74  * Monitor Data
75  *
76  * The timecode format used for debugging and data recording includes
77  * data helpful in diagnosing problems with the IRIG signal and codec
78  * connections. The driver produces one line for each timecode in the
79  * following format:
80  *
81  * 00 00 98 23 19:26:52 2782 143 0.694 10 0.3 66.5 3094572411.00027
82  *
83  * If clockstats is enabled, the most recent line is written to the
84  * clockstats file every 64 s. If verbose recording is enabled (fudge
85  * flag 4) each line is written as generated.
86  *
87  * The first field containes the error flags in hex, where the hex bits
88  * are interpreted as below. This is followed by the year of century,
89  * day of year and time of day. Note that the time of day is for the
90  * previous minute, not the current time. The status indicator and year
91  * are not produced by some IRIG devices and appear as zeros. Following
92  * these fields are the carrier amplitude (0-3000), codec gain (0-255),
93  * modulation index (0-1), time constant (4-10), carrier phase error
94  * +-.5) and carrier frequency error (PPM). The last field is the on-
95  * time timestamp in NTP format.
96  *
97  * The error flags are defined as follows in hex:
98  *
99  * x01	Low signal. The carrier amplitude is less than 100 units. This
100  *	is usually the result of no signal or wrong input port.
101  * x02	Frequency error. The codec frequency error is greater than 250
102  *	PPM. This may be due to wrong signal format or (rarely)
103  *	defective codec.
104  * x04	Modulation error. The IRIG modulation index is less than 0.5.
105  *	This is usually the result of an overdriven codec, wrong signal
106  *	format or wrong input port.
107  * x08	Frame synch error. The decoder frame does not match the IRIG
108  *	frame. This is usually the result of an overdriven codec, wrong
109  *	signal format or noisy IRIG signal. It may also be the result of
110  *	an IRIG signature check which indicates a failure of the IRIG
111  *	signal synchronization source.
112  * x10	Data bit error. The data bit length is out of tolerance. This is
113  *	usually the result of an overdriven codec, wrong signal format
114  *	or noisy IRIG signal.
115  * x20	Seconds numbering discrepancy. The decoder second does not match
116  *	the IRIG second. This is usually the result of an overdriven
117  *	codec, wrong signal format or noisy IRIG signal.
118  * x40	Codec error (overrun). The machine is not fast enough to keep up
119  *	with the codec.
120  * x80	Device status error (Spectracom).
121  *
122  *
123  * Once upon a time, an UltrSPARC 30 and Solaris 2.7 kept the clock
124  * within a few tens of microseconds relative to the IRIG-B signal.
125  * Accuracy with IRIG-E was about ten times worse. Unfortunately, Sun
126  * broke the 2.7 audio driver in 2.8, which has a 10-ms sawtooth
127  * modulation.
128  *
129  * Unlike other drivers, which can have multiple instantiations, this
130  * one supports only one. It does not seem likely that more than one
131  * audio codec would be useful in a single machine. More than one would
132  * probably chew up too much CPU time anyway.
133  *
134  * Fudge factors
135  *
136  * Fudge flag4 causes the dubugging output described above to be
137  * recorded in the clockstats file. Fudge flag2 selects the audio input
138  * port, where 0 is the mike port (default) and 1 is the line-in port.
139  * It does not seem useful to select the compact disc player port. Fudge
140  * flag3 enables audio monitoring of the input signal. For this purpose,
141  * the monitor gain is set t a default value. Fudgetime2 is used as a
142  * frequency vernier for broken codec sample frequency.
143  *
144  * Alarm codes
145  *
146  * CEVNT_BADTIME	invalid date or time
147  * CEVNT_TIMEOUT	no IRIG data since last poll
148  */
149 /*
150  * Interface definitions
151  */
152 #define	DEVICE_AUDIO	"/dev/audio" /* audio device name */
153 #define	PRECISION	(-17)	/* precision assumed (about 10 us) */
154 #define	REFID		"IRIG"	/* reference ID */
155 #define	DESCRIPTION	"Generic IRIG Audio Driver" /* WRU */
156 #define	AUDIO_BUFSIZ	320	/* audio buffer size (40 ms) */
157 #define SECOND		8000	/* nominal sample rate (Hz) */
158 #define BAUD		80	/* samples per baud interval */
159 #define OFFSET		128	/* companded sample offset */
160 #define SIZE		256	/* decompanding table size */
161 #define CYCLE		8	/* samples per bit */
162 #define SUBFLD		10	/* bits per frame */
163 #define FIELD		100	/* bits per second */
164 #define MINTC		2	/* min PLL time constant */
165 #define MAXTC		10	/* max PLL time constant max */
166 #define	MAXAMP		3000.	/* maximum signal amplitude */
167 #define	MINAMP		2000.	/* minimum signal amplitude */
168 #define DRPOUT		100.	/* dropout signal amplitude */
169 #define MODMIN		0.5	/* minimum modulation index */
170 #define MAXFREQ		(250e-6 * SECOND) /* freq tolerance (.025%) */
171 
172 /*
173  * The on-time synchronization point is the positive-going zero crossing
174  * of the first cycle of the second. The IIR baseband filter phase delay
175  * is 1.03 ms for IRIG-B and 3.47 ms for IRIG-E. The fudge value 2.68 ms
176  * due to the codec and other causes was determined by calibrating to a
177  * PPS signal from a GPS receiver.
178  *
179  * The results with a 2.4-GHz P4 running FreeBSD 6.1 are generally
180  * within .02 ms short-term with .02 ms jitter. The processor load due
181  * to the driver is 0.51 percent.
182  */
183 #define IRIG_B	((1.03 + 2.68) / 1000)	/* IRIG-B system delay (s) */
184 #define IRIG_E	((3.47 + 2.68) / 1000)	/* IRIG-E system delay (s) */
185 
186 /*
187  * Data bit definitions
188  */
189 #define BIT0		0	/* zero */
190 #define BIT1		1	/* one */
191 #define BITP		2	/* position identifier */
192 
193 /*
194  * Error flags
195  */
196 #define IRIG_ERR_AMP	0x01	/* low carrier amplitude */
197 #define IRIG_ERR_FREQ	0x02	/* frequency tolerance exceeded */
198 #define IRIG_ERR_MOD	0x04	/* low modulation index */
199 #define IRIG_ERR_SYNCH	0x08	/* frame synch error */
200 #define IRIG_ERR_DECODE	0x10	/* frame decoding error */
201 #define IRIG_ERR_CHECK	0x20	/* second numbering discrepancy */
202 #define IRIG_ERR_ERROR	0x40	/* codec error (overrun) */
203 #define IRIG_ERR_SIGERR	0x80	/* IRIG status error (Spectracom) */
204 
205 static	char	hexchar[] = "0123456789abcdef";
206 
207 /*
208  * IRIG unit control structure
209  */
210 struct irigunit {
211 	u_char	timecode[2 * SUBFLD + 1]; /* timecode string */
212 	l_fp	timestamp;	/* audio sample timestamp */
213 	l_fp	tick;		/* audio sample increment */
214 	l_fp	refstamp;	/* reference timestamp */
215 	l_fp	chrstamp;	/* baud timestamp */
216 	l_fp	prvstamp;	/* previous baud timestamp */
217 	double	integ[BAUD];	/* baud integrator */
218 	double	phase, freq;	/* logical clock phase and frequency */
219 	double	zxing;		/* phase detector integrator */
220 	double	yxing;		/* cycle phase */
221 	double	exing;		/* envelope phase */
222 	double	modndx;		/* modulation index */
223 	double	irig_b;		/* IRIG-B signal amplitude */
224 	double	irig_e;		/* IRIG-E signal amplitude */
225 	int	errflg;		/* error flags */
226 	/*
227 	 * Audio codec variables
228 	 */
229 	double	comp[SIZE];	/* decompanding table */
230 	double	signal;		/* peak signal for AGC */
231 	int	port;		/* codec port */
232 	int	gain;		/* codec gain */
233 	int	mongain;	/* codec monitor gain */
234 	int	seccnt;		/* second interval counter */
235 
236 	/*
237 	 * RF variables
238 	 */
239 	double	bpf[9];		/* IRIG-B filter shift register */
240 	double	lpf[5];		/* IRIG-E filter shift register */
241 	double	envmin, envmax;	/* envelope min and max */
242 	double	slice;		/* envelope slice level */
243 	double	intmin, intmax;	/* integrated envelope min and max */
244 	double	maxsignal;	/* integrated peak amplitude */
245 	double	noise;		/* integrated noise amplitude */
246 	double	lastenv[CYCLE];	/* last cycle amplitudes */
247 	double	lastint[CYCLE];	/* last integrated cycle amplitudes */
248 	double	lastsig;	/* last carrier sample */
249 	double	fdelay;		/* filter delay */
250 	int	decim;		/* sample decimation factor */
251 	int	envphase;	/* envelope phase */
252 	int	envptr;		/* envelope phase pointer */
253 	int	envsw;		/* envelope state */
254 	int	envxing;	/* envelope slice crossing */
255 	int	tc;		/* time constant */
256 	int	tcount;		/* time constant counter */
257 	int	badcnt;		/* decimation interval counter */
258 
259 	/*
260 	 * Decoder variables
261 	 */
262 	int	pulse;		/* cycle counter */
263 	int	cycles;		/* carrier cycles */
264 	int	dcycles;	/* data cycles */
265 	int	lastbit;	/* last code element */
266 	int	second;		/* previous second */
267 	int	bitcnt;		/* bit count in frame */
268 	int	frmcnt;		/* bit count in second */
269 	int	xptr;		/* timecode pointer */
270 	int	bits;		/* demodulated bits */
271 };
272 
273 /*
274  * Function prototypes
275  */
276 static	int	irig_start	(int, struct peer *);
277 static	void	irig_shutdown	(int, struct peer *);
278 static	void	irig_receive	(struct recvbuf *);
279 static	void	irig_poll	(int, struct peer *);
280 
281 /*
282  * More function prototypes
283  */
284 static	void	irig_base	(struct peer *, double);
285 static	void	irig_rf		(struct peer *, double);
286 static	void	irig_baud	(struct peer *, int);
287 static	void	irig_decode	(struct peer *, int);
288 static	void	irig_gain	(struct peer *);
289 
290 /*
291  * Transfer vector
292  */
293 struct	refclock refclock_irig = {
294 	irig_start,		/* start up driver */
295 	irig_shutdown,		/* shut down driver */
296 	irig_poll,		/* transmit poll message */
297 	noentry,		/* not used (old irig_control) */
298 	noentry,		/* initialize driver (not used) */
299 	noentry,		/* not used (old irig_buginfo) */
300 	NOFLAGS			/* not used */
301 };
302 
303 
304 /*
305  * irig_start - open the devices and initialize data for processing
306  */
307 static int
308 irig_start(
309 	int	unit,		/* instance number (used for PCM) */
310 	struct peer *peer	/* peer structure pointer */
311 	)
312 {
313 	struct refclockproc *pp;
314 	struct irigunit *up;
315 
316 	/*
317 	 * Local variables
318 	 */
319 	int	fd;		/* file descriptor */
320 	int	i;		/* index */
321 	double	step;		/* codec adjustment */
322 
323 	/*
324 	 * Open audio device
325 	 */
326 	fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit);
327 	if (fd < 0)
328 		return (0);
329 #ifdef DEBUG
330 	if (debug)
331 		audio_show();
332 #endif
333 
334 	/*
335 	 * Allocate and initialize unit structure
336 	 */
337 	if (!(up = (struct irigunit *)
338 	      emalloc(sizeof(struct irigunit)))) {
339 		(void) close(fd);
340 		return (0);
341 	}
342 	memset((char *)up, 0, sizeof(struct irigunit));
343 	pp = peer->procptr;
344 	pp->unitptr = (caddr_t)up;
345 	pp->io.clock_recv = irig_receive;
346 	pp->io.srcclock = (caddr_t)peer;
347 	pp->io.datalen = 0;
348 	pp->io.fd = fd;
349 	if (!io_addclock(&pp->io)) {
350 		(void)close(fd);
351 		free(up);
352 		return (0);
353 	}
354 
355 	/*
356 	 * Initialize miscellaneous variables
357 	 */
358 	peer->precision = PRECISION;
359 	pp->clockdesc = DESCRIPTION;
360 	memcpy((char *)&pp->refid, REFID, 4);
361 	up->tc = MINTC;
362 	up->decim = 1;
363 	up->gain = 127;
364 
365 	/*
366 	 * The companded samples are encoded sign-magnitude. The table
367 	 * contains all the 256 values in the interest of speed.
368 	 */
369 	up->comp[0] = up->comp[OFFSET] = 0.;
370 	up->comp[1] = 1; up->comp[OFFSET + 1] = -1.;
371 	up->comp[2] = 3; up->comp[OFFSET + 2] = -3.;
372 	step = 2.;
373 	for (i = 3; i < OFFSET; i++) {
374 		up->comp[i] = up->comp[i - 1] + step;
375 		up->comp[OFFSET + i] = -up->comp[i];
376                 if (i % 16 == 0)
377 			step *= 2.;
378 	}
379 	DTOLFP(1. / SECOND, &up->tick);
380 	return (1);
381 }
382 
383 
384 /*
385  * irig_shutdown - shut down the clock
386  */
387 static void
388 irig_shutdown(
389 	int	unit,		/* instance number (not used) */
390 	struct peer *peer	/* peer structure pointer */
391 	)
392 {
393 	struct refclockproc *pp;
394 	struct irigunit *up;
395 
396 	pp = peer->procptr;
397 	up = (struct irigunit *)pp->unitptr;
398 	io_closeclock(&pp->io);
399 	free(up);
400 }
401 
402 
403 /*
404  * irig_receive - receive data from the audio device
405  *
406  * This routine reads input samples and adjusts the logical clock to
407  * track the irig clock by dropping or duplicating codec samples.
408  */
409 static void
410 irig_receive(
411 	struct recvbuf *rbufp	/* receive buffer structure pointer */
412 	)
413 {
414 	struct peer *peer;
415 	struct refclockproc *pp;
416 	struct irigunit *up;
417 
418 	/*
419 	 * Local variables
420 	 */
421 	double	sample;		/* codec sample */
422 	u_char	*dpt;		/* buffer pointer */
423 	int	bufcnt;		/* buffer counter */
424 	l_fp	ltemp;		/* l_fp temp */
425 
426 	peer = (struct peer *)rbufp->recv_srcclock;
427 	pp = peer->procptr;
428 	up = (struct irigunit *)pp->unitptr;
429 
430 	/*
431 	 * Main loop - read until there ain't no more. Note codec
432 	 * samples are bit-inverted.
433 	 */
434 	DTOLFP((double)rbufp->recv_length / SECOND, &ltemp);
435 	L_SUB(&rbufp->recv_time, &ltemp);
436 	up->timestamp = rbufp->recv_time;
437 	dpt = rbufp->recv_buffer;
438 	for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
439 		sample = up->comp[~*dpt++ & 0xff];
440 
441 		/*
442 		 * Variable frequency oscillator. The codec oscillator
443 		 * runs at the nominal rate of 8000 samples per second,
444 		 * or 125 us per sample. A frequency change of one unit
445 		 * results in either duplicating or deleting one sample
446 		 * per second, which results in a frequency change of
447 		 * 125 PPM.
448 		 */
449 		up->phase += (up->freq + clock_codec) / SECOND;
450 		up->phase += pp->fudgetime2 / 1e6;
451 		if (up->phase >= .5) {
452 			up->phase -= 1.;
453 		} else if (up->phase < -.5) {
454 			up->phase += 1.;
455 			irig_rf(peer, sample);
456 			irig_rf(peer, sample);
457 		} else {
458 			irig_rf(peer, sample);
459 		}
460 		L_ADD(&up->timestamp, &up->tick);
461 		sample = fabs(sample);
462 		if (sample > up->signal)
463 			up->signal = sample;
464 			up->signal += (sample - up->signal) /
465 			    1000;
466 
467 		/*
468 		 * Once each second, determine the IRIG format and gain.
469 		 */
470 		up->seccnt = (up->seccnt + 1) % SECOND;
471 		if (up->seccnt == 0) {
472 			if (up->irig_b > up->irig_e) {
473 				up->decim = 1;
474 				up->fdelay = IRIG_B;
475 			} else {
476 				up->decim = 10;
477 				up->fdelay = IRIG_E;
478 			}
479 			up->irig_b = up->irig_e = 0;
480 			irig_gain(peer);
481 
482 		}
483 	}
484 
485 	/*
486 	 * Set the input port and monitor gain for the next buffer.
487 	 */
488 	if (pp->sloppyclockflag & CLK_FLAG2)
489 		up->port = 2;
490 	else
491 		up->port = 1;
492 	if (pp->sloppyclockflag & CLK_FLAG3)
493 		up->mongain = MONGAIN;
494 	else
495 		up->mongain = 0;
496 }
497 
498 
499 /*
500  * irig_rf - RF processing
501  *
502  * This routine filters the RF signal using a bandass filter for IRIG-B
503  * and a lowpass filter for IRIG-E. In case of IRIG-E, the samples are
504  * decimated by a factor of ten. Note that the codec filters function as
505  * roofing filters to attenuate both the high and low ends of the
506  * passband. IIR filter coefficients were determined using Matlab Signal
507  * Processing Toolkit.
508  */
509 static void
510 irig_rf(
511 	struct peer *peer,	/* peer structure pointer */
512 	double	sample		/* current signal sample */
513 	)
514 {
515 	struct refclockproc *pp;
516 	struct irigunit *up;
517 
518 	/*
519 	 * Local variables
520 	 */
521 	double	irig_b, irig_e;	/* irig filter outputs */
522 
523 	pp = peer->procptr;
524 	up = (struct irigunit *)pp->unitptr;
525 
526 	/*
527 	 * IRIG-B filter. Matlab 4th-order IIR elliptic, 800-1200 Hz
528 	 * bandpass, 0.3 dB passband ripple, -50 dB stopband ripple,
529 	 * phase delay 1.03 ms.
530 	 */
531 	irig_b = (up->bpf[8] = up->bpf[7]) * 6.505491e-001;
532 	irig_b += (up->bpf[7] = up->bpf[6]) * -3.875180e+000;
533 	irig_b += (up->bpf[6] = up->bpf[5]) * 1.151180e+001;
534 	irig_b += (up->bpf[5] = up->bpf[4]) * -2.141264e+001;
535 	irig_b += (up->bpf[4] = up->bpf[3]) * 2.712837e+001;
536 	irig_b += (up->bpf[3] = up->bpf[2]) * -2.384486e+001;
537 	irig_b += (up->bpf[2] = up->bpf[1]) * 1.427663e+001;
538 	irig_b += (up->bpf[1] = up->bpf[0]) * -5.352734e+000;
539 	up->bpf[0] = sample - irig_b;
540 	irig_b = up->bpf[0] * 4.952157e-003
541 	    + up->bpf[1] * -2.055878e-002
542 	    + up->bpf[2] * 4.401413e-002
543 	    + up->bpf[3] * -6.558851e-002
544 	    + up->bpf[4] * 7.462108e-002
545 	    + up->bpf[5] * -6.558851e-002
546 	    + up->bpf[6] * 4.401413e-002
547 	    + up->bpf[7] * -2.055878e-002
548 	    + up->bpf[8] * 4.952157e-003;
549 	up->irig_b += irig_b * irig_b;
550 
551 	/*
552 	 * IRIG-E filter. Matlab 4th-order IIR elliptic, 130-Hz lowpass,
553 	 * 0.3 dB passband ripple, -50 dB stopband ripple, phase delay
554 	 * 3.47 ms.
555 	 */
556 	irig_e = (up->lpf[4] = up->lpf[3]) * 8.694604e-001;
557 	irig_e += (up->lpf[3] = up->lpf[2]) * -3.589893e+000;
558 	irig_e += (up->lpf[2] = up->lpf[1]) * 5.570154e+000;
559 	irig_e += (up->lpf[1] = up->lpf[0]) * -3.849667e+000;
560 	up->lpf[0] = sample - irig_e;
561 	irig_e = up->lpf[0] * 3.215696e-003
562 	    + up->lpf[1] * -1.174951e-002
563 	    + up->lpf[2] * 1.712074e-002
564 	    + up->lpf[3] * -1.174951e-002
565 	    + up->lpf[4] * 3.215696e-003;
566 	up->irig_e += irig_e * irig_e;
567 
568 	/*
569 	 * Decimate by a factor of either 1 (IRIG-B) or 10 (IRIG-E).
570 	 */
571 	up->badcnt = (up->badcnt + 1) % up->decim;
572 	if (up->badcnt == 0) {
573 		if (up->decim == 1)
574 			irig_base(peer, irig_b);
575 		else
576 			irig_base(peer, irig_e);
577 	}
578 }
579 
580 /*
581  * irig_base - baseband processing
582  *
583  * This routine processes the baseband signal and demodulates the AM
584  * carrier using a synchronous detector. It then synchronizes to the
585  * data frame at the baud rate and decodes the width-modulated data
586  * pulses.
587  */
588 static void
589 irig_base(
590 	struct peer *peer,	/* peer structure pointer */
591 	double	sample		/* current signal sample */
592 	)
593 {
594 	struct refclockproc *pp;
595 	struct irigunit *up;
596 
597 	/*
598 	 * Local variables
599 	 */
600 	double	lope;		/* integrator output */
601 	double	env;		/* envelope detector output */
602 	double	dtemp;
603 	int	carphase;	/* carrier phase */
604 
605 	pp = peer->procptr;
606 	up = (struct irigunit *)pp->unitptr;
607 
608 	/*
609 	 * Synchronous baud integrator. Corresponding samples of current
610 	 * and past baud intervals are integrated to refine the envelope
611 	 * amplitude and phase estimate. We keep one cycle (1 ms) of the
612 	 * raw data and one baud (10 ms) of the integrated data.
613 	 */
614 	up->envphase = (up->envphase + 1) % BAUD;
615 	up->integ[up->envphase] += (sample - up->integ[up->envphase]) /
616 	    (5 * up->tc);
617 	lope = up->integ[up->envphase];
618 	carphase = up->envphase % CYCLE;
619 	up->lastenv[carphase] = sample;
620 	up->lastint[carphase] = lope;
621 
622 	/*
623 	 * Phase detector. Find the negative-going zero crossing
624 	 * relative to sample 4 in the 8-sample sycle. A phase change of
625 	 * 360 degrees produces an output change of one unit.
626 	 */
627 	if (up->lastsig > 0 && lope <= 0)
628 		up->zxing += (double)(carphase - 4) / CYCLE;
629 	up->lastsig = lope;
630 
631 	/*
632 	 * End of the baud. Update signal/noise estimates and PLL
633 	 * phase, frequency and time constant.
634 	 */
635 	if (up->envphase == 0) {
636 		up->maxsignal = up->intmax; up->noise = up->intmin;
637 		up->intmin = 1e6; up->intmax = -1e6;
638 		if (up->maxsignal < DRPOUT)
639 			up->errflg |= IRIG_ERR_AMP;
640 		if (up->maxsignal > 0)
641 			up->modndx = (up->maxsignal - up->noise) /
642 			    up->maxsignal;
643  		else
644 			up->modndx = 0;
645 		if (up->modndx < MODMIN)
646 			up->errflg |= IRIG_ERR_MOD;
647 		if (up->errflg & (IRIG_ERR_AMP | IRIG_ERR_FREQ |
648 		   IRIG_ERR_MOD | IRIG_ERR_SYNCH)) {
649 			up->tc = MINTC;
650 			up->tcount = 0;
651 		}
652 
653 		/*
654 		 * Update PLL phase and frequency. The PLL time constant
655 		 * is set initially to stabilize the frequency within a
656 		 * minute or two, then increases to the maximum. The
657 		 * frequency is clamped so that the PLL capture range
658 		 * cannot be exceeded.
659 		 */
660 		dtemp = up->zxing * up->decim / BAUD;
661 		up->yxing = dtemp;
662 		up->zxing = 0.;
663 		up->phase += dtemp / up->tc;
664 		up->freq += dtemp / (4. * up->tc * up->tc);
665 		if (up->freq > MAXFREQ) {
666 			up->freq = MAXFREQ;
667 			up->errflg |= IRIG_ERR_FREQ;
668 		} else if (up->freq < -MAXFREQ) {
669 			up->freq = -MAXFREQ;
670 			up->errflg |= IRIG_ERR_FREQ;
671 		}
672 	}
673 
674 	/*
675 	 * Synchronous demodulator. There are eight samples in the cycle
676 	 * and ten cycles in the baud. Since the PLL has aligned the
677 	 * negative-going zero crossing at sample 4, the maximum
678 	 * amplitude is at sample 2 and minimum at sample 6. The
679 	 * beginning of the data pulse is determined from the integrated
680 	 * samples, while the end of the pulse is determined from the
681 	 * raw samples. The raw data bits are demodulated relative to
682 	 * the slice level and left-shifted in the decoding register.
683 	 */
684 	if (carphase != 7)
685 		return;
686 
687 	lope = (up->lastint[2] - up->lastint[6]) / 2.;
688 	if (lope > up->intmax)
689 		up->intmax = lope;
690 	if (lope < up->intmin)
691 		up->intmin = lope;
692 
693 	/*
694 	 * Pulse code demodulator and reference timestamp. The decoder
695 	 * looks for a sequence of ten bits; the first two bits must be
696 	 * one, the last two bits must be zero. Frame synch is asserted
697 	 * when three correct frames have been found.
698 	 */
699 	up->pulse = (up->pulse + 1) % 10;
700 	up->cycles <<= 1;
701 	if (lope >= (up->maxsignal + up->noise) / 2.)
702 		up->cycles |= 1;
703 	if ((up->cycles & 0x303c0f03) == 0x300c0300) {
704 		if (up->pulse != 0)
705 			up->errflg |= IRIG_ERR_SYNCH;
706 		up->pulse = 0;
707 	}
708 
709 	/*
710 	 * Assemble the baud and max/min to get the slice level for the
711 	 * next baud. The slice level is based on the maximum over the
712 	 * first two bits and the minimum over the last two bits, with
713 	 * the slice level halfway between the maximum and minimum.
714 	 */
715 	env = (up->lastenv[2] - up->lastenv[6]) / 2.;
716 	up->dcycles <<= 1;
717 	if (env >= up->slice)
718 		up->dcycles |= 1;
719 	switch(up->pulse) {
720 
721 	case 0:
722 		irig_baud(peer, up->dcycles);
723 		if (env < up->envmin)
724 			up->envmin = env;
725 		up->slice = (up->envmax + up->envmin) / 2;
726 		up->envmin = 1e6; up->envmax = -1e6;
727 		break;
728 
729 	case 1:
730 		up->envmax = env;
731 		break;
732 
733 	case 2:
734 		if (env > up->envmax)
735 			up->envmax = env;
736 		break;
737 
738 	case 9:
739 		up->envmin = env;
740 		break;
741 	}
742 }
743 
744 /*
745  * irig_baud - update the PLL and decode the pulse-width signal
746  */
747 static void
748 irig_baud(
749 	struct peer *peer,	/* peer structure pointer */
750 	int	bits		/* decoded bits */
751 	)
752 {
753 	struct refclockproc *pp;
754 	struct irigunit *up;
755 	double	dtemp;
756 	l_fp	ltemp;
757 
758         pp = peer->procptr;
759 	up = (struct irigunit *)pp->unitptr;
760 
761 	/*
762 	 * The PLL time constant starts out small, in order to
763 	 * sustain a frequency tolerance of 250 PPM. It
764 	 * gradually increases as the loop settles down. Note
765 	 * that small wiggles are not believed, unless they
766 	 * persist for lots of samples.
767 	 */
768 	up->exing = -up->yxing;
769 	if (fabs(up->envxing - up->envphase) <= 1) {
770 		up->tcount++;
771 		if (up->tcount > 20 * up->tc) {
772 			up->tc++;
773 			if (up->tc > MAXTC)
774 				up->tc = MAXTC;
775 			up->tcount = 0;
776 			up->envxing = up->envphase;
777 		} else {
778 			up->exing -= up->envxing - up->envphase;
779 		}
780 	} else {
781 		up->tcount = 0;
782 		up->envxing = up->envphase;
783 	}
784 
785 	/*
786 	 * Strike the baud timestamp as the positive zero crossing of
787 	 * the first bit, accounting for the codec delay and filter
788 	 * delay.
789 	 */
790 	up->prvstamp = up->chrstamp;
791 	dtemp = up->decim * (up->exing / SECOND) + up->fdelay;
792 	DTOLFP(dtemp, &ltemp);
793 	up->chrstamp = up->timestamp;
794 	L_SUB(&up->chrstamp, &ltemp);
795 
796 	/*
797 	 * The data bits are collected in ten-bit bauds. The first two
798 	 * bits are not used. The resulting patterns represent runs of
799 	 * 0-1 bits (0), 2-4 bits (1) and 5-7 bits (PI). The remaining
800 	 * 8-bit run represents a soft error and is treated as 0.
801 	 */
802 	switch (up->dcycles & 0xff) {
803 
804 	case 0x00:		/* 0-1 bits (0) */
805 	case 0x80:
806 		irig_decode(peer, BIT0);
807 		break;
808 
809 	case 0xc0:		/* 2-4 bits (1) */
810 	case 0xe0:
811 	case 0xf0:
812 		irig_decode(peer, BIT1);
813 		break;
814 
815 	case 0xf8:		/* (5-7 bits (PI) */
816 	case 0xfc:
817 	case 0xfe:
818 		irig_decode(peer, BITP);
819 		break;
820 
821 	default:		/* 8 bits (error) */
822 		irig_decode(peer, BIT0);
823 		up->errflg |= IRIG_ERR_DECODE;
824 	}
825 }
826 
827 
828 /*
829  * irig_decode - decode the data
830  *
831  * This routine assembles bauds into digits, digits into frames and
832  * frames into the timecode fields. Bits can have values of zero, one
833  * or position identifier. There are four bits per digit, ten digits per
834  * frame and ten frames per second.
835  */
836 static void
837 irig_decode(
838 	struct	peer *peer,	/* peer structure pointer */
839 	int	bit		/* data bit (0, 1 or 2) */
840 	)
841 {
842 	struct refclockproc *pp;
843 	struct irigunit *up;
844 
845 	/*
846 	 * Local variables
847 	 */
848 	int	syncdig;	/* sync digit (Spectracom) */
849 	char	sbs[6];		/* binary seconds since 0h */
850 	char	spare[2];	/* mulligan digits */
851 	int	temp;
852 
853         pp = peer->procptr;
854 	up = (struct irigunit *)pp->unitptr;
855 
856 	/*
857 	 * Assemble frame bits.
858 	 */
859 	up->bits >>= 1;
860 	if (bit == BIT1) {
861 		up->bits |= 0x200;
862 	} else if (bit == BITP && up->lastbit == BITP) {
863 
864 		/*
865 		 * Frame sync - two adjacent position identifiers, which
866 		 * mark the beginning of the second. The reference time
867 		 * is the beginning of the second position identifier,
868 		 * so copy the character timestamp to the reference
869 		 * timestamp.
870 		 */
871 		if (up->frmcnt != 1)
872 			up->errflg |= IRIG_ERR_SYNCH;
873 		up->frmcnt = 1;
874 		up->refstamp = up->prvstamp;
875 	}
876 	up->lastbit = bit;
877 	if (up->frmcnt % SUBFLD == 0) {
878 
879 		/*
880 		 * End of frame. Encode two hexadecimal digits in
881 		 * little-endian timecode field. Note frame 1 is shifted
882 		 * right one bit to account for the marker PI.
883 		 */
884 		temp = up->bits;
885 		if (up->frmcnt == 10)
886 			temp >>= 1;
887 		if (up->xptr >= 2) {
888 			up->timecode[--up->xptr] = hexchar[temp & 0xf];
889 			up->timecode[--up->xptr] = hexchar[(temp >> 5) &
890 			    0xf];
891 		}
892 		if (up->frmcnt == 0) {
893 
894 			/*
895 			 * End of second. Decode the timecode and wind
896 			 * the clock. Not all IRIG generators have the
897 			 * year; if so, it is nonzero after year 2000.
898 			 * Not all have the hardware status bit; if so,
899 			 * it is lit when the source is okay and dim
900 			 * when bad. We watch this only if the year is
901 			 * nonzero. Not all are configured for signature
902 			 * control. If so, all BCD digits are set to
903 			 * zero if the source is bad. In this case the
904 			 * refclock_process() will reject the timecode
905 			 * as invalid.
906 			 */
907 			up->xptr = 2 * SUBFLD;
908 			if (sscanf((char *)up->timecode,
909 			   "%6s%2d%1d%2s%3d%2d%2d%2d", sbs, &pp->year,
910 			    &syncdig, spare, &pp->day, &pp->hour,
911 			    &pp->minute, &pp->second) != 8)
912 				pp->leap = LEAP_NOTINSYNC;
913 			else
914 				pp->leap = LEAP_NOWARNING;
915 			up->second = (up->second + up->decim) % 60;
916 
917 			/*
918 			 * Raise an alarm if the day field is zero,
919 			 * which happens when signature control is
920 			 * enabled and the device has lost
921 			 * synchronization. Raise an alarm if the year
922 			 * field is nonzero and the sync indicator is
923 			 * zero, which happens when a Spectracom radio
924 			 * has lost synchronization. Raise an alarm if
925 			 * the expected second does not agree with the
926 			 * decoded second, which happens with a garbled
927 			 * IRIG signal. We are very particular.
928 			 */
929 			if (pp->day == 0 || (pp->year != 0 && syncdig ==
930 			    0))
931 				up->errflg |= IRIG_ERR_SIGERR;
932 			if (pp->second != up->second)
933 				up->errflg |= IRIG_ERR_CHECK;
934 			up->second = pp->second;
935 
936 			/*
937 			 * Wind the clock only if there are no errors
938 			 * and the time constant has reached the
939 			 * maximum.
940 			 */
941 			if (up->errflg == 0 && up->tc == MAXTC) {
942 				pp->lastref = pp->lastrec;
943 				pp->lastrec = up->refstamp;
944 				if (!refclock_process(pp))
945 					refclock_report(peer,
946 					    CEVNT_BADTIME);
947 			}
948 			sprintf(pp->a_lastcode,
949 			    "%02x %02d %03d %02d:%02d:%02d %4.0f %3d %6.3f %2d %6.2f %6.1f %s",
950 			    up->errflg, pp->year, pp->day,
951 			    pp->hour, pp->minute, pp->second,
952 			    up->maxsignal, up->gain, up->modndx,
953 			    up->tc, up->exing * 1e6 / SECOND, up->freq *
954 			    1e6 / SECOND, ulfptoa(&pp->lastrec, 6));
955 			pp->lencode = strlen(pp->a_lastcode);
956 			up->errflg = 0;
957 			if (pp->sloppyclockflag & CLK_FLAG4) {
958 				record_clock_stats(&peer->srcadr,
959 				    pp->a_lastcode);
960 #ifdef DEBUG
961 				if (debug)
962 					printf("irig %s\n",
963 					    pp->a_lastcode);
964 #endif /* DEBUG */
965 			}
966 		}
967 	}
968 	up->frmcnt = (up->frmcnt + 1) % FIELD;
969 }
970 
971 
972 /*
973  * irig_poll - called by the transmit procedure
974  *
975  * This routine sweeps up the timecode updates since the last poll. For
976  * IRIG-B there should be at least 60 updates; for IRIG-E there should
977  * be at least 6. If nothing is heard, a timeout event is declared.
978  */
979 static void
980 irig_poll(
981 	int	unit,		/* instance number (not used) */
982 	struct peer *peer	/* peer structure pointer */
983 	)
984 {
985 	struct refclockproc *pp;
986 	struct irigunit *up;
987 
988 	pp = peer->procptr;
989 	up = (struct irigunit *)pp->unitptr;
990 
991 	if (pp->coderecv == pp->codeproc) {
992 		refclock_report(peer, CEVNT_TIMEOUT);
993 		return;
994 
995 	}
996 	refclock_receive(peer);
997 	if (!(pp->sloppyclockflag & CLK_FLAG4)) {
998 		record_clock_stats(&peer->srcadr, pp->a_lastcode);
999 #ifdef DEBUG
1000 		if (debug)
1001 			printf("irig %s\n", pp->a_lastcode);
1002 #endif /* DEBUG */
1003 	}
1004 	pp->polls++;
1005 
1006 }
1007 
1008 
1009 /*
1010  * irig_gain - adjust codec gain
1011  *
1012  * This routine is called at the end of each second. It uses the AGC to
1013  * bradket the maximum signal level between MINAMP and MAXAMP to avoid
1014  * hunting. The routine also jiggles the input port and selectively
1015  * mutes the monitor.
1016  */
1017 static void
1018 irig_gain(
1019 	struct peer *peer	/* peer structure pointer */
1020 	)
1021 {
1022 	struct refclockproc *pp;
1023 	struct irigunit *up;
1024 
1025 	pp = peer->procptr;
1026 	up = (struct irigunit *)pp->unitptr;
1027 
1028 	/*
1029 	 * Apparently, the codec uses only the high order bits of the
1030 	 * gain control field. Thus, it may take awhile for changes to
1031 	 * wiggle the hardware bits.
1032 	 */
1033 	if (up->maxsignal < MINAMP) {
1034 		up->gain += 4;
1035 		if (up->gain > MAXGAIN)
1036 			up->gain = MAXGAIN;
1037 	} else if (up->maxsignal > MAXAMP) {
1038 		up->gain -= 4;
1039 		if (up->gain < 0)
1040 			up->gain = 0;
1041 	}
1042 	audio_gain(up->gain, up->mongain, up->port);
1043 }
1044 
1045 
1046 #else
1047 int refclock_irig_bs;
1048 #endif /* REFCLOCK */
1049