+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdlib.h>
+
+double frac(double x)
+{
+ return x - (int)x;
+}
+
+double filter(double x)
+{
+ static double l = 0.0;
+ l = 0.001 * x + 0.999 * l;
+ return l;
+}
+
+int main(int argc, char **argv)
+{
+ FILE *in1, *in2, *skew;
+
+ // open input1 (reference)
+ if (strcmp(argv[1], "-") == 0) {
+ in1 = stdin;
+ } else {
+ in1 = fopen(argv[1], "rb");
+ if (in1 == NULL) {
+ perror(argv[1]);
+ exit(1);
+ }
+ }
+
+ // open input2
+ if (strcmp(argv[2], "-") == 0) {
+ in2 = stdin;
+ } else {
+ in2 = fopen(argv[2], "rb");
+ if (in2 == NULL) {
+ perror(argv[2]);
+ exit(1);
+ }
+ }
+
+ // open (estimated) skew
+ if (strcmp(argv[3], "-") == 0) {
+ skew = stderr;
+ } else {
+ skew = fopen(argv[3], "wb");
+ if (skew == NULL) {
+ perror(argv[3]);
+ exit(1);
+ }
+ }
+
+ short prev_sample, sample = 0;
+ int in_pos = -1;
+ double p = -1.0;
+ double speed = 1.0;
+
+ while (!feof(in1) && !feof(in2)) {
+ // read sample from reference
+ if (fread(&sample, sizeof(short), 1, in1) != 1) {
+ exit(0);
+ }
+
+ // read sample from in2 (at current speed)
+ p += speed;
+ while ((int)(ceil(p)) > in_pos) {
+ prev_sample = sample;
+ if (fread(&sample, sizeof(short), 1, in2) != 1) {
+ exit(0);
+ }
+ ++in_pos;
+ }
+
+ // linear interpolation (works well since delta_p varies so slowly)
+ double t = frac(p);
+ short intp_sample = prev_sample * (1.0 - t) + sample * t;
+
+ // subtract the two samples and pull it through a filter. we assume
+ // the sines are normalized for now (and that there's no dc skew)
+ double offset = -0.001 * filter(intp_sample - sample);
+ speed += offset;
+ fwrite(&speed, sizeof(double), 1, skew);
+
+ printf("%f (offset=%f)\n", speed, offset);
+ }
+}