dist/ntp/ntpd/refclock_irig.c

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