2 * refclock_chu - clock driver for Canadian CHU time/frequency station
8 #if defined(REFCLOCK) && defined(CLOCK_CHU)
12 #include "ntp_refclock.h"
13 #include "ntp_calendar.h"
14 #include "ntp_stdlib.h"
22 #endif /* HAVE_AUDIO */
24 #define ICOM 1 /* undefine to suppress ICOM code */
31 * Audio CHU demodulator/decoder
33 * This driver synchronizes the computer time using data encoded in
34 * radio transmissions from Canadian time/frequency station CHU in
35 * Ottawa, Ontario. Transmissions are made continuously on 3330 kHz,
36 * 7335 kHz and 14670 kHz in upper sideband, compatible AM mode. An
37 * ordinary shortwave receiver can be tuned manually to one of these
38 * frequencies or, in the case of ICOM receivers, the receiver can be
39 * tuned automatically using this program as propagation conditions
40 * change throughout the day and night.
42 * The driver receives, demodulates and decodes the radio signals when
43 * connected to the audio codec of a suported workstation hardware and
44 * operating system. These include Solaris, SunOS, FreeBSD, NetBSD and
45 * Linux. In this implementation, only one audio driver and codec can be
46 * supported on a single machine.
48 * The driver can be compiled to use a Bell 103 compatible modem or
49 * modem chip to receive the radio signal and demodulate the data.
50 * Alternatively, the driver can be compiled to use the audio codec of
51 * the Sun workstation or another with compatible audio drivers. In the
52 * latter case, the driver implements the modem using DSP routines, so
53 * the radio can be connected directly to either the microphone on line
54 * input port. In either case, the driver decodes the data using a
55 * maximum likelihood technique which exploits the considerable degree
56 * of redundancy available to maximize accuracy and minimize errors.
58 * The CHU time broadcast includes an audio signal compatible with the
59 * Bell 103 modem standard (mark = 2225 Hz, space = 2025 Hz). It consist
60 * of nine, ten-character bursts transmitted at 300 bps and beginning
61 * each second from second 31 to second 39 of the minute. Each character
62 * consists of eight data bits plus one start bit and two stop bits to
63 * encode two hex digits. The burst data consist of five characters (ten
64 * hex digits) followed by a repeat of these characters. In format A,
65 * the characters are repeated in the same polarity; in format B, the
66 * characters are repeated in the opposite polarity.
68 * Format A bursts are sent at seconds 32 through 39 of the minute in
71 * 6dddhhmmss6dddhhmmss
73 * The first ten digits encode a frame marker (6) followed by the day
74 * (ddd), hour (hh in UTC), minute (mm) and the second (ss). Since
75 * format A bursts are sent during the third decade of seconds the tens
76 * digit of ss is always 3. The driver uses this to determine correct
77 * burst synchronization. These digits are then repeated with the same
80 * Format B bursts are sent at second 31 of the minute in hex digits
82 * xdyyyyttaaxdyyyyttaa
84 * The first ten digits encode a code (x described below) followed by
85 * the DUT1 (d in deciseconds), Gregorian year (yyyy), difference TAI -
86 * UTC (tt) and daylight time indicator (aa) peculiar to Canada. These
87 * digits are then repeated with inverted polarity.
91 * 1 Sign of DUT (0 = +)
92 * 2 Leap second warning. One second will be added.
93 * 4 Leap second warning. One second will be subtracted.
94 * 8 Even parity bit for this nibble.
96 * By design, the last stop bit of the last character in the burst
97 * coincides with 0.5 second. Since characters have 11 bits and are
98 * transmitted at 300 bps, the last stop bit of the first character
99 * coincides with 0.5 - 10 * 11/300 = 0.133 second. Depending on the
100 * UART, character interrupts can vary somewhere between the beginning
101 * of bit 9 and end of bit 11. These eccentricities can be corrected
102 * along with the radio propagation delay using fudge time 1.
106 * The timecode format used for debugging and data recording includes
107 * data helpful in diagnosing problems with the radio signal and serial
108 * connections. With debugging enabled (-d on the ntpd command line),
109 * the driver produces one line for each burst in two formats
110 * corresponding to format A and B. Following is format A:
114 * where n is the number of characters in the burst (0-11), b the burst
115 * distance (0-40), f the field alignment (-1, 0, 1), s the
116 * synchronization distance (0-16), m the burst number (2-9) and code
117 * the burst characters as received. Note that the hex digits in each
118 * character are reversed, so the burst
120 * 10 38 0 16 9 06851292930685129293
122 * is interpreted as containing 11 characters with burst distance 38,
123 * field alignment 0, synchronization distance 16 and burst number 9.
124 * The nibble-swapped timecode shows day 58, hour 21, minute 29 and
127 * When the audio driver is compiled, format A is preceded by
128 * the current gain (0-255) and relative signal level (0-9999). The
129 * receiver folume control should be set so that the gain is somewhere
130 * near the middle of the range 0-255, which results in a signal level
133 * Following is format B:
137 * where n is the number of characters in the burst (0-11), b the burst
138 * distance (0-40), s the synchronization distance (0-40) and code the
139 * burst characters as received. Note that the hex digits in each
140 * character are reversed and the last ten digits inverted, so the burst
142 * 11 40 1091891300ef6e76ecff
144 * is interpreted as containing 11 characters with burst distance 40.
145 * The nibble-swapped timecode shows DUT1 +0.1 second, year 1998 and TAI
148 * In addition to the above, the reference timecode is updated and
149 * written to the clockstats file and debug score after the last burst
150 * received in the minute. The format is
152 * qq yyyy ddd hh:mm:ss nn dd tt
154 * where qq are the error flags, as described below, yyyy is the year,
155 * ddd the day, hh:mm:ss the time of day, nn the number of format A
156 * bursts received during the previous minute, dd the decoding distance
157 * and tt the number of timestamps. The error flags are cleared after
162 * For accuracies better than the low millisceconds, fudge time1 can be
163 * set to the radio propagation delay from CHU to the receiver. This can
164 * be done conviently using the minimuf program.
166 * Fudge flag4 causes the dubugging output described above to be
167 * recorded in the clockstats file. When the audio driver is compiled,
168 * fudge flag2 selects the audio input port, where 0 is the mike port
169 * (default) and 1 is the line-in port. It does not seem useful to
170 * select the compact disc player port. Fudge flag3 enables audio
171 * monitoring of the input signal. For this purpose, the monitor gain is
172 * set to a default value.
174 * The audio codec code is normally compiled in the driver if the
175 * architecture supports it (HAVE_AUDIO defined), but is used only if
176 * the link /dev/chu_audio is defined and valid. The serial port code is
177 * always compiled in the driver, but is used only if the autdio codec
178 * is not available and the link /dev/chu%d is defined and valid.
180 * The ICOM code is normally compiled in the driver if selected (ICOM
181 * defined), but is used only if the link /dev/icom%d is defined and
182 * valid and the mode keyword on the server configuration command
183 * specifies a nonzero mode (ICOM ID select code). The C-IV speed is
184 * 9600 bps if the high order 0x80 bit of the mode is zero and 1200 bps
185 * if one. The C-IV trace is turned on if the debug level is greater
189 * Interface definitions
191 #define SPEED232 B300 /* uart speed (300 baud) */
192 #define PRECISION (-10) /* precision assumed (about 1 ms) */
193 #define REFID "CHU" /* reference ID */
194 #define DEVICE "/dev/chu%d" /* device name and unit */
195 #define SPEED232 B300 /* UART speed (300 baud) */
197 #define TUNE .001 /* offset for narrow filter (kHz) */
198 #define DWELL 5 /* minutes in a probe cycle */
199 #define NCHAN 3 /* number of channels */
200 #define ISTAGE 3 /* number of integrator stages */
205 * Audio demodulator definitions
207 #define SECOND 8000 /* nominal sample rate (Hz) */
208 #define BAUD 300 /* modulation rate (bps) */
209 #define OFFSET 128 /* companded sample offset */
210 #define SIZE 256 /* decompanding table size */
211 #define MAXAMP 6000. /* maximum signal level */
212 #define MAXCLP 100 /* max clips above reference per s */
213 #define LIMIT 1000. /* soft limiter threshold */
214 #define AGAIN 6. /* baseband gain */
215 #define LAG 10 /* discriminator lag */
216 #define DEVICE_AUDIO "/dev/audio" /* device name */
217 #define DESCRIPTION "CHU Audio/Modem Receiver" /* WRU */
218 #define AUDIO_BUFSIZ 240 /* audio buffer size (30 ms) */
220 #define DESCRIPTION "CHU Modem Receiver" /* WRU */
221 #endif /* HAVE_AUDIO */
224 * Decoder definitions
226 #define CHAR (11. / 300.) /* character time (s) */
227 #define FUDGE .185 /* offset to first stop bit (s) */
228 #define BURST 11 /* max characters per burst */
229 #define MINCHAR 9 /* min characters per burst */
230 #define MINDIST 28 /* min burst distance (of 40) */
231 #define MINBURST 4 /* min bursts in minute */
232 #define MINSYNC 8 /* min sync distance (of 16) */
233 #define MINSTAMP 20 /* min timestamps (of 60) */
234 #define METRIC 50. /* min channel metric */
235 #define PANIC 1440 /* panic timeout (m) */
236 #define HOLD 30 /* reach hold (m) */
239 * Hex extension codes (>= 16)
241 #define HEX_MISS 16 /* miss _ */
242 #define HEX_SOFT 17 /* soft error * */
243 #define HEX_HARD 18 /* hard error = */
246 * Status bits (status)
248 #define RUNT 0x0001 /* runt burst */
249 #define NOISE 0x0002 /* noise burst */
250 #define BFRAME 0x0004 /* invalid format B frame sync */
251 #define BFORMAT 0x0008 /* invalid format B data */
252 #define AFRAME 0x0010 /* invalid format A frame sync */
253 #define AFORMAT 0x0020 /* invalid format A data */
254 #define DECODE 0x0040 /* invalid data decode */
255 #define STAMP 0x0080 /* too few timestamps */
256 #define AVALID 0x0100 /* valid A frame */
257 #define BVALID 0x0200 /* valid B frame */
258 #define INSYNC 0x0400 /* clock synchronized */
261 * Alarm status bits (alarm)
263 * These alarms are set at the end of a minute in which at least one
264 * burst was received. SYNERR is raised if the AFRAME or BFRAME status
265 * bits are set during the minute, FMTERR is raised if the AFORMAT or
266 * BFORMAT status bits are set, DECERR is raised if the DECODE status
267 * bit is set and TSPERR is raised if the STAMP status bit is set.
269 #define SYNERR 0x01 /* frame sync error */
270 #define FMTERR 0x02 /* data format error */
271 #define DECERR 0x04 /* data decoding error */
272 #define TSPERR 0x08 /* insufficient data */
276 * Maximum likelihood UART structure. There are eight of these
277 * corresponding to the number of phases.
280 double shift[12]; /* mark register */
281 double es_max, es_min; /* max/min envelope signals */
282 double dist; /* sample distance */
283 int uart; /* decoded character */
285 #endif /* HAVE_AUDIO */
289 * CHU station structure. There are three of these corresponding to the
293 double integ[ISTAGE]; /* circular integrator */
294 double metric; /* integrator sum */
295 int iptr; /* integrator pointer */
296 int probe; /* dwells since last probe */
301 * CHU unit control structure
304 u_char decode[20][16]; /* maximum likelihood decoding matrix */
305 l_fp cstamp[BURST]; /* character timestamps */
306 l_fp tstamp[MAXSTAGE]; /* timestamp samples */
307 l_fp timestamp; /* current buffer timestamp */
308 l_fp laststamp; /* last buffer timestamp */
309 l_fp charstamp; /* character time as a l_fp */
310 int errflg; /* error flags */
311 int status; /* status bits */
312 char ident[5]; /* station ID and channel */
314 int fd_icom; /* ICOM file descriptor */
315 int chan; /* data channel */
316 int achan; /* active channel */
317 int dwell; /* dwell cycle */
318 struct xmtr xmtr[NCHAN]; /* station metric */
322 * Character burst variables
324 int cbuf[BURST]; /* character buffer */
325 int ntstamp; /* number of timestamp samples */
326 int ndx; /* buffer start index */
327 int prevsec; /* previous burst second */
328 int burdist; /* burst distance */
329 int syndist; /* sync distance */
330 int burstcnt; /* format A bursts this minute */
335 int leap; /* leap/dut code */
336 int dut; /* UTC1 correction */
337 int tai; /* TAI - UTC correction */
338 int dst; /* Canadian DST code */
342 * Audio codec variables
344 int fd_audio; /* audio port file descriptor */
345 double comp[SIZE]; /* decompanding table */
346 int port; /* codec port */
347 int gain; /* codec gain */
348 int mongain; /* codec monitor gain */
349 int clipcnt; /* sample clip count */
350 int seccnt; /* second interval counter */
355 l_fp tick; /* audio sample increment */
356 double bpf[9]; /* IIR bandpass filter */
357 double disc[LAG]; /* discriminator shift register */
358 double lpf[27]; /* FIR lowpass filter */
359 double monitor; /* audio monitor */
360 double maxsignal; /* signal level */
361 int discptr; /* discriminator pointer */
364 * Maximum likelihood UART variables
366 double baud; /* baud interval */
367 struct surv surv[8]; /* UART survivor structures */
368 int decptr; /* decode pointer */
369 int dbrk; /* holdoff counter */
370 #endif /* HAVE_AUDIO */
374 * Function prototypes
376 static int chu_start P((int, struct peer *));
377 static void chu_shutdown P((int, struct peer *));
378 static void chu_receive P((struct recvbuf *));
379 static void chu_poll P((int, struct peer *));
382 * More function prototypes
384 static void chu_decode P((struct peer *, int));
385 static void chu_burst P((struct peer *));
386 static void chu_clear P((struct peer *));
387 static void chu_a P((struct peer *, int));
388 static void chu_b P((struct peer *, int));
389 static int chu_dist P((int, int));
390 static double chu_major P((struct peer *));
392 static void chu_uart P((struct surv *, double));
393 static void chu_rf P((struct peer *, double));
394 static void chu_gain P((struct peer *));
395 static void chu_audio_receive P((struct recvbuf *rbufp));
396 #endif /* HAVE_AUDIO */
398 static int chu_newchan P((struct peer *, double));
400 static void chu_serial_receive P((struct recvbuf *rbufp));
405 static char hexchar[] = "0123456789abcdef_*=";
409 * Note the tuned frequencies are 1 kHz higher than the carrier. CHU
410 * transmits on USB with carrier so we can use AM and the narrow SSB
413 static double qsy[NCHAN] = {3.330, 7.335, 14.670}; /* freq (MHz) */
419 struct refclock refclock_chu = {
420 chu_start, /* start up driver */
421 chu_shutdown, /* shut down driver */
422 chu_poll, /* transmit poll message */
423 noentry, /* not used (old chu_control) */
424 noentry, /* initialize driver (not used) */
425 noentry, /* not used (old chu_buginfo) */
426 NOFLAGS /* not used */
431 * chu_start - open the devices and initialize data for processing
435 int unit, /* instance number (not used) */
436 struct peer *peer /* peer structure pointer */
440 struct refclockproc *pp;
441 char device[20]; /* device name */
442 int fd; /* file descriptor */
447 int fd_audio; /* audio port file descriptor */
449 double step; /* codec adjustment */
454 fd_audio = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit);
456 if (fd_audio > 0 && debug)
461 * Open serial port in raw mode.
466 sprintf(device, DEVICE, unit);
467 fd = refclock_open(device, SPEED232, LDISC_RAW);
469 #else /* HAVE_AUDIO */
472 * Open serial port in raw mode.
474 sprintf(device, DEVICE, unit);
475 fd = refclock_open(device, SPEED232, LDISC_RAW);
476 #endif /* HAVE_AUDIO */
481 * Allocate and initialize unit structure
483 if (!(up = (struct chuunit *)
484 emalloc(sizeof(struct chuunit)))) {
488 memset((char *)up, 0, sizeof(struct chuunit));
490 pp->unitptr = (caddr_t)up;
491 pp->io.clock_recv = chu_receive;
492 pp->io.srcclock = (caddr_t)peer;
495 if (!io_addclock(&pp->io)) {
502 * Initialize miscellaneous variables
504 peer->precision = PRECISION;
505 pp->clockdesc = DESCRIPTION;
506 strcpy(up->ident, "CHU");
507 memcpy(&peer->refid, up->ident, 4);
508 DTOLFP(CHAR, &up->charstamp);
512 * The companded samples are encoded sign-magnitude. The table
513 * contains all the 256 values in the interest of speed. We do
514 * this even if the audio codec is not available. C'est la lazy.
516 up->fd_audio = fd_audio;
518 up->comp[0] = up->comp[OFFSET] = 0.;
519 up->comp[1] = 1; up->comp[OFFSET + 1] = -1.;
520 up->comp[2] = 3; up->comp[OFFSET + 2] = -3.;
522 for (i = 3; i < OFFSET; i++) {
523 up->comp[i] = up->comp[i - 1] + step;
524 up->comp[OFFSET + i] = -up->comp[i];
528 DTOLFP(1. / SECOND, &up->tick);
529 #endif /* HAVE_AUDIO */
537 if (peer->ttl & 0x80)
538 up->fd_icom = icom_init("/dev/icom", B1200,
541 up->fd_icom = icom_init("/dev/icom", B9600,
544 if (up->fd_icom > 0) {
545 if (chu_newchan(peer, 0) != 0) {
546 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
548 "icom: radio not found");
549 up->errflg = CEVNT_FAULT;
553 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT)
555 "icom: autotune enabled");
564 * chu_shutdown - shut down the clock
568 int unit, /* instance number (not used) */
569 struct peer *peer /* peer structure pointer */
573 struct refclockproc *pp;
576 up = (struct chuunit *)pp->unitptr;
580 io_closeclock(&pp->io);
590 * chu_receive - receive data from the audio or serial device
594 struct recvbuf *rbufp /* receive buffer structure pointer */
599 struct refclockproc *pp;
602 peer = (struct peer *)rbufp->recv_srcclock;
604 up = (struct chuunit *)pp->unitptr;
607 * If the audio codec is warmed up, the buffer contains codec
608 * samples which need to be demodulated and decoded into CHU
609 * characters using the software UART. Otherwise, the buffer
610 * contains CHU characters from the serial port, so the software
611 * UART is bypassed. In this case the CPU will probably run a
612 * few degrees cooler.
614 if (up->fd_audio > 0)
615 chu_audio_receive(rbufp);
617 chu_serial_receive(rbufp);
619 chu_serial_receive(rbufp);
620 #endif /* HAVE_AUDIO */
626 * chu_audio_receive - receive data from the audio device
630 struct recvbuf *rbufp /* receive buffer structure pointer */
634 struct refclockproc *pp;
637 double sample; /* codec sample */
638 u_char *dpt; /* buffer pointer */
639 int bufcnt; /* buffer counter */
640 l_fp ltemp; /* l_fp temp */
642 peer = (struct peer *)rbufp->recv_srcclock;
644 up = (struct chuunit *)pp->unitptr;
647 * Main loop - read until there ain't no more. Note codec
648 * samples are bit-inverted.
650 DTOLFP((double)rbufp->recv_length / SECOND, <emp);
651 L_SUB(&rbufp->recv_time, <emp);
652 up->timestamp = rbufp->recv_time;
653 dpt = rbufp->recv_buffer;
654 for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
655 sample = up->comp[~*dpt++ & 0xff];
658 * Clip noise spikes greater than MAXAMP. If no clips,
659 * increase the gain a tad; if the clips are too high,
662 if (sample > MAXAMP) {
665 } else if (sample < -MAXAMP) {
669 chu_rf(peer, sample);
670 L_ADD(&up->timestamp, &up->tick);
673 * Once each second ride gain.
675 up->seccnt = (up->seccnt + 1) % SECOND;
676 if (up->seccnt == 0) {
677 pp->second = (pp->second + 1) % 60;
683 * Set the input port and monitor gain for the next buffer.
685 if (pp->sloppyclockflag & CLK_FLAG2)
689 if (pp->sloppyclockflag & CLK_FLAG3)
690 up->mongain = MONGAIN;
697 * chu_rf - filter and demodulate the FSK signal
699 * This routine implements a 300-baud Bell 103 modem with mark 2225 Hz
700 * and space 2025 Hz. It uses a bandpass filter followed by a soft
701 * limiter, FM discriminator and lowpass filter. A maximum likelihood
702 * decoder samples the baseband signal at eight times the baud rate and
703 * detects the start bit of each character.
705 * The filters are built for speed, which explains the rather clumsy
706 * code. Hopefully, the compiler will efficiently implement the move-
707 * and-muiltiply-and-add operations.
711 struct peer *peer, /* peer structure pointer */
712 double sample /* analog sample */
715 struct refclockproc *pp;
722 double signal; /* bandpass signal */
723 double limit; /* limiter signal */
724 double disc; /* discriminator signal */
725 double lpf; /* lowpass signal */
726 double span; /* UART signal span */
727 double dist; /* UART signal distance */
731 up = (struct chuunit *)pp->unitptr;
734 * Bandpass filter. 4th-order elliptic, 500-Hz bandpass centered
735 * at 2125 Hz. Passband ripple 0.3 dB, stopband ripple 50 dB.
737 signal = (up->bpf[8] = up->bpf[7]) * 5.844676e-01;
738 signal += (up->bpf[7] = up->bpf[6]) * 4.884860e-01;
739 signal += (up->bpf[6] = up->bpf[5]) * 2.704384e+00;
740 signal += (up->bpf[5] = up->bpf[4]) * 1.645032e+00;
741 signal += (up->bpf[4] = up->bpf[3]) * 4.644557e+00;
742 signal += (up->bpf[3] = up->bpf[2]) * 1.879165e+00;
743 signal += (up->bpf[2] = up->bpf[1]) * 3.522634e+00;
744 signal += (up->bpf[1] = up->bpf[0]) * 7.315738e-01;
745 up->bpf[0] = sample - signal;
746 signal = up->bpf[0] * 6.176213e-03
747 + up->bpf[1] * 3.156599e-03
748 + up->bpf[2] * 7.567487e-03
749 + up->bpf[3] * 4.344580e-03
750 + up->bpf[4] * 1.190128e-02
751 + up->bpf[5] * 4.344580e-03
752 + up->bpf[6] * 7.567487e-03
753 + up->bpf[7] * 3.156599e-03
754 + up->bpf[8] * 6.176213e-03;
756 up->monitor = signal / 4.; /* note monitor after filter */
759 * Soft limiter/discriminator. The 11-sample discriminator lag
760 * interval corresponds to three cycles of 2125 Hz, which
761 * requires the sample frequency to be 2125 * 11 / 3 = 7791.7
762 * Hz. The discriminator output varies +-0.5 interval for input
763 * frequency 2025-2225 Hz. However, we don't get to sample at
764 * this frequency, so the discriminator output is biased. Life
770 else if (limit < -LIMIT)
772 disc = up->disc[up->discptr] * -limit;
773 up->disc[up->discptr] = limit;
774 up->discptr = (up->discptr + 1 ) % LAG;
781 * Lowpass filter. Raised cosine, Ts = 1 / 300, beta = 0.1.
783 lpf = (up->lpf[26] = up->lpf[25]) * 2.538771e-02;
784 lpf += (up->lpf[25] = up->lpf[24]) * 1.084671e-01;
785 lpf += (up->lpf[24] = up->lpf[23]) * 2.003159e-01;
786 lpf += (up->lpf[23] = up->lpf[22]) * 2.985303e-01;
787 lpf += (up->lpf[22] = up->lpf[21]) * 4.003697e-01;
788 lpf += (up->lpf[21] = up->lpf[20]) * 5.028552e-01;
789 lpf += (up->lpf[20] = up->lpf[19]) * 6.028795e-01;
790 lpf += (up->lpf[19] = up->lpf[18]) * 6.973249e-01;
791 lpf += (up->lpf[18] = up->lpf[17]) * 7.831828e-01;
792 lpf += (up->lpf[17] = up->lpf[16]) * 8.576717e-01;
793 lpf += (up->lpf[16] = up->lpf[15]) * 9.183463e-01;
794 lpf += (up->lpf[15] = up->lpf[14]) * 9.631951e-01;
795 lpf += (up->lpf[14] = up->lpf[13]) * 9.907208e-01;
796 lpf += (up->lpf[13] = up->lpf[12]) * 1.000000e+00;
797 lpf += (up->lpf[12] = up->lpf[11]) * 9.907208e-01;
798 lpf += (up->lpf[11] = up->lpf[10]) * 9.631951e-01;
799 lpf += (up->lpf[10] = up->lpf[9]) * 9.183463e-01;
800 lpf += (up->lpf[9] = up->lpf[8]) * 8.576717e-01;
801 lpf += (up->lpf[8] = up->lpf[7]) * 7.831828e-01;
802 lpf += (up->lpf[7] = up->lpf[6]) * 6.973249e-01;
803 lpf += (up->lpf[6] = up->lpf[5]) * 6.028795e-01;
804 lpf += (up->lpf[5] = up->lpf[4]) * 5.028552e-01;
805 lpf += (up->lpf[4] = up->lpf[3]) * 4.003697e-01;
806 lpf += (up->lpf[3] = up->lpf[2]) * 2.985303e-01;
807 lpf += (up->lpf[2] = up->lpf[1]) * 2.003159e-01;
808 lpf += (up->lpf[1] = up->lpf[0]) * 1.084671e-01;
809 lpf += up->lpf[0] = disc * 2.538771e-02;
812 * Maximum likelihood decoder. The UART updates each of the
813 * eight survivors and determines the span, slice level and
814 * tentative decoded character. Valid 11-bit characters are
815 * framed so that bit 1 and bit 11 (stop bits) are mark and bit
816 * 2 (start bit) is space. When a valid character is found, the
817 * survivor with maximum distance determines the final decoded
820 up->baud += 1. / SECOND;
821 if (up->baud > 1. / (BAUD * 8.)) {
822 up->baud -= 1. / (BAUD * 8.);
823 sp = &up->surv[up->decptr];
824 span = sp->es_max - sp->es_min;
825 up->maxsignal += (span - up->maxsignal) / 80.;
828 } else if ((sp->uart & 0x403) == 0x401 && span > 1000.)
832 for (i = 0; i < 8; i++) {
833 if (up->surv[i].dist > dist) {
834 dist = up->surv[i].dist;
838 chu_decode(peer, (up->surv[j].uart >> 2) &
842 up->decptr = (up->decptr + 1) % 8;
843 chu_uart(sp, -lpf * AGAIN);
849 * chu_uart - maximum likelihood UART
851 * This routine updates a shift register holding the last 11 envelope
852 * samples. It then computes the slice level and span over these samples
853 * and determines the tentative data bits and distance. The calling
854 * program selects over the last eight survivors the one with maximum
855 * distance to determine the decoded character.
859 struct surv *sp, /* survivor structure pointer */
860 double sample /* baseband signal */
863 double es_max, es_min; /* max/min envelope */
864 double slice; /* slice level */
865 double dist; /* distance */
870 * Save the sample and shift right. At the same time, measure
871 * the maximum and minimum over all eleven samples.
875 sp->shift[0] = sample;
876 for (i = 11; i > 0; i--) {
877 sp->shift[i] = sp->shift[i - 1];
878 if (sp->shift[i] > es_max)
879 es_max = sp->shift[i];
880 if (sp->shift[i] < es_min)
881 es_min = sp->shift[i];
885 * Determine the slice level midway beteen the maximum and
886 * minimum and the span as the maximum less the minimum. Compute
887 * the distance on the assumption the first and last bits must
888 * be mark, the second space and the rest either mark or space.
890 slice = (es_max + es_min) / 2.;
893 for (i = 1; i < 12; i++) {
895 dtemp = sp->shift[i];
898 if (i == 1 || i == 11) {
899 dist += dtemp - es_min;
900 } else if (i == 10) {
901 dist += es_max - dtemp;
904 dist += dtemp - es_min;
906 dist += es_max - dtemp;
911 sp->dist = dist / (11 * (es_max - es_min));
913 #endif /* HAVE_AUDIO */
917 * chu_serial_receive - receive data from the serial device
921 struct recvbuf *rbufp /* receive buffer structure pointer */
925 struct refclockproc *pp;
928 u_char *dpt; /* receive buffer pointer */
930 peer = (struct peer *)rbufp->recv_srcclock;
932 up = (struct chuunit *)pp->unitptr;
935 * Initialize pointers and read the timecode and timestamp.
937 up->timestamp = rbufp->recv_time;
938 dpt = (u_char *)&rbufp->recv_space;
939 chu_decode(peer, *dpt);
944 * chu_decode - decode the character data
948 struct peer *peer, /* peer structure pointer */
949 int hexhex /* data character */
952 struct refclockproc *pp;
955 l_fp tstmp; /* timestamp temp */
959 up = (struct chuunit *)pp->unitptr;
962 * If the interval since the last character is greater than the
963 * longest burst, process the last burst and start a new one. If
964 * the interval is less than this but greater than two
965 * characters, consider this a noise burst and reject it.
967 tstmp = up->timestamp;
968 if (L_ISZERO(&up->laststamp))
969 up->laststamp = up->timestamp;
970 L_SUB(&tstmp, &up->laststamp);
971 up->laststamp = up->timestamp;
972 LFPTOD(&tstmp, dtemp);
973 if (dtemp > BURST * CHAR) {
976 } else if (dtemp > 2.5 * CHAR) {
981 * Append the character to the current burst and append the
982 * timestamp to the timestamp list.
984 if (up->ndx < BURST) {
985 up->cbuf[up->ndx] = hexhex & 0xff;
986 up->cstamp[up->ndx] = up->timestamp;
994 * chu_burst - search for valid burst format
1002 struct refclockproc *pp;
1007 up = (struct chuunit *)pp->unitptr;
1010 * Correlate a block of five characters with the next block of
1011 * five characters. The burst distance is defined as the number
1012 * of bits that match in the two blocks for format A and that
1013 * match the inverse for format B.
1015 if (up->ndx < MINCHAR) {
1020 for (i = 0; i < 5 && i < up->ndx - 5; i++)
1021 up->burdist += chu_dist(up->cbuf[i], up->cbuf[i + 5]);
1024 * If the burst distance is at least MINDIST, this must be a
1025 * format A burst; if the value is not greater than -MINDIST, it
1026 * must be a format B burst. If the B burst is perfect, we
1027 * believe it; otherwise, it is a noise burst and of no use to
1030 if (up->burdist >= MINDIST) {
1031 chu_a(peer, up->ndx);
1032 } else if (up->burdist <= -MINDIST) {
1033 chu_b(peer, up->ndx);
1035 up->status |= NOISE;
1040 * If this is a valid burst, wait a guard time of ten seconds to
1041 * allow for more bursts, then arm the poll update routine to
1042 * process the minute. Don't do this if this is called from the
1043 * timer interrupt routine.
1045 if (peer->outdate != current_time)
1046 peer->nextdate = current_time + 10;
1051 * chu_b - decode format B burst
1059 struct refclockproc *pp;
1062 u_char code[11]; /* decoded timecode */
1063 char tbuf[80]; /* trace buffer */
1064 l_fp offset; /* timestamp offset */
1068 up = (struct chuunit *)pp->unitptr;
1071 * In a format B burst, a character is considered valid only if
1072 * the first occurrence matches the last occurrence. The burst
1073 * is considered valid only if all characters are valid; that
1074 * is, only if the distance is 40. Note that once a valid frame
1075 * has been found errors are ignored.
1077 sprintf(tbuf, "chuB %04x %2d %2d ", up->status, nchar,
1079 for (i = 0; i < nchar; i++)
1080 sprintf(&tbuf[strlen(tbuf)], "%02x", up->cbuf[i]);
1081 if (pp->sloppyclockflag & CLK_FLAG4)
1082 record_clock_stats(&peer->srcadr, tbuf);
1085 printf("%s\n", tbuf);
1087 if (up->burdist > -40) {
1088 up->status |= BFRAME;
1091 up->status |= BVALID;
1094 * Convert the burst data to internal format. If this succeeds,
1095 * save the timestamps for later.
1097 for (i = 0; i < 5; i++) {
1098 code[2 * i] = hexchar[up->cbuf[i] & 0xf];
1099 code[2 * i + 1] = hexchar[(up->cbuf[i] >>
1102 if (sscanf((char *)code, "%1x%1d%4d%2d%2x", &up->leap, &up->dut,
1103 &pp->year, &up->tai, &up->dst) != 5) {
1104 up->status |= BFORMAT;
1111 for (i = 0; i < nchar && i < 10; i++) {
1112 up->tstamp[up->ntstamp] = up->cstamp[i];
1113 L_SUB(&up->tstamp[up->ntstamp], &offset);
1114 L_ADD(&offset, &up->charstamp);
1115 if (up->ntstamp < MAXSTAGE - 1)
1122 * chu_a - decode format A burst
1130 struct refclockproc *pp;
1133 char tbuf[80]; /* trace buffer */
1134 l_fp offset; /* timestamp offset */
1135 int val; /* distance */
1140 up = (struct chuunit *)pp->unitptr;
1143 * Determine correct burst phase. There are three cases
1144 * corresponding to in-phase, one character early or one
1145 * character late. These cases are distinguished by the position
1146 * of the framing digits x6 at positions 0 and 5 and x3 at
1147 * positions 4 and 9. The correct phase is when the distance
1148 * relative to the framing digits is maximum. The burst is valid
1149 * only if the maximum distance is at least MINSYNC.
1151 up->syndist = k = 0;
1153 for (i = -1; i < 2; i++) {
1154 temp = up->cbuf[i + 4] & 0xf;
1156 temp |= (up->cbuf[i] & 0xf) << 4;
1157 val = chu_dist(temp, 0x63);
1158 temp = (up->cbuf[i + 5] & 0xf) << 4;
1160 temp |= up->cbuf[i + 9] & 0xf;
1161 val += chu_dist(temp, 0x63);
1162 if (val > up->syndist) {
1167 temp = (up->cbuf[k + 4] >> 4) & 0xf;
1168 if (temp > 9 || k + 9 >= nchar || temp != ((up->cbuf[k + 9] >>
1173 sprintf(tbuf, "chuA %04x %4.0f %2d %2d %2d %2d %1d ",
1174 up->status, up->maxsignal, nchar, up->burdist, k,
1177 sprintf(tbuf, "chuA %04x %2d %2d %2d %2d %1d ",
1178 up->status, nchar, up->burdist, k, up->syndist,
1182 sprintf(tbuf, "chuA %04x %2d %2d %2d %2d %1d ", up->status,
1183 nchar, up->burdist, k, up->syndist, temp);
1184 #endif /* HAVE_AUDIO */
1185 for (i = 0; i < nchar; i++)
1186 sprintf(&tbuf[strlen(tbuf)], "%02x",
1188 if (pp->sloppyclockflag & CLK_FLAG4)
1189 record_clock_stats(&peer->srcadr, tbuf);
1192 printf("%s\n", tbuf);
1194 if (up->syndist < MINSYNC) {
1195 up->status |= AFRAME;
1200 * A valid burst requires the first seconds number to match the
1201 * last seconds number. If so, the burst timestamps are
1202 * corrected to the current minute and saved for later
1203 * processing. In addition, the seconds decode is advanced from
1204 * the previous burst to the current one.
1207 pp->second = 30 + temp;
1208 offset.l_ui = 30 + temp;
1212 offset = up->charstamp;
1215 for (; i < nchar && i < k + 10; i++) {
1216 up->tstamp[up->ntstamp] = up->cstamp[i];
1217 L_SUB(&up->tstamp[up->ntstamp], &offset);
1218 L_ADD(&offset, &up->charstamp);
1219 if (up->ntstamp < MAXSTAGE - 1)
1222 while (temp > up->prevsec) {
1223 for (j = 15; j > 0; j--) {
1224 up->decode[9][j] = up->decode[9][j - 1];
1226 up->decode[19][j - 1];
1228 up->decode[9][j] = up->decode[19][j] = 0;
1233 for (j = 0; j < nchar; j++) {
1234 if (i < 0 || i > 18) {
1238 up->decode[i][up->cbuf[j] & 0xf]++;
1240 up->decode[i][(up->cbuf[j] >> 4) & 0xf]++;
1243 up->status |= AVALID;
1249 * chu_poll - called by the transmit procedure
1254 struct peer *peer /* peer structure pointer */
1257 struct refclockproc *pp;
1260 char synchar, qual, leapchar;
1265 up = (struct chuunit *)pp->unitptr;
1266 if (pp->coderecv == pp->codeproc)
1267 up->errflg = CEVNT_TIMEOUT;
1272 * If once in sync and the radio has not been heard for awhile
1273 * (30 m), it is no longer reachable. If not heard in a long
1274 * while (one day), turn out the lights and start from scratch.
1276 minset = ((current_time - peer->update) + 30) / 60;
1277 if (up->status & INSYNC) {
1280 else if (minset <= HOLD)
1285 * Process the last burst, if still in the burst buffer.
1286 * Don't mess with anything if nothing has been heard. If the
1287 * minute contains a valid A frame and valid B frame, assume
1288 * synchronized; however, believe the time only if within metric
1289 * threshold. Note the quality indicator is only for
1290 * diagnostics; the data are used only if in sync and above
1294 if (up->burstcnt == 0) {
1296 chu_newchan(peer, 0);
1300 dtemp = chu_major(peer);
1302 if (up->status & (BFRAME | AFRAME))
1304 if (up->status & (BFORMAT | AFORMAT))
1306 if (up->status & DECODE)
1308 if (up->status & STAMP)
1310 if (up->status & AVALID && up->status & BVALID)
1311 up->status |= INSYNC;
1312 synchar = leapchar = ' ';
1313 if (!(up->status & INSYNC)) {
1314 pp->leap = LEAP_NOTINSYNC;
1316 } else if (up->leap & 0x2) {
1317 pp->leap = LEAP_ADDSECOND;
1319 } else if (up->leap & 0x4) {
1320 pp->leap = LEAP_DELSECOND;
1323 pp->leap = LEAP_NOWARNING;
1327 sprintf(pp->a_lastcode,
1328 "%c%1X %04d %3d %02d:%02d:%02d %c%x %+d %d %d %s %.0f %d",
1329 synchar, qual, pp->year, pp->day, pp->hour,
1330 pp->minute, pp->second, leapchar, up->dst, up->dut,
1331 minset, up->gain, up->ident, dtemp, up->ntstamp);
1333 sprintf(pp->a_lastcode,
1334 "%c%1X %04d %3d %02d:%02d:%02d %c%x %+d %d %s %.0f %d",
1335 synchar, qual, pp->year, pp->day, pp->hour,
1336 pp->minute, pp->second, leapchar, up->dst, up->dut,
1337 minset, up->ident, dtemp, up->ntstamp);
1339 sprintf(pp->a_lastcode,
1340 "%c%1X %04d %3d %02d:%02d:%02d %c%x %+d %d %s %.0f %d",
1341 synchar, qual, pp->year, pp->day, pp->hour, pp->minute,
1342 pp->second, leapchar, up->dst, up->dut, minset, up->ident,
1343 dtemp, up->ntstamp);
1344 #endif /* HAVE_AUDIO */
1345 pp->lencode = strlen(pp->a_lastcode);
1348 * If in sync and the signal metric is above threshold, the
1349 * timecode is ipso fatso valid and can be selected to
1350 * discipline the clock. Be sure not to leave stray timestamps
1351 * around if signals are too weak or the clock time is invalid.
1353 if (up->status & INSYNC && dtemp > METRIC) {
1354 if (!clocktime(pp->day, pp->hour, pp->minute, 0, GMT,
1355 up->tstamp[0].l_ui, &pp->yearstart, &offset.l_ui)) {
1356 up->errflg = CEVNT_BADTIME;
1359 for (i = 0; i < up->ntstamp; i++)
1360 refclock_process_offset(pp, offset,
1361 up->tstamp[i], FUDGE +
1363 pp->lastref = up->timestamp;
1364 refclock_receive(peer);
1366 record_clock_stats(&peer->srcadr, pp->a_lastcode);
1367 } else if (pp->sloppyclockflag & CLK_FLAG4) {
1368 record_clock_stats(&peer->srcadr, pp->a_lastcode);
1372 printf("chu: timecode %d %s\n", pp->lencode,
1376 chu_newchan(peer, dtemp);
1380 refclock_report(peer, up->errflg);
1386 * chu_major - majority decoder
1390 struct peer *peer /* peer structure pointer */
1393 struct refclockproc *pp;
1396 u_char code[11]; /* decoded timecode */
1397 int mindist; /* minimum distance */
1398 int val1, val2; /* maximum distance */
1399 int synchar; /* stray cat */
1404 up = (struct chuunit *)pp->unitptr;
1407 * Majority decoder. Each burst encodes two replications at each
1408 * digit position in the timecode. Each row of the decoding
1409 * matrix encodes the number of occurrences of each digit found
1410 * at the corresponding position. The maximum over all
1411 * occurrences at each position is the distance for this
1412 * position and the corresponding digit is the maximum
1413 * likelihood candidate. If the distance is zero, assume a miss
1414 * '_'; if the distance is not more than half the total number
1415 * of occurrences, assume a soft error '*'; if two different
1416 * digits with the same distance are found, assume a hard error
1417 * '='. These will later cause a format error when the timecode
1418 * is interpreted. The decoding distance is defined as the
1419 * minimum distance over the first nine digits. The tenth digit
1420 * varies over the seconds, so we don't count it.
1423 for (i = 0; i < 9; i++) {
1426 for (j = 0; j < 16; j++) {
1427 temp = up->decode[i][j] + up->decode[i + 10][j];
1436 else if (val1 == val2)
1438 else if (val1 <= up->burstcnt)
1444 code[i] = hexchar[code[i]];
1449 * A valid timecode requires a minimum distance at least half
1450 * the total number of occurrences. A valid timecode also
1451 * requires at least 20 valid timestamps.
1453 if (up->burstcnt < MINBURST || mindist < up->burstcnt)
1454 up->status |= DECODE;
1455 if (up->ntstamp < MINSTAMP)
1456 up->status |= STAMP;
1459 * Compute the timecode timestamp from the days, hours and
1460 * minutes of the timecode. Use clocktime() for the aggregate
1461 * minutes and the minute offset computed from the burst
1462 * seconds. Note that this code relies on the filesystem time
1463 * for the years and does not use the years of the timecode.
1465 if (sscanf((char *)code, "%1x%3d%2d%2d", &synchar, &pp->day,
1466 &pp->hour, &pp->minute) != 4) {
1467 up->status |= AFORMAT;
1470 if (up->status & (DECODE | STAMP)) {
1471 up->errflg = CEVNT_BADREPLY;
1474 return (mindist * 100. / (2. * up->burstcnt));
1479 * chu_clear - clear decoding matrix
1483 struct peer *peer /* peer structure pointer */
1486 struct refclockproc *pp;
1491 up = (struct chuunit *)pp->unitptr;
1494 * Clear stuff for the minute.
1496 up->ndx = up->prevsec = 0;
1497 up->burstcnt = up->ntstamp = 0;
1498 up->status &= INSYNC;
1499 for (i = 0; i < 20; i++) {
1500 for (j = 0; j < 16; j++)
1501 up->decode[i][j] = 0;
1507 * chu_newchan - called once per minute to find the best channel;
1508 * returns zero on success, nonzero if ICOM error.
1517 struct refclockproc *pp;
1519 char tbuf[80]; /* trace buffer */
1525 up = (struct chuunit *)pp->unitptr;
1528 * The radio can be tuned to three channels: 0 (3330 kHz), 1
1529 * (7335 kHz) and 2 (14670 kHz). There are five one-minute
1530 * dwells in each cycle. During the first dwell the radio is
1531 * tuned to one of three probe channels; during the remaining
1532 * four dwells the radio is tuned to the data channel. The probe
1533 * channel is selects as the least recently used. At the end of
1534 * each dwell the channel metrics are measured and the highest
1535 * one is selected as the data channel.
1537 if (up->fd_icom <= 0)
1540 sp = &up->xmtr[up->achan];
1541 sp->metric -= sp->integ[sp->iptr];
1542 sp->integ[sp->iptr] = met;
1543 sp->metric += sp->integ[sp->iptr];
1544 sp->iptr = (sp->iptr + 1) % ISTAGE;
1547 for (i = 0; i < NCHAN; i++) {
1548 up->xmtr[i].probe++;
1550 up->xmtr[i].probe = 0;
1551 if (up->xmtr[i].metric < metric)
1553 metric = up->xmtr[i].metric;
1556 if (j != up->chan && metric > 0) {
1558 sprintf(tbuf, "chu: QSY to %.3f MHz metric %.0f",
1559 qsy[up->chan], metric);
1560 if (pp->sloppyclockflag & CLK_FLAG4)
1561 record_clock_stats(&peer->srcadr, tbuf);
1564 printf("%s\n", tbuf);
1569 * Start the next dwell. We speed up the initial sync a little.
1570 * If not in sync and no bursts were heard the previous dwell,
1571 * restart the probe.
1574 if (up->burstcnt == 0 && !(up->status & INSYNC))
1579 "chu: at %ld dwell %d achan %d metric %.0f chan %d\n",
1580 current_time, up->dwell, up->achan, sp->metric,
1583 if (up->dwell == 0) {
1585 for (i = 0; i < NCHAN; i++) {
1586 if (up->xmtr[i].probe < rval)
1588 rval = up->xmtr[i].probe;
1591 rval = icom_freq(up->fd_icom, peer->ttl & 0x7f,
1592 qsy[up->achan] + TUNE);
1595 printf("chu: at %ld probe channel %d\n",
1596 current_time, up->achan);
1599 if (up->achan != up->chan) {
1600 rval = icom_freq(up->fd_icom, peer->ttl & 0x7f,
1601 qsy[up->chan] + TUNE);
1602 up->achan = up->chan;
1605 sprintf(up->ident, "CHU%d", up->achan);
1606 memcpy(&peer->refid, up->ident, 4);
1607 up->dwell = (up->dwell + 1) % DWELL;
1613 * chu_dist - determine the distance of two octet arguments
1617 int x, /* an octet of bits */
1618 int y /* another octet of bits */
1621 int val; /* bit count */
1626 * The distance is determined as the weight of the exclusive OR
1627 * of the two arguments. The weight is determined by the number
1628 * of one bits in the result. Each one bit increases the weight,
1629 * while each zero bit decreases it.
1633 for (i = 0; i < 8; i++) {
1634 if ((temp & 0x1) == 0)
1646 * chu_gain - adjust codec gain
1648 * This routine is called once each second. If the signal envelope
1649 * amplitude is too low, the codec gain is bumped up by four units; if
1650 * too high, it is bumped down. The decoder is relatively insensitive to
1651 * amplitude, so this crudity works just fine. The input port is set and
1652 * the error flag is cleared, mostly to be ornery.
1656 struct peer *peer /* peer structure pointer */
1659 struct refclockproc *pp;
1663 up = (struct chuunit *)pp->unitptr;
1666 * Apparently, the codec uses only the high order bits of the
1667 * gain control field. Thus, it may take awhile for changes to
1668 * wiggle the hardware bits.
1670 if (up->clipcnt == 0) {
1672 if (up->gain > MAXGAIN)
1674 } else if (up->clipcnt > MAXCLP) {
1679 audio_gain(up->gain, up->mongain, up->port);
1682 #endif /* HAVE_AUDIO */
1686 int refclock_chu_bs;
1687 #endif /* REFCLOCK */