xref: /freebsd/contrib/ntp/ntpd/refclock_chu.c (revision daf1cffce2e07931f27c6c6998652e90df6ba87e)
1 /*
2  * refclock_chu - clock driver for Canadian CHU time/frequency station
3  */
4 #ifdef HAVE_CONFIG_H
5 #include <config.h>
6 #endif
7 
8 #if defined(REFCLOCK) && defined(CLOCK_CHU)
9 
10 #include <stdio.h>
11 #include <ctype.h>
12 #include <sys/time.h>
13 #include <time.h>
14 #include <math.h>
15 
16 #include "ntpd.h"
17 #include "ntp_io.h"
18 #include "ntp_refclock.h"
19 #include "ntp_calendar.h"
20 #include "ntp_stdlib.h"
21 #ifdef AUDIO_CHU
22 #include "audio.h"
23 #endif /* AUDIO_CHU */
24 
25 #define ICOM 	1		/* undefine to suppress ICOM code */
26 
27 #ifdef ICOM
28 #include "icom.h"
29 #endif /* ICOM */
30 
31 /*
32  * Audio CHU demodulator/decoder
33  *
34  * This driver synchronizes the computer time using data encoded in
35  * radio transmissions from Canadian time/frequency station CHU in
36  * Ottawa, Ontario. Transmissions are made continuously on 3330 kHz,
37  * 7335 kHz and 14670 kHz in upper sideband, compatible AM mode. An
38  * ordinary shortwave receiver can be tuned manually to one of these
39  * frequencies or, in the case of ICOM receivers, the receiver can be
40  * tuned automatically using this program as propagation conditions
41  * change throughout the day and night.
42  *
43  * The driver receives, demodulates and decodes the radio signals when
44  * connected to the audio codec of a Sun workstation running SunOS or
45  * Solaris, and with a little help, other workstations with similar
46  * codecs or sound cards. In this implementation, only one audio driver
47  * and codec can be supported on a single machine.
48  *
49  * The driver can be compiled to use a Bell 103 compatible modem or
50  * modem chip to receive the radio signal and demodulate the data.
51  * Alternatively, the driver can be compiled to use the audio codec of
52  * the Sun workstation or another with compatible audio drivers. In the
53  * latter case, the driver implements the modem using DSP routines, so
54  * the radio can be connected directly to either the microphone on line
55  * input port. In either case, the driver decodes the data using a
56  * maximum likelihood technique which exploits the considerable degree
57  * of redundancy available to maximize accuracy and minimize errors.
58  *
59  * The CHU time broadcast includes an audio signal compatible with the
60  * Bell 103 modem standard (mark = 2225 Hz, space = 2025 Hz). It consist
61  * of nine, ten-character bursts transmitted at 300 bps and beginning
62  * each second from second 31 to second 39 of the minute. Each character
63  * consists of eight data bits plus one start bit and two stop bits to
64  * encode two hex digits. The burst data consist of five characters (ten
65  * hex digits) followed by a repeat of these characters. In format A,
66  * the characters are repeated in the same polarity; in format B, the
67  * characters are repeated in the opposite polarity.
68  *
69  * Format A bursts are sent at seconds 32 through 39 of the minute in
70  * hex digits
71  *
72  *	6dddhhmmss6dddhhmmss
73  *
74  * The first ten digits encode a frame marker (6) followed by the day
75  * (ddd), hour (hh in UTC), minute (mm) and the second (ss). Since
76  * format A bursts are sent during the third decade of seconds the tens
77  * digit of ss is always 3. The driver uses this to determine correct
78  * burst synchronization. These digits are then repeated with the same
79  * polarity.
80  *
81  * Format B bursts are sent at second 31 of the minute in hex digits
82  *
83  *	xdyyyyttaaxdyyyyttaa
84  *
85  * The first ten digits encode a code (x described below) followed by
86  * the DUT1 (d in deciseconds), Gregorian year (yyyy), difference TAI -
87  * UTC (tt) and daylight time indicator (aa) peculiar to Canada. These
88  * digits are then repeated with inverted polarity.
89  *
90  * The x is coded
91  *
92  * 1 Sign of DUT (0 = +)
93  * 2 Leap second warning. One second will be added.
94  * 4 Leap second warning. One second will be subtracted.
95  * 8 Even parity bit for this nibble.
96  *
97  * By design, the last stop bit of the last character in the burst
98  * coincides with 0.5 second. Since characters have 11 bits and are
99  * transmitted at 300 bps, the last stop bit of the first character
100  * coincides with 0.5 - 10 * 11/300 = 0.133 second. Depending on the
101  * UART, character interrupts can vary somewhere between the beginning
102  * of bit 9 and end of bit 11. These eccentricities can be corrected
103  * along with the radio propagation delay using fudge time 1.
104  *
105  * Debugging aids
106  *
107  * The timecode format used for debugging and data recording includes
108  * data helpful in diagnosing problems with the radio signal and serial
109  * connections. With debugging enabled (-d -d -d on the ntpd command
110  * line), the driver produces one line for each burst in two formats
111  * corresponding to format A and B. Following is format A:
112  *
113  *	n b f s m code
114  *
115  * where n is the number of characters in the burst (0-11), b the burst
116  * distance (0-40), f the field alignment (-1, 0, 1), s the
117  * synchronization distance (0-16), m the burst number (2-9) and code
118  * the burst characters as received. Note that the hex digits in each
119  * character are reversed, so the burst
120  *
121  *	10 38 0 16 9 06851292930685129293
122  *
123  * is interpreted as containing 11 characters with burst distance 38,
124  * field alignment 0, synchronization distance 16 and burst number 9.
125  * The nibble-swapped timecode shows day 58, hour 21, minute 29 and
126  * second 39.
127  *
128  * When the audio driver is compiled, format A is preceded by
129  * the current gain (0-255) and relative signal level (0-9999). The
130  * receiver folume control should be set so that the gain is somewhere
131  * near the middle of the range 0-255, which results in a signal level
132  * near 1000.
133  *
134  * Following is format B:
135  *
136  *	n b s code
137  *
138  * where n is the number of characters in the burst (0-11), b the burst
139  * distance (0-40), s the synchronization distance (0-40) and code the
140  * burst characters as received. Note that the hex digits in each
141  * character are reversed and the last ten digits inverted, so the burst
142  *
143  *	11 40 1091891300ef6e76ecff
144  *
145  * is interpreted as containing 11 characters with burst distance 40.
146  * The nibble-swapped timecode shows DUT1 +0.1 second, year 1998 and TAI
147  * - UTC 31 seconds.
148  *
149  * In addition to the above, the reference timecode is updated and
150  * written to the clockstats file and debug score after the last burst
151  * received in the minute. The format is
152  *
153  *	qq yyyy ddd hh:mm:ss nn dd tt
154  *
155  * where qq are the error flags, as described below, yyyy is the year,
156  * ddd the day, hh:mm:ss the time of day, nn the number of format A
157  * bursts received during the previous minute, dd the decoding distance
158  * and tt the number of timestamps. The error flags are cleared after
159  * every update.
160  *
161  * Fudge factors
162  *
163  * For accuracies better than the low millisceconds, fudge time1 can be
164  * set to the radio propagation delay from CHU to the receiver. This can
165  * be done conviently using the minimuf program. When the modem driver
166  * is compiled, fudge flag3 enables the ppsclock line discipline. Fudge
167  * flag4 causes the dubugging output described above to be recorded in
168  * the clockstats file.
169  *
170  * When the audio driver is compiled, fudge flag2 selects the audio
171  * input port, where 0 is the mike port (default) and 1 is the line-in
172  * port. It does not seem useful to select the compact disc player port.
173  * Fudge flag3 enables audio monitoring of the input signal. For this
174  * purpose, the speaker volume must be set before the driver is started.
175  *
176  * The ICOM code is normally compiled in the driver. It isn't used,
177  * unless the mode keyword on the server configuration command specifies
178  * a nonzero ICOM ID select code. The C-IV trace is turned on if the
179  * debug level is greater than one.
180  */
181 /*
182  * Interface definitions
183  */
184 #define	SPEED232	B300	/* uart speed (300 baud) */
185 #define	PRECISION	(-10)	/* precision assumed (about 1 ms) */
186 #define	REFID		"CHU"	/* reference ID */
187 #ifdef ICOM
188 #define DWELL		5	/* minutes before qsy */
189 #define NCHAN		3	/* number of channels */
190 #endif /* ICOM */
191 #ifdef AUDIO_CHU
192 #define	DESCRIPTION	"CHU Modem Receiver" /* WRU */
193 
194 /*
195  * Audio demodulator definitions
196  */
197 #define SECOND		8000	/* nominal sample rate (Hz) */
198 #define BAUD		300	/* modulation rate (bps) */
199 #define OFFSET		128	/* companded sample offset */
200 #define SIZE		256	/* decompanding table size */
201 #define	MAXSIG		6000.	/* maximum signal level */
202 #define LIMIT		1000.	/* soft limiter threshold */
203 #define AGAIN		6.	/* baseband gain */
204 #define LAG		10	/* discriminator lag */
205 #else
206 #define	DEVICE		"/dev/chu%d" /* device name and unit */
207 #define	SPEED232	B300	/* UART speed (300 baud) */
208 #define	DESCRIPTION	"CHU Audio Receiver" /* WRU */
209 #endif /* AUDIO_CHU */
210 
211 /*
212  * Decoder definitions
213  */
214 #define CHAR		(11. / 300.) /* character time (s) */
215 #define	FUDGE		.185	/* offset to first stop bit (s) */
216 #define BURST		11	/* max characters per burst */
217 #define MINCHAR		9	/* min characters per burst */
218 #define MINDIST		28	/* min burst distance (of 40)  */
219 #define MINSYNC		8	/* min sync distance (of 16) */
220 #define MINSTAMP	20	/* min timestamps (of 60) */
221 #define PANIC		(4 * 1440) /* panic restart */
222 
223 /*
224  * Hex extension codes (>= 16)
225  */
226 #define HEX_MISS	16	/* miss */
227 #define HEX_SOFT	17	/* soft error */
228 #define HEX_HARD	18	/* hard error */
229 
230 /*
231  * Status bits (status)
232  */
233 #define RUNT		0x0001	/* runt burst */
234 #define NOISE		0x0002	/* noise burst */
235 #define BFRAME		0x0004	/* invalid format B frame sync */
236 #define BFORMAT		0x0008	/* invalid format B data */
237 #define AFRAME		0x0010	/* invalid format A frame sync */
238 #define AFORMAT		0x0020	/* invalid format A data */
239 #define DECODE		0x0040	/* invalid data decode */
240 #define STAMP		0x0080	/* too few timestamps */
241 #define INYEAR		0x0100	/* valid B frame */
242 #define INSYNC		0x0200	/* clock synchronized */
243 
244 /*
245  * Alarm status bits (alarm)
246  *
247  * These alarms are set at the end of a minute in which at least one
248  * burst was received. SYNERR is raised if the AFRAME or BFRAME status
249  * bits are set during the minute, FMTERR is raised if the AFORMAT or
250  * BFORMAT status bits are set, DECERR is raised if the DECODE status
251  * bit is set and TSPERR is raised if the STAMP status bit is set.
252  */
253 #define SYNERR		0x01	/* frame sync error */
254 #define FMTERR		0x02	/* data format error */
255 #define DECERR		0x04	/* data decoding error */
256 #define TSPERR		0x08	/* insufficient data */
257 
258 #ifdef AUDIO_CHU
259 struct surv {
260 	double	shift[12];	/* mark register */
261 	double	es_max, es_min;	/* max/min envelope signals */
262 	double	dist;		/* sample distance */
263 	int	uart;		/* decoded character */
264 };
265 #endif /* AUDIO_CHU */
266 
267 /*
268  * CHU unit control structure
269  */
270 struct chuunit {
271 	u_char	decode[20][16];	/* maximum likelihood decoding matrix */
272 	l_fp	cstamp[BURST];	/* character timestamps */
273 	l_fp	tstamp[MAXSTAGE]; /* timestamp samples */
274 	l_fp	timestamp;	/* current buffer timestamp */
275 	l_fp	laststamp;	/* last buffer timestamp */
276 	l_fp	charstamp;	/* character time as a l_fp */
277 	int	errflg;		/* error flags */
278 	int	status;		/* status bits */
279 	int	bufptr;		/* buffer index pointer */
280 	char	ident[10];	/* transmitter frequency */
281 #ifdef ICOM
282 	int	chan;		/* frequency identifier */
283 	int	dwell;		/* dwell minutes at current frequency */
284 	int	fd_icom;	/* ICOM file descriptor */
285 #endif /* ICOM */
286 
287 	/*
288 	 * Character burst variables
289 	 */
290 	int	cbuf[BURST];	/* character buffer */
291 	int	ntstamp;	/* number of timestamp samples */
292 	int	ndx;		/* buffer start index */
293 	int	prevsec;	/* previous burst second */
294 	int	burdist;	/* burst distance */
295 	int	mindist;	/* minimum distance */
296 	int	syndist;	/* sync distance */
297 	int	burstcnt;	/* format A bursts this minute */
298 
299 	/*
300 	 * Format particulars
301 	 */
302 	int	leap;		/* leap/dut code */
303 	int	dut;		/* UTC1 correction */
304 	int	tai;		/* TAI - UTC correction */
305 	int	dst;		/* Canadian DST code */
306 
307 #ifdef AUDIO_CHU
308 	/*
309 	 * Audio codec variables
310 	 */
311 	double	comp[SIZE];	/* decompanding table */
312 	int	port;		/* codec port */
313 	int	gain;		/* codec gain */
314 	int	bufcnt;		/* samples in buffer */
315 	int	clipcnt;	/* sample clip count */
316 	int	seccnt;		/* second interval counter */
317 
318 	/*
319 	 * Modem variables
320 	 */
321 	l_fp	tick;		/* audio sample increment */
322 	double	bpf[9];		/* IIR bandpass filter */
323 	double	disc[LAG];	/* discriminator shift register */
324 	double	lpf[27];	/* FIR lowpass filter */
325 	double	monitor;	/* audio monitor */
326 	double	maxsignal;	/* signal level */
327 	int	discptr;	/* discriminator pointer */
328 
329 	/*
330 	 * Maximum likelihood UART variables
331 	 */
332 	double	baud;		/* baud interval */
333 	struct surv surv[8];	/* UART survivor structures */
334 	int	decptr;		/* decode pointer */
335 	int	dbrk;		/* holdoff counter */
336 #endif /* AUDIO_CHU */
337 };
338 
339 /*
340  * Function prototypes
341  */
342 static	int	chu_start	P((int, struct peer *));
343 static	void	chu_shutdown	P((int, struct peer *));
344 static	void	chu_receive	P((struct recvbuf *));
345 static	void	chu_poll	P((int, struct peer *));
346 
347 /*
348  * More function prototypes
349  */
350 static	void	chu_decode	P((struct peer *, int));
351 static	void	chu_burst	P((struct peer *));
352 static	void	chu_clear	P((struct peer *));
353 static	void	chu_a		P((struct peer *, int));
354 static	void	chu_b		P((struct peer *, int));
355 static	int	chu_dist	P((int, int));
356 static	int	chu_major	P((struct peer *));
357 #ifdef AUDIO_CHU
358 static	void	chu_uart	P((struct surv *, double));
359 static	void	chu_rf		P((struct peer *, double));
360 static	void	chu_gain	P((struct peer *));
361 #endif /* AUDIO_CHU */
362 
363 /*
364  * Global variables
365  */
366 static char hexchar[] = "0123456789abcdef_-=";
367 #ifdef ICOM
368 static double qsy[NCHAN] = {3.33, 7.335, 14.67}; /* frequencies (MHz) */
369 #endif /* ICOM */
370 
371 /*
372  * Transfer vector
373  */
374 struct	refclock refclock_chu = {
375 	chu_start,		/* start up driver */
376 	chu_shutdown,		/* shut down driver */
377 	chu_poll,		/* transmit poll message */
378 	noentry,		/* not used (old chu_control) */
379 	noentry,		/* initialize driver (not used) */
380 	noentry,		/* not used (old chu_buginfo) */
381 	NOFLAGS			/* not used */
382 };
383 
384 
385 /*
386  * chu_start - open the devices and initialize data for processing
387  */
388 static int
389 chu_start(
390 	int	unit,		/* instance number (not used) */
391 	struct peer *peer	/* peer structure pointer */
392 	)
393 {
394 	struct chuunit *up;
395 	struct refclockproc *pp;
396 	int	fd;		/* file descriptor */
397 #ifdef ICOM
398 	char	tbuf[80];	/* trace buffer */
399 	int	temp;
400 #endif /* ICOM */
401 #ifdef AUDIO_CHU
402 	int	i;		/* index */
403 	double	step;		/* codec adjustment */
404 
405 	/*
406 	 * Open audio device
407 	 */
408 	fd = audio_init();
409 	if (fd < 0)
410 		return (0);
411 #ifdef DEBUG
412 	if (debug)
413 		audio_show();
414 #endif
415 #else
416 	char device[20];	/* device name */
417 
418 	/*
419 	 * Open serial port in raw mode.
420 	 */
421 	(void)sprintf(device, DEVICE, unit);
422 	if (!(fd = refclock_open(device, SPEED232, LDISC_RAW))) {
423 		return (0);
424 	}
425 #endif /* AUDIO_CHU */
426 
427 	/*
428 	 * Allocate and initialize unit structure
429 	 */
430 	if (!(up = (struct chuunit *)
431 	      emalloc(sizeof(struct chuunit)))) {
432 		(void) close(fd);
433 		return (0);
434 	}
435 	memset((char *)up, 0, sizeof(struct chuunit));
436 	pp = peer->procptr;
437 	pp->unitptr = (caddr_t)up;
438 	pp->io.clock_recv = chu_receive;
439 	pp->io.srcclock = (caddr_t)peer;
440 	pp->io.datalen = 0;
441 	pp->io.fd = fd;
442 	if (!io_addclock(&pp->io)) {
443 		(void)close(fd);
444 		free(up);
445 		return (0);
446 	}
447 
448 	/*
449 	 * Initialize miscellaneous variables
450 	 */
451 	peer->precision = PRECISION;
452 	pp->clockdesc = DESCRIPTION;
453 	memcpy((char *)&pp->refid, REFID, 4);
454 	DTOLFP(CHAR, &up->charstamp);
455 #ifdef AUDIO_CHU
456 	up->gain = 127;
457 
458 	/*
459 	 * The companded samples are encoded sign-magnitude. The table
460 	 * contains all the 256 values in the interest of speed.
461 	 */
462 	up->comp[0] = up->comp[OFFSET] = 0.;
463 	up->comp[1] = 1; up->comp[OFFSET + 1] = -1.;
464 	up->comp[2] = 3; up->comp[OFFSET + 2] = -3.;
465 	step = 2.;
466 	for (i = 3; i < OFFSET; i++) {
467 		up->comp[i] = up->comp[i - 1] + step;
468 		up->comp[OFFSET + i] = -up->comp[i];
469                 if (i % 16 == 0)
470                 	step *= 2.;
471 	}
472 	DTOLFP(1. / SECOND, &up->tick);
473 #endif /* AUDIO_CHU */
474 	strcpy(up->ident, "X");
475 #ifdef ICOM
476 	temp = 0;
477 #ifdef DEBUG
478 	if (debug > 1)
479 		temp = P_TRACE;
480 #endif
481 	if (peer->ttl > 0) {
482 		if (peer->ttl & 0x80)
483 			up->fd_icom = icom_init("/dev/icom", B1200,
484 			    temp);
485 		else
486 			up->fd_icom = icom_init("/dev/icom", B9600,
487 			    temp);
488 	}
489 	if (up->fd_icom > 0) {
490 		if (icom_freq(up->fd_icom, peer->ttl & 0x7f,
491 		    qsy[up->chan]) < 0) {
492 			NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
493 			    msyslog(LOG_ERR,
494 			    "ICOM bus error; autotune disabled");
495 			up->errflg = CEVNT_FAULT;
496 			close(up->fd_icom);
497 			up->fd_icom = 0;
498 		} else {
499 			sprintf(up->ident, "%.1f", qsy[up->chan]);
500 			sprintf(tbuf, "chu: QSY to %s MHz", up->ident);
501 			record_clock_stats(&peer->srcadr, tbuf);
502 #ifdef DEBUG
503 			if (debug)
504 				printf("%s\n", tbuf);
505 #endif
506 		}
507 	}
508 #endif /* ICOM */
509 	return (1);
510 }
511 
512 
513 /*
514  * chu_shutdown - shut down the clock
515  */
516 static void
517 chu_shutdown(
518 	int	unit,		/* instance number (not used) */
519 	struct peer *peer	/* peer structure pointer */
520 	)
521 {
522 	struct chuunit *up;
523 	struct refclockproc *pp;
524 
525 	pp = peer->procptr;
526 	up = (struct chuunit *)pp->unitptr;
527 	io_closeclock(&pp->io);
528 	if (up->fd_icom > 0)
529 		close(up->fd_icom);
530 	free(up);
531 }
532 
533 #ifdef AUDIO_CHU
534 
535 /*
536  * chu_receive - receive data from the audio device
537  */
538 static void
539 chu_receive(
540 	struct recvbuf *rbufp	/* receive buffer structure pointer */
541 	)
542 {
543 	struct chuunit *up;
544 	struct refclockproc *pp;
545 	struct peer *peer;
546 
547 	double	sample;		/* codec sample */
548 	u_char	*dpt;		/* buffer pointer */
549 	l_fp	ltemp;		/* l_fp temp */
550 	int	isneg;		/* parity flag */
551 	double	dtemp;
552 	int	i, j;
553 
554 	peer = (struct peer *)rbufp->recv_srcclock;
555 	pp = peer->procptr;
556 	up = (struct chuunit *)pp->unitptr;
557 
558 	/*
559 	 * Main loop - read until there ain't no more. Note codec
560 	 * samples are bit-inverted.
561 	 */
562 	up->timestamp = rbufp->recv_time;
563 	up->bufcnt = rbufp->recv_length;
564 	DTOLFP(up->bufcnt * 1. / SECOND, &ltemp);
565 	L_SUB(&up->timestamp, &ltemp);
566 	dpt = (u_char *)&rbufp->recv_space;
567 	for (up->bufptr = 0; up->bufptr < up->bufcnt; up->bufptr++) {
568 		sample = up->comp[~*dpt & 0xff];
569 
570 		/*
571 		 * Clip noise spikes greater than MAXSIG. If no clips,
572 		 * increase the gain a tad; if the clips are too high,
573 		 * decrease a tad.
574 		 */
575 		if (sample > MAXSIG) {
576 			sample = MAXSIG;
577 			up->clipcnt++;
578 		} else if (sample < -MAXSIG) {
579 			sample = -MAXSIG;
580 			up->clipcnt++;
581 		}
582 		up->seccnt = (up->seccnt + 1) % SECOND;
583 		if (up->seccnt == 0) {
584 			if (pp->sloppyclockflag & CLK_FLAG2)
585 				up->port = 2;
586 			else
587 				up->port = 1;
588 			chu_gain(peer);
589 		}
590 		chu_rf(peer, sample);
591 
592 		/*
593 		 * During development, it is handy to have an audio
594 		 * monitor that can be switched to various signals. This
595 		 * code converts the linear signal left in up->monitor
596 		 * to codec format. If we can get the grass out of this
597 		 * thing and improve modem performance, this expensive
598 		 * code will be permanently nixed.
599 		 */
600 		isneg = 0;
601 		dtemp = up->monitor;
602 		if (sample < 0) {
603 			isneg = 1;
604 			dtemp-= dtemp;
605 		}
606 		i = 0;
607 		j = OFFSET >> 1;
608 		while (j != 0) {
609 			if (dtemp > up->comp[i])
610 				i += j;
611 			else if (dtemp < up->comp[i])
612 				i -= j;
613 			else
614 				break;
615 			j >>= 1;
616 		}
617 		if (isneg)
618 			*dpt = ~(i + OFFSET);
619 		else
620 			*dpt = ~i;
621 		dpt++;
622 		L_ADD(&up->timestamp, &up->tick);
623 	}
624 
625 	/*
626 	 * Squawk to the monitor speaker if enabled.
627 	 */
628 	if (pp->sloppyclockflag & CLK_FLAG3)
629 		if (write(pp->io.fd, (u_char *)&rbufp->recv_space,
630 		    (u_int)up->bufcnt) < 0)
631 			perror("chu:");
632 }
633 
634 
635 /*
636  * chu_rf - filter and demodulate the FSK signal
637  *
638  * This routine implements a 300-baud Bell 103 modem with mark 2225 Hz
639  * and space 2025 Hz. It uses a bandpass filter followed by a soft
640  * limiter, FM discriminator and lowpass filter. A maximum likelihood
641  * decoder samples the baseband signal at eight times the baud rate and
642  * detects the start bit of each character.
643  *
644  * The filters are built for speed, which explains the rather clumsy
645  * code. Hopefully, the compiler will efficiently implement the move-
646  * and-muiltiply-and-add operations.
647  */
648 static void
649 chu_rf(
650 	struct peer *peer,	/* peer structure pointer */
651 	double	sample		/* analog sample */
652 	)
653 {
654 	struct refclockproc *pp;
655 	struct chuunit *up;
656 	struct surv *sp;
657 
658 	/*
659 	 * Local variables
660 	 */
661 	double	signal;		/* bandpass signal */
662 	double	limit;		/* limiter signal */
663 	double	disc;		/* discriminator signal */
664 	double	lpf;		/* lowpass signal */
665 	double	span;		/* UART signal span */
666 	double	dist;		/* UART signal distance */
667 	int	i, j;
668 
669 	pp = peer->procptr;
670 	up = (struct chuunit *)pp->unitptr;
671 
672 	/*
673 	 * Bandpass filter. 4th-order elliptic, 500-Hz bandpass centered
674 	 * at 2125 Hz. Passband ripple 0.3 dB, stopband ripple 50 dB.
675 	 */
676 	signal = (up->bpf[8] = up->bpf[7]) * 5.844676e-01;
677 	signal += (up->bpf[7] = up->bpf[6]) * 4.884860e-01;
678 	signal += (up->bpf[6] = up->bpf[5]) * 2.704384e+00;
679 	signal += (up->bpf[5] = up->bpf[4]) * 1.645032e+00;
680 	signal += (up->bpf[4] = up->bpf[3]) * 4.644557e+00;
681 	signal += (up->bpf[3] = up->bpf[2]) * 1.879165e+00;
682 	signal += (up->bpf[2] = up->bpf[1]) * 3.522634e+00;
683 	signal += (up->bpf[1] = up->bpf[0]) * 7.315738e-01;
684 	up->bpf[0] = sample - signal;
685 	signal = up->bpf[0] * 6.176213e-03
686 	    + up->bpf[1] * 3.156599e-03
687 	    + up->bpf[2] * 7.567487e-03
688 	    + up->bpf[3] * 4.344580e-03
689 	    + up->bpf[4] * 1.190128e-02
690 	    + up->bpf[5] * 4.344580e-03
691 	    + up->bpf[6] * 7.567487e-03
692 	    + up->bpf[7] * 3.156599e-03
693 	    + up->bpf[8] * 6.176213e-03;
694 
695 	up->monitor = signal / 4.;	/* note monitor after filter */
696 
697 	/*
698 	 * Soft limiter/discriminator. The 11-sample discriminator lag
699 	 * interval corresponds to three cycles of 2125 Hz, which
700 	 * requires the sample frequency to be 2125 * 11 / 3 = 7791.7
701 	 * Hz. The discriminator output varies +-0.5 interval for input
702 	 * frequency 2025-2225 Hz. However, we don't get to sample at
703 	 * this frequency, so the discriminator output is biased. Life
704 	 * at 8000 Hz sucks.
705 	 */
706 	limit = signal;
707 	if (limit > LIMIT)
708 		limit = LIMIT;
709 	else if (limit < -LIMIT)
710 		limit = -LIMIT;
711 	disc = up->disc[up->discptr] * -limit;
712 	up->disc[up->discptr] = limit;
713 	up->discptr = (up->discptr + 1 ) % LAG;
714 	if (disc >= 0)
715 		disc = sqrt(disc);
716 	else
717 		disc = -sqrt(-disc);
718 
719 	/*
720 	 * Lowpass filter. Raised cosine, Ts = 1 / 300, beta = 0.1.
721 	 */
722 	lpf = (up->lpf[26] = up->lpf[25]) * 2.538771e-02;
723 	lpf += (up->lpf[25] = up->lpf[24]) * 1.084671e-01;
724 	lpf += (up->lpf[24] = up->lpf[23]) * 2.003159e-01;
725 	lpf += (up->lpf[23] = up->lpf[22]) * 2.985303e-01;
726 	lpf += (up->lpf[22] = up->lpf[21]) * 4.003697e-01;
727 	lpf += (up->lpf[21] = up->lpf[20]) * 5.028552e-01;
728 	lpf += (up->lpf[20] = up->lpf[19]) * 6.028795e-01;
729 	lpf += (up->lpf[19] = up->lpf[18]) * 6.973249e-01;
730 	lpf += (up->lpf[18] = up->lpf[17]) * 7.831828e-01;
731 	lpf += (up->lpf[17] = up->lpf[16]) * 8.576717e-01;
732 	lpf += (up->lpf[16] = up->lpf[15]) * 9.183463e-01;
733 	lpf += (up->lpf[15] = up->lpf[14]) * 9.631951e-01;
734 	lpf += (up->lpf[14] = up->lpf[13]) * 9.907208e-01;
735 	lpf += (up->lpf[13] = up->lpf[12]) * 1.000000e+00;
736 	lpf += (up->lpf[12] = up->lpf[11]) * 9.907208e-01;
737 	lpf += (up->lpf[11] = up->lpf[10]) * 9.631951e-01;
738 	lpf += (up->lpf[10] = up->lpf[9]) * 9.183463e-01;
739 	lpf += (up->lpf[9] = up->lpf[8]) * 8.576717e-01;
740 	lpf += (up->lpf[8] = up->lpf[7]) * 7.831828e-01;
741 	lpf += (up->lpf[7] = up->lpf[6]) * 6.973249e-01;
742 	lpf += (up->lpf[6] = up->lpf[5]) * 6.028795e-01;
743 	lpf += (up->lpf[5] = up->lpf[4]) * 5.028552e-01;
744 	lpf += (up->lpf[4] = up->lpf[3]) * 4.003697e-01;
745 	lpf += (up->lpf[3] = up->lpf[2]) * 2.985303e-01;
746 	lpf += (up->lpf[2] = up->lpf[1]) * 2.003159e-01;
747 	lpf += (up->lpf[1] = up->lpf[0]) * 1.084671e-01;
748 	lpf += up->lpf[0] = disc * 2.538771e-02;
749 
750 	/*
751 	 * Maximum likelihood decoder. The UART updates each of the
752 	 * eight survivors and determines the span, slice level and
753 	 * tentative decoded character. Valid 11-bit characters are
754 	 * framed so that bit 1 and bit 11 (stop bits) are mark and bit
755 	 * 2 (start bit) is space. When a valid character is found, the
756 	 * survivor with maximum distance determines the final decoded
757 	 * character.
758 	 */
759 	up->baud += 1. / SECOND;
760 	if (up->baud > 1. / (BAUD * 8.)) {
761 		up->baud -= 1. / (BAUD * 8.);
762 		sp = &up->surv[up->decptr];
763 		span = sp->es_max - sp->es_min;
764 		up->maxsignal += (span - up->maxsignal) / 80.;
765 		if (up->dbrk > 0) {
766 			up->dbrk--;
767 		} else if ((sp->uart & 0x403) == 0x401 && span > 1000.)
768 		    {
769 			dist = 0;
770 			j = 0;
771 			for (i = 0; i < 8; i++) {
772 				if (up->surv[i].dist > dist) {
773 					dist = up->surv[i].dist;
774 					j = i;
775 				}
776 			}
777 			chu_decode(peer, (up->surv[j].uart >> 2) &
778 			    0xff);
779 			up->dbrk = 80;
780 		}
781 		up->decptr = (up->decptr + 1) % 8;
782 		chu_uart(sp, -lpf * AGAIN);
783 	}
784 }
785 
786 
787 /*
788  * chu_uart - maximum likelihood UART
789  *
790  * This routine updates a shift register holding the last 11 envelope
791  * samples. It then computes the slice level and span over these samples
792  * and determines the tentative data bits and distance. The calling
793  * program selects over the last eight survivors the one with maximum
794  * distance to determine the decoded character.
795  */
796 static void
797 chu_uart(
798 	struct surv *sp,	/* survivor structure pointer */
799 	double	sample		/* baseband signal */
800 	)
801 {
802 	double	es_max, es_min;	/* max/min envelope */
803 	double	slice;		/* slice level */
804 	double	dist;		/* distance */
805 	double	dtemp;
806 	int	i;
807 
808 	/*
809 	 * Save the sample and shift right. At the same time, measure
810 	 * the maximum and minimum over all eleven samples.
811 	 */
812 	es_max = -1e6;
813 	es_min = 1e6;
814 	sp->shift[0] = sample;
815 	for (i = 11; i > 0; i--) {
816 		sp->shift[i] = sp->shift[i - 1];
817 		if (sp->shift[i] > es_max)
818 			es_max = sp->shift[i];
819 		if (sp->shift[i] < es_min)
820 			es_min = sp->shift[i];
821 	}
822 
823 	/*
824 	 * Determine the slice level midway beteen the maximum and
825 	 * minimum and the span as the maximum less the minimum. Compute
826 	 * the distance on the assumption the first and last bits must
827 	 * be mark, the second space and the rest either mark or space.
828 	 */
829 	slice = (es_max + es_min) / 2.;
830 	dist = 0;
831 	sp->uart = 0;
832 	for (i = 1; i < 12; i++) {
833 		sp->uart <<= 1;
834 		dtemp = sp->shift[i];
835 		if (dtemp > slice)
836 			sp->uart |= 0x1;
837 		if (i == 1 || i == 11) {
838 			dist += dtemp - es_min;
839 		} else if (i == 10) {
840 			dist += es_max - dtemp;
841 		} else {
842 			if (dtemp > slice)
843 				dist += dtemp - es_min;
844 			else
845 				dist += es_max - dtemp;
846 		}
847 	}
848 	sp->es_max = es_max;
849 	sp->es_min = es_min;
850 	sp->dist = dist / (11 * (es_max - es_min));
851 }
852 
853 
854 #else /* AUDIO_CHU */
855 /*
856  * chu_receive - receive data from the serial interface
857  */
858 static void
859 chu_receive(
860 	struct recvbuf *rbufp	/* receive buffer structure pointer */
861 	)
862 {
863 	struct chuunit *up;
864 	struct refclockproc *pp;
865 	struct peer *peer;
866 
867 	u_char	*dpt;		/* receive buffer pointer */
868 
869 	peer = (struct peer *)rbufp->recv_srcclock;
870 	pp = peer->procptr;
871 	up = (struct chuunit *)pp->unitptr;
872 
873 	/*
874 	 * Initialize pointers and read the timecode and timestamp.
875 	 */
876 	up->timestamp = rbufp->recv_time;
877 	dpt = (u_char *)&rbufp->recv_space;
878 	chu_decode(peer, *dpt);
879 }
880 #endif /* AUDIO_CHU */
881 
882 
883 /*
884  * chu_decode - decode the data
885  */
886 static void
887 chu_decode(
888 	struct peer *peer,	/* peer structure pointer */
889 	int	hexhex		/* data character */
890 	)
891 {
892 	struct refclockproc *pp;
893 	struct chuunit *up;
894 
895 	l_fp	tstmp;		/* timestamp temp */
896 	double	dtemp;
897 
898 	pp = peer->procptr;
899 	up = (struct chuunit *)pp->unitptr;
900 
901 	/*
902 	 * If the interval since the last character is greater than the
903 	 * longest burst, process the last burst and start a new one. If
904 	 * the interval is less than this but greater than two
905 	 * characters, consider this a noise burst and reject it.
906 	 */
907 	tstmp = up->timestamp;
908 	if (L_ISZERO(&up->laststamp))
909 		up->laststamp = up->timestamp;
910 	L_SUB(&tstmp, &up->laststamp);
911 	up->laststamp = up->timestamp;
912 	LFPTOD(&tstmp, dtemp);
913 	if (dtemp > BURST * CHAR) {
914 		chu_burst(peer);
915 		up->ndx = 0;
916 	} else if (dtemp > 2.5 * CHAR) {
917 		up->ndx = 0;
918 	}
919 
920 	/*
921 	 * Append the character to the current burst and append the
922 	 * timestamp to the timestamp list.
923 	 */
924 	if (up->ndx < BURST) {
925 		up->cbuf[up->ndx] = hexhex & 0xff;
926 		up->cstamp[up->ndx] = up->timestamp;
927 		up->ndx++;
928 
929 	}
930 }
931 
932 
933 /*
934  * chu_burst - search for valid burst format
935  */
936 static void
937 chu_burst(
938 	struct peer *peer
939 	)
940 {
941 	struct chuunit *up;
942 	struct refclockproc *pp;
943 
944 	int	i;
945 
946 	pp = peer->procptr;
947 	up = (struct chuunit *)pp->unitptr;
948 
949 	/*
950 	 * Correlate a block of five characters with the next block of
951 	 * five characters. The burst distance is defined as the number
952 	 * of bits that match in the two blocks for format A and that
953 	 * match the inverse for format B.
954 	 */
955 	if (up->ndx < MINCHAR) {
956 		up->status |= RUNT;
957 		return;
958 	}
959 	up->burdist = 0;
960 	for (i = 0; i < 5 && i < up->ndx - 5; i++)
961 		up->burdist += chu_dist(up->cbuf[i], up->cbuf[i + 5]);
962 
963 	/*
964 	 * If the burst distance is at least MINDIST, this must be a
965 	 * format A burst; if the value is not greater than -MINDIST, it
966 	 * must be a format B burst. If the B burst is perfect, we
967 	 * believe it; otherwise, it is a noise burst and of no use to
968 	 * anybody.
969 	 */
970 	if (up->burdist >= MINDIST) {
971 		chu_a(peer, up->ndx);
972 	} else if (up->burdist <= -MINDIST) {
973 		chu_b(peer, up->ndx);
974 	} else {
975 		up->status |= NOISE;
976 		return;
977 	}
978 
979 	/*
980 	 * If this is a valid burst, wait a guard time of ten seconds to
981 	 * allow for more bursts, then arm the poll update routine to
982 	 * process the minute. Don't do this if this is called from the
983 	 * timer interrupt routine.
984 	 */
985 	if (peer->outdate != current_time)
986 		peer->nextdate = current_time + 10;
987 }
988 
989 
990 /*
991  * chu_b - decode format B burst
992  */
993 static void
994 chu_b(
995 	struct peer *peer,
996 	int	nchar
997 	)
998 {
999 	struct	refclockproc *pp;
1000 	struct	chuunit *up;
1001 
1002 	u_char	code[11];	/* decoded timecode */
1003 	char	tbuf[80];	/* trace buffer */
1004 	l_fp	offset;		/* timestamp offset */
1005 	int	i;
1006 
1007 	pp = peer->procptr;
1008 	up = (struct chuunit *)pp->unitptr;
1009 
1010 	/*
1011 	 * In a format B burst, a character is considered valid only if
1012 	 * the first occurrence matches the last occurrence. The burst
1013 	 * is considered valid only if all characters are valid; that
1014 	 * is, only if the distance is 40.
1015 	 */
1016 	sprintf(tbuf, "chuB %04x %2d %2d ", up->status, nchar,
1017 	    -up->burdist);
1018 	for (i = 0; i < nchar; i++)
1019 		sprintf(&tbuf[strlen(tbuf)], "%02x",
1020 		    up->cbuf[i]);
1021 	if (pp->sloppyclockflag & CLK_FLAG4)
1022 		record_clock_stats(&peer->srcadr, tbuf);
1023 #ifdef DEBUG
1024 	if (debug)
1025 		printf("%s\n", tbuf);
1026 #endif
1027 	if (up->burdist > -40) {
1028 		up->status |= BFRAME;
1029 		return;
1030 	}
1031 	up->status |= INYEAR;
1032 
1033 	/*
1034 	 * Convert the burst data to internal format. If this succeeds,
1035 	 * save the timestamps for later.
1036 	 */
1037 	for (i = 0; i < 5; i++) {
1038 		code[2 * i] = hexchar[up->cbuf[i] & 0xf];
1039 		code[2 * i + 1] = hexchar[(up->cbuf[i] >>
1040 		    4) & 0xf];
1041 	}
1042 	if (sscanf((char *)code, "%1x%1d%4d%2d%2x", &up->leap, &up->dut,
1043 	    &pp->year, &up->tai, &up->dst) != 5) {
1044 		up->status |= BFORMAT;
1045 		return;
1046 	}
1047 	if (up->leap & 0x8)
1048 		up->dut = -up->dut;
1049 	offset.l_ui = 31;
1050 	offset.l_f = 0;
1051 	for (i = 0; i < nchar && i < 10; i++) {
1052 		up->tstamp[up->ntstamp] = up->cstamp[i];
1053 		L_SUB(&up->tstamp[up->ntstamp], &offset);
1054 		L_ADD(&offset, &up->charstamp);
1055 		if (up->ntstamp < MAXSTAGE)
1056 			up->ntstamp++;
1057 	}
1058 }
1059 
1060 
1061 /*
1062  * chu_a - decode format A burst
1063  */
1064 static void
1065 chu_a(
1066 	struct peer *peer,
1067 	int nchar
1068 	)
1069 {
1070 	struct refclockproc *pp;
1071 	struct chuunit *up;
1072 
1073 	char	tbuf[80];	/* trace buffer */
1074 	l_fp	offset;		/* timestamp offset */
1075 	int	val;		/* distance */
1076 	int	temp;
1077 	int	i, j, k;
1078 
1079 	pp = peer->procptr;
1080 	up = (struct chuunit *)pp->unitptr;
1081 
1082 	/*
1083 	 * Determine correct burst phase. There are three cases
1084 	 * corresponding to in-phase, one character early or one
1085 	 * character late. These cases are distinguished by the position
1086 	 * of the framing digits x6 at positions 0 and 5 and x3 at
1087 	 * positions 4 and 9. The correct phase is when the distance
1088 	 * relative to the framing digits is maximum. The burst is valid
1089 	 * only if the maximum distance is at least MINSYNC.
1090 	 */
1091 	up->syndist = k = 0;
1092 	val = -16;
1093 	for (i = -1; i < 2; i++) {
1094 		temp = up->cbuf[i + 4] & 0xf;
1095 		if (i >= 0)
1096 			temp |= (up->cbuf[i] & 0xf) << 4;
1097 		val = chu_dist(temp, 0x63);
1098 		temp = (up->cbuf[i + 5] & 0xf) << 4;
1099 		if (i + 9 < nchar)
1100 			temp |= up->cbuf[i + 9] & 0xf;
1101 		val += chu_dist(temp, 0x63);
1102 		if (val > up->syndist) {
1103 			up->syndist = val;
1104 			k = i;
1105 		}
1106 	}
1107 	temp = (up->cbuf[k + 4] >> 4) & 0xf;
1108 	if (temp > 9 || k + 9 >= nchar || temp != ((up->cbuf[k + 9] >>
1109 	    4) & 0xf))
1110 		temp = 0;
1111 #ifdef AUDIO_CHU
1112 	sprintf(tbuf, "chuA %04x %4.0f %2d %2d %2d %2d %1d ",
1113 	    up->status, up->maxsignal, nchar, up->burdist, k,
1114 	    up->syndist, temp);
1115 #else
1116 	sprintf(tbuf, "chuA %04x %2d %2d %2d %2d %1d ", up->status,
1117 	    nchar, up->burdist, k, up->syndist, temp);
1118 #endif /* AUDIO_CHU */
1119 	for (i = 0; i < nchar; i++)
1120 		sprintf(&tbuf[strlen(tbuf)], "%02x",
1121 		    up->cbuf[i]);
1122 	if (pp->sloppyclockflag & CLK_FLAG4)
1123 		record_clock_stats(&peer->srcadr, tbuf);
1124 #ifdef DEBUG
1125 	if (debug)
1126 		printf("%s\n", tbuf);
1127 #endif
1128 	if (up->syndist < MINSYNC) {
1129 		up->status |= AFRAME;
1130 		return;
1131 	}
1132 
1133 	/*
1134 	 * A valid burst requires the first seconds number to match the
1135 	 * last seconds number. If so, the burst timestamps are
1136 	 * corrected to the current minute and saved for later
1137 	 * processing. In addition, the seconds decode is advanced from
1138 	 * the previous burst to the current one.
1139 	 */
1140 	if (temp != 0) {
1141 		offset.l_ui = 30 + temp;
1142 		offset.l_f = 0;
1143 		i = 0;
1144 		if (k < 0)
1145 			offset = up->charstamp;
1146 		else if (k > 0)
1147 			i = 1;
1148 		for (; i < nchar && i < k + 10; i++) {
1149 			up->tstamp[up->ntstamp] = up->cstamp[i];
1150 			L_SUB(&up->tstamp[up->ntstamp], &offset);
1151 			L_ADD(&offset, &up->charstamp);
1152 			if (up->ntstamp < MAXSTAGE)
1153 				up->ntstamp++;
1154 		}
1155 		while (temp > up->prevsec) {
1156 			for (j = 15; j > 0; j--) {
1157 				up->decode[9][j] = up->decode[9][j - 1];
1158 				up->decode[19][j] =
1159 				    up->decode[19][j - 1];
1160 			}
1161 			up->decode[9][j] = up->decode[19][j] = 0;
1162 			up->prevsec++;
1163 		}
1164 	}
1165 	i = -(2 * k);
1166 	for (j = 0; j < nchar; j++) {
1167 		if (i < 0 || i > 19) {
1168 			i += 2;
1169 			continue;
1170 		}
1171 		up->decode[i][up->cbuf[j] & 0xf]++;
1172 		i++;
1173 		up->decode[i][(up->cbuf[j] >> 4) & 0xf]++;
1174 		i++;
1175 	}
1176 	up->burstcnt++;
1177 }
1178 
1179 
1180 /*
1181  * chu_poll - called by the transmit procedure
1182  */
1183 static void
1184 chu_poll(
1185 	int unit,
1186 	struct peer *peer	/* peer structure pointer */
1187 	)
1188 {
1189 	struct refclockproc *pp;
1190 	struct chuunit *up;
1191 	char	synchar, qual, leapchar;
1192 	int	minset;
1193 	int	temp;
1194 #ifdef ICOM
1195 	char	tbuf[80];	/* trace buffer */
1196 #endif /* ICOM */
1197 	pp = peer->procptr;
1198 	up = (struct chuunit *)pp->unitptr;
1199 	if (pp->coderecv == pp->codeproc)
1200 		up->errflg = CEVNT_TIMEOUT;
1201 	else
1202 		pp->polls++;
1203 	minset = ((current_time - peer->update) + 30) / 60;
1204 	if (up->status & INSYNC) {
1205 		if (minset > PANIC)
1206 			up->status = 0;
1207 		else
1208 			peer->reach |= 1;
1209 	}
1210 
1211 	/*
1212 	 * Process the last burst, if still in the burst buffer.
1213 	 * Don't mess with anything if nothing has been heard.
1214 	 */
1215 	chu_burst(peer);
1216 #ifdef ICOM
1217 	if (up->burstcnt > 2) {
1218 		up->dwell = 0;
1219 	} else if (up->dwell < DWELL) {
1220 		up->dwell++;
1221 	} else if (up->fd_icom > 0) {
1222 		up->dwell = 0;
1223 		up->chan = (up->chan + 1) % NCHAN;
1224 		icom_freq(up->fd_icom, peer->ttl & 0x7f, qsy[up->chan]);
1225 		sprintf(up->ident, "%.3f", qsy[up->chan]);
1226 		sprintf(tbuf, "chu: QSY to %s MHz", up->ident);
1227 		record_clock_stats(&peer->srcadr, tbuf);
1228 #ifdef DEBUG
1229 		if (debug)
1230 			printf("%s\n", tbuf);
1231 #endif
1232 	}
1233 #endif /* ICOM */
1234 	if (up->burstcnt == 0)
1235 		return;
1236 	temp = chu_major(peer);
1237 	if (up->status & INYEAR)
1238 		up->status |= INSYNC;
1239 	qual = 0;
1240 	if (up->status & (BFRAME | AFRAME))
1241 		qual |= SYNERR;
1242 	if (up->status & (BFORMAT | AFORMAT))
1243 		qual |= FMTERR;
1244 	if (up->status & DECODE)
1245 		qual |= DECERR;
1246 	if (up->status & STAMP)
1247 		qual |= TSPERR;
1248 	synchar = leapchar = ' ';
1249 	if (!(up->status & INSYNC)) {
1250 		pp->leap = LEAP_NOTINSYNC;
1251 		synchar = '?';
1252 	} else if (up->leap & 0x2) {
1253 		pp->leap = LEAP_ADDSECOND;
1254 		leapchar = 'L';
1255 	} else {
1256 		pp->leap = LEAP_NOWARNING;
1257 	}
1258 #ifdef AUDIO_CHU
1259 	sprintf(pp->a_lastcode,
1260 	    "%c%1X %4d %3d %02d:%02d:%02d.000 %c%x %+d %d %d %s %d %d %d %d",
1261 	    synchar, qual, pp->year, pp->day, pp->hour, pp->minute,
1262 	    pp->second, leapchar, up->dst, up->dut, minset, up->gain,
1263 	    up->ident, up->tai, up->burstcnt, up->mindist, up->ntstamp);
1264 #else
1265 	sprintf(pp->a_lastcode,
1266 	    "%c%1X %4d %3d %02d:%02d:%02d.000 %c%x %+d %d %s %d %d %d %d",
1267 	    synchar, qual, pp->year, pp->day, pp->hour, pp->minute,
1268 	    pp->second, leapchar, up->dst, up->dut, minset,
1269 	    up->ident, up->tai, up->burstcnt, up->mindist, up->ntstamp);
1270 #endif /* AUDIO_CHU */
1271 	pp->lencode = strlen(pp->a_lastcode);
1272 
1273 	/*
1274 	 * If timestamps have been stuffed, the timecode is ipso fatso
1275 	 * correct and can be selected to discipline the clock.
1276 	 */
1277 	if (temp > 0) {
1278 		record_clock_stats(&peer->srcadr, pp->a_lastcode);
1279 		refclock_receive(peer);
1280 	} else if (pp->sloppyclockflag & CLK_FLAG4) {
1281 		record_clock_stats(&peer->srcadr, pp->a_lastcode);
1282 	}
1283 #ifdef DEBUG
1284 	if (debug)
1285 		printf("chu: timecode %d %s\n", pp->lencode,
1286 		    pp->a_lastcode);
1287 #endif
1288 	chu_clear(peer);
1289 	if (up->errflg)
1290 		refclock_report(peer, up->errflg);
1291 	up->errflg = 0;
1292 }
1293 
1294 
1295 /*
1296  * chu_major - majority decoder
1297  */
1298 static int
1299 chu_major(
1300 	struct peer *peer	/* peer structure pointer */
1301 	)
1302 {
1303 	struct refclockproc *pp;
1304 	struct chuunit *up;
1305 
1306 	u_char	code[11];	/* decoded timecode */
1307 	l_fp	toffset, offset; /* l_fp temps */
1308 	int	val1, val2;	/* maximum distance */
1309 	int	synchar;	/* stray cat */
1310 	double	dtemp;
1311 	int	temp;
1312 	int	i, j, k;
1313 
1314 	pp = peer->procptr;
1315 	up = (struct chuunit *)pp->unitptr;
1316 
1317 	/*
1318 	 * Majority decoder. Each burst encodes two replications at each
1319 	 * digit position in the timecode. Each row of the decoding
1320 	 * matrix encodes the number of occurences of each digit found
1321 	 * at the corresponding position. The maximum over all
1322 	 * occurences at each position is the distance for this position
1323 	 * and the corresponding digit is the maximumn likelihood
1324 	 * candidate. If the distance is zero, assume a miss '_'; if the
1325 	 * distance is not more than half the total number of
1326 	 * occurences, assume a soft error '-'; if two different digits
1327 	 * with the same distance are found, assume a hard error '='.
1328 	 * These will later cause a format error when the timecode is
1329 	 * interpreted. The decoding distance is defined as the minimum
1330 	 * distance over the first nine digits. The tenth digit varies
1331 	 * over the seconds, so we don't count it.
1332 	 */
1333 	up->mindist = 16;
1334 	for (i = 0; i < 9; i++) {
1335 		val1 = val2 = 0;
1336 		k = 0;
1337 		for (j = 0; j < 16; j++) {
1338 			temp = up->decode[i][j] + up->decode[i + 10][j];
1339 			if (temp > val1) {
1340 				val2 = val1;
1341 				val1 = temp;
1342 				k = j;
1343 			}
1344 		}
1345 		if (val1 == 0)
1346 			code[i] = HEX_MISS;
1347 		else if (val1 == val2)
1348 			code[i] = HEX_HARD;
1349 		else if (val1 <= up->burstcnt)
1350 			code[i] = HEX_SOFT;
1351 		else
1352 			code[i] = k;
1353 		if (val1 < up->mindist)
1354 			up->mindist = val1;
1355 		code[i] = hexchar[code[i]];
1356 	}
1357 	code[i] = 0;
1358 
1359 	/*
1360 	 * A valid timecode requires at least three bursts and a
1361 	 * decoding distance greater than half the total number of
1362 	 * occurences. A valid timecode also requires at least 20 valid
1363 	 * timestamps.
1364 	 */
1365 	if (up->burstcnt < 3 || up->mindist <= up->burstcnt)
1366 		up->status |= DECODE;
1367 	if (up->ntstamp < MINSTAMP)
1368 		up->status |= STAMP;
1369 
1370 	/*
1371 	 * Compute the timecode timestamp from the days, hours and
1372 	 * minutes of the timecode. Use clocktime() for the aggregate
1373 	 * minutes and the minute offset computed from the burst
1374 	 * seconds. Note that this code relies on the filesystem time
1375 	 * for the years and does not use the years of the timecode.
1376 	 */
1377 	if (sscanf((char *)code, "%1x%3d%2d%2d", &synchar, &pp->day,
1378 	    &pp->hour, &pp->minute) != 4) {
1379 		up->status |= AFORMAT;
1380 		return (0);
1381 	}
1382 	if (up->status & (DECODE | STAMP)) {
1383 		up->errflg = CEVNT_BADREPLY;
1384 		return (0);
1385 	}
1386 	L_CLR(&offset);
1387 	if (!clocktime(pp->day, pp->hour, pp->minute, 0, GMT,
1388 	    up->tstamp[0].l_ui, &pp->yearstart, &offset.l_ui)) {
1389 		up->errflg = CEVNT_BADTIME;
1390 		return (0);
1391 	}
1392 	pp->lastref = offset;
1393 	pp->variance = 0;
1394 	for (i = 0; i < up->ntstamp; i++) {
1395 		toffset = offset;
1396 		L_SUB(&toffset, &up->tstamp[i]);
1397 		LFPTOD(&toffset, dtemp);
1398 		SAMPLE(dtemp + FUDGE + pp->fudgetime1);
1399 	}
1400 	return (i);
1401 }
1402 
1403 
1404 /*
1405  * chu_clear - clear decoding matrix
1406  */
1407 static void
1408 chu_clear(
1409 	struct peer *peer	/* peer structure pointer */
1410 	)
1411 {
1412 	struct refclockproc *pp;
1413 	struct chuunit *up;
1414 	int	i, j;
1415 
1416 	pp = peer->procptr;
1417 	up = (struct chuunit *)pp->unitptr;
1418 
1419 	/*
1420 	 * Clear stuff for the minute.
1421 	 */
1422 	up->ndx = up->prevsec = 0;
1423 	up->burstcnt = up->mindist = up->ntstamp = 0;
1424 	up->status &= INSYNC | INYEAR;
1425 	up->burstcnt = 0;
1426 	for (i = 0; i < 20; i++) {
1427 		for (j = 0; j < 16; j++)
1428 			up->decode[i][j] = 0;
1429 	}
1430 }
1431 
1432 
1433 /*
1434  * chu_dist - determine the distance of two octet arguments
1435  */
1436 static int
1437 chu_dist(
1438 	int	x,		/* an octet of bits */
1439 	int	y		/* another octet of bits */
1440 	)
1441 {
1442 	int	val;		/* bit count */
1443 	int	temp;
1444 	int	i;
1445 
1446 	/*
1447 	 * The distance is determined as the weight of the exclusive OR
1448 	 * of the two arguments. The weight is determined by the number
1449 	 * of one bits in the result. Each one bit increases the weight,
1450 	 * while each zero bit decreases it.
1451 	 */
1452 	temp = x ^ y;
1453 	val = 0;
1454 	for (i = 0; i < 8; i++) {
1455 		if ((temp & 0x1) == 0)
1456 			val++;
1457 		else
1458 			val--;
1459 		temp >>= 1;
1460 	}
1461 	return (val);
1462 }
1463 
1464 
1465 #ifdef AUDIO_CHU
1466 /*
1467  * chu_gain - adjust codec gain
1468  *
1469  * This routine is called once each second. If the signal envelope
1470  * amplitude is too low, the codec gain is bumped up by four units; if
1471  * too high, it is bumped down. The decoder is relatively insensitive to
1472  * amplitude, so this crudity works just fine. The input port is set and
1473  * the error flag is cleared, mostly to be ornery.
1474  */
1475 static void
1476 chu_gain(
1477 	struct peer *peer	/* peer structure pointer */
1478 	)
1479 {
1480 	struct refclockproc *pp;
1481 	struct chuunit *up;
1482 
1483 	pp = peer->procptr;
1484 	up = (struct chuunit *)pp->unitptr;
1485 
1486 	/*
1487 	 * Apparently, the codec uses only the high order bits of the
1488 	 * gain control field. Thus, it may take awhile for changes to
1489 	 * wiggle the hardware bits.
1490 	 */
1491 	if (up->clipcnt == 0) {
1492 		up->gain += 4;
1493 		if (up->gain > 255)
1494 			up->gain = 255;
1495 	} else if (up->clipcnt > SECOND / 100) {
1496 		up->gain -= 4;
1497 		if (up->gain < 0)
1498 			up->gain = 0;
1499 	}
1500 	audio_gain(up->gain, up->port);
1501 	up->clipcnt = 0;
1502 }
1503 #endif /* AUDIO_CHU */
1504 
1505 
1506 #else
1507 int refclock_chu_bs;
1508 #endif /* REFCLOCK */
1509