/contrib/ntp/ntpd/refclock_wwv.c

https://bitbucket.org/freebsd/freebsd-head/ · C · 2709 lines · 1594 code · 204 blank · 911 comment · 264 complexity · d5531729a944f68437d2d9bd4047b8da MD5 · raw file

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