/contrib/ntp/ntpd/refclock_wwv.c
C | 2709 lines | 1594 code | 204 blank | 911 comment | 264 complexity | d5531729a944f68437d2d9bd4047b8da MD5 | raw file
Possible License(s): MPL-2.0-no-copyleft-exception, BSD-3-Clause, LGPL-2.0, LGPL-2.1, BSD-2-Clause, 0BSD, JSON, AGPL-1.0, GPL-2.0
Large files files are truncated, but you can click here to view the full file
1/* 2 * refclock_wwv - clock driver for NIST WWV/H time/frequency station 3 */ 4#ifdef HAVE_CONFIG_H 5#include <config.h> 6#endif 7 8#if defined(REFCLOCK) && defined(CLOCK_WWV) 9 10#include "ntpd.h" 11#include "ntp_io.h" 12#include "ntp_refclock.h" 13#include "ntp_calendar.h" 14#include "ntp_stdlib.h" 15#include "audio.h" 16 17#include <stdio.h> 18#include <ctype.h> 19#include <math.h> 20#ifdef HAVE_SYS_IOCTL_H 21# include <sys/ioctl.h> 22#endif /* HAVE_SYS_IOCTL_H */ 23 24#define ICOM 1 25 26#ifdef ICOM 27#include "icom.h" 28#endif /* ICOM */ 29 30/* 31 * Audio WWV/H demodulator/decoder 32 * 33 * This driver synchronizes the computer time using data encoded in 34 * radio transmissions from NIST time/frequency stations WWV in Boulder, 35 * CO, and WWVH in Kauai, HI. Transmissions are made continuously on 36 * 2.5, 5, 10 and 15 MHz from WWV and WWVH, and 20 MHz from WWV. An 37 * ordinary AM 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 weasons, both day and night. 41 * 42 * The driver receives, demodulates and decodes the radio signals when 43 * connected to the audio codec of a workstation running Solaris, SunOS 44 * FreeBSD or Linux, and with a little help, other workstations with 45 * similar codecs or sound cards. In this implementation, only one audio 46 * driver and codec can be supported on a single machine. 47 * 48 * The demodulation and decoding algorithms used in this driver are 49 * based on those developed for the TAPR DSP93 development board and the 50 * TI 320C25 digital signal processor described in: Mills, D.L. A 51 * precision radio clock for WWV transmissions. Electrical Engineering 52 * Report 97-8-1, University of Delaware, August 1997, 25 pp., available 53 * from www.eecis.udel.edu/~mills/reports.html. The algorithms described 54 * in this report have been modified somewhat to improve performance 55 * under weak signal conditions and to provide an automatic station 56 * identification feature. 57 * 58 * The ICOM code is normally compiled in the driver. It isn't used, 59 * unless the mode keyword on the server configuration command specifies 60 * a nonzero ICOM ID select code. The C-IV trace is turned on if the 61 * debug level is greater than one. 62 * 63 * Fudge factors 64 * 65 * Fudge flag4 causes the dubugging output described above to be 66 * recorded in the clockstats file. Fudge flag2 selects the audio input 67 * port, where 0 is the mike port (default) and 1 is the line-in port. 68 * It does not seem useful to select the compact disc player port. Fudge 69 * flag3 enables audio monitoring of the input signal. For this purpose, 70 * the monitor gain is set to a default value. 71 */ 72/* 73 * General definitions. These ordinarily do not need to be changed. 74 */ 75#define DEVICE_AUDIO "/dev/audio" /* audio device name */ 76#define AUDIO_BUFSIZ 320 /* audio buffer size (50 ms) */ 77#define PRECISION (-10) /* precision assumed (about 1 ms) */ 78#define DESCRIPTION "WWV/H Audio Demodulator/Decoder" /* WRU */ 79#define SECOND 8000 /* second epoch (sample rate) (Hz) */ 80#define MINUTE (SECOND * 60) /* minute epoch */ 81#define OFFSET 128 /* companded sample offset */ 82#define SIZE 256 /* decompanding table size */ 83#define MAXAMP 6000. /* max signal level reference */ 84#define MAXCLP 100 /* max clips above reference per s */ 85#define MAXSNR 40. /* max SNR reference */ 86#define MAXFREQ 1.5 /* max frequency tolerance (187 PPM) */ 87#define DATCYC 170 /* data filter cycles */ 88#define DATSIZ (DATCYC * MS) /* data filter size */ 89#define SYNCYC 800 /* minute filter cycles */ 90#define SYNSIZ (SYNCYC * MS) /* minute filter size */ 91#define TCKCYC 5 /* tick filter cycles */ 92#define TCKSIZ (TCKCYC * MS) /* tick filter size */ 93#define NCHAN 5 /* number of radio channels */ 94#define AUDIO_PHI 5e-6 /* dispersion growth factor */ 95 96/* 97 * Tunable parameters. The DGAIN parameter can be changed to fit the 98 * audio response of the radio at 100 Hz. The WWV/WWVH data subcarrier 99 * is transmitted at about 20 percent percent modulation; the matched 100 * filter boosts it by a factor of 17 and the receiver response does 101 * what it does. The compromise value works for ICOM radios. If the 102 * radio is not tunable, the DCHAN parameter can be changed to fit the 103 * expected best propagation frequency: higher if further from the 104 * transmitter, lower if nearer. The compromise value works for the US 105 * right coast. The FREQ_OFFSET parameter can be used as a frequency 106 * vernier to correct codec requency if greater than MAXFREQ. 107 */ 108#define DCHAN 3 /* default radio channel (15 Mhz) */ 109#define DGAIN 5. /* subcarrier gain */ 110#define FREQ_OFFSET 0. /* codec frequency correction (PPM) */ 111 112/* 113 * General purpose status bits (status) 114 * 115 * SELV and/or SELH are set when WWV or WWVH have been heard and cleared 116 * on signal loss. SSYNC is set when the second sync pulse has been 117 * acquired and cleared by signal loss. MSYNC is set when the minute 118 * sync pulse has been acquired. DSYNC is set when the units digit has 119 * has reached the threshold and INSYNC is set when all nine digits have 120 * reached the threshold. The MSYNC, DSYNC and INSYNC bits are cleared 121 * only by timeout, upon which the driver starts over from scratch. 122 * 123 * DGATE is lit if the data bit amplitude or SNR is below thresholds and 124 * BGATE is lit if the pulse width amplitude or SNR is below thresolds. 125 * LEPSEC is set during the last minute of the leap day. At the end of 126 * this minute the driver inserts second 60 in the seconds state machine 127 * and the minute sync slips a second. 128 */ 129#define MSYNC 0x0001 /* minute epoch sync */ 130#define SSYNC 0x0002 /* second epoch sync */ 131#define DSYNC 0x0004 /* minute units sync */ 132#define INSYNC 0x0008 /* clock synchronized */ 133#define FGATE 0x0010 /* frequency gate */ 134#define DGATE 0x0020 /* data pulse amplitude error */ 135#define BGATE 0x0040 /* data pulse width error */ 136#define LEPSEC 0x1000 /* leap minute */ 137 138/* 139 * Station scoreboard bits 140 * 141 * These are used to establish the signal quality for each of the five 142 * frequencies and two stations. 143 */ 144#define SELV 0x0100 /* WWV station select */ 145#define SELH 0x0200 /* WWVH station select */ 146 147/* 148 * Alarm status bits (alarm) 149 * 150 * These bits indicate various alarm conditions, which are decoded to 151 * form the quality character included in the timecode. 152 */ 153#define CMPERR 1 /* digit or misc bit compare error */ 154#define LOWERR 2 /* low bit or digit amplitude or SNR */ 155#define NINERR 4 /* less than nine digits in minute */ 156#define SYNERR 8 /* not tracking second sync */ 157 158/* 159 * Watchcat timeouts (watch) 160 * 161 * If these timeouts expire, the status bits are mashed to zero and the 162 * driver starts from scratch. Suitably more refined procedures may be 163 * developed in future. All these are in minutes. 164 */ 165#define ACQSN 6 /* station acquisition timeout */ 166#define DATA 15 /* unit minutes timeout */ 167#define SYNCH 40 /* station sync timeout */ 168#define PANIC (2 * 1440) /* panic timeout */ 169 170/* 171 * Thresholds. These establish the minimum signal level, minimum SNR and 172 * maximum jitter thresholds which establish the error and false alarm 173 * rates of the driver. The values defined here may be on the 174 * adventurous side in the interest of the highest sensitivity. 175 */ 176#define MTHR 13. /* minute sync gate (percent) */ 177#define TTHR 50. /* minute sync threshold (percent) */ 178#define AWND 20 /* minute sync jitter threshold (ms) */ 179#define ATHR 2500. /* QRZ minute sync threshold */ 180#define ASNR 20. /* QRZ minute sync SNR threshold (dB) */ 181#define QTHR 2500. /* QSY minute sync threshold */ 182#define QSNR 20. /* QSY minute sync SNR threshold (dB) */ 183#define STHR 2500. /* second sync threshold */ 184#define SSNR 15. /* second sync SNR threshold (dB) */ 185#define SCMP 10 /* second sync compare threshold */ 186#define DTHR 1000. /* bit threshold */ 187#define DSNR 10. /* bit SNR threshold (dB) */ 188#define AMIN 3 /* min bit count */ 189#define AMAX 6 /* max bit count */ 190#define BTHR 1000. /* digit threshold */ 191#define BSNR 3. /* digit likelihood threshold (dB) */ 192#define BCMP 3 /* digit compare threshold */ 193#define MAXERR 40 /* maximum error alarm */ 194 195/* 196 * Tone frequency definitions. The increments are for 4.5-deg sine 197 * table. 198 */ 199#define MS (SECOND / 1000) /* samples per millisecond */ 200#define IN100 ((100 * 80) / SECOND) /* 100 Hz increment */ 201#define IN1000 ((1000 * 80) / SECOND) /* 1000 Hz increment */ 202#define IN1200 ((1200 * 80) / SECOND) /* 1200 Hz increment */ 203 204/* 205 * Acquisition and tracking time constants 206 */ 207#define MINAVG 8 /* min averaging time */ 208#define MAXAVG 1024 /* max averaging time */ 209#define FCONST 3 /* frequency time constant */ 210#define TCONST 16 /* data bit/digit time constant */ 211 212/* 213 * Miscellaneous status bits (misc) 214 * 215 * These bits correspond to designated bits in the WWV/H timecode. The 216 * bit probabilities are exponentially averaged over several minutes and 217 * processed by a integrator and threshold. 218 */ 219#define DUT1 0x01 /* 56 DUT .1 */ 220#define DUT2 0x02 /* 57 DUT .2 */ 221#define DUT4 0x04 /* 58 DUT .4 */ 222#define DUTS 0x08 /* 50 DUT sign */ 223#define DST1 0x10 /* 55 DST1 leap warning */ 224#define DST2 0x20 /* 2 DST2 DST1 delayed one day */ 225#define SECWAR 0x40 /* 3 leap second warning */ 226 227/* 228 * The on-time synchronization point for the driver is the second epoch 229 * sync pulse produced by the FIR matched filters. As the 5-ms delay of 230 * these filters is compensated, the program delay is 1.1 ms due to the 231 * 600-Hz IIR bandpass filter. The measured receiver delay is 4.7 ms and 232 * the codec delay less than 0.2 ms. The additional propagation delay 233 * specific to each receiver location can be programmed in the fudge 234 * time1 and time2 values for WWV and WWVH, respectively. 235 */ 236#define PDELAY (.0011 + .0047 + .0002) /* net system delay (s) */ 237 238/* 239 * Table of sine values at 4.5-degree increments. This is used by the 240 * synchronous matched filter demodulators. 241 */ 242double sintab[] = { 243 0.000000e+00, 7.845910e-02, 1.564345e-01, 2.334454e-01, /* 0-3 */ 244 3.090170e-01, 3.826834e-01, 4.539905e-01, 5.224986e-01, /* 4-7 */ 245 5.877853e-01, 6.494480e-01, 7.071068e-01, 7.604060e-01, /* 8-11 */ 246 8.090170e-01, 8.526402e-01, 8.910065e-01, 9.238795e-01, /* 12-15 */ 247 9.510565e-01, 9.723699e-01, 9.876883e-01, 9.969173e-01, /* 16-19 */ 248 1.000000e+00, 9.969173e-01, 9.876883e-01, 9.723699e-01, /* 20-23 */ 249 9.510565e-01, 9.238795e-01, 8.910065e-01, 8.526402e-01, /* 24-27 */ 250 8.090170e-01, 7.604060e-01, 7.071068e-01, 6.494480e-01, /* 28-31 */ 251 5.877853e-01, 5.224986e-01, 4.539905e-01, 3.826834e-01, /* 32-35 */ 252 3.090170e-01, 2.334454e-01, 1.564345e-01, 7.845910e-02, /* 36-39 */ 253-0.000000e+00, -7.845910e-02, -1.564345e-01, -2.334454e-01, /* 40-43 */ 254-3.090170e-01, -3.826834e-01, -4.539905e-01, -5.224986e-01, /* 44-47 */ 255-5.877853e-01, -6.494480e-01, -7.071068e-01, -7.604060e-01, /* 48-51 */ 256-8.090170e-01, -8.526402e-01, -8.910065e-01, -9.238795e-01, /* 52-55 */ 257-9.510565e-01, -9.723699e-01, -9.876883e-01, -9.969173e-01, /* 56-59 */ 258-1.000000e+00, -9.969173e-01, -9.876883e-01, -9.723699e-01, /* 60-63 */ 259-9.510565e-01, -9.238795e-01, -8.910065e-01, -8.526402e-01, /* 64-67 */ 260-8.090170e-01, -7.604060e-01, -7.071068e-01, -6.494480e-01, /* 68-71 */ 261-5.877853e-01, -5.224986e-01, -4.539905e-01, -3.826834e-01, /* 72-75 */ 262-3.090170e-01, -2.334454e-01, -1.564345e-01, -7.845910e-02, /* 76-79 */ 263 0.000000e+00}; /* 80 */ 264 265/* 266 * Decoder operations at the end of each second are driven by a state 267 * machine. The transition matrix consists of a dispatch table indexed 268 * by second number. Each entry in the table contains a case switch 269 * number and argument. 270 */ 271struct progx { 272 int sw; /* case switch number */ 273 int arg; /* argument */ 274}; 275 276/* 277 * Case switch numbers 278 */ 279#define IDLE 0 /* no operation */ 280#define COEF 1 /* BCD bit */ 281#define COEF1 2 /* BCD bit for minute unit */ 282#define COEF2 3 /* BCD bit not used */ 283#define DECIM9 4 /* BCD digit 0-9 */ 284#define DECIM6 5 /* BCD digit 0-6 */ 285#define DECIM3 6 /* BCD digit 0-3 */ 286#define DECIM2 7 /* BCD digit 0-2 */ 287#define MSCBIT 8 /* miscellaneous bit */ 288#define MSC20 9 /* miscellaneous bit */ 289#define MSC21 10 /* QSY probe channel */ 290#define MIN1 11 /* latch time */ 291#define MIN2 12 /* leap second */ 292#define SYNC2 13 /* latch minute sync pulse */ 293#define SYNC3 14 /* latch data pulse */ 294 295/* 296 * Offsets in decoding matrix 297 */ 298#define MN 0 /* minute digits (2) */ 299#define HR 2 /* hour digits (2) */ 300#define DA 4 /* day digits (3) */ 301#define YR 7 /* year digits (2) */ 302 303struct progx progx[] = { 304 {SYNC2, 0}, /* 0 latch minute sync pulse */ 305 {SYNC3, 0}, /* 1 latch data pulse */ 306 {MSCBIT, DST2}, /* 2 dst2 */ 307 {MSCBIT, SECWAR}, /* 3 lw */ 308 {COEF, 0}, /* 4 1 year units */ 309 {COEF, 1}, /* 5 2 */ 310 {COEF, 2}, /* 6 4 */ 311 {COEF, 3}, /* 7 8 */ 312 {DECIM9, YR}, /* 8 */ 313 {IDLE, 0}, /* 9 p1 */ 314 {COEF1, 0}, /* 10 1 minute units */ 315 {COEF1, 1}, /* 11 2 */ 316 {COEF1, 2}, /* 12 4 */ 317 {COEF1, 3}, /* 13 8 */ 318 {DECIM9, MN}, /* 14 */ 319 {COEF, 0}, /* 15 10 minute tens */ 320 {COEF, 1}, /* 16 20 */ 321 {COEF, 2}, /* 17 40 */ 322 {COEF2, 3}, /* 18 80 (not used) */ 323 {DECIM6, MN + 1}, /* 19 p2 */ 324 {COEF, 0}, /* 20 1 hour units */ 325 {COEF, 1}, /* 21 2 */ 326 {COEF, 2}, /* 22 4 */ 327 {COEF, 3}, /* 23 8 */ 328 {DECIM9, HR}, /* 24 */ 329 {COEF, 0}, /* 25 10 hour tens */ 330 {COEF, 1}, /* 26 20 */ 331 {COEF2, 2}, /* 27 40 (not used) */ 332 {COEF2, 3}, /* 28 80 (not used) */ 333 {DECIM2, HR + 1}, /* 29 p3 */ 334 {COEF, 0}, /* 30 1 day units */ 335 {COEF, 1}, /* 31 2 */ 336 {COEF, 2}, /* 32 4 */ 337 {COEF, 3}, /* 33 8 */ 338 {DECIM9, DA}, /* 34 */ 339 {COEF, 0}, /* 35 10 day tens */ 340 {COEF, 1}, /* 36 20 */ 341 {COEF, 2}, /* 37 40 */ 342 {COEF, 3}, /* 38 80 */ 343 {DECIM9, DA + 1}, /* 39 p4 */ 344 {COEF, 0}, /* 40 100 day hundreds */ 345 {COEF, 1}, /* 41 200 */ 346 {COEF2, 2}, /* 42 400 (not used) */ 347 {COEF2, 3}, /* 43 800 (not used) */ 348 {DECIM3, DA + 2}, /* 44 */ 349 {IDLE, 0}, /* 45 */ 350 {IDLE, 0}, /* 46 */ 351 {IDLE, 0}, /* 47 */ 352 {IDLE, 0}, /* 48 */ 353 {IDLE, 0}, /* 49 p5 */ 354 {MSCBIT, DUTS}, /* 50 dut+- */ 355 {COEF, 0}, /* 51 10 year tens */ 356 {COEF, 1}, /* 52 20 */ 357 {COEF, 2}, /* 53 40 */ 358 {COEF, 3}, /* 54 80 */ 359 {MSC20, DST1}, /* 55 dst1 */ 360 {MSCBIT, DUT1}, /* 56 0.1 dut */ 361 {MSCBIT, DUT2}, /* 57 0.2 */ 362 {MSC21, DUT4}, /* 58 0.4 QSY probe channel */ 363 {MIN1, 0}, /* 59 p6 latch time */ 364 {MIN2, 0} /* 60 leap second */ 365}; 366 367/* 368 * BCD coefficients for maximum likelihood digit decode 369 */ 370#define P15 1. /* max positive number */ 371#define N15 -1. /* max negative number */ 372 373/* 374 * Digits 0-9 375 */ 376#define P9 (P15 / 4) /* mark (+1) */ 377#define N9 (N15 / 4) /* space (-1) */ 378 379double bcd9[][4] = { 380 {N9, N9, N9, N9}, /* 0 */ 381 {P9, N9, N9, N9}, /* 1 */ 382 {N9, P9, N9, N9}, /* 2 */ 383 {P9, P9, N9, N9}, /* 3 */ 384 {N9, N9, P9, N9}, /* 4 */ 385 {P9, N9, P9, N9}, /* 5 */ 386 {N9, P9, P9, N9}, /* 6 */ 387 {P9, P9, P9, N9}, /* 7 */ 388 {N9, N9, N9, P9}, /* 8 */ 389 {P9, N9, N9, P9}, /* 9 */ 390 {0, 0, 0, 0} /* backstop */ 391}; 392 393/* 394 * Digits 0-6 (minute tens) 395 */ 396#define P6 (P15 / 3) /* mark (+1) */ 397#define N6 (N15 / 3) /* space (-1) */ 398 399double bcd6[][4] = { 400 {N6, N6, N6, 0}, /* 0 */ 401 {P6, N6, N6, 0}, /* 1 */ 402 {N6, P6, N6, 0}, /* 2 */ 403 {P6, P6, N6, 0}, /* 3 */ 404 {N6, N6, P6, 0}, /* 4 */ 405 {P6, N6, P6, 0}, /* 5 */ 406 {N6, P6, P6, 0}, /* 6 */ 407 {0, 0, 0, 0} /* backstop */ 408}; 409 410/* 411 * Digits 0-3 (day hundreds) 412 */ 413#define P3 (P15 / 2) /* mark (+1) */ 414#define N3 (N15 / 2) /* space (-1) */ 415 416double bcd3[][4] = { 417 {N3, N3, 0, 0}, /* 0 */ 418 {P3, N3, 0, 0}, /* 1 */ 419 {N3, P3, 0, 0}, /* 2 */ 420 {P3, P3, 0, 0}, /* 3 */ 421 {0, 0, 0, 0} /* backstop */ 422}; 423 424/* 425 * Digits 0-2 (hour tens) 426 */ 427#define P2 (P15 / 2) /* mark (+1) */ 428#define N2 (N15 / 2) /* space (-1) */ 429 430double bcd2[][4] = { 431 {N2, N2, 0, 0}, /* 0 */ 432 {P2, N2, 0, 0}, /* 1 */ 433 {N2, P2, 0, 0}, /* 2 */ 434 {0, 0, 0, 0} /* backstop */ 435}; 436 437/* 438 * DST decode (DST2 DST1) for prettyprint 439 */ 440char dstcod[] = { 441 'S', /* 00 standard time */ 442 'I', /* 01 set clock ahead at 0200 local */ 443 'O', /* 10 set clock back at 0200 local */ 444 'D' /* 11 daylight time */ 445}; 446 447/* 448 * The decoding matrix consists of nine row vectors, one for each digit 449 * of the timecode. The digits are stored from least to most significant 450 * order. The maximum likelihood timecode is formed from the digits 451 * corresponding to the maximum likelihood values reading in the 452 * opposite order: yy ddd hh:mm. 453 */ 454struct decvec { 455 int radix; /* radix (3, 4, 6, 10) */ 456 int digit; /* current clock digit */ 457 int mldigit; /* maximum likelihood digit */ 458 int count; /* match count */ 459 double digprb; /* max digit probability */ 460 double digsnr; /* likelihood function (dB) */ 461 double like[10]; /* likelihood integrator 0-9 */ 462}; 463 464/* 465 * The station structure (sp) is used to acquire the minute pulse from 466 * WWV and/or WWVH. These stations are distinguished by the frequency 467 * used for the second and minute sync pulses, 1000 Hz for WWV and 1200 468 * Hz for WWVH. Other than frequency, the format is the same. 469 */ 470struct sync { 471 double epoch; /* accumulated epoch differences */ 472 double maxeng; /* sync max energy */ 473 double noieng; /* sync noise energy */ 474 long pos; /* max amplitude position */ 475 long lastpos; /* last max position */ 476 long mepoch; /* minute synch epoch */ 477 478 double amp; /* sync signal */ 479 double syneng; /* sync signal max */ 480 double synmax; /* sync signal max latched at 0 s */ 481 double synsnr; /* sync signal SNR */ 482 double metric; /* signal quality metric */ 483 int reach; /* reachability register */ 484 int count; /* bit counter */ 485 int select; /* select bits */ 486 char refid[5]; /* reference identifier */ 487}; 488 489/* 490 * The channel structure (cp) is used to mitigate between channels. 491 */ 492struct chan { 493 int gain; /* audio gain */ 494 struct sync wwv; /* wwv station */ 495 struct sync wwvh; /* wwvh station */ 496}; 497 498/* 499 * WWV unit control structure (up) 500 */ 501struct wwvunit { 502 l_fp timestamp; /* audio sample timestamp */ 503 l_fp tick; /* audio sample increment */ 504 double phase, freq; /* logical clock phase and frequency */ 505 double monitor; /* audio monitor point */ 506#ifdef ICOM 507 int fd_icom; /* ICOM file descriptor */ 508#endif /* ICOM */ 509 int errflg; /* error flags */ 510 int watch; /* watchcat */ 511 512 /* 513 * Audio codec variables 514 */ 515 double comp[SIZE]; /* decompanding table */ 516 int port; /* codec port */ 517 int gain; /* codec gain */ 518 int mongain; /* codec monitor gain */ 519 int clipcnt; /* sample clipped count */ 520 521 /* 522 * Variables used to establish basic system timing 523 */ 524 int avgint; /* master time constant */ 525 int yepoch; /* sync epoch */ 526 int repoch; /* buffered sync epoch */ 527 double epomax; /* second sync amplitude */ 528 double eposnr; /* second sync SNR */ 529 double irig; /* data I channel amplitude */ 530 double qrig; /* data Q channel amplitude */ 531 int datapt; /* 100 Hz ramp */ 532 double datpha; /* 100 Hz VFO control */ 533 int rphase; /* second sample counter */ 534 long mphase; /* minute sample counter */ 535 536 /* 537 * Variables used to mitigate which channel to use 538 */ 539 struct chan mitig[NCHAN]; /* channel data */ 540 struct sync *sptr; /* station pointer */ 541 int dchan; /* data channel */ 542 int schan; /* probe channel */ 543 int achan; /* active channel */ 544 545 /* 546 * Variables used by the clock state machine 547 */ 548 struct decvec decvec[9]; /* decoding matrix */ 549 int rsec; /* seconds counter */ 550 int digcnt; /* count of digits synchronized */ 551 552 /* 553 * Variables used to estimate signal levels and bit/digit 554 * probabilities 555 */ 556 double datsig; /* data signal max */ 557 double datsnr; /* data signal SNR (dB) */ 558 559 /* 560 * Variables used to establish status and alarm conditions 561 */ 562 int status; /* status bits */ 563 int alarm; /* alarm flashers */ 564 int misc; /* miscellaneous timecode bits */ 565 int errcnt; /* data bit error counter */ 566}; 567 568/* 569 * Function prototypes 570 */ 571static int wwv_start P((int, struct peer *)); 572static void wwv_shutdown P((int, struct peer *)); 573static void wwv_receive P((struct recvbuf *)); 574static void wwv_poll P((int, struct peer *)); 575 576/* 577 * More function prototypes 578 */ 579static void wwv_epoch P((struct peer *)); 580static void wwv_rf P((struct peer *, double)); 581static void wwv_endpoc P((struct peer *, int)); 582static void wwv_rsec P((struct peer *, double)); 583static void wwv_qrz P((struct peer *, struct sync *, int)); 584static void wwv_corr4 P((struct peer *, struct decvec *, 585 double [], double [][4])); 586static void wwv_gain P((struct peer *)); 587static void wwv_tsec P((struct peer *)); 588static int timecode P((struct wwvunit *, char *)); 589static double wwv_snr P((double, double)); 590static int carry P((struct decvec *)); 591static int wwv_newchan P((struct peer *)); 592static void wwv_newgame P((struct peer *)); 593static double wwv_metric P((struct sync *)); 594static void wwv_clock P((struct peer *)); 595#ifdef ICOM 596static int wwv_qsy P((struct peer *, int)); 597#endif /* ICOM */ 598 599static double qsy[NCHAN] = {2.5, 5, 10, 15, 20}; /* frequencies (MHz) */ 600 601/* 602 * Transfer vector 603 */ 604struct refclock refclock_wwv = { 605 wwv_start, /* start up driver */ 606 wwv_shutdown, /* shut down driver */ 607 wwv_poll, /* transmit poll message */ 608 noentry, /* not used (old wwv_control) */ 609 noentry, /* initialize driver (not used) */ 610 noentry, /* not used (old wwv_buginfo) */ 611 NOFLAGS /* not used */ 612}; 613 614 615/* 616 * wwv_start - open the devices and initialize data for processing 617 */ 618static int 619wwv_start( 620 int unit, /* instance number (used by PCM) */ 621 struct peer *peer /* peer structure pointer */ 622 ) 623{ 624 struct refclockproc *pp; 625 struct wwvunit *up; 626#ifdef ICOM 627 int temp; 628#endif /* ICOM */ 629 630 /* 631 * Local variables 632 */ 633 int fd; /* file descriptor */ 634 int i; /* index */ 635 double step; /* codec adjustment */ 636 637 /* 638 * Open audio device 639 */ 640 fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit); 641 if (fd < 0) 642 return (0); 643#ifdef DEBUG 644 if (debug) 645 audio_show(); 646#endif /* DEBUG */ 647 648 /* 649 * Allocate and initialize unit structure 650 */ 651 if (!(up = (struct wwvunit *)emalloc(sizeof(struct wwvunit)))) { 652 close(fd); 653 return (0); 654 } 655 memset(up, 0, sizeof(struct wwvunit)); 656 pp = peer->procptr; 657 pp->unitptr = (caddr_t)up; 658 pp->io.clock_recv = wwv_receive; 659 pp->io.srcclock = (caddr_t)peer; 660 pp->io.datalen = 0; 661 pp->io.fd = fd; 662 if (!io_addclock(&pp->io)) { 663 close(fd); 664 free(up); 665 return (0); 666 } 667 668 /* 669 * Initialize miscellaneous variables 670 */ 671 peer->precision = PRECISION; 672 pp->clockdesc = DESCRIPTION; 673 674 /* 675 * The companded samples are encoded sign-magnitude. The table 676 * contains all the 256 values in the interest of speed. 677 */ 678 up->comp[0] = up->comp[OFFSET] = 0.; 679 up->comp[1] = 1.; up->comp[OFFSET + 1] = -1.; 680 up->comp[2] = 3.; up->comp[OFFSET + 2] = -3.; 681 step = 2.; 682 for (i = 3; i < OFFSET; i++) { 683 up->comp[i] = up->comp[i - 1] + step; 684 up->comp[OFFSET + i] = -up->comp[i]; 685 if (i % 16 == 0) 686 step *= 2.; 687 } 688 DTOLFP(1. / SECOND, &up->tick); 689 690 /* 691 * Initialize the decoding matrix with the radix for each digit 692 * position. 693 */ 694 up->decvec[MN].radix = 10; /* minutes */ 695 up->decvec[MN + 1].radix = 6; 696 up->decvec[HR].radix = 10; /* hours */ 697 up->decvec[HR + 1].radix = 3; 698 up->decvec[DA].radix = 10; /* days */ 699 up->decvec[DA + 1].radix = 10; 700 up->decvec[DA + 2].radix = 4; 701 up->decvec[YR].radix = 10; /* years */ 702 up->decvec[YR + 1].radix = 10; 703 704#ifdef ICOM 705 /* 706 * Initialize autotune if available. Note that the ICOM select 707 * code must be less than 128, so the high order bit can be used 708 * to select the line speed 0 (9600 bps) or 1 (1200 bps). 709 */ 710 temp = 0; 711#ifdef DEBUG 712 if (debug > 1) 713 temp = P_TRACE; 714#endif /* DEBUG */ 715 if (peer->ttl != 0) { 716 if (peer->ttl & 0x80) 717 up->fd_icom = icom_init("/dev/icom", B1200, 718 temp); 719 else 720 up->fd_icom = icom_init("/dev/icom", B9600, 721 temp); 722 if (up->fd_icom < 0) { 723 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT) 724 msyslog(LOG_NOTICE, 725 "icom: %m"); 726 up->errflg = CEVNT_FAULT; 727 } 728 } 729 if (up->fd_icom > 0) { 730 if (wwv_qsy(peer, DCHAN) != 0) { 731 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT) 732 msyslog(LOG_NOTICE, 733 "icom: radio not found"); 734 up->errflg = CEVNT_FAULT; 735 close(up->fd_icom); 736 up->fd_icom = 0; 737 } else { 738 NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT) 739 msyslog(LOG_NOTICE, 740 "icom: autotune enabled"); 741 } 742 } 743#endif /* ICOM */ 744 745 /* 746 * Let the games begin. 747 */ 748 wwv_newgame(peer); 749 return (1); 750} 751 752 753/* 754 * wwv_shutdown - shut down the clock 755 */ 756static void 757wwv_shutdown( 758 int unit, /* instance number (not used) */ 759 struct peer *peer /* peer structure pointer */ 760 ) 761{ 762 struct refclockproc *pp; 763 struct wwvunit *up; 764 765 pp = peer->procptr; 766 up = (struct wwvunit *)pp->unitptr; 767 if (up == NULL) 768 return; 769 770 io_closeclock(&pp->io); 771#ifdef ICOM 772 if (up->fd_icom > 0) 773 close(up->fd_icom); 774#endif /* ICOM */ 775 free(up); 776} 777 778 779/* 780 * wwv_receive - receive data from the audio device 781 * 782 * This routine reads input samples and adjusts the logical clock to 783 * track the A/D sample clock by dropping or duplicating codec samples. 784 * It also controls the A/D signal level with an AGC loop to mimimize 785 * quantization noise and avoid overload. 786 */ 787static void 788wwv_receive( 789 struct recvbuf *rbufp /* receive buffer structure pointer */ 790 ) 791{ 792 struct peer *peer; 793 struct refclockproc *pp; 794 struct wwvunit *up; 795 796 /* 797 * Local variables 798 */ 799 double sample; /* codec sample */ 800 u_char *dpt; /* buffer pointer */ 801 int bufcnt; /* buffer counter */ 802 l_fp ltemp; 803 804 peer = (struct peer *)rbufp->recv_srcclock; 805 pp = peer->procptr; 806 up = (struct wwvunit *)pp->unitptr; 807 808 /* 809 * Main loop - read until there ain't no more. Note codec 810 * samples are bit-inverted. 811 */ 812 DTOLFP((double)rbufp->recv_length / SECOND, <emp); 813 L_SUB(&rbufp->recv_time, <emp); 814 up->timestamp = rbufp->recv_time; 815 dpt = rbufp->recv_buffer; 816 for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) { 817 sample = up->comp[~*dpt++ & 0xff]; 818 819 /* 820 * Clip noise spikes greater than MAXAMP (6000) and 821 * record the number of clips to be used later by the 822 * AGC. 823 */ 824 if (sample > MAXAMP) { 825 sample = MAXAMP; 826 up->clipcnt++; 827 } else if (sample < -MAXAMP) { 828 sample = -MAXAMP; 829 up->clipcnt++; 830 } 831 832 /* 833 * Variable frequency oscillator. The codec oscillator 834 * runs at the nominal rate of 8000 samples per second, 835 * or 125 us per sample. A frequency change of one unit 836 * results in either duplicating or deleting one sample 837 * per second, which results in a frequency change of 838 * 125 PPM. 839 */ 840 up->phase += up->freq / SECOND; 841 up->phase += FREQ_OFFSET / 1e6; 842 if (up->phase >= .5) { 843 up->phase -= 1.; 844 } else if (up->phase < -.5) { 845 up->phase += 1.; 846 wwv_rf(peer, sample); 847 wwv_rf(peer, sample); 848 } else { 849 wwv_rf(peer, sample); 850 } 851 L_ADD(&up->timestamp, &up->tick); 852 } 853 854 /* 855 * Set the input port and monitor gain for the next buffer. 856 */ 857 if (pp->sloppyclockflag & CLK_FLAG2) 858 up->port = 2; 859 else 860 up->port = 1; 861 if (pp->sloppyclockflag & CLK_FLAG3) 862 up->mongain = MONGAIN; 863 else 864 up->mongain = 0; 865} 866 867 868/* 869 * wwv_poll - called by the transmit procedure 870 * 871 * This routine keeps track of status. If no offset samples have been 872 * processed during a poll interval, a timeout event is declared. If 873 * errors have have occurred during the interval, they are reported as 874 * well. 875 */ 876static void 877wwv_poll( 878 int unit, /* instance number (not used) */ 879 struct peer *peer /* peer structure pointer */ 880 ) 881{ 882 struct refclockproc *pp; 883 struct wwvunit *up; 884 885 pp = peer->procptr; 886 up = (struct wwvunit *)pp->unitptr; 887 if (pp->coderecv == pp->codeproc) 888 up->errflg = CEVNT_TIMEOUT; 889 if (up->errflg) 890 refclock_report(peer, up->errflg); 891 up->errflg = 0; 892 pp->polls++; 893} 894 895 896/* 897 * wwv_rf - process signals and demodulate to baseband 898 * 899 * This routine grooms and filters decompanded raw audio samples. The 900 * output signal is the 100-Hz filtered baseband data signal in 901 * quadrature phase. The routine also determines the minute synch epoch, 902 * as well as certain signal maxima, minima and related values. 903 * 904 * There are two 1-s ramps used by this program. Both count the 8000 905 * logical clock samples spanning exactly one second. The epoch ramp 906 * counts the samples starting at an arbitrary time. The rphase ramp 907 * counts the samples starting at the 5-ms second sync pulse found 908 * during the epoch ramp. 909 * 910 * There are two 1-m ramps used by this program. The mphase ramp counts 911 * the 480,000 logical clock samples spanning exactly one minute and 912 * starting at an arbitrary time. The rsec ramp counts the 60 seconds of 913 * the minute starting at the 800-ms minute sync pulse found during the 914 * mphase ramp. The rsec ramp drives the seconds state machine to 915 * determine the bits and digits of the timecode. 916 * 917 * Demodulation operations are based on three synthesized quadrature 918 * sinusoids: 100 Hz for the data signal, 1000 Hz for the WWV sync 919 * signal and 1200 Hz for the WWVH sync signal. These drive synchronous 920 * matched filters for the data signal (170 ms at 100 Hz), WWV minute 921 * sync signal (800 ms at 1000 Hz) and WWVH minute sync signal (800 ms 922 * at 1200 Hz). Two additional matched filters are switched in 923 * as required for the WWV second sync signal (5 cycles at 1000 Hz) and 924 * WWVH second sync signal (6 cycles at 1200 Hz). 925 */ 926static void 927wwv_rf( 928 struct peer *peer, /* peerstructure pointer */ 929 double isig /* input signal */ 930 ) 931{ 932 struct refclockproc *pp; 933 struct wwvunit *up; 934 struct sync *sp, *rp; 935 936 static double lpf[5]; /* 150-Hz lpf delay line */ 937 double data; /* lpf output */ 938 static double bpf[9]; /* 1000/1200-Hz bpf delay line */ 939 double syncx; /* bpf output */ 940 static double mf[41]; /* 1000/1200-Hz mf delay line */ 941 double mfsync; /* mf output */ 942 943 static int iptr; /* data channel pointer */ 944 static double ibuf[DATSIZ]; /* data I channel delay line */ 945 static double qbuf[DATSIZ]; /* data Q channel delay line */ 946 947 static int jptr; /* sync channel pointer */ 948 static int kptr; /* tick channel pointer */ 949 950 static int csinptr; /* wwv channel phase */ 951 static double cibuf[SYNSIZ]; /* wwv I channel delay line */ 952 static double cqbuf[SYNSIZ]; /* wwv Q channel delay line */ 953 static double ciamp; /* wwv I channel amplitude */ 954 static double cqamp; /* wwv Q channel amplitude */ 955 956 static double csibuf[TCKSIZ]; /* wwv I tick delay line */ 957 static double csqbuf[TCKSIZ]; /* wwv Q tick delay line */ 958 static double csiamp; /* wwv I tick amplitude */ 959 static double csqamp; /* wwv Q tick amplitude */ 960 961 static int hsinptr; /* wwvh channel phase */ 962 static double hibuf[SYNSIZ]; /* wwvh I channel delay line */ 963 static double hqbuf[SYNSIZ]; /* wwvh Q channel delay line */ 964 static double hiamp; /* wwvh I channel amplitude */ 965 static double hqamp; /* wwvh Q channel amplitude */ 966 967 static double hsibuf[TCKSIZ]; /* wwvh I tick delay line */ 968 static double hsqbuf[TCKSIZ]; /* wwvh Q tick delay line */ 969 static double hsiamp; /* wwvh I tick amplitude */ 970 static double hsqamp; /* wwvh Q tick amplitude */ 971 972 static double epobuf[SECOND]; /* second sync comb filter */ 973 static double epomax, nxtmax; /* second sync amplitude buffer */ 974 static int epopos; /* epoch second sync position buffer */ 975 976 static int iniflg; /* initialization flag */ 977 int pdelay; /* propagation delay (samples) */ 978 int epoch; /* comb filter index */ 979 double dtemp; 980 int i; 981 982 pp = peer->procptr; 983 up = (struct wwvunit *)pp->unitptr; 984 985 if (!iniflg) { 986 iniflg = 1; 987 memset((char *)lpf, 0, sizeof(lpf)); 988 memset((char *)bpf, 0, sizeof(bpf)); 989 memset((char *)mf, 0, sizeof(mf)); 990 memset((char *)ibuf, 0, sizeof(ibuf)); 991 memset((char *)qbuf, 0, sizeof(qbuf)); 992 memset((char *)cibuf, 0, sizeof(cibuf)); 993 memset((char *)cqbuf, 0, sizeof(cqbuf)); 994 memset((char *)csibuf, 0, sizeof(csibuf)); 995 memset((char *)csqbuf, 0, sizeof(csqbuf)); 996 memset((char *)hibuf, 0, sizeof(hibuf)); 997 memset((char *)hqbuf, 0, sizeof(hqbuf)); 998 memset((char *)hsibuf, 0, sizeof(hsibuf)); 999 memset((char *)hsqbuf, 0, sizeof(hsqbuf)); 1000 memset((char *)epobuf, 0, sizeof(epobuf)); 1001 } 1002 1003 /* 1004 * Baseband data demodulation. The 100-Hz subcarrier is 1005 * extracted using a 150-Hz IIR lowpass filter. This attenuates 1006 * the 1000/1200-Hz sync signals, as well as the 440-Hz and 1007 * 600-Hz tones and most of the noise and voice modulation 1008 * components. 1009 * 1010 * The subcarrier is transmitted 10 dB down from the carrier. 1011 * The DGAIN parameter can be adjusted for this and to 1012 * compensate for the radio audio response at 100 Hz. 1013 * 1014 * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB 1015 * passband ripple, -50 dB stopband ripple. 1016 */ 1017 data = (lpf[4] = lpf[3]) * 8.360961e-01; 1018 data += (lpf[3] = lpf[2]) * -3.481740e+00; 1019 data += (lpf[2] = lpf[1]) * 5.452988e+00; 1020 data += (lpf[1] = lpf[0]) * -3.807229e+00; 1021 lpf[0] = isig * DGAIN - data; 1022 data = lpf[0] * 3.281435e-03 1023 + lpf[1] * -1.149947e-02 1024 + lpf[2] * 1.654858e-02 1025 + lpf[3] * -1.149947e-02 1026 + lpf[4] * 3.281435e-03; 1027 1028 /* 1029 * The 100-Hz data signal is demodulated using a pair of 1030 * quadrature multipliers, matched filters and a phase lock 1031 * loop. The I and Q quadrature data signals are produced by 1032 * multiplying the filtered signal by 100-Hz sine and cosine 1033 * signals, respectively. The signals are processed by 170-ms 1034 * synchronous matched filters to produce the amplitude and 1035 * phase signals used by the demodulator. The signals are scaled 1036 * to produce unit energy at the maximum value. 1037 */ 1038 i = up->datapt; 1039 up->datapt = (up->datapt + IN100) % 80; 1040 dtemp = sintab[i] * data / (MS / 2. * DATCYC); 1041 up->irig -= ibuf[iptr]; 1042 ibuf[iptr] = dtemp; 1043 up->irig += dtemp; 1044 1045 i = (i + 20) % 80; 1046 dtemp = sintab[i] * data / (MS / 2. * DATCYC); 1047 up->qrig -= qbuf[iptr]; 1048 qbuf[iptr] = dtemp; 1049 up->qrig += dtemp; 1050 iptr = (iptr + 1) % DATSIZ; 1051 1052 /* 1053 * Baseband sync demodulation. The 1000/1200 sync signals are 1054 * extracted using a 600-Hz IIR bandpass filter. This removes 1055 * the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz 1056 * tones and most of the noise and voice modulation components. 1057 * 1058 * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB 1059 * passband ripple, -50 dB stopband ripple. 1060 */ 1061 syncx = (bpf[8] = bpf[7]) * 4.897278e-01; 1062 syncx += (bpf[7] = bpf[6]) * -2.765914e+00; 1063 syncx += (bpf[6] = bpf[5]) * 8.110921e+00; 1064 syncx += (bpf[5] = bpf[4]) * -1.517732e+01; 1065 syncx += (bpf[4] = bpf[3]) * 1.975197e+01; 1066 syncx += (bpf[3] = bpf[2]) * -1.814365e+01; 1067 syncx += (bpf[2] = bpf[1]) * 1.159783e+01; 1068 syncx += (bpf[1] = bpf[0]) * -4.735040e+00; 1069 bpf[0] = isig - syncx; 1070 syncx = bpf[0] * 8.203628e-03 1071 + bpf[1] * -2.375732e-02 1072 + bpf[2] * 3.353214e-02 1073 + bpf[3] * -4.080258e-02 1074 + bpf[4] * 4.605479e-02 1075 + bpf[5] * -4.080258e-02 1076 + bpf[6] * 3.353214e-02 1077 + bpf[7] * -2.375732e-02 1078 + bpf[8] * 8.203628e-03; 1079 1080 /* 1081 * The 1000/1200 sync signals are demodulated using a pair of 1082 * quadrature multipliers and matched filters. However, 1083 * synchronous demodulation at these frequencies is impractical, 1084 * so only the signal amplitude is used. The I and Q quadrature 1085 * sync signals are produced by multiplying the filtered signal 1086 * by 1000-Hz (WWV) and 1200-Hz (WWVH) sine and cosine signals, 1087 * respectively. The WWV and WWVH signals are processed by 800- 1088 * ms synchronous matched filters and combined to produce the 1089 * minute sync signal and detect which one (or both) the WWV or 1090 * WWVH signal is present. The WWV and WWVH signals are also 1091 * processed by 5-ms synchronous matched filters and combined to 1092 * produce the second sync signal. The signals are scaled to 1093 * produce unit energy at the maximum value. 1094 * 1095 * Note the master timing ramps, which run continuously. The 1096 * minute counter (mphase) counts the samples in the minute, 1097 * while the second counter (epoch) counts the samples in the 1098 * second. 1099 */ 1100 up->mphase = (up->mphase + 1) % MINUTE; 1101 epoch = up->mphase % SECOND; 1102 1103 /* 1104 * WWV 1105 */ 1106 i = csinptr; 1107 csinptr = (csinptr + IN1000) % 80; 1108 1109 dtemp = sintab[i] * syncx / (MS / 2.); 1110 ciamp -= cibuf[jptr]; 1111 cibuf[jptr] = dtemp; 1112 ciamp += dtemp; 1113 csiamp -= csibuf[kptr]; 1114 csibuf[kptr] = dtemp; 1115 csiamp += dtemp; 1116 1117 i = (i + 20) % 80; 1118 dtemp = sintab[i] * syncx / (MS / 2.); 1119 cqamp -= cqbuf[jptr]; 1120 cqbuf[jptr] = dtemp; 1121 cqamp += dtemp; 1122 csqamp -= csqbuf[kptr]; 1123 csqbuf[kptr] = dtemp; 1124 csqamp += dtemp; 1125 1126 sp = &up->mitig[up->achan].wwv; 1127 sp->amp = sqrt(ciamp * ciamp + cqamp * cqamp) / SYNCYC; 1128 if (!(up->status & MSYNC)) 1129 wwv_qrz(peer, sp, (int)(pp->fudgetime1 * SECOND)); 1130 1131 /* 1132 * WWVH 1133 */ 1134 i = hsinptr; 1135 hsinptr = (hsinptr + IN1200) % 80; 1136 1137 dtemp = sintab[i] * syncx / (MS / 2.); 1138 hiamp -= hibuf[jptr]; 1139 hibuf[jptr] = dtemp; 1140 hiamp += dtemp; 1141 hsiamp -= hsibuf[kptr]; 1142 hsibuf[kptr] = dtemp; 1143 hsiamp += dtemp; 1144 1145 i = (i + 20) % 80; 1146 dtemp = sintab[i] * syncx / (MS / 2.); 1147 hqamp -= hqbuf[jptr]; 1148 hqbuf[jptr] = dtemp; 1149 hqamp += dtemp; 1150 hsqamp -= hsqbuf[kptr]; 1151 hsqbuf[kptr] = dtemp; 1152 hsqamp += dtemp; 1153 1154 rp = &up->mitig[up->achan].wwvh; 1155 rp->amp = sqrt(hiamp * hiamp + hqamp * hqamp) / SYNCYC; 1156 if (!(up->status & MSYNC)) 1157 wwv_qrz(peer, rp, (int)(pp->fudgetime2 * SECOND)); 1158 jptr = (jptr + 1) % SYNSIZ; 1159 kptr = (kptr + 1) % TCKSIZ; 1160 1161 /* 1162 * The following section is called once per minute. It does 1163 * housekeeping and timeout functions and empties the dustbins. 1164 */ 1165 if (up->mphase == 0) { 1166 up->watch++; 1167 if (!(up->status & MSYNC)) { 1168 1169 /* 1170 * If minute sync has not been acquired before 1171 * ACQSN timeout (6 min), or if no signal is 1172 * heard, the program cycles to the next 1173 * frequency and tries again. 1174 */ 1175 if (!wwv_newchan(peer)) 1176 up->watch = 0; 1177#ifdef ICOM 1178 if (up->fd_icom > 0) 1179 wwv_qsy(peer, up->dchan); 1180#endif /* ICOM */ 1181 } else { 1182 1183 /* 1184 * If the leap bit is set, set the minute epoch 1185 * back one second so the station processes 1186 * don't miss a beat. 1187 */ 1188 if (up->status & LEPSEC) { 1189 up->mphase -= SECOND; 1190 if (up->mphase < 0) 1191 up->mphase += MINUTE; 1192 } 1193 } 1194 } 1195 1196 /* 1197 * When the channel metric reaches threshold and the second 1198 * counter matches the minute epoch within the second, the 1199 * driver has synchronized to the station. The second number is 1200 * the remaining seconds until the next minute epoch, while the 1201 * sync epoch is zero. Watch out for the first second; if 1202 * already synchronized to the second, the buffered sync epoch 1203 * must be set. 1204 * 1205 * Note the guard interval is 200 ms; if for some reason the 1206 * clock drifts more than that, it might wind up in the wrong 1207 * second. If the maximum frequency error is not more than about 1208 * 1 PPM, the clock can go as much as two days while still in 1209 * the same second. 1210 */ 1211 if (up->status & MSYNC) { 1212 wwv_epoch(peer); 1213 } else if (up->sptr != NULL) { 1214 sp = up->sptr; 1215 if (sp->metric >= TTHR && epoch == sp->mepoch % SECOND) { 1216 up->rsec = (60 - sp->mepoch / SECOND) % 60; 1217 up->rphase = 0; 1218 up->status |= MSYNC; 1219 up->watch = 0; 1220 if (!(up->status & SSYNC)) 1221 up->repoch = up->yepoch = epoch; 1222 else 1223 up->repoch = up->yepoch; 1224 1225 } 1226 } 1227 1228 /* 1229 * The second sync pulse is extracted using 5-ms (40 sample) FIR 1230 * matched filters at 1000 Hz for WWV or 1200 Hz for WWVH. This 1231 * pulse is used for the most precise synchronization, since if 1232 * provides a resolution of one sample (125 us). The filters run 1233 * only if the station has been reliably determined. 1234 */ 1235 if (up->status & SELV) { 1236 pdelay = (int)(pp->fudgetime1 * SECOND); 1237 mfsync = sqrt(csiamp * csiamp + csqamp * csqamp) / 1238 TCKCYC; 1239 } else if (up->status & SELH) { 1240 pdelay = (int)(pp->fudgetime2 * SECOND); 1241 mfsync = sqrt(hsiamp * hsiamp + hsqamp * hsqamp) / 1242 TCKCYC; 1243 } else { 1244 pdelay = 0; 1245 mfsync = 0; 1246 } 1247 1248 /* 1249 * Enhance the seconds sync pulse using a 1-s (8000-sample) comb 1250 * filter. Correct for the FIR matched filter delay, which is 5 1251 * ms for both the WWV and WWVH filters, and also for the 1252 * propagation delay. Once each second look for second sync. If 1253 * not in minute sync, fiddle the codec gain. Note the SNR is 1254 * computed from the maximum sample and the envelope of the 1255 * sample 6 ms before it, so if we slip more than a cycle the 1256 * SNR should plummet. The signal is scaled to produce unit 1257 * energy at the maximum value. 1258 */ 1259 dtemp = (epobuf[epoch] += (mfsync - epobuf[epoch]) / 1260 up->avgint); 1261 if (dtemp > epomax) { 1262 int j; 1263 1264 epomax = dtemp; 1265 epopos = epoch; 1266 j = epoch - 6 * MS; 1267 if (j < 0) 1268 j += SECOND; 1269 nxtmax = fabs(epobuf[j]); 1270 } 1271 if (epoch == 0) { 1272 up->epomax = epomax; 1273 up->eposnr = wwv_snr(epomax, nxtmax); 1274 epopos -= pdelay + TCKCYC * MS; 1275 if (epopos < 0) 1276 epopos += SECOND; 1277 wwv_endpoc(peer, epopos); 1278 if (!(up->status & SSYNC)) 1279 up->alarm |= SYNERR; 1280 epomax = 0; 1281 if (!(up->status & MSYNC)) 1282 wwv_gain(peer); 1283 } 1284} 1285 1286 1287/* 1288 * wwv_qrz - identify and acquire WWV/WWVH minute sync pulse 1289 * 1290 * This routine implements a virtual station process used to acquire 1291 * minute sync and to mitigate among the ten frequency and station 1292 * combinations. During minute sync acquisition the process probes each 1293 * frequency and station in turn for the minute pulse, which 1294 * involves searching through the entire 480,000-sample minute. The 1295 * process finds the maximum signal and RMS noise plus signal. Then, the 1296 * actual noise is determined by subtracting the energy of the matched 1297 * filter. 1298 * 1299 * Students of radar receiver technology will discover this algorithm 1300 * amounts to a range-gate discriminator. A valid pulse must have peak 1301 * amplitude at least QTHR (2500) and SNR at least QSNR (20) dB and the 1302 * difference between the current and previous epoch must be less than 1303 * AWND (20 ms). Note that the discriminator peak occurs about 800 ms 1304 * into the second, so the timing is retarded to the previous second 1305 * epoch. 1306 */ 1307static void 1308wwv_qrz( 1309 struct peer *peer, /* peer structure pointer */ 1310 struct sync *sp, /* sync channel structure */ 1311 int pdelay /* propagation delay (samples) */ 1312 ) 1313{ 1314 struct refclockproc *pp; 1315 struct wwvunit *up; 1316 char tbuf[80]; /* monitor buffer */ 1317 long epoch; 1318 1319 pp = peer->procptr; 1320 up = (struct wwvunit *)pp->unitptr; 1321 1322 /* 1323 * Find the sample with peak amplitude, which defines the minute 1324 * epoch. Accumulate all samples to determine the total noise 1325 * energy. 1326 */ 1327 epoch = up->mphase - pdelay - SYNSIZ; 1328 if (epoch < 0) 1329 epoch += MINUTE; 1330 if (sp->amp > sp->maxeng) { 1331 sp->maxeng = sp->amp; 1332 sp->pos = epoch; 1333 } 1334 sp->noieng += sp->amp; 1335 1336 /* 1337 * At the end of the minute, determine the epoch of the minute 1338 * sync pulse, as well as the difference between the current and 1339 * previous epoches due to the intrinsic frequency error plus 1340 * jitter. When calculating the SNR, subtract the pulse energy 1341 * from the total noise energy and then normalize. 1342 */ 1343 if (up->mphase == 0) { 1344 sp->synmax = sp->maxeng; 1345 sp->synsnr = wwv_snr(sp->synmax, (sp->noieng - 1346 sp->synmax) / MINUTE); 1347 if (sp->count == 0) 1348 sp->lastpos = sp->pos; 1349 epoch = (sp->pos - sp->lastpos) % MINUTE; 1350 sp->reach <<= 1; 1351 if (sp->reach & (1 << AMAX)) 1352 sp->count--; 1353 if (sp->synmax > ATHR && sp->synsnr > ASNR) { 1354 if (abs(epoch) < AWND * MS) { 1355 sp->reach |= 1; 1356 sp->count++; 1357 sp->mepoch = sp->lastpos = sp->pos; 1358 } else if (sp->count == 1) { 1359 sp->lastpos = sp->pos; 1360 } 1361 } 1362 if (up->watch > ACQSN) 1363 sp->metric = 0; 1364 else 1365 sp->metric = wwv_metric(sp); 1366 if (pp->sloppyclockflag & CLK_FLAG4) { 1367 sprintf(tbuf, 1368 "wwv8 %04x %3d %s %04x %.0f %.0f/%.1f %4ld %4ld", 1369 up->status, up->gain, sp->refid, 1370 sp->reach & 0xffff, sp->metric, sp->synmax, 1371 sp->synsnr, sp->pos % SECOND, epoch); 1372 record_clock_stats(&peer->srcadr, tbuf); 1373#ifdef DEBUG 1374 if (debug) 1375 printf("%s\n", tbuf); 1376#endif /* DEBUG */ 1377 } 1378 sp->maxeng = sp->noieng = 0; 1379 } 1380} 1381 1382 1383/* 1384 * wwv_endpoc - identify and acquire second sync pulse 1385 * 1386 * This routine is called at the end of the second sync interval. It 1387 * determines the second sync epoch position within the second and 1388 * disciplines the sample clock using a frequency-lock loop (FLL). 1389 * 1390 * Second sync is determined in the RF input routine as the maximum 1391 * over all 8000 samples in the second comb filter. To assure accurate 1392 * and reliable time and frequency discipline, this routine performs a 1393 * great deal of heavy-handed heuristic data filtering and grooming. 1394 */ 1395static void 1396wwv_endpoc( 1397 struct peer *peer, /* peer structure pointer */ 1398 int epopos /* epoch max position */ 1399 ) 1400{ 1401 struct refclockproc *pp; 1402 struct wwvunit *up; 1403 static int epoch_mf[3]; /* epoch median filter */ 1404 static int tepoch; /* current second epoch */ 1405 static int xepoch; /* last second epoch */ 1406 static int zepoch; /* last run epoch */ 1407 static int zcount; /* last run end time */ 1408 static int scount; /* seconds counter */ 1409 static int syncnt; /* run length counter */ 1410 static int maxrun; /* longest run length */ 1411 static int mepoch; /* longest run end epoch */ 1412 static int mcount; /* longest run end time */ 1413 static int avgcnt; /* averaging interval counter */ 1414 static int avginc; /* averaging ratchet */ 1415 static int iniflg; /* initialization flag */ 1416 char tbuf[80]; /* monitor buffer */ 1417 double dtemp; 1418 int tmp2; 1419 1420 pp = peer->procptr; 1421 up = (struct wwvunit *)pp->unitptr; 1422 if (!iniflg) { 1423 iniflg = 1; 1424 memset((char *)epoch_mf, 0, sizeof(epoch_mf)); 1425 } 1426 1427 /* 1428 * If the signal amplitude or SNR fall below thresholds, dim the 1429 * second sync lamp and wait for hotter ions. If no stations are 1430 * heard, we are either in a probe cycle or the ions are really 1431 * cold. 1432 */ 1433 scount++; 1434 if (up->epomax < STHR || up->eposnr < SSNR) { 1435 up->status &= ~(SSYNC | FGATE); 1436 avgcnt = syncnt = maxrun = 0; 1437 return; 1438 } 1439 if (!(up->status & (SELV | SELH))) 1440 return; 1441 1442 /* 1443 * A three-stage median filter is used to help denoise the 1444 * second sync pulse. The median sample becomes the candidate 1445 * epoch. 1446 */ 1447 epoch_mf[2] = epoch_mf[1]; 1448 epoch_mf[1] = epoch_mf[0]; 1449 epoch_mf[0] = epopos; 1450 if (epoch_mf[0] > epoch_mf[1]) { 1451 if (epoch_mf[1] > epoch_mf[2]) 1452 tepoch = epoch_mf[1]; /* 0 1 2 */ 1453 else if (epoch_mf[2] > epoch_mf[0]) 1454 tepoch = epoch_mf[0]; /* 2 0 1 */ 1455 else 1456 tepoch = epoch_mf[2]; /* 0 2 1 */ 1457 } else { 1458 if (epoch_mf[1] < epoch_mf[2]) 1459 tepoch = epoch_mf[1]; /* 2 1 0 */ 1460 else if (epoch_mf[2] < epoch_mf[0]) 1461 tepoch = epoch_mf[0]; /* 1 0 2 */ 1462 else 1463 tepoch = epoch_mf[2]; /* 1 2 0 */ 1464 } 1465 1466 1467 /* 1468 * If the epoch candidate is the same as the last one, increment 1469 * the run counter. If not, save the length, epoch and end 1470 * time of the current run for use later and reset the counter. 1471 * The epoch is considered valid if the run is at least SCMP 1472 * (10) s, the minute is synchronized and the interval since the 1473 * last epoch is not greater than the averaging interval. Thus, 1474 * after a long absence, the program will wait a full averaging 1475 * interval while the comb filter charges up and noise 1476 * dissapates.. 1477 */ 1478 tmp2 = (tepoch - xepoch) % SECOND; 1479 if (tmp2 == 0) { 1480 syncnt++; 1481 if (syncnt > SCMP && up->status & MSYNC && (up->status & 1482 FGATE || scount - zcount <= up->avgint)) { 1483 up->status |= SSYNC; 1484 up->yepoch = tepoch; 1485 } 1486 } else if (syncnt >= maxrun) { 1487 maxrun = syncnt; 1488 mcount = scount; 1489 mepoch = xepoch; 1490 syncnt = 0; 1491 } 1492 if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status & MSYNC)) 1493 { 1494 sprintf(tbuf, 1495 "wwv1 %04x %3d %4d %5.0f %5.1f %5d %4d %4d %4d", 1496 up->status, up->gain, tepoch, up->epomax, 1497 up->eposnr, tmp2, avgcnt, syncnt, 1498 maxrun); 1499 record_clock_stats(&peer->srcadr, tbuf); 1500#ifdef DEBUG 1501 if (debug) 1502 printf("%s\n", tbuf); 1503#endif /* DEBUG */ 1504 } 1505 avgcnt++; 1506 if (avgcnt < up->avgint) { 1507 xepoch = tepoch; 1508 return; 1509 } 1510 1511 /* 1512 * The sample clock frequency is disciplined using a first-order 1513 * feedback loop with time constant consistent with the Allan 1514 * intercept of typical computer clocks. During each averaging 1515 * interval the candidate epoch at the end of the longest run is 1516 * determined. If the longest run is zero, all epoches in the 1517 * interval are different, so the candidate epoch is the current 1518 * epoch. The frequency update is computed from the candidate 1519 * epoch difference (125-us units) and time difference (seconds) 1520 * between updates. 1521 */ 1522 if (syncnt >= maxrun) { 1523 maxrun = syncnt; 1524 mcount = scount; 1525 mepoch = xepoch; 1526 } 1527 xepoch = tepoch; 1528 if (maxrun == 0) { 1529 mepoch = tepoch; 1530 mcount = scount; 1531 } 1532 1533 /* 1534 * The master clock runs at the codec sample frequency of 8000 1535 * Hz, so the intrinsic time resolution is 125 us. The frequency 1536 * resolution ranges from 18 PPM at the minimum averaging 1537 * interval of 8 s to 0.12 PPM at the maximum interval of 1024 1538 * s. An offset update is determined at the end of the longest 1539 * run in each averaging interval. The frequency adjustment is 1540 * computed from the difference between offset updates and the 1541 * interval between them. 1542 * 1543 * The maximum frequency adjustment ranges from 187 PPM at the 1544 * minimum interval to 1.5 PPM at the maximum. If the adjustment 1545 * exceeds the maximum, the update is discarded and the 1546 * hysteresis counter is decremented. Otherwise, the frequency 1547 * is incremented by the adjustment, but clamped to the maximum 1548 * 187.5 PPM. If the update is less than half the maximum, the 1549 * hysteresis counter is incremented. If the counter increments 1550 * to +3, the averaging interval is doubled and the counter set 1551 * to zero; if it decrements to -3, the interval is halved and 1552 * the counter set to zero. 1553 */ 1554 dtemp = (mepoch - zepoch) % SECOND; 1555 if (up->status & FGATE) { 1556 if (abs(dtemp) < MAXFREQ * MINAVG) { 1557 up->freq += (dtemp / 2.) / ((mcount - zcount) * 1558 FCONST); 1559 if (up->freq > MAXFREQ) 1560 up->freq = MAXFREQ; 1561 else if (up->freq < -MAXFREQ) 1562 up->freq = -MAXFREQ; 1563 if (abs(dtemp) < MAXFREQ * MINAVG / 2.) { 1564 if (avginc < 3) { 1565 avginc++; 1566 } else { 1567 if (up->avgint < MAXAVG) { 1568 up->avgint <<= 1; 1569 avginc = 0; 1570 } 1571 } 1572 } 1573 } else { 1574 if (avginc > -3) { 1575 avginc--; 1576 } else { 1577 if (up->avgint > MINAVG) { 1578 up->avgint >>= 1; 1579 avginc = 0; 1580 } 1581 } 1582 } 1583 } 1584 if (pp->sloppyclockflag & CLK_FLAG4) { 1585 sprintf(tbuf, 1586 "wwv2 %04x %5.0f %5.1f %5d %4d %4d %4d %4.0f %7.2f", 1587 up->status, up->epomax, up->eposnr, mepoch, 1588 up->avgint, maxrun, mcount - zcount, dtemp, 1589 up->freq * 1e6 / SECOND); 1590 r…
Large files files are truncated, but you can click here to view the full file