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